CIVET-Chimp: An automated pipeline for MRI-based cortical surface extraction in chimpanzees.

1. #!/bin/bash
3. # Created by argbash-init v2.8.0
4. # Rearrange the order of options below according to what you would like to see in the help message.
5. # ARG_OPTIONAL_SINGLE([exclude],[e],[Mask file defining regions to exclude from classifcation, region is still corrected])
6. # ARG_OPTIONAL_SINGLE([config],[c],[Path to an alternative config file defining priors to use, use "auto" to use automatic template selection])
7. # ARG_OPTIONAL_SINGLE([logfile],[l],[Path to file to log all output])
8. # ARG_OPTIONAL_BOOLEAN([standalone],[s],[Script is run standalone so save all outputs])
9. # ARG_OPTIONAL_BOOLEAN([autocrop],[a],[Crop the final output to 10 mm around the head determined by headmask from modelspace])
10. # ARG_OPTIONAL_SINGLE([max-iterations],[],[Maximum number of iterations to run],[10])
11. # ARG_OPTIONAL_SINGLE([convergence-threshold],[],[Coeffcient of variation limit between two bias field estimates],[0.01])
12. # ARG_OPTIONAL_SINGLE([classification-prior-weight],[],[How much weight is given to prior classification proabilities during iteration],[0.25])
13. # ARG_OPTIONAL_BOOLEAN([debug],[],[Debug mode, increase verbosity further, don't cleanup])
14. # ARG_VERBOSE([v])
15. # ARG_POSITIONAL_SINGLE([input],[T1w scan to be corrected])
16. # ARG_POSITIONAL_SINGLE([output],[Output filename for corrected T1w (also used as basename for other outputs)])
17. # ARGBASH_SET_INDENT([ ])
18. # ARGBASH_SET_DELIM([ =])
19. # ARG_OPTION_STACKING([getopt])
20. # ARG_RESTRICT_VALUES([no-local-options])
21. # ARG_DEFAULTS_POS([])
22. # ARG_HELP([iterativeN4_multispectral.sh is script which performs iterative inhomogeneity (bias field) correction and classification on T1w (and optionally T2w/PDw) MRI scans])
23. # ARGBASH_GO()
24. # needed because of Argbash --> m4_ignore([
25. ### START OF CODE GENERATED BY Argbash v2.8.1 one line above ###
26. # Argbash is a bash code generator used to get arguments parsing right.
27. # Argbash is FREE SOFTWARE, see https://argbash.io for more info
30. die()
31. {
32. local _ret=$2
33. test -n "$_ret" || _ret=1
34. test "$_PRINT_HELP" = yes && print_help >&2
35. echo "$1" >&2
36. exit ${_ret}
37. }
40. evaluate_strictness()
41. {
42. [[ "$2" =~ ^-(-(exclude|config|logfile|standalone|autocrop|max-iterations|convergence-threshold|classification-prior-weight|debug|verbose|input|output|help)$|[eclsavh]) ]] && die "You have passed '$2' as a value of argument '$1', which makes it look like that you have omitted the actual value, since '$2' is an option accepted by this script. This is considered a fatal error."
43. }
46. begins_with_short_option()
47. {
48. local first_option all_short_options='eclsavh'
49. first_option="${1:0:1}"
50. test "$all_short_options" = "${all_short_options/$first_option/}" && return 1 || return 0
51. }
53. # THE DEFAULTS INITIALIZATION - POSITIONALS
54. _positionals=()
55. _arg_input=
56. _arg_output=
57. # THE DEFAULTS INITIALIZATION - OPTIONALS
58. _arg_exclude=
59. _arg_config=
60. _arg_logfile=
61. _arg_standalone="off"
62. _arg_autocrop="off"
63. _arg_max_iterations="10"
64. _arg_convergence_threshold="0.01"
65. _arg_classification_prior_weight="0.25"
66. _arg_debug="off"
67. _arg_verbose=0
70. print_help()
71. {
72. printf '%s\n' "iterativeN4_multispectral.sh is script which performs iterative inhomogeneity (bias field) correction and classification on T1w (and optionally T2w/PDw) MRI scans"
73. printf 'Usage: %s [-e|--exclude <arg>] [-c|--config <arg>] [-l|--logfile <arg>] [-s|--(no-)standalone] [-a|--(no-)autocrop] [--max-iterations <arg>] [--convergence-threshold <arg>] [--classification-prior-weight <arg>] [--(no-)debug] [-v|--verbose] [-h|--help] <input> <output>\n' "$0"
74. printf '\t%s\n' "<input>: T1w scan to be corrected"
75. printf '\t%s\n' "<output>: Output filename for corrected T1w (also used as basename for other outputs)"
76. printf '\t%s\n' "-e, --exclude: Mask file defining regions to exclude from classifcation, region is still corrected (no default)"
77. printf '\t%s\n' "-c, --config: Path to an alternative config file defining priors to use, use \"auto\" to use automatic template selection (no default)"
78. printf '\t%s\n' "-l, --logfile: Path to file to log all output (no default)"
79. printf '\t%s\n' "-s, --standalone, --no-standalone: Script is run standalone so save all outputs (off by default)"
80. printf '\t%s\n' "-a, --autocrop, --no-autocrop: Crop the final output to 10 mm around the head determined by headmask from modelspace (off by default)"
81. printf '\t%s\n' "--max-iterations: Maximum number of iterations to run (default: '10')"
82. printf '\t%s\n' "--convergence-threshold: Coeffcient of variation limit between two bias field estimates (default: '0.01')"
83. printf '\t%s\n' "--classification-prior-weight: How much weight is given to prior classification proabilities during iteration (default: '0.25')"
84. printf '\t%s\n' "--debug, --no-debug: Debug mode, increase verbosity further, don't cleanup (off by default)"
85. printf '\t%s\n' "-v, --verbose: Set verbose output (can be specified multiple times to increase the effect)"
86. printf '\t%s\n' "-h, --help: Prints help"
87. }
90. parse_commandline()
91. {
92. _positionals_count=0
93. while test $# -gt 0
94. do
95. _key="$1"
96. case "$_key" in
97. -e|--exclude)
98. test $# -lt 2 && die "Missing value for the optional argument '$_key'." 1
99. _arg_exclude="$2"
100. shift
101. evaluate_strictness "$_key" "$_arg_exclude"
102. ;;
103. --exclude=*)
104. _arg_exclude="${_key##--exclude=}"
105. evaluate_strictness "$_key" "$_arg_exclude"
106. ;;
107. -e*)
108. _arg_exclude="${_key##-e}"
109. evaluate_strictness "$_key" "$_arg_exclude"
110. ;;
111. -c|--config)
112. test $# -lt 2 && die "Missing value for the optional argument '$_key'." 1
113. _arg_config="$2"
114. shift
115. evaluate_strictness "$_key" "$_arg_config"
116. ;;
117. --config=*)
118. _arg_config="${_key##--config=}"
119. evaluate_strictness "$_key" "$_arg_config"
120. ;;
121. -c*)
122. _arg_config="${_key##-c}"
123. evaluate_strictness "$_key" "$_arg_config"
124. ;;
125. -l|--logfile)
126. test $# -lt 2 && die "Missing value for the optional argument '$_key'." 1
127. _arg_logfile="$2"
128. shift
129. evaluate_strictness "$_key" "$_arg_logfile"
130. ;;
131. --logfile=*)
132. _arg_logfile="${_key##--logfile=}"
133. evaluate_strictness "$_key" "$_arg_logfile"
134. ;;
135. -l*)
136. _arg_logfile="${_key##-l}"
137. evaluate_strictness "$_key" "$_arg_logfile"
138. ;;
139. -s|--no-standalone|--standalone)
140. _arg_standalone="on"
141. test "${1:0:5}" = "--no-" && _arg_standalone="off"
142. ;;
143. -s*)
144. _arg_standalone="on"
145. _next="${_key##-s}"
146. if test -n "$_next" -a "$_next" != "$_key"
147. then
148. { begins_with_short_option "$_next" && shift && set -- "-s" "-${_next}" "$@"; } || die "The short option '$_key' can't be decomposed to ${_key:0:2} and -${_key:2}, because ${_key:0:2} doesn't accept value and '-${_key:2:1}' doesn't correspond to a short option."
149. fi
150. ;;
151. -a|--no-autocrop|--autocrop)
152. _arg_autocrop="on"
153. test "${1:0:5}" = "--no-" && _arg_autocrop="off"
154. ;;
155. -a*)
156. _arg_autocrop="on"
157. _next="${_key##-a}"
158. if test -n "$_next" -a "$_next" != "$_key"
159. then
160. { begins_with_short_option "$_next" && shift && set -- "-a" "-${_next}" "$@"; } || die "The short option '$_key' can't be decomposed to ${_key:0:2} and -${_key:2}, because ${_key:0:2} doesn't accept value and '-${_key:2:1}' doesn't correspond to a short option."
161. fi
162. ;;
163. --max-iterations)
164. test $# -lt 2 && die "Missing value for the optional argument '$_key'." 1
165. _arg_max_iterations="$2"
166. shift
167. evaluate_strictness "$_key" "$_arg_max_iterations"
168. ;;
169. --max-iterations=*)
170. _arg_max_iterations="${_key##--max-iterations=}"
171. evaluate_strictness "$_key" "$_arg_max_iterations"
172. ;;
173. --convergence-threshold)
174. test $# -lt 2 && die "Missing value for the optional argument '$_key'." 1
175. _arg_convergence_threshold="$2"
176. shift
177. evaluate_strictness "$_key" "$_arg_convergence_threshold"
178. ;;
179. --convergence-threshold=*)
180. _arg_convergence_threshold="${_key##--convergence-threshold=}"
181. evaluate_strictness "$_key" "$_arg_convergence_threshold"
182. ;;
183. --classification-prior-weight)
184. test $# -lt 2 && die "Missing value for the optional argument '$_key'." 1
185. _arg_classification_prior_weight="$2"
186. shift
187. evaluate_strictness "$_key" "$_arg_classification_prior_weight"
188. ;;
189. --classification-prior-weight=*)
190. _arg_classification_prior_weight="${_key##--classification-prior-weight=}"
191. evaluate_strictness "$_key" "$_arg_classification_prior_weight"
192. ;;
193. --no-debug|--debug)
194. _arg_debug="on"
195. test "${1:0:5}" = "--no-" && _arg_debug="off"
196. ;;
197. -v|--verbose)
198. _arg_verbose=$((_arg_verbose + 1))
199. ;;
200. -v*)
201. _arg_verbose=$((_arg_verbose + 1))
202. _next="${_key##-v}"
203. if test -n "$_next" -a "$_next" != "$_key"
204. then
205. { begins_with_short_option "$_next" && shift && set -- "-v" "-${_next}" "$@"; } || die "The short option '$_key' can't be decomposed to ${_key:0:2} and -${_key:2}, because ${_key:0:2} doesn't accept value and '-${_key:2:1}' doesn't correspond to a short option."
206. fi
207. ;;
208. -h|--help)
209. print_help
210. exit 0
211. ;;
212. -h*)
213. print_help
214. exit 0
215. ;;
216. *)
217. _last_positional="$1"
218. _positionals+=("$_last_positional")
219. _positionals_count=$((_positionals_count + 1))
220. ;;
221. esac
222. shift
223. done
224. }
227. handle_passed_args_count()
228. {
229. local _required_args_string="'input' and 'output'"
230. test "${_positionals_count}" -ge 2 || _PRINT_HELP=yes die "FATAL ERROR: Not enough positional arguments - we require exactly 2 (namely: $_required_args_string), but got only ${_positionals_count}." 1
231. test "${_positionals_count}" -le 2 || _PRINT_HELP=yes die "FATAL ERROR: There were spurious positional arguments --- we expect exactly 2 (namely: $_required_args_string), but got ${_positionals_count} (the last one was: '${_last_positional}')." 1
232. }
235. assign_positional_args()
236. {
237. local _positional_name _shift_for=$1
238. _positional_names="_arg_input _arg_output "
240. shift "$_shift_for"
241. for _positional_name in ${_positional_names}
242. do
243. test $# -gt 0 || break
244. eval "$_positional_name=\${1}" || die "Error during argument parsing, possibly an Argbash bug." 1
245. shift
246. done
247. }
249. parse_commandline "$@"
250. handle_passed_args_count
251. assign_positional_args 1 "${_positionals[@]}"
253. # OTHER STUFF GENERATED BY Argbash
255. ### END OF CODE GENERATED BY Argbash (sortof) ### ])
256. # [ <-- needed because of Argbash
258. set -euoE pipefail
260. #Special trick to redirect all output within script into logfile
261. #https://unix.stackexchange.com/questions/145651/using-exec-and-tee-to-redirect-logs-to-stdout-and-a-log-file-in-the-same-time
262. if [[ -n ${_arg_logfile} ]]; then
263. exec > >(tee -ia ${_arg_logfile})
264. exec 2> >(tee -ia ${_arg_logfile} >&2)
265. fi
267. #If debug, print timestamps and every command run
268. if [[ ${_arg_debug} == "on" ]]; then
269. set -xT
270. set -o functrace
271. PS4='+\t '
272. fi
274. #If verbose (or debug) turn on verbose outputs for commands
275. if [[ ${_arg_verbose} -ge 1 || ${_arg_debug} == "on" ]]; then
276. N4_VERBOSE=1
277. fi
279. #Create temporary directory for work
280. tmpdir=$(mktemp -d)
282. #Setup exit trap for cleanup, don't do if debug
283. function finish() {
284. if [[ ${_arg_debug} == "off" ]]; then
285. rm -rf "${tmpdir}"
286. fi
288. }
289. trap finish EXIT
291. #Add handler for failure to show where things went wrong
292. failure() {
293. local lineno=$1
294. local msg=$2
295. echo "Failed at $lineno: $msg"
296. }
297. trap 'failure ${LINENO} "$BASH_COMMAND"' ERR
299. #Set local parallelism inherited from QBATCH
300. export ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS=${THREADS_PER_COMMAND:-$(nproc)}
301. export OMP_NUM_THREADS=${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS}
303. ### DEFAULT PRIORS ###
305. #BeAST configuration
306. BEASTLIBRARY_DIR="${QUARANTINE_PATH}/resources/BEaST_libraries/combined"
307. BEAST_CONFIG=${BEASTLIBRARY_DIR}/default.1mm.conf
309. #mni_icbm152_nlin_sym_09c priors as default
310. REGISTRATIONMODEL="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c.mnc"
311. REGISTRATIONBRAINMASK="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c_mask.mnc"
312. WMPRIOR="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_wm_tal_nlin_sym_09c.mnc"
313. GMPRIOR="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_gm_tal_nlin_sym_09c.mnc"
314. CSFPRIOR="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_csf_tal_nlin_sym_09c.mnc"
316. #Files used to define MNI space
317. RESAMPLEMODEL="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c.mnc"
318. RESAMPLEMODELBRAINMASK="${QUARANTINE_PATH}/resources/mni_icbm152_nlin_sym_09c_minc2/mni_icbm152_t1_tal_nlin_sym_09c_mask.mnc"
320. # Check config files, eventually argbash will do this
321. if [[ -n ${_arg_config} && ${_arg_config} != "auto" ]]; then
322. if [[ -r ${_arg_config} ]]; then
323. source ${_arg_config}
324. else
325. echo "iterativeN4_multispectral.sh ERROR: config file does not exist or is not readable" && exit 2
326. fi
327. fi
329. if [[ ! -d ${BEASTLIBRARY_DIR} ]]; then
330. echo "iterativeN4_multispectral.sh ERROR: ${BEASTLIBRARY_DIR} does not exist"
331. fi
333. for prior in ${REGISTRATIONMODEL} ${REGISTRATIONBRAINMASK} ${WMPRIOR} ${GMPRIOR} ${CSFPRIOR} ${RESAMPLEMODEL} ${RESAMPLEMODELBRAINMASK} ${BEAST_CONFIG}; do
334. if [[ ! -s ${prior} ]]; then
335. echo "iterativeN4_multispectral.sh ERROR: File ${prior} does not exist or is zero size" && exit 3
336. fi
337. done
339. #Setup internal variables
340. output=${_arg_output}
341. originput=${_arg_input}
343. #Internal resampled input used for processing
344. input=${tmpdir}/t1.mnc
346. function outlier_mask() {
347. #Generate an outlier mask which combines the vessel segmentation >4.75 * MAD of WM
348. local outlier_input=$1
349. local outlier_mask=$2
350. local outlier_output=$3
352. local median
353. local mad
355. median=$(mincstats -quiet -median -mask ${outlier_mask} -mask_binvalue 1 ${outlier_input})
357. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -expression "abs(A[0]-${median})" ${outlier_input} ${tmpdir}/${n}/madmap.mnc
358. mad=$(mincstats -quiet -median -mask ${outlier_mask} -mask_binvalue 1 ${tmpdir}/${n}/madmap.mnc)
360. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -unsigned -byte -expression "(((0.6745*(A[0]-${median}))/${mad})<4.75)&&(A[1]<45)?1:0" \
361. ${outlier_input} ${tmpdir}/vessels.mnc ${outlier_output}
363. }
365. function renorm() {
366. #Renormalize image 0.1%-(GM/WM mean)-99.9% to 0-32767-65535 using a classification mask
367. #Compute the percentiles using the GM/WM mask
368. #Achieved via solving a linear system to get a 2nd order polynomial remapping
369. #of the intensity values
371. local renorm_input=$1
372. local renorm_classification=$2
373. local usebrainmask="${3:-}"
375. local wmbinvalue
376. local gmbinvalue
378. if [[ -n ${usebrainmask} ]]; then
379. wmbinvalue=1
380. gmbinvalue=1
381. else
382. wmbinvalue=3
383. gmbinvalue=2
384. fi
386. #Compute the percentiles and median values of GM and WM
387. valuelow=$(mincstats -quiet -mask ${tmpdir}/headmask.mnc -mask_binvalue 1 -pctT 1 ${renorm_input})
388. valuewm=$(mincstats -quiet -median -mask ${renorm_classification} -mask_binvalue ${wmbinvalue} ${renorm_input})
389. valuegm=$(mincstats -quiet -median -mask ${renorm_classification} -mask_binvalue ${gmbinvalue} ${renorm_input})
390. valuehigh=$(mincstats -quiet -mask ${tmpdir}/headmask.mnc -mask_binvalue 1 -pctT 99 ${renorm_input})
392. #Solve the linear system of a quadratic polynomial mapping the input values to 0-32767-65535
393. mapping=($(python -c "import numpy as np; print(np.array2string(np.linalg.solve(np.array([[1, ${valuelow}, ${valuelow}**2], [1, ((${valuewm}+${valuegm})/2.0), ((${valuewm}+${valuegm})/2.0)**2], [1, ${valuehigh}, ${valuehigh}**2]]),np.array([0,32767,65535])),separator= ' ')[1:-1])"))
395. #Apply the mpapping
396. minccalc -quiet ${N4_VERBOSE:+-verbose} -short -unsigned -expression "clamp(A[0]^2*${mapping[2]} + A[0]*${mapping[1]} + ${mapping[0]},0,65535)" \
397. ${renorm_input} $(dirname ${renorm_input})/$(basename ${renorm_input} .mnc).norm.mnc
399. mv -f $(dirname ${renorm_input})/$(basename ${renorm_input} .mnc).norm.mnc ${renorm_input}
400. }
402. #Function used to do bias field correction
403. function do_N4_correct() {
404. #input fov mask weight output bias shrink classifymask
405. local n4input=$1
406. local n4initmask=$2
407. local n4brainmask=$3
408. local n4weight=$4
409. local n4corrected=$5
410. local n4bias=$6
411. local n4shrink=$7
412. local n4classifymask=$8
414. #Estimate bias field
415. N4BiasFieldCorrection ${N4_VERBOSE:+--verbose} -d 3 -s ${n4shrink} -w ${n4weight} -x ${n4initmask} \
416. -b [ 200 ] -c [ 300x300x300x300,1e-5 ] --histogram-sharpening [ 0.05,0.01,200 ] \
417. -i ${tmpdir}/${n}/t1.mnc \
418. -o [ ${n4corrected},${tmpdir}/${n}/bias2.mnc ] -r 0
420. ImageMath 3 ${tmpdir}/${n}/bias2.mnc / ${tmpdir}/${n}/bias2.mnc $(mincstats -quiet -mean ${tmpdir}/${n}/bias2.mnc)
421. ImageMath 3 ${n4bias} m ${tmpdir}/prebias.mnc ${tmpdir}/${n}/bias2.mnc
422. ImageMath 3 ${n4bias} / ${n4bias} $(mincstats -quiet -mean ${n4bias})
423. ImageMath 3 ${n4corrected} / ${n4input} ${n4bias}
424. cp -f ${n4bias} ${tmpdir}/prebias.mnc
426. renorm ${n4corrected} ${n4classifymask}
427. }
429. function iterative_precorrect() {
431. local pctTlow
432. local pctThigh
434. #Foreground/background via multi-level otsu
435. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/weight1.mnc Otsu 4 ${tmpdir}/nonzero.mnc
436. ThresholdImage 3 ${tmpdir}/${n}/weight1.mnc ${tmpdir}/${n}/weight1.mnc 2 Inf 1 0
437. ImageMath 3 ${tmpdir}/${n}/weight1.mnc GetLargestComponent ${tmpdir}/${n}/weight1.mnc
438. iMath 3 ${tmpdir}/${n}/weight1.mnc MC ${tmpdir}/${n}/weight1.mnc 2 1 ball 1
439. ImageMath 3 ${tmpdir}/${n}/weight1.mnc FillHoles ${tmpdir}/${n}/weight1.mnc 2
440. cp -f ${tmpdir}/${n}/weight1.mnc ${tmpdir}/${n}/mask1.mnc
441. ImageMath 3 ${tmpdir}/${n}/weight1.mnc m ${tmpdir}/${n}/weight1.mnc ${tmpdir}/nonzero.mnc
442. ImageMath 3 ${tmpdir}/${n}/weight1.mnc GetLargestComponent ${tmpdir}/${n}/weight1.mnc
444. #Exclude Hotspots
445. minccalc -quiet ${N4_VERBOSE:+-verbose} \
446. -expression "A[0]<$(mincstats -quiet -mask ${tmpdir}/${n}/mask1.mnc -mask_binvalue 1 -pctT 99.9 ${tmpdir}/${n}/t1.mnc)?A[1]:0" \
447. ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/weight1.mnc ${tmpdir}/${n}/weighttemp.mnc
448. mv -f ${tmpdir}/${n}/weighttemp.mnc ${tmpdir}/${n}/weight1.mnc
450. #First round of correction
451. N4BiasFieldCorrection -d 3 -i ${tmpdir}/${n}/t1.mnc -b [ 200 ] -c [ 50x50x50x50,0 ] \
452. -w ${tmpdir}/${n}/weight1.mnc -o [ ${tmpdir}/${n}/t1.mnc,${tmpdir}/${n}/bias.mnc ] -s 4 --verbose \
453. --histogram-sharpening [ 0.15,0.01,200 ] -r 0 -x ${tmpdir}/initmask.mnc
455. ImageMath 3 ${tmpdir}/${n}/bias.mnc / ${tmpdir}/${n}/bias.mnc \
456. $(mincstats -quiet -mean ${tmpdir}/${n}/bias.mnc)
457. ImageMath 3 ${tmpdir}/${n}/t1.mnc / ${input} ${tmpdir}/${n}/bias.mnc
459. #Renormalize intensity
460. pctTlow=$(mincstats -quiet -mask ${tmpdir}/${n}/mask1.mnc -mask_binvalue 1 -pctT 0.1 ${tmpdir}/${n}/t1.mnc)
461. pctThigh=$(mincstats -quiet -mask ${tmpdir}/${n}/mask1.mnc -mask_binvalue 1 -pctT 99.9 ${tmpdir}/${n}/t1.mnc)
462. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -expression "clamp(clamp(A[0]-${pctTlow},0,65535)/(${pctThigh}-${pctTlow})*65535,0,65535)" \
463. ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/t1.norm.mnc
464. mv -f ${tmpdir}/${n}/t1.norm.mnc ${tmpdir}/${n}/t1.mnc
466. #Second round Foreground/background via multi-level otsu
467. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/weight2.mnc Otsu 4 ${tmpdir}/nonzero.mnc
468. ThresholdImage 3 ${tmpdir}/${n}/weight2.mnc ${tmpdir}/${n}/weight2.mnc 2 Inf 1 0
469. ImageMath 3 ${tmpdir}/${n}/weight2.mnc GetLargestComponent ${tmpdir}/${n}/weight2.mnc
470. iMath 3 ${tmpdir}/${n}/weight2.mnc MC ${tmpdir}/${n}/weight2.mnc 3 1 ball 1
471. ImageMath 3 ${tmpdir}/${n}/weight2.mnc FillHoles ${tmpdir}/${n}/weight2.mnc 2
472. cp -f ${tmpdir}/${n}/weight2.mnc ${tmpdir}/${n}/mask2.mnc
473. cp -f ${tmpdir}/${n}/mask2.mnc ${tmpdir}/fgmask.mnc
474. ImageMath 3 ${tmpdir}/${n}/weight2.mnc m ${tmpdir}/${n}/weight2.mnc ${tmpdir}/nonzero.mnc
475. ImageMath 3 ${tmpdir}/${n}/weight2.mnc GetLargestComponent ${tmpdir}/${n}/weight2.mnc
477. pctTlow=$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -pctT 0.1 ${input})
478. pctThigh=$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -pctT 99.9 ${input})
479. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -expression "clamp(clamp(A[0]-${pctTlow},0,65535)/(${pctThigh}-${pctTlow})*65535,0,65535)" \
480. ${input} ${tmpdir}/${n}/t1.mnc
481. cp ${tmpdir}/${n}/t1.mnc ${tmpdir}/t1.renorm.mnc
482. input=${tmpdir}/t1.renorm.mnc
484. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -unsigned -byte -expression 'A[0]>1.01?1:0' ${input} ${tmpdir}/nonzero.mnc
485. ImageMath 3 ${tmpdir}/${n}/weight2.mnc m ${tmpdir}/${n}/weight2.mnc ${tmpdir}/nonzero.mnc
487. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} \
488. -expression "A[0]<$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -pctT 99.5 ${tmpdir}/${n}/t1.mnc)?A[1]:0" \
489. ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/weight2.mnc ${tmpdir}/${n}/weighttemp.mnc
490. mv -f ${tmpdir}/${n}/weighttemp.mnc ${tmpdir}/${n}/weight2.mnc
492. N4BiasFieldCorrection -d 3 -i ${tmpdir}/${n}/t1.mnc -b [ 200 ] -c [ 50x50x50x50,0 ] \
493. -w ${tmpdir}/${n}/weight2.mnc -o [ ${tmpdir}/${n}/t1.mnc,${tmpdir}/${n}/bias.mnc ] -s 4 --verbose \
494. --histogram-sharpening [ 0.15,0.01,200 ] -r 0 -x ${tmpdir}/initmask.mnc
496. ImageMath 3 ${tmpdir}/${n}/bias.mnc / ${tmpdir}/${n}/bias.mnc \
497. $(mincstats -quiet -mean ${tmpdir}/${n}/bias.mnc)
498. ImageMath 3 ${tmpdir}/${n}/t1.mnc / ${input} ${tmpdir}/${n}/bias.mnc
500. pctTlow=$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -pctT 0.1 ${tmpdir}/${n}/t1.mnc)
501. pctThigh=$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -pctT 99.9 ${tmpdir}/${n}/t1.mnc)
502. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -expression "clamp(clamp(A[0]-${pctTlow},0,65535)/(${pctThigh}-${pctTlow})*65535,0,65535)" \
503. ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/t1.norm.mnc
504. mv -f ${tmpdir}/${n}/t1.norm.mnc ${tmpdir}/${n}/t1.mnc
506. cp -f ${tmpdir}/${n}/bias.mnc ${tmpdir}/${n}/prebias.mnc
508. itk_vesselness --scales 8 --rescale ${tmpdir}/${n}/t1.mnc ${tmpdir}/vessels.mnc
510. i=0
512. #Iterative correction using tissue masking, some badly biased scans can't be
513. #corrected in one-shot
514. while true; do
516. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/otsu.mnc Otsu 4 ${tmpdir}/${n}/mask$((2 + i)).mnc
517. ThresholdImage 3 ${tmpdir}/${n}/otsu.mnc ${tmpdir}/${n}/otsu.mnc 2 Inf 1 0
519. ImageMath 3 ${tmpdir}/${n}/mask$((3 + i)).mnc GetLargestComponent ${tmpdir}/${n}/otsu.mnc
520. iMath 3 ${tmpdir}/${n}/mask$((3 + i)).mnc MC ${tmpdir}/${n}/mask$((3 + i)).mnc 8 1 ball 1
521. ImageMath 3 ${tmpdir}/${n}/mask$((3 + i)).mnc FillHoles ${tmpdir}/${n}/mask$((3 + i)).mnc 2
522. cp -f ${tmpdir}/${n}/mask$((3 + i)).mnc ${tmpdir}/fgmask.mnc
524. if [[ $i == 0 ]]; then
525. ImageMath 3 ${tmpdir}/${n}/weight$((3 + i)).mnc + ${tmpdir}/${n}/otsu.mnc ${tmpdir}/${n}/mask$((2 + i)).mnc
526. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} \
527. -expression "(A[0]<45)&&(A[1]<$(mincstats -quiet -mask ${tmpdir}/${n}/mask$((2 + i)).mnc -mask_binvalue 1 -pctT 99.5 ${tmpdir}/${n}/t1.mnc))?A[2]:0" \
528. ${tmpdir}/vessels.mnc ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/otsu.mnc ${tmpdir}/${n}/weight$((3 + i)).mnc
529. else
530. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} \
531. -expression "(A[0]<45)&&(A[1]<$(mincstats -quiet -mask ${tmpdir}/${n}/mask$((2 + i)).mnc -mask_binvalue 1 -pctT 99.5 ${tmpdir}/${n}/t1.mnc))?A[2]:0" \
532. ${tmpdir}/vessels.mnc ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/otsu.mnc ${tmpdir}/${n}/weight$((3 + i)).mnc
533. iMath 3 ${tmpdir}/${n}/weight$((3 + i)).mnc ME ${tmpdir}/${n}/otsu.mnc 1 1 ball 1
534. ImageMath 3 ${tmpdir}/${n}/weight$((3 + i)).mnc GetLargestComponent ${tmpdir}/${n}/weight$((3 + i)).mnc
535. iMath 3 ${tmpdir}/${n}/weight$((3 + i)).mnc MD ${tmpdir}/${n}/weight$((3 + i)).mnc 1 1 ball 1
536. fi
538. N4BiasFieldCorrection -d 3 -i ${tmpdir}/${n}/t1.mnc -b [ 200 ] -c [ 300x300x300x300,1e-4 ] \
539. -w ${tmpdir}/${n}/weight$((3 + i)).mnc -o [ ${tmpdir}/${n}/t1.mnc,${tmpdir}/${n}/bias.mnc ] -s 2 --verbose \
540. --histogram-sharpening [ 0.05,0.01,200 ] -r 0 -x ${tmpdir}/initmask.mnc
542. ImageMath 3 ${tmpdir}/${n}/bias.mnc / ${tmpdir}/${n}/bias.mnc $(mincstats -quiet -mean ${tmpdir}/${n}/bias.mnc)
543. ImageMath 3 ${tmpdir}/${n}/prebias.mnc m ${tmpdir}/${n}/prebias.mnc ${tmpdir}/${n}/bias.mnc
544. ImageMath 3 ${tmpdir}/${n}/prebias.mnc / ${tmpdir}/${n}/prebias.mnc $(mincstats -quiet -mean ${tmpdir}/${n}/prebias.mnc)
545. ImageMath 3 ${tmpdir}/${n}/t1.mnc / ${input} ${tmpdir}/${n}/prebias.mnc
547. ((++i))
548. [[ ( ${i} -le 2 ) ]] || break
550. pctTlow=$(mincstats -quiet -mask ${tmpdir}/${n}/mask$((2 + i)).mnc -mask_binvalue 1 -pctT 0.1 ${tmpdir}/${n}/t1.mnc)
551. pctThigh=$(mincstats -quiet -mask ${tmpdir}/${n}/mask$((2 + i)).mnc -mask_binvalue 1 -pctT 99.9 ${tmpdir}/${n}/t1.mnc)
552. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -expression "clamp(clamp(A[0]-${pctTlow},0,65535)/(${pctThigh}-${pctTlow})*65535,0,65535)" \
553. ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/t1.norm.mnc
554. mv -f ${tmpdir}/${n}/t1.norm.mnc ${tmpdir}/${n}/t1.mnc
556. done
558. pctThigh=$(mincstats -quiet -mask ${tmpdir}/${n}/mask5.mnc -mask_binvalue 1 -pctT 99.9 ${tmpdir}/${n}/t1.mnc)
559. pctTlow=$(mincstats -quiet -mask ${tmpdir}/${n}/mask5.mnc -mask_binvalue 1 -pctT 0.1 ${tmpdir}/${n}/t1.mnc)
560. minccalc -quiet ${N4_VERBOSE:+-verbose} -short -unsigned -expression "clamp(clamp(A[0]-${pctTlow},0,65535)/(${pctThigh}-${pctTlow})*65535,0,65535)" \
561. ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/corrected.mnc
563. minc_anlm ${N4_VERBOSE:+--verbose} --clobber --mt ${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS} ${tmpdir}/${n}/corrected.mnc ${tmpdir}/${n}/t1.mnc
564. }
567. function classify_to_mask() {
568. #Convert classify image into a mask
569. #Mostly a clone of the supersteps of the antsBrainExtraction supersteps
570. ThresholdImage 3 ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/gm.mnc 2 2 1 0
571. ThresholdImage 3 ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/wm.mnc 3 3 1 0
573. ImageMath 3 ${tmpdir}/${n}/gm.mnc GetLargestComponent ${tmpdir}/${n}/gm.mnc
574. ImageMath 3 ${tmpdir}/${n}/wm.mnc GetLargestComponent ${tmpdir}/${n}/wm.mnc
576. ImageMath 3 ${tmpdir}/${n}/gm.mnc FillHoles ${tmpdir}/${n}/gm.mnc 2
578. ImageMath 3 ${tmpdir}/${n}/classifymask.mnc addtozero ${tmpdir}/${n}/gm.mnc ${tmpdir}/${n}/wm.mnc
580. iMath 3 ${tmpdir}/${n}/classifymask.mnc ME ${tmpdir}/${n}/classifymask.mnc 1 1 ball 1
581. ImageMath 3 ${tmpdir}/${n}/classifymask.mnc GetLargestComponent ${tmpdir}/${n}/classifymask.mnc
582. iMath 3 ${tmpdir}/${n}/classifymask.mnc MD ${tmpdir}/${n}/classifymask.mnc 2 1 ball 1
583. iMath 3 ${tmpdir}/bmask_E.mnc ME ${tmpdir}/masks/mnimask.mnc 10 1 ball 1
584. ImageMath 3 ${tmpdir}/${n}/classifymask.mnc addtozero ${tmpdir}/${n}/classifymask.mnc ${tmpdir}/bmask_E.mnc
585. ImageMath 3 ${tmpdir}/${n}/classifymask.mnc FillHoles ${tmpdir}/${n}/classifymask.mnc 2
586. }
588. function make_qc() {
589. #Generate a standardized view of the final correct brain in MNI space, with classification overlayed
590. #Create animated version if img2webp is available
591. mkdir -p ${tmpdir}/qc
593. #Resample into MNI space for all the inputs
594. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 ${MNI_XFM:+-t ${MNI_XFM}} -t ${tmpdir}/mni0_GenericAffine.xfm \
595. -i ${tmpdir}/${n}/classify.mnc -o ${tmpdir}/qc/classify.mnc -r ${RESAMPLEMODEL} -n GenericLabel
596. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 ${MNI_XFM:+-t ${MNI_XFM}} -t ${tmpdir}/mni0_GenericAffine.xfm \
597. -i ${tmpdir}/corrected.mnc -o ${tmpdir}/qc/corrected.mnc -r ${RESAMPLEMODEL} -n BSpline[5]
598. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 ${MNI_XFM:+-t ${MNI_XFM}} -t ${tmpdir}/mni0_GenericAffine.xfm \
599. -i ${tmpdir}/origqcref.mnc -o ${tmpdir}/qc/orig.mnc -r ${RESAMPLEMODEL} -n BSpline[5]
600. mincmath -clobber -quiet ${N4_VERBOSE:+-verbose} -clamp -const2 0 65535 ${tmpdir}/qc/corrected.mnc ${tmpdir}/qc/corrected.clamp.mnc
601. mv -f ${tmpdir}/qc/corrected.clamp.mnc ${tmpdir}/qc/corrected.mnc
602. mincmath -clobber -quiet ${N4_VERBOSE:+-verbose} -clamp -const2 0 65535 ${tmpdir}/qc/orig.mnc ${tmpdir}/qc/orig.clamp.mnc
603. mv -f ${tmpdir}/qc/orig.clamp.mnc ${tmpdir}/qc/orig.mnc
605. #Create the bounding box for create_verify_image
606. mincresample -clobber -quiet ${N4_VERBOSE:+-verbose} $(mincbbox -mincresample ${tmpdir}/qc/classify.mnc) ${tmpdir}/qc/classify.mnc ${tmpdir}/qc/label-crop.mnc
607. minccalc -quiet ${N4_VERBOSE:+-verbose} -unsigned -byte -expression '1' ${tmpdir}/qc/label-crop.mnc ${tmpdir}/qc/bounding.mnc
609. #Trasverse
610. create_verify_image -range_floor 0 ${tmpdir}/qc/trans_classify.rgb \
611. -width 1920 -autocols 10 -autocol_planes t \
612. -bounding_volume ${tmpdir}/qc/bounding.mnc \
613. -row ${tmpdir}/qc/corrected.mnc color:gray:0:65535 \
614. volume_overlay:${tmpdir}/qc/classify.mnc:0.4
616. create_verify_image -range_floor 0 ${tmpdir}/qc/trans_corrected.rgb \
617. -width 1920 -autocols 10 -autocol_planes t \
618. -bounding_volume ${tmpdir}/qc/bounding.mnc \
619. -row ${tmpdir}/qc/corrected.mnc color:spect:0:65535
621. create_verify_image -range_floor 0 ${tmpdir}/qc/trans_corrected_gray.rgb \
622. -width 1920 -autocols 10 -autocol_planes t \
623. -bounding_volume ${tmpdir}/qc/bounding.mnc \
624. -row ${tmpdir}/qc/corrected.mnc color:gray:0:65535
626. create_verify_image -range_floor 0 ${tmpdir}/qc/trans_orig.rgb \
627. -width 1920 -autocols 10 -autocol_planes t \
628. -bounding_volume ${tmpdir}/qc/bounding.mnc \
629. -row ${tmpdir}/qc/orig.mnc color:spect:0:65535
631. #Sagital
632. create_verify_image -range_floor 0 ${tmpdir}/qc/sag_classify.rgb \
633. -width 1920 -autocols 10 -autocol_planes s \
634. -bounding_volume ${tmpdir}/qc/bounding.mnc \
635. -row ${tmpdir}/qc/corrected.mnc color:gray:0:65535 \
636. volume_overlay:${tmpdir}/qc/classify.mnc:0.4
638. create_verify_image -range_floor 0 ${tmpdir}/qc/sag_corrected.rgb \
639. -width 1920 -autocols 10 -autocol_planes s \
640. -bounding_volume ${tmpdir}/qc/bounding.mnc \
641. -row ${tmpdir}/qc/corrected.mnc color:spect:0:65535
643. create_verify_image -range_floor 0 ${tmpdir}/qc/sag_corrected_gray.rgb \
644. -width 1920 -autocols 10 -autocol_planes s \
645. -bounding_volume ${tmpdir}/qc/bounding.mnc \
646. -row ${tmpdir}/qc/corrected.mnc color:gray:0:65535
648. create_verify_image -range_floor 0 ${tmpdir}/qc/sag_orig.rgb \
649. -width 1920 -autocols 10 -autocol_planes s \
650. -bounding_volume ${tmpdir}/qc/bounding.mnc \
651. -row ${tmpdir}/qc/orig.mnc color:spect:0:65535
653. #Coronal
654. create_verify_image -range_floor 0 ${tmpdir}/qc/cor_classify.rgb \
655. -width 1920 -autocols 10 -autocol_planes c \
656. -bounding_volume ${tmpdir}/qc/bounding.mnc \
657. -row ${tmpdir}/qc/corrected.mnc color:gray:0:65535 \
658. volume_overlay:${tmpdir}/qc/classify.mnc:0.4
660. create_verify_image -range_floor 0 ${tmpdir}/qc/cor_corrected.rgb \
661. -width 1920 -autocols 10 -autocol_planes c \
662. -bounding_volume ${tmpdir}/qc/bounding.mnc \
663. -row ${tmpdir}/qc/corrected.mnc color:spect:0:65535
665. create_verify_image -range_floor 0 ${tmpdir}/qc/cor_corrected_gray.rgb \
666. -width 1920 -autocols 10 -autocol_planes c \
667. -bounding_volume ${tmpdir}/qc/bounding.mnc \
668. -row ${tmpdir}/qc/corrected.mnc color:gray:0:65535
670. create_verify_image -range_floor 0 ${tmpdir}/qc/cor_orig.rgb \
671. -width 1920 -autocols 10 -autocol_planes c \
672. -bounding_volume ${tmpdir}/qc/bounding.mnc \
673. -row ${tmpdir}/qc/orig.mnc color:spect:0:65535
675. convert -background black -strip -append \
676. ${tmpdir}/qc/cor_corrected.rgb \
677. ${tmpdir}/qc/cor_classify.rgb \
678. ${tmpdir}/qc/sag_corrected.rgb \
679. ${tmpdir}/qc/sag_classify.rgb \
680. ${tmpdir}/qc/trans_corrected.rgb \
681. ${tmpdir}/qc/trans_classify.rgb \
682. ${tmpdir}/qc/corrected.mpc
684. convert -background black -strip -append \
685. ${tmpdir}/qc/cor_orig.rgb \
686. ${tmpdir}/qc/cor_corrected_gray.rgb \
687. ${tmpdir}/qc/sag_orig.rgb \
688. ${tmpdir}/qc/sag_corrected_gray.rgb \
689. ${tmpdir}/qc/trans_orig.rgb \
690. ${tmpdir}/qc/trans_corrected_gray.rgb \
691. ${tmpdir}/qc/orig.mpc
693. #Save static QC jpg
694. convert -background black -strip -interlace Plane -sampling-factor 4:2:0 -quality "85%" \
695. ${tmpdir}/qc/corrected.mpc $(dirname ${output})/$(basename ${output} .mnc).jpg
697. #If webp software is available animate a before/after image
698. if command -v img2webp; then
699. convert -background black ${tmpdir}/qc/corrected.mpc ${tmpdir}/qc/corrected.png
700. convert -background black ${tmpdir}/qc/orig.mpc ${tmpdir}/qc/orig.png
701. img2webp -d 750 -lossy -min_size ${tmpdir}/qc/corrected.png ${tmpdir}/qc/orig.png -o $(dirname ${output})/$(basename ${output} .mnc).webp || true
702. fi
703. }
705. function test_templates() {
706. #Automatic template selection for most similar template for use as prior
707. #Loop over the configs/auto config files and choose the best one based on ants CC
708. mkdir -p ${tmpdir}/test_templates
710. for configfile in $(dirname "$(readlink -f "$0")")/configs/auto/*cfg; do
712. source ${configfile}
714. antsRegistration ${N4_VERBOSE:+--verbose} -d 3 --float 1 --minc \
715. --output [ ${tmpdir}/test_templates/$(basename ${configfile} .cfg),${tmpdir}/test_templates/$(basename ${configfile} .cfg).mnc ] \
716. --use-histogram-matching 1 \
717. --initial-moving-transform [ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1 ] \
718. --transform Translation[ 0.1 ] \
719. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
720. --convergence [ 500x500x500x500x500x500x500x500,1e-6,10 ] \
721. --shrink-factors 6x6x6x6x6x6x6x6 \
722. --smoothing-sigmas 6.35574237559x5.93006674681x5.50423435717x5.07820577132x4.65192708599x4.22532260674x3.79828256043x3.37064139994mm \
723. --masks [ NOMASK,NOMASK ] \
724. --transform Rigid[ 0.1 ] \
725. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
726. --convergence [ 500x500x500x500x500x500x500,1e-6,10 ] \
727. --shrink-factors 6x6x6x6x6x6x5 \
728. --smoothing-sigmas 4.65192708599x4.22532260674x3.79828256043x3.37064139994x2.94213702015x2.51232776601x2.08040503813mm \
729. --masks [ NOMASK,NOMASK ] \
730. --transform Similarity[ 0.1 ] \
731. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
732. --convergence [ 500x500x500x500x450x150,1e-6,10 ] \
733. --shrink-factors 6x6x5x4x3x2 \
734. --smoothing-sigmas 2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
735. --masks [ NOMASK,NOMASK ] \
736. --transform Similarity[ 0.1 ] \
737. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
738. --convergence [ 500x500x500x500x450x150,1e-6,10 ] \
739. --shrink-factors 6x6x5x4x3x2 \
740. --smoothing-sigmas 2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
741. --masks [ ${REGISTRATIONBRAINMASK},NOMASK ] \
742. --transform Affine[ 0.1 ] \
743. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,64,None ] \
744. --convergence [ 500x450x150x50,1e-6,10 ] \
745. --shrink-factors 4x3x2x1 \
746. --smoothing-sigmas 1.64470459404x1.20112240879x0.735534255037x0.0mm \
747. --masks [ ${REGISTRATIONBRAINMASK},NOMASK ]
749. echo ${configfile},$(MeasureImageSimilarity -d 3 -m CC[${REGISTRATIONMODEL},${tmpdir}/test_templates/$(basename ${configfile} .cfg).mnc,1,4] \
750. -x ${REGISTRATIONBRAINMASK}) >> ${tmpdir}/test_templates/results.csv
751. done
753. #Prep and load winner template
754. unset MNI_XFM
755. echo "Choosing template $(sort -k2 -g -t, ${tmpdir}/test_templates/results.csv | cut -d"," -f 1 | head -1)"
756. source $(sort -k2 -g -t, ${tmpdir}/test_templates/results.csv | cut -d"," -f 1 | head -1)
758. #Store template registration for later use
759. cp -f ${tmpdir}/test_templates/$(basename $(sort -k2 -g -t, ${tmpdir}/test_templates/results.csv | cut -d"," -f 1 | head -1) .cfg)0_GenericAffine.xfm ${tmpdir}/template_bootstrap.xfm
761. if [[ ${_arg_debug} == "off" ]]; then
762. rm -rf ${tmpdir}/test_templates
763. fi
764. }
766. ##########START OF SCRIPT#############
767. #Forceably convert to MINC2, and clamp range to avoid negative numbers, rescale to 0-65535
768. mincconvert -2 ${originput} ${tmpdir}/originput.mnc
769. #Rescale initial data into entirely positive range (fix for completely negative data)
770. ImageMath 3 ${tmpdir}/originput.mnc RescaleImage ${tmpdir}/originput.mnc 0 65535
771. #Very mild range clamp for very hot voxels
772. mincmath -quiet ${N4_VERBOSE:+-verbose} -clamp \
773. -const2 $(mincstats -quiet -floor 1e-12 -pctT 0.1 ${tmpdir}/originput.mnc) \
774. $(mincstats -quiet -floor 1e-12 -pctT 99.9 ${tmpdir}/originput.mnc) \
775. ${tmpdir}/originput.mnc ${tmpdir}/originput.clamp.mnc
776. ImageMath 3 ${tmpdir}/originput.clamp.mnc RescaleImage ${tmpdir}/originput.clamp.mnc 0 65535
777. mincresample -quiet ${N4_VERBOSE:+-verbose} -like ${tmpdir}/originput.mnc -keep -unsigned -short \
778. ${tmpdir}/originput.clamp.mnc ${tmpdir}/originput.clamp.resample.mnc
779. mv -f ${tmpdir}/originput.clamp.resample.mnc ${tmpdir}/originput.mnc
780. rm -f ${tmpdir}/originput.clamp.mnc
781. originput=${tmpdir}/originput.mnc
782. cp -f ${originput} ${tmpdir}/origqcref.mnc
784. #Isotropize, and normalize intensity range, this is the file that will be processed in the pipeline
785. #Need smoothing for downsampling to avoid aliasing
786. #Ideas stolen from https://discourse.itk.org/t/resampling-to-isotropic-signal-processing-theory/1403
787. isostep=1.0
788. inputres=$(python -c "print('\n'.join([str(abs(x)) for x in [float(x) for x in \"$(PrintHeader ${originput} 1)\".split(\"x\")]]))")
789. blurs=""
791. for dim in ${inputres}; do
792. if [[ $(python -c "print(${dim}>(${isostep}-1e-6))") == True ]]; then
793. blurs+=1e-12x
794. else
795. blurs+=$(python -c "import math; print(math.sqrt((${isostep}**2.0 - ${dim}**2.0)/(2.0*math.sqrt(2.0*math.log(2.0)))**2.0))")x
796. fi
797. done
799. SmoothImage 3 ${originput} "${blurs%?}" ${tmpdir}/smoothed.mnc 1 0
800. ResampleImage 3 ${tmpdir}/smoothed.mnc ${input} ${isostep}x${isostep}x${isostep} 0 4
801. mincmath -quiet ${N4_VERBOSE:+-verbose} -clamp -const2 0 $(mincstats -max -quiet ${input}) ${input} ${tmpdir}/input.clamp.mnc
802. ImageMath 3 ${input} RescaleImage ${tmpdir}/input.clamp.mnc 0 65535
803. rm -f ${tmpdir}/input.clamp.mnc
804. ImageMath 3 ${input} PadImage ${input} 20
806. #Generate a global nonzero mask to always exclude pure background
807. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -unsigned -byte -expression 'A[0]>1.01?1:0' ${input} ${tmpdir}/nonzero.mnc
809. #If exclusion mask exists, negate it to produce a multiplicative exlcusion mask, resample to internal resolution
810. if [[ -n ${_arg_exclude} ]]; then
811. ImageMath 3 ${tmpdir}/exclude.mnc Neg ${_arg_exclude}
812. excludemask=${tmpdir}/exclude.mnc
813. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${excludemask} -r ${input} -n GenericLabel -o ${excludemask}
814. else
815. excludemask=""
816. fi
818. mkdir -p ${tmpdir}/masks
820. ################################################################################
821. #Round 0
822. #Iterative estimation of a mask with multilevel otsu to find foreground-background
823. #Also found forground mask and use it combined to template FOV registration to
824. #trim the FOV to generate a headmask
825. ################################################################################
826. n=0
828. mkdir -p ${tmpdir}/${n}
830. minc_anlm ${N4_VERBOSE:+--verbose} --mt ${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS} ${input} ${tmpdir}/${n}/t1.mnc
832. iterative_precorrect
834. #If the "auto" template method is selected, do the registrations and CC
835. #estimate to find best matching model
836. if [[ ${_arg_config} == "auto" ]]; then
837. test_templates
838. fi
840. #Register to model to resample back a FOV mask
841. if [[ -s ${tmpdir}/template_bootstrap.xfm ]]; then
842. cp -f ${tmpdir}/template_bootstrap.xfm ${tmpdir}/${n}/mni0_GenericAffine.xfm
843. else
844. antsRegistration ${N4_VERBOSE:+--verbose} -d 3 --float 1 --minc \
845. --output [ ${tmpdir}/${n}/mni ] \
846. --use-histogram-matching 1 \
847. --initial-moving-transform [ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1 ] \
848. --transform Translation[ 0.1 ] \
849. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
850. --convergence [ 500x500x500x500x500x500x500x500,1e-6,10 ] \
851. --shrink-factors 6x6x6x6x6x6x6x6 \
852. --smoothing-sigmas 6.35574237559x5.93006674681x5.50423435717x5.07820577132x4.65192708599x4.22532260674x3.79828256043x3.37064139994mm \
853. --masks [ NOMASK,NOMASK ] \
854. --transform Rigid[ 0.1 ] \
855. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
856. --convergence [ 500x500x500x500x500x500x500,1e-6,10 ] \
857. --shrink-factors 6x6x6x6x6x6x5 \
858. --smoothing-sigmas 4.65192708599x4.22532260674x3.79828256043x3.37064139994x2.94213702015x2.51232776601x2.08040503813mm \
859. --masks [ NOMASK,NOMASK ] \
860. --transform Similarity[ 0.1 ] \
861. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
862. --convergence [ 500x500x500x500x450x150,1e-6,10 ] \
863. --shrink-factors 6x6x5x4x3x2 \
864. --smoothing-sigmas 2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
865. --masks [ NOMASK,NOMASK ] \
866. --transform Similarity[ 0.1 ] \
867. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
868. --convergence [ 500x500x500x500x450x150,1e-6,10 ] \
869. --shrink-factors 6x6x5x4x3x2 \
870. --smoothing-sigmas 2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
871. --masks [ ${REGISTRATIONBRAINMASK},NOMASK ] \
872. --transform Affine[ 0.1 ] \
873. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,64,None ] \
874. --convergence [ 500x450x150x0,1e-6,10 ] \
875. --shrink-factors 4x3x2x1 \
876. --smoothing-sigmas 1.64470459404x1.20112240879x0.735534255037x0.0mm \
877. --masks [ ${REGISTRATIONBRAINMASK},NOMASK ]
878. fi
880. #Make a fov mask from all 1's of the
881. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -unsigned -byte -expression '1' ${REGISTRATIONMODEL} ${tmpdir}/modelfovmask.mnc
882. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/modelfovmask.mnc \
883. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -o ${tmpdir}/headmask.mnc -r ${tmpdir}/${n}/t1.mnc -n GenericLabel
885. #Headmask is intersection of filled otsu foreground mask and FOV from model
886. ImageMath 3 ${tmpdir}/headmask.mnc m ${tmpdir}/headmask.mnc ${tmpdir}/fgmask.mnc
887. cp -f ${tmpdir}/fgmask.mnc ${tmpdir}/fgmask_orig.mnc
889. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} \
890. -expression "(A[0]<45)&&(A[1]<$(mincstats -quiet -mask ${tmpdir}/headmask.mnc -mask_binvalue 1 -pctT 99.5 ${tmpdir}/${n}/t1.mnc))?A[2]:0" \
891. ${tmpdir}/vessels.mnc ${tmpdir}/${n}/t1.mnc ${tmpdir}/headmask.mnc ${tmpdir}/${n}/otsu.mnc
893. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/otsu.mnc Otsu 4 ${tmpdir}/${n}/otsu.mnc
894. ThresholdImage 3 ${tmpdir}/${n}/otsu.mnc ${tmpdir}/${n}/weight6.mnc 2 Inf 1 0
896. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/otsu.mnc Otsu 4 ${tmpdir}/${n}/weight6.mnc
897. ThresholdImage 3 ${tmpdir}/${n}/otsu.mnc ${tmpdir}/${n}/weight6.mnc 2 Inf 1 0
899. iMath 3 ${tmpdir}/${n}/weight6.mnc ME ${tmpdir}/${n}/weight6.mnc 1 1 ball 1
900. ImageMath 3 ${tmpdir}/${n}/weight6.mnc GetLargestComponent ${tmpdir}/${n}/weight6.mnc
901. iMath 3 ${tmpdir}/${n}/weight6.mnc MD ${tmpdir}/${n}/weight6.mnc 1 1 ball 1
903. #Use exclude mask if provided
904. if [[ -n ${excludemask} ]]; then
905. ImageMath 3 ${tmpdir}/${n}/weight6.mnc m ${tmpdir}/${n}/weight6.mnc ${excludemask}
906. fi
908. N4BiasFieldCorrection -d 3 -i ${tmpdir}/${n}/t1.mnc -b [ 200 ] -c [ 300x300x300x300,1e-4 ] \
909. -w ${tmpdir}/${n}/weight6.mnc -o [ ${tmpdir}/${n}/t1.mnc,${tmpdir}/${n}/bias.mnc ] -s 2 --verbose \
910. --histogram-sharpening [ 0.05,0.01,200 ] -r 0
912. ImageMath 3 ${tmpdir}/${n}/bias.mnc / ${tmpdir}/${n}/bias.mnc $(mincstats -quiet -mean ${tmpdir}/${n}/bias.mnc)
913. ImageMath 3 ${tmpdir}/${n}/prebias.mnc m ${tmpdir}/${n}/prebias.mnc ${tmpdir}/${n}/bias.mnc
914. ImageMath 3 ${tmpdir}/${n}/prebias.mnc / ${tmpdir}/${n}/prebias.mnc $(mincstats -quiet -mean ${tmpdir}/${n}/prebias.mnc)
915. ImageMath 3 ${tmpdir}/${n}/t1.mnc / ${input} ${tmpdir}/${n}/prebias.mnc
917. renorm ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/weight6.mnc usemask
918. cp ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/corrected.mnc
920. #Resample headmask into subject space, zero background and recrop
921. ImageMath 3 ${input} PadImage ${input} 50
922. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/headmask.mnc -o ${tmpdir}/headmask.mnc -r ${input} -n GenericLabel
923. ImageMath 3 ${input} m ${input} ${tmpdir}/headmask.mnc
924. ExtractRegionFromImageByMask 3 ${input} ${tmpdir}/input.crop.mnc ${tmpdir}/headmask.mnc 1 10
925. mv -f ${tmpdir}/input.crop.mnc ${input}
926. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/headmask.mnc -o ${tmpdir}/headmask.mnc -r ${input} -n GenericLabel
927. minccalc -clobber -quiet ${N4_VERBOSE:+-verbose} -unsigned -byte -expression 'A[0]>1.01?1:0' ${input} ${tmpdir}/nonzero.mnc
929. #Backup the original bias field estimate, in case cropping is not done, so we can correct the neck tissues
930. cp -f ${tmpdir}/${n}/prebias.mnc ${tmpdir}/bias_orig.mnc
932. #Need to fill the bias field with 1's in case we're padding the image
933. mincresample -clobber -quiet ${N4_VERBOSE:+-verbose} -fill -fillvalue 1 -like ${input} ${tmpdir}/${n}/prebias.mnc ${tmpdir}/${n}/bias_resample.mnc
934. cp -f ${tmpdir}/${n}/prebias.mnc ${tmpdir}/prebias.mnc
935. mincresample -clobber -quiet ${N4_VERBOSE:+-verbose} -fill -fillvalue 1 -like ${input} ${tmpdir}/prebias.mnc ${tmpdir}/prebias_resample.mnc
936. mv -f ${tmpdir}/${n}/bias_resample.mnc ${tmpdir}/${n}/bias.mnc
937. mv -f ${tmpdir}/prebias_resample.mnc ${tmpdir}/prebias.mnc
939. #Resample the exlude mask into the new recropped space
940. if [[ -n ${_arg_exclude} ]]; then
941. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${excludemask} -r ${input} -n GenericLabel -o ${excludemask}
942. fi
944. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/${n}/corrected.mnc -o ${tmpdir}/${n}/corrected.mnc -r ${input}
945. ImageMath 3 ${tmpdir}/${n}/corrected.mnc m ${tmpdir}/${n}/corrected.mnc ${tmpdir}/headmask.mnc
947. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -unsigned -byte -expression '1' ${input} ${tmpdir}/initmask.mnc
949. ################################################################################
950. #Round 1, N4 with estimate weight mask using affine registered GM/WM/CSF priors
951. ################################################################################
952. ((++n))
953. mkdir -p ${tmpdir}/${n}
955. minc_anlm ${N4_VERBOSE:+--verbose} --mt ${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS} ${tmpdir}/$((n - 1))/corrected.mnc ${tmpdir}/${n}/t1.mnc
956. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/masks/tissuemask.mnc Otsu 4 ${tmpdir}/headmask.mnc
957. ThresholdImage 3 ${tmpdir}/masks/tissuemask.mnc ${tmpdir}/masks/tissuemask.mnc 2 Inf 1 0
958. itk_vesselness --clobber --scales 8 --rescale ${tmpdir}/${n}/t1.mnc ${tmpdir}/vessels.mnc
961. antsRegistration ${N4_VERBOSE:+--verbose} -d 3 --float 1 --minc \
962. --output [ ${tmpdir}/${n}/mni ] \
963. --use-histogram-matching 1 \
964. --initial-moving-transform ${tmpdir}/$((n - 1))/mni0_GenericAffine.xfm \
965. --transform Similarity[ 0.1 ] \
966. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,32,None ] \
967. --convergence [ 500x500x500x500x450x150,1e-6,10 ] \
968. --shrink-factors 6x6x5x4x3x2 \
969. --smoothing-sigmas 2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037mm \
970. --masks [ ${REGISTRATIONBRAINMASK},NOMASK ] \
971. --transform Affine[ 0.1 ] \
972. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,64,None ] \
973. --convergence [ 500x450x150x50,1e-6,10 ] \
974. --shrink-factors 4x3x2x1 \
975. --smoothing-sigmas 1.64470459404x1.20112240879x0.735534255037x0.0mm \
976. --masks [ ${REGISTRATIONBRAINMASK},NOMASK ]
978. unset reg_initalization
980. #Make MNI-space copy of brain for BeAST
981. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/${n}/t1.mnc \
982. ${MNI_XFM:+-t ${MNI_XFM}} -t ${tmpdir}/${n}/mni0_GenericAffine.xfm -n BSpline[ 5 ] -o ${tmpdir}/${n}/mni.mnc -r ${RESAMPLEMODEL}
984. #BSpline[ 5 ] does weird things to intensity, clip back to positive range
985. mincmath -quiet ${N4_VERBOSE:+-verbose} -clamp -const2 0 65535 ${tmpdir}/${n}/mni.mnc ${tmpdir}/${n}/mni.clamp.mnc
986. mv -f ${tmpdir}/${n}/mni.clamp.mnc ${tmpdir}/${n}/mni.mnc
988. #Shrink the MNI mask for the first intensity matching
989. iMath 3 ${tmpdir}/${n}/shrinkmask.mnc ME ${RESAMPLEMODELBRAINMASK} 2 1 ball 1
991. #Intensity normalize
992. volume_pol ${N4_VERBOSE:+--verbose} --order 1 --min 0 --max 100 --noclamp \
993. --source_mask ${tmpdir}/${n}/shrinkmask.mnc --target_mask ${RESAMPLEMODELBRAINMASK} \
994. ${tmpdir}/${n}/mni.mnc ${RESAMPLEMODEL} ${tmpdir}/${n}/mni.norm.mnc
996. #Run a quick beast to get a brain mask
997. mincbeast ${N4_VERBOSE:+-verbose} -sparse -v2 -double -fill -median -same_res -flip -conf ${BEAST_CONFIG} \
998. ${BEASTLIBRARY_DIR} ${tmpdir}/${n}/mni.norm.mnc ${tmpdir}/${n}/beastmask.mnc
1000. #Resample beast mask and MNI mask to native space
1001. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -r ${tmpdir}/${n}/t1.mnc \
1002. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] ${MNI_XFM:+-t [${MNI_XFM},1]} \
1003. -i ${tmpdir}/${n}/beastmask.mnc -o ${tmpdir}/${n}/bmask.mnc -n GenericLabel
1004. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -r ${tmpdir}/${n}/t1.mnc \
1005. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -i ${REGISTRATIONBRAINMASK} \
1006. -o ${tmpdir}/${n}/mnimask.mnc -n GenericLabel
1008. #BeAST Failure mode of a chunk of almost unattached voxels, try to remove
1009. iMath 3 ${tmpdir}/${n}/bmask.mnc ME ${tmpdir}/${n}/bmask.mnc 1 1 ball 1
1010. ImageMath 3 ${tmpdir}/${n}/bmask.mnc GetLargestComponent ${tmpdir}/${n}/bmask.mnc
1011. iMath 3 ${tmpdir}/${n}/bmask.mnc MD ${tmpdir}/${n}/bmask.mnc 1 1 ball 1
1013. cp -f ${tmpdir}/${n}/bmask.mnc ${tmpdir}/masks/bmask.mnc
1014. cp -f ${tmpdir}/${n}/mnimask.mnc ${tmpdir}/masks/affinemask.mnc
1016. ImageMath 3 ${tmpdir}/${n}/mask.mnc addtozero ${tmpdir}/masks/bmask.mnc ${tmpdir}/masks/affinemask.mnc
1017. ImageMath 3 ${tmpdir}/${n}/mask.mnc GetLargestComponent ${tmpdir}/${n}/mask.mnc
1019. iMath 3 ${tmpdir}/${n}/mask_D.mnc MD ${tmpdir}/${n}/mask.mnc 2 1 ball 1
1021. #Resample MNI Priors to Native space for classification
1022. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${WMPRIOR} \
1023. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/SegmentationPrior3.mnc -n Linear
1024. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${GMPRIOR} \
1025. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/SegmentationPrior2.mnc -n Linear
1026. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${CSFPRIOR} \
1027. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/SegmentationPrior1.mnc -n Linear
1029. if [[ -n ${excludemask} ]]; then
1030. ImageMath 3 ${tmpdir}/${n}/mask_D.mnc m ${tmpdir}/${n}/mask_D.mnc ${excludemask}
1031. fi
1033. #Estimate outlier to exclude from classification
1034. outlier_mask ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/bmask.mnc ${tmpdir}/${n}/hotmask.mnc
1035. ImageMath 3 ${tmpdir}/${n}/mask_D.mnc m ${tmpdir}/${n}/mask_D.mnc ${tmpdir}/${n}/hotmask.mnc
1037. #Classify brain
1038. Atropos ${N4_VERBOSE:+--verbose} -d 3 -x ${tmpdir}/${n}/mask_D.mnc -c [ 5,0.005 ] -a ${tmpdir}/${n}/t1.mnc -s 1x2 -s 2x3 \
1039. -i PriorProbabilityImages[ 3,${tmpdir}/${n}/SegmentationPrior%d.mnc,0.1 ] -k Gaussian -m [ 0.1,1x1x1 ] \
1040. -o ${tmpdir}/${n}/classify.mnc -r 1 -p Socrates[ 0 ] --winsorize-outliers BoxPlot
1042. #Convert classification to a brain mask and brain tissue mask
1043. ThresholdImage 3 ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/2.mnc 2 2 1 0
1044. ThresholdImage 3 ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/3.mnc 3 3 1 0
1045. ImageMath 3 ${tmpdir}/${n}/2.mnc GetLargestComponent ${tmpdir}/${n}/2.mnc
1046. ImageMath 3 ${tmpdir}/${n}/3.mnc GetLargestComponent ${tmpdir}/${n}/3.mnc
1048. ImageMath 3 ${tmpdir}/${n}/weight.mnc addtozero ${tmpdir}/${n}/2.mnc ${tmpdir}/${n}/3.mnc
1049. iMath 3 ${tmpdir}/${n}/weight.mnc ME ${tmpdir}/${n}/weight.mnc 1 1 ball 1
1050. ImageMath 3 ${tmpdir}/${n}/weight.mnc GetLargestComponent ${tmpdir}/${n}/weight.mnc
1051. iMath 3 ${tmpdir}/${n}/weight.mnc MD ${tmpdir}/${n}/weight.mnc 2 1 ball 1
1053. iMath 3 ${tmpdir}/${n}/mask2.mnc MC ${tmpdir}/${n}/weight.mnc 5 1 ball 1
1054. ImageMath 3 ${tmpdir}/${n}/mask2.mnc FillHoles ${tmpdir}/${n}/mask2.mnc 2
1055. cp -f ${tmpdir}/${n}/mask2.mnc ${tmpdir}/masks/classifymask${n}.mnc
1057. #User provided exclusion mask
1058. if [[ -n ${excludemask} ]]; then
1059. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${excludemask}
1060. fi
1062. #Remove outliers round 2
1063. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/hotmask.mnc
1064. ImageMath 3 ${tmpdir}/${n}/weight.mnc GetLargestComponent ${tmpdir}/${n}/weight.mnc
1066. ImageMath 3 ${tmpdir}/${n}/classify.mnc m ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/mask2.mnc
1068. #Always exclude 0 from correction
1069. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -unsigned -byte -expression 'A[0]>1.01?1:0' ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/nonzero.mnc
1070. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/nonzero.mnc
1071. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/nonzero.mnc
1073. do_N4_correct ${input} ${tmpdir}/initmask.mnc ${tmpdir}/${n}/mask2.mnc ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/corrected.mnc ${tmpdir}/${n}/bias.mnc 2 ${tmpdir}/${n}/classify.mnc
1075. #Calculate coeffcient of variation between this round bias field and prior round
1076. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -zero -expression 'A[0]/A[1]' ${tmpdir}/$((n - 1))/bias.mnc ${tmpdir}/${n}/bias.mnc ${tmpdir}/${n}/ratio.mnc
1077. python -c "print(float(\"$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -stddev ${tmpdir}/${n}/ratio.mnc)\") / float(\"$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -mean ${tmpdir}/${n}/ratio.mnc)\"))" >>${tmpdir}/convergence.txt
1079. if [[ ${_arg_debug} == "off" ]]; then
1080. rm -rf ${tmpdir}/$((n - 1))
1081. fi
1083. ################################################################################
1084. #Round 2, N4 with classification from nonlinear registered priors
1085. ################################################################################
1086. ((++n))
1087. mkdir -p ${tmpdir}/${n}
1089. minc_anlm ${N4_VERBOSE:+--verbose} --mt ${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS} ${tmpdir}/$((n - 1))/corrected.mnc ${tmpdir}/${n}/t1.mnc
1090. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/masks/tissuemask.mnc Otsu 4 ${tmpdir}/headmask.mnc
1091. ThresholdImage 3 ${tmpdir}/masks/tissuemask.mnc ${tmpdir}/masks/tissuemask.mnc 2 Inf 1 0
1093. #Affine register to MNI space, tweak registration
1094. antsRegistration ${N4_VERBOSE:+--verbose} -d 3 --float 1 --minc \
1095. --output [ ${tmpdir}/${n}/mni ] \
1096. --use-histogram-matching 1 \
1097. --initial-moving-transform ${tmpdir}/$((n - 1))/mni0_GenericAffine.xfm \
1098. --transform Affine[ 0.05 ] \
1099. --metric Mattes[ ${REGISTRATIONMODEL},${tmpdir}/${n}/t1.mnc,1,64,None ] \
1100. --convergence [ 500x450x150x50,1e-6,10 ] \
1101. --shrink-factors 4x3x2x1 \
1102. --smoothing-sigmas 1.64470459404x1.20112240879x0.735534255037x0.0mm \
1103. --masks [ ${REGISTRATIONBRAINMASK},${tmpdir}/$((n - 1))/mask2.mnc ]
1105. cp -f ${tmpdir}/$((n - 1))/mask2.mnc ${tmpdir}/${n}/mask.mnc
1106. iMath 3 ${tmpdir}/${n}/extractmask.mnc MD ${tmpdir}/${n}/mask.mnc 1 1 ball 1
1108. ImageMath 3 ${tmpdir}/${n}/t1.extracted.mnc m ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/extractmask.mnc
1109. ImageMath 3 ${tmpdir}/extractmodel.mnc m ${REGISTRATIONMODEL} ${REGISTRATIONBRAINMASK}
1111. #Non linearly register priors
1112. #We use the extracted images because subjects with different distance between
1113. #brain and skull consistently fail
1114. antsRegistration ${N4_VERBOSE:+--verbose} -d 3 --float 1 --minc \
1115. --output [ ${tmpdir}/${n}/nonlin ] \
1116. --initial-moving-transform ${tmpdir}/${n}/mni0_GenericAffine.xfm \
1117. --use-histogram-matching 1 \
1118. --transform SyN[ 0.1,3,0 ] \
1119. --metric CC[ ${tmpdir}/extractmodel.mnc,${tmpdir}/${n}/t1.extracted.mnc,1,4 ] \
1120. --convergence [ 500x500x500x500x500x500x500x500x500x500x0x0x0,1e-6,10 ] \
1121. --shrink-factors 5x5x5x5x5x5x5x5x5x4x3x2x1 \
1122. --smoothing-sigmas 5.50423435717x5.07820577132x4.65192708599x4.22532260674x3.79828256043x3.37064139994x2.94213702015x2.51232776601x2.08040503813x1.64470459404x1.20112240879x0.735534255037x0.0mm \
1123. --masks [ NOMASK,NOMASK ] \
1124. --transform SyN[ 0.1,3,0 ] \
1125. --metric CC[ ${tmpdir}/extractmodel.mnc,${tmpdir}/${n}/t1.extracted.mnc,1,2 ] \
1126. --convergence [ 500x500x225x225x0,1e-6,10 ] \
1127. --shrink-factors 5x4x3x2x1 \
1128. --smoothing-sigmas 2.08040503813x1.64470459404x1.20112240879x0.735534255037x0.0mm \
1129. --masks [ ${REGISTRATIONBRAINMASK},${tmpdir}/${n}/extractmask.mnc ]
1131. #Save MNI space registration for QC later
1132. cp -f ${tmpdir}/${n}/mni0_GenericAffine.xfm ${tmpdir}/mni0_GenericAffine.xfm
1134. #Resample MNI Priors to Native space for classification
1135. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${WMPRIOR} \
1136. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -t ${tmpdir}/${n}/nonlin1_inverse_NL.xfm \
1137. -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/SegmentationPrior3.mnc -n Linear
1138. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${GMPRIOR} \
1139. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -t ${tmpdir}/${n}/nonlin1_inverse_NL.xfm \
1140. -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/SegmentationPrior2.mnc -n Linear
1141. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${CSFPRIOR} \
1142. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -t ${tmpdir}/${n}/nonlin1_inverse_NL.xfm \
1143. -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/SegmentationPrior1.mnc -n Linear
1145. #Resample back to subject space
1146. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${REGISTRATIONBRAINMASK} \
1147. -t [ ${tmpdir}/${n}/mni0_GenericAffine.xfm,1 ] -t ${tmpdir}/${n}/nonlin1_inverse_NL.xfm -r ${tmpdir}/${n}/t1.mnc -o ${tmpdir}/${n}/mnimask.mnc -n GenericLabel
1148. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -unsigned -byte -expression '(A[0]>=0.25||A[1]>=0.25)?1:0' \
1149. ${tmpdir}/${n}/SegmentationPrior3.mnc ${tmpdir}/${n}/SegmentationPrior2.mnc ${tmpdir}/${n}/mniprobmask.mnc
1151. #Make MNI mask a combination of both the MNI mask and the tissue probabilty mask
1152. ImageMath 3 ${tmpdir}/${n}/mnimask.mnc addtozero ${tmpdir}/${n}/mnimask.mnc ${tmpdir}/${n}/mniprobmask.mnc
1154. #Last time we generate MNI mask, save it outside iterations
1155. cp -f ${tmpdir}/${n}/mnimask.mnc ${tmpdir}/masks/mnimask.mnc
1157. #Vote a consensus mask from prior masking estimates
1158. ImageMath 3 ${tmpdir}/${n}/mask.mnc MajorityVoting ${tmpdir}/masks/*mnc
1159. iMath 3 ${tmpdir}/${n}/mask.mnc MC ${tmpdir}/${n}/mask.mnc 1 1 ball 1
1161. #Expand the mask a bit
1162. iMath 3 ${tmpdir}/${n}/mask_D.mnc MD ${tmpdir}/${n}/mask.mnc 1 1 ball 1
1164. if [[ -n ${excludemask} ]]; then
1165. ImageMath 3 ${tmpdir}/${n}/mask_D.mnc m ${tmpdir}/${n}/mask_D.mnc ${excludemask}
1166. fi
1168. #Find outliers to exclude from classification
1169. ThresholdImage 3 ${tmpdir}/$((n - 1))/classify.mnc ${tmpdir}/${n}/outlier_wm.mnc 3 3 1 0
1170. ImageMath 3 ${tmpdir}/${n}/outlier_wm.mnc GetLargestComponent ${tmpdir}/${n}/outlier_wm.mnc
1171. outlier_mask ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/outlier_wm.mnc ${tmpdir}/${n}/hotmask.mnc
1172. ImageMath 3 ${tmpdir}/${n}/mask_D.mnc m ${tmpdir}/${n}/mask_D.mnc ${tmpdir}/${n}/hotmask.mnc
1174. #Do an initial classification using the MNI priors
1175. Atropos ${N4_VERBOSE:+--verbose} -d 3 -x ${tmpdir}/${n}/mask_D.mnc -c [ 5,0.005 ] -a ${tmpdir}/${n}/t1.mnc -s 1x2 -s 2x3 \
1176. -i PriorProbabilityImages[ 3,${tmpdir}/${n}/SegmentationPrior%d.mnc,${_arg_classification_prior_weight} ] -k Gaussian -m [ 0.1,1x1x1 ] \
1177. -o [ ${tmpdir}/${n}/classify.mnc,${tmpdir}/${n}/SegmentationPosteriors%d.mnc ] -r 1 -p Aristotle[ 0 ] --winsorize-outliers BoxPlot \
1178. -l [ 0.69314718055994530942,1 ]
1180. #Convert classification to the mask
1181. classify_to_mask
1183. cp ${tmpdir}/${n}/classifymask.mnc ${tmpdir}/masks/classifymask${n}.mnc
1185. ImageMath 3 ${tmpdir}/${n}/mask2.mnc MajorityVoting ${tmpdir}/masks/*mnc
1186. iMath 3 ${tmpdir}/${n}/mask2.mnc MC ${tmpdir}/${n}/mask2.mnc 1 1 ball 1
1188. #Combine GM and WM proabability images into a N4 mask,
1189. ImageMath 3 ${tmpdir}/${n}/weight.mnc PureTissueN4WeightMask ${tmpdir}/${n}/SegmentationPosteriors2.mnc ${tmpdir}/${n}/SegmentationPosteriors3.mnc
1190. ImageMath 3 ${tmpdir}/${n}/weight.mnc RescaleImage ${tmpdir}/${n}/weight.mnc 0 1
1192. ImageMath 3 ${tmpdir}/${n}/weightmask.mnc GetLargestComponent ${tmpdir}/${n}/weight.mnc
1193. iMath 3 ${tmpdir}/${n}/weightmask.mnc ME ${tmpdir}/${n}/weightmask.mnc 1 1 ball 1
1194. ImageMath 3 ${tmpdir}/${n}/weightmask.mnc GetLargestComponent ${tmpdir}/${n}/weightmask.mnc
1195. iMath 3 ${tmpdir}/${n}/weightmask.mnc MD ${tmpdir}/${n}/weightmask.mnc 1 1 ball 1
1196. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/weightmask.mnc
1198. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/hotmask.mnc
1200. #Clip the classification weight and posteriors
1201. for item in ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/SegmentationPosteriors1.mnc ${tmpdir}/${n}/SegmentationPosteriors2.mnc ${tmpdir}/${n}/SegmentationPosteriors3.mnc; do
1202. ImageMath 3 ${item} m ${item} ${tmpdir}/${n}/mask2.mnc
1203. done
1205. if [[ -n ${excludemask} ]]; then
1206. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${excludemask}
1207. fi
1209. #Always exclude 0 from correction
1210. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -unsigned -byte -expression 'A[0]>1.01?1:0' ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/nonzero.mnc
1211. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/nonzero.mnc
1212. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/nonzero.mnc
1214. do_N4_correct ${input} ${tmpdir}/initmask.mnc ${tmpdir}/${n}/mask2.mnc ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/corrected.mnc ${tmpdir}/${n}/bias.mnc 2 ${tmpdir}/${n}/classify.mnc
1216. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -zero -expression 'A[0]/A[1]' ${tmpdir}/$((n - 1))/bias.mnc ${tmpdir}/${n}/bias.mnc ${tmpdir}/${n}/ratio.mnc
1217. python -c "print(float(\"$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -stddev ${tmpdir}/${n}/ratio.mnc)\") / float(\"$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -mean ${tmpdir}/${n}/ratio.mnc)\"))" >>${tmpdir}/convergence.txt
1219. if [[ ${_arg_debug} == "off" ]]; then
1220. rm -rf ${tmpdir}/$((n - 1))
1221. fi
1223. ################################################################################
1224. #Remaining rounds, N4 with segmentation posteriors bootstrapped from prior run until convergence
1225. ################################################################################
1226. while true; do
1227. ((++n))
1228. mkdir -p ${tmpdir}/${n}
1230. minc_anlm ${N4_VERBOSE:+--verbose} --mt ${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS} ${tmpdir}/$((n - 1))/corrected.mnc ${tmpdir}/${n}/t1.mnc
1231. ThresholdImage 3 ${tmpdir}/${n}/t1.mnc ${tmpdir}/masks/tissuemask.mnc Otsu 4 ${tmpdir}/headmask.mnc
1232. ThresholdImage 3 ${tmpdir}/masks/tissuemask.mnc ${tmpdir}/masks/tissuemask.mnc 2 Inf 1 0
1234. cp -f ${tmpdir}/$((n - 1))/mask2.mnc ${tmpdir}/${n}/mask.mnc
1235. iMath 3 ${tmpdir}/${n}/mask_D.mnc MD ${tmpdir}/${n}/mask.mnc 1 1 ball 1
1237. if [[ -n ${excludemask} ]]; then
1238. ImageMath 3 ${tmpdir}/${n}/mask_D.mnc m ${tmpdir}/${n}/mask_D.mnc ${excludemask}
1239. fi
1241. ThresholdImage 3 ${tmpdir}/$((n - 1))/classify.mnc ${tmpdir}/${n}/outlier_wm.mnc 3 3 1 0
1242. ImageMath 3 ${tmpdir}/${n}/outlier_wm.mnc GetLargestComponent ${tmpdir}/${n}/outlier_wm.mnc
1243. outlier_mask ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/outlier_wm.mnc ${tmpdir}/${n}/hotmask.mnc
1244. ImageMath 3 ${tmpdir}/${n}/mask_D.mnc m ${tmpdir}/${n}/mask_D.mnc ${tmpdir}/${n}/hotmask.mnc
1246. #Do a classification using the last round posteriors, remove outliers
1247. Atropos ${N4_VERBOSE:+--verbose} -d 3 -x ${tmpdir}/${n}/mask_D.mnc -c [ 5,0.005 ] -a ${tmpdir}/${n}/t1.mnc -s 1x2 -s 2x3 \
1248. -i PriorProbabilityImages[ 3,${tmpdir}/$((n - 1))/SegmentationPosteriors%d.mnc,0.5 ] -k Gaussian -m [ 0.1,1x1x1 ] \
1249. -o [ ${tmpdir}/${n}/classify.mnc,${tmpdir}/${n}/SegmentationPosteriors%d.mnc ] -r 1 -p Aristotle[ 1 ] --winsorize-outliers BoxPlot \
1250. -l [ 0.69314718055994530942,1 ]
1252. classify_to_mask
1254. cp -f ${tmpdir}/${n}/classifymask.mnc ${tmpdir}/masks/classifymask${n}.mnc
1256. ImageMath 3 ${tmpdir}/${n}/mask2.mnc MajorityVoting ${tmpdir}/masks/*mnc
1257. iMath 3 ${tmpdir}/${n}/mask2.mnc MC ${tmpdir}/${n}/mask2.mnc 1 1 ball 1
1259. #Combine GM and WM probably images into a N4 mask,
1260. ImageMath 3 ${tmpdir}/${n}/weight.mnc PureTissueN4WeightMask ${tmpdir}/${n}/SegmentationPosteriors2.mnc ${tmpdir}/${n}/SegmentationPosteriors3.mnc
1261. ImageMath 3 ${tmpdir}/${n}/weight.mnc RescaleImage ${tmpdir}/${n}/weight.mnc 0 1
1263. ImageMath 3 ${tmpdir}/${n}/weightmask.mnc GetLargestComponent ${tmpdir}/${n}/weight.mnc
1264. iMath 3 ${tmpdir}/${n}/weightmask.mnc ME ${tmpdir}/${n}/weightmask.mnc 1 1 ball 1
1265. ImageMath 3 ${tmpdir}/${n}/weightmask.mnc GetLargestComponent ${tmpdir}/${n}/weightmask.mnc
1266. iMath 3 ${tmpdir}/${n}/weightmask.mnc MD ${tmpdir}/${n}/weightmask.mnc 1 1 ball 1
1267. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/weightmask.mnc
1269. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/hotmask.mnc
1271. #Clip the classification weight and posteriors
1272. for item in ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/classify.mnc ${tmpdir}/${n}/SegmentationPosteriors1.mnc ${tmpdir}/${n}/SegmentationPosteriors2.mnc ${tmpdir}/${n}/SegmentationPosteriors3.mnc; do
1273. ImageMath 3 ${item} m ${item} ${tmpdir}/${n}/mask2.mnc
1274. done
1276. if [[ -n ${excludemask} ]]; then
1277. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${excludemask}
1278. fi
1280. #Always exclude 0 from correction
1281. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -unsigned -byte -expression 'A[0]>1.01?1:0' ${tmpdir}/${n}/t1.mnc ${tmpdir}/${n}/nonzero.mnc
1282. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/nonzero.mnc
1283. ImageMath 3 ${tmpdir}/${n}/weight.mnc m ${tmpdir}/${n}/weight.mnc ${tmpdir}/nonzero.mnc
1285. do_N4_correct ${input} ${tmpdir}/initmask.mnc ${tmpdir}/${n}/mask2.mnc ${tmpdir}/${n}/weight.mnc ${tmpdir}/${n}/corrected.mnc ${tmpdir}/${n}/bias.mnc 2 ${tmpdir}/${n}/classify.mnc
1287. #Compute coeffcient of variation
1288. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -zero -expression 'A[0]/A[1]' ${tmpdir}/$((n - 1))/bias.mnc ${tmpdir}/${n}/bias.mnc ${tmpdir}/${n}/ratio.mnc
1289. python -c "print(float(\"$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -stddev ${tmpdir}/${n}/ratio.mnc)\") / float(\"$(mincstats -quiet -mask ${tmpdir}/${n}/mask2.mnc -mask_binvalue 1 -mean ${tmpdir}/${n}/ratio.mnc)\"))" >>${tmpdir}/convergence.txt
1291. if [[ ${_arg_debug} == "off" ]]; then
1292. rm -rf ${tmpdir}/$((n - 1))
1293. fi
1295. # Break if greater than max iterations or less than convergence threshold
1296. [[ (${n} -lt ${_arg_max_iterations}) && ($(python -c "print($(tail -1 ${tmpdir}/convergence.txt) > ${_arg_convergence_threshold})") == "True") ]] || break
1298. done
1300. echo "--------------------"
1301. echo "Convergence results:"
1302. cat ${tmpdir}/convergence.txt
1303. echo "--------------------"
1305. #If cropping is enabled, recrop the originput file and resample the mask again
1306. if [[ ${_arg_autocrop} == "on" ]]; then
1307. ImageMath 3 ${originput} PadImage ${originput} 50
1308. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/headmask.mnc -o ${tmpdir}/finalheadmask.mnc -r ${originput} -n GenericLabel
1309. ImageMath 3 ${originput} m ${originput} ${tmpdir}/finalheadmask.mnc
1310. ExtractRegionFromImageByMask 3 ${originput} ${tmpdir}/originput.crop.mnc ${tmpdir}/finalheadmask.mnc 1 10
1311. mv -f ${tmpdir}/originput.crop.mnc ${originput}
1312. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -r ${originput} -i ${tmpdir}/headmask.mnc \
1313. -o ${tmpdir}/headmask.mnc -n GenericLabel
1314. cp -f ${tmpdir}/headmask.mnc ${tmpdir}/fgmask.mnc
1315. else
1316. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -r ${originput} -i ${tmpdir}/fgmask_orig.mnc \
1317. -o ${tmpdir}/fgmask.mnc -n GenericLabel
1318. ImageMath 3 ${originput} m ${originput} ${tmpdir}/fgmask.mnc
1319. fi
1321. #Resample final results into original space and correct original input file
1322. n4input=${originput}
1323. n4corrected=${tmpdir}/corrected.mnc
1324. n4classifymask=${tmpdir}/finalclassify.mnc
1326. #Reconstruct a final bias field
1327. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -r ${originput} \
1328. -n BSpline[5] -i ${tmpdir}/${n}/bias.mnc -o ${tmpdir}/finalbias.mnc
1330. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -r ${originput} \
1331. -n BSpline[5] -i ${tmpdir}/bias_orig.mnc -o ${tmpdir}/bias_orig.mnc
1333. ImageMath 3 ${tmpdir}/finalbias.mnc addtozero ${tmpdir}/finalbias.mnc ${tmpdir}/bias_orig.mnc
1334. ImageMath 3 ${tmpdir}/finalbias.mnc addtozero ${tmpdir}/finalbias.mnc 1
1336. ImageMath 3 ${n4corrected} / ${originput} ${tmpdir}/finalbias.mnc
1338. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/${n}/mask2.mnc -o ${tmpdir}/finalmask.mnc -r ${n4corrected} -n GenericLabel
1339. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/masks/bmask.mnc -o ${tmpdir}/finalbmask.mnc -r ${n4corrected} -n GenericLabel
1340. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/${n}/classifymask.mnc -o ${tmpdir}/finalclassifymask.mnc -r ${n4corrected} -n GenericLabel
1341. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/masks/mnimask.mnc -o ${tmpdir}/finalmnimask.mnc -r ${n4corrected} -n GenericLabel
1342. antsApplyTransforms ${N4_VERBOSE:+--verbose} -d 3 -i ${tmpdir}/${n}/classify.mnc -o ${tmpdir}/finalclassify.mnc -r ${n4corrected} -n GenericLabel
1344. valuelow=$(mincstats -quiet -mask ${tmpdir}/fgmask.mnc -mask_binvalue 1 -pctT 0.1 ${n4corrected})
1345. valuewm=$(mincstats -quiet -median -mask ${n4classifymask} -mask_binvalue 3 ${n4corrected})
1346. valuegm=$(mincstats -quiet -median -mask ${n4classifymask} -mask_binvalue 2 ${n4corrected})
1347. valuehigh=$(mincstats -quiet -mask ${tmpdir}/fgmask.mnc -mask_binvalue 1 -pctT 99.9 ${n4corrected})
1349. mapping=($(python -c "import numpy as np; print(np.array2string(np.linalg.solve(np.array([[1, ${valuelow}, ${valuelow}**2], [1, ((${valuewm}+${valuegm})/2.0), ((${valuewm}+${valuegm})/2.0)**2], [1, ${valuehigh}, ${valuehigh}**2]]),np.array([0,32767,65535])),separator= ' ')[1:-1])"))
1351. minccalc -quiet ${N4_VERBOSE:+-verbose} -short -unsigned -expression "clamp(A[0]^2*${mapping[2]} + A[0]*${mapping[1]} + ${mapping[0]},0,65535)" \
1352. ${n4corrected} $(dirname ${n4corrected})/$(basename ${n4corrected} .mnc).norm.mnc
1354. mv -f $(dirname ${n4corrected})/$(basename ${n4corrected} .mnc).norm.mnc ${n4corrected}
1356. cp -f ${tmpdir}/corrected.mnc ${output}
1358. #Output final classification files if standalone
1359. if [[ ${_arg_standalone} == "on" || ${_arg_debug} == "on" ]]; then
1360. make_qc
1361. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -byte -unsigned ${tmpdir}/finalbmask.mnc $(dirname ${output})/$(basename ${output} .mnc).beastmask.mnc
1362. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -byte -unsigned ${tmpdir}/finalmnimask.mnc $(dirname ${output})/$(basename ${output} .mnc).mnimask.mnc
1363. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -byte -unsigned ${tmpdir}/finalclassify.mnc $(dirname $output)/$(basename ${output} .mnc).classify.mnc
1364. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -byte -unsigned ${tmpdir}/finalmask.mnc $(dirname $output)/$(basename ${output} .mnc).mask.mnc
1365. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -byte -unsigned ${tmpdir}/finalclassifymask.mnc $(dirname $output)/$(basename ${output} .mnc).classifymask.mnc
1366. minccalc -quiet ${N4_VERBOSE:+-verbose} -clobber -short -unsigned -expression 'A[0]*A[1]' ${output} ${tmpdir}/finalmask.mnc ${tmpdir}/output.extracted.mnc
1367. ExtractRegionFromImageByMask 3 ${tmpdir}/output.extracted.mnc ${tmpdir}/output.extracted.crop.mnc ${tmpdir}/finalmask.mnc 1 10
1368. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -short -unsigned ${tmpdir}/output.extracted.crop.mnc $(dirname $output)/$(basename ${output} .mnc).extracted.mnc
1370. minc_anlm ${N4_VERBOSE:+--verbose} --mt ${ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS} ${tmpdir}/corrected.mnc ${tmpdir}/corrected.denoise.mnc
1371. mincreshape -quiet ${N4_VERBOSE:+-verbose} -clobber -short -unsigned ${tmpdir}/corrected.denoise.mnc $(dirname $output)/$(basename ${output} .mnc).denoise.mnc
1372. fi
1374. if [[ ${_arg_debug} == "off" ]]; then
1375. rm -rf ${tmpdir}
1376. fi
1378. # ] <-- needed because of Argbash

iterativeN4_multispectral.sh at commit 70dbcef, under other · at the source