ENH: Ingest ITKBoneEnhancement into Modules/Filtering/BoneEnhancement

Modules/Filtering/BoneEnhancement/CMakeLists.txt

@@ -0,0 +1,3 @@
+project(BoneEnhancement)
+
+itk_module_impl()

Modules/Filtering/BoneEnhancement/include/itkDescoteauxEigenToMeasureImageFilter.h

@@ -0,0 +1,139 @@
+/*=========================================================================
+ *
+ *  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
+ *
+ *         https://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 itkDescoteauxEigenToMeasureImageFilter_h
+#define itkDescoteauxEigenToMeasureImageFilter_h
+
+#include "itkEigenToMeasureImageFilter.h"
+#include "itkMath.h"
+
+namespace itk
+{
+/** \class DescoteauxEigenToMeasureImageFilter
+ * \brief Convert eigenvalues into a measure of sheetness according to the method of Descoteaux et al.
+ *
+ * Converts a 3D fixed array of eigenvalues into a measure of sheetness according to the method
+ * of Descoteaux et al. The parameters of the filter should be set using
+ * DescoteauxEigentoScalarParameterEstimationImageFilter.
+ *
+ * Computes the following equation for eigenvalues in a three dimensional fixed array:
+ *  \f{eqnarray*}{
+ *      R_{sheet} &=&  \frac{|\lambda_2|}{\lambda_3|} \\
+ *      R_{blob} &=&  \frac{|2 |\lambda_3| - |\lambda_2| - |\lambda_1| |}{|\lambda_3|} \\
+ *      R_{noise} &=&  \sqrt{|\lambda_1|^2 + |\lambda_2|^2 + \lambda_3|^2} \\
+ *      s &=& \exp\left(- \frac{R_{sheet}^2}{\alpha^2} \right) \left(1 - \exp\left(- \frac{R_{blob}^2}{\beta^2} \right)
+ * \right) \left(1 - \exp\left(- \frac{R_{noise}^2}{c^2} \right) \right) \f}
+ *
+ * Note that if \f$ \lambda_3 > 0 \f$, \f$ s = 0 \f$.
+ *
+ * \sa DescoteauxEigenToMeasureParameterEstimationFilter
+ * \sa EigenToMeasureImageFilter
+ * \sa MultiScaleHessianEnhancementImageFilter
+ *
+ * \author: Bryce Besler
+ * \ingroup BoneEnhancement
+ */
+template <typename TInputImage, typename TOutputImage>
+class DescoteauxEigenToMeasureImageFilter : public EigenToMeasureImageFilter<TInputImage, TOutputImage>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(DescoteauxEigenToMeasureImageFilter);
+
+  /** Standard Self typedef */
+  using Self = DescoteauxEigenToMeasureImageFilter;
+  using Superclass = EigenToMeasureImageFilter<TInputImage, TOutputImage>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** Input typedefs */
+  using InputImageType = typename Superclass::InputImageType;
+  using InputImagePixelType = typename Superclass::InputImagePixelType;
+  using InputImagePointer = typename Superclass::InputImagePointer;
+  using InputImageConstPointer = typename Superclass::InputImageConstPointer;
+  using InputImageRegionType = typename Superclass::InputImageRegionType;
+
+  /** Output typedefs */
+  using OutputImageType = typename Superclass::OutputImageType;
+  using OutputImagePointer = typename Superclass::OutputImagePointer;
+  using OutputImageRegionType = typename Superclass::OutputImageRegionType;
+  using OutputImagePixelType = typename Superclass::OutputImagePixelType;
+
+  /** Parameter typedefs */
+  using RealType = typename Superclass::RealType;
+  using ParameterArrayType = typename Superclass::ParameterArrayType;
+  using ParameterDecoratedType = typename Superclass::ParameterDecoratedType;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Runtime information support. */
+  itkOverrideGetNameOfClassMacro(DescoteauxEigenToMeasureImageFilter);
+
+  /** Explicitely state the eigenvalues are ordered by magnitude for this filter */
+  typename Superclass::EigenValueOrderEnum
+  GetEigenValueOrder() const override
+  {
+    return Superclass::EigenValueOrderEnum::OrderByMagnitude;
+  }
+
+  /** Setter/Getter methods for setting Direction */
+  itkSetMacro(EnhanceType, RealType);
+  itkGetConstMacro(EnhanceType, RealType);
+  void
+  SetEnhanceBrightObjects()
+  {
+    SetEnhanceType(-1.0);
+  }
+  void
+  SetEnhanceDarkObjects()
+  {
+    SetEnhanceType(1.0);
+  }
+
+#ifdef ITK_USE_CONCEPT_CHECKING
+  // Begin concept checking
+  itkConceptMacro(InputHaveDimension3Check, (Concept::SameDimension<TInputImage::ImageDimension, 3u>));
+  itkConceptMacro(OutputHaveDimension3Check, (Concept::SameDimension<TOutputImage::ImageDimension, 3u>));
+  itkConceptMacro(InputFixedArrayHasDimension3Check, (Concept::SameDimension<TInputImage::PixelType::Dimension, 3u>));
+  // End concept checking
+#endif
+protected:
+  DescoteauxEigenToMeasureImageFilter();
+  ~DescoteauxEigenToMeasureImageFilter() override = default;
+
+  OutputImagePixelType
+  ProcessPixel(const InputImagePixelType & pixel) override;
+
+  /** Check the input has the right number of parameters. */
+  void
+  BeforeThreadedGenerateData() override;
+
+  void
+  PrintSelf(std::ostream & os, Indent indent) const override;
+
+private:
+  /* Member variables */
+  RealType m_EnhanceType;
+}; // end class
+} /* end namespace itk */
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkDescoteauxEigenToMeasureImageFilter.hxx"
+#endif
+
+#endif /* itkDescoteauxEigenToMeasureImageFilter_h */

