# DTIlib

Autores: guilhermesaraiva e gustavopinheiro 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]
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)
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]
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
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:
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)
8
9 y, x = vector_field(lambda y,x: x**3, lambda y,x: -y)
14
15 y, x = vector_field(lambda y,x: x**3, lambda y,x: -y**3)
20
21 y, x = vector_field(lambda y,x: sin(2*x)+sin(2*y), lambda y,x: sin(2*x)+sin(2*y))
26
27 y, x = vector_field(lambda y,x: np.tanh(y), lambda y,x: np.tanh(x))
32
33 y, x = vector_field(lambda y,x: np.tanh(x), lambda y,x: -np.tanh(y))
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))
44
45 y, x = vector_field(lambda y,x: y, lambda y,x: x)
50
51 y, x = vector_field(lambda y,x: x, lambda y,x: -y)
56
57 y, x = vector_field(lambda y,x: x+1+y, lambda y,x: x-1-y)
62
63
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,:]
73
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,:]
83
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,:]
93
94
```1 pass