ENH: GradientBasedAngleOfIncidentImageFilter

ITKKWStyleOverwrite.txt

@@ -3,3 +3,4 @@ include/itkBlockMatchingBayesianRegularizationDisplacementCalculator\.h Comments
 include/itkBlockMatchingMetricImageToDisplacementCalculator\.h Comments Disable
 include/itkBlockMatchingOptimizingInterpolationDisplacementCalculator\.h Comments Disable
 include/itkBlockMatchingParabolicInterpolationDisplacementCalculator\.h Comments Disable
+test/itkAcousticImpulseResponseImageFilterTest\.cxx Namespace Disable

include/itkAcousticImpulseResponseImageFilter.h

@@ -0,0 +1,135 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  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 itkAcousticImpulseResponseImageFilter_h
+#define itkAcousticImpulseResponseImageFilter_h
+
+#include "itkCastImageFilter.h"
+#include "itkCovariantVector.h"
+#include "itkImageToImageFilter.h"
+
+namespace itk
+{
+/** \class AcousticImpulseResponseImageFilter
+ *
+ * \brief Compute the acoustic impulse response along the beam direction.
+ *
+ * This filter computes the acoustic pressure impulse response along a beam
+ * direction.
+ *
+ * \f[
+ * T( \mathbf{x} ) =
+ *   \cos^n \theta \left( \frac{|\nabla Z( \mathbf{x} )|}
+ *     {2 * Z( \mathbf{x} )} \right )
+ * \f]
+ *
+ * where:
+ *
+ * \f{eqnarray*}
+ *   T( \mathbf{x} ) &=& \mbox{acoustic pressure impulse response}
+ *   n             &=& \mbox{specifies angle dependence}
+ *   Z             &=& \mbox{acoustic impedance}
+ *   \cos \theta   &=& \mbox{angle of incidence}
+ * \f}
+ *
+ * This filter requires two inputs.  The first input is the input acoustic
+ * impedance image.  The second input is an angle of incidence image.
+ * \f$n\f$ is specified with the \c AngleDependence parameter, and
+ * defaults to one.
+ *
+ * It is possible to specify the filter used to calculate the gradient
+ * magnitude with \c SetGradientMagnitudeFilter.
+ *
+ * \ingroup ImageToImageFilter
+ * \ingroup Ultrasound
+ */
+template< typename TInputImage, typename TOutputImage,
+  typename TOperatorValue = float >
+class AcousticImpulseResponseImageFilter
+  : public ImageToImageFilter< TInputImage, TOutputImage >
+{
+public:
+  /** Standard class typedefs. */
+  typedef AcousticImpulseResponseImageFilter                  Self;
+  typedef ImageToImageFilter< TInputImage, TOutputImage >     Superclass;
+  typedef SmartPointer< Self >                                Pointer;
+  typedef SmartPointer< const Self >                          ConstPointer;
+
+  /** Method for creation through the object factory. */
+  itkNewMacro( Self );
+
+  /** Run-time type information ( and related methods ). */
+  itkTypeMacro( AcousticImpulseResponseImageFilter, ImageToImageFilter );
+
+  /** Some convenient typedefs. */
+  typedef TInputImage                           InputImageType;
+  typedef typename InputImageType::PointType    OriginType;
+  typedef TOutputImage                          OutputImageType;
+  typedef typename OutputImageType::RegionType  OutputImageRegionType;
+
+  itkStaticConstMacro( ImageDimension, unsigned int,
+    InputImageType::ImageDimension );
+
+  typedef TOperatorValue
+    OperatorValueType;
+  typedef Image< OperatorValueType, ImageDimension >
+    OperatorImageType;
+  typedef ImageToImageFilter< OperatorImageType, OperatorImageType >
+    GradientMagnitudeFilterType;
+
+  /** Set/Get the angle dependence term \c n in the class documentation. */
+  itkSetMacro( AngleDependence, double );
+  itkGetConstMacro( AngleDependence, double );
+
+  /** Set/Get the filter used to compute the gradient magnitude. */
+  itkSetObjectMacro( GradientMagnitudeFilter, GradientMagnitudeFilterType );
+  itkGetObjectMacro( GradientMagnitudeFilter, GradientMagnitudeFilterType );
+
+protected:
+  AcousticImpulseResponseImageFilter( void );
+  virtual ~AcousticImpulseResponseImageFilter( void ) {}
+
+  void PrintSelf( std::ostream & os, Indent indent ) const override;
+
+  void BeforeThreadedGenerateData( void ) override;
+
+  void ThreadedGenerateData( const OutputImageRegionType &
+    outputRegionForThread, ThreadIdType threadId ) override;
+
+private:
+  // purposely not implemented
+  AcousticImpulseResponseImageFilter( const Self & );
+  // purposely not implemented
+  void operator=( const Self & );
+
+  typename GradientMagnitudeFilterType::Pointer m_GradientMagnitudeFilter;
+
+  double m_AngleDependence;
+
+  typedef CastImageFilter< InputImageType, OperatorImageType >
+    CastImageFilterType;
+  typename CastImageFilterType::Pointer m_CastImageFilter;
+
+}; // End class AcousticImpulseResponseImageFilter
+
+} // End namespace itk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "itkAcousticImpulseResponseImageFilter.hxx"
+#endif
+
+#endif // End !defined( __itkAcousticImpulseResponseImageFilter_h )

include/itkAcousticImpulseResponseImageFilter.hxx