Modules/Filtering/BoneEnhancement/include/itkDescoteauxEigenToMeasureImageFilter.hxx

@@ -0,0 +1,98 @@
+/*=========================================================================
+ *
+ *  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
+ *
+ *         https://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 itkDescoteauxEigenToMeasureImageFilter_hxx
+#define itkDescoteauxEigenToMeasureImageFilter_hxx
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIterator.h"
+
+namespace itk
+{
+template <typename TInputImage, typename TOutputImage>
+DescoteauxEigenToMeasureImageFilter<TInputImage, TOutputImage>::DescoteauxEigenToMeasureImageFilter()
+  : m_EnhanceType(-1.0)
+{}
+
+template <typename TInputImage, typename TOutputImage>
+void
+DescoteauxEigenToMeasureImageFilter<TInputImage, TOutputImage>::BeforeThreadedGenerateData()
+{
+  ParameterArrayType parameters = this->GetParametersInput()->Get();
+  if (parameters.GetSize() != 3)
+  {
+    itkExceptionMacro(<< "Parameters must have size 3. Given array of size " << parameters.GetSize());
+  }
+}
+
+template <typename TInputImage, typename TOutputImage>
+typename DescoteauxEigenToMeasureImageFilter<TInputImage, TOutputImage>::OutputImagePixelType
+DescoteauxEigenToMeasureImageFilter<TInputImage, TOutputImage>::ProcessPixel(const InputImagePixelType & pixel)
+{
+  /* Grab parameters */
+  ParameterArrayType parameters = this->GetParametersInput()->Get();
+  RealType           alpha = parameters[0];
+  RealType           beta = parameters[1];
+  RealType           c = parameters[2];
+
+  /* Grab pixel values */
+  double sheetness = 0.0;
+  auto   a1 = static_cast<double>(pixel[0]);
+  auto   a2 = static_cast<double>(pixel[1]);
+  auto   a3 = static_cast<double>(pixel[2]);
+  double l1 = itk::Math::abs(a1);
+  double l2 = itk::Math::abs(a2);
+  double l3 = itk::Math::abs(a3);
+
+  /* Deal with l3 > 0 */
+  if (m_EnhanceType * a3 < 0)
+  {
+    return static_cast<OutputImagePixelType>(0.0);
+  }
+
+  /* Avoid divisions by zero (or close to zero) */
+  if (l3 < Math::eps)
+  {
+    return static_cast<OutputImagePixelType>(0.0);
+  }
+
+  /* Compute measures */
+  const double Rsheet = l2 / l3;
+  const double Rblob = itk::Math::abs(2 * l3 - l2 - l1) / l3;
+  const double Rnoise = sqrt(l1 * l1 + l2 * l2 + l3 * l3);
+
+  /* Multiply together to get sheetness */
+  sheetness = 1.0;
+  sheetness *= std::exp(-(Rsheet * Rsheet) / (2 * alpha * alpha));
+  sheetness *= (1.0 - std::exp(-(Rblob * Rblob) / (2 * beta * beta)));
+  sheetness *= (1.0 - std::exp(-(Rnoise * Rnoise) / (2 * c * c)));
+
+  return static_cast<OutputImagePixelType>(sheetness);
+}
+
+template <typename TInputImage, typename TOutputImage>
+void
+DescoteauxEigenToMeasureImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os, Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+  os << indent << "Direction: " << GetEnhanceType() << std::endl;
+}
+
+} // namespace itk
+
+#endif /* itkDescoteauxEigenToMeasureImageFilter_hxx */

Modules/Filtering/BoneEnhancement/include/itkDescoteauxEigenToMeasureParameterEstimationFilter.h

