ENH: Adding Spectra1DNormalizeImageFilter

The regression image consists of all ones,
since we normalize by the same image we use as input.

Co-authored-by: Tom Birdsong <tom.birdsong@kitware.com>

include/itkSpectra1DNormalizeImageFilter.h

@@ -0,0 +1,77 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkSpectra1DNormalizeImageFilter_h
+#define itkSpectra1DNormalizeImageFilter_h
+
+#include "itkImageToImageFilter.h"
+#include "itkFrequencyDomain1DFilterFunction.h"
+
+namespace itk
+{
+/** \class Spectra1DNormalizeImageFilter
+ * \brief Normalize (divide) a spectral image by reference line spectra.
+ *
+ * Single transducer line is the fastest-varying dimension.
+ *
+ * \ingroup Ultrasound
+ */
+template <typename TInputImage, typename TReferenceImage>
+class ITK_TEMPLATE_EXPORT Spectra1DNormalizeImageFilter : public ImageToImageFilter<TInputImage, TInputImage>
+{
+public:
+  ITK_DISALLOW_COPY_AND_ASSIGN(Spectra1DNormalizeImageFilter);
+
+  /** Standard class type alias. */
+  using InputImageType = TInputImage;
+  using OutputImageType = TInputImage;
+  using OutputImageRegionType = typename OutputImageType::RegionType;
+  using ReferenceImageType = TReferenceImage;
+
+  itkStaticConstMacro(ImageDimension, unsigned int, InputImageType::ImageDimension);
+
+  using Self = Spectra1DNormalizeImageFilter;
+  using Superclass = ImageToImageFilter<InputImageType, OutputImageType>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  itkTypeMacro(Spectra1DNormalizeImageFilter, ImageToImageFilter);
+  itkNewMacro(Self);
+
+  void
+  SetReferenceImage(ReferenceImageType * referenceImage)
+  {
+    this->SetInput("ReferenceImage", referenceImage);
+  }
+
+protected:
+  Spectra1DNormalizeImageFilter() { this->AddRequiredInputName("ReferenceImage", 1); }
+  ~Spectra1DNormalizeImageFilter() override = default;
+
+  void
+  GenerateInputRequestedRegion() override;
+
+  void
+  DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
+};
+} // namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkSpectra1DNormalizeImageFilter.hxx"
+#endif
+
+#endif // itkSpectra1DNormalizeImageFilter_h

include/itkSpectra1DNormalizeImageFilter.hxx

