Doxygen Source Code Documentation
Intracranial.c File Reference
Go to the source code of this file.
Defines | |
| #define | AR(i, j) ar[(i)+(j)*nx] |
| #define | MAXLAT 90 |
| #define | MAXLNG 180 |
| #define | AIZE |
Functions | |
| int | verify_inputs () |
| void | draw_2dfiller (int nx, int ny, int ix, int jy, byte *ar) |
| void | center_of_mass (int *cx, int *cy, int *cz) |
| void | segment_x_slices (short *cv, int cx, Boolean axial_slice) |
| void | segment_y_slices (short *cv, int cy, Boolean axial_slice) |
| void | segment_z_slices (short *cv, int cz, Boolean axial_slice) |
| void | segment_envelope (short *cv, int cx, int cy, int cz) |
| void | segment_volume (short *cv) |
| int | connectivity_tests (short *cv) |
Define Documentation
|
|
|
|
|
|
|
|
|
|
|
|
Function Documentation
|
||||||||||||||||
|
Definition at line 163 of file Intracranial.c. References DSET_BRICK_ARRAY, DSET_NX, DSET_NY, DSET_NZ, and nz. Referenced by segment_volume().
00165 {
00166 int nx, ny, nz, nxy, nxyz; /* dataset dimensions in voxels */
00167 int ix, jy, kz, ixyz; /* voxel indices */
00168 short * anat_data = NULL; /* data from anatomical image */
00169 float f, sum, sumx, sumy, sumz; /* sums for calculating means */
00170
00171
00172 /*----- Progress report -----*/
00173 if (! quiet) printf ("\nCalculating center of mass \n");
00174
00175
00176 /*----- Initialize local variables -----*/
00177 anat_data = (short *) DSET_BRICK_ARRAY(anat,0) ;
00178 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00179 nxy = nx * ny; nxyz = nxy * nz;
00180
00181
00182 /*----- Sum over all voxels -----*/
00183 sum = 0.0; sumx = 0.0; sumy = 0.0; sumz = 0.0;
00184 for (kz = 0; kz < nz; kz++)
00185 for (jy = 0; jy < ny; jy++)
00186 for (ix = 0; ix < nx; ix++)
00187 {
00188 ixyz = ix + jy*nx + kz*nxy;
00189 f = anat_data[ixyz];
00190 sum += f; sumx += ix*f; sumy += jy*f; sumz += kz*f;
00191 }
00192
00193
00194 /*----- Calculate center of mass -----*/
00195 *cx = sumx / sum;
00196 *cy = sumy / sum;
00197 *cz = sumz / sum;
00198
00199 if (! quiet) printf ("Center of mass: ix = %d jy = %d kz = %d \n",
00200 *cx, *cy, *cz);
00201
00202 }
|
|
|
Definition at line 988 of file Intracranial.c. References DSET_NX, DSET_NY, DSET_NZ, free, IJK_TO_THREE, malloc, MTEST, nz, and THREE_TO_IJK. Referenced by SEGMENT_auto().
00990 {
00991 int nx, ny, nz, nxy, nxyz; /* voxel counters */
00992 int ix, jy, kz, ixyz; /* voxel indices */
00993 byte * dv = NULL; /* volume for intermediate data */
00994
00995
00996 /*----- Progress report -----*/
00997 if (! quiet) printf ("\nPerforming voxel connectivity tests \n");
00998
00999
01000 /*----- Initialize local variables -----*/
01001 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
01002 nxy = nx*ny; nxyz = nxy*nz;
01003
01004
01005 /*----- Initialize voxel classification indicators -----*/
01006 dv = (byte *) malloc (nxyz * sizeof(byte));
01007 MTEST (dv); if (dv == NULL) return (0);
01008 for (ixyz = 0; ixyz < nxyz; ixyz++)
01009 dv[ixyz] = 0;
01010
01011
01012 /*----- Determine which voxels will leave the target structure -----*/
01013 if (max_conn_int > -1)
01014 {
01015 for (ixyz = 0; ixyz < nxyz; ixyz++)
01016 {
01017 if (! cv[ixyz])
01018 {
01019 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
01020 if (ix > 1) dv[THREE_TO_IJK(ix-1,jy,kz,nx,nxy)]++;
01021 if (ix < nx-1) dv[THREE_TO_IJK(ix+1,jy,kz,nx,nxy)]++;
01022 if (jy > 1) dv[THREE_TO_IJK(ix,jy-1,kz,nx,nxy)]++;
01023 if (jy < ny-1) dv[THREE_TO_IJK(ix,jy+1,kz,nx,nxy)]++;
01024 if (kz > 1) dv[THREE_TO_IJK(ix,jy,kz-1,nx,nxy)]++;
01025 if (kz < nz-1) dv[THREE_TO_IJK(ix,jy,kz+1,nx,nxy)]++;
01026 }
01027 }
01028
01029 for (ixyz = 0; ixyz < nxyz; ixyz++)
01030 {
01031 if (dv[ixyz] <= max_conn_int) cv[ixyz] = 1;
01032 dv[ixyz] = 0;
01033 }
01034 }
01035
01036
01037 /*----- Determine which voxels will enter the target structure -----*/
01038 if (min_conn_int < 7)
01039 {
01040 for (ixyz = 0; ixyz < nxyz; ixyz++)
01041 {
01042 if (! cv[ixyz])
01043 {
01044 IJK_TO_THREE (ixyz, ix, jy, kz, nx, nxy);
01045 if (ix > 1) dv[THREE_TO_IJK(ix-1,jy,kz,nx,nxy)]++;
01046 if (ix < nx-1) dv[THREE_TO_IJK(ix+1,jy,kz,nx,nxy)]++;
01047 if (jy > 1) dv[THREE_TO_IJK(ix,jy-1,kz,nx,nxy)]++;
01048 if (jy < ny-1) dv[THREE_TO_IJK(ix,jy+1,kz,nx,nxy)]++;
01049 if (kz > 1) dv[THREE_TO_IJK(ix,jy,kz-1,nx,nxy)]++;
01050 if (kz < nz-1) dv[THREE_TO_IJK(ix,jy,kz+1,nx,nxy)]++;
01051 }
01052 }
01053
01054 for (ixyz = 0; ixyz < nxyz; ixyz++)
01055 if (dv[ixyz] >= min_conn_int) cv[ixyz] = 0;
01056 }
01057
01058
01059 /*----- Deallocate memory -----*/
01060 free (dv); dv = NULL;
01061
01062
01063 return (1);
01064 }
|
|
||||||||||||||||||||||||
|
Definition at line 118 of file Intracranial.c. Referenced by segment_x_slices(), segment_y_slices(), and segment_z_slices().
00119 {
00120 int ii,jj , ip,jp , num ;
00121
00122 #define AR(i,j) ar[(i)+(j)*nx]
00123
00124
00125 /*----- Test for early termination -----*/
00126 if (AR(ix,jy) != 0) return;
00127
00128
00129 /* fill out in cross from 1st point */
00130
00131 ip = ix ; jp = jy ; AR(ip,jp) = 2 ;
00132
00133 for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ) AR(ii,jp) = 2;
00134 for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ) AR(ii,jp) = 2;
00135 for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ) AR(ip,jj) = 2;
00136 for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ) AR(ip,jj) = 2;
00137
00138 /* brute force repetition of the cross technique */
00139
00140 do {
00141 num = 0 ;
00142 for( jp=0 ; jp < ny ; jp++ ){
00143 for( ip=0 ; ip < nx ; ip++ ){
00144 if( AR(ip,jp) == 2 ){
00145 for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ){ AR(ii,jp) = 2; num++; }
00146 for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ){ AR(ii,jp) = 2; num++; }
00147 for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ){ AR(ip,jj) = 2; num++; }
00148 for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ){ AR(ip,jj) = 2; num++; }
00149 }
00150 }
00151 }
00152 } while( num > 0 ) ;
00153
00154 return ;
00155 }
|
|
||||||||||||||||||||
|
Definition at line 760 of file Intracranial.c. References DSET_NX, DSET_NY, DSET_NZ, free, i, malloc, nr, nz, and r. Referenced by segment_volume().
00765 {
00766 #define MAXLAT 90
00767 #define MAXLNG 180
00768
00769 int nx, ny, nz, nxy, nxz, nyz, nxyz; /* dataset dimensions in voxels */
00770 int ix, jy, kz, ixy, ixz, iyz, ixyz; /* voxel indices */
00771
00772 /* 24 Mar 2004: on Mac OS X, the compiler generates incorrect code
00773 for these large arrays when optimization is turned on;
00774 instead, use malloc() to get them if AIZE is not defined */
00775
00776 #undef AIZE
00777 #ifndef DARWIN
00778 # define AIZE
00779 #endif
00780
00781 #ifdef AIZE
00782 float radius[MAXLAT][MAXLNG]; /* max. radius vector to brain voxel */
00783 float smradius[MAXLAT][MAXLNG]; /* array of smoothed radius vectors */
00784 #else
00785 float **radius , **smradius ;
00786 #endif
00787
00788 int ilat, jlat; /* radius array latitude index */
00789 int ilng, jlng; /* radius array longitude index */
00790 float deltalat, deltalng; /* increments in latitude and longitude */
00791 float lat, clat, slat; /* voxel latitude */
00792 float lng, clng, slng; /* voxel longitude */
00793 int ir, i, j; /* indices */
00794 float rxy; /* radius to voxel in xy-plane */
00795 float r; /* radius to voxel from center of brain */
00796 int nr; /* maximum test radius */
00797 float maxr; /* maximum radius to voxel inside the brain */
00798
00799
00800 /*----- Progress report -----*/
00801 if (! quiet) printf ("\nCreating envelope for segmented image \n");
00802
00803 #ifndef AIZE
00804 radius = (float **) malloc(sizeof(float)*MAXLAT) ;
00805 smradius = (float **) malloc(sizeof(float)*MAXLAT) ;
00806 for (ilat = 0; ilat < MAXLAT; ilat++){
00807 radius[ilat] = (float *)malloc(sizeof(float)*MAXLNG) ;
00808 smradius[ilat] = (float *)malloc(sizeof(float)*MAXLNG) ;
00809 }
00810 #endif
00811
00812
00813 /*----- Initialize local variables -----*/
00814 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00815 nxy = nx*ny; nxz = nx*nz; nyz = ny*nz; nxyz = nxy*nz;
00816 deltalat = PI / MAXLAT;
00817 deltalng = 2.0 * PI / MAXLNG;
00818 nr = nx;
00819 if (ny > nr) nr = ny;
00820 if (nz > nr) nr = nz;
00821
00822
00823
00824 /*----- Calculate the radius vector -----*/
00825 for (ilat = 0; ilat < MAXLAT; ilat++)
00826 {
00827 lat = ilat * deltalat - PI/2;
00828 slat = sin(lat);
00829 clat = cos(lat);
00830
00831 for (ilng = 0; ilng < MAXLNG; ilng++)
00832 {
00833 lng = ilng * deltalng;
00834 clng = cos(lng);
00835 slng = sin(lng);
00836
00837 radius[ilat][ilng] = 0.0;
00838 for (ir = 0; ir < nr; ir++)
00839 {
00840 ix = cx + ir*clat*clng;
00841 jy = cy + ir*clat*slng;
00842 kz = cz + ir * slat;
00843
00844 if ((ix >= 0) && (ix < nx) && (jy >= 0) && (jy < ny)
00845 && (kz >= 0) && (kz < nz))
00846 {
00847 ixyz = ix + jy*nx + kz*nxy;
00848 if (! cv[ixyz]) radius[ilat][ilng] = (float) ir;
00849 }
00850 }
00851 }
00852 }
00853
00854 /*----- Smooth the radius vectors -----*/
00855 for (ilat = 0; ilat < MAXLAT; ilat++)
00856 {
00857 for (ilng = 0; ilng < MAXLNG; ilng++)
00858 {
00859 smradius[ilat][ilng] = 0.0;
00860
00861 for (i = -1; i <= 1; i++)
00862 {
00863 jlat = (ilat + i + MAXLAT) % MAXLAT;
00864
00865 if (ilat == 0) jlat = 0;
00866 if (ilat == MAXLAT-1) jlat = MAXLAT-1;
00867
00868 for (j = -2; j <= 2; j++)
00869 {
00870 jlng = (ilng + j + MAXLNG) % MAXLNG;
00871 smradius[ilat][ilng] += radius[jlat][jlng] / 15.0;
00872 }
00873 }
00874 }
00875 }
00876
00877 /*----- Find maximum radius to a brain voxel -----*/
00878 maxr = 0.0;
00879 for (ilat = 0; ilat < MAXLAT; ilat++)
00880 for (ilng = 0; ilng < MAXLNG; ilng++)
00881 if (maxr < smradius[ilat][ilng]) maxr = smradius[ilat][ilng];
00882
00883 /*----- Smooth the brain mask -----*/
00884 for (ix = 0; ix < nx; ix++){
00885 for (jy = 0; jy < ny; jy++)
00886 {
00887 rxy = hypot ((float) (ix-cx), (float) (jy-cy));
00888 lng = atan2 (jy-cy, ix-cx);
00889 ilng = (int) ((lng/deltalng) + MAXLNG) % MAXLNG;
00890
00891 for (kz = 0; kz < nz; kz++)
00892 {
00893 ixyz = ix + jy*nx + kz*nxy;
00894
00895 r = hypot (rxy, (float) (kz-cz));
00896
00897 if (r < maxr)
00898 {
00899 lat = atan2 (kz-cz, rxy);
00900 ilat = (int) (((lat+PI/2)/deltalat) + MAXLAT) % MAXLAT;
00901
00902 if (r < smradius[ilat][ilng]) cv[ixyz] = 0;
00903 }
00904
00905 }
00906 }
00907 }
00908
00909 #ifndef AIZE
00910 for (ilat = 0; ilat < MAXLAT; ilat++){
00911 free( radius[ilat]) ; free(smradius[ilat]) ;
00912 }
00913 free(smradius); free(radius);
00914 #endif
00915 }
|
|
|
Definition at line 924 of file Intracranial.c. References center_of_mass(), DSET_NX, DSET_NY, DSET_NZ, segment_envelope(), segment_x_slices(), segment_y_slices(), segment_z_slices(), SI_error(), THD_dataxes::xxorient, THD_dataxes::yyorient, and THD_dataxes::zzorient. Referenced by SEGMENT_auto().
00928 {
00929 THD_dataxes * daxes; /* dataset axes */
00930 char xxorient, yyorient, zzorient; /* dataset axes orientations */
00931 int nxyz; /* dataset dimension in voxels */
00932 int ixyz; /* voxel indices */
00933 int cx, cy, cz; /* location of center of mass */
00934
00935
00936 /*----- Initialize local variables -----*/
00937 daxes = anat->daxes;
00938 xxorient = ORIENT_typestr[daxes->xxorient][0];
00939 yyorient = ORIENT_typestr[daxes->yyorient][0];
00940 zzorient = ORIENT_typestr[daxes->zzorient][0];
00941 nxyz = DSET_NX(anat) * DSET_NY(anat) * DSET_NZ(anat);
00942
00943
00944 /*----- Calculate center of mass of brain image -----*/
00945 center_of_mass (&cx, &cy, &cz);
00946
00947
00948 /*----- Initialize mask for non-brain voxels -----*/
00949 for (ixyz = 0; ixyz < nxyz; ixyz++)
00950 cv[ixyz] = 1;
00951
00952
00953 /*----- Segment the axial slices -----*/
00954 if ((xxorient == 'S') || (xxorient == 'I')) segment_x_slices (cv,cx,1);
00955 else if ((yyorient == 'S') || (yyorient == 'I')) segment_y_slices (cv,cy,1);
00956 else if ((zzorient == 'S') || (zzorient == 'I')) segment_z_slices (cv,cz,1);
00957 else SI_error ("Unable to determine dataset orientation");
00958
00959
00960 /*----- Segment the sagittal slices -----*/
00961 if ((xxorient == 'L') || (xxorient == 'R')) segment_x_slices (cv,cx,0);
00962 else if ((yyorient == 'L') || (yyorient == 'R')) segment_y_slices (cv,cy,0);
00963 else if ((zzorient == 'L') || (zzorient == 'R')) segment_z_slices (cv,cz,0);
00964 else SI_error ("Unable to determine dataset orientation");
00965
00966
00967 /*----- Segment the coronal slices -----*/
00968 if ((xxorient == 'A') || (xxorient == 'P')) segment_x_slices (cv,cx,0);
00969 else if ((yyorient == 'A') || (yyorient == 'P')) segment_y_slices (cv,cy,0);
00970 else if ((zzorient == 'A') || (zzorient == 'P')) segment_z_slices (cv,cz,0);
00971 else SI_error ("Unable to determine dataset orientation");
00972
00973
00974 /*----- Create envelope for segmented image -----*/
00975 if (! nosmooth)
00976 segment_envelope (cv, cx, cy, cz);
00977
00978
00979 return;
00980 }
|
|
||||||||||||||||
|
Definition at line 211 of file Intracranial.c. References draw_2dfiller(), DSET_BRICK_ARRAY, DSET_NX, DSET_NY, DSET_NZ, free, malloc, nz, and rand_uniform(). Referenced by segment_volume().
00217 {
00218 const int MAXNPTS = 1000; /* max. number of random initial test points */
00219 const int MAXFOUND = 10; /* max. number of initial search points */
00220
00221 int nx, ny, nz, nxy, nxz, nyz, nxyz; /* dataset dimensions in voxels */
00222 int ix, jy, kz, ixy, ixz, iyz, ixyz; /* voxel indices */
00223
00224 int m; /* slice index */
00225 short * anat_data = NULL; /* data from anatomical image */
00226 float z; /* voxel gray-scale intensity */
00227 int nmpts, found; /* random search counters */
00228 byte * slice = NULL; /* current slice brain mask */
00229 byte * pslice = NULL; /* previous slice non-brain mask */
00230
00231 int count; /* count of brain voxels in current slice */
00232 int maxcount = 0; /* maximum number of brain voxels in a slice */
00233 float threshold = 0.05; /* stop if count falls below 5% of maximum */
00234
00235
00236 /*----- Progress report -----*/
00237 if (! quiet) printf ("\nSegmenting slices perpendicular to x-axis \n");
00238
00239
00240 /*----- Initialize local variables -----*/
00241 anat_data = (short *) DSET_BRICK_ARRAY(anat,0) ;
00242 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00243 nxy = nx*ny; nxz = nx*nz; nyz = ny*nz; nxyz = nxy*nz;
00244
00245
00246 /*----- Allocate memory -----*/
00247 slice = (byte *) malloc (nyz * sizeof(byte));
00248 pslice = (byte *) malloc (nyz * sizeof(byte));
00249
00250
00251 /*----- Loop over x-slices -----*/
00252 m = 0;
00253 while (m <= nx)
00254 {
00255 /*----- Start from center slice -----*/
00256 if (m <= cx) ix = cx - m;
00257 else ix = m-1;
00258
00259
00260 /*----- Initialize mask -----*/
00261 if (ix == cx)
00262 {
00263 if (axial_slice)
00264 for (iyz = 0; iyz < nyz; iyz++) pslice[iyz] = 0;
00265 else
00266 {
00267 for (kz = 0; kz < nz; kz++)
00268 for (jy = 0; jy < ny; jy++)
00269 {
00270 iyz = jy + kz*ny;
00271 ixyz = ix + jy*nx + kz*nxy;
00272 if (cv[ixyz]) pslice[iyz] = 2;
00273 else pslice[iyz] = 0;
00274 }
00275 }
00276 }
00277
00278
00279 /*----- Examine each voxel for possible inclusion in brain -----*/
00280 for (kz = 0; kz < nz; kz++)
00281 for (jy = 0; jy < ny; jy++)
00282 {
00283 iyz = jy + kz*ny;
00284 ixyz = ix + jy*nx + kz*nxy;
00285 z = anat_data[ixyz];
00286
00287 if (pslice[iyz] == 2) slice[iyz] = 1;
00288 else if (z < min_val_float) slice[iyz] = 1;
00289 else if (z > max_val_float) slice[iyz] = 1;
00290 else slice[iyz] = 0;
00291 }
00292
00293
00294 /*----- Fill brain voxels from center out -----*/
00295 draw_2dfiller (ny, nz, ny/2, nz/2, slice);
00296
00297
00298 /*----- Use additional random starting points for robustness -----*/
00299 nmpts = 0; found = 0;
00300 while ((nmpts < MAXNPTS) && (found < MAXFOUND))
00301 {
00302 nmpts++;
00303 if ((ix == cx) && axial_slice)
00304 {
00305 jy = (3*ny/8) + (int) (rand_uniform(0.0,1.0) * ny/4);
00306 kz = (3*nz/8) + (int) (rand_uniform(0.0,1.0) * nz/4);
00307 }
00308 else
00309 {
00310 jy = (int) (rand_uniform(0.0,1.0) * ny);
00311 kz = (int) (rand_uniform(0.0,1.0) * nz);
00312 }
00313 if (slice[jy+kz*ny] != 1) found++;
00314 draw_2dfiller (ny, nz, jy, kz, slice);
00315 }
00316
00317
00318 /*----- Filled brain voxels become barriers -----*/
00319 for (iyz = 0; iyz < nyz; iyz++)
00320 if (slice[iyz] == 2) pslice[iyz] = 1;
00321 else pslice[iyz] = 0;
00322
00323
00324 /*----- Now fill non-brain voxels from outside in -----*/
00325 draw_2dfiller (ny, nz, 0, 0, pslice);
00326 draw_2dfiller (ny, nz, ny-1, 0, pslice);
00327 draw_2dfiller (ny, nz, 0, nz-1, pslice);
00328 draw_2dfiller (ny, nz, ny-1, nz-1, pslice);
00329
00330
00331 /*----- Save results for this slice -----*/
00332 count = 0;
00333 for (kz = 0; kz < nz; kz++)
00334 for (jy = 0; jy < ny; jy++)
00335 {
00336 iyz = jy + kz*ny;
00337 ixyz = ix + jy*nx + kz*nxy;
00338 if (pslice[iyz] != 2)
00339 {
00340 cv[ixyz] = 0;
00341 count++;
00342 }
00343 }
00344
00345
00346 /*----- Slightly increase size of brain mask for next slice -----*/
00347 for (kz = 1; kz < nz-1; kz++)
00348 for (jy = 1; jy < ny-1; jy++)
00349 {
00350 iyz = jy + kz*ny;
00351 ixyz = ix + jy*nx + kz*nxy;
00352 if (!cv[ixyz])
00353 {
00354 pslice[iyz] = 0;
00355 pslice[iyz+1] = 0;
00356 pslice[iyz-1] = 0;
00357 pslice[iyz-ny] = 0;
00358 pslice[iyz+ny] = 0;
00359 }
00360 }
00361
00362
00363 /*----- Examine count of brain voxels for this slice -----*/
00364 if (count > maxcount) maxcount = count;
00365 if (count < threshold * maxcount)
00366 {
00367 if (m <= cx) m = cx;
00368 else m = nx;
00369 }
00370
00371
00372 /*----- Increment slice index -----*/
00373 m++;
00374
00375 } /* Loop over x-slices */
00376
00377
00378 /*----- Release memory -----*/
00379 free (pslice); pslice = NULL;
00380 free (slice); slice = NULL;
00381
00382
00383 return;
00384
00385 }
|
|
||||||||||||||||
|
Definition at line 394 of file Intracranial.c. References draw_2dfiller(), DSET_BRICK_ARRAY, DSET_NX, DSET_NY, DSET_NZ, free, malloc, nz, and rand_uniform(). Referenced by segment_volume().
00400 {
00401 const int MAXNPTS = 1000; /* max. number of random initial test points */
00402 const int MAXFOUND = 10; /* max. number of initial search points */
00403
00404 int nx, ny, nz, nxy, nxz, nyz, nxyz; /* dataset dimensions in voxels */
00405 int ix, jy, kz, ixy, ixz, iyz, ixyz; /* voxel indices */
00406
00407 int m; /* slice index */
00408 short * anat_data = NULL; /* data from anatomical image */
00409 float z; /* voxel gray-scale intensity */
00410 int nmpts, found; /* random search counters */
00411 byte * slice = NULL; /* current slice brain mask */
00412 byte * pslice = NULL; /* previous slice non-brain mask */
00413
00414 int count; /* count of brain voxels in current slice */
00415 int maxcount = 0; /* maximum number of brain voxels in a slice */
00416 float threshold = 0.05; /* stop if count falls below 5% of maximum */
00417
00418
00419 /*----- Progress report -----*/
00420 if (! quiet) printf ("\nSegmenting slices perpendicular to y-axis \n");
00421
00422
00423 /*----- Initialize local variables -----*/
00424 anat_data = (short *) DSET_BRICK_ARRAY(anat,0) ;
00425 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00426 nxy = nx*ny; nxz = nx*nz; nyz = ny*nz; nxyz = nxy*nz;
00427
00428
00429 /*----- Allocate memory -----*/
00430 slice = (byte *) malloc (nxz * sizeof(byte));
00431 pslice = (byte *) malloc (nxz * sizeof(byte));
00432
00433
00434 /*----- Loop over y-slices -----*/
00435 m = 0;
00436 while (m <= ny)
00437 {
00438 /*----- Start from center slice -----*/
00439 if (m <= cy) jy = cy - m;
00440 else jy = m-1;
00441
00442
00443 /*----- Initialize mask -----*/
00444 if (jy == cy)
00445 {
00446 if (axial_slice)
00447 for (ixz = 0; ixz < nxz; ixz++) pslice[ixz] = 0;
00448 else
00449 {
00450 for (kz = 0; kz < nz; kz++)
00451 for (ix = 0; ix < nx; ix++)
00452 {
00453 ixz = ix + kz*nx;
00454 ixyz = ix + jy*nx + kz*nxy;
00455 if (cv[ixyz]) pslice[ixz] = 2;
00456 else pslice[ixz] = 0;
00457 }
00458 }
00459 }
00460
00461
00462 /*----- Examine each voxel for possible inclusion in brain -----*/
00463 for (kz = 0; kz < nz; kz++)
00464 for (ix = 0; ix < nx; ix++)
00465 {
00466 ixz = ix + kz*nx;
00467 ixyz = ix + jy*nx + kz*nxy;
00468 z = anat_data[ixyz];
00469
00470 if (pslice[ixz] == 2) slice[ixz] = 1;
00471 else if (z < min_val_float) slice[ixz] = 1;
00472 else if (z > max_val_float) slice[ixz] = 1;
00473 else slice[ixz] = 0;
00474 }
00475
00476
00477 /*----- Fill brain voxels from center out -----*/
00478 draw_2dfiller (nx, nz, nx/2, nz/2, slice);
00479
00480
00481 /*----- Use additional random starting points for robustness -----*/
00482 nmpts = 0; found = 0;
00483 while ((nmpts < MAXNPTS) && (found < MAXFOUND))
00484 {
00485 nmpts++;
00486 if ((jy == cy) && axial_slice)
00487 {
00488 ix = (3*nx/8) + (int) (rand_uniform(0.0,1.0) * nx/4);
00489 kz = (3*nz/8) + (int) (rand_uniform(0.0,1.0) * nz/4);
00490 }
00491 else
00492 {
00493 ix = (int) (rand_uniform(0.0,1.0) * nx);
00494 kz = (int) (rand_uniform(0.0,1.0) * nz);
00495 }
00496 if (slice[ix+kz*nx] != 1) found++;
00497 draw_2dfiller (nx, nz, ix, kz, slice);
00498 }
00499
00500
00501 /*----- Filled brain voxels become barriers -----*/
00502 for (ixz = 0; ixz < nxz; ixz++)
00503 if (slice[ixz] == 2) pslice[ixz] = 1;
00504 else pslice[ixz] = 0;
00505
00506
00507 /*----- Now fill non-brain voxels from outside in -----*/
00508 draw_2dfiller (nx, nz, 0, 0, pslice);
00509 draw_2dfiller (nx, nz, nx-1, 0, pslice);
00510 draw_2dfiller (nx, nz, 0, nz-1, pslice);
00511 draw_2dfiller (nx, nz, nx-1, nz-1, pslice);
00512
00513
00514 /*----- Save results for this slice -----*/
00515 count = 0;
00516 for (kz = 0; kz < nz; kz++)
00517 for (ix = 0; ix < nx; ix++)
00518 {
00519 ixz = ix + kz*nx;
00520 ixyz = ix + jy*nx + kz*nxy;
00521 if (pslice[ixz] != 2)
00522 {
00523 cv[ixyz] = 0;
00524 count++;
00525 }
00526 }
00527
00528
00529 /*----- Slightly increase size of brain mask for next slice -----*/
00530 for (kz = 1; kz < nz-1; kz++)
00531 for (ix = 1; ix < nx-1; ix++)
00532 {
00533 ixz = ix + kz*nx;
00534 ixyz = ix + jy*nx + kz*nxy;
00535 if (!cv[ixyz])
00536 {
00537 pslice[ixz] = 0;
00538 pslice[ixz+1] = 0;
00539 pslice[ixz-1] = 0;
00540 pslice[ixz-nx] = 0;
00541 pslice[ixz+nx] = 0;
00542 }
00543 }
00544
00545
00546 /*----- Examine count of brain voxels for this slice -----*/
00547 if (count > maxcount) maxcount = count;
00548 if (count < threshold * maxcount)
00549 {
00550 if (m <= cy) m = cy;
00551 else m = ny;
00552 }
00553
00554
00555 /*----- Increment slice index -----*/
00556 m++;
00557
00558 } /* Loop over y-slices */
00559
00560
00561 /*----- Release memory -----*/
00562 free (pslice); pslice = NULL;
00563 free (slice); slice = NULL;
00564
00565
00566 return;
00567
00568 }
|
|
||||||||||||||||
|
Definition at line 577 of file Intracranial.c. References draw_2dfiller(), DSET_BRICK_ARRAY, DSET_NX, DSET_NY, DSET_NZ, free, malloc, nz, and rand_uniform(). Referenced by segment_volume().
00583 {
00584 const int MAXNPTS = 1000; /* max. number of random initial test points */
00585 const int MAXFOUND = 10; /* max. number of initial search points */
00586
00587 int nx, ny, nz, nxy, nxz, nyz, nxyz; /* dataset dimensions in voxels */
00588 int ix, jy, kz, ixy, ixz, iyz, ixyz; /* voxel indices */
00589
00590 int m; /* slice index */
00591 short * anat_data = NULL; /* data from anatomical image */
00592 float z; /* voxel gray-scale intensity */
00593 int nmpts, found; /* random search counters */
00594 byte * slice = NULL; /* current slice brain mask */
00595 byte * pslice = NULL; /* previous slice non-brain mask */
00596
00597 int count; /* count of brain voxels in current slice */
00598 int maxcount = 0; /* maximum number of brain voxels in a slice */
00599 float threshold = 0.05; /* stop if count falls below 5% of maximum */
00600
00601
00602 /*----- Progress report -----*/
00603 if (! quiet) printf ("\nSegmenting slices perpendicular to z-axis \n");
00604
00605
00606 /*----- Initialize local variables -----*/
00607 anat_data = (short *) DSET_BRICK_ARRAY(anat,0) ;
00608 nx = DSET_NX(anat); ny = DSET_NY(anat); nz = DSET_NZ(anat);
00609 nxy = nx*ny; nxz = nx*nz; nyz = ny*nz; nxyz = nxy*nz;
00610
00611
00612 /*----- Allocate memory -----*/
00613 slice = (byte *) malloc (nxy * sizeof(byte));
00614 pslice = (byte *) malloc (nxy * sizeof(byte));
00615
00616
00617 /*----- Loop over z-slices -----*/
00618 m = 0;
00619 while (m <= nz)
00620 {
00621 /*----- Start from center slice -----*/
00622 if (m <= cz) kz = cz - m;
00623 else kz = m-1;
00624
00625
00626 /*----- Initialize mask -----*/
00627 if (kz == cz)
00628 {
00629 if (axial_slice)
00630 for (ixy = 0; ixy < nxy; ixy++) pslice[ixy] = 0;
00631 else
00632 {
00633 for (ix = 0; ix < nx; ix++)
00634 for (jy = 0; jy < ny; jy++)
00635 {
00636 ixy = ix + jy*nx;
00637 ixyz = ix + jy*nx + kz*nxy;
00638 if (cv[ixyz]) pslice[ixy] = 2;
00639 else pslice[ixy] = 0;
00640 }
00641 }
00642 }
00643
00644
00645 /*----- Examine each voxel for possible inclusion in brain -----*/
00646 for (jy = 0; jy < ny; jy++)
00647 for (ix = 0; ix < nx; ix++)
00648 {
00649 ixy = ix + jy*nx;
00650 ixyz = ix + jy*nx + kz*nxy;
00651 z = anat_data[ixyz];
00652
00653 if (pslice[ixy] == 2) slice[ixy] = 1;
00654 else if (z < min_val_float) slice[ixy] = 1;
00655 else if (z > max_val_float) slice[ixy] = 1;
00656 else slice[ixy] = 0;
00657 }
00658
00659
00660 /*----- Fill brain voxels from center out -----*/
00661 draw_2dfiller (nx, ny, nx/2, ny/2, slice);
00662
00663
00664 /*----- Use additional random starting points for robustness -----*/
00665 nmpts = 0; found = 0;
00666 while ((nmpts < MAXNPTS) && (found < MAXFOUND))
00667 {
00668 nmpts++;
00669 if ((kz == cz) && axial_slice)
00670 {
00671 ix = (3*nx/8) + (int) (rand_uniform(0.0,1.0) * nx/4);
00672 jy = (3*ny/8) + (int) (rand_uniform(0.0,1.0) * ny/4);
00673 }
00674 else
00675 {
00676 ix = (int) (rand_uniform(0.0,1.0) * nx);
00677 jy = (int) (rand_uniform(0.0,1.0) * ny);
00678 }
00679 if (slice[ix+jy*nx] != 1) found++;
00680 draw_2dfiller (nx, ny, ix, jy, slice);
00681 }
00682
00683
00684 /*----- Filled brain voxels become barriers -----*/
00685 for (ixy = 0; ixy < nxy; ixy++)
00686 if (slice[ixy] == 2) pslice[ixy] = 1;
00687 else pslice[ixy] = 0;
00688
00689
00690 /*----- Now fill non-brain voxels from outside in -----*/
00691 draw_2dfiller (nx, ny, 0, 0, pslice);
00692 draw_2dfiller (nx, ny, nx-1, 0, pslice);
00693 draw_2dfiller (nx, ny, 0, ny-1, pslice);
00694 draw_2dfiller (nx, ny, nx-1, ny-1, pslice);
00695
00696
00697 /*----- Save results for this slice -----*/
00698 count = 0;
00699 for (jy = 0; jy < ny; jy++)
00700 for (ix = 0; ix < nx; ix++)
00701 {
00702 ixy = ix + jy*nx;
00703 ixyz = ix + jy*nx + kz*nxy;
00704 if (pslice[ixy] != 2)
00705 {
00706 cv[ixyz] = 0;
00707 count++;
00708 }
00709 }
00710
00711
00712 /*----- Slightly increase size of brain mask for next slice -----*/
00713 for (jy = 1; jy < ny-1; jy++)
00714 for (ix = 1; ix < nx-1; ix++)
00715 {
00716 ixy = ix + jy*nx;
00717 ixyz = ix + jy*nx + kz*nxy;
00718 if (!cv[ixyz])
00719 {
00720 pslice[ixy] = 0;
00721 pslice[ixy+1] = 0;
00722 pslice[ixy-1] = 0;
00723 pslice[ixy-nx] = 0;
00724 pslice[ixy+nx] = 0;
00725 }
00726 }
00727
00728
00729 /*----- Examine count of brain voxels for this slice -----*/
00730 if (count > maxcount) maxcount = count;
00731 if (count < threshold * maxcount)
00732 {
00733 if (m <= cz) m = cz;
00734 else m = nz;
00735 }
00736
00737
00738 /*----- Increment slice index -----*/
00739 m++;
00740
00741 } /* Loop over z-slices */
00742
00743
00744 /*----- Release memory -----*/
00745 free (pslice); pslice = NULL;
00746 free (slice); slice = NULL;
00747
00748
00749 return;
00750
00751 }
|
|
|
Definition at line 28 of file Intracranial.c. References check_one_output_file(), DSET_BRICK_TYPE, EQUIV_DATAXES, ISANAT, and SI_error(). Referenced by initialize_program().
00029 {
00030 char message[MAX_STRING_LENGTH]; /* error message */
00031
00032
00033 /*----- Check for required datasets -----*/
00034 if (anat == NULL)
00035 {
00036 sprintf (message, "No anatomical image");
00037 SI_error (message);
00038 return (0);
00039 }
00040 if (! ISANAT(anat))
00041 {
00042 sprintf (message, "Input dataset must be anatomical image");
00043 SI_error (message);
00044 return (0);
00045 }
00046 if (DSET_BRICK_TYPE(anat,0) != MRI_short)
00047 {
00048 sprintf (message, "Anatomical image must have data type: short integer");
00049 SI_error (message);
00050 return (0);
00051 }
00052
00053 #ifdef PLUG_INTRACRANIAL
00054 if (dset == NULL)
00055 {
00056 sprintf (message, "No output dataset");
00057 SI_error (message);
00058 return (0);
00059 }
00060 if (DSET_BRICK_TYPE(dset,0) != MRI_short)
00061 {
00062 sprintf (message, "Output dataset must have data type: short integer");
00063 SI_error (message);
00064 return (0);
00065 }
00066 #else
00067 if (prefix_filename == NULL)
00068 {
00069 sprintf (message, "Must specify output prefix filename");
00070 SI_error (message);
00071 return (0);
00072 }
00073
00074 /*----- Check whether output file already exists -----*/
00075 check_one_output_file (prefix_filename);
00076 #endif
00077
00078
00079 /*----- Check compatibility of datasets -----*/
00080 #ifdef PLUG_INTRACRANIAL
00081 if (! EQUIV_DATAXES (anat->daxes, dset->daxes))
00082 {
00083 sprintf (message,
00084 "Anatomical image and output dataset are incompatible");
00085 SI_error (message);
00086 return (0);
00087 }
00088 #endif
00089
00090
00091 /*----- Check allowed range of intensity values -----*/
00092 if (min_val_float >= max_val_float)
00093 {
00094 sprintf (message, "min_val >= max_val ???");
00095 SI_error (message);
00096 return (0);
00097 }
00098
00099 return (1);
00100
00101 }
|
