DTIlib

Autores: guilhermesaraiva e gustavopinheiro
Data: 14/06/2015

Biblioteca para análise de imagens DTI.

  1 import numpy as np
  2 import cv2
  3 import ia636 as ia
  4 
  5 
  6 def aumenta(I, kh, kw):
  7 
  8     (h,w) = I.shape
  9     H = np.zeros((kh*h,kw*w))
 10     H[::kh,::kw] = I
 11     k = np.ones((kh,kw))
 12     H = ia.iapconv(H, k)
 13     return H
 14 
 15 def resize(I, sx, sy=0):
 16     if sy == 0:
 17         sy = sx
 18     if I.ndim == 3:
 19         I = np.rollaxis(I, 0, 3)
 20         I = cv2.resize(I, dsize=(0,0), fx=sx, fy=sy, interpolation=cv2.INTER_NEAREST)
 21         return np.rollaxis(I, 2)
 22     return cv2.resize(I, dsize=(0,0), fx=sx, fy=sy, interpolation=cv2.INTER_NEAREST)
 23 
 24 def normaliza(I):
 25     Imin = np.amin(I)
 26     Imax = np.amax(I)
 27     return (255 * (I-Imin) / (Imax-Imin)).astype(np.uint8)
 28 
 29 def normaliza_cor(I):
 30     Ir = np.zeros(I.shape)
 31     Ig = np.zeros(I.shape)
 32     Ib = np.zeros(I.shape)
 33     Imax = np.amax(np.abs(I))
 34     Ipos = I>0
 35     Ineg = I<0
 36     Ir[Ipos] = ( 255 * I[Ipos]).astype(np.uint8)
 37     Ib[Ineg] = (-255 * I[Ineg]).astype(np.uint8)
 38     Ir = np.expand_dims(Ir, 0)
 39     Ig = np.expand_dims(Ig, 0)
 40     Ib = np.expand_dims(Ib, 0)
 41     return np.concatenate((Ir,Ig,Ib))
 42 
 43 def interpola(w, axis):
 44 
 45     W = (w[:,:,:,:-1]+w[:,:,:,1:])/2
 46 
 47     return W
 48 
 49 def teste(v, axis):
 50 
 51     #produto escalar
 52     V = np.sum(v[:,:,:,:-1]*v[:,:,:,1:], axis=0)
 53 
 54     #interpolacao do campo vetorial e ponderacao pelo produto escalar
 55     V = V*interpola(v,1)
 56 
 57     return V
 58 
 59 
 60 def set_direction(v, direction):
 61     direction = direction.reshape(3,1,1,1)
 62     dot_product = np.sum(direction*v, axis=0)
 63     s = np.sign(dot_product)
 64     s[s==0] = 1
 65     return s*v
 66 
 67 def set_direction2(v, direction):
 68     direction = direction/(np.linalg.norm(direction))
 69     direction = direction.reshape(3,1,1,1)
 70     dot_product = np.sum(direction*v, axis=0)
 71     s = dot_product
 72     return s*v
 73 
 74 def divergente(vy, vx = None, ks = 3):
 75 
 76     if ks == 2:
 77         k1 = np.array([[1, -1]])
 78 
 79         if vx == None:
 80 
 81             vx = ia.iaconv(vy,k1)[:,:-1]
 82             vy = ia.iaconv(vy,k1.T)[:-1,:]
 83 
 84         dx1 = ia.iaconv(vx,k1)[:,:-1]
 85         dy1 = ia.iaconv(vy,k1.T)[:-1,:]
 86 
 87 
 88         d1=dx1+dy1
 89 
 90         return d1[1:-1,1:-1]
 91 
 92     else:
 93         k2 = np.array([[1, 0, -1]])
 94 
 95         if vx == None:
 96 
 97             vx = ia.iaconv(vy,k2)[:,1:-1]
 98             vy = ia.iaconv(vy,k2.T)[1:-1,:]
 99 
