Doxygen Source Code Documentation
thd_orient.c File Reference
#include "afni_warp.h"Go to the source code of this file.
Defines | |
| #define | NPER 10 |
| #define | BAR(i, j, k) bar[(i)+(j)*nx+(k)*nxy] |
Functions | |
| int_triple | THD_orient_guess (MRI_IMAGE *imm) |
| int | main (int argc, char *argv[]) |
Variables | |
| int | nper = NPER |
Define Documentation
|
|
|
|
|
Definition at line 3 of file thd_orient.c. Referenced by THD_orient_guess(). |
Function Documentation
|
||||||||||||
|
find asymmetry measures in 1/2 spaces perp to L-R * Definition at line 246 of file thd_orient.c. References argc, THD_3dim_dataset::daxes, DSET_BRICK, DSET_delete, DSET_DX, DSET_DY, DSET_DZ, DSET_load, MRI_IMAGE::dx, MRI_IMAGE::dy, MRI_IMAGE::dz, nper, THD_open_dataset(), THD_orient_guess(), THD_dataxes::xxorient, THD_dataxes::yyorient, and THD_dataxes::zzorient.
00247 {
00248 THD_3dim_dataset *dset ;
00249 MRI_IMAGE *im ;
00250 int iarg=1 ;
00251
00252 if( strcmp(argv[1],"-nper") == 0 ){
00253 nper = strtol( argv[2] , NULL , 10 ) ;
00254 iarg = 3 ;
00255 }
00256
00257 for( ; iarg < argc ; iarg++ ){
00258 printf("======= dataset %s nper=%d ========\n",argv[iarg],nper) ;
00259 dset = THD_open_dataset(argv[iarg]) ;
00260 if( dset == NULL ) continue ;
00261 DSET_load(dset) ;
00262 im = DSET_BRICK(dset,0) ;
00263 im->dx = DSET_DX(dset) ;
00264 im->dy = DSET_DY(dset) ;
00265 im->dz = DSET_DZ(dset) ;
00266 (void) THD_orient_guess( im ) ;
00267
00268 printf( "Data Axes Orientation:\n"
00269 " first (x) = %s\n"
00270 " second (y) = %s\n"
00271 " third (z) = %s\n" ,
00272 ORIENT_typestr[dset->daxes->xxorient] ,
00273 ORIENT_typestr[dset->daxes->yyorient] ,
00274 ORIENT_typestr[dset->daxes->zzorient] ) ;
00275 DSET_delete(dset) ;
00276 }
00277
00278 exit(0) ;
00279 }
|
|
|
Definition at line 7 of file thd_orient.c. References abs, MRI_IMAGE::dx, MRI_IMAGE::dy, MRI_IMAGE::dz, free, MRI_IMAGE::kind, malloc, MRI_BYTE_PTR, MRI_FLOAT_PTR, mri_free(), mri_new(), mri_percents(), MRI_SHORT_PTR, NPER, nper, MRI_IMAGE::nx, MRI_IMAGE::ny, MRI_IMAGE::nz, nz, THD_cliplevel(), THD_countmask(), and THD_mask_clust(). Referenced by main().
00008 {
00009 int nvox , ii , nx,ny,nxy,nz , ix,jy,kz , icm,jcm,kcm , nbar ;
00010 byte *bar , bp,bm ;
00011 float xcm , ycm , zcm , ff , dx,dy,dz ;
00012 float xx,yy,zz ;
00013 int ni,nj,nk , itop,jtop,ktop , im,ip , jm,jp , km,kp ;
00014 float axx,ayy,azz , clip , qx,qy,qz , arr[3] ;
00015 int d_LR , d_AP , d_IS ;
00016
00017 int_triple it = {-1,-1,-1} ;
00018
00019 /*-- check for bad input --*/
00020
00021 if( imm == NULL || imm->nx < 5 || imm->ny < 5 || imm->nz < 5 ) return it ;
00022
00023 nx = imm->nx; ny = imm->ny; nz = imm->nz; nxy = nx*ny; nvox = nx*ny*nz;
00024
00025 dx = fabs(imm->dx) ; if( dx == 0.0 ) dx = 1.0 ;
00026 dy = fabs(imm->dy) ; if( dy == 0.0 ) dy = 1.0 ;
00027 dz = fabs(imm->dz) ; if( dz == 0.0 ) dz = 1.0 ;
00028
00029 /*-- make mask of NPER levels --*/
00030
00031 bar = (byte *) malloc( sizeof(byte) * nvox ) ;
00032 clip = THD_cliplevel( imm , 0.5 ) ;
00033
00034 /* start with a binary mask */
00035
00036 switch( imm->kind ){
00037 case MRI_float:{
00038 float *ar = MRI_FLOAT_PTR(imm) ;
00039 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (ar[ii] >= clip);
00040 }
00041 break ;
00042 case MRI_short:{
00043 short *ar = MRI_SHORT_PTR(imm) ;
00044 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (ar[ii] >= clip);
00045 }
00046 break ;
00047 case MRI_byte:{
00048 byte *ar = MRI_BYTE_PTR(imm) ;
00049 for( ii=0 ; ii < nvox ; ii++ ) bar[ii] = (ar[ii] >= clip);
00050 }
00051 break ;
00052 }
00053
00054 nbar = THD_countmask(nvox,bar) ;
00055 printf("%d voxels in initial binary mask\n",nbar) ;
00056 if( nbar == 0 ){ free(bar); return it; } /* bad */
00057
00058 THD_mask_clust( nx,ny,nz , bar ) ; /* take biggest cluster */
00059
00060 nbar = THD_countmask(nvox,bar) ;
00061 printf("%d voxels in final binary mask\n",nbar) ;
00062
00063 #ifdef NPER
00064 if( nper > 1 ){
00065 float per[NPER+1] ; MRI_IMAGE *qim ; int jj ;
00066 qim = mri_new( nbar , 1 , imm->kind ) ;
00067 switch(imm->kind){
00068 case MRI_float:{
00069 float *ar=MRI_FLOAT_PTR(imm) , *qar=MRI_FLOAT_PTR(qim) ;
00070 for( jj=ii=0 ; ii < nvox ; ii++ ) if( bar[ii] ) qar[jj++] = ar[ii] ;
00071 }
00072 break ;
00073 case MRI_short:{
00074 short *ar=MRI_SHORT_PTR(imm) , *qar=MRI_SHORT_PTR(qim) ;
00075 for( jj=ii=0 ; ii < nvox ; ii++ ) if( bar[ii] ) qar[jj++] = ar[ii] ;
00076 }
00077 break ;
00078 case MRI_byte:{
00079 byte *ar=MRI_BYTE_PTR(imm) , *qar=MRI_BYTE_PTR(qim) ;
00080 for( jj=ii=0 ; ii < nvox ; ii++ ) if( bar[ii] ) qar[jj++] = ar[ii] ;
00081 }
00082 break ;
00083 }
00084 printf("call mri_percents\n") ;
00085 mri_percents( qim , nper , per ) ; /* compute nper break points */
00086 mri_free(qim) ;
00087 printf("per:") ;
00088 for( ii=0 ; ii <= nper ; ii++ ) printf(" %g",per[ii]) ;
00089 printf("\n") ;
00090 switch( imm->kind ){
00091 case MRI_float:{
00092 float *ar = MRI_FLOAT_PTR(imm) , val ;
00093 for( ii=0 ; ii < nvox ; ii++ ){
00094 if( bar[ii] ){
00095 val = ar[ii] ;
00096 for( jj=1 ; jj <= nper && val >= per[jj] ; jj++ ) ; /*spin*/
00097 bar[ii] = jj ;
00098 }
00099 }
00100 }
00101 break ;
00102 case MRI_short:{
00103 short *ar = MRI_SHORT_PTR(imm) , val ;
00104 for( ii=0 ; ii < nvox ; ii++ ){
00105 if( bar[ii] ){
00106 val = ar[ii] ;
00107 for( jj=1 ; jj <= nper && val >= per[jj] ; jj++ ) ; /*spin*/
00108 bar[ii] = jj ;
00109 }
00110 }
00111 }
00112 break ;
00113 case MRI_byte:{
00114 byte *ar = MRI_BYTE_PTR(imm) , val ;
00115 for( ii=0 ; ii < nvox ; ii++ ){
00116 if( bar[ii] ){
00117 val = ar[ii] ;
00118 for( jj=1 ; jj <= nper && val >= per[jj] ; jj++ ) ; /*spin*/
00119 bar[ii] = jj ;
00120 }
00121 }
00122 }
00123 break ;
00124 }
00125 }
00126 #endif /* NPER */
00127
00128 /* find center of mass of mask */
00129
00130 xcm = ycm = zcm = ff = 0.0 ;
00131 for( ii=0 ; ii < nvox ; ii++ ){
00132 if( bar[ii] ){
00133 ix = (ii % nx) ; xx = ix*dx ; xcm += xx*bar[ii] ;
00134 jy = (ii / nx) % ny ; yy = jy*dy ; ycm += yy*bar[ii] ;
00135 kz = (ii /nxy) ; zz = kz*dz ; zcm += zz*bar[ii] ;
00136 ff += bar[ii] ;
00137 }
00138 }
00139 xcm /= ff ; ycm /= ff ; zcm /= ff ;
00140
00141 icm = rint(xcm/dx) ;
00142 itop = 2*icm ; if( itop >= nx ) itop = nx-1 ;
00143 ni = itop-icm ;
00144
00145 jcm = rint(ycm/dy) ;
00146 jtop = 2*jcm ; if( jtop >= ny ) jtop = ny-1 ;
00147 nj = jtop-jcm ;
00148
00149 kcm = rint(zcm/dz) ;
00150 ktop = 2*kcm ; if( ktop >= nz ) ktop = nz-1 ;
00151 nk = ktop-kcm ;
00152
00153 printf("Mask count = %d\n"
00154 "icm = %d jcm = %d kcm = %d\n"
00155 "ni = %d nj = %d nk = %d\n",
00156 (int)ff , icm,jcm,kcm , ni,nj,nk ) ;
00157
00158 /** compute asymmetry measures about CM voxel **/
00159
00160 #define BAR(i,j,k) bar[(i)+(j)*nx+(k)*nxy]
00161
00162 axx = 0.0 ;
00163 for( ix=1 ; ix <= ni ; ix++ ){
00164 im = icm-ix ; ip = icm+ix ;
00165 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
00166 for( jy=jcm-nj ; jy <= jcm+nj ; jy++ )
00167 axx += abs(BAR(ip,jy,kz) - BAR(im,jy,kz)) ;
00168 }
00169 }
00170 axx /= (ni*nj*nk) ; printf("axx = %g\n",axx) ;
00171
00172 ayy = 0.0 ;
00173 for( jy=1 ; jy <= nj ; jy++ ){
00174 jm = jcm-jy ; jp = jcm+jy ;
00175 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
00176 for( ix=icm-ni ; ix <= icm+ni ; ix++ )
00177 ayy += abs(BAR(ix,jp,kz) - BAR(ix,jm,kz)) ;
00178 }
00179 }
00180 ayy /= (ni*nj*nk) ; printf("ayy = %g\n",ayy) ;
00181
00182 azz = 0.0 ;
00183 for( kz=1 ; kz <= nk ; kz++ ){
00184 km = kcm-kz ; kp = kcm+kz ;
00185 for( jy=jcm-nj ; jy <= jcm+nj ; jy++ ){
00186 for( ix=icm-ni ; ix <= icm+ni ; ix++ )
00187 azz += abs(BAR(ix,jy,kp) - BAR(ix,jy,km)) ;
00188 }
00189 }
00190 azz /= (ni*nj*nk) ; printf("azz = %g\n",azz) ;
00191
00192 /** least asymettric is L-R direction **/
00193
00194 if( axx < ayy ){
00195 if( axx < azz ) d_LR = 1 ;
00196 else d_LR = 3 ;
00197 } else {
00198 if( ayy < azz ) d_LR = 2 ;
00199 else d_LR = 3 ;
00200 }
00201 printf("axis %d is L-R\n",d_LR) ;
00202
00203 arr[0] = axx ; arr[1] = ayy ; arr[2] = azz ; ff = arr[d_LR-1] ;
00204 arr[0] /= ff ;
00205 arr[1] /= ff ;
00206 arr[2] /= ff ;
00207 printf("a ratios = %g %g %g\n",arr[0],arr[1],arr[2]) ;
00208
00209 /** find asymmetry measures in 1/2 spaces perp to L-R **/
00210
00211 switch( d_LR ){
00212
00213 case 3:{ /* L-R is z-axis */
00214 float axx_jp=0.0, axx_jm=0.0, ayy_ip=0.0, ayy_im=0.0 ;
00215 for( ix=1 ; ix <= ni ; ix++ ){
00216 im = icm-ix ; ip = icm+ix ;
00217 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
00218 for( jy=1 ; jy <= nj ; jy++ ){
00219 axx_jp += abs(BAR(ip,jcm+jy,kz) - BAR(im,jcm+jy,kz)) ;
00220 axx_jm += abs(BAR(ip,jcm-jy,kz) - BAR(im,jcm-jy,kz)) ;
00221 }
00222 }
00223 }
00224 for( jy=1 ; jy <= nj ; jy++ ){
00225 jm = jcm-jy ; jp = jcm+jy ;
00226 for( kz=kcm-nk ; kz <= kcm+nk ; kz++ ){
00227 for( ix=1 ; ix <= ni ; ix++ ){
00228 ayy_ip += abs(BAR(icm+ix,jp,kz) - BAR(icm+ix,jm,kz)) ;
00229 ayy_im += abs(BAR(icm-ix,jp,kz) - BAR(icm-ix,jm,kz)) ;
00230 }
00231 }
00232 }
00233 axx_jp /= (ni*nj*nk) ; axx_jm /= (ni*nj*nk) ;
00234 ayy_ip /= (ni*nj*nk) ; ayy_im /= (ni*nj*nk) ;
00235
00236 printf("axx_jp=%g axx_jm=%g ayy_ip=%g ayy_im=%g\n",
00237 axx_jp,axx_jm , ayy_ip,ayy_im ) ;
00238 } /* end of case 3 */
00239 break ;
00240
00241 }
00242
00243 return it ;
00244 }
|
Variable Documentation
|
|
Definition at line 5 of file thd_orient.c. Referenced by main(), and THD_orient_guess(). |