@@ -0,0 +1,126 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+#ifndef itkSpectra1DNormalizeImageFilter_hxx
+#define itkSpectra1DNormalizeImageFilter_hxx
+
+
+#include "itkImageScanlineIterator.h"
+#include "itkTotalProgressReporter.h"
+
+namespace itk
+{
+template <typename TInputImage, typename TReferenceImage>
+void
+Spectra1DNormalizeImageFilter<TInputImage, TReferenceImage>::GenerateInputRequestedRegion()
+{
+  // call the superclass' implementation of this method
+  Superclass::GenerateInputRequestedRegion();
+
+  // ImageToImageFilter expects all inputs to be of TInputImage type
+  // so we don't use this->GetInput(), instead we use ProcessObject directly
+  auto * in1 = reinterpret_cast<const ReferenceImageType *>(this->ProcessObject::GetInput(1));
+  auto * reference = const_cast<ReferenceImageType *>(in1);
+  // instead of fiddling with cropping the requested region to reference line,
+  // we simply request the entire reference line
+  reference->SetRequestedRegionToLargestPossibleRegion();
+}
+
+#define DeclareVectorDivideOperator(TypeA, TypeB)                      \
+  template <typename TScalarA, typename TScalarB, unsigned VDimension> \
+  TypeA ElementWiseDivide(const TypeA & a, const TypeB & b)            \
+  {                                                                    \
+    TypeA result{ a };                                                 \
+    for (unsigned i = 0; i < VDimension; ++i)                          \
+    {                                                                  \
+      if (b[i] != 0) /* else implicitly assume b[i] == 1.0 */          \
+      {                                                                \
+        result[i] /= b[i];                                             \
+      }                                                                \
+    }                                                                  \
+    return result;                                                     \
+  }                                                                    \
+  ITK_NOOP_STATEMENT
+
+// TypeA nor TypeB can have a comma in their specification
+// because commas are used to separate macro arguments
+#define ITK_COMMA ,
+
+DeclareVectorDivideOperator(Vector<TScalarA ITK_COMMA VDimension>, VariableLengthVector<TScalarB>);
+DeclareVectorDivideOperator(VariableLengthVector<TScalarA>, Vector<TScalarB ITK_COMMA VDimension>);
+DeclareVectorDivideOperator(Vector<TScalarA ITK_COMMA VDimension>, Vector<TScalarB ITK_COMMA VDimension>);
+
+// VI->VI case needs to be separate, because it does not have compile-time Dimension
+template <typename TScalarA, typename TScalarB>
+VariableLengthVector<TScalarA>
+ElementWiseDivide(const VariableLengthVector<TScalarA> & a, const VariableLengthVector<TScalarB> & b)
+{
+  VariableLengthVector<TScalarA> result{ a };
+
+  itkAssertInDebugAndIgnoreInReleaseMacro(a.GetNumberOfElements() == b.GetNumberOfElements());
+  const unsigned aDimension = a.GetNumberOfElements();
+  for (unsigned i = 0; i < aDimension; ++i)
+  {
+    if (b[i] != 0) /* else implicitly assume b[i] == 1.0 */
+    {
+      result[i] /= b[i];
+    }
+  }
+  return result;
+}
+
+
+template <typename TInputImage, typename TReferenceImage>
+void
+Spectra1DNormalizeImageFilter<TInputImage, TReferenceImage>::DynamicThreadedGenerateData(
+  const OutputImageRegionType & outputRegionForThread)
+{
+  const InputImageType * input = this->GetInput();
+  const auto *           reference = reinterpret_cast<const ReferenceImageType *>(this->ProcessObject::GetInput(1));
+  OutputImageType *      output = this->GetOutput();
+
+  TotalProgressReporter progress(this, output->GetRequestedRegion().GetNumberOfPixels());
+
+  MultiThreaderBase * multiThreader = this->GetMultiThreader();
+  multiThreader->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
+
+  typename ReferenceImageType::IndexType ind1{ 0 }; // initialize all indices to zero
+
+  ImageScanlineConstIterator<InputImageType> iIt(input, outputRegionForThread);
+  ImageScanlineIterator<OutputImageType>     oIt(output, outputRegionForThread);
+
+  while (!iIt.IsAtEnd())
+  {
+    while (!iIt.IsAtEndOfLine())
+    {
+      ind1[0] = iIt.GetIndex()[0]; // set the index along the depth dimension
+
+      oIt.Set(ElementWiseDivide(iIt.Get(), reference->GetPixel(ind1)));
+
+      ++iIt;
+      ++oIt;
+    }
+
+    iIt.NextLine();
+    oIt.NextLine();
+    progress.Completed(outputRegionForThread.GetSize()[0]);
+  }
+}
+
+} // end namespace itk
+
+#endif // itkSpectra1DNormalizeImageFilter_hxx

test/Baseline/itkSpectra1DNormalizeImageFilterTest.mha.sha512

@@ -0,0 +1 @@
+9b257ffe8fb43feae9acf99b154a3235b34e2247869bb399f70bf77d4b4b7bace4c699f3032ba39c81019e339f39f8158d522e5b2e96e28f7079965b55cba6d0

test/Baseline/itkSpectra1DNormalizeImageFilterTestBaseline.mha.sha512

@@ -0,0 +1 @@
+a5e945abada50dcfd3cda04bec80f48fe4de108b4f99b31fd74bc77090df0d1b07f40076d6461812773ef67284519ccc7e7da601b2ee4429a144f4b76e318ed6

test/CMakeLists.txt

@@ -19,6 +19,7 @@ set(UltrasoundTests
   itkSpectra1DSupportWindowImageFilterTest.cxx
   itkSpectra1DSupportWindowToMaskImageFilterTest.cxx
   itkSpectra1DAveragingImageFilterTest.cxx
+  itkSpectra1DNormalizeImageFilterTest.cxx
   itkTimeGainCompensationImageFilterTest.cxx
   itkInverse1DFFTImageFilterTest.cxx
   itkComplexToComplex1DFFTImageFilterTest.cxx
@@ -185,6 +186,16 @@ itk_add_test(NAME itkSpectra1DAveragingImageFilterTest
     DATA{Input/spectra_pixel_image.mha}
     DATA{Input/spectra_pixel_image.mha}
     )
+itk_add_test(NAME itkSpectra1DNormalizeImageFilterTest
+  COMMAND UltrasoundTestDriver
+    --compare
+    DATA{Baseline/itkSpectra1DNormalizeImageFilterTestBaseline.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkSpectra1DNormalizeImageFilterTest.mha
+  itkSpectra1DNormalizeImageFilterTest
+    DATA{Input/spectra_pixel_image.mha}
+    DATA{Input/spectra_pixel_image.mha}
+    ${ITK_TEST_OUTPUT_DIR}/itkSpectra1DNormalizeImageFilterTest.mha
+    )
 itk_add_test(NAME itkTimeGainCompensationImageFilterTest
   COMMAND UltrasoundTestDriver
   --compare