@@ -0,0 +1,134 @@
+/*=========================================================================
+ *
+ *  Copyright Insight Software Consortium
+ *
+ *  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 itkAcousticImpulseResponseImageFilter_hxx
+#define itkAcousticImpulseResponseImageFilter_hxx
+
+#include "itkAcousticImpulseResponseImageFilter.h"
+#include "itkGradientMagnitudeImageFilter.h"
+
+namespace itk
+{
+
+template< typename TInputImage, class TOutputImage, class TOperatorValue >
+AcousticImpulseResponseImageFilter< TInputImage, TOutputImage,
+  TOperatorValue >
+::AcousticImpulseResponseImageFilter( void )
+  : m_AngleDependence( 1.0 )
+{
+  this->SetNumberOfRequiredInputs( 2 );
+
+  typedef GradientMagnitudeImageFilter< OperatorImageType,
+    OperatorImageType > DefaultGradientMagnitudeFilterType;
+  this->m_GradientMagnitudeFilter = DefaultGradientMagnitudeFilterType::
+    New();
+
+  this->m_CastImageFilter = CastImageFilterType::New();
+  //Use the ITKv4 Threading Model (call ThreadedGenerateData instead of DynamicThreadedGenerateData)
+  this->DynamicMultiThreadingOff();
+}
+
+
+template< typename TInputImage, class TOutputImage, class TOperatorValue >
+void
+AcousticImpulseResponseImageFilter< TInputImage, TOutputImage,
+  TOperatorValue >
+::BeforeThreadedGenerateData( void )
+{
+  this->m_CastImageFilter->SetInput( this->GetInput() );
+  this->m_GradientMagnitudeFilter->SetNumberOfWorkUnits(
+    this->GetNumberOfWorkUnits() );
+  this->m_GradientMagnitudeFilter->SetInput(
+    this->m_CastImageFilter->GetOutput() );
+  this->m_GradientMagnitudeFilter->Update();
+}
+
+
+template< typename TInputImage, class TOutputImage, class TOperatorValue >
+void
+AcousticImpulseResponseImageFilter< TInputImage, TOutputImage,
+  TOperatorValue >
+::ThreadedGenerateData( const OutputImageRegionType & outputRegionForThread,
+  ThreadIdType threadId )
+{
+  const InputImageType * input = this->GetInput( 0 );
+  const OperatorImageType * angleOfIncidence = this->GetInput( 1 );
+  const OperatorImageType * gradientMagnitude =
+    this->m_GradientMagnitudeFilter->GetOutput();
+  OutputImageType * output = this->GetOutput();
+
+  typedef ImageRegionConstIterator< InputImageType > InputIteratorType;
+  typedef ImageRegionIterator< OutputImageType >     OutputIteratorType;
+
+  typedef ImageRegionConstIterator< OperatorImageType >
+    AngleOfIncidenceIteratorType;
+  typedef ImageRegionConstIterator< OperatorImageType >
+    GradientMagnitudeIteratorType;
+
+  InputIteratorType inputIt( input, outputRegionForThread );
+  AngleOfIncidenceIteratorType angleOfIncidenceIt( angleOfIncidence,
+    outputRegionForThread );
+  GradientMagnitudeIteratorType gradientMagnitudeIt( gradientMagnitude,
+    outputRegionForThread );
+  OutputIteratorType outputIt( output, outputRegionForThread );
+
+  ProgressReporter progress( this, threadId,
+    outputRegionForThread.GetNumberOfPixels() );
+
+  if( this->m_AngleDependence == 1.0 )  // avoid the pow call
+    {
+    for( inputIt.GoToBegin(), angleOfIncidenceIt.GoToBegin(),
+      gradientMagnitudeIt.GoToBegin(), outputIt.GoToBegin();
+      !outputIt.IsAtEnd();
+      ++inputIt, ++angleOfIncidenceIt, ++gradientMagnitudeIt, ++outputIt )
+      {
+      outputIt.Set( angleOfIncidenceIt.Get() *
+        gradientMagnitudeIt.Get() / ( 2.0 * inputIt.Get() ) );
+      }
+    }
+  else
+    {
+    for( inputIt.GoToBegin(), angleOfIncidenceIt.GoToBegin(),
+      gradientMagnitudeIt.GoToBegin(), outputIt.GoToBegin();
+      !outputIt.IsAtEnd();
+      ++inputIt, ++angleOfIncidenceIt, ++gradientMagnitudeIt, ++outputIt )
+      {
+      outputIt.Set( static_cast< typename OutputImageType::PixelType >(
+        std::pow( static_cast< OperatorValueType >(
+            angleOfIncidenceIt.Get() ),
+        static_cast< OperatorValueType >( this->m_AngleDependence *
+        gradientMagnitudeIt.Get() / ( 2.0 * inputIt.Get() ) ) ) ) );
+      }
+    }
+}
+
+
+template< typename TInputImage, typename TOutputImage, typename TOperatorValue >
+void
+AcousticImpulseResponseImageFilter< TInputImage, TOutputImage,
+  TOperatorValue >
+::PrintSelf( std::ostream & os, Indent indent ) const
+{
+  Superclass::PrintSelf( os, indent );
+  os << indent << "AngleDependence: "
+     << this->m_AngleDependence
+     << std::endl;
+}
+
+} // End namespace itk
+
+#endif // End !defined( __itkAcousticImpulseResponseImageFilter_hxx )