[ITK-users] Metadata Information in FFT images can be misleading.

[Pablo Hernández]

Hey there,

I am contributing with the IsotropicWavelet external module<https://github.com/phcerdan/ITKIsotropicWavelets>, and I am facing some doubts about the meaning of image metadata in frequency-domain images.

Current behavior performing a FFTforward is to copy the input image metadata: spacing,origin, and direction to the output, even though the output is in the frequency domain or dual space. I guess it is better copy it than to lose it, but doesn't mean that the metadata is meaningful in the frequency domain. Spacing information can led the user to think that between each pixel holding a frequency value, that spacing is like a frequency resolution, but it is not.

Some refresh of the lingo to help:
dual space: f ~ 1/T (where T can represent time or space)
Units of frequency are: Hertz = 1 / [ T ], or Radians (Rad = 2pi Hz)
Spacing (or Sampling in the lingo of digital signaling processing) in a spatial domain image is associated with the resolution of the image, and represents the pixel width in physical units.
The Origin is an array holding some relative units to world coordinates that are meaningful for the experimenter.

What does spacing mean in a frequency domain image? The frequency resolution after an FFT is related to the size of the FFT (the size of the image) and the sampling rate of the original image (more here<http://electronics.stackexchange.com/questions/12407/what-is-the-relation-between-fft-length-and-frequency-resolution>).

Freq_bin_resolution = Freq_sampling / N. Where N is the size of the image, and Freq_sampling = 1 / Spacing. All these variables are arrays, with size equal to the dimension of the image.

And the origin, in the case of the output of an FFT, depends on the layout of that particular FFT algorithm. VNL and FFTW share the same layout, where, for example, the zero frequency bin is stored the first index {{0}} , (also note that the physical units of that first index is always 0.0 Hz, regardless of origin, or spacing of the original image).

If we set the Freq_sampling to {{1.0}} (corresponding to a Spacing of {{1.0}}, then the freq resolution will be {{1/N}} (depending on the size!), and the range of frequencies will always be: (-0.5, 0.5] Hz, or (-pi, pi] rads (independent of the size).

So right now, in the wavelet module that works in the frequency domain, and does some shrinkage in this domain, I have chosen to ignore all this metadata, and let the user recover it if he/she needs it, but it might be worth to think about this.

For example, in a shrinkage by a factor of 2, the width of a pixel representing a frequency, doubles. F_bin_resolution =1/(N/2) = 2/N, and the range of allowed frequencies keeps the same: (-Freq_sampling/2.0, Freq_sampling/2.0] or (-0.5, 0.5] or (-pi,pi] if Freq_sampling =1.0

This is all good, but hard to grasp when the metadata can be misleading. If you have any thoughts on how to improve this or be clearer to the user, happy to hear, meanwhile I think I will keep ignoring any metadata associated with the spatial domain when dealing with frequency domain images.

Cheers,
Pablo

[Matt McCormick]

Hi Pablo,

Thanks for discussing this issue and the excellent overview of the
topic. The meaning of metadata in the frequency domain could
definitely a source of confusion and bugs, and we would benefit from a
clarity and consensus on the subject.

For practical reasons of pipeline use and existing algorithm
application, re-use of itk::Image as much about seems beneficial.

What do you think about adding documentation to
itk::ImageBase::GetSpacing() [1]? We could state that the Spacing is
always in physical units. If the image is in the frequency domain,
then Spacing is equal to 1 / Frequency Sampling and the Frequency Bin
Resolution equals 1 / (Frequency Sampling * N), where N is the Size of
the LargestPossibleImageRegion (but not in the half-Hermitian
storage).

Then, the frequency shrinkage filter would increase its output Spacing
accordingly so metadata is updated throughout an analysis pipeline.

Thanks,
Matt

[1] https://itk.org/Doxygen/html/classitk_1_1ImageBase.html#aaadef7c0a9627cf22b0fbf2de6de913c

On Tue, Mar 14, 2017 at 4:54 AM, Pablo Hernández

Post by Pablo Hernández
Hey there,
I am contributing with the IsotropicWavelet external module, and I am facing
some doubts about the meaning of image metadata in frequency-domain images.
Current behavior performing a FFTforward is to copy the input image
metadata: spacing,origin, and direction to the output, even though the
output is in the frequency domain or dual space. I guess it is better copy
it than to lose it, but doesn't mean that the metadata is meaningful in the
frequency domain. Spacing information can led the user to think that between
each pixel holding a frequency value, that spacing is like a frequency
resolution, but it is not.
dual space: f ~ 1/T (where T can represent time or space)
Units of frequency are: Hertz = 1 / [ T ], or Radians (Rad = 2pi Hz)
Spacing (or Sampling in the lingo of digital signaling processing) in a
spatial domain image is associated with the resolution of the image, and
represents the pixel width in physical units.
The Origin is an array holding some relative units to world coordinates that
are meaningful for the experimenter.
What does spacing mean in a frequency domain image? The frequency resolution
after an FFT is related to the size of the FFT (the size of the image) and
the sampling rate of the original image (more here).
Freq_bin_resolution = Freq_sampling / N. Where N is the size of the image,
and Freq_sampling = 1 / Spacing. All these variables are arrays, with size
equal to the dimension of the image.
And the origin, in the case of the output of an FFT, depends on the layout
of that particular FFT algorithm. VNL and FFTW share the same layout, where,
for example, the zero frequency bin is stored the first index {{0}} , (also
note that the physical units of that first index is always 0.0 Hz,
regardless of origin, or spacing of the original image).
If we set the Freq_sampling to {{1.0}} (corresponding to a Spacing of
{{1.0}}, then the freq resolution will be {{1/N}} (depending on the size!),
and the range of frequencies will always be: (-0.5, 0.5] Hz, or (-pi, pi]
rads (independent of the size).
So right now, in the wavelet module that works in the frequency domain, and
does some shrinkage in this domain, I have chosen to ignore all this
metadata, and let the user recover it if he/she needs it, but it might be
worth to think about this.
For example, in a shrinkage by a factor of 2, the width of a pixel
representing a frequency, doubles. F_bin_resolution =1/(N/2) = 2/N, and the
range of allowed frequencies keeps the same: (-Freq_sampling/2.0,
Freq_sampling/2.0] or (-0.5, 0.5] or (-pi,pi] if Freq_sampling =1.0
This is all good, but hard to grasp when the metadata can be misleading. If
you have any thoughts on how to improve this or be clearer to the user,
happy to hear, meanwhile I think I will keep ignoring any metadata
associated with the spatial domain when dealing with frequency domain
images.
Cheers,
Pablo
_____________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

_____________________________________
Powered by www.kitware.com

Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html

Kitware offers ITK Training Courses, for more information visit:
http://www.kitware.com/products/protraining.php

Please keep messages on-topic and check the ITK FAQ at:
http://www.itk.org/Wiki/ITK_FAQ

Follow this link to subscribe/unsubscribe:
http://public.

[Pablo Hernández]

Thanks Matt, I agree that the more appropriate thing to do would be to document better, but maybe only in the ForwardFFT filter instead of ImageBase?

I can imagine experiments where the input images are in the frequency domain: photo detectors in scattering ( from this search [1], these images [2] for example). In this case origin and spacing metadata make sense as they are, even though the experimenters know the physical units of the images are in hertz, not in meters.

The problem with the metadata arises exclusively doing the FFT. If we modify it, we will lose the input image metadata. If we don't, the metadata of the output requires extra knowledge and read the documentation to be meaningful.

Also the Origin is tricky... the first bin, or index {0}, of the output of an FFT depends on the frequency layout of the algorithm, usually in the "standard" layout, it corresponds to zero frequency [3]. So the Origin(Frequency) will be 0.0 if we understand the origin with the value of the first index..., but if the FFT is shifted, the first index would be now a large, negative frequency.

In the ExternalModule, I am using a FrequencyIterator, derived from a regular ImageRegionIterator. It adds a GetFrequency() function, providing an abstraction to let the user stop worrying about the specifics of the frequency layout.
There are 3 different layouts, Regular, FFTLayout, and ShiftedFFTLayout.
- Regular is a normal iterator, with GetFrequency() calling Get(), for the case I commented when the image is taken in the frequency domain, but the user needs to use the IsotropicWavelet module on it (which requires frequency iterators).
- FFTLayout is the standard FFT layout, vnl, fftw, and I think python uses it as well. This is the default type.
- ShiftedFFTLayout, the layout after shifting the FFT.

This allows to iterate an image from the frequency domain, and get the frequency value, the frequency bin, as well as actual value of the current pixel.
And also they are a good place to store frequency 'metadata' information such as FrequencyOrigin, and FrequencySpacing, that depends on the frequency layout: See [4]

Cheers,
Pablo

[1] https://www.google.co.nz/search?q=photo+detector+scattering+result&source=lnms&tbm=isch&sa=X&ved=0ahUKEwjK2NbzntzSAhUEmZQKHZenBU4Q_AUIBigB&biw=1920&bih=947#imgrc=OA28hyPkRDniVM:

[2] http://www.azom.com/article.aspx?ArticleID=12502

[3] https://docs.scipy.org/doc/numpy/reference/routines.fft.html

[4] https://github.com/phcerdan/ITKIsotropicWavelets/blob/master/include/itkFrequencyFFTLayoutImageRegionConstIteratorWithIndex.h#L288

On Fri, Mar 17, 2017 at 8:43 AM, ***@kitware.com<mailto:***@kitware.com> <***@kitware.com<mailto:***@kitware.com>> wrote:
Hi Pablo,

Thanks for discussing this issue and the excellent overview of the
topic. The meaning of metadata in the frequency domain could
definitely a source of confusion and bugs, and we would benefit from a
clarity and consensus on the subject.

For practical reasons of pipeline use and existing algorithm
application, re-use of itk::Image as much about seems beneficial.

What do you think about adding documentation to
itk::ImageBase::GetSpacing() [1]? We could state that the Spacing is
always in physical units. If the image is in the frequency domain,
then Spacing is equal to 1 / Frequency Sampling and the Frequency Bin
Resolution equals 1 / (Frequency Sampling * N), where N is the Size of
the LargestPossibleImageRegion (but not in the half-Hermitian
storage).

Then, the frequency shrinkage filter would increase its output Spacing
accordingly so metadata is updated throughout an analysis pipeline.

Thanks,
Matt

[1] https://itk.org/Doxygen/html/classitk_1_1ImageBase.html#aaadef7c0a9627cf22b0fbf2de6de913c

On Tue, Mar 14, 2017 at 4:54 AM, Pablo Hernández

Post by Pablo Hernández
Hey there,
I am contributing with the IsotropicWavelet external module, and I am facing
some doubts about the meaning of image metadata in frequency-domain images.
Current behavior performing a FFTforward is to copy the input image
metadata: spacing,origin, and direction to the output, even though the
output is in the frequency domain or dual space. I guess it is better copy
it than to lose it, but doesn't mean that the metadata is meaningful in the
frequency domain. Spacing information can led the user to think that between
each pixel holding a frequency value, that spacing is like a frequency
resolution, but it is not.
dual space: f ~ 1/T (where T can represent time or space)
Units of frequency are: Hertz = 1 / [ T ], or Radians (Rad = 2pi Hz)
Spacing (or Sampling in the lingo of digital signaling processing) in a
spatial domain image is associated with the resolution of the image, and
represents the pixel width in physical units.
The Origin is an array holding some relative units to world coordinates that
are meaningful for the experimenter.
What does spacing mean in a frequency domain image? The frequency resolution
after an FFT is related to the size of the FFT (the size of the image) and
the sampling rate of the original image (more here).
Freq_bin_resolution = Freq_sampling / N. Where N is the size of the image,
and Freq_sampling = 1 / Spacing. All these variables are arrays, with size
equal to the dimension of the image.
And the origin, in the case of the output of an FFT, depends on the layout
of that particular FFT algorithm. VNL and FFTW share the same layout, where,
for example, the zero frequency bin is stored the first index {{0}} , (also
note that the physical units of that first index is always 0.0 Hz,
regardless of origin, or spacing of the original image).
If we set the Freq_sampling to {{1.0}} (corresponding to a Spacing of
{{1.0}}, then the freq resolution will be {{1/N}} (depending on the size!),
and the range of frequencies will always be: (-0.5, 0.5] Hz, or (-pi, pi]
rads (independent of the size).
So right now, in the wavelet module that works in the frequency domain, and
does some shrinkage in this domain, I have chosen to ignore all this
metadata, and let the user recover it if he/she needs it, but it might be
worth to think about this.
For example, in a shrinkage by a factor of 2, the width of a pixel
representing a frequency, doubles. F_bin_resolution =1/(N/2) = 2/N, and the
range of allowed frequencies keeps the same: (-Freq_sampling/2.0,
Freq_sampling/2.0] or (-0.5, 0.5] or (-pi,pi] if Freq_sampling =1.0
This is all good, but hard to grasp when the metadata can be misleading. If
you have any thoughts on how to improve this or be clearer to the user,
happy to hear, meanwhile I think I will keep ignoring any metadata
associated with the spatial domain when dealing with frequency domain
images.
Cheers,
Pablo
_____________________________________
Powered by www.kitware.com<http://www.kitware.com>
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Matt McCormick]

On Mon, May 1, 2017 at 10:41 PM, Pablo Hernández <

Post by Pablo Hernández
Thanks Matt, I agree that the more appropriate thing to do would be to
document better, but maybe only in the ForwardFFT filter instead of
ImageBase?

Yes, ForwardFFTImageFilter is a good location, too.


I can imagine experiments where the input images are in the frequency

Post by Pablo Hernández
domain: photo detectors in scattering ( from this search [1], these images
[2] for example). In this case origin and spacing metadata make sense as
they are, even though the experimenters know the physical units of the
images are in hertz, not in meters.

Clarification and a good interface is also important for a common
application encountered in medical imaging: magnetic resonance imaging
collects samples in k-space.

Post by Pablo Hernández
The problem with the metadata arises exclusively doing the FFT. If we
modify it, we will lose the input image metadata. If we don't, the metadata
of the output requires extra knowledge and read the documentation to be
meaningful.
Also the Origin is tricky... the first bin, or index {0}, of the output of
an FFT depends on the frequency layout of the algorithm, usually in the
"standard" layout, it corresponds to zero frequency [3]. So the
Origin(Frequency) will be 0.0 if we understand the origin with the value of
the first index..., but if the FFT is shifted, the first index would be now
a large, negative frequency.
In the ExternalModule, I am using a FrequencyIterator, derived from a
regular ImageRegionIterator. It adds a GetFrequency() function, providing
an abstraction to let the user stop worrying about the specifics of the
frequency layout.
There are 3 different layouts, Regular, FFTLayout, and ShiftedFFTLayout.
- Regular is a normal iterator, with GetFrequency() calling Get(), for the
case I commented when the image is taken in the frequency domain, but the
user needs to use the IsotropicWavelet module on it (which requires
frequency iterators).
- FFTLayout is the standard FFT layout, vnl, fftw, and I think python uses
it as well. This is the default type.
- ShiftedFFTLayout, the layout after shifting the FFT.
This allows to iterate an image from the frequency domain, and get the
frequency value, the frequency bin, as well as actual value of the current
pixel.
And also they are a good place to store frequency 'metadata' information
such as FrequencyOrigin, and FrequencySpacing, that depends on the
frequency layout: See [4]

These iterators with the explicit methods are excellent! The naming makes
their meaning clear, and they are also practically useful for processing
FFT data. Great idea.


Cheers,
Matt

[Pablo Hernández]

Thanks Matt, I agree that the more appropriate thing to do would be to document better, but maybe only in the ForwardFFT filter instead of ImageBase?

I can imagine experiments where the input images are in the frequency domain: photo detectors in scattering ( from this search [1], these images [2] for example). In this case origin and spacing metadata make sense as they are, even though the experimenters know the physical units of the images are in hertz, not in meters.

The problem with the metadata arises exclusively doing the FFT. If we modify it, we will lose the input image metadata. If we don't, the metadata of the output requires extra knowledge and read the documentation to be meaningful.

Also the Origin is tricky... the first bin, or index {0}, of the output of an FFT depends on the frequency layout of the algorithm, usually in the "standard" layout, it corresponds to zero frequency [3]. So the Origin(Frequency) will be 0.0 if we understand the origin with the value of the first index..., but if the FFT is shifted, the first index would be now a large, negative frequency.

In the ExternalModule, I am using a FrequencyIterator, derived from a regular ImageRegionIterator. It adds a GetFrequency() function, providing an abstraction to let the user stop worrying about the specifics of the frequency layout.
There are 3 different layouts, Regular, FFTLayout, and ShiftedFFTLayout.
- Regular is a normal iterator, with GetFrequency() calling Get(), for the case I commented when the image is taken in the frequency domain, but the user needs to use the IsotropicWavelet module on it (which requires frequency iterators).
- FFTLayout is the standard FFT layout, vnl, fftw, and I think python uses it as well. This is the default type.
- ShiftedFFTLayout, the layout after shifting the FFT.

This allows to iterate an image from the frequency domain, and get the frequency value, the frequency bin, as well as actual value of the current pixel.
And also they are a good place to store frequency 'metadata' information such as FrequencyOrigin, and FrequencySpacing, that depends on the frequency layout: See [4]

Cheers,
Pablo

[1] https://www.google.co.nz/search?q=photo+detector+scattering+result&source=lnms&tbm=isch&sa=X&ved=0ahUKEwjK2NbzntzSAhUEmZQKHZenBU4Q_AUIBigB&biw=1920&bih=947#imgrc=OA28hyPkRDniVM:

[2] http://www.azom.com/article.aspx?ArticleID=12502

[3] https://docs.scipy.org/doc/numpy/reference/routines.fft.html

[4] https://github.com/phcerdan/ITKIsotropicWavelets/blob/master/include/itkFrequencyFFTLayoutImageRegionConstIteratorWithIndex.h#L288

On Fri, Mar 17, 2017 at 8:43 AM, ***@kitware.com<mailto:***@kitware.com> <***@kitware.com<mailto:***@kitware.com>> wrote:
Hi Pablo,

Thanks for discussing this issue and the excellent overview of the
topic. The meaning of metadata in the frequency domain could
definitely a source of confusion and bugs, and we would benefit from a
clarity and consensus on the subject.

For practical reasons of pipeline use and existing algorithm
application, re-use of itk::Image as much about seems beneficial.

What do you think about adding documentation to
itk::ImageBase::GetSpacing() [1]? We could state that the Spacing is
always in physical units. If the image is in the frequency domain,
then Spacing is equal to 1 / Frequency Sampling and the Frequency Bin
Resolution equals 1 / (Frequency Sampling * N), where N is the Size of
the LargestPossibleImageRegion (but not in the half-Hermitian
storage).

Then, the frequency shrinkage filter would increase its output Spacing
accordingly so metadata is updated throughout an analysis pipeline.

Thanks,
Matt

[1] https://itk.org/Doxygen/html/classitk_1_1ImageBase.html#aaadef7c0a9627cf22b0fbf2de6de913c

On Tue, Mar 14, 2017 at 4:54 AM, Pablo Hernández

Post by Pablo Hernández
Hey there,
I am contributing with the IsotropicWavelet external module, and I am facing
some doubts about the meaning of image metadata in frequency-domain images.
Current behavior performing a FFTforward is to copy the input image
metadata: spacing,origin, and direction to the output, even though the
output is in the frequency domain or dual space. I guess it is better copy
it than to lose it, but doesn't mean that the metadata is meaningful in the
frequency domain. Spacing information can led the user to think that between
each pixel holding a frequency value, that spacing is like a frequency
resolution, but it is not.
dual space: f ~ 1/T (where T can represent time or space)
Units of frequency are: Hertz = 1 / [ T ], or Radians (Rad = 2pi Hz)
Spacing (or Sampling in the lingo of digital signaling processing) in a
spatial domain image is associated with the resolution of the image, and
represents the pixel width in physical units.
The Origin is an array holding some relative units to world coordinates that
are meaningful for the experimenter.
What does spacing mean in a frequency domain image? The frequency resolution
after an FFT is related to the size of the FFT (the size of the image) and
the sampling rate of the original image (more here).
Freq_bin_resolution = Freq_sampling / N. Where N is the size of the image,
and Freq_sampling = 1 / Spacing. All these variables are arrays, with size
equal to the dimension of the image.
And the origin, in the case of the output of an FFT, depends on the layout
of that particular FFT algorithm. VNL and FFTW share the same layout, where,
for example, the zero frequency bin is stored the first index {{0}} , (also
note that the physical units of that first index is always 0.0 Hz,
regardless of origin, or spacing of the original image).
If we set the Freq_sampling to {{1.0}} (corresponding to a Spacing of
{{1.0}}, then the freq resolution will be {{1/N}} (depending on the size!),
and the range of frequencies will always be: (-0.5, 0.5] Hz, or (-pi, pi]
rads (independent of the size).
So right now, in the wavelet module that works in the frequency domain, and
does some shrinkage in this domain, I have chosen to ignore all this
metadata, and let the user recover it if he/she needs it, but it might be
worth to think about this.
For example, in a shrinkage by a factor of 2, the width of a pixel
representing a frequency, doubles. F_bin_resolution =1/(N/2) = 2/N, and the
range of allowed frequencies keeps the same: (-Freq_sampling/2.0,
Freq_sampling/2.0] or (-0.5, 0.5] or (-pi,pi] if Freq_sampling =1.0
This is all good, but hard to grasp when the metadata can be misleading. If
you have any thoughts on how to improve this or be clearer to the user,
happy to hear, meanwhile I think I will keep ignoring any metadata
associated with the spatial domain when dealing with frequency domain
images.
Cheers,
Pablo
_____________________________________
Powered by www.kitware.com<http://www.kitware.com>
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Pablo Hernández]

Post by Pablo Hernández

Post by Pablo Hernández
- Regular is a normal iterator, with GetFrequency() calling Get(), for the case I commented when the image is taken in the frequency domain, but the user needs to use the IsotropicWavelet module on it (which requires frequency iterators)

Sorry, here, the iterator method GetFrequency() does not calls Get(), it returns an array similar to using TransformIndexToPhysicalPoint on GetIndex() (but ignoring Direction metadata) [1]

[1] https://github.com/phcerdan/ITKIsotropicWavelets/blob/master/include/itkFrequencyImageRegionConstIteratorWithIndex.h#L135

On Tue, May 2, 2017 at 2:41 PM, Pablo Hernández <***@outlook.com<mailto:***@outlook.com>> wrote:
Thanks Matt, I agree that the more appropriate thing to do would be to document better, but maybe only in the ForwardFFT filter instead of ImageBase?

I can imagine experiments where the input images are in the frequency domain: photo detectors in scattering ( from this search [1], these images [2] for example). In this case origin and spacing metadata make sense as they are, even though the experimenters know the physical units of the images are in hertz, not in meters.

The problem with the metadata arises exclusively doing the FFT. If we modify it, we will lose the input image metadata. If we don't, the metadata of the output requires extra knowledge and read the documentation to be meaningful.

Also the Origin is tricky... the first bin, or index {0}, of the output of an FFT depends on the frequency layout of the algorithm, usually in the "standard" layout, it corresponds to zero frequency [3]. So the Origin(Frequency) will be 0.0 if we understand the origin with the value of the first index..., but if the FFT is shifted, the first index would be now a large, negative frequency.

In the ExternalModule, I am using a FrequencyIterator, derived from a regular ImageRegionIterator. It adds a GetFrequency() function, providing an abstraction to let the user stop worrying about the specifics of the frequency layout.
There are 3 different layouts, Regular, FFTLayout, and ShiftedFFTLayout.
- Regular is a normal iterator, with GetFrequency() calling Get(), for the case I commented when the image is taken in the frequency domain, but the user needs to use the IsotropicWavelet module on it (which requires frequency iterators).
- FFTLayout is the standard FFT layout, vnl, fftw, and I think python uses it as well. This is the default type.
- ShiftedFFTLayout, the layout after shifting the FFT.

This allows to iterate an image from the frequency domain, and get the frequency value, the frequency bin, as well as actual value of the current pixel.
And also they are a good place to store frequency 'metadata' information such as FrequencyOrigin, and FrequencySpacing, that depends on the frequency layout: See [4]

Cheers,
Pablo

[1] https://www.google.co.nz/search?q=photo+detector+scattering+result&source=lnms&tbm=isch&sa=X&ved=0ahUKEwjK2NbzntzSAhUEmZQKHZenBU4Q_AUIBigB&biw=1920&bih=947#imgrc=OA28hyPkRDniVM:

[2] http://www.azom.com/article.aspx?ArticleID=12502

[3] https://docs.scipy.org/doc/numpy/reference/routines.fft.html

[4] https://github.com/phcerdan/ITKIsotropicWavelets/blob/master/include/itkFrequencyFFTLayoutImageRegionConstIteratorWithIndex.h#L288

On Fri, Mar 17, 2017 at 8:43 AM, ***@kitware.com<mailto:***@kitware.com> <***@kitware.com<mailto:***@kitware.com>> wrote:
Hi Pablo,

Thanks for discussing this issue and the excellent overview of the
topic. The meaning of metadata in the frequency domain could
definitely a source of confusion and bugs, and we would benefit from a
clarity and consensus on the subject.

For practical reasons of pipeline use and existing algorithm
application, re-use of itk::Image as much about seems beneficial.

What do you think about adding documentation to
itk::ImageBase::GetSpacing() [1]? We could state that the Spacing is
always in physical units. If the image is in the frequency domain,
then Spacing is equal to 1 / Frequency Sampling and the Frequency Bin
Resolution equals 1 / (Frequency Sampling * N), where N is the Size of
the LargestPossibleImageRegion (but not in the half-Hermitian
storage).

Then, the frequency shrinkage filter would increase its output Spacing
accordingly so metadata is updated throughout an analysis pipeline.

Thanks,
Matt

[1] https://itk.org/Doxygen/html/classitk_1_1ImageBase.html#aaadef7c0a9627cf22b0fbf2de6de913c

On Tue, Mar 14, 2017 at 4:54 AM, Pablo Hernández

Post by Pablo Hernández
Hey there,
I am contributing with the IsotropicWavelet external module, and I am facing
some doubts about the meaning of image metadata in frequency-domain images.
Current behavior performing a FFTforward is to copy the input image
metadata: spacing,origin, and direction to the output, even though the
output is in the frequency domain or dual space. I guess it is better copy
it than to lose it, but doesn't mean that the metadata is meaningful in the
frequency domain. Spacing information can led the user to think that between
each pixel holding a frequency value, that spacing is like a frequency
resolution, but it is not.
dual space: f ~ 1/T (where T can represent time or space)
Units of frequency are: Hertz = 1 / [ T ], or Radians (Rad = 2pi Hz)
Spacing (or Sampling in the lingo of digital signaling processing) in a
spatial domain image is associated with the resolution of the image, and
represents the pixel width in physical units.
The Origin is an array holding some relative units to world coordinates that
are meaningful for the experimenter.
What does spacing mean in a frequency domain image? The frequency resolution
after an FFT is related to the size of the FFT (the size of the image) and
the sampling rate of the original image (more here).
Freq_bin_resolution = Freq_sampling / N. Where N is the size of the image,
and Freq_sampling = 1 / Spacing. All these variables are arrays, with size
equal to the dimension of the image.
And the origin, in the case of the output of an FFT, depends on the layout
of that particular FFT algorithm. VNL and FFTW share the same layout, where,
for example, the zero frequency bin is stored the first index {{0}} , (also
note that the physical units of that first index is always 0.0 Hz,
regardless of origin, or spacing of the original image).
If we set the Freq_sampling to {{1.0}} (corresponding to a Spacing of
{{1.0}}, then the freq resolution will be {{1/N}} (depending on the size!),
and the range of frequencies will always be: (-0.5, 0.5] Hz, or (-pi, pi]
rads (independent of the size).
So right now, in the wavelet module that works in the frequency domain, and
does some shrinkage in this domain, I have chosen to ignore all this
metadata, and let the user recover it if he/she needs it, but it might be
worth to think about this.
For example, in a shrinkage by a factor of 2, the width of a pixel
representing a frequency, doubles. F_bin_resolution =1/(N/2) = 2/N, and the
range of allowed frequencies keeps the same: (-Freq_sampling/2.0,
Freq_sampling/2.0] or (-0.5, 0.5] or (-pi,pi] if Freq_sampling =1.0
This is all good, but hard to grasp when the metadata can be misleading. If
you have any thoughts on how to improve this or be clearer to the user,
happy to hear, meanwhile I think I will keep ignoring any metadata
associated with the spatial domain when dealing with frequency domain
images.
Cheers,
Pablo
_____________________________________
Powered by www.kitware.com<http://www.kitware.com>
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users

[Pablo Hernández]

Post by Pablo Hernández

Post by Pablo Hernández
- Regular is a normal iterator, with GetFrequency() calling Get(), for the case I commented when the image is taken in the frequency domain, but the user needs to use the IsotropicWavelet module on it (which requires frequency iterators)

Sorry, here, the iterator method GetFrequency() does not calls Get(), it returns an array similar to using TransformIndexToPhysicalPoint on GetIndex() (but ignoring Direction metadata) [1]

[1] https://github.com/phcerdan/ITKIsotropicWavelets/blob/master/include/itkFrequencyImageRegionConstIteratorWithIndex.h#L135

On Tue, May 2, 2017 at 2:41 PM, Pablo Hernández <***@outlook.com<mailto:***@outlook.com>> wrote:
Thanks Matt, I agree that the more appropriate thing to do would be to document better, but maybe only in the ForwardFFT filter instead of ImageBase?

I can imagine experiments where the input images are in the frequency domain: photo detectors in scattering ( from this search [1], these images [2] for example). In this case origin and spacing metadata make sense as they are, even though the experimenters know the physical units of the images are in hertz, not in meters.

The problem with the metadata arises exclusively doing the FFT. If we modify it, we will lose the input image metadata. If we don't, the metadata of the output requires extra knowledge and read the documentation to be meaningful.

Also the Origin is tricky... the first bin, or index {0}, of the output of an FFT depends on the frequency layout of the algorithm, usually in the "standard" layout, it corresponds to zero frequency [3]. So the Origin(Frequency) will be 0.0 if we understand the origin with the value of the first index..., but if the FFT is shifted, the first index would be now a large, negative frequency.

In the ExternalModule, I am using a FrequencyIterator, derived from a regular ImageRegionIterator. It adds a GetFrequency() function, providing an abstraction to let the user stop worrying about the specifics of the frequency layout.
There are 3 different layouts, Regular, FFTLayout, and ShiftedFFTLayout.
- Regular is a normal iterator, with GetFrequency() calling Get(), for the case I commented when the image is taken in the frequency domain, but the user needs to use the IsotropicWavelet module on it (which requires frequency iterators).
- FFTLayout is the standard FFT layout, vnl, fftw, and I think python uses it as well. This is the default type.
- ShiftedFFTLayout, the layout after shifting the FFT.

This allows to iterate an image from the frequency domain, and get the frequency value, the frequency bin, as well as actual value of the current pixel.
And also they are a good place to store frequency 'metadata' information such as FrequencyOrigin, and FrequencySpacing, that depends on the frequency layout: See [4]

Cheers,
Pablo

[1] https://www.google.co.nz/search?q=photo+detector+scattering+result&source=lnms&tbm=isch&sa=X&ved=0ahUKEwjK2NbzntzSAhUEmZQKHZenBU4Q_AUIBigB&biw=1920&bih=947#imgrc=OA28hyPkRDniVM:

[2] http://www.azom.com/article.aspx?ArticleID=12502

[3] https://docs.scipy.org/doc/numpy/reference/routines.fft.html

[4] https://github.com/phcerdan/ITKIsotropicWavelets/blob/master/include/itkFrequencyFFTLayoutImageRegionConstIteratorWithIndex.h#L288

On Fri, Mar 17, 2017 at 8:43 AM, ***@kitware.com<mailto:***@kitware.com> <***@kitware.com<mailto:***@kitware.com>> wrote:
Hi Pablo,

Thanks for discussing this issue and the excellent overview of the
topic. The meaning of metadata in the frequency domain could
definitely a source of confusion and bugs, and we would benefit from a
clarity and consensus on the subject.

For practical reasons of pipeline use and existing algorithm
application, re-use of itk::Image as much about seems beneficial.

What do you think about adding documentation to
itk::ImageBase::GetSpacing() [1]? We could state that the Spacing is
always in physical units. If the image is in the frequency domain,
then Spacing is equal to 1 / Frequency Sampling and the Frequency Bin
Resolution equals 1 / (Frequency Sampling * N), where N is the Size of
the LargestPossibleImageRegion (but not in the half-Hermitian
storage).

Then, the frequency shrinkage filter would increase its output Spacing
accordingly so metadata is updated throughout an analysis pipeline.

Thanks,
Matt

[1] https://itk.org/Doxygen/html/classitk_1_1ImageBase.html#aaadef7c0a9627cf22b0fbf2de6de913c

On Tue, Mar 14, 2017 at 4:54 AM, Pablo Hernández

Post by Pablo Hernández
Hey there,
I am contributing with the IsotropicWavelet external module, and I am facing
some doubts about the meaning of image metadata in frequency-domain images.
Current behavior performing a FFTforward is to copy the input image
metadata: spacing,origin, and direction to the output, even though the
output is in the frequency domain or dual space. I guess it is better copy
it than to lose it, but doesn't mean that the metadata is meaningful in the
frequency domain. Spacing information can led the user to think that between
each pixel holding a frequency value, that spacing is like a frequency
resolution, but it is not.
dual space: f ~ 1/T (where T can represent time or space)
Units of frequency are: Hertz = 1 / [ T ], or Radians (Rad = 2pi Hz)
Spacing (or Sampling in the lingo of digital signaling processing) in a
spatial domain image is associated with the resolution of the image, and
represents the pixel width in physical units.
The Origin is an array holding some relative units to world coordinates that
are meaningful for the experimenter.
What does spacing mean in a frequency domain image? The frequency resolution
after an FFT is related to the size of the FFT (the size of the image) and
the sampling rate of the original image (more here).
Freq_bin_resolution = Freq_sampling / N. Where N is the size of the image,
and Freq_sampling = 1 / Spacing. All these variables are arrays, with size
equal to the dimension of the image.
And the origin, in the case of the output of an FFT, depends on the layout
of that particular FFT algorithm. VNL and FFTW share the same layout, where,
for example, the zero frequency bin is stored the first index {{0}} , (also
note that the physical units of that first index is always 0.0 Hz,
regardless of origin, or spacing of the original image).
If we set the Freq_sampling to {{1.0}} (corresponding to a Spacing of
{{1.0}}, then the freq resolution will be {{1/N}} (depending on the size!),
and the range of frequencies will always be: (-0.5, 0.5] Hz, or (-pi, pi]
rads (independent of the size).
So right now, in the wavelet module that works in the frequency domain, and
does some shrinkage in this domain, I have chosen to ignore all this
metadata, and let the user recover it if he/she needs it, but it might be
worth to think about this.
For example, in a shrinkage by a factor of 2, the width of a pixel
representing a frequency, doubles. F_bin_resolution =1/(N/2) = 2/N, and the
range of allowed frequencies keeps the same: (-Freq_sampling/2.0,
Freq_sampling/2.0] or (-0.5, 0.5] or (-pi,pi] if Freq_sampling =1.0
This is all good, but hard to grasp when the metadata can be misleading. If
you have any thoughts on how to improve this or be clearer to the user,
happy to hear, meanwhile I think I will keep ignoring any metadata
associated with the spatial domain when dealing with frequency domain
images.
Cheers,
Pablo
_____________________________________
Powered by www.kitware.com<http://www.kitware.com>
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://www.kitware.com/products/protraining.php
http://www.itk.org/Wiki/ITK_FAQ
http://public.kitware.com/mailman/listinfo/insight-users