1#!/bin/tcsh 2# @measure_erosion_thick 3# erode a dataset repeatedly to separate datasets until nothing left 4# add up all the erosions at the end to get erosion_levels depth dataset 5# and normalize to get erosionnorm dataset 6# 12 Dec 2019 - DRG, addition of centers 7set ver = 1.01 8 9# h_view help instead of regular help (opens editor with -help output) 10@global_parse `basename $0` "$*" ; if ($status) exit 0 11 12# initializations 13set maskset = "" 14set surfset = "" 15set thickdir = "erosion_thickdir" 16set maxthick = "6" 17set resample = "auto" 18set surfsmooth = "8" 19set smoothmm = "" 20set depthsearch = "3" 21set inval = "" 22set outval = "" 23set maskval = "" 24set keep_temp_files = "" 25set balls_only = "" 26set heat_method = "HEAT_07" 27 28if ($# < 1) then 29# help 30HELP: 31 echo "@measure_erosion_thick - compute thickness of mask using erosion method" 32 echo "usage:" 33 echo "@measure_erosion_thick -maskset maskset -surfset surfacedset.gii -outdir thickdir" 34 echo 35 echo "where maskset is the dataset to find thickness" 36 echo " using the largest non-zero value in the mask." 37 echo " If dataset has values -2,-1 and 1 for different regions, this script" 38 echo " calculates the thickness only for voxels with a value of 1" 39 echo "surfset is a surface to use to find normals into the volume" 40 echo "output is in directory thickdir. If not specified, erosion_thickdir is used" 41 echo 42 echo "This script finds thickness by eroding voxels away, using facing voxels" 43 echo "a layer at a time." 44 echo 45 echo "Because of limitations in the growth of the spheres used in this method," 46 echo "it is recommended to use oversampled data, particularly when using 1mm data" 47 echo "See -resample option below" 48 echo 49 echo "Main options:" 50 echo " -maskset mydset mask dataset for input" 51 echo " -surfset mydset.gii surface dataset onto which to map thickness" 52 echo " (probably a pial/gray matter surface)" 53 echo " -outdir thickdir output directory" 54 echo 55 echo "Other options:" 56 echo 57 echo " -resample mm resample input to mm in millimeters (put a number here)" 58 echo ' set this to half a voxel or \"auto\".' 59 echo " No resampling is done by default" 60 echo " Resampling is highly recommended for most 1mm data" 61 echo " -surfsmooth mm smooth surface map of thickness by mm millimeters" 62 echo " default is 8 mm" 63 echo " -smoothmm mm smooth volume by mm FWHM in mask" 64 echo " default is 2*voxelsize of mask or resampled mask" 65 echo " -maxthick mm search for maximum thickness value of mm millimeters" 66 echo " default is 6 mm" 67 echo " -depthsearch mm map to surface by looking for max along mm millimeter" 68 echo " normal vectors. default is 3 mm" 69 echo " -keep_temp_files do not delete the intermediate files (for testing)" 70 echo " -surfsmooth_method heattype heat method used for smoothing surfaces" 71 echo " default is HEAT_07 but HEAT_05 is also useful for models" 72 echo 73 echo "Output:" 74 echo " erosion_depth.nii.gz - depth dataset" 75 echo " erosion_thick.nii.gz - volumetric thickness dataset" 76 echo " erosion_thick_smooth.nii.gz - smoothed volumetric thickness dataset" 77 echo " erosion_thick.niml.dset - unsmoothed thickness mapped to surface nodes" 78 echo " erosion_thick_smooth_nn_mm.niml.dset - smoothed thickness mapped to surface nodes" 79 echo 80 echo " Other datasets included in output:" 81 echo " maskset.nii.gz, maskset_rs.nii.gz - mask and optional resampled mask" 82 echo " anat.gii - surface representation of mask volume" 83 echo " quick.spec - simple specification file for surface to use with suma commands" 84 echo 85 echo "See related scripts and programs for computing thickness:" 86 echo " @measure_in2out, @measure_bb_thick and SurfMeasures" 87 exit 88endif 89 90 91# process user options 92set ac = 1 93while ( $ac <= $#argv ) 94 # maybe just a desperate plea for help 95 if ( ( "$argv[$ac]" == "-help" ) || ( "$argv[$ac]" == "-HELP" ) || \ 96 ( "$argv[$ac]" == "--help" ) || ( "$argv[$ac]" == "-h" ) ) then 97 goto HELP 98 else if ( "$argv[$ac]" == "-ver" ) then 99 echo "@measure_erosion_thick version: $ver" 100 exit 0 101 # get the basics - input dataset, surface and output directory 102 # only the input is really required, but the other two are nice to have 103 else if ( "$argv[$ac]" == "-maskset" ) then 104 @ ac ++ 105 if ( $ac > $#argv ) then 106 echo "** missing parameter for option '-maskset'" 107 exit 1 108 endif 109 set maskset = $argv[$ac] 110 @ ac ++; continue; 111 else if ( "$argv[$ac]" == "-surfset" ) then 112 @ ac ++ 113 if ( $ac > $#argv ) then 114 echo "** missing parameter for option '-surfset'" 115 exit 1 116 endif 117 set surfset = $argv[$ac] 118 @ ac ++; continue; 119 else if ( "$argv[$ac]" == "-outdir" ) then 120 @ ac ++ 121 if ( $ac > $#argv ) then 122 echo "** missing parameter for option '-outdir'" 123 exit 1 124 endif 125 set thickdir = $argv[$ac] 126 @ ac ++; continue; 127 128 # now tweak the options 129 else if ( "$argv[$ac]" == "-resample" ) then 130 @ ac ++ 131 if ( $ac > $#argv ) then 132 echo "** missing parameter for option '-resample'" 133 exit 1 134 endif 135 set resample = $argv[$ac] 136 if ("$resample" == "auto") then 137 echo "resampling set to auto" 138 else 139 set resample = `ccalc -expr "$resample"` 140 if ("$status" != "0") then 141 echo "resample set to $argv[$ac] is not valid" 142 exit 1 143 endif 144 echo "resample voxel size to $resample mm" 145 endif 146 @ ac ++; continue; 147 else if ( "$argv[$ac]" == "-maxthick" ) then 148 @ ac ++ 149 if ( $ac > $#argv ) then 150 echo "** missing parameter for option '-maxthick'" 151 exit 1 152 endif 153 set maxthick = $argv[$ac] 154 set maxthick = `ccalc -expr "$maxthick"` 155 if ("$status" != "0") then 156 echo "maxthick set to $argv[$ac] is not valid" 157 exit 1 158 endif 159 @ ac ++; continue; 160 else if ( "$argv[$ac]" == "-surfsmooth" ) then 161 @ ac ++ 162 if ( $ac > $#argv ) then 163 echo "** missing parameter for option '-surfsmooth'" 164 exit 1 165 endif 166 set surfsmooth = $argv[$ac] 167 set surfsmooth = `ccalc -expr "$surfsmooth"` 168 if ("$status" != "0") then 169 echo "surfsmooth set to $argv[$ac] is not valid" 170 exit 1 171 endif 172 @ ac ++; continue; 173 else if ( "$argv[$ac]" == "-surfsmooth_method" ) then 174 @ ac ++ 175 if ( $ac > $#argv ) then 176 echo "** missing parameter for option '-surfsmooth'" 177 exit 1 178 endif 179 set heat_method = $argv[$ac] 180 @ ac ++; continue; 181 else if ( "$argv[$ac]" == "-smoothmm" ) then 182 @ ac ++ 183 if ( $ac > $#argv ) then 184 echo "** missing parameter for option '-smoothmm'" 185 exit 1 186 endif 187 set smoothmm = $argv[$ac] 188 set smoothmm = `ccalc -expr "$smoothmm"` 189 if ("$status" != "0") then 190 echo "smoothmm set to $argv[$ac] is not valid" 191 exit 1 192 endif 193 @ ac ++; continue; 194 else if ( "$argv[$ac]" == "-depthsearch" ) then 195 @ ac ++ 196 if ( $ac > $#argv ) then 197 echo "** missing parameter for option '-depthsearch'" 198 exit 1 199 endif 200 set depthsearch = $argv[$ac] 201 set depthsearch = `ccalc -expr "$depthsearch"` 202 if ("$status" != "0") then 203 echo "depthsearch set to $argv[$ac] is not valid" 204 exit 1 205 endif 206 @ ac ++; continue; 207 else if ( "$argv[$ac]" == "-keep_temp_files" ) then 208 set keep_temp_files = "keep" 209 @ ac ++; continue; 210 else if ( "$argv[$ac]" == "-ignore_unknown_options" ) then 211 set ignore_unknown_options = "yes" 212 @ ac ++; continue; 213 214 # unknown options - exit 215 else 216 if ("$ignore_unknown_options" != "yes") then 217 echo "** unknown option $argv[$ac]" 218 exit 1 219 endif 220 endif 221 @ ac ++ 222end 223 224if ("$maskset" == "") then 225 echo 'No input dataset. Please provide "-maskset mydset.nii" ' 226 exit 1 227endif 228 229# get voxel size - minimum of 2D square at least in-plane (should use cubic voxels though) 230# but this could potentially work for 2D images too ("thick" 2D images) 231set v3 = `3dinfo -d3 $maskset` 232set voxsize = `ccalc "min(abs($v3[1]),abs($v3[2]))"` 233# starting radius in mm 234set inirad = $voxsize 235 236# maximum radius in mm (not absolutely necessary for this method) 237if ("$maxthick" == "") then 238 set maxrad = 4 239else 240 set maxrad = `ccalc -expr "$maxthick/2+0.01"` 241endif 242 243# do everything in a new directory for simplicity 244mkdir $thickdir 245 246# put original dataset into output directory with specific name 247# this is easier and allows for mask range selectors, sub-brick selectors, and NIFTI 2483dTcat -prefix $thickdir/allmaskset.nii.gz -overwrite $maskset 249if ($surfset != "") then 250 ConvertSurface -i $surfset -o $thickdir/anat.gii 251endif 252 253cd $thickdir 254set maskset = maskset_rs.nii.gz 255 256# goto CENTERS 257 258# get rid of any previous results if directory already exists 259\rm -f __tt_*_erode_* 260 261set maskset = maskset.nii.gz 262# if there are multiple values in the mask dataset 263# then assume the lower values are inside and outside masks for now 264set vrange = `3dBrickStat -non-zero -min -max allmaskset.nii.gz` 265if ($vrange[1] != $vrange[2]) then 266 # can do in-out, but use a simpler mask for the other methods 267 3dTcat -overwrite -prefix maskset.nii.gz allmaskset.nii.gz"<${vrange[2]}>" 268else 269 # just a single value mask - can't do in-out calculation 270 rm maskset.nii.gz 271 3dcopy allmaskset.nii.gz maskset.nii.gz 272endif 273 274# check if resampling is desired 275if ("$resample" != "") then 276 # for the "auto" way, use 1/2 voxel size 277 if ("$resample" == "auto") then 278 set voxsize = `ccalc -expr "${voxsize}/2"` 279 else 280 set voxsize = $resample 281 endif 282 3dresample -rmode NN -prefix maskset_rs.nii.gz -overwrite \ 283 -bound_type SLAB -dxyz $voxsize $voxsize $voxsize -input $maskset 284 set maskset = maskset_rs.nii.gz 285endif 286 287set iter = `printf "%2.2d" 0` 288set iterp1 = `printf "%2.2d" 1` 289 290# put copy of original dset in 1st iteration as mask and for easier bookkeeping 2913dcalc -a $maskset -prefix __tt_erode_$iter.nii.gz -overwrite -expr 'step(a)' 292 293# get voxel size - minimum of 2D square at least in-plane (should use cubic voxels though) 294# but this will for 2D images too (maybe) 295set v3 = `3dinfo -d3 $maskset` 296set voxsize = `ccalc "min(abs($v3[1]),abs($v3[2]))"` 297set nvox = 1 298 299# iterate until no more voxels in mask 300while ( $nvox != 0 ) 301 # erode facing voxels with 3dcalc 302 3dcalc -a __tt_erode_${iter}.nii.gz \ 303 -b a+i -c a-i -d a+j -e a-j -f a+k -g a-k \ 304 -expr 'and(a,b,c,d,e,f,g)' \ 305 -prefix __tt_erode_${iterp1}.nii.gz -overwrite 306 307 set iter = $iterp1 308 set iterp1 = `ccalc -form "%2.2d" "$iterp1+1"` 309 310# @ operator is weird with 08, that is converted octal and doesn't know what to do... 311# so used ccalc instead to add 1 above 312# @ iterp1++ 313# set iterp1 = `printf "%2.2d" $iterp1` 314 315 # check for maximum thickness 316 set too_thick = `ccalc "step($iterp1*$voxsize-$maxthick)"` 317 if ($too_thick == "1") then 318 set nvox = 0 319 else 320 # count the voxels left 321 set nvox = `3dBrickStat -count -non-zero __tt_erode_${iter}.nii.gz` 322 endif 323end 324rm -rf erosion_depth_vox.nii.gz erosion_depth.nii.gz erosion_thick.nii.gz erosion_norm.nii.gz 325# add all the erosions together for a continuous map of the erosion 326# each iteration left masks always with a value of 1, so how many 1's across 327# erosion iterations is the level. Inside values get more 1's 328# up to 99 erosions allowed here with the ?? 3293dMean -sum -overwrite -prefix erosion_depth_vox.nii.gz __tt_erode_??.nii.gz 330 331# compute normalized erosion, to range from 0-1 3323dcalc -a erosion_depth_vox.nii.gz -expr "a / ($iterp1-1)" \ 333 -overwrite -prefix erosion_norm.nii.gz 334 335# convert erosion level to thickness by voxelsize 336# unfortunately there's a problem converting depth to thickness 337# for every layer there can be two possible thicknesses 338# that will give rise to that depth. 339# These examples are 1D versions with x's for each voxel 340# and a number for the maximum depth 341# x - 1 342# xx - 1 343# xxx - 2 344# xxxx - 2 345# xxxxx - 3 346# xxxxxx - 3 347# xxxxxxx - 4 348# xxxxxxxx - 4 349# 1->1.5v, 2->3.5v, 3->5.5v, 4-> 7.5 3503dcalc -a erosion_depth_vox.nii.gz -expr "step(a)*max((2*a-0.5)*$voxsize, $voxsize)" \ 351 -overwrite -prefix erosion_depth.nii.gz -datum float 352 353# Convert depth to thickness by putting balls down (same as in ball and box method) 354 355# now put spheres back down at each coordinate starting at the maximum radius 356# write out all the coordinates in the mask with the current thickness 3573dmaskdump -mask $maskset -nozero erosion_depth.nii.gz > sphere_centers_pre.1D 358 359# sort coordinates from smallest to largest thickness 360sort -n -k 4 sphere_centers_pre.1D > sphere_centers_sorted.1D 361# calculate radius again - need this for 3dUndump 3621deval -a sphere_centers_sorted.1D'[3]' -expr "(a-${voxsize})/2" > radii.1D 363# 1D file is "i j k thickness radius" sorted by radius 3641dcat sphere_centers_sorted.1D radii.1D > sphere_centers.1D 365 366# fill spheres around centers with 1/2 thickness radii and values of the thickness 367# minimum value in mask should be voxsize 3683dUndump -master $maskset -overwrite -prefix erosion_thick.nii.gz \ 369 -datum float -mask $maskset \ 370 -srad 1.0 -dval $voxsize -ijk sphere_centers.1D 371 372# smooth volume in mask 373if ("$smoothmm" == "") then 374 set smoothmm = `ccalc "$voxsize*2"` 375endif 3763dBlurInMask -overwrite -FWHM $smoothmm -prefix erosion_thick_smooth.nii.gz \ 377 -mask $maskset -input erosion_thick.nii.gz 378 379# map the thickness to a surface if we have one 380if ($surfset != "") then 381 quickspec -tn gii anat.gii -spec quick.spec 382 # how many steps along normal vector for surface mapping 383 set norm_step_size = `ccalc -int "1.5*$depthsearch/$voxsize"` 384 385 # use the erosionlevels dset and search for the max going along a normal 386 # from the surface. This requires looking a little deeper, and maybe too deep 387 3dVol2Surf -spec quick.spec -surf_A anat.gii -sv $maskset \ 388 -grid_parent erosion_depth.nii.gz -map_func max \ 389 -f_steps $norm_step_size -f_index nodes \ 390 -use_norms -norm_len -$depthsearch \ 391 -cmask "-a $maskset -expr astep(a,0.000000)" \ 392 -oob_value $voxsize -oom_value $voxsize \ 393 -out_niml erosion_thick.niml.dset -overwrite 394 395 # smooth data on surface 396 # (HEAT_05 method is more robust from quick tests than HEAT_07 method 397 # recommended in help for SurfSmooth, but HEAT_07 gives smoother 398 # results that look reasonable) 399 SurfSmooth -spec quick.spec -surf_A anat.gii -met HEAT_07 \ 400 -input erosion_thick.niml.dset -fwhm $surfsmooth \ 401 -output erosion_thick_smooth.niml.dset -overwrite 402endif 403 404CENTERS: 405# find center 406# first center of mass of original dataset 407set CM = `3dCM -mask $maskset $maskset` 408 409# find second kind of center - deepest, then closest to CM 410# maximum erosion level 411set maxdepth = `3dBrickStat -max -slow erosion_depth.nii.gz` 412# calculate distance for maximum depth voxels from CM map 4133dcalc -RAI -a erosion_depth.nii.gz -datum float -overwrite \ 414 -expr "within(a,$maxdepth-0.01,$maxdepth+0.01)*(sqrt(($CM[1]-x)^2 + ($CM[2]-y)^2 + \ 415 ($CM[3]-z)^2) + 0.0001)" \ 416 -prefix max_erosion_distance.nii.gz -overwrite 417# smallest distance in deepest erosion level 418set closest_erosion = `3dBrickStat -min -slow \ 419 -mask max_erosion_distance.nii.gz max_erosion_distance.nii.gz` 420# *** Erosion Center - get xyz of closest erosion point 421# allow some tolerance to get the distance 422set lowval = `ccalc "$closest_erosion - 0.0001"` 423set highval = `ccalc "$closest_erosion + 0.0001"` 424# get xyz value 425set erosioncm_xyz = `3dmaskdump -xyz -noijk \ 426 -mask max_erosion_distance.nii.gz \ 427 -mrange $lowval $highval max_erosion_distance.nii.gz| head -1` 428 429# find third kind of center - closest to CM in mask 430# calculate distance for maximum depth voxels from CM map 4313dcalc -RAI -a erosion_depth.nii.gz -datum float -overwrite \ 432 -expr "step(a)*(sqrt(($CM[1]-x)^2 + ($CM[2]-y)^2 + \ 433 ($CM[3]-z)^2) + 0.0001)" \ 434 -prefix cm_distance.nii.gz -overwrite 435# smallest distance in all of mask (not deepest this time) 436set closest_cm = `3dBrickStat -min -slow \ 437 -mask cm_distance.nii.gz cm_distance.nii.gz` 438# *** Closest to Center of mass - get xyz of closest voxel in mask 439# allow some tolerance to get the distance 440set lowval = `ccalc "$closest_cm - 0.0001"` 441set highval = `ccalc "$closest_cm + 0.0001"` 442# get xyz value 443set closestcm_xyz = `3dmaskdump -xyz -noijk \ 444 -mask cm_distance.nii.gz \ 445 -mrange $lowval $highval cm_distance.nii.gz| head -1` 446 447# find 4th kind of center, deeper near the voxel that is close to CM 448# at closestcm_xyz, get deepest voxel nearby 449# get thickness at closestcm_xyz voxel 450set thickclosestcm = `3dmaskdump -xyz -noijk \ 451 -mask cm_distance.nii.gz \ 452 -mrange $lowval $highval erosion_thick_smooth.nii.gz | head -1` 453# look a little more than halfway into volume (thickness/2) 454set deeprad = `ccalc "${thickclosestcm[4]}*0.51"` 455 456echo ${closestcm_xyz[1-3]} 3 > __tt_closest_center.txt 4573dresample -orient RAI -inset erosion_depth.nii.gz \ 458 -prefix __tt_erosion_depth_RAI.nii.gz -overwrite 459# make a small sphere that goes a little into the mask from the edge 4603dUndump -srad $deeprad -xyz -master __tt_erosion_depth_RAI.nii.gz -datum byte \ 461 -prefix __tt_sphere_closest_center.nii.gz -overwrite __tt_closest_center.txt 462# maximum depth nearby 463set deepclose = `3dBrickStat -max -slow \ 464 -mask __tt_sphere_closest_center.nii.gz __tt_erosion_depth_RAI.nii.gz` 465# *** find deepest - get xyz of deepest voxel in very small sphere mask 466# allow some tolerance to get the value at this depth 467set lowval = `ccalc "$deepclose - 0.0001"` 468# calc distance from max depth in nearby sphere to CM as final deciding factor 4693dcalc -RAI \ 470 -a __tt_sphere_closest_center.nii.gz -b __tt_erosion_depth_RAI.nii.gz \ 471 -datum float -overwrite \ 472 -expr "step(a)*step(b-$lowval)*(sqrt(($CM[1]-x)^2 + ($CM[2]-y)^2 + \ 473 ($CM[3]-z)^2) + 0.0001)" \ 474 -prefix deepclosest_cm_distance.nii.gz -overwrite 475# smallest distance in all of mask (not deepest this time) 476set deepclosest_cm = `3dBrickStat -min -slow \ 477 -mask deepclosest_cm_distance.nii.gz deepclosest_cm_distance.nii.gz` 478# *** Closest to Center of mass - get xyz of closest voxel in mask 479# allow some tolerance to get the distance 480set lowval = `ccalc "$deepclosest_cm - 0.0001"` 481set highval = `ccalc "$deepclosest_cm + 0.0001"` 482# get xyz value 483set deepclose_xyz = `3dmaskdump -xyz -noijk \ 484 -mask deepclosest_cm_distance.nii.gz \ 485 -mrange $lowval $highval deepclosest_cm_distance.nii.gz| head -1` 486 487# put various centers into stats 488echo "CM::$CM" > erosion_stats.txt 489echo "ErosionCM::${erosioncm_xyz[1-3]}" >> erosion_stats.txt 490echo "Erosion_dist_to_CM::$closest_erosion" >> erosion_stats.txt 491 492echo "ClosestCM::${closestcm_xyz[1-3]}" >> erosion_stats.txt 493echo "Closest_dist_to_CM::$closest_cm" >> erosion_stats.txt 494 495echo "DeepClose::${deepclose_xyz[1-3]}" >> erosion_stats.txt 496 497# make a dataset to show the centers 498# make temporary RAI file 4993dresample -overwrite -orient RAI -prefix __tt_mask_RAI.nii.gz -input $maskset 500# put centers in a text file to then put as spheres in a dataset 501echo $CM 1 > centers.txt 502echo ${erosioncm_xyz[1-3]} 2 >> centers.txt 503echo ${closestcm_xyz[1-3]} 3 >> centers.txt 504echo ${deepclose_xyz[1-3]} 4 >> centers.txt 505 506# get voxel size of mask 507set voxsize = `@GetAfniRes -min $maskset` 508# set a sphere size to put the voxel down 509set srad = `ccalc "8*$voxsize"` 510 511# make volume with volumetric spheres 5123dUndump -srad $srad -xyz -master __tt_mask_RAI.nii.gz -datum byte \ 513 -prefix __tt_spheres_coords.nii.gz -overwrite centers.txt 5143dAutobox -noclust -prefix spheres_coords.nii.gz -overwrite __tt_spheres_coords.nii.gz 5153drefit -cmap INT_CMAP spheres_coords.nii.gz 516 517# unless the user is a hoarder, trash the intermediate files 518if ("$keep_temp_files" != "keep") then 519# cleanup all the erosion iterations and the voxel level erosions 520 \rm -f __tt_* erosion_depth_vox.nii.gz sphere_centers*.1D *.smrec 521endif 522 523echo "#\!/bin/tcsh" > @show_thick 524echo "setenv SUMA_Sym_I_Range NO" >> @show_thick 525echo "suma -i anat.gii -drive_com " '"-com surf_cont -view_surf_cont y \' >> @show_thick 526echo ' -com surf_cont -load_dset erosion_thick_smooth.niml.dset \' >> @show_thick 527echo ' -com surf_cont -Dim 0.4 -com surf_cont -I_range 0 5 \' >> @show_thick 528echo ' -com viewer_cont -key b" &' >> @show_thick 529echo "To show surface with thickness map, use this command:" 530echo " cd $thickdir; tcsh @show_thick" 531echo "To see volume datasets of thickness, view erosion_thick_smooth.nii.gz in the afni GUI" 532