1#!/bin/tcsh 2# h_view help instead of regular help (opens editor with -help output) 3@global_parse `basename $0` "$*" ; if ($status) exit 0 4 5# initializations 6set maskset = "" 7set surfset = "" 8set thickdir = "bb_thickdir" 9set increment = "" 10set maxthick = "" 11set resample = "auto" 12set surfsmooth = "6" 13set smoothmm = "" 14set depthsearch = "3" 15set inval = "" 16set outval = "" 17set maskval = "" 18set keep_temp_files = "" 19set balls_only = "" 20set heat_method = "HEAT_07" 21 22if ($# < 1) then 23HELP: 24 echo "@measure_bb_thick - compute thickness of mask using ball and box method" 25 echo "usage:" 26 echo "@measure_bb_thick -maskset maskset -surfset surfacedset.gii -outdir thickdir" 27 echo 28 echo "where maskset is the dataset to find thickness" 29 echo " using the largest non-zero value in the mask." 30 echo " If dataset has values -2,-1 and 1 for different regions, this script" 31 echo " calculates the thickness only for voxels with a value of 1" 32 echo "surfset is a surface to use to find normals into the volume" 33 echo "output is in directory thickdir. If not specified, bb_thickdir is used" 34 echo 35 echo 'This script finds thickness by finding the largest sphere or cube that fits' 36 echo 'within the mask at each voxel. The cube search has little effect on' 37 echo 'surface mapping of thickness, affecting only some edges in the volume.' 38 echo 'If one is primarily interested in the surface mapping, then consider' 39 echo 'the -balls_only to skip the cube search.' 40 echo 41 echo "Because of limitations in the growth of the spheres used in this method," 42 echo "it is recommended to use oversampled data, particularly when using 1mm data" 43 echo "See -resample option below" 44 echo 45 echo "Main options:" 46 echo " -maskset mydset mask dataset for input" 47 echo " -surfset mydset.gii surface dataset onto which to map thickness" 48 echo " (probably a pial/gray matter surface)" 49 echo " -outdir thickdir output directory" 50 echo 51 echo "Other options:" 52 echo 53 echo " -resample mm resample input to mm in millimeters (put a number here)" 54 echo ' set this to half a voxel or \"auto\".' 55 echo " No resampling is done by default" 56 echo " Resampling is highly recommended for most 1mm data" 57 echo " -increment mm test thickness at increments of sub-voxel distance" 58 echo " default is 1/4 voxel minimum distance (in-plane)" 59 echo " -surfsmooth mm smooth surface map of thickness by mm millimeters" 60 echo " default is 6 mm" 61 echo " -smoothmm mm smooth volume by mm FWHM in mask" 62 echo " default is 2*voxelsize of mask or resampled mask" 63 echo " -maxthick mm search for maximum thickness value of mm millimeters" 64 echo " default is 6 mm" 65 echo " -depthsearch mm map to surface by looking for max along mm millimeter" 66 echo " normal vectors. default is 3 mm" 67 echo " -keep_temp_files do not delete the intermediate files (for testing)" 68 echo " -balls_only calculate only with spheres and skip boxes" 69 echo " -surfsmooth_method heattype heat method used for smoothing surfaces" 70 echo " default is HEAT_07 but HEAT_05 is also useful for models" 71 echo 72 echo "Output:" 73 echo " maxfill.nii.gz - thickness/depth dataset" 74 echo " bb_thick.nii.gz - volumetric thickness dataset" 75 echo " bb_thick_smooth.nii.gz - smoothed volumetric thickness dataset" 76 echo " bb_thick.niml.dset - unsmoothed thickness mapped to surface nodes" 77 echo " bb_thick_smooth.niml.dset - smoothed thickness mapped to surface nodes" 78 echo 79 echo " Other datasets included in output:" 80 echo " maskset.nii.gz, maskset_rs.nii.gz - mask and optional resampled mask" 81 echo " anat.gii - surface representation of mask volume" 82 echo " quick.spec - simple specification file for surface to use with suma commands" 83 echo 84 echo "See related scripts and programs for computing thickness:" 85 echo " @measure_in2out, @measure_erosion_thick and SurfMeasures" 86 exit 87endif 88 89# process user options 90set ac = 1 91while ( $ac <= $#argv ) 92 # maybe just a desperate plea for help 93 if ( ( "$argv[$ac]" == "-help" ) || ( "$argv[$ac]" == "-HELP" ) || \ 94 ( "$argv[$ac]" == "--help" ) || ( "$argv[$ac]" == "-h" ) ) then 95 goto HELP 96 # get the basics - input dataset, surface and output directory 97 # only the input is really required, but the other two are nice to have 98 else if ( "$argv[$ac]" == "-maskset" ) then 99 @ ac ++ 100 if ( $ac > $#argv ) then 101 echo "** missing parameter for option '-maskset'" 102 exit 1 103 endif 104 set maskset = $argv[$ac] 105 @ ac ++; continue; 106 else if ( "$argv[$ac]" == "-surfset" ) then 107 @ ac ++ 108 if ( $ac > $#argv ) then 109 echo "** missing parameter for option '-surfset'" 110 exit 1 111 endif 112 set surfset = $argv[$ac] 113 @ ac ++; continue; 114 else if ( "$argv[$ac]" == "-outdir" ) then 115 @ ac ++ 116 if ( $ac > $#argv ) then 117 echo "** missing parameter for option '-outdir'" 118 exit 1 119 endif 120 set thickdir = $argv[$ac] 121 @ ac ++; continue; 122 123 # now tweak the options 124 else if ( "$argv[$ac]" == "-resample" ) then 125 @ ac ++ 126 if ( $ac > $#argv ) then 127 echo "** missing parameter for option '-resample'" 128 exit 1 129 endif 130 set resample = $argv[$ac] 131 if ("$resample" == "auto") then 132 echo "resampling set to auto" 133 else 134 set resample = `ccalc -expr "$resample"` 135 if ("$status" != "0") then 136 echo "resample set to $argv[$ac] is not valid" 137 exit 1 138 endif 139 echo "resample voxel size to $resample mm" 140 endif 141 @ ac ++; continue; 142 else if ( "$argv[$ac]" == "-increment" ) then 143 @ ac ++ 144 if ( $ac > $#argv ) then 145 echo "** missing parameter for option '-increment'" 146 exit 1 147 endif 148 set increment = $argv[$ac] 149 set increment = `ccalc -expr "$increment"` 150 if ("$status" != "0") then 151 echo "increment set to $argv[$ac] is not valid" 152 exit 1 153 endif 154 @ ac ++; continue; 155 else if ( "$argv[$ac]" == "-maxthick" ) then 156 @ ac ++ 157 if ( $ac > $#argv ) then 158 echo "** missing parameter for option '-maxthick'" 159 exit 1 160 endif 161 set maxthick = $argv[$ac] 162 set maxthick = `ccalc -expr "$maxthick"` 163 if ("$status" != "0") then 164 echo "maxthick set to $argv[$ac] is not valid" 165 exit 1 166 endif 167 @ ac ++; continue; 168 else if ( "$argv[$ac]" == "-surfsmooth" ) then 169 @ ac ++ 170 if ( $ac > $#argv ) then 171 echo "** missing parameter for option '-surfsmooth'" 172 exit 1 173 endif 174 set surfsmooth = $argv[$ac] 175 set surfsmooth = `ccalc -expr "$surfsmooth"` 176 if ("$status" != "0") then 177 echo "surfsmooth set to $argv[$ac] is not valid" 178 exit 1 179 endif 180 @ ac ++; continue; 181 else if ( "$argv[$ac]" == "-surfsmooth_method" ) then 182 @ ac ++ 183 if ( $ac > $#argv ) then 184 echo "** missing parameter for option '-surfsmooth'" 185 exit 1 186 endif 187 set heat_method = $argv[$ac] 188 @ ac ++; continue; 189 else if ( "$argv[$ac]" == "-smoothmm" ) then 190 @ ac ++ 191 if ( $ac > $#argv ) then 192 echo "** missing parameter for option '-smoothmm'" 193 exit 1 194 endif 195 set smoothmm = $argv[$ac] 196 set smoothmm = `ccalc -expr "$smoothmm"` 197 if ("$status" != "0") then 198 echo "smoothmm set to $argv[$ac] is not valid" 199 exit 1 200 endif 201 @ ac ++; continue; 202 else if ( "$argv[$ac]" == "-depthsearch" ) then 203 @ ac ++ 204 if ( $ac > $#argv ) then 205 echo "** missing parameter for option '-depthsearch'" 206 exit 1 207 endif 208 set depthsearch = $argv[$ac] 209 set depthsearch = `ccalc -expr "$depthsearch"` 210 if ("$status" != "0") then 211 echo "depthsearch set to $argv[$ac] is not valid" 212 exit 1 213 endif 214 @ ac ++; continue; 215 else if ( "$argv[$ac]" == "-keep_temp_files" ) then 216 set keep_temp_files = "keep" 217 @ ac ++; continue; 218 else if ( "$argv[$ac]" == "-balls_only" ) then 219 set balls_only = "true" 220 @ ac ++; continue; 221 else if ( "$argv[$ac]" == "-ignore_unknown_options" ) then 222 set ignore_unknown_options = "yes" 223 @ ac ++; continue; 224 225 # unknown options - exit 226 else 227 if ("$ignore_unknown_options" != "yes") then 228 echo "** unknown option $argv[$ac]" 229 exit 1 230 endif 231 endif 232 @ ac ++ 233end 234 235if ("$maskset" == "") then 236 echo 'No input dataset. Please provide "-maskset mydset.nii" ' 237 exit 1 238endif 239 240# get voxel size - minimum of 2D square at least in-plane (should use cubic voxels though) 241# but this could potentially work for 2D images too ("thick" 2D images) 242set v3 = `3dinfo -d3 $maskset` 243set voxsize = `ccalc "min(abs($v3[1]),abs($v3[2]))"` 244# starting radius in mm 245set inirad = $voxsize 246 247# maximum radius in mm (not absolutely necessary for this method) 248if ("$maxthick" == "") then 249 set maxrad = 6 250else 251 set maxrad = `ccalc -expr "$maxthick/2+0.01"` 252endif 253 254# do everything in a new directory for simplicity 255mkdir $thickdir 256 257# put original dataset into output directory with specific name 258# this is easier and allows for mask range selectors, sub-brick selectors, and NIFTI 2593dTcat -prefix $thickdir/allmaskset.nii.gz -overwrite $maskset 260if ($surfset != "") then 261 ConvertSurface -i $surfset -o $thickdir/anat.gii 262endif 263 264cd $thickdir 265# get rid of any previous results if directory already exists 266rm fill_r*.nii.gz fill_c*.nii.gz 267set maskset = maskset.nii.gz 268# if there are multiple values in the mask dataset 269# then assume the lower values are inside and outside masks for now 270set vrange = `3dBrickStat -non-zero -min -max allmaskset.nii.gz` 271if ($vrange[1] != $vrange[2]) then 272 # can do in-out, but use a simpler mask for the other methods 273 3dTcat -overwrite -prefix maskset.nii.gz allmaskset.nii.gz"<${vrange[2]}>" 274else 275 # just a single value mask - can't do in-out calculation 276 rm maskset.nii.gz 277 3dcopy allmaskset.nii.gz maskset.nii.gz 278endif 279 280# check if resampling is desired 281if ("$resample" != "") then 282 # for the "auto" way, use 1/2 voxel size 283 if ("$resample" == "auto") then 284 set voxsize = `ccalc -expr "${voxsize}/2"` 285 else 286 set voxsize = $resample 287 endif 288 3dresample -rmode NN -prefix maskset_rs.nii.gz -overwrite \ 289 -bound_type SLAB -dxyz $voxsize $voxsize $voxsize -input $maskset 290 set maskset = maskset_rs.nii.gz 291endif 292 293# increment radial search in mm 294# if nothing set, use 1/4 voxels 295if ("$increment" == "") then 296 set radinc = $voxsize 297 298 set r1 = `ccalc "0.25*$voxsize"` 299 set r2 = `ccalc "0.5*$voxsize"` 300 set r3 = `ccalc "0.75*$voxsize"` 301 #set radoffs = ( 0.00 $edge_r $corner_r ) 302 set radoffs = ( 0.00 $r1 $r2 $r3 ) 303else 304 set radinc = $increment 305 set radoffs = ( 0.00 ) 306endif 307 308set rad = $inirad 309set radi = 0 310set donext = 1 311set radoff = 0 312# try putting spheres down at different sizes to see what's the largest one that fits 313while ($donext ) 314# took out radius offsets because thickness measurements are still 315# limited by voxel increments even if large neighborhoods could still be included 316# the smaller neighborhood still gets assigned the same thickness as the larger 317# one, and that voxel only gets assigned at the center voxel. 318 foreach radoff ($radoffs) 319 set rad = `ccalc -expr "$inirad + $radinc * $radi + $radoff"` 320 set donext = `ccalc -int "step($maxrad -$rad)"` 321 if ($donext != "0") then 322 # thickness is double radius, add voxelsize to allow for full size past centers 323 # may want to adjust for edge,corner too 324 # set rthick = `ccalc "int($rad/$voxsize+0.001)*2*$voxsize+$voxsize"` 325 set rthick = `ccalc "$rad*2+$voxsize"` 326 3dLocalstat -nbhd "SPHERE($rad)" -stat filled -mask $maskset -fillvalue $rthick \ 327 -prefix fill_r$rad.nii.gz -overwrite $maskset 328 # all spheres MUST be filled. 329 # If the spheres are too big, then no spheres are found, and we can stop growing. 330 set anytherer = `3dBrickStat -max fill_r$rad.nii.gz` 331 set anythere = `ccalc -int "step($anytherer)"` 332 if ($anythere == 0) then 333 echo maximum thickness is $rthick 334 set maxfound = $rad 335 set donext = 0 336 goto finish_sphere_search 337 endif 338 endif 339 end 340 @ radi ++ 341end 342finish_sphere_search: 343 344set rad = $inirad 345set donext = 1 346# skip cubes if balls only 347if ("$balls_only" == "true") then 348 set donext = 0 349endif 350 351# Now try putting cubes down instead at different sizes to see what's the largest one that fits 352while ($donext ) 353 # calculate approx. thickness from radius - roughly equivalent to edge 354 # consider the sqrt(2) or sqrt(3) as a factor here for diagonal 355# set rthick = `ccalc "int($rad/$voxsize+0.001)*2*$voxsize+$voxsize"` 356 set rthick = `ccalc "$rad*2+$voxsize"` 357 3dLocalstat -nbhd "RECT($rad,$rad,$rad)" -stat filled -mask $maskset -fillvalue $rthick \ 358 -prefix fill_c$rad.nii.gz -overwrite $maskset 359 set anytherer = `3dBrickStat -max fill_c$rad.nii.gz` 360 set anythere = `ccalc -int "step($anytherer)"` 361 if ($anythere) then 362 set rad = `ccalc "$rad + $radinc"` 363 set donext = `ccalc -int "step($maxrad -$rad)"` 364 else 365 echo maximum thickness of cubes is $rthick 366 set maxfoundc = $rad 367 set donext = 0 368 goto finish_cube_search 369 endif 370end 371finish_cube_search: 372 373# find the biggest sphere that fits at any voxel 374# max of spheres 375if ("$balls_only" == "true") then 376 3dMean -max -prefix tmaxfill.nii.gz -overwrite fill_r*.nii.gz 377else 378 3dMean -max -prefix _tt_maxfill_spheresonly.nii.gz -overwrite fill_r*.nii.gz 379 3dcalc -a _tt_maxfill_spheresonly.nii.gz -b $maskset \ 380 -prefix maxfill_spheresonly.nii.gz -overwrite \ 381 -expr "step(b)*max(a,$voxsize)" 382 383 # max of cubes 384 3dMean -max -prefix _tt_maxfill_cubesonly.nii.gz -overwrite fill_c*.nii.gz 385 3dcalc -a _tt_maxfill_cubesonly.nii.gz -b $maskset \ 386 -prefix maxfill_cubesonly.nii.gz -overwrite \ 387 -expr "step(b)*max(a,$voxsize)" 388 # max of spheres and cubes together 389 3dMean -max -prefix tmaxfill.nii.gz -overwrite maxfill_spheresonly.nii.gz maxfill_cubesonly.nii.gz 390endif 391 392# all voxels should have at least a thickness of a voxel size 393# (don't need a 3dLocalstat for that) 3943dcalc -a tmaxfill.nii.gz -b $maskset -expr "step(b)*max(a,$voxsize)" \ 395 -prefix maxfill.nii.gz -overwrite 396 397 398# now put spheres back down at each coordinate starting at the maximum radius 399# write out all the coordinates in the mask with the current thickness 4003dmaskdump -mask maxfill.nii.gz -nozero maxfill.nii.gz > sphere_centers_pre.1D 401 402# sort coordinates from smallest to largest thickness 403sort -n -k 4 sphere_centers_pre.1D > sphere_centers_sorted.1D 404# calculate radius again - need this for 3dUndump 4051deval -a sphere_centers_sorted.1D'[3]' -expr "(a-${voxsize})/2" > radii.1D 406# 1D file is "i j k thickness radius" sorted by radius 4071dcat sphere_centers_sorted.1D radii.1D > sphere_centers.1D 408 409if ("$balls_only" == "true") then 410 set thickvol = bb_thick.nii.gz 411else 412 set thickvol = bb_spheresonly_thick.nii.gz 413endif 414 415# fill spheres around centers with 1/2 thickness radii and values of the thickness 416# minimum value in mask should be voxsize 4173dUndump -master $maskset -overwrite -prefix $thickvol \ 418 -datum float -mask $maskset \ 419 -srad 1.0 -dval $voxsize -ijk sphere_centers.1D 420 421# do similar process for boxes 422if ("$balls_only" != "true") then 423 # put boxes back down at each coordinate starting at the maximum radius 424 # write out all the coordinates in the mask with the current thickness 425 3dmaskdump -mask maxfill_cubesonly.nii.gz -nozero maxfill_cubesonly.nii.gz > cube_centers_pre.1D 426 427 # sort coordinates from smallest to largest thickness 428 sort -n -k 4 cube_centers_pre.1D > cube_centers_sorted.1D 429 # calculate radius again - need this for 3dUndump 430 1deval -a cube_centers_sorted.1D'[3]' -expr "(a-${voxsize})/2" > radii.1D 431 # 1D file is "i j k thickness radius" sorted by radius 432 1dcat cube_centers_sorted.1D radii.1D > cube_centers.1D 433 434 # fill cubes around centers with 1/2 thickness radii and values of the thickness 435 # minimum value in mask should be voxsize 436 3dUndump -master $maskset -overwrite -prefix bb_cubesonly_thick.nii.gz \ 437 -datum float -mask $maskset -cubes \ 438 -srad 1.0 -dval $voxsize -ijk cube_centers.1D 439 440 # big boxes beat out little balls to fill corners! 441 3dMean -max -overwrite -prefix bb_thick.nii.gz bb_cubesonly_thick.nii.gz bb_spheresonly_thick.nii.gz 442endif 443 444 445# smooth these out in the volume (mean, median, mode all work, but mean gives smoothest result here) 446#3dLocalstat -stat mean -nbhd "SPHERE(1)" -mask $maskset \ 447# -prefix bb_thick_smooth.nii.gz -overwrite bb_thick.nii.gz 448 449# smooth volume in mask 450if ("$smoothmm" == "") then 451 set smoothmm = `ccalc "$voxsize*2"` 452endif 4533dBlurInMask -overwrite -FWHM $smoothmm -prefix bb_thick_smooth.nii.gz -mask $maskset -input bb_thick.nii.gz 454 455# do some surface stuff if we have a surface to work with 456# could make one with IsoSurface, but if it's already done, we can use that instead 457if ($surfset != "") then 458 quickspec -tn gii anat.gii -spec quick.spec 459 # try some surface-volume interaction methods also to extract the cortical thickness 460 # one way is to look for the mean just a little on the interior 461 # because it has already been filled by spheres, do no have to look too deep 462# the surface representation from the thickness volume is blotchier than from the 463# thickness/depth volume, so I've commented out these representations to 464# avoid user confusion with too many output datasets 465# Could do this optionally in the future 466 #3dVol2Surf -spec quick.spec -surf_A anat.gii -sv $maskset \ 467 #-grid_parent bb_thick_smooth.nii.gz -gp_index 0 -map_func max \ 468 #-f_steps 10 -f_index nodes \ 469 #-cmask "-a $maskset<$maskval> -expr astep(a,0.000000)" \ 470 #-oob_value $voxsize -oom_value $voxsize \ 471 #-use_norms -norm_len -$depthsearch \ 472 #-dnode -1 -out_niml bb_thick_max.niml.dset -overwrite 473 474 ## smooth data on surface 475 ## (HEAT_05 method is more robust from quick tests than HEAT_07 method 476 ## recommended in help for SurfSmooth) 477 #SurfSmooth -spec quick.spec -surf_A anat.gii -met HEAT_05 \ 478 #-input bb_thick_max.niml.dset -fwhm $surfsmooth \ 479 #-output bb_thick_max_smooth_${surfsmooth}_mm.niml.dset -overwrite 480 481# as long as we don't go too deep into the volume, this works well 482 # another way is to use the max radii dset and search for the max going along a normal 483 # from the surface. This requires looking a little deeper, and maybe too deep 484 485 # how many steps along normal vector for surface mapping 486 set norm_step_size = `ccalc -int "1.5*$depthsearch/$voxsize"` 487 488 3dVol2Surf -spec quick.spec -surf_A anat.gii -sv $maskset \ 489 -grid_parent maxfill.nii.gz -gp_index 0 -map_func max \ 490 -f_steps $norm_step_size -f_index nodes \ 491 -cmask "-a maxfill.nii.gz -expr astep(a,0.000000)" \ 492 -oob_value $voxsize -oom_value $voxsize \ 493 -use_norms -norm_len -$depthsearch \ 494 -dnode -1 -out_niml bb_thick.niml.dset -overwrite 495 # smooth data on surface 496 # (HEAT_05 method is more robust from quick tests than HEAT_07 method 497 # recommended in help for SurfSmooth but HEAT_07 gives smoother 498 # results that look reasonable) 499 SurfSmooth -spec quick.spec -surf_A anat.gii -met $heat_method \ 500 -input bb_thick.niml.dset -fwhm $surfsmooth \ 501 -output bb_thick_smooth.niml.dset -overwrite 502endif 503 504# unless the user is a hoarder, trash the intermediate files 505if ("$keep_temp_files" != "keep") then 506 rm fill_r*.nii.gz fill_c*.nii.gz sphere_centers*.1D allmaskset.nii.gz \ 507 radii.1D tmaxfill*.nii.gz maxfill_*.nii.gz *.smrec \ 508 bb_cubesonly_thick.nii.gz bb_spheresonly_thick.nii.gz _tt*.nii.gz 509endif 510 511echo "#\!/bin/tcsh" > @show_thick 512echo "setenv SUMA_Sym_I_Range NO" >> @show_thick 513echo "suma -i anat.gii -drive_com " '"-com surf_cont -view_surf_cont y \' >> @show_thick 514echo ' -com surf_cont -load_dset bb_thick_smooth.niml.dset \' >> @show_thick 515echo ' -com surf_cont -Dim 0.4 -com surf_cont -I_range 0 5 \' >> @show_thick 516echo ' -com viewer_cont -key b" &' >> @show_thick 517echo "To show surface with thickness map, use this command:" 518echo " cd $thickdir; tcsh @show_thick" 519echo "To see volume datasets of thickness, view bb_thick_smooth.nii.gz in the afni GUI" 520