From 75ef8206592485118161caeb5fd460459a3650df Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Sun, 10 Nov 2024 22:44:30 +0100 Subject: [PATCH 1/7] draft to normalize weights --- ants/core/__init__.py | 1 + ants/core/ants_image_io.py | 23 ++++++++++++++- ants/registration/build_template.py | 46 +++++++++++++++++++++++------ 3 files changed, 60 insertions(+), 10 deletions(-) diff --git a/ants/core/__init__.py b/ants/core/__init__.py index 665d08f7..531a7782 100644 --- a/ants/core/__init__.py +++ b/ants/core/__init__.py @@ -3,6 +3,7 @@ image_read, dicom_read, image_write, + ones_like, make_image, from_numpy, from_numpy_like, diff --git a/ants/core/ants_image_io.py b/ants/core/ants_image_io.py index 88210318..6549c752 100644 --- a/ants/core/ants_image_io.py +++ b/ants/core/ants_image_io.py @@ -527,4 +527,25 @@ def new_image_like(image, data): has_components=image.has_components) def from_numpy_like(data, image): - return new_image_like(image, data) \ No newline at end of file + return new_image_like(image, data) + + +def ones_like(image: ants.ANTsImage) -> ants.ANTsImage: + """ + Return an image of ones with the same shape and info as a given image. + + Arguments + --------- + image : ANTsImage + Image to get image shape and info from. + + Returns + ------- + ANTsImage + Image of ones with reference header information + """ + ones = image.clone() + ones.view().fill(1) + return ones + + diff --git a/ants/registration/build_template.py b/ants/registration/build_template.py index 35f8b6c3..1bd33dd7 100644 --- a/ants/registration/build_template.py +++ b/ants/registration/build_template.py @@ -6,6 +6,7 @@ import ants + def build_template( initial_template=None, image_list=None, @@ -14,7 +15,7 @@ def build_template( blending_weight=0.75, weights=None, useNoRigid=True, - **kwargs + **kwargs, ): """ Estimate an optimal template from an input image_list @@ -70,30 +71,44 @@ def build_template( weights = [x / sum(weights) for x in weights] if initial_template is None: initial_template = image_list[0] * 0 + wimg = initial_template.clone("float") for i in range(len(image_list)): temp = image_list[i] * weights[i] temp = ants.resample_image_to_target(temp, initial_template) initial_template = initial_template + temp + wtemp = image_list[i].clone("float") + wtemp.view().fill(weights[i]) + wtemp = ants.resample_image_to_target(wtemp, initial_template) + wimg = wimg + wtemp + nonzero = wimg.view() != 0 + initial_template.view()[nonzero] = initial_template.view()[nonzero] / wimg.view()[nonzero] + xavg = initial_template.clone() for i in range(iterations): affinelist = [] for k in range(len(image_list)): w1 = ants.registration( - xavg, image_list[k], type_of_transform=type_of_transform, **kwargs + xavg, image_list[k] * weights[k], type_of_transform=type_of_transform, **kwargs ) L = len(w1["fwdtransforms"]) # affine is the last one - affinelist.append(w1["fwdtransforms"][L-1]) + affinelist.append(w1["fwdtransforms"][L - 1]) if k == 0: if L == 2: wavg = ants.image_read(w1["fwdtransforms"][0]) * weights[k] - xavgNew = w1["warpedmovout"] * weights[k] + xavgNew = w1["warpedmovout"] + wimgNew = ants.apply_transforms(xavg, ants.ones_like(image_list[k]) * weights[k], transformlist=w1["fwdtransforms"]) else: if L == 2: wavg = wavg + ants.image_read(w1["fwdtransforms"][0]) * weights[k] - xavgNew = xavgNew + w1["warpedmovout"] * weights[k] + xavgNew = xavgNew + w1["warpedmovout"] + wimgNew = wimgNew + ants.apply_transforms(xavg, ants.ones_like(image_list[k]) * weights[k], transformlist=w1["fwdtransforms"]) + + # normalize + nonzero = wimgNew.view() != 0 + xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimgNew.view()[nonzero] if useNoRigid: avgaffine = ants.average_affine_transform_no_rigid(affinelist) @@ -108,13 +123,26 @@ def build_template( wavg = wavg * wscl # apply affine to the nonlinear? # need to save the average - wavgA = ants.apply_transforms(fixed = xavgNew, moving = wavg, imagetype=1, transformlist=afffn, whichtoinvert=[1]) + wavgA = ants.apply_transforms( + fixed=xavgNew, + moving=wavg, + imagetype=1, + transformlist=afffn, + whichtoinvert=[1], + ) wavgfn = mktemp(suffix=".nii.gz") ants.image_write(wavgA, wavgfn) - xavg = ants.apply_transforms(fixed=xavgNew, moving=xavgNew, transformlist=[wavgfn, afffn], whichtoinvert=[0, 1]) + xavg = ants.apply_transforms( + fixed=xavgNew, + moving=xavgNew, + transformlist=[wavgfn, afffn], + whichtoinvert=[0, 1], + ) else: - xavg = ants.apply_transforms(fixed=xavgNew, moving=xavgNew, transformlist=[afffn], whichtoinvert=[1]) - + xavg = ants.apply_transforms( + fixed=xavgNew, moving=xavgNew, transformlist=[afffn], whichtoinvert=[1] + ) + os.remove(afffn) if blending_weight is not None: xavg = xavg * blending_weight + ants.iMath(xavg, "Sharpen") * ( From ab4fdf3ebccc711ec424f8e31c4be38ed9418f8a Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Mon, 11 Nov 2024 13:37:48 +0100 Subject: [PATCH 2/7] simplify normalization --- ants/core/__init__.py | 4 +- ants/registration/build_template.py | 60 ++++++++++++++--------------- 2 files changed, 31 insertions(+), 33 deletions(-) diff --git a/ants/core/__init__.py b/ants/core/__init__.py index 531a7782..6d1d459f 100644 --- a/ants/core/__init__.py +++ b/ants/core/__init__.py @@ -3,11 +3,11 @@ image_read, dicom_read, image_write, - ones_like, make_image, from_numpy, from_numpy_like, - new_image_like) + new_image_like, + ones_like) from .ants_image import (ANTsImage, copy_image_info, set_origin, diff --git a/ants/registration/build_template.py b/ants/registration/build_template.py index 1bd33dd7..ae1c9b74 100644 --- a/ants/registration/build_template.py +++ b/ants/registration/build_template.py @@ -2,11 +2,11 @@ import numpy as np import os +import shutil from tempfile import mktemp import ants - def build_template( initial_template=None, image_list=None, @@ -15,7 +15,8 @@ def build_template( blending_weight=0.75, weights=None, useNoRigid=True, - **kwargs, + output_dir=None, + **kwargs ): """ Estimate an optimal template from an input image_list @@ -45,6 +46,10 @@ def build_template( useNoRigid : boolean equivalent of -y in the script. Template update step will not use the rigid component if this is True. + + output_dir : path + directory name where intermediate transforms are written + kwargs : keyword args extra arguments passed to ants registration @@ -61,6 +66,12 @@ def build_template( >>> timage = ants.build_template( image_list = ( image, image2, image3 ) ).resample_image( (45,45)) >>> timagew = ants.build_template( image_list = ( image, image2, image3 ), weights = (5,1,1) ) """ + work_dir = mktemp() if output_dir is None else output_dir + + def make_outprefix(k: int): + os.makedirs(os.path.join(work_dir, f"img{k:04d}"), exist_ok=True) + return os.path.join(work_dir, f"img{k:04d}", "out") + if "type_of_transform" not in kwargs: type_of_transform = "SyN" else: @@ -76,11 +87,8 @@ def build_template( temp = image_list[i] * weights[i] temp = ants.resample_image_to_target(temp, initial_template) initial_template = initial_template + temp - - wtemp = image_list[i].clone("float") - wtemp.view().fill(weights[i]) - wtemp = ants.resample_image_to_target(wtemp, initial_template) - wimg = wimg + wtemp + wtemp = ants.resample_image_to_target(ants.ones_like(image_list[i]), initial_template) + wimg = wimg + wtemp * weights[i] nonzero = wimg.view() != 0 initial_template.view()[nonzero] = initial_template.view()[nonzero] / wimg.view()[nonzero] @@ -89,26 +97,26 @@ def build_template( affinelist = [] for k in range(len(image_list)): w1 = ants.registration( - xavg, image_list[k] * weights[k], type_of_transform=type_of_transform, **kwargs + xavg, image_list[k], type_of_transform=type_of_transform, outprefix=make_outprefix(k), **kwargs ) L = len(w1["fwdtransforms"]) # affine is the last one - affinelist.append(w1["fwdtransforms"][L - 1]) + affinelist.append(w1["fwdtransforms"][L-1]) if k == 0: if L == 2: wavg = ants.image_read(w1["fwdtransforms"][0]) * weights[k] - xavgNew = w1["warpedmovout"] + xavgNew = w1["warpedmovout"] * weights[k] wimgNew = ants.apply_transforms(xavg, ants.ones_like(image_list[k]) * weights[k], transformlist=w1["fwdtransforms"]) else: if L == 2: wavg = wavg + ants.image_read(w1["fwdtransforms"][0]) * weights[k] - xavgNew = xavgNew + w1["warpedmovout"] + xavgNew = xavgNew + w1["warpedmovout"] * weights[k] wimgNew = wimgNew + ants.apply_transforms(xavg, ants.ones_like(image_list[k]) * weights[k], transformlist=w1["fwdtransforms"]) - # normalize - nonzero = wimgNew.view() != 0 - xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimgNew.view()[nonzero] + # normalize + nonzero = wimgNew.view() != 0 + xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimgNew.view()[nonzero] if useNoRigid: avgaffine = ants.average_affine_transform_no_rigid(affinelist) @@ -123,30 +131,20 @@ def build_template( wavg = wavg * wscl # apply affine to the nonlinear? # need to save the average - wavgA = ants.apply_transforms( - fixed=xavgNew, - moving=wavg, - imagetype=1, - transformlist=afffn, - whichtoinvert=[1], - ) + wavgA = ants.apply_transforms(fixed = xavgNew, moving = wavg, imagetype=1, transformlist=afffn, whichtoinvert=[1]) wavgfn = mktemp(suffix=".nii.gz") ants.image_write(wavgA, wavgfn) - xavg = ants.apply_transforms( - fixed=xavgNew, - moving=xavgNew, - transformlist=[wavgfn, afffn], - whichtoinvert=[0, 1], - ) + xavg = ants.apply_transforms(fixed=xavgNew, moving=xavgNew, transformlist=[wavgfn, afffn], whichtoinvert=[0, 1]) + os.remove(wavgfn) else: - xavg = ants.apply_transforms( - fixed=xavgNew, moving=xavgNew, transformlist=[afffn], whichtoinvert=[1] - ) - + xavg = ants.apply_transforms(fixed=xavgNew, moving=xavgNew, transformlist=[afffn], whichtoinvert=[1]) + os.remove(afffn) if blending_weight is not None: xavg = xavg * blending_weight + ants.iMath(xavg, "Sharpen") * ( 1.0 - blending_weight ) + if output_dir is None: + shutil.rmtree(work_dir) return xavg From d955daa35e193d6eafd87803d29bfb7c6a8d1635 Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Mon, 11 Nov 2024 13:52:39 +0100 Subject: [PATCH 3/7] fix doc action --- ants/core/ants_image_io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ants/core/ants_image_io.py b/ants/core/ants_image_io.py index 6549c752..cafbe7fd 100644 --- a/ants/core/ants_image_io.py +++ b/ants/core/ants_image_io.py @@ -530,7 +530,7 @@ def from_numpy_like(data, image): return new_image_like(image, data) -def ones_like(image: ants.ANTsImage) -> ants.ANTsImage: +def ones_like(image): """ Return an image of ones with the same shape and info as a given image. From 10184842883dbbe5b6c76672ab9b5f4d60b693b5 Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Mon, 11 Nov 2024 21:12:12 +0100 Subject: [PATCH 4/7] remove temp file changes, normalize wavg --- ants/registration/build_template.py | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/ants/registration/build_template.py b/ants/registration/build_template.py index ae1c9b74..7746b2c9 100644 --- a/ants/registration/build_template.py +++ b/ants/registration/build_template.py @@ -2,7 +2,6 @@ import numpy as np import os -import shutil from tempfile import mktemp import ants @@ -15,7 +14,6 @@ def build_template( blending_weight=0.75, weights=None, useNoRigid=True, - output_dir=None, **kwargs ): """ @@ -47,9 +45,6 @@ def build_template( equivalent of -y in the script. Template update step will not use the rigid component if this is True. - output_dir : path - directory name where intermediate transforms are written - kwargs : keyword args extra arguments passed to ants registration @@ -66,12 +61,6 @@ def build_template( >>> timage = ants.build_template( image_list = ( image, image2, image3 ) ).resample_image( (45,45)) >>> timagew = ants.build_template( image_list = ( image, image2, image3 ), weights = (5,1,1) ) """ - work_dir = mktemp() if output_dir is None else output_dir - - def make_outprefix(k: int): - os.makedirs(os.path.join(work_dir, f"img{k:04d}"), exist_ok=True) - return os.path.join(work_dir, f"img{k:04d}", "out") - if "type_of_transform" not in kwargs: type_of_transform = "SyN" else: @@ -97,7 +86,7 @@ def make_outprefix(k: int): affinelist = [] for k in range(len(image_list)): w1 = ants.registration( - xavg, image_list[k], type_of_transform=type_of_transform, outprefix=make_outprefix(k), **kwargs + xavg, image_list[k], type_of_transform=type_of_transform, **kwargs ) L = len(w1["fwdtransforms"]) # affine is the last one @@ -117,6 +106,8 @@ def make_outprefix(k: int): # normalize nonzero = wimgNew.view() != 0 xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimgNew.view()[nonzero] + if L == 2: + wavg.view()[nonzero] = wavg.view()[nonzero] / wimgNew.view()[nonzero] if useNoRigid: avgaffine = ants.average_affine_transform_no_rigid(affinelist) @@ -135,7 +126,6 @@ def make_outprefix(k: int): wavgfn = mktemp(suffix=".nii.gz") ants.image_write(wavgA, wavgfn) xavg = ants.apply_transforms(fixed=xavgNew, moving=xavgNew, transformlist=[wavgfn, afffn], whichtoinvert=[0, 1]) - os.remove(wavgfn) else: xavg = ants.apply_transforms(fixed=xavgNew, moving=xavgNew, transformlist=[afffn], whichtoinvert=[1]) @@ -145,6 +135,4 @@ def make_outprefix(k: int): 1.0 - blending_weight ) - if output_dir is None: - shutil.rmtree(work_dir) return xavg From d94761446bbc6c67284174a73f669fe27c5d1d0c Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Mon, 11 Nov 2024 21:22:22 +0100 Subject: [PATCH 5/7] make normalization optional --- ants/registration/build_template.py | 30 +++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/ants/registration/build_template.py b/ants/registration/build_template.py index 7746b2c9..a6ed5f89 100644 --- a/ants/registration/build_template.py +++ b/ants/registration/build_template.py @@ -14,6 +14,7 @@ def build_template( blending_weight=0.75, weights=None, useNoRigid=True, + normalize=False, **kwargs ): """ @@ -44,6 +45,10 @@ def build_template( useNoRigid : boolean equivalent of -y in the script. Template update step will not use the rigid component if this is True. + + normalize : boolean + if this is True, the intensity contribution from each input + image is renormalized on a per-pixel level. kwargs : keyword args extra arguments passed to ants registration @@ -76,10 +81,12 @@ def build_template( temp = image_list[i] * weights[i] temp = ants.resample_image_to_target(temp, initial_template) initial_template = initial_template + temp - wtemp = ants.resample_image_to_target(ants.ones_like(image_list[i]), initial_template) - wimg = wimg + wtemp * weights[i] - nonzero = wimg.view() != 0 - initial_template.view()[nonzero] = initial_template.view()[nonzero] / wimg.view()[nonzero] + if normalize: + for i in range(len(image_list)): + wtemp = ants.resample_image_to_target(ants.ones_like(image_list[i]), wimg) + wimg = wimg + wtemp * weights[i] + nonzero = wimg.view() != 0 + initial_template.view()[nonzero] = initial_template.view()[nonzero] / wimg.view()[nonzero] xavg = initial_template.clone() for i in range(iterations): @@ -96,18 +103,21 @@ def build_template( if L == 2: wavg = ants.image_read(w1["fwdtransforms"][0]) * weights[k] xavgNew = w1["warpedmovout"] * weights[k] - wimgNew = ants.apply_transforms(xavg, ants.ones_like(image_list[k]) * weights[k], transformlist=w1["fwdtransforms"]) + if normalize: + wimg = ants.apply_transforms(xavg, ants.ones_like(image_list[k]), transformlist=w1["fwdtransforms"]) * weights[k] else: if L == 2: wavg = wavg + ants.image_read(w1["fwdtransforms"][0]) * weights[k] xavgNew = xavgNew + w1["warpedmovout"] * weights[k] - wimgNew = wimgNew + ants.apply_transforms(xavg, ants.ones_like(image_list[k]) * weights[k], transformlist=w1["fwdtransforms"]) + if normalize: + wimg = wimg + ants.apply_transforms(xavg, ants.ones_like(image_list[k]), transformlist=w1["fwdtransforms"]) * weights[k] # normalize - nonzero = wimgNew.view() != 0 - xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimgNew.view()[nonzero] - if L == 2: - wavg.view()[nonzero] = wavg.view()[nonzero] / wimgNew.view()[nonzero] + if normalize: + nonzero = wimg.view() != 0 + xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimg.view()[nonzero] + if L == 2: + wavg.view()[nonzero] = wavg.view()[nonzero] / wimg.view()[nonzero] if useNoRigid: avgaffine = ants.average_affine_transform_no_rigid(affinelist) From a0e61546cb903719832f089918880428dac8414b Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Mon, 11 Nov 2024 21:36:00 +0100 Subject: [PATCH 6/7] move wimg to 'if normalize scope' --- ants/registration/build_template.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ants/registration/build_template.py b/ants/registration/build_template.py index a6ed5f89..c4f063a0 100644 --- a/ants/registration/build_template.py +++ b/ants/registration/build_template.py @@ -76,12 +76,12 @@ def build_template( weights = [x / sum(weights) for x in weights] if initial_template is None: initial_template = image_list[0] * 0 - wimg = initial_template.clone("float") for i in range(len(image_list)): temp = image_list[i] * weights[i] temp = ants.resample_image_to_target(temp, initial_template) initial_template = initial_template + temp if normalize: + wimg = initial_template.clone("float") for i in range(len(image_list)): wtemp = ants.resample_image_to_target(ants.ones_like(image_list[i]), wimg) wimg = wimg + wtemp * weights[i] From 334671ce380cb2966acec246ecd923ad077a9128 Mon Sep 17 00:00:00 2001 From: Bryn Lloyd <12702862+dyollb@users.noreply.github.com> Date: Mon, 11 Nov 2024 21:37:15 +0100 Subject: [PATCH 7/7] remove redundant comment --- ants/registration/build_template.py | 1 - 1 file changed, 1 deletion(-) diff --git a/ants/registration/build_template.py b/ants/registration/build_template.py index c4f063a0..99ae15a9 100644 --- a/ants/registration/build_template.py +++ b/ants/registration/build_template.py @@ -112,7 +112,6 @@ def build_template( if normalize: wimg = wimg + ants.apply_transforms(xavg, ants.ones_like(image_list[k]), transformlist=w1["fwdtransforms"]) * weights[k] - # normalize if normalize: nonzero = wimg.view() != 0 xavgNew.view()[nonzero] = xavgNew.view()[nonzero] / wimg.view()[nonzero]