@@ -0,0 +1,138 @@
+/*=========================================================================
+ *
+ *  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
+ *
+ *         https://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 itkDescoteauxEigenToMeasureParameterEstimationFilter_h
+#define itkDescoteauxEigenToMeasureParameterEstimationFilter_h
+
+#include "itkMath.h"
+#include "itkEigenToMeasureParameterEstimationFilter.h"
+#include <mutex>
+
+namespace itk
+{
+/** \class DescoteauxEigenToMeasureParameterEstimationFilter
+ * \brief Automatic parameter estimation as defined by Descoteaux et al.
+ *
+ * The default parameters are:
+ *   \f{eqnarray*}{
+ *      \alpha &=& 0.5 \\
+ *      \beta &=& 0.5 \\
+ *      \gamma &=& 0.5 max\( Frobenius norm \)
+ *   \f}
+ *
+ * Where the Frobenius norm for a real, symmetric matrix is given by
+ * the square root of the sum of squares of the eigenvalues.
+ *
+ * The parameters are estimated over the whole volume unless a mask is given.
+ * If a mask is given, parameters are evaluated only where IsInside returns
+ * true.
+ *
+ * \sa DescoteauxEigenToMeasureImageFilter
+ * \sa EigenToMeasureParameterEstimationFilter
+ * \sa MultiScaleHessianEnhancementImageFilter
+ *
+ * \author: Bryce Besler
+ * \ingroup BoneEnhancement
+ */
+template <typename TInputImage, typename TOutputImage = TInputImage>
+class DescoteauxEigenToMeasureParameterEstimationFilter
+  : public EigenToMeasureParameterEstimationFilter<TInputImage, TOutputImage>
+{
+public:
+  ITK_DISALLOW_COPY_AND_MOVE(DescoteauxEigenToMeasureParameterEstimationFilter);
+
+  /** Standard Self typedef */
+  using Self = DescoteauxEigenToMeasureParameterEstimationFilter;
+  using Superclass = EigenToMeasureParameterEstimationFilter<TInputImage, TOutputImage>;
+  using Pointer = SmartPointer<Self>;
+  using ConstPointer = SmartPointer<const Self>;
+
+  /** Input typedefs */
+  using InputImageType = typename Superclass::InputImageType;
+  using InputImagePointer = typename Superclass::InputImagePointer;
+  using InputImageConstPointer = typename Superclass::InputImageConstPointer;
+  using InputImageRegionType = typename Superclass::InputImageRegionType;
+  using InputImagePixelType = typename Superclass::InputImagePixelType;
+  using PixelValueType = typename Superclass::PixelValueType;
+
+  /** Output typedefs */
+  using OutputImageType = typename Superclass::OutputImageType;
+  using OutputImageRegionType = typename Superclass::OutputImageRegionType;
+  using OutputImagePixelType = typename Superclass::OutputImagePixelType;
+
+  /** Input Mask typedefs. */
+  using MaskSpatialObjectType = typename Superclass::MaskSpatialObjectType;
+  using MaskSpatialObjectTypeConstPointer = typename Superclass::MaskSpatialObjectTypeConstPointer;
+
+  /** Parameter typedefs */
+  using RealType = typename Superclass::RealType;
+  using ParameterArrayType = typename Superclass::ParameterArrayType;
+  using ParameterDecoratedType = typename Superclass::ParameterDecoratedType;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  /** Runtime information support. */
+  itkOverrideGetNameOfClassMacro(DescoteauxEigenToMeasureParameterEstimationFilter);
+
+  /** Setter/Getter methods for setting FrobeniusNormWeight */
+  itkSetMacro(FrobeniusNormWeight, RealType);
+  itkGetConstMacro(FrobeniusNormWeight, RealType);
+
+#ifdef ITK_USE_CONCEPT_CHECKING
+  // Begin concept checking
+  itkConceptMacro(InputHaveDimension3Check, (Concept::SameDimension<TInputImage::ImageDimension, 3u>));
+  itkConceptMacro(InputFixedArrayHasDimension3Check, (Concept::SameDimension<TInputImage::PixelType::Dimension, 3u>));
+  // End concept checking
+#endif
+protected:
+  DescoteauxEigenToMeasureParameterEstimationFilter();
+  ~DescoteauxEigenToMeasureParameterEstimationFilter() override = default;
+
+  /** Initialize some accumulators before the threads run. */
+  void
+  BeforeThreadedGenerateData() override;
+
+  /** Do final mean and variance computation from data accumulated in threads. */
+  void
+  AfterThreadedGenerateData() override;
+
+  /** Multi-thread version GenerateData. */
+  void
+  DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread) override;
+
+  inline RealType
+  CalculateFrobeniusNorm(const InputImagePixelType & pixel) const;
+
+  void
+  PrintSelf(std::ostream & os, Indent indent) const override;
+
+private:
+  /* Member variables */
+  RealType m_FrobeniusNormWeight;
+  RealType m_MaxFrobeniusNorm;
+
+  std::mutex m_Mutex;
+}; // end class
+} // namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#  include "itkDescoteauxEigenToMeasureParameterEstimationFilter.hxx"
+#endif
+
+#endif /* itkDescoteauxEigenToMeasureParameterEstimationFilter_h */