# drafts2_gustavopinheiro

Autor: gustavopinheiro 28/05/2015

# divergente e rotacional

teste para o projeto

```  1 import ia636 as ia
2 import numpy as np
3 import cv2
4 from DTIlib import *
5 #import dtimp
6 #from nibabel import nifti1
7 #from glob import glob
8
10 #print L.shape
11
12
13 evt_path = find_attachment_file('evt.npy')
14 evl_path = find_attachment_file('evl.npy')
15 fa_path = find_attachment_file('fa.npy')
16
17 #evt_path = find_attachment_file('evt005.npy')
18 #evl_path = find_attachment_file('evl005.npy')
19 #fa_path = find_attachment_file('fa005.npy')
20
21
25
26 # Altera sentido do campo vetorial
27 #s = np.sign(evt[1,:,:,:])
28 #s[s==0] = 1
29 #evt = s * evt
30 # 0: Plano sagital
31 # 1: Plano transversal
32 # 2: Plano coronal
33 #direction = np.array([1,0,0])
34 #direction = np.array([0,1,0])
35 direction = np.array([0,0,1])
36 #direction = np.array([1,20,0])
37 #direction = np.array([tan(3),1,0])
38
39 evt = set_direction(evt, direction)
40 #evt = set_direction2(evt, direction)
41
42 print 'evt.shape =', evt.shape
43 print 'evl.shape =', evl.shape
44 print 'fa.shape =', fa.shape
45
46 print 'np.max(evt[0,:,:,:]) =', np.max(evt[0,:,:,:])
47 print 'np.max(evt[1,:,:,:]) =', np.max(evt[1,:,:,:])
48 print 'np.max(evt[2,:,:,:]) =', np.max(evt[2,:,:,:])
49 print 'np.min(evt[0,:,:,:]) =', np.min(evt[0,:,:,:])
50 print 'np.min(evt[1,:,:,:]) =', np.min(evt[1,:,:,:])
51 print 'np.min(evt[2,:,:,:]) =', np.min(evt[2,:,:,:])
52 print 'np.max(evl) =', np.max(evl)
53 print 'np.min(evl) =', np.min(evl)
54 print 'np.max(fa) =', np.max(fa)
55 print 'np.min(fa) =', np.min(fa)
56
60
64
68
69 C, D, H, W = evt.shape
70 vf_evl = evl * evt
71 vf_fa = fa * evt
72 vf_fa_evl = vf_evl * fa
73 div_evt = divergente3(tuple(evt))
74 div_evl = divergente3(tuple(vf_evl))
75 div_fa = divergente3(tuple(vf_fa))
76 rot_evt = rotacional3(tuple(evt))
77 rot_evl = rotacional3(tuple(vf_evl))
78 rot_fa = rotacional3(tuple(vf_fa))
79
80 print 'np.min(vf_evl) =', np.min(vf_evl)
81 print 'np.max(vf_evl) =', np.max(vf_evl)
82 print 'np.min(vf_fa) =', np.min(vf_fa)
83 print 'np.max(vf_fa) =', np.max(vf_fa)
84 print 'np.min(vf_fa_evl) =', np.min(vf_fa_evl)
85 print 'np.max(vf_fa_evl) =', np.max(vf_fa_evl)
86 print 'np.min(div_evt) =', np.min(div_evt)
87 print 'np.max(div_evt) =', np.max(div_evt)
88 print 'np.min(div_evl) =', np.min(div_evl)
89 print 'np.max(div_evl) =', np.max(div_evl)
90 print 'np.min(div_fa) =', np.min(div_fa)
91 print 'np.max(div_fa) =', np.max(div_fa)
92 print 'np.min(rot_evt) =', np.min(rot_evt)
93 print 'np.max(rot_evt) =', np.max(rot_evt)
94 print 'np.min(rot_evl) =', np.min(rot_evl)
95 print 'np.max(rot_evl) =', np.max(rot_evl)
96 print 'np.min(rot_fa) =', np.min(rot_fa)
97 print 'np.max(rot_fa) =', np.max(rot_fa)
98
99 # Normaliza rotacionais (para visualização das direções apenas)
100 #rot_evl = rot_evl/(rot_evl[0]**2 + rot_evl[1]**2 + rot_evl[2]**2)**0.5
101 #rot_fa = rot_fa/(rot_fa[0]**2 + rot_fa[1]**2 + rot_fa[2]**2)**0.5
102
103 # Normaliza globalmente para exibição
104 div_evt_n = normaliza(div_evt/np.max(div_evt))
105 div_evl_n = normaliza(div_evl/np.max(div_evl))
106 #div_fa_n = normaliza_cor(div_fa/np.max(div_fa))
107 div_fa_n = normaliza(div_fa/np.max(div_fa))       #grayscale
108
109 #print div_fa_n[:,:,0,:].shape
110 rot0_evt_n = normaliza(rot_evt[0]/np.max(rot_evt[0]))
111 rot1_evt_n = normaliza(rot_evt[1]/np.max(rot_evt[1]))
112 rot2_evt_n = normaliza(rot_evt[2]/np.max(rot_evt[2]))
113 rot0_evl_n = normaliza(rot_evl[0]/np.max(rot_evl[0]))
114 rot1_evl_n = normaliza(rot_evl[1]/np.max(rot_evl[1]))
115 rot2_evl_n = normaliza(rot_evl[2]/np.max(rot_evl[2]))
116 rot0_fa_n = normaliza(rot_fa[0]/np.max(rot_fa[0]))
117 rot1_fa_n = normaliza(rot_fa[1]/np.max(rot_fa[1]))
118 rot2_fa_n = normaliza(rot_fa[2]/np.max(rot_fa[2]))
119
120 #k=1
121 #adshow(show_vector_field(evt[0,:,k,:], evt[2,:,k,:], step=1, maxlen=24, rescale=16), 'projeção do campo vetorial unitário (plano transversal, fatia %s)' % k)
122
123 for k in np.arange(H):
124     adshow(show_vector_field(evt[0,:,k,:], evt[2,:,k,:], step=1, maxlen=24, rescale=16), 'projeção do campo vetorial unitário (plano transversal, fatia %s)' % k)
125     adshow(show_vector_field(vf_evl[0,:,k,:], vf_evl[2,:,k,:], step=1, maxlen=24, rescale=16), 'proj. do c. v. ponderado pelos autovalores (plano transversal, fatia %s)' % k)
126     adshow(show_vector_field(vf_fa[0,:,k,:], vf_fa[2,:,k,:], step=1, maxlen=48, rescale=32), 'proj. do c. v. ponderado pela anisotropia fracional (plano transversal, fatia %s)' % k)
127     adshow(resize(div_evt_n[:,k,:],8), 'divergente do campo vetorial unitário (plano transversal, fatia %s)' % k)
128     adshow(resize(div_evl_n[:,k,:],8), 'divergente do c. v. ponderado pelos autovalores (plano transversal, fatia %s)' % k)
129     #adshow(resize(div_fa_n[:,:,k,:],8), 'divergente do c. v. ponderado pela anisotropia fracional (plano transversal, fatia %s)' % k)             #colorido
130     adshow(resize(div_fa_n[:,k,:],8), 'divergente do c. v. ponderado pela anisotropia fracional (plano transversal, fatia %s)' % k)                #grayscale
131     adshow(show_vector_field(rot_evt[0][:,k,:], rot_evt[2][:,k,:], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do campo vetorial unitário (plano transversal, fatia %s)' % k)
132     adshow(show_vector_field(rot_evl[0][:,k,:], rot_evl[2][:,k,:], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do c. v. ponderado pelos autovalores (plano transversal, fatia %s)' % k)
133     adshow(show_vector_field(rot_fa[0][:,k,:], rot_fa[2][:,k,:], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do c. v. ponderado pela anisotropia fracional (plano transversal, fatia %s)' % k)
134     adshow(resize(rot1_evt_n[:,k,:],8), 'comp. perpendicular do rot. do campo vetorial unitário (plano transversal, fatia %s)' % k)
135     adshow(resize(rot1_evl_n[:,k,:],8), 'comp. perpendicular do rot. do c. v. ponderado pelos autovalores (plano transversal, fatia %s)' % k)
136     adshow(resize(rot1_fa_n[:,k,:],8), 'comp. perpendicular do rot. do c. v. ponderado pela anisotropia fracional (plano transversal, fatia %s)' % k)
137
138 for k in np.arange(D):
139     adshow(show_vector_field(evt[1,k,:,:], evt[2,k,:,:], step=1, maxlen=24, rescale=16), 'projeção do campo vetorial unitário (plano sagital, fatia %s)' % k)
140     adshow(show_vector_field(vf_evl[1,k,:,:], vf_evl[2,k,:,:], step=1, maxlen=24, rescale=16), 'proj. do c. v. ponderado pelos autovalores (plano sagital, fatia %s)' % k)
141     adshow(show_vector_field(vf_fa[1,k,:,:], vf_fa[2,k,:,:], step=1, maxlen=24, rescale=16), 'proj. do c. v. ponderado pela anisotropia fracional (plano sagital, fatia %s)' % k)
142     adshow(resize(div_evt_n[k,:,:],8), 'divergente do campo vetorial unitário (plano sagital, fatia %s)' % k)
143     adshow(resize(div_evl_n[k,:,:],8), 'divergente do c. v. ponderado pelos autovalores (plano sagital, fatia %s)' % k)
144     #adshow(resize(div_fa_n[:,k,:,:],8), 'divergente do c. v. ponderado pela anisotropia fracional (plano sagital, fatia %s)' % k)           #colorido
145     adshow(resize(div_fa_n[k,:,:],8), 'divergente do c. v. ponderado pela anisotropia fracional (plano sagital, fatia %s)' % k)              #grayscale
146     adshow(show_vector_field(rot_evt[1][k,:,:], rot_evt[2][k,:,:], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do campo vetorial unitário (plano sagital, fatia %s)' % k)
147     adshow(show_vector_field(rot_evl[1][k,:,:], rot_evl[2][k,:,:], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do c. v. ponderado pelos autovalores (plano sagital, fatia %s)' % k)
148     adshow(show_vector_field(rot_fa[1][k,:,:], rot_fa[2][k,:,:], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do c. v. ponderado pela anisotropia fracional (plano sagital, fatia %s)' % k)
149     adshow(resize(rot0_evt_n[k,:,:],8), 'comp. perpendicular do rot. do campo vetorial unitário (plano sagital, fatia %s)' % k)
150     adshow(resize(rot0_evl_n[k,:,:],8), 'comp. perpendicular do rot. do c. v. ponderado pelos autovalores (plano sagital, fatia %s)' % k)
151     adshow(resize(rot0_fa_n[k,:,:],8), 'comp. perpendicular do rot. do c. v. ponderado pela anisotropia fracional (plano sagital, fatia %s)' % k)
152
153 for k in np.arange(W):
154     adshow(show_vector_field(evt[0,:,:,k], evt[1,:,:,k], step=1, maxlen=24, rescale=16), 'projeção do campo vetorial unitário (plano coronal, fatia %s)' % k)
155     adshow(show_vector_field(vf_evl[0,:,:,k], vf_evl[1,:,:,k], step=1, maxlen=24, rescale=16), 'proj. do c. v. ponderado pelos autovalores (plano coronal, fatia %s)' % k)
156     adshow(show_vector_field(vf_fa[0,:,:,k], vf_fa[1,:,:,k], step=1, maxlen=24, rescale=16), 'proj. do c. v. ponderado pela anisotropia fracional (plano coronal, fatia %s)' % k)
157     adshow(resize(div_evt_n[:,:,k],8), 'divergente do campo vetorial unitário (plano coronal, fatia %s)' % k)
158     adshow(resize(div_evl_n[:,:,k],8), 'divergente do c. v. ponderado pelos autovalores (plano coronal, fatia %s)' % k)
159     #adshow(resize(div_fa_n[:,:,:,k],8), 'divergente do c. v. ponderado pela anisotropia fracional (plano coronal, fatia %s)' % k)           #colorido
160     adshow(resize(div_fa_n[:,:,k],8), 'divergente do c. v. ponderado pela anisotropia fracional (plano coronal, fatia %s)' % k)            #grayscale
161     adshow(show_vector_field(rot_evt[0][:,:,k], rot_evt[1][:,:,k], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do campo vetorial unitário (plano coronal, fatia %s)' % k)
162     adshow(show_vector_field(rot_evl[0][:,:,k], rot_evl[1][:,:,k], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do c. v. ponderado pelos autovalores (plano coronal, fatia %s)' % k)
163     adshow(show_vector_field(rot_fa[0][:,:,k], rot_fa[1][:,:,k], step=1, maxlen=12, rescale=8, color=(255,0,0), bgcolor=(0,0,0)), 'proj. do rotacional do c. v. ponderado pela anisotropia fracional (plano coronal, fatia %s)' % k)
164     adshow(resize(rot2_evt_n[:,:,k],8), 'comp. perpendicular do rot. do campo vetorial unitário (plano coronal, fatia %s)' % k)
165     adshow(resize(rot2_evl_n[:,:,k],8), 'comp. perpendicular do rot. do c. v. ponderado pelos autovalores (plano coronal, fatia %s)' % k)
166     adshow(resize(rot2_fa_n[:,:,k],8), 'comp. perpendicular do rot. do c. v. ponderado pela anisotropia fracional (plano coronal, fatia %s)' % k)
```
```evt.shape = (3, 21, 17, 81)
evl.shape = (21, 17, 81)
fa.shape = (21, 17, 81)
np.max(evt[0,:,:,:]) = 0.999995708466
np.max(evt[1,:,:,:]) = 0.999978601933
np.max(evt[2,:,:,:]) = 0.99995881319
np.min(evt[0,:,:,:]) = -0.999998629093
np.min(evt[1,:,:,:]) = -0.999996364117
np.min(evt[2,:,:,:]) = 1.01282703326e-05
np.max(evl) = 0.00522658
np.min(evl) = 0.0
np.max(fa) = 1.0
np.min(fa) = 0.0
np.min(vf_evl) = -0.0044590559602
np.max(vf_evl) = 0.00517595795979
np.min(vf_fa) = -0.999494969845
np.max(vf_fa) = 0.996011972427
np.min(vf_fa_evl) = -0.00242391554166
np.max(vf_fa_evl) = 0.00241512827758
np.min(div_evt) = -3.98587960005
np.max(div_evt) = 3.98248390853
np.min(div_evl) = -0.0110773401067
np.max(div_evl) = 0.0144738547319
np.min(div_fa) = -2.63336205914
np.max(div_fa) = 2.39728021897
np.min(rot_evt) = -3.69270610809
np.max(rot_evt) = 3.69779676199
np.min(rot_evl) = -0.0111799504302
np.max(rot_evl) = 0.0132315950596
np.min(rot_fa) = -2.35710590064
np.max(rot_fa) = 2.32674971867
```

teste para o projeto

``` 1 import ia636 as ia
2 import numpy as np
3
5
6 #adshow(ia.ianormalize(divergente3((np.repeat(np.expand_dims(f,0), 5, 0),np.repeat(np.expand_dims(f,0), 5, 0),np.repeat(np.expand_dims(f,0), 5, 0)),2)))
11
12 #z, y, x = vector_field3(lambda z,y,x: sin(4*x)+sin(4*y)+sin(4*z), lambda z,y,x: sin(4*x)+sin(4*y)+sin(4*z), lambda z,y,x: sin(4*x)+sin(4*y)+sin(4*z))
13 #div = divergente3((z,y,x))
14 #rotz, roty, rotx = rotacional((z,y,x))
16 #adshow(ia.ianormalize(rotx[40,:,:]), 'Rotacional na direção x')
17 #adshow(ia.ianormalize(roty[40,:,:]), 'Rotacional na direção y')
18 #adshow(ia.ianormalize(rotz[40,:,:]), 'Rotacional na direção z')
19
21 #f = cv2.GaussianBlur(f, (0,0), 20, borderType=cv2.BORDER_REPLICATE)
22 #k2 = np.array([[1, 0, -1]])
23 #x = ia.iaconv(f,k2)[:,1:-1]
24 #y = ia.iaconv(f,k2.T)[1:-1,:]