From 0544321a7247767351c78fe951ddf52864721ff8 Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Thu, 28 Mar 2024 17:27:03 +0000 Subject: [PATCH 1/7] Added minor tweak to exclude ve data from standard differencing --- oxasl/image.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/oxasl/image.py b/oxasl/image.py index e13e45d..92a20a6 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -808,9 +808,10 @@ def perf_weighted(self, name=None): :param name: Optional name for returned image. Defaults to original name with suffix ``_pwi`` :return: 3D fsl.data.image.Image. Not an AslImage as timing information lost """ + print('T.O. test') if not name: name = self.name + "_pwi" - if self.iaf != "mp": + if self.iaf != "mp" and self.iaf != "ve": mean_diffdata = self.diff().mean_across_repeats().data if mean_diffdata.ndim > 3: mean_diffdata = np.mean(mean_diffdata, axis=-1) From bc7b7d3a62030ead3774dc695220471ec78e1d9e Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Tue, 2 Apr 2024 16:52:57 +0100 Subject: [PATCH 2/7] Implemented standard deviation across time to generate an approximate PWI for vessel-encoded data for use in e.g. registration --- oxasl/image.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/oxasl/image.py b/oxasl/image.py index 92a20a6..cab0ebc 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -808,14 +808,10 @@ def perf_weighted(self, name=None): :param name: Optional name for returned image. Defaults to original name with suffix ``_pwi`` :return: 3D fsl.data.image.Image. Not an AslImage as timing information lost """ - print('T.O. test') + if not name: name = self.name + "_pwi" - if self.iaf != "mp" and self.iaf != "ve": - mean_diffdata = self.diff().mean_across_repeats().data - if mean_diffdata.ndim > 3: - mean_diffdata = np.mean(mean_diffdata, axis=-1) - else: + if self.iaf == "mp": # Special case for multiphase data - we cannot difference all # time points but we can take the difference between maximal in/out phase # using FFT. We do not make this part of diff() because it is not a 'proper' @@ -824,6 +820,24 @@ def perf_weighted(self, name=None): meandata = self.mean_across_repeats(diff=False) fft = np.fft.fft(meandata.data, axis=-1) mean_diffdata = np.abs(fft[..., 1]) + elif self.iaf == "ve": + # Special case for vessel-encoded data. As above, we do not do "proper" differencing + # here, but aim for a workaround to generate a PWI e.g. for registration. In this case, + # without knowing the specifics of the vessel-encoding procedure, we make the (reasonable) + # assumption that each vessel is encoded reasonably efficiently across the vessel-encoding + # cycles, in which case the standard deviation across cycles will give a reasonable + # representation of the perfusion signal, irrespective of the order in which the vessels + # have been labelled/controlled. + print('T.O. debugging: running ve differencing...') + meandata = self.mean_across_repeats(diff=False) + print('T.O. debugging: meandata.data.shape = ',meandata.data.shape) + mean_diffdata = np.std(meandata.data, axis=-1) + print('T.O. debugging: mean_diffdata.shape = ',mean_diffdata.shape) + else: # Standard differencing + mean_diffdata = self.diff().mean_across_repeats().data + if mean_diffdata.ndim > 3: + mean_diffdata = np.mean(mean_diffdata, axis=-1) + ret = Image(image=mean_diffdata, name=name, header=self.header) return ret From fb321eec5d0d825bd1fa62f6f61498ad8440121b Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Wed, 3 Apr 2024 13:06:59 +0100 Subject: [PATCH 3/7] Previous approach doesn't work for multi-PLD data as SD across volumes is dominated by static tissue. Trying to copy now an approach from the init-loc option of oxasl_ve which averages over PLDs first --- oxasl/image.py | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/oxasl/image.py b/oxasl/image.py index cab0ebc..c16b627 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -824,15 +824,39 @@ def perf_weighted(self, name=None): # Special case for vessel-encoded data. As above, we do not do "proper" differencing # here, but aim for a workaround to generate a PWI e.g. for registration. In this case, # without knowing the specifics of the vessel-encoding procedure, we make the (reasonable) - # assumption that each vessel is encoded reasonably efficiently across the vessel-encoding + # assumption that each vessel is encoded reasonably efficiently encoded across the vessel-encoding # cycles, in which case the standard deviation across cycles will give a reasonable # representation of the perfusion signal, irrespective of the order in which the vessels # have been labelled/controlled. print('T.O. debugging: running ve differencing...') - meandata = self.mean_across_repeats(diff=False) + #meandata = self.mean_across_repeats(diff=False) + #print('T.O. debugging: meandata.data.shape = ',meandata.data.shape) + #mean_diffdata = np.std(meandata.data, axis=-1) + #print('T.O. debugging: mean_diffdata.shape = ',mean_diffdata.shape) + + # Above doesn't work with multi-PLD data since the SD is dominated by static tissue signal changes + # across PLDs + # Try a similar approach to the init-loc functionality which averages across PLDs first + + # First average over repeats + meandata = self.mean_across_repeats(diff=False).reorder(out_order="lrt") print('T.O. debugging: meandata.data.shape = ',meandata.data.shape) - mean_diffdata = np.std(meandata.data, axis=-1) - print('T.O. debugging: mean_diffdata.shape = ',mean_diffdata.shape) + + # Now average over TIs + if self.ntis > 1: + mean_diffdata = np.zeros(list(meandata.data.shape[:3]) + [self.ntc], dtype=np.float32) + print('T.O. debugging: Initialised mean_diffdata.shape = ',mean_diffdata.shape) + for idx in range(self.ntis): + mean_diffdata += meandata[..., idx*self.ntc:(idx+1)*self.ntc] + mean_diffdata /= self.ntis + print('T.O. debugging: After averaging mean_diffdata.shape = ',mean_diffdata.shape) + else: + mean_diffdata = meandata.data + + # Now take the standard deviation across vessel-encodings + mean_diffdata = np.std(mean_diffdata,axis=-1) + print('T.O. debugging: After SD mean_diffdata.shape = ',mean_diffdata.shape) + else: # Standard differencing mean_diffdata = self.diff().mean_across_repeats().data if mean_diffdata.ndim > 3: From 54ebe4d6885018dd7553fc297e012d1b60d125a9 Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Wed, 3 Apr 2024 15:55:11 +0100 Subject: [PATCH 4/7] Restructured to perform multi-PLD averaging for multi-phase and vessel-encoded data --- oxasl/image.py | 72 +++++++++++++++++++++++++++----------------------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/oxasl/image.py b/oxasl/image.py index c16b627..2e48826 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -811,39 +811,23 @@ def perf_weighted(self, name=None): if not name: name = self.name + "_pwi" - if self.iaf == "mp": - # Special case for multiphase data - we cannot difference all - # time points but we can take the difference between maximal in/out phase - # using FFT. We do not make this part of diff() because it is not a 'proper' - # differencing and is just a workaround for generating a PWI, e.g. for - # registration - meandata = self.mean_across_repeats(diff=False) - fft = np.fft.fft(meandata.data, axis=-1) - mean_diffdata = np.abs(fft[..., 1]) - elif self.iaf == "ve": - # Special case for vessel-encoded data. As above, we do not do "proper" differencing - # here, but aim for a workaround to generate a PWI e.g. for registration. In this case, - # without knowing the specifics of the vessel-encoding procedure, we make the (reasonable) - # assumption that each vessel is encoded reasonably efficiently encoded across the vessel-encoding - # cycles, in which case the standard deviation across cycles will give a reasonable - # representation of the perfusion signal, irrespective of the order in which the vessels - # have been labelled/controlled. - print('T.O. debugging: running ve differencing...') - #meandata = self.mean_across_repeats(diff=False) - #print('T.O. debugging: meandata.data.shape = ',meandata.data.shape) - #mean_diffdata = np.std(meandata.data, axis=-1) - #print('T.O. debugging: mean_diffdata.shape = ',mean_diffdata.shape) - - # Above doesn't work with multi-PLD data since the SD is dominated by static tissue signal changes - # across PLDs - # Try a similar approach to the init-loc functionality which averages across PLDs first - - # First average over repeats + + if self.iaf == ("mp" or "ve"): + # Special case for multiphase and vessel-encoded data - we cannot + # difference all time points but we can take the difference between + # maximal in/out phase using FFT for multiphase data or SD across + # encoding cycles for vessel-encoded data. We do not make this part + # of diff() because it is not a 'proper' differencing and is just a + # workaround for generating a PWI, e.g. for registration + + # First take the mean across repeats, leaving TIs (PLDs) as the last dimension meandata = self.mean_across_repeats(diff=False).reorder(out_order="lrt") print('T.O. debugging: meandata.data.shape = ',meandata.data.shape) - - # Now average over TIs + + # If there are multiple TIs (PLDs), then average across these first + # NB. Uses a similar approach to the init-loc functionality in oxasl_ve if self.ntis > 1: + print('T.O. debugging: Averaging across TIs...') mean_diffdata = np.zeros(list(meandata.data.shape[:3]) + [self.ntc], dtype=np.float32) print('T.O. debugging: Initialised mean_diffdata.shape = ',mean_diffdata.shape) for idx in range(self.ntis): @@ -853,9 +837,31 @@ def perf_weighted(self, name=None): else: mean_diffdata = meandata.data - # Now take the standard deviation across vessel-encodings - mean_diffdata = np.std(mean_diffdata,axis=-1) - print('T.O. debugging: After SD mean_diffdata.shape = ',mean_diffdata.shape) + # Now run specific processing for multiphase or vessel-encoded data + if self.iaf == "mp": + print('T.O. debugging: running FFT analysis for multiphase data...') + # Take the FFT across the multiphase dimension + fft = np.fft.fft(mean_diffdata, axis=-1) + + # The first (non-DC) component represents the sinusoidal + # variation across cycles which gives us an approximation to the + # perfusion signal + mean_diffdata = np.abs(fft[..., 1]) + + elif self.iaf == "ve": + # Special case for vessel-encoded data. In this case, without + # knowing the specifics of the vessel-encoding procedure, we + # make the (reasonable) assumption that each vessel is + # reasonably efficiently encoded across the vessel-encoding + # cycles, in which case the standard deviation across cycles + # will give a reasonable representation of the perfusion signal, + # irrespective of the order in which the vessels have been + # labelled/controlled. + print('T.O. debugging: running ve differencing...') + + # Now take the standard deviation across vessel-encodings + mean_diffdata = np.std(mean_diffdata,axis=-1) + print('T.O. debugging: After SD mean_diffdata.shape = ',mean_diffdata.shape) else: # Standard differencing mean_diffdata = self.diff().mean_across_repeats().data From 5cd3cbd53b47b1a80d18ffd61a4865f0e303551d Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Wed, 3 Apr 2024 16:05:43 +0100 Subject: [PATCH 5/7] Added some debug messages --- oxasl/image.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/oxasl/image.py b/oxasl/image.py index 2e48826..1d22825 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -812,6 +812,11 @@ def perf_weighted(self, name=None): if not name: name = self.name + "_pwi" + print('T.O. debugging: self.iaf:',self.iaf) + print('T.O. debugging: self.iaf == "mp":',self.iaf == "mp") + print('T.O. debugging: self.iaf == "ve":',self.iaf == "ve") + print('T.O. debugging: self.iaf == ("mp" or "ve"): ',self.iaf == ("mp" or "ve")) + if self.iaf == ("mp" or "ve"): # Special case for multiphase and vessel-encoded data - we cannot # difference all time points but we can take the difference between From 5ca628b9a8f0cd35ee721d57b2b7e2ec6375bf5e Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Wed, 3 Apr 2024 16:08:05 +0100 Subject: [PATCH 6/7] More debugging... --- oxasl/image.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/oxasl/image.py b/oxasl/image.py index 1d22825..e303b21 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -815,9 +815,9 @@ def perf_weighted(self, name=None): print('T.O. debugging: self.iaf:',self.iaf) print('T.O. debugging: self.iaf == "mp":',self.iaf == "mp") print('T.O. debugging: self.iaf == "ve":',self.iaf == "ve") - print('T.O. debugging: self.iaf == ("mp" or "ve"): ',self.iaf == ("mp" or "ve")) - - if self.iaf == ("mp" or "ve"): + print('T.O. debugging: (self.iaf == "mp") or (self.iaf == "ve"): ',(self.iaf == "mp") or (self.iaf == "ve")) + + if (self.iaf == "mp") or (self.iaf == "ve"): # Special case for multiphase and vessel-encoded data - we cannot # difference all time points but we can take the difference between # maximal in/out phase using FFT for multiphase data or SD across From 072cdeff81c4d5c9fade45b42729a647198651b7 Mon Sep 17 00:00:00 2001 From: Tom Okell Date: Wed, 3 Apr 2024 16:14:46 +0100 Subject: [PATCH 7/7] Variable renaming, removing debug messages --- oxasl/image.py | 26 ++++++-------------------- 1 file changed, 6 insertions(+), 20 deletions(-) diff --git a/oxasl/image.py b/oxasl/image.py index e303b21..8e912e2 100755 --- a/oxasl/image.py +++ b/oxasl/image.py @@ -812,11 +812,6 @@ def perf_weighted(self, name=None): if not name: name = self.name + "_pwi" - print('T.O. debugging: self.iaf:',self.iaf) - print('T.O. debugging: self.iaf == "mp":',self.iaf == "mp") - print('T.O. debugging: self.iaf == "ve":',self.iaf == "ve") - print('T.O. debugging: (self.iaf == "mp") or (self.iaf == "ve"): ',(self.iaf == "mp") or (self.iaf == "ve")) - if (self.iaf == "mp") or (self.iaf == "ve"): # Special case for multiphase and vessel-encoded data - we cannot # difference all time points but we can take the difference between @@ -827,26 +822,21 @@ def perf_weighted(self, name=None): # First take the mean across repeats, leaving TIs (PLDs) as the last dimension meandata = self.mean_across_repeats(diff=False).reorder(out_order="lrt") - print('T.O. debugging: meandata.data.shape = ',meandata.data.shape) # If there are multiple TIs (PLDs), then average across these first # NB. Uses a similar approach to the init-loc functionality in oxasl_ve if self.ntis > 1: - print('T.O. debugging: Averaging across TIs...') - mean_diffdata = np.zeros(list(meandata.data.shape[:3]) + [self.ntc], dtype=np.float32) - print('T.O. debugging: Initialised mean_diffdata.shape = ',mean_diffdata.shape) + mean_acrosstis = np.zeros(list(meandata.data.shape[:3]) + [self.ntc], dtype=np.float32) for idx in range(self.ntis): - mean_diffdata += meandata[..., idx*self.ntc:(idx+1)*self.ntc] - mean_diffdata /= self.ntis - print('T.O. debugging: After averaging mean_diffdata.shape = ',mean_diffdata.shape) + mean_acrosstis += meandata[..., idx*self.ntc:(idx+1)*self.ntc] + mean_acrosstis /= self.ntis else: - mean_diffdata = meandata.data + mean_acrosstis = meandata.data # Now run specific processing for multiphase or vessel-encoded data if self.iaf == "mp": - print('T.O. debugging: running FFT analysis for multiphase data...') # Take the FFT across the multiphase dimension - fft = np.fft.fft(mean_diffdata, axis=-1) + fft = np.fft.fft(mean_acrosstis, axis=-1) # The first (non-DC) component represents the sinusoidal # variation across cycles which gives us an approximation to the @@ -862,11 +852,7 @@ def perf_weighted(self, name=None): # will give a reasonable representation of the perfusion signal, # irrespective of the order in which the vessels have been # labelled/controlled. - print('T.O. debugging: running ve differencing...') - - # Now take the standard deviation across vessel-encodings - mean_diffdata = np.std(mean_diffdata,axis=-1) - print('T.O. debugging: After SD mean_diffdata.shape = ',mean_diffdata.shape) + mean_diffdata = np.std(mean_acrosstis,axis=-1) else: # Standard differencing mean_diffdata = self.diff().mean_across_repeats().data