100         dx2 = ia.iaconv(vx,k2)[:,1:-1]
101         dy2 = ia.iaconv(vy,k2.T)[1:-1,:]
102 
103         d2=dx2+dy2
104 
105         return d2[1:-1,1:-1]
106 
107     return dd2[1:-1,1:-1]
108 
109 
110 def rotacional(vy, vx = None, ks = 3):
111 
112     if ks == 2:
113         k1 = np.array([[1, -1]])
114 
115         if vx == None:
116 
117             vx = ia.iaconv(vy,k1)[:,:-1]
118             vy = ia.iaconv(vy,k1.T)[:-1,:]
119 
120         dxy = ia.iaconv(vx,k1.T)[:-1,:]
121         dyx = ia.iaconv(vy,k1)[:,:-1]
122 
123 
124         r=dyx - dxy
125 
126         return r[1:-1,1:-1]
127 
128     else:
129         k2 = np.array([[1, 0, -1]])
130 
131         if vx == None:
132 
133             vx = ia.iaconv(vy,k2)[:,1:-1]
134             vy = ia.iaconv(vy,k2.T)[1:-1,:]
135 
136         dxy = ia.iaconv(vx,k2.T)[1:-1,:]
137         dyx = ia.iaconv(vy,k2)[:,1:-1]
138 
139 
140         r=dyx - dxy
141 
142         return r[1:-1,1:-1]
143 
144     return r[1:-1,1:-1]
145 
146 def vector_field(fy, fx, domain=(-2.,2.,-2.,2.), size=(512,512)):
147     ymin, ymax, xmin, xmax = domain
148     H,W = size
149     x = np.linspace(xmin, xmax, W)
150     y = np.linspace(ymin, ymax, H)
151     xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
152     fx_v = np.vectorize(fx)
153     fy_v = np.vectorize(fy)
154     return (fy_v(yv,xv), fx_v(yv,xv))
155 
156 def show_magnitude(y, x):
157     return ia.ianormalize((x**2+y**2)**0.5)
158 
159 def show_vector_field(y, x, step=32, maxlen=24, rescale=1, showPoints=False, color=(0,0,255), bgcolor=(255,255,255), pointsColor=(255,0,0), precision=3):
160     H, W = x.shape
161     img = np.array(bgcolor, np.uint8).reshape((1,1,3)) * np.ones((rescale*H, rescale*W, 3), np.uint8)
162     mask_y, mask_x = np.mgrid[step/2 : H : step, step/2 : W : step]
163     ssx = x[mask_y, mask_x]
164     ssy = y[mask_y, mask_x]
165     mag = np.sqrt(ssx**2+ssy**2)
166     magmax = np.max(mag)
167     cx = ssx / magmax
168     cy = ssy / magmax
169     dx = ssx / mag
170     dy = ssy / mag
171     dx[np.isnan(dx)] = 0
172     dy[np.isnan(dy)] = 0
173     p = 2**precision
174     pt1_x = (p*(rescale*mask_x - 0.5*maxlen*cx)).astype(int)
175     pt1_y = (p*(rescale*mask_y - 0.5*maxlen*cy)).astype(int)
176     pt2_x = (p*(rescale*mask_x + 0.5*maxlen*cx)).astype(int)
177     pt2_y = (p*(rescale*mask_y + 0.5*maxlen*cy)).astype(int)
178     pt3_x = ((pt2_x - p*maxlen*(0.2*dx-0.1*dy))).astype(int)
179     pt3_y = ((pt2_y - p*maxlen*(0.2*dy+0.1*dx))).astype(int)
180     pt4_x = ((pt2_x - p*maxlen*(0.2*dx+0.1*dy))).astype(int)
181     pt4_y = ((pt2_y - p*maxlen*(0.2*dy-0.1*dx))).astype(int)
182     for i in np.arange(mask_y.shape[0]):
183         for j in np.arange(mask_x.shape[1]):
184             cv2.line(img, (pt1_x[i,j], pt1_y[i,j]), (pt2_x[i,j], pt2_y[i,j]), color, lineType=cv2.CV_AA, shift=precision)
185             cv2.line(img, (pt2_x[i,j], pt2_y[i,j]), (pt3_x[i,j], pt3_y[i,j]), color, lineType=cv2.CV_AA, shift=precision)
186             cv2.line(img, (pt2_x[i,j], pt2_y[i,j]), (pt4_x[i,j], pt4_y[i,j]), color, lineType=cv2.CV_AA, shift=precision)
187 
188     if showPoints:
189         from scipy.ndimage.measurements import minimum_position
190         mag2 = x**2+y**2
191         min_pos = minimum_position(mag2)
192         if type(min_pos) == tuple:
193             min_pos = [min_pos]
194         for i,j in min_pos:
195             cv2.circle(img, (rescale*j, rescale*i), 3, pointsColor, thickness=-1)
196 
197     img = np.swapaxes(np.swapaxes(img,0,2),1,2)
198 
199     return img
200 
201 
202 def show_vector_field_image(y, x, f, step=32, color=(0,0,255), maxlen=24, showPoints=True, pointsColor=(255,0,0), precision=3):
203     H, W = x.shape
204     img = np.empty((H,W,3), np.uint8)
205 
206     img[:,:,:] = f.reshape(H,W,1)
207 
208     mask_x, mask_y = np.mgrid[step/2 : step/2+H : step, step/2 : step/2+W : step]
209     ssx = x[mask_y, mask_x]
210     ssy = y[mask_y, mask_x]
211     mag = np.sqrt(ssx**2+ssy**2)
212     magmax = np.max(mag)
213     cx = ssx / magmax
214     cy = ssy / magmax
215     dx = ssx / mag
216     dy = ssy / mag
217     dx[np.isnan(dx)] = 0
218     dy[np.isnan(dy)] = 0
219     p = 2**precision
220     pt1_x = (p*(mask_x - 0.5*maxlen*cx)).astype(int)
221     pt1_y = (p*(mask_y - 0.5*maxlen*cy)).astype(int)
222     pt2_x = (p*(mask_x + 0.5*maxlen*cx)).astype(int)
223     pt2_y = (p*(mask_y + 0.5*maxlen*cy)).astype(int)
224     pt3_x = ((pt2_x - p*maxlen*(0.2*dx-0.1*dy))).astype(int)
225     pt3_y = ((pt2_y - p*maxlen*(0.2*dy+0.1*dx))).astype(int)
226     pt4_x = ((pt2_x - p*maxlen*(0.2*dx+0.1*dy))).astype(int)
227     pt4_y = ((pt2_y - p*maxlen*(0.2*dy-0.1*dx))).astype(int)
228     p = 2**precision
229     for i in np.arange(mask_y.shape[0]):
230         for j in np.arange(mask_x.shape[1]):
231             cv2.line(img, (pt1_x[i,j], pt1_y[i,j]), (pt2_x[i,j], pt2_y[i,j]), color, lineType=cv2.CV_AA, shift=precision)
232             cv2.line(img, (pt2_x[i,j], pt2_y[i,j]), (pt3_x[i,j], pt3_y[i,j]), color, lineType=cv2.CV_AA, shift=precision)
233             cv2.line(img, (pt2_x[i,j], pt2_y[i,j]), (pt4_x[i,j], pt4_y[i,j]), color, lineType=cv2.CV_AA, shift=precision)
234 
235     if showPoints:
236         from scipy.ndimage.measurements import minimum_position
237         mag2 = x**2+y**2
238         min_pos = minimum_position(mag2)
239         if type(min_pos) == tuple:
240             min_pos = [min_pos]
241         for i,j in min_pos:
242             cv2.circle(img, (j,i), 3, pointsColor, thickness=-1)
243 
244     img = np.swapaxes(np.swapaxes(img,0,2),1,2)
245 
246 
247     return img
  1 import ia636 as ia
  2 import numpy as np
  3 import cv2
  4 
  5 def vector_field3(fz, fy, fx, domain=(-2.,2.,-2.,2.,-2.,2.), size=(80,80,80)):
  6     zmin, zmax, ymin, ymax, xmin, xmax = domain
  7     D,H,W = size
  8     x = np.linspace(xmin, xmax, W)
  9     y = np.linspace(ymin, ymax, H)
 10     z = np.linspace(zmin, zmax, D)
 11     xv, yv, zv = np.meshgrid(x, y, z, sparse=False, indexing='xy')
 12     fx_v = np.vectorize(fx)
 13     fy_v = np.vectorize(fy)
 14     fz_v = np.vectorize(fz)
 15     return (fz_v(zv,yv,xv), fy_v(zv,yv,xv), fx_v(zv,yv,xv))
 16 
 17 def gradiente3(v, ks = 3):
 18     if v.ndim == 2:
 19         v = np.expand_dims(v, 0)
 20     if ks == 2:
 21         k = np.array([[[1, -1]]])
 22         vx = ia.iaconv(v,k)[:,:,:-1]
 23         vy = ia.iaconv(v,k.swapaxes(1,2))[:,:-1,:]
 24         vz = ia.iaconv(v,k.swapaxes(0,2))[:-1,:,:]
 25     else:
 26         k = np.array([[[1, 0, -1]]])
 27         vx = ia.iaconv(v,k)[:,:,1:-1]
 28         vy = ia.iaconv(v,k.swapaxes(1,2))[:,1:-1,:]
 29         vz = ia.iaconv(v,k.swapaxes(0,2))[1:-1,:,:]
 30     return (vz, vy, vx)
 31 
 32 
 33 def divergente3(v, ks = 3):
 34 
 35     if type(v) != tuple:
 36         v = (v,)
 37 
 38     n = len(v)
 39 
 40     # Se a imagem de entrada for um campo escalar
 41     if n == 1:
 42         v = gradiente3(v[0], ks)
 43         n = 3
 44 
 45     # CASO BIDIMENSIONAL
 46     if n == 2:
 47         vy, vx = v
 48         vx = np.expand_dims(vx, 0)
 49         vy = np.expand_dims(vy, 0)
 50         vz = np.zeros(vx.shape, vx.dtype)
 51         v = (vz, vy, vx)
 52         n = 3
 53 
 54     # CASO TRIDIMENSIONAL
 55     if n == 3:
 56 
 57         vz, vy, vx = v
 58 
 59         if ks == 2:
 60             k = np.array([[[1, -1]]])
 61 
 62             dx = ia.iaconv(vx,k)[:,:,:-1]
 63             dy = ia.iaconv(vy,k.swapaxes(1,2))[:,:-1,:]
 64             dz = ia.iaconv(vz,k.swapaxes(0,2))[:-1,:,:]
 65 
 66         else:
 67             k = np.array([[[1, 0, -1]]])
 68 
 69             dx = ia.iaconv(vx,k)[:,:,1:-1]
 70             dy = ia.iaconv(vy,k.swapaxes(1,2))[:,1:-1,:]
 71             dz = ia.iaconv(vz,k.swapaxes(0,2))[1:-1,:,:]
 72 
 73         d=dx+dy+dz
 74 
 75         return d
 76 
 77 def rotacional3(v, ks = 3):
 78 
 79     if type(v) != tuple:
 80         v = (v,)
 81 
 82     n = len(v)
 83 
 84     # CASO BIDIMENSIONAL
 85     if n == 2:
 86         vy, vx = v
 87         vx = np.expand_dims(vx, 0)
 88         vy = np.expand_dims(vy, 0)
 89         vz = np.zeros(vx.shape, vx.dtype)
 90         v = (vz, vy, vx)
 91         n = 3
 92 
 93     # CASO TRIDIMENSIONAL
 94     if n == 3:
 95 
 96         vz, vy, vx = v
 97 
 98         if ks == 2:
 99             kx = np.array([[[1, -1]]])
100             ky = kx.swapaxes(1,2)
101             kz = kx.swapaxes(0,2)
102 
103             dxy = ia.iaconv(vx,ky)[:,:-1,:]
104             dxz = ia.iaconv(vx,kz)[:-1,:,:]
105             dyx = ia.iaconv(vy,kx)[:,:,:-1]
106             dyz = ia.iaconv(vy,kz)[:-1,:,:]
107             dzx = ia.iaconv(vz,kx)[:,:,:-1]
108             dzy = ia.iaconv(vz,ky)[:,:-1,:]
109 
110         else:
111             kx = np.array([[[1, 0, -1]]])
112             ky = kx.swapaxes(1,2)
113             kz = kx.swapaxes(0,2)
114 
115             dxy = ia.iaconv(vx,ky)[:,1:-1,:]
116             dxz = ia.iaconv(vx,kz)[1:-1,:,:]
117             dyx = ia.iaconv(vy,kx)[:,:,1:-1]
118             dyz = ia.iaconv(vy,kz)[1:-1,:,:]
119             dzx = ia.iaconv(vz,kx)[:,:,1:-1]
120             dzy = ia.iaconv(vz,ky)[:,1:-1,:]
121 
122         rx = dzy - dyz
123         ry = dxz - dzx
124         rz = dyx - dxy
125 
126         return (rz, ry, rx)
  1 from DTIlib import *
  2 
  3 y, x = vector_field(lambda y,x: y, lambda y,x: x**3)
  4 adshow(show_vector_field(y,x), 'Campo vetorial')
  5 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
  6 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
  7 adshow(show_magnitude(y,x), 'Magnitude do campo')
  8 
  9 y, x = vector_field(lambda y,x: x**3, lambda y,x: -y)
 10 adshow(show_vector_field(y,x), 'Campo vetorial')
 11 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 12 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 13 adshow(show_magnitude(y,x), 'Magnitude do campo')
 14 
 15 y, x = vector_field(lambda y,x: x**3, lambda y,x: -y**3)
 16 adshow(show_vector_field(y,x), 'Campo vetorial')
 17 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 18 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 19 adshow(show_magnitude(y,x), 'Magnitude do campo')
 20 
 21 y, x = vector_field(lambda y,x: sin(2*x)+sin(2*y), lambda y,x: sin(2*x)+sin(2*y))
 22 adshow(show_vector_field(y,x), 'Campo vetorial')
 23 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 24 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 25 adshow(show_magnitude(y,x), 'Magnitude do campo')
 26 
 27 y, x = vector_field(lambda y,x: np.tanh(y), lambda y,x: np.tanh(x))
 28 adshow(show_vector_field(y,x), 'Campo vetorial')
 29 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 30 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 31 adshow(show_magnitude(y,x), 'Magnitude do campo')
 32 
 33 y, x = vector_field(lambda y,x: np.tanh(x), lambda y,x: -np.tanh(y))
 34 adshow(show_vector_field(y,x), 'Campo vetorial')
 35 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 36 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 37 adshow(show_magnitude(y,x), 'Magnitude do campo')
 38 
 39 y, x = vector_field(lambda y,x: np.tanh(y)+np.tanh(x+1), lambda y,x: np.tanh(x-1)-np.tanh(y))
 40 adshow(show_vector_field(y,x), 'Campo vetorial')
 41 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 42 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 43 adshow(show_magnitude(y,x), 'Magnitude do campo')
 44 
 45 y, x = vector_field(lambda y,x: y, lambda y,x: x)
 46 adshow(show_vector_field(y,x), 'Campo vetorial')
 47 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 48 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 49 adshow(show_magnitude(y,x), 'Magnitude do campo')
 50 
 51 y, x = vector_field(lambda y,x: x, lambda y,x: -y)
 52 adshow(show_vector_field(y,x), 'Campo vetorial')
 53 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 54 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 55 adshow(show_magnitude(y,x), 'Magnitude do campo')
 56 
 57 y, x = vector_field(lambda y,x: x+1+y, lambda y,x: x-1-y)
 58 adshow(show_vector_field(y,x), 'Campo vetorial')
 59 adshow(ia.ianormalize(divergente(y,x)), 'Divergente')
 60 adshow(ia.ianormalize(rotacional(y,x)), 'Rotacional')
 61 adshow(show_magnitude(y,x), 'Magnitude do campo')
 62 
 63 
 64 f = adreadgray('lenina.tif')
 65 f = cv2.GaussianBlur(f, (0,0), 20, borderType=cv2.BORDER_REPLICATE)
 66 k2 = np.array([[1, 0, -1]])
 67 x = ia.iaconv(f,k2)[:,1:-1]
 68 y = ia.iaconv(f,k2.T)[1:-1,:]
 69 adshow(show_vector_field_image(y,x,f,8), 'Campo vetorial')
 70 adshow(ia.ianormalize(divergente(f)[1:-1,1:-1]), 'Divergente')
 71 adshow(ia.ianormalize(rotacional(f)[1:-1,1:-1]), 'Rotacional')
 72 adshow(show_magnitude(y,x), 'Magnitude do campo')
 73 
 74 f = adreadgray('lenina.tif')
 75 f = cv2.GaussianBlur(f, (0,0), 10, borderType=cv2.BORDER_REPLICATE)
 76 k2 = np.array([[1, 0, -1]])
 77 x = ia.iaconv(f,k2)[:,1:-1]
 78 y = ia.iaconv(f,k2.T)[1:-1,:]
 79 adshow(show_vector_field_image(y,x,f,8), 'Campo vetorial')
 80 adshow(ia.ianormalize(divergente(f)[1:-1,1:-1]), 'Divergente')
 81 adshow(ia.ianormalize(rotacional(f)[1:-1,1:-1]), 'Rotacional')
 82 adshow(show_magnitude(y,x), 'Magnitude do campo')
 83 
 84 f = adreadgray('cameraman.tif')
 85 f = cv2.GaussianBlur(f, (0,0), 20, borderType=cv2.BORDER_REPLICATE)
 86 k2 = np.array([[1, 0, -1]])
 87 x = ia.iaconv(f,k2)[:,1:-1]
 88 y = ia.iaconv(f,k2.T)[1:-1,:]
 89 adshow(show_vector_field_image(y,x,f,8), 'Campo vetorial')
 90 adshow(ia.ianormalize(divergente(f)[1:-1,1:-1]), 'Divergente')
 91 adshow(ia.ianormalize(rotacional(f)[1:-1,1:-1]), 'Rotacional')
 92 adshow(show_magnitude(y,x), 'Magnitude do campo')
 93 
 94 
 95 f = adreadgray('cameraman.tif')
 96 f = cv2.GaussianBlur(f, (0,0), 5, borderType=cv2.BORDER_REPLICATE)
 97 k2 = np.array([[1, 0, -1]])
 98 x = ia.iaconv(f,k2)[:,1:-1]
 99 y = ia.iaconv(f,k2.T)[1:-1,:]
100 adshow(show_vector_field_image(y,x,f,8), 'Campo vetorial')
101 adshow(ia.ianormalize(divergente(f)[1:-1,1:-1]), 'Divergente')
102 adshow(ia.ianormalize(rotacional(f)[1:-1,1:-1]), 'Rotacional')
103 adshow(show_magnitude(y,x), 'Magnitude do campo')

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

Campo vetorial

Divergente

Rotacional

Magnitude do campo

1 pass