Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Template initialization to improve convergence #1227

Merged
merged 5 commits into from
Aug 24, 2021
Merged

Template initialization to improve convergence #1227

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
bc1ca3f
ENH: Better template initialization from population averages
cookpa Aug 18, 2021
6900651
ENH: Check all expected warp images get created
cookpa Aug 18, 2021
4f2d4e3
DOC: Explain new features, and some older ones
cookpa Aug 18, 2021
f3390d4
ENH: Also check smoothing levels consistent with iteration levels
cookpa Aug 19, 2021
14853e4
ENH: Allow user to control averaging method in initial stages
cookpa Aug 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Examples/AverageImages.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ int AverageImages1(unsigned int argc, char *argv[])
reader->SetFileName(argv[bigimage]);
reader->Update();
typename ImageType::Pointer averageimage = reader->GetOutput();
std::cout << " Setting physcal space of output average image based on largest image " << std::endl;
std::cout << " Setting physical space of output average image based on largest image " << std::endl;
unsigned int vectorlength = reader->GetImageIO()->GetNumberOfComponents();
std::cout << " Averaging " << numberofimages << " images with dim = " << ImageDimension << " vector components "
<< vectorlength << std::endl;
Expand Down
85 changes: 70 additions & 15 deletions Scripts/antsMultivariateTemplateConstruction.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ Optional arguments:
-A sharpening applied to template at each iteration (default 1)
0 = none
1 = Laplacian
2 = Unsharp mask
2 = Unsharp mask

-c: Control for parallel computation (default 1) -- 0 == run serially, 1 == SGE qsub,
2 == use PEXEC (localhost), 3 == Apple XGrid, 4 == PBS qsub, 5 == SLURM
Expand Down Expand Up @@ -128,7 +128,15 @@ Optional arguments:

-s: Type of similarity metric used for registration.

-t: Type of transformation model used for registration.
-t: Type of transformation model used for registration. Supported options are case sensitive.
GR : Greedy SyN (default for scalar data)
GR_Constrained : Greedy SyN with regularization on the total deformation (default for time series)
EL : Elastic
EX : Exponential
DD : Greedy exponential, diffemorphic-demons-style optimization
SY : LDDMM-style SyN with symmetric time-dependent gradient estimation
LDDMM : Like SY, but with asymmetric time-dependent gradient estimation
S2 : Like SY, but with no time-dependent gradient estimation

-x: XGrid arguments (e.g., -x "-p password -h controlhost")

Expand Down Expand Up @@ -237,9 +245,15 @@ function shapeupdatetotemplate() {
echo " shapeupdatetotemplate---voxel-wise averaging of the warped images to the current template"
echo "--------------------------------------------------------------------------------------"

imagelist=(`ls ${outputname}template${whichtemplate}*WarpedToTemplate.nii.gz`)
if [[ ${#imagelist[@]} -ne ${IMAGESPERMODALITY} ]] ; then
echo "ERROR shapeupdatedtotemplate - imagelist length is ${#imagelist[@]}, expected ${IMAGESPERMODALITY}"
exit 1
fi

case $summarizemethod in
0) #mean
${ANTSPATH}/AverageImages $dim ${template} 0 ${templatename}${whichtemplate}*WarpedToTemplate.nii.gz
${ANTSPATH}/AverageImages $dim ${template} 0 ${templatename}${whichtemplate}*WarpedToTemplate.nii.gz
;;
1) #mean of normalized images
${ANTSPATH}/AverageImages $dim ${template} 2 ${templatename}${whichtemplate}*WarpedToTemplate.nii.gz
Expand All @@ -259,20 +273,20 @@ function shapeupdatetotemplate() {
echo " shapeupdatetotemplate---sharpening of the new template"
echo "--------------------------------------------------------------------------------------"

case $sharpenmethod in
case $sharpenmethod in
0)
echo "Sharpening method none"
;;
1)
echo "Laplacian sharpening"
${ANTSPATH}/ImageMath $dim $template Sharpen $template
${ANTSPATH}/ImageMath $dim $template Sharpen $template
;;
2)
echo "Unsharp mask sharpening"
${ANTSPATH}/ImageMath $dim $template UnsharpMask $template 0.5 1 0 0
;;
esac

if [[ $whichtemplate -eq 0 ]] ;
then
echo
Expand Down Expand Up @@ -452,7 +466,7 @@ while getopts "A:a:b:c:d:g:h:i:j:k:m:n:o:p:s:r:t:w:x:y:z:" OPT
exit 0
;;
A) # Sharpening method
SHARPENMETHOD=$OPTARG
SHARPENMETHOD=$OPTARG
;;
a) # summarizing statistic
STATSMETHOD=$OPTARG
Expand Down Expand Up @@ -540,13 +554,18 @@ elif [[ $nargs -lt 6 ]]
Usage >&2
fi

OUTPUT_DIR=${OUTPUTNAME%\/*}
OUTPUT_DIR=`dirname ${OUTPUTNAME}`
if [[ ! -d $OUTPUT_DIR ]];
then
echo "The output directory \"$OUTPUT_DIR\" does not exist. Making it."
mkdir -p $OUTPUT_DIR
fi

# Intermediate template output. Keep the template for each iteration and also the average warp if defined.
# Useful for debugging and monitoring convergence
intermediateTemplateDir=${OUTPUT_DIR}/intermediateTemplates
mkdir -p $intermediateTemplateDir

if [[ $DOQSUB -eq 1 || $DOQSUB -eq 4 ]];
then
qq=`which qsub`
Expand Down Expand Up @@ -811,14 +830,17 @@ if [[ $NUMBEROFMODALITIES -gt 1 ]];
echo "--------------------------------------------------------------------------------------"
fi

# Useful to check the right number of images exist for various ops
IMAGESPERMODALITY=$(( ${#IMAGESETARRAY[@]} / ${NUMBEROFMODALITIES} ))

# check for initial template images
for (( i = 0; i < $NUMBEROFMODALITIES; i++ ))
do
setCurrentImageSet $i

if [[ -n "${REGTEMPLATES[$i]}" ]];
then
if [[ ! -r "${REGTEMPLATES[$i]}" ]];
if [[ ! -r "${REGTEMPLATES[$i]}" ]];
then
echo "Initial template {REGTEMPLATES[$i]} cannot be read"
exit 1
Expand All @@ -837,13 +859,34 @@ for (( i = 0; i < $NUMBEROFMODALITIES; i++ ))
echo "--------------------------------------------------------------------------------------"
# Normalize but don't sharpen at this stage
${ANTSPATH}/AverageImages $DIM ${TEMPLATES[$i]} 2 ${CURRENTIMAGESET[@]}
# Quickly align COM of input images to average, and then recompute average
IMAGECOMSET=()
for (( j = 0; j < ${#CURRENTIMAGESET[@]}; j+=1 ))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One other implementation note here, I found for best performance I often had to do two rounds of this.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you have an example of this? I can't see how it would help except to shift the location of the image set inside the bounding box. Which may be desirable (or not if it shifts the wrong way)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

After checking things over, my tests show I don't actually need two rounds, however, padding/recropping the template is important as otherwise you can get a non-centered FOV.

Compare, left no pad/recrop vs right:
https://imgur.com/a/ded9IPL

Pad/recrop code:
https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction/blob/master/modelbuild.sh#L486-L490

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that's a good idea, but this whole PR started because a user wanted to control the size of their template image. Also, unless -y 0 is used, the template will re-center itself to minimize the average translation.

do
IMGbase=`basename ${CURRENTIMAGESET[$j]}`
BASENAME=` echo ${IMGbase} | cut -d '.' -f 1 `
COM="${OUTPUT_DIR}/initialCOM${i}_${j}_${IMGbase}"
COMTRANSFORM="${OUTPUT_DIR}/initialCOM${i}_${j}_${BASENAME}.mat"
antsAI -d 3 --convergence 0 --verbose 1 -m Mattes[${TEMPLATES[$i]},${CURRENTIMAGESET[$j]},32,None] -o ${COMTRANSFORM} -t AlignCentersOfMass
antsApplyTransforms -d 3 -r ${TEMPLATES[$i]} -i ${CURRENTIMAGESET[$j]} -t ${COMTRANSFORM} -o ${COM} --verbose
rm -f $COMTRANSFORM
IMAGECOMSET[${#IMAGECOMSET[@]}]=$COM
done
${ANTSPATH}/AverageImages $DIM ${TEMPLATES[$i]} 2 ${IMAGECOMSET[@]}
cookpa marked this conversation as resolved.
Show resolved Hide resolved
# Clean up
rm -f ${IMAGECOMSET[@]}
fi

if [[ ! -s ${TEMPLATES[$i]} ]];
then
echo "Your template : $TEMPLATES[$i] was not created. This indicates trouble! You may want to check correctness of your input parameters. exiting."
exit
fi

# Back up template
intermediateTemplateBase=`basename ${TEMPLATES[$i]}`
cp ${TEMPLATES[$i]} ${intermediateTemplateDir}/initial_${intermediateTemplateBase}

done

# remove old job bash scripts
Expand Down Expand Up @@ -1034,9 +1077,12 @@ if [[ "$RIGID" -eq 1 ]];
done
echo
echo "${ANTSPATH}/AverageImages $DIM ${TEMPLATES[$j]} 2 ${IMAGERIGIDSET[@]}"

# Don't sharpen after rigid alignment
${ANTSPATH}/AverageImages $DIM ${TEMPLATES[$j]} 2 ${IMAGERIGIDSET[@]}

# Don't sharpen after rigid alignment
${ANTSPATH}/AverageImages $DIM ${TEMPLATES[$j]} 2 ${IMAGERIGIDSET[@]}
intermediateTemplateBase=`basename ${TEMPLATES[$j]}`
cp ${TEMPLATES[$j]} ${intermediateTemplateDir}/initialRigid_${intermediateTemplateBase}

done

# cleanup and save output in seperate folder
Expand Down Expand Up @@ -1172,7 +1218,7 @@ while [[ $i -lt ${ITERATIONLIMIT} ]];
rm -f ${outdir}/job*.sh
# Used to save time by only running coarse registration for the first couple of iterations
# This may also help convergence, but because there's no way to turn it off, it makes it harder
# to refine templates with multiple calls to this script.
# to refine templates with multiple calls to this script.
# If you uncomment this, replace MAXITERATIONS with ITERATIONS in the call to ants below
#
# # For the first couple of iterations, use high-level registration only
Expand Down Expand Up @@ -1396,10 +1442,19 @@ while [[ $i -lt ${ITERATIONLIMIT} ]];
exit 1;
fi
fi

for (( j = 0; j < $NUMBEROFMODALITIES; j++ ))
do
do
shapeupdatetotemplate ${DIM} ${TEMPLATES[$j]} ${TEMPLATENAME} ${OUTPUTNAME} ${GRADIENTSTEP} ${STATSMETHOD} ${SHARPENMETHOD} ${j}
done
intermediateTemplateBase=`basename ${TEMPLATES[$j]}`
cp ${TEMPLATES[$j]} ${intermediateTemplateDir}/${TRANSFORMATIONTYPE}_iteration${i}_${intermediateTemplateBase}
done

if [[ -f "${TEMPLATENAME}0warp.nii.gz" ]]
then
cp ${TEMPLATENAME}0warp.nii.gz ${intermediateTemplateDir}/${TRANSFORMATIONTYPE}_iteration${i}_shapeUpdateWarp.nii.gz
fi

if [[ BACKUP_EACH_ITERATION -eq 1 ]];
then
echo
Expand Down
Loading