1 #include "afni_warp.h"
2
3 #define NPER 10
4
5 static int nper = NPER ;
6
THD_orient_guess(MRI_IMAGE * imm)7 int_triple THD_orient_guess( MRI_IMAGE *imm )
8 {
9 int nvox , ii , nx,ny,nxy,nz , ix,jy,kz , icm,jcm,kcm , nbar ;
10 byte *bar , bp,bm ;
11 float xcm , ycm , zcm , ff , dx,dy,dz ;
12 float xx,yy,zz ;
13 int ni,nj,nk , itop,jtop,ktop , im,ip , jm,jp , km,kp ;
14 float axx,ayy,azz , clip , qx,qy,qz , arr[3] ;
15 int d_LR , d_AP , d_IS ;
16
17 int_triple it = {-1,-1,-1} ;
18
19 /*-- check for bad input --*/
20
21 if( imm == NULL || imm->nx < 5 || imm->ny < 5 || imm->nz < 5 ) return it ;
22
23 nx = imm->nx; ny = imm->ny; nz = imm->nz; nxy = nx*ny; nvox = nx*ny*nz;
24
25 dx = fabs(imm->dx) ; if( dx == 0.0 ) dx = 1.0 ;
26 dy = fabs(imm->dy) ; if( dy == 0.0 ) dy = 1.0 ;
27 dz = fabs(imm->dz) ; if( dz == 0.0 ) dz = 1.0 ;
28
29 /*-- make mask of NPER levels --*/
30
31 bar = (byte *) malloc( sizeof(byte) * nvox ) ;
32 clip = THD_cliplevel( imm , 0.5 ) ;
33
34 /* start with a binary mask */
35
36 switch( imm->kind ){
37 case MRI_float:{
38 float *ar = MRI_FLOAT_PTR(imm) ;
39 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (ar[ii] >= clip);
40 }
41 break ;
42 case MRI_short:{
43 short *ar = MRI_SHORT_PTR(imm) ;
44 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (ar[ii] >= clip);
45 }
46 break ;
47 case MRI_byte:{
48 byte *ar = MRI_BYTE_PTR(imm) ;
49 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (ar[ii] >= clip);
50 }
51 break ;
52 }
53
54 nbar = THD_countmask(nvox,bar) ;
55 printf("%d voxels in initial binary mask\n",nbar) ;
56 if( nbar == 0 ){ free(bar); return it; } /* bad */
57
58 THD_mask_clust( nx,ny,nz , bar ) ; /* take biggest cluster */
59
60 nbar = THD_countmask(nvox,bar) ;
61 printf("%d voxels in final binary mask\n",nbar) ;
62
63 #ifdef NPER
64 if( nper > 1 ){
65 float per[NPER+1] ; MRI_IMAGE *qim ; int jj ;
66 qim = mri_new( nbar , 1 , imm->kind ) ;
67 switch(imm->kind){
68 case MRI_float:{
69 float *ar=MRI_FLOAT_PTR(imm) , *qar=MRI_FLOAT_PTR(qim) ;
70 for( jj=ii=0 ; ii < nvox ; ii++ ) if( bar[ii] ) qar[jj++] = ar[ii] ;
71 }
72 break ;
73 case MRI_short:{
74 short *ar=MRI_SHORT_PTR(imm) , *qar=MRI_SHORT_PTR(qim) ;
75 for( jj=ii=0 ; ii < nvox ; ii++ ) if( bar[ii] ) qar[jj++] = ar[ii] ;
76 }
77 break ;
78 case MRI_byte:{
79 byte *ar=MRI_BYTE_PTR(imm) , *qar=MRI_BYTE_PTR(qim) ;
80 for( jj=ii=0 ; ii < nvox ; ii++ ) if( bar[ii] ) qar[jj++] = ar[ii] ;
81 }
82 break ;
83 }
84 printf("call mri_percents\n") ;
85 mri_percents( qim , nper , per ) ; /* compute nper break points */
86 mri_free(qim) ;
87 printf("per:") ;
88 for( ii=0 ; ii <= nper ; ii++ ) printf(" %g",per[ii]) ;
89 printf("\n") ;
90 switch( imm->kind ){
91 case MRI_float:{
92 float *ar = MRI_FLOAT_PTR(imm) , val ;
93 for( ii=0 ; ii < nvox ; ii++ ){
94 if( bar[ii] ){
95 val = ar[ii] ;
96 for( jj=1 ; jj <= nper && val >= per[jj] ; jj++ ) ; /*spin*/
97 bar[ii] = jj ;
98 }
99 }
100 }
101 break ;
102 case MRI_short:{
103 short *ar = MRI_SHORT_PTR(imm) , val ;
104 for( ii=0 ; ii < nvox ; ii++ ){
105 if( bar[ii] ){
106 val = ar[ii] ;
107 for( jj=1 ; jj <= nper && val >= per[jj] ; jj++ ) ; /*spin*/
108 bar[ii] = jj ;
109 }
110 }
111 }
112 break ;
113 case MRI_byte:{
114 byte *ar = MRI_BYTE_PTR(imm) , val ;
115 for( ii=0 ; ii < nvox ; ii++ ){
116 if( bar[ii] ){
117 val = ar[ii] ;
118 for( jj=1 ; jj <= nper && val >= per[jj] ; jj++ ) ; /*spin*/
119 bar[ii] = jj ;
120 }
121 }
122 }
123 break ;
124 }
125 }
126 #endif /* NPER */
127
128 /* find center of mass of mask */
129
130 xcm = ycm = zcm = ff = 0.0 ;
131 for( ii=0 ; ii < nvox ; ii++ ){
132 if( bar[ii] ){
133 ix = (ii % nx) ; xx = ix*dx ; xcm += xx*bar[ii] ;
134 jy = (ii / nx) % ny ; yy = jy*dy ; ycm += yy*bar[ii] ;
135 kz = (ii /nxy) ; zz = kz*dz ; zcm += zz*bar[ii] ;
136 ff += bar[ii] ;
137 }
138 }
139 xcm /= ff ; ycm /= ff ; zcm /= ff ;
140
141 icm = rint(xcm/dx) ;
142 itop = 2*icm ; if( itop >= nx ) itop = nx-1 ;
143 ni = itop-icm ;
144
145 jcm = rint(ycm/dy) ;
146 jtop = 2*jcm ; if( jtop >= ny ) jtop = ny-1 ;
147 nj = jtop-jcm ;
148
149 kcm = rint(zcm/dz) ;
150 ktop = 2*kcm ; if( ktop >= nz ) ktop = nz-1 ;
151 nk = ktop-kcm ;
152
153 printf("Mask count = %d\n"
154 "icm = %d jcm = %d kcm = %d\n"
155 "ni = %d nj = %d nk = %d\n",
156 (int)ff , icm,jcm,kcm , ni,nj,nk ) ;
157
158 /** compute asymmetry measures about CM voxel **/
159
160 #define BAR(i,j,k) bar[(i)+(j)*nx+(k)*nxy]
161
162 axx = 0.0 ;
163 for( ix=1 ; ix <= ni ; ix++ ){
164 im = icm-ix ; ip = icm+ix ;
165 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
166 for( jy=jcm-nj ; jy <= jcm+nj ; jy++ )
167 axx += abs(BAR(ip,jy,kz) - BAR(im,jy,kz)) ;
168 }
169 }
170 axx /= (ni*nj*nk) ; printf("axx = %g\n",axx) ;
171
172 ayy = 0.0 ;
173 for( jy=1 ; jy <= nj ; jy++ ){
174 jm = jcm-jy ; jp = jcm+jy ;
175 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
176 for( ix=icm-ni ; ix <= icm+ni ; ix++ )
177 ayy += abs(BAR(ix,jp,kz) - BAR(ix,jm,kz)) ;
178 }
179 }
180 ayy /= (ni*nj*nk) ; printf("ayy = %g\n",ayy) ;
181
182 azz = 0.0 ;
183 for( kz=1 ; kz <= nk ; kz++ ){
184 km = kcm-kz ; kp = kcm+kz ;
185 for( jy=jcm-nj ; jy <= jcm+nj ; jy++ ){
186 for( ix=icm-ni ; ix <= icm+ni ; ix++ )
187 azz += abs(BAR(ix,jy,kp) - BAR(ix,jy,km)) ;
188 }
189 }
190 azz /= (ni*nj*nk) ; printf("azz = %g\n",azz) ;
191
192 /** least asymettric is L-R direction **/
193
194 if( axx < ayy ){
195 if( axx < azz ) d_LR = 1 ;
196 else d_LR = 3 ;
197 } else {
198 if( ayy < azz ) d_LR = 2 ;
199 else d_LR = 3 ;
200 }
201 printf("axis %d is L-R\n",d_LR) ;
202
203 arr[0] = axx ; arr[1] = ayy ; arr[2] = azz ; ff = arr[d_LR-1] ;
204 arr[0] /= ff ;
205 arr[1] /= ff ;
206 arr[2] /= ff ;
207 printf("a ratios = %g %g %g\n",arr[0],arr[1],arr[2]) ;
208
209 /** find asymmetry measures in 1/2 spaces perp to L-R **/
210
211 switch( d_LR ){
212
213 case 3:{ /* L-R is z-axis */
214 float axx_jp=0.0, axx_jm=0.0, ayy_ip=0.0, ayy_im=0.0 ;
215 for( ix=1 ; ix <= ni ; ix++ ){
216 im = icm-ix ; ip = icm+ix ;
217 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
218 for( jy=1 ; jy <= nj ; jy++ ){
219 axx_jp += abs(BAR(ip,jcm+jy,kz) - BAR(im,jcm+jy,kz)) ;
220 axx_jm += abs(BAR(ip,jcm-jy,kz) - BAR(im,jcm-jy,kz)) ;
221 }
222 }
223 }
224 for( jy=1 ; jy <= nj ; jy++ ){
225 jm = jcm-jy ; jp = jcm+jy ;
226 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
227 for( ix=1 ; ix <= ni ; ix++ ){
228 ayy_ip += abs(BAR(icm+ix,jp,kz) - BAR(icm+ix,jm,kz)) ;
229 ayy_im += abs(BAR(icm-ix,jp,kz) - BAR(icm-ix,jm,kz)) ;
230 }
231 }
232 }
233 axx_jp /= (ni*nj*nk) ; axx_jm /= (ni*nj*nk) ;
234 ayy_ip /= (ni*nj*nk) ; ayy_im /= (ni*nj*nk) ;
235
236 printf("axx_jp=%g axx_jm=%g ayy_ip=%g ayy_im=%g\n",
237 axx_jp,axx_jm , ayy_ip,ayy_im ) ;
238 } /* end of case 3 */
239 break ;
240
241 }
242
243 return it ;
244 }
245
main(int argc,char * argv[])246 int main( int argc , char *argv[] )
247 {
248 THD_3dim_dataset *dset ;
249 MRI_IMAGE *im ;
250 int iarg=1 ;
251
252 if( strcmp(argv[1],"-nper") == 0 ){
253 nper = strtol( argv[2] , NULL , 10 ) ;
254 iarg = 3 ;
255 }
256
257 for( ; iarg < argc ; iarg++ ){
258 printf("======= dataset %s nper=%d ========\n",argv[iarg],nper) ;
259 dset = THD_open_dataset(argv[iarg]) ;
260 if( dset == NULL ) continue ;
261 DSET_load(dset) ;
262 im = DSET_BRICK(dset,0) ;
263 im->dx = DSET_DX(dset) ;
264 im->dy = DSET_DY(dset) ;
265 im->dz = DSET_DZ(dset) ;
266 (void) THD_orient_guess( im ) ;
267
268 printf( "Data Axes Orientation:\n"
269 " first (x) = %s\n"
270 " second (y) = %s\n"
271 " third (z) = %s\n" ,
272 ORIENT_typestr[dset->daxes->xxorient] ,
273 ORIENT_typestr[dset->daxes->yyorient] ,
274 ORIENT_typestr[dset->daxes->zzorient] ) ;
275 DSET_delete(dset) ;
276 }
277
278 exit(0) ;
279 }
280