diff --git a/CMakeLists.txt b/CMakeLists.txt index a4a62d69..4ca1edcf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,6 +4,17 @@ project(LesionSizingToolkit) if(NOT ITK_SOURCE_DIR) find_package(ITK REQUIRED) list(APPEND CMAKE_MODULE_PATH ${ITK_CMAKE_DIR}) +endif() + +# ITK_VERSION_TWEAK is the commit date (YYYYMMDD) of ITK builds after v6.0b02; it is empty before that. +if(ITK_VERSION_TWEAK AND ITK_VERSION_TWEAK GREATER_EQUAL 20260227) + add_compile_definitions(LST_ITK_HAS_ITERATOR_CTAD) +endif() +if(ITK_VERSION_TWEAK AND ITK_VERSION_TWEAK GREATER_EQUAL 20260508) + add_compile_definitions(LST_ITK_HAS_PRINT_NUMERIC_TRAIT) +endif() + +if(NOT ITK_SOURCE_DIR) include(ITKModuleExternal) else() itk_module_impl() diff --git a/Data/Baseline/CannyEdgeDetectionRecursiveGaussianImageFilterTest3.png b/Data/Baseline/CannyEdgeDetectionRecursiveGaussianImageFilterTest3.png new file mode 100644 index 00000000..1cea13e6 Binary files /dev/null and b/Data/Baseline/CannyEdgeDetectionRecursiveGaussianImageFilterTest3.png differ diff --git a/Data/Input/cthead1.png b/Data/Input/cthead1.png new file mode 100644 index 00000000..8c40498a Binary files /dev/null and b/Data/Input/cthead1.png differ diff --git a/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.h b/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.h index 0c1df5b2..dfa6ad86 100644 --- a/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.h +++ b/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.h @@ -18,42 +18,34 @@ #ifndef itkCannyEdgeDetectionRecursiveGaussianImageFilter_h #define itkCannyEdgeDetectionRecursiveGaussianImageFilter_h -#include "itkImageToImageFilter.h" -#include "itkImage.h" -#include "itkFixedArray.h" #include "itkConstNeighborhoodIterator.h" #include "itkSmoothingRecursiveGaussianImageFilter.h" #include "itkMultiplyImageFilter.h" -#include "itkZeroFluxNeumannBoundaryCondition.h" #include "itkMultiThreaderBase.h" #include "itkDerivativeOperator.h" #include "itkSparseFieldLayer.h" #include "itkObjectStore.h" - +#include "itkMath.h" namespace itk { -using MultiThreader = MultiThreaderBase; - - -template -class ListNode +template +class ITK_TEMPLATE_EXPORT ListNode { public: - TValueType m_Value; + TValue m_Value; ListNode * Next; ListNode * Previous; }; - -/** \class CannyEdgeDetectionRecursiveGaussianImageFilter +/** + * \class CannyEdgeDetectionRecursiveGaussianImageFilter + * \brief This filter is an implementation of a Canny edge detector for + * scalar-valued images. * - * This filter is an implementation of a Canny edge detector for scalar-valued - * images. Based on John Canny's paper "A Computational Approach to Edge - * Detection"(IEEE Transactions on Pattern Analysis and Machine Intelligence, - * Vol. PAMI-8, No.6, November 1986), there are four major steps used in the - * edge-detection scheme: + * Based on John Canny's paper \cite canny1986 there are four major steps + * used in the edge-detection scheme: * (1) Smooth the input image with Gaussian filter. * (2) Calculate the second directional derivatives of the smoothed image. * (3) Non-Maximum Suppression: the zero-crossings of 2nd derivative are found, @@ -85,16 +77,19 @@ class ListNode * \todo Edge-linking will be added when an itk connected component labeling * algorithm is available. * - * \ingroup LesionSizingToolkit - * * \sa SmoothingRecursiveGaussianImageFilter * \sa ZeroCrossingImageFilter - * \sa ThresholdImageFilter */ + * \sa ThresholdImageFilter + * \ingroup ITKImageFeature + * \ingroup LesionSizingToolkit + */ template class ITK_TEMPLATE_EXPORT CannyEdgeDetectionRecursiveGaussianImageFilter : public ImageToImageFilter { public: + ITK_DISALLOW_COPY_AND_MOVE(CannyEdgeDetectionRecursiveGaussianImageFilter); + /** Standard "Self" & Superclass type alias. */ using Self = CannyEdgeDetectionRecursiveGaussianImageFilter; using Superclass = ImageToImageFilter; @@ -111,10 +106,7 @@ class ITK_TEMPLATE_EXPORT CannyEdgeDetectionRecursiveGaussianImageFilter using InputImagePixelType = typename TInputImage::PixelType; using OutputImagePixelType = typename TOutputImage::PixelType; using IndexType = typename TInputImage::IndexType; - - using GaussianImageFilterType = SmoothingRecursiveGaussianImageFilter; - using ScalarRealType = typename GaussianImageFilterType::ScalarRealType; - using SigmaArrayType = typename GaussianImageFilterType::SigmaArrayType; + using SizeValueType = itk::SizeValueType; /** The default boundary condition is used unless overridden *in the Evaluate() method. */ @@ -135,17 +127,23 @@ class ITK_TEMPLATE_EXPORT CannyEdgeDetectionRecursiveGaussianImageFilter /** Typedef to describe the output image region type. */ using OutputImageRegionType = typename TOutputImage::RegionType; + using InputImageRegionType = typename TInputImage::RegionType; - /** Run-time type information (and related methods). */ + /** \see LightObject::GetNameOfClass() */ itkOverrideGetNameOfClassMacro(CannyEdgeDetectionRecursiveGaussianImageFilter); - /** ImageDimension constant */ + /** ImageDimension constant. */ static constexpr unsigned int ImageDimension = TInputImage::ImageDimension; static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; - /** Typedef of double containers */ + /** Typedef of double containers. */ using ArrayType = FixedArray; + using ScalarRealType = + typename SmoothingRecursiveGaussianImageFilter::ScalarRealType; + using SigmaArrayType = + typename SmoothingRecursiveGaussianImageFilter::SigmaArrayType; + /** Smoothing parameters for the Gaussian filter. */ void SetSigmaArray(const SigmaArrayType & sigmas); @@ -156,198 +154,100 @@ class ITK_TEMPLATE_EXPORT CannyEdgeDetectionRecursiveGaussianImageFilter ScalarRealType GetSigma() const; - /* Set the Threshold value for detected edges. */ - void - SetThreshold(const OutputImagePixelType th) - { - this->m_Threshold = th; - this->m_UpperThreshold = m_Threshold; - this->m_LowerThreshold = m_Threshold / 2.0; - itkLegacyReplaceBodyMacro(SetThreshold, 2.2, SetUpperThreshold); - } - - OutputImagePixelType - GetThreshold(OutputImagePixelType itkNotUsed(th)) - { - itkLegacyReplaceBodyMacro(GetThreshold, 2.2, GetUpperThreshold); - return this->m_Threshold; - } - - ///* Set the Threshold value for detected edges. */ + /** \brief Set the upper threshold value for detected edges. + */ itkSetMacro(UpperThreshold, OutputImagePixelType); - itkGetMacro(UpperThreshold, OutputImagePixelType); + itkGetConstMacro(UpperThreshold, OutputImagePixelType); + /** \brief Set the lower threshold value for detected edges. + */ itkSetMacro(LowerThreshold, OutputImagePixelType); - itkGetMacro(LowerThreshold, OutputImagePixelType); - - /* Set the Thresholdvalue for detected edges. */ - itkSetMacro(OutsideValue, OutputImagePixelType); - itkGetMacro(OutsideValue, OutputImagePixelType); + itkGetConstMacro(LowerThreshold, OutputImagePixelType); OutputImageType * - GetNonMaximumSuppressionImage() const + GetNonMaximumSuppressionImage() { return this->m_MultiplyImageFilter->GetOutput(); } - /** CannyEdgeDetectionRecursiveGaussianImageFilter needs a larger input requested - * region than the output requested region ( derivative operators, etc). - * As such, CannyEdgeDetectionRecursiveGaussianImageFilter needs to provide an implementation - * for GenerateInputRequestedRegion() in order to inform the - * pipeline execution model. - * - * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ - void - GenerateInputRequestedRegion() noexcept(false) override; - -#ifdef ITK_USE_CONCEPT_CHECKING - /** Begin concept checking */ itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); itkConceptMacro(OutputHasNumericTraitsCheck, (Concept::HasNumericTraits)); itkConceptMacro(SameDimensionCheck, (Concept::SameDimension)); itkConceptMacro(InputIsFloatingPointCheck, (Concept::IsFloatingPoint)); itkConceptMacro(OutputIsFloatingPointCheck, (Concept::IsFloatingPoint)); - /** End concept checking */ -#endif protected: CannyEdgeDetectionRecursiveGaussianImageFilter(); - CannyEdgeDetectionRecursiveGaussianImageFilter(const Self &) {} void PrintSelf(std::ostream & os, Indent indent) const override; void GenerateData() override; + using GaussianImageFilterType = SmoothingRecursiveGaussianImageFilter; using MultiplyImageFilterType = MultiplyImageFilter; private: ~CannyEdgeDetectionRecursiveGaussianImageFilter() override = default; - /** Thread-Data Structure */ - struct CannyThreadStruct - { - CannyEdgeDetectionRecursiveGaussianImageFilter * Filter; - }; - - /** This allocate storage for m_UpdateBuffer, m_UpdateBuffer1 */ + /** Allocate storage for update buffers used during calculation of multiple steps. */ void AllocateUpdateBuffer(); - /** Implement hysteresis thresholding */ + /** Implement hysteresis thresholding. */ void HysteresisThresholding(); - /** Edge linking funciton */ + /** Edge linking function. */ void - FollowEdge(IndexType index); - - /** Check if the index is in bounds or not */ - bool - InBounds(IndexType index); - + FollowEdge(IndexType index, const OutputImageType * multiplyImageFilterOutput); /** Calculate the second derivative of the smoothed image, it writes the - * result to m_UpdateBuffer using the ThreadedCompute2ndDerivative() method - * and multithreading mechanism. */ - void - Compute2ndDerivative(); - - /** - * Split the input into "num" pieces, returning region "i" as - * "splitRegion". This method is called "num" times to return non-overlapping - * regions. The method returns the number of pieces that the input - * can be split into by the routine. i.e. return value is less than or equal - * to "num". - * \sa ImageSource - */ - // virtual - // int SplitUpdateContainer(int i, int num, ThreadRegionType& splitRegion); - - /** Does the actual work of calculating of the 2nd derivative over a region - * supplied by the multithreading mechanism. - * - * \sa Compute2ndDerivative - * \sa Compute2ndDerivativeThreaderCallBack */ + * result to the update buffer */ void - ThreadedCompute2ndDerivative(const OutputImageRegionType & outputRegionForThread, int threadId); + ThreadedCompute2ndDerivative(const OutputImageRegionType & outputRegionForThread); - /** This callback method uses ImageSource::SplitRequestedRegion to acquire an - * output region that it passes to ThreadedCompute2ndDerivative for - * processing. */ - static ITK_THREAD_RETURN_TYPE - Compute2ndDerivativeThreaderCallback(void * arg); - /** This methos is used to calculate the 2nd derivative for + /** This method is used to calculate the 2nd derivative for * non-boundary pixels. It is called by the ThreadedCompute2ndDerivative * method. */ OutputImagePixelType ComputeCannyEdge(const NeighborhoodType & it, void * globalData); /** Calculate the gradient of the second derivative of the smoothed image, - * it writes the result to m_UpdateBuffer1 using the - * ThreadedCompute2ndDerivativePos() method and multithreading mechanism. - */ - void - Compute2ndDerivativePos(); - - /** Does the actual work of calculating of the 2nd derivative over a region - * supplied by the multithreading mechanism. - * - * \sa Compute2ndDerivativePos - * \sa Compute3ndDerivativePosThreaderCallBack */ + * it writes the result to m_UpdateBuffer1 */ void - ThreadedCompute2ndDerivativePos(const OutputImageRegionType & outputRegionForThread, int threadId); - - /**This callback method uses ImageSource::SplitRequestedRegion to acquire an - * output region that it passes to ThreadedCompute2ndDerivative for - * processing. */ - static ITK_THREAD_RETURN_TYPE - Compute2ndDerivativePosThreaderCallback(void * arg); - - /** Standard deviation of the gaussian used for smoothing */ - SigmaArrayType m_Sigma; + ThreadedCompute2ndDerivativePos(const OutputImageRegionType & outputRegionForThread); - /** The maximum error of the gaussian blurring kernel in each dimensional - * direction. */ - ArrayType m_MaximumError; + SigmaArrayType m_Sigma{}; - /** Upper threshold value for identifying edges. */ - OutputImagePixelType m_UpperThreshold; // should be float here? + OutputImagePixelType m_UpperThreshold{}; // should be float here? + OutputImagePixelType m_LowerThreshold{}; // should be float here? - /** Lower threshold value for identifying edges. */ - OutputImagePixelType m_LowerThreshold; // should be float here? + typename OutputImageType::Pointer m_UpdateBuffer1{}; - /** Threshold value for identifying edges. */ - OutputImagePixelType m_Threshold; - - /** "Background" value for use in thresholding. */ - OutputImagePixelType m_OutsideValue; - - /** Update buffers used during calculation of multiple steps */ - typename OutputImageType::Pointer m_UpdateBuffer1; - - /** Gaussian filter to smooth the input image */ - typename GaussianImageFilterType::Pointer m_GaussianFilter; + /** Gaussian filter to smooth the input image. */ + typename GaussianImageFilterType::Pointer m_GaussianFilter{}; /** Multiply image filter to multiply with the zero crossings of the second - * derivative. */ - typename MultiplyImageFilterType::Pointer m_MultiplyImageFilter; + * derivative. */ + typename MultiplyImageFilterType::Pointer m_MultiplyImageFilter{}; - /** Function objects that are used in the inner loops of derivatiVex - calculations. */ - DerivativeOperator m_ComputeCannyEdge1stDerivativeOper; - DerivativeOperator m_ComputeCannyEdge2ndDerivativeOper; + /** Function objects that are used in the inner loops of derivative + * calculations. */ + DerivativeOperator m_ComputeCannyEdge1stDerivativeOper{}; + DerivativeOperator m_ComputeCannyEdge2ndDerivativeOper{}; - std::slice m_ComputeCannyEdgeSlice[ImageDimension]; + std::slice m_ComputeCannyEdgeSlice[ImageDimension]{}; - unsigned long m_Stride[ImageDimension]; - unsigned long m_Center; + SizeValueType m_Stride[ImageDimension]{}; + SizeValueType m_Center{}; - typename ListNodeStorageType::Pointer m_NodeStore; - ListPointerType m_NodeList; -}; + typename ListNodeStorageType::Pointer m_NodeStore{}; + ListPointerType m_NodeList{}; + OutputImageType * m_OutputImage{}; +}; } // end of namespace itk #ifndef ITK_MANUAL_INSTANTIATION diff --git a/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.hxx b/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.hxx index cab667ea..f2195e8f 100644 --- a/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.hxx +++ b/include/itkCannyEdgeDetectionRecursiveGaussianImageFilter.hxx @@ -24,49 +24,46 @@ #include "itkProgressReporter.h" #include "itkGradientMagnitudeImageFilter.h" #include "itkImageRegionIteratorWithIndex.h" -#include +#include "itkMath.h" +#include "itkProgressTransformer.h" +#ifdef LST_ITK_HAS_PRINT_NUMERIC_TRAIT +# include "itkPrintHelper.h" +#endif + namespace itk { - template CannyEdgeDetectionRecursiveGaussianImageFilter::CannyEdgeDetectionRecursiveGaussianImageFilter() + : m_UpperThreshold(OutputImagePixelType{}) + , m_LowerThreshold(OutputImagePixelType{}) + , m_UpdateBuffer1(OutputImageType::New()) + , m_GaussianFilter(GaussianImageFilterType::New()) + , m_MultiplyImageFilter(MultiplyImageFilterType::New()) { - unsigned int i; - m_Sigma.Fill(1.0); - m_OutsideValue = NumericTraits::Zero; - m_Threshold = NumericTraits::Zero; - m_UpperThreshold = NumericTraits::Zero; - m_LowerThreshold = NumericTraits::Zero; - - m_GaussianFilter = GaussianImageFilterType::New(); - m_MultiplyImageFilter = MultiplyImageFilterType::New(); - m_UpdateBuffer1 = OutputImageType::New(); - // Set up neighborhood slices for all the dimensions. - typename Neighborhood::RadiusType r; - r.Fill(1); + constexpr auto r = MakeFilled::RadiusType>(1); - // Dummy neighborhood used to set up the slices. + // Dummy neighborhood used to set up the slices Neighborhood it; it.SetRadius(r); // Slice the neighborhood m_Center = it.Size() / 2; - for (i = 0; i < ImageDimension; ++i) + for (unsigned int i = 0; i < ImageDimension; ++i) { m_Stride[i] = it.GetStride(i); } - for (i = 0; i < ImageDimension; ++i) + for (unsigned int i = 0; i < ImageDimension; ++i) { m_ComputeCannyEdgeSlice[i] = std::slice(m_Center - m_Stride[i], 3, m_Stride[i]); } - // Allocate the derivative operator. + // Allocate the derivative operator m_ComputeCannyEdge1stDerivativeOper.SetDirection(0); m_ComputeCannyEdge1stDerivativeOper.SetOrder(1); m_ComputeCannyEdge1stDerivativeOper.CreateDirectional(); @@ -78,15 +75,16 @@ CannyEdgeDetectionRecursiveGaussianImageFilter void CannyEdgeDetectionRecursiveGaussianImageFilter::AllocateUpdateBuffer() { - // The update buffer looks just like the input. - - typename TInputImage::ConstPointer input = this->GetInput(); + // The update buffer looks just like the input + const typename TInputImage::ConstPointer input = this->GetInput(); m_UpdateBuffer1->CopyInformation(input); m_UpdateBuffer1->SetRequestedRegion(input->GetRequestedRegion()); @@ -94,97 +92,38 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Alloc m_UpdateBuffer1->Allocate(); } -template -void -CannyEdgeDetectionRecursiveGaussianImageFilter::GenerateInputRequestedRegion() noexcept( - false) -{ - // call the superclass' implementation of this method - Superclass::GenerateInputRequestedRegion(); - return; - // get pointers to the input and output - typename Superclass::InputImagePointer inputPtr = const_cast(this->GetInput()); - typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); - - if (!inputPtr || !outputPtr) - { - return; - } - - // Set the kernel size. - unsigned long radius = 1; - - // get a copy of the input requested region (should equal the output - // requested region) - typename TInputImage::RegionType inputRequestedRegion; - inputRequestedRegion = inputPtr->GetRequestedRegion(); - - // pad the input requested region by the operator radius - inputRequestedRegion.PadByRadius(radius); - - // crop the input requested region at the input's largest possible region - if (inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion())) - { - inputPtr->SetRequestedRegion(inputRequestedRegion); - return; - } - else - { - // Couldn't crop the region (requested region is outside the largest - // possible region). Throw an exception. - - // store what we tried to request (prior to trying to crop) - inputPtr->SetRequestedRegion(inputRequestedRegion); - - // build an exception - InvalidRequestedRegionError e(__FILE__, __LINE__); - std::ostringstream msg; - msg << this->GetNameOfClass() << "::GenerateInputRequestedRegion()"; - e.SetLocation(msg.str().c_str()); - e.SetDescription("Requested region is (at least partially) outside the largest possible region."); - e.SetDataObject(inputPtr); - throw e; - } -} - template void CannyEdgeDetectionRecursiveGaussianImageFilter::ThreadedCompute2ndDerivative( - const OutputImageRegionType & outputRegionForThread, - int threadId) + const OutputImageRegionType & outputRegionForThread) { ZeroFluxNeumannBoundaryCondition nbc; - ImageRegionIterator it; - void * globalData = nullptr; - // Here input is the result from the gaussian filter - // output is the update buffer. - typename OutputImageType::Pointer input = m_GaussianFilter->GetOutput(); - typename OutputImageType::Pointer output = this->GetOutput(); + // Here input is the result from the gaussian filter output is the update + // buffer + const typename OutputImageType::Pointer input = m_GaussianFilter->GetOutput(); - // set iterator radius - Size radius; - radius.Fill(1); + // Set iterator radius + constexpr auto radius = Size::Filled(1); // Find the data-set boundary "faces" - typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList; - NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; - faceList = bC(input, outputRegionForThread, radius); - - typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType::iterator fit; + NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; + const typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList = + bC(input, outputRegionForThread, radius); - // support progress methods/callbacks - ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels(), 100, 0.0f, 0.5f); - - // Process the non-boundady region and then each of the boundary faces. + // Process the non-boundary region and then each of the boundary faces. // These are N-d regions which border the edge of the buffer. - for (fit = faceList.begin(); fit != faceList.end(); ++fit) + for (const auto & face : faceList) { - NeighborhoodType bit(radius, input, *fit); + NeighborhoodType bit(radius, input, face); - it = ImageRegionIterator(output, *fit); +#ifdef LST_ITK_HAS_ITERATOR_CTAD + ImageRegionIterator it(this->m_OutputImage, face); +#else + ImageRegionIterator it(this->m_OutputImage, face); +#endif bit.OverrideBoundaryCondition(&nbc); bit.GoToBegin(); @@ -193,42 +132,39 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Threa it.Value() = ComputeCannyEdge(bit, globalData); ++bit; ++it; - progress.CompletedPixel(); } } } template -typename CannyEdgeDetectionRecursiveGaussianImageFilter::OutputImagePixelType +auto CannyEdgeDetectionRecursiveGaussianImageFilter::ComputeCannyEdge( const NeighborhoodType & it, - void * itkNotUsed(globalData)) + void * itkNotUsed(globalData)) -> OutputImagePixelType { - unsigned int i, j; - NeighborhoodInnerProduct innerProduct; + + const NeighborhoodInnerProduct innerProduct; OutputImagePixelType dx[ImageDimension]; OutputImagePixelType dxx[ImageDimension]; OutputImagePixelType dxy[ImageDimension * (ImageDimension - 1) / 2]; - OutputImagePixelType deriv; - OutputImagePixelType gradMag; // double alpha = 0.01; // Calculate 1st & 2nd order derivative - for (i = 0; i < ImageDimension; i++) + for (unsigned int i = 0; i < ImageDimension; ++i) { dx[i] = innerProduct(m_ComputeCannyEdgeSlice[i], it, m_ComputeCannyEdge1stDerivativeOper); dxx[i] = innerProduct(m_ComputeCannyEdgeSlice[i], it, m_ComputeCannyEdge2ndDerivativeOper); } - deriv = NumericTraits::Zero; - int k = 0; + OutputImagePixelType deriv{}; + int k = 0; // Calculate the 2nd derivative - for (i = 0; i < ImageDimension - 1; i++) + for (unsigned int i = 0; i < ImageDimension - 1; ++i) { - for (j = i + 1; j < ImageDimension; j++) + for (unsigned int j = i + 1; j < ImageDimension; ++j) { dxy[k] = 0.25 * it.GetPixel(m_Center - m_Stride[i] - m_Stride[j]) - 0.25 * it.GetPixel(m_Center - m_Stride[i] + m_Stride[j]) - @@ -236,12 +172,12 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Compu 0.25 * it.GetPixel(m_Center + m_Stride[i] + m_Stride[j]); deriv += 2.0 * dx[i] * dx[j] * dxy[k]; - k++; + ++k; } } - gradMag = 0.0001; // alpha * alpha; - for (i = 0; i < ImageDimension; i++) + auto gradMag = static_cast(0.0001); // alpha * alpha; + for (unsigned int i = 0; i < ImageDimension; ++i) { deriv += dx[i] * dx[i] * dxx[i]; gradMag += dx[i] * dx[i]; @@ -252,91 +188,67 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Compu return deriv; } -// Calculate the second derivative -template -void -CannyEdgeDetectionRecursiveGaussianImageFilter::Compute2ndDerivative() -{ - CannyThreadStruct str; - str.Filter = this; - - this->GetMultiThreader()->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); - this->GetMultiThreader()->SetSingleMethod(this->Compute2ndDerivativeThreaderCallback, &str); - - this->GetMultiThreader()->SingleMethodExecute(); -} - -template -ITK_THREAD_RETURN_TYPE -CannyEdgeDetectionRecursiveGaussianImageFilter::Compute2ndDerivativeThreaderCallback( - void * arg) -{ - CannyThreadStruct * str; - - int total, threadId, threadCount; - - threadId = ((MultiThreaderBase::WorkUnitInfo *)(arg))->WorkUnitID; - threadCount = ((MultiThreaderBase::WorkUnitInfo *)(arg))->NumberOfWorkUnits; - - str = (CannyThreadStruct *)(((MultiThreaderBase::WorkUnitInfo *)(arg))->UserData); - - // Execute the actual method with appropriate output region - // first find out how many pieces extent can be split into. - // Using the SplitRequestedRegion method from itk::ImageSource. - OutputImageRegionType splitRegion; - total = str->Filter->SplitRequestedRegion(threadId, threadCount, splitRegion); - - if (threadId < total) - { - str->Filter->ThreadedCompute2ndDerivative(splitRegion, threadId); - } - - return ITK_THREAD_RETURN_DEFAULT_VALUE; -} template void CannyEdgeDetectionRecursiveGaussianImageFilter::GenerateData() { - // Allocate the output - this->GetOutput()->SetBufferedRegion(this->GetOutput()->GetRequestedRegion()); - this->GetOutput()->Allocate(); - - typename InputImageType::ConstPointer input = this->GetInput(); + this->UpdateProgress(0.0f); - typename ZeroCrossingImageFilter::Pointer zeroCrossFilter = - ZeroCrossingImageFilter::New(); + // Use grafting of the input and output of this filter to isolate + // the mini-pipeline and other modifications from the pipeline. + auto input = InputImageType::New(); + input->Graft(const_cast(this->GetInput())); - typename GradientMagnitudeImageFilter::Pointer gradMag = - GradientMagnitudeImageFilter::New(); + // Allocate the output, and graft + Superclass::AllocateOutputs(); + auto output = OutputImageType::New(); + output->Graft(this->GetOutput()); + this->m_OutputImage = output; - typename MultiplyImageFilter::Pointer multFilter = - MultiplyImageFilter::New(); + auto zeroCrossFilter = ZeroCrossingImageFilter::New(); this->AllocateUpdateBuffer(); - // 1.Apply the Gaussian Filter to the input image.------- + // 1.Apply the Gaussian Filter to the input image m_GaussianFilter->SetSigmaArray(this->m_Sigma); m_GaussianFilter->SetNormalizeAcrossScale(true); m_GaussianFilter->SetInput(input); + // Modify to force execution, due to grafting complications + m_GaussianFilter->Modified(); m_GaussianFilter->Update(); + this->UpdateProgress(0.01f); - // 2. Calculate 2nd order directional derivative------- + ProgressTransformer progress1(0.01f, 0.45f, this); + // 2. Calculate 2nd order directional derivative // Calculate the 2nd order directional derivative of the smoothed image. // The output of this filter will be used to store the directional // derivative. - this->Compute2ndDerivative(); - - this->Compute2ndDerivativePos(); - - // 3. Non-maximum suppression---------- + this->GetMultiThreader()->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); + this->GetMultiThreader()->template ParallelizeImageRegion( + this->GetOutput()->GetRequestedRegion(), + [this](const OutputImageRegionType & outputRegionForThread) { + this->ThreadedCompute2ndDerivative(outputRegionForThread); + }, + progress1.GetProcessObject()); + + ProgressTransformer progress2(0.45f, 0.9f, this); + this->GetMultiThreader()->template ParallelizeImageRegion( + this->GetOutput()->GetRequestedRegion(), + [this](const OutputImageRegionType & outputRegionForThread) { + this->ThreadedCompute2ndDerivativePos(outputRegionForThread); + }, + progress2.GetProcessObject()); + + // 3. Non-maximum suppression // Calculate the zero crossings of the 2nd directional derivative and write // the result to output buffer. - zeroCrossFilter->SetInput(this->GetOutput()); + zeroCrossFilter->SetInput(this->m_OutputImage); zeroCrossFilter->Update(); + this->UpdateProgress(0.92f); - // 4. Hysteresis Thresholding--------- + // 4. Hysteresis Thresholding // First get all the edges corresponding to zerocrossings m_MultiplyImageFilter->SetInput1(m_UpdateBuffer1); @@ -346,9 +258,15 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Gener // which is no longer needed, into the m_MultiplyImageFilter. m_MultiplyImageFilter->GraftOutput(m_GaussianFilter->GetOutput()); m_MultiplyImageFilter->Update(); + this->UpdateProgress(0.95f); - // Then do the double threshoulding upon the edge reponses + // Then do the double thresholding upon the edge responses this->HysteresisThresholding(); + this->UpdateProgress(0.99f); + + this->GraftOutput(output); + this->m_OutputImage = nullptr; + this->UpdateProgress(1.0f); } template @@ -358,33 +276,41 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Hyste // This is the Zero crossings of the Second derivative multiplied with the // gradients of the image. HysteresisThresholding of this image should give // the Canny output. - typename OutputImageType::Pointer input = m_MultiplyImageFilter->GetOutput(); - float value; + const typename OutputImageType::Pointer input = m_MultiplyImageFilter->GetOutput(); + float value = NAN; - ListNodeType * node; + ListNodeType * node = nullptr; + // fix me +#ifdef LST_ITK_HAS_ITERATOR_CTAD + ImageRegionIterator oit(input, input->GetRequestedRegion()); +#else ImageRegionIterator oit(input, input->GetRequestedRegion()); +#endif - oit.GoToBegin(); - - ImageRegionIterator uit(this->GetOutput(), this->GetOutput()->GetRequestedRegion()); - uit.GoToBegin(); +#ifdef LST_ITK_HAS_ITERATOR_CTAD + ImageRegionIterator uit(this->m_OutputImage, this->m_OutputImage->GetRequestedRegion()); +#else + ImageRegionIterator uit(this->m_OutputImage, this->m_OutputImage->GetRequestedRegion()); +#endif while (!uit.IsAtEnd()) { - uit.Value() = 0; + uit.Value() = OutputImagePixelType{}; ++uit; } + const OutputImageType * multiplyImageFilterOutput = this->m_MultiplyImageFilter->GetOutput(); while (!oit.IsAtEnd()) { value = oit.Value(); if (value > m_UpperThreshold) { + const auto computedIndex = oit.ComputeIndex(); node = m_NodeStore->Borrow(); - node->m_Value = oit.GetIndex(); + node->m_Value = computedIndex; m_NodeList->PushFront(node); - FollowEdge(oit.GetIndex()); + FollowEdge(computedIndex, multiplyImageFilterOutput); } ++oit; @@ -393,121 +319,97 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Hyste template void -CannyEdgeDetectionRecursiveGaussianImageFilter::FollowEdge(IndexType index) +CannyEdgeDetectionRecursiveGaussianImageFilter::FollowEdge( + IndexType index, + const OutputImageType * multiplyImageFilterOutput) { // This is the Zero crossings of the Second derivative multiplied with the // gradients of the image. HysteresisThresholding of this image should give // the Canny output. - typename OutputImageType::Pointer input = m_MultiplyImageFilter->GetOutput(); + const InputImageRegionType inputRegion = multiplyImageFilterOutput->GetRequestedRegion(); - IndexType nIndex; - IndexType cIndex; - ListNodeType * node; + // Assign iterator radius + constexpr auto radius = Size::Filled(1); - // assign iterator radius - Size radius; - radius.Fill(1); - - ConstNeighborhoodIterator oit(radius, input, input->GetRequestedRegion()); - ImageRegionIteratorWithIndex uit(this->GetOutput(), this->GetOutput()->GetRequestedRegion()); + ConstNeighborhoodIterator oit( + radius, multiplyImageFilterOutput, multiplyImageFilterOutput->GetRequestedRegion()); + ImageRegionIteratorWithIndex uit(this->m_OutputImage, this->m_OutputImage->GetRequestedRegion()); uit.SetIndex(index); - if (uit.Get() == 1) + if (Math::ExactlyEquals(uit.Get(), NumericTraits::OneValue())) { + // Remove the node if we are not going to follow it! + // + // Pop the front node from the list and read its index value. + auto node = m_NodeList->Front(); // get a pointer to the first node + m_NodeList->PopFront(); // unlink the front node + m_NodeStore->Return(node); // return the memory for reuse return; } - int nSize = m_Center * 2 + 1; + const int nSize = m_Center * 2 + 1; while (!m_NodeList->Empty()) { // Pop the front node from the list and read its index value. - node = m_NodeList->Front(); // get a pointer to the first node - cIndex = node->m_Value; // read the value of the first node - m_NodeList->PopFront(); // unlink the front node - m_NodeStore->Return(node); // return the memory for reuse + auto node = m_NodeList->Front(); // Get a pointer to the first node + const auto cIndex = node->m_Value; // Read the value of the first node + m_NodeList->PopFront(); // Unlink the front node + m_NodeStore->Return(node); // Return the memory for reuse // Move iterators to the correct index position. oit.SetLocation(cIndex); uit.SetIndex(cIndex); uit.Value() = 1; - // Search the neighbors for new indicies to add to the list. - for (int i = 0; i < nSize; i++) + // Search the neighbors for new indices to add to the list. + for (int i = 0; i < nSize; ++i) { - nIndex = oit.GetIndex(i); + const auto nIndex = oit.GetIndex(i); uit.SetIndex(nIndex); - if (InBounds(nIndex)) + if (inputRegion.IsInside(nIndex)) { - if (oit.GetPixel(i) > m_LowerThreshold && uit.Value() != 1) + if (oit.GetPixel(i) > m_LowerThreshold && + Math::NotExactlyEquals(uit.Value(), NumericTraits::OneValue())) { - node = m_NodeStore->Borrow(); // get a new node struct - node->m_Value = nIndex; // set its value - m_NodeList->PushFront(node); // add the new node to the list + node = m_NodeStore->Borrow(); // Get a new node struct + node->m_Value = nIndex; // Set its value + m_NodeList->PushFront(node); // Add the new node to the list uit.SetIndex(nIndex); - uit.Value() = 1; + uit.Value() = NumericTraits::OneValue(); } } } } } -template -bool -CannyEdgeDetectionRecursiveGaussianImageFilter::InBounds(IndexType index) -{ - typename InputImageType::ConstPointer input = this->GetInput(); - typename InputImageType::SizeType sz; - sz = (input->GetRequestedRegion()).GetSize(); - - for (unsigned int i = 0; i < ImageDimension; i++) - { - if (index[i] < 0 || index[i] >= static_cast(sz[i])) - { - return false; - } - } - return true; -} - template void CannyEdgeDetectionRecursiveGaussianImageFilter::ThreadedCompute2ndDerivativePos( - const OutputImageRegionType & outputRegionForThread, - int threadId) + const OutputImageRegionType & outputRegionForThread) { - ZeroFluxNeumannBoundaryCondition nbc; ConstNeighborhoodIterator bit; ConstNeighborhoodIterator bit1; - ImageRegionIterator it; - // Here input is the result from the gaussian filter // input1 is the 2nd derivative result // output is the gradient of 2nd derivative - typename OutputImageType::Pointer input1 = this->GetOutput(); - typename OutputImageType::Pointer input = m_GaussianFilter->GetOutput(); - - typename InputImageType::Pointer output = m_UpdateBuffer1; + const typename OutputImageType::Pointer input1 = this->m_OutputImage; + const typename OutputImageType::Pointer input = m_GaussianFilter->GetOutput(); + const typename InputImageType::Pointer output = m_UpdateBuffer1; - // set iterator radius - Size radius; - radius.Fill(1); + // Set iterator radius + constexpr auto radius = Size::Filled(1); // Find the data-set boundary "faces" - typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList; - NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; - faceList = bC(input, outputRegionForThread, radius); + NeighborhoodAlgorithm::ImageBoundaryFacesCalculator bC; + const typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType faceList = + bC(input, outputRegionForThread, radius); - typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator::FaceListType::iterator fit; - - // support progress methods/callbacks - ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels(), 100, 0.5f, 0.5f); - - InputImagePixelType zero = NumericTraits::Zero; + constexpr InputImagePixelType zero{}; OutputImagePixelType dx[ImageDimension]; OutputImagePixelType dx1[ImageDimension]; @@ -520,13 +422,17 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Threa // Process the non-boundary region and then each of the boundary faces. // These are N-d regions which border the edge of the buffer. - NeighborhoodInnerProduct IP; + const NeighborhoodInnerProduct IP; - for (fit = faceList.begin(); fit != faceList.end(); ++fit) + for (const auto & face : faceList) { - bit = ConstNeighborhoodIterator(radius, input, *fit); - bit1 = ConstNeighborhoodIterator(radius, input1, *fit); - it = ImageRegionIterator(output, *fit); + bit = ConstNeighborhoodIterator(radius, input, face); + bit1 = ConstNeighborhoodIterator(radius, input1, face); +#ifdef LST_ITK_HAS_ITERATOR_CTAD + ImageRegionIterator it(output, face); +#else + ImageRegionIterator it(output, face); +#endif bit.OverrideBoundaryCondition(&nbc); bit.GoToBegin(); bit1.GoToBegin(); @@ -534,10 +440,9 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Threa while (!bit.IsAtEnd()) { - gradMag = 0.0001; - for (unsigned int i = 0; i < ImageDimension; i++) + for (unsigned int i = 0; i < ImageDimension; ++i) { dx[i] = IP(m_ComputeCannyEdgeSlice[i], bit, m_ComputeCannyEdge1stDerivativeOper); gradMag += dx[i] * dx[i]; @@ -545,17 +450,14 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Threa dx1[i] = IP(m_ComputeCannyEdgeSlice[i], bit1, m_ComputeCannyEdge1stDerivativeOper); } - gradMag = std::sqrt((double)gradMag); + gradMag = std::sqrt(static_cast(gradMag)); derivPos = zero; - for (unsigned int i = 0; i < ImageDimension; i++) + for (unsigned int i = 0; i < ImageDimension; ++i) { - // First calculate the directional derivative - directional[i] = dx[i] / gradMag; - // calculate gradient of 2nd derivative - + // Calculate gradient of 2nd derivative derivPos += dx1[i] * directional[i]; } @@ -564,54 +466,10 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Threa ++bit; ++bit1; ++it; - progress.CompletedPixel(); } } } -// Calculate the second derivative -template -void -CannyEdgeDetectionRecursiveGaussianImageFilter::Compute2ndDerivativePos() -{ - CannyThreadStruct str; - str.Filter = this; - - this->GetMultiThreader()->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits()); - this->GetMultiThreader()->SetSingleMethod(this->Compute2ndDerivativePosThreaderCallback, &str); - - this->GetMultiThreader()->SingleMethodExecute(); -} - -template -ITK_THREAD_RETURN_TYPE -CannyEdgeDetectionRecursiveGaussianImageFilter::Compute2ndDerivativePosThreaderCallback( - void * arg) -{ - CannyThreadStruct * str; - - int total, threadId, threadCount; - - threadId = ((MultiThreaderBase::WorkUnitInfo *)(arg))->WorkUnitID; - threadCount = ((MultiThreaderBase::WorkUnitInfo *)(arg))->NumberOfWorkUnits; - - str = (CannyThreadStruct *)(((MultiThreaderBase::WorkUnitInfo *)(arg))->UserData); - - // Execute the actual method with appropriate output region - // first find out how many pieces extent can be split into. - // Using the SplitRequestedRegion method from itk::ImageSource. - - OutputImageRegionType splitRegion; - total = str->Filter->SplitRequestedRegion(threadId, threadCount, splitRegion); - - if (threadId < total) - { - str->Filter->ThreadedCompute2ndDerivativePos(splitRegion, threadId); - } - - return ITK_THREAD_RETURN_DEFAULT_VALUE; -} - // Set value of Sigma (isotropic) template @@ -664,26 +522,22 @@ CannyEdgeDetectionRecursiveGaussianImageFilter::Print Superclass::PrintSelf(os, indent); os << "Sigma: " << m_Sigma << std::endl; - os << indent << "Threshold: " << static_cast::PrintType>(m_Threshold) - << std::endl; +#ifdef LST_ITK_HAS_PRINT_NUMERIC_TRAIT + print_helper::PrintNumericTrait(os, indent, "UpperThreshold", m_UpperThreshold); + print_helper::PrintNumericTrait(os, indent, "LowerThreshold", m_LowerThreshold); +#else os << indent << "UpperThreshold: " << static_cast::PrintType>(m_UpperThreshold) << std::endl; os << indent << "LowerThreshold: " << static_cast::PrintType>(m_LowerThreshold) << std::endl; - os << indent - << "OutsideValue: " << static_cast::PrintType>(m_OutsideValue) - << std::endl; - os << "Center: " << m_Center << std::endl; - os << "Stride: " << m_Stride << std::endl; - os << "Gaussian Filter: " << std::endl; - m_GaussianFilter->Print(os, indent.GetNextIndent()); - os << "Multiply image Filter: " << std::endl; - m_MultiplyImageFilter->Print(os, indent.GetNextIndent()); - os << "UpdateBuffer1: " << std::endl; - m_UpdateBuffer1->Print(os, indent.GetNextIndent()); +#endif + os << indent << "Center: " << m_Center << std::endl; + os << indent << "Stride: " << m_Stride << std::endl; + itkPrintSelfObjectMacro(GaussianFilter); + itkPrintSelfObjectMacro(MultiplyImageFilter); + itkPrintSelfObjectMacro(UpdateBuffer1); } - } // namespace itk #endif diff --git a/include/itkCannyEdgesDistanceAdvectionFieldFeatureGenerator.hxx b/include/itkCannyEdgesDistanceAdvectionFieldFeatureGenerator.hxx index 86ef1bb3..19200c54 100644 --- a/include/itkCannyEdgesDistanceAdvectionFieldFeatureGenerator.hxx +++ b/include/itkCannyEdgesDistanceAdvectionFieldFeatureGenerator.hxx @@ -109,7 +109,6 @@ CannyEdgesDistanceAdvectionFieldFeatureGenerator::GenerateData() this->m_CannyFilter->SetSigma(this->m_Sigma); this->m_CannyFilter->SetUpperThreshold(this->m_UpperThreshold); this->m_CannyFilter->SetLowerThreshold(this->m_LowerThreshold); - this->m_CannyFilter->SetOutsideValue(NumericTraits::Zero); this->m_DistanceMapFilter->Update(); diff --git a/include/itkCannyEdgesDistanceFeatureGenerator.hxx b/include/itkCannyEdgesDistanceFeatureGenerator.hxx index 0ad470c2..35b47281 100644 --- a/include/itkCannyEdgesDistanceFeatureGenerator.hxx +++ b/include/itkCannyEdgesDistanceFeatureGenerator.hxx @@ -106,7 +106,6 @@ CannyEdgesDistanceFeatureGenerator::GenerateData() this->m_CannyFilter->SetSigmaArray(this->m_Sigma); this->m_CannyFilter->SetUpperThreshold(this->m_UpperThreshold); this->m_CannyFilter->SetLowerThreshold(this->m_LowerThreshold); - this->m_CannyFilter->SetOutsideValue(NumericTraits::Zero); this->m_DistanceMapFilter->Update(); diff --git a/include/itkCannyEdgesFeatureGenerator.hxx b/include/itkCannyEdgesFeatureGenerator.hxx index fce5d680..022bae98 100644 --- a/include/itkCannyEdgesFeatureGenerator.hxx +++ b/include/itkCannyEdgesFeatureGenerator.hxx @@ -112,7 +112,6 @@ CannyEdgesFeatureGenerator::GenerateData() this->m_CannyFilter->SetSigmaArray(this->m_Sigma); this->m_CannyFilter->SetUpperThreshold(this->m_UpperThreshold); this->m_CannyFilter->SetLowerThreshold(this->m_LowerThreshold); - this->m_CannyFilter->SetOutsideValue(NumericTraits::Zero); this->m_RescaleFilter->Update(); diff --git a/src/itkCannyEdgeDetectionImageFilter2.cxx b/src/itkCannyEdgeDetectionImageFilter2.cxx index bd0a6a76..96892292 100644 --- a/src/itkCannyEdgeDetectionImageFilter2.cxx +++ b/src/itkCannyEdgeDetectionImageFilter2.cxx @@ -99,7 +99,6 @@ main(int argc, char * argv[]) cannyFilter->SetSigma(sigma); cannyFilter->SetUpperThreshold(upperThreshold); cannyFilter->SetLowerThreshold(lowerThreshold); - cannyFilter->SetOutsideValue(255); writer->SetInput(cannyFilter->GetOutput()); writer->UseCompressionOn(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f4f4b325..82c6242e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,6 +1,14 @@ itk_module_test() set(LesionSizingToolkitTests itkBinaryThresholdFeatureGeneratorTest1.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7.cxx +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8.cxx itkCannyEdgesDistanceAdvectionFieldFeatureGeneratorTest1.cxx itkCannyEdgesDistanceFeatureGeneratorTest1.cxx itkCannyEdgesFeatureGeneratorTest1.cxx @@ -265,6 +273,55 @@ itk_add_test(NAME itkCannyEdgesDistanceFeatureGeneratorTest1 75 # Lower hysteresis threshold ) +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1 + COMMAND LesionSizingToolkitTestDriver itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1 + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2 + COMMAND LesionSizingToolkitTestDriver + --compare ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest2_A.png + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest2_B.png + itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2 + ${TEST_DATA_ROOT}/Input/cthead1.png + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest2_A.png + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest2_B.png + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3 + COMMAND LesionSizingToolkitTestDriver + --compare ${TEST_DATA_ROOT}/Baseline/CannyEdgeDetectionRecursiveGaussianImageFilterTest3.png + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest3.png + itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3 + ${TEST_DATA_ROOT}/Input/cthead1.png + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest3.png + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4 + COMMAND LesionSizingToolkitTestDriver + --compare ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest4_A.mha + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest4_B.mha + itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4 + ${TEST_DATA_ROOT}/Input/PartSolidLesionCropped.mha + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest4_A.mha + ${TEMP}/CannyEdgeDetectionRecursiveGaussianImageFilterTest4_B.mha + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5 + COMMAND LesionSizingToolkitTestDriver itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5 + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6 + COMMAND LesionSizingToolkitTestDriver itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6 + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7 + COMMAND LesionSizingToolkitTestDriver itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7 + ) + +itk_add_test(NAME itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8 + COMMAND LesionSizingToolkitTestDriver itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8 + ) + itk_add_test(NAME itkCannyEdgesFeatureGeneratorTest1 COMMAND LesionSizingToolkitTestDriver itkCannyEdgesFeatureGeneratorTest1 ${TEST_DATA_ROOT}/Input/PartSolidLesionCropped.mha diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1.cxx new file mode 100644 index 00000000..104c56fd --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1.cxx @@ -0,0 +1,232 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkCommand.h" +#include "itkImage.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkTestingMacros.h" + +#include +#include + +namespace +{ + +template +typename itk::Image::Pointer +CreateBoxImage(itk::SizeValueType sideLength, itk::IndexValueType boxLower, itk::IndexValueType boxUpper) +{ + using ImageType = itk::Image; + typename ImageType::SizeType size; + size.Fill(sideLength); + + auto image = ImageType::New(); + image->SetRegions(typename ImageType::RegionType(size)); + image->Allocate(true); + + itk::ImageRegionIteratorWithIndex it(image, image->GetBufferedRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + const auto index = it.GetIndex(); + bool inside = true; + for (unsigned int d = 0; d < VDimension; ++d) + { + inside = inside && index[d] >= boxLower && index[d] < boxUpper; + } + it.Set(inside ? 100.0f : 0.0f); + } + return image; +} + +// True when the index lies within `band` voxels of the surface of the box. +template +bool +IsNearBoxSurface(const itk::Index & index, + itk::IndexValueType boxLower, + itk::IndexValueType boxUpper, + itk::IndexValueType band) +{ + bool withinOuter = true; + bool withinInner = true; + for (unsigned int d = 0; d < VDimension; ++d) + { + withinOuter = withinOuter && index[d] >= boxLower - band && index[d] < boxUpper + band; + withinInner = withinInner && index[d] >= boxLower + band && index[d] < boxUpper - band; + } + return withinOuter && !withinInner; +} + +class ProgressRecorder : public itk::Command +{ +public: + using Self = ProgressRecorder; + using Superclass = itk::Command; + using Pointer = itk::SmartPointer; + itkNewMacro(Self); + + void + Execute(itk::Object * caller, const itk::EventObject & event) override + { + this->Execute(static_cast(caller), event); + } + + void + Execute(const itk::Object * caller, const itk::EventObject &) override + { + const auto * process = dynamic_cast(caller); + m_Values.push_back(process->GetProgress()); + } + + std::vector m_Values; + +protected: + ProgressRecorder() = default; +}; + +template +int +RunTest(itk::SizeValueType sideLength, itk::IndexValueType boxLower, itk::IndexValueType boxUpper) +{ + using ImageType = itk::Image; + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + + constexpr itk::IndexValueType band = 3; + + auto input = CreateBoxImage(sideLength, boxLower, boxUpper); + + auto filter = FilterType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, CannyEdgeDetectionRecursiveGaussianImageFilter, ImageToImageFilter); + + filter->SetInput(input); + filter->SetSigma(1.0); + filter->SetUpperThreshold(20.0f); + filter->SetLowerThreshold(10.0f); + ITK_TEST_SET_GET_VALUE(20.0f, filter->GetUpperThreshold()); + ITK_TEST_SET_GET_VALUE(10.0f, filter->GetLowerThreshold()); + + auto progress = ProgressRecorder::New(); + filter->AddObserver(itk::ProgressEvent(), progress); + + // The result must not depend on the number of work units. + typename ImageType::Pointer reference; + for (const unsigned int workUnits : { 1u, 2u, 3u, 8u }) + { + filter->SetNumberOfWorkUnits(workUnits); + progress->m_Values.clear(); + filter->Modified(); + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + + if (progress->m_Values.empty()) + { + std::cerr << "No progress events with " << workUnits << " work units." << std::endl; + return EXIT_FAILURE; + } + const auto [minProgress, maxProgress] = std::minmax_element(progress->m_Values.begin(), progress->m_Values.end()); + if (*minProgress < 0.0f || *maxProgress > 1.0f) + { + std::cerr << "Progress outside [0, 1]: [" << *minProgress << ", " << *maxProgress << "]" << std::endl; + return EXIT_FAILURE; + } + if (*maxProgress < 0.9f) + { + std::cerr << "Progress never approached completion: max " << *maxProgress << std::endl; + return EXIT_FAILURE; + } + + typename ImageType::Pointer output = filter->GetOutput(); + output->DisconnectPipeline(); + + if (reference.IsNull()) + { + reference = output; + continue; + } + itk::ImageRegionConstIteratorWithIndex refIt(reference, reference->GetBufferedRegion()); + itk::ImageRegionConstIteratorWithIndex outIt(output, output->GetBufferedRegion()); + for (; !refIt.IsAtEnd(); ++refIt, ++outIt) + { + if (refIt.Get() != outIt.Get()) + { + std::cerr << "Output differs between 1 and " << workUnits << " work units at " << refIt.GetIndex() << std::endl; + return EXIT_FAILURE; + } + } + } + + // Edges are binary, present, and lie on the surface of the box. + itk::SizeValueType edgeCount = 0; + itk::ImageRegionConstIteratorWithIndex it(reference, reference->GetBufferedRegion()); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + const float value = it.Get(); + if (value != 0.0f && value != 1.0f) + { + std::cerr << "Non-binary edge value " << value << " at " << it.GetIndex() << std::endl; + return EXIT_FAILURE; + } + if (value == 1.0f) + { + ++edgeCount; + if (!IsNearBoxSurface(it.GetIndex(), boxLower, boxUpper, band)) + { + std::cerr << "Edge pixel away from the box surface at " << it.GetIndex() << std::endl; + return EXIT_FAILURE; + } + } + } + if (edgeCount == 0) + { + std::cerr << "No edge pixels detected." << std::endl; + return EXIT_FAILURE; + } + + // A smaller box must move the edges, which guards against a constant output. + auto smallerInput = CreateBoxImage(sideLength, boxLower + 4, boxUpper - 4); + filter->SetInput(smallerInput); + filter->SetNumberOfWorkUnits(4); + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + itk::ImageRegionConstIteratorWithIndex smallIt(filter->GetOutput(), + filter->GetOutput()->GetBufferedRegion()); + for (smallIt.GoToBegin(); !smallIt.IsAtEnd(); ++smallIt) + { + if (smallIt.Get() == 1.0f && !IsNearBoxSurface(smallIt.GetIndex(), boxLower + 4, boxUpper - 4, band)) + { + std::cerr << "Stale edge pixel after changing the input at " << smallIt.GetIndex() << std::endl; + return EXIT_FAILURE; + } + } + + return EXIT_SUCCESS; +} + +} // namespace + +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest1(int, char *[]) +{ + if (RunTest<2>(64, 16, 48) != EXIT_SUCCESS) + { + return EXIT_FAILURE; + } + if (RunTest<3>(40, 10, 30) != EXIT_SUCCESS) + { + return EXIT_FAILURE; + } + + std::cout << "Test finished." << std::endl; + return EXIT_SUCCESS; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2.cxx new file mode 100644 index 00000000..4e2b681c --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2.cxx @@ -0,0 +1,92 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkRescaleIntensityImageFilter.h" +#include "itkTestingMacros.h" + +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest2(int argc, char * argv[]) +{ + if (argc < 4) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << " InputImage OutputImage1 OutputImage2" + << std::endl; + return EXIT_FAILURE; + } + + constexpr unsigned int Dimension{ 2 }; + using InputImage = itk::Image; + using OutputImage = itk::Image; + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + + auto reader = itk::ImageFileReader::New(); + reader->SetFileName(argv[1]); + + auto filter = FilterType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, CannyEdgeDetectionRecursiveGaussianImageFilter, ImageToImageFilter); + + filter->SetInput(reader->GetOutput()); + + FilterType::OutputImagePixelType upperThreshold = 25; + filter->SetUpperThreshold(upperThreshold); + ITK_TEST_SET_GET_VALUE(upperThreshold, filter->GetUpperThreshold()); + + FilterType::OutputImagePixelType lowerThreshold = 10; + filter->SetLowerThreshold(lowerThreshold); + ITK_TEST_SET_GET_VALUE(lowerThreshold, filter->GetLowerThreshold()); + + constexpr FilterType::ScalarRealType sigma = 1.0; + filter->SetSigma(sigma); + ITK_TEST_SET_GET_VALUE(sigma, filter->GetSigma()); + + auto rescale = itk::RescaleIntensityImageFilter::New(); + rescale->SetInput(filter->GetOutput()); + rescale->SetOutputMinimum(0); + rescale->SetOutputMaximum(255); + + auto writer = itk::ImageFileWriter::New(); + writer->SetInput(rescale->GetOutput()); + writer->SetFileName(argv[2]); + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + // Different thresholds must re-execute the pipeline. + upperThreshold = 20; + filter->SetUpperThreshold(upperThreshold); + ITK_TEST_SET_GET_VALUE(upperThreshold, filter->GetUpperThreshold()); + + lowerThreshold = 5; + filter->SetLowerThreshold(lowerThreshold); + ITK_TEST_SET_GET_VALUE(lowerThreshold, filter->GetLowerThreshold()); + + ITK_TRY_EXPECT_NO_EXCEPTION(rescale->Update()); + + // Restoring the thresholds must reproduce the first result. + upperThreshold = 25; + filter->SetUpperThreshold(upperThreshold); + ITK_TEST_SET_GET_VALUE(upperThreshold, filter->GetUpperThreshold()); + + lowerThreshold = 10; + filter->SetLowerThreshold(lowerThreshold); + ITK_TEST_SET_GET_VALUE(lowerThreshold, filter->GetLowerThreshold()); + + writer->SetFileName(argv[3]); + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + return EXIT_SUCCESS; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3.cxx new file mode 100644 index 00000000..b0c4d7e9 --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3.cxx @@ -0,0 +1,67 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkRescaleIntensityImageFilter.h" +#include "itkSimpleFilterWatcher.h" +#include "itkTestingMacros.h" + +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest3(int argc, char * argv[]) +{ + if (argc < 3) + { + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << " InputImage OutputImage" << std::endl; + return EXIT_FAILURE; + } + + constexpr unsigned int Dimension{ 2 }; + using InputImage = itk::Image; + using OutputImage = itk::Image; + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + + auto reader = itk::ImageFileReader::New(); + reader->SetFileName(argv[1]); + + auto filter = FilterType::New(); + const itk::SimpleFilterWatcher watcher(filter); + filter->SetInput(reader->GetOutput()); + + constexpr FilterType::OutputImagePixelType upperThreshold{ 30 }; + filter->SetUpperThreshold(upperThreshold); + ITK_TEST_SET_GET_VALUE(upperThreshold, filter->GetUpperThreshold()); + + constexpr FilterType::OutputImagePixelType lowerThreshold{ 15 }; + filter->SetLowerThreshold(lowerThreshold); + ITK_TEST_SET_GET_VALUE(lowerThreshold, filter->GetLowerThreshold()); + + const FilterType::SigmaArrayType sigma(1.0); + filter->SetSigmaArray(sigma); + ITK_TEST_SET_GET_VALUE(sigma, filter->GetSigmaArray()); + + auto rescale = itk::RescaleIntensityImageFilter::New(); + rescale->SetInput(filter->GetOutput()); + rescale->SetOutputMinimum(0); + rescale->SetOutputMaximum(255); + + auto writer = itk::ImageFileWriter::New(); + writer->SetInput(rescale->GetOutput()); + writer->SetFileName(argv[2]); + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + return EXIT_SUCCESS; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4.cxx new file mode 100644 index 00000000..10ef1021 --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4.cxx @@ -0,0 +1,70 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImageFileReader.h" +#include "itkImageFileWriter.h" +#include "itkCastImageFilter.h" +#include "itkTestingMacros.h" + +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest4(int argc, char * argv[]) +{ + if (argc < 4) + { + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) + << " InputImage OutputImage1 OutputImage2 [sigma upperThreshold lowerThreshold]" << std::endl; + return EXIT_FAILURE; + } + + constexpr unsigned int Dimension{ 3 }; + using InputImage = itk::Image; + using InternalImage = itk::Image; + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + + const double sigma = argc > 4 ? std::stod(argv[4]) : 0.7; + const float upperThreshold = argc > 5 ? std::stof(argv[5]) : 150.0f; + const float lowerThreshold = argc > 6 ? std::stof(argv[6]) : 75.0f; + + auto reader = itk::ImageFileReader::New(); + reader->SetFileName(argv[1]); + + auto cast = itk::CastImageFilter::New(); + cast->SetInput(reader->GetOutput()); + + auto filter = FilterType::New(); + filter->SetInput(cast->GetOutput()); + filter->SetSigma(sigma); + filter->SetUpperThreshold(upperThreshold); + filter->SetLowerThreshold(lowerThreshold); + + auto writer = itk::ImageFileWriter::New(); + writer->SetInput(filter->GetOutput()); + writer->SetFileName(argv[2]); + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + // Looser thresholds must produce a different, re-executed result. + filter->SetUpperThreshold(upperThreshold / 2.0f); + filter->SetLowerThreshold(lowerThreshold / 2.0f); + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + + // Restoring the thresholds must reproduce the first result. + filter->SetUpperThreshold(upperThreshold); + filter->SetLowerThreshold(lowerThreshold); + writer->SetFileName(argv[3]); + ITK_TRY_EXPECT_NO_EXCEPTION(writer->Update()); + + return EXIT_SUCCESS; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5.cxx new file mode 100644 index 00000000..59b45e02 --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5.cxx @@ -0,0 +1,91 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImage.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkTestingMacros.h" + +namespace +{ +using ImageType = itk::Image; + +// A box whose edge strength varies along its contour, so that hysteresis +// linking decides which contour pixels survive. +ImageType::Pointer +CreateImage(const ImageType::IndexType & start) +{ + ImageType::SizeType size; + size.Fill(64); + + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType(start, size)); + image->Allocate(true); + + itk::ImageRegionIteratorWithIndex it(image, image->GetBufferedRegion()); + for (; !it.IsAtEnd(); ++it) + { + const auto offset = it.GetIndex() - start; + const bool inside = offset[0] >= 16 && offset[0] < 48 && offset[1] >= 16 && offset[1] < 48; + it.Set(inside ? 10.0f + 90.0f * offset[0] / 64.0f + 40.0f * offset[1] / 64.0f : 0.0f); + } + return image; +} + +ImageType::Pointer +Detect(const ImageType::IndexType & start) +{ + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + auto filter = FilterType::New(); + filter->SetInput(CreateImage(start)); + filter->SetSigma(1.0); + filter->SetUpperThreshold(32.0f); + filter->SetLowerThreshold(8.0f); + filter->Update(); + ImageType::Pointer output = filter->GetOutput(); + output->DisconnectPipeline(); + return output; +} +} // namespace + +// The detected edges must not depend on where the image region starts. +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest5(int, char *[]) +{ + ImageType::IndexType origin; + origin.Fill(0); + const auto reference = Detect(origin); + + int status = EXIT_SUCCESS; + for (const auto & start : { ImageType::IndexType{ { 10, -7 } }, ImageType::IndexType{ { -30, 40 } } }) + { + const auto shifted = Detect(start); + itk::ImageRegionConstIterator refIt(reference, reference->GetBufferedRegion()); + itk::ImageRegionConstIterator shiftedIt(shifted, shifted->GetBufferedRegion()); + itk::SizeValueType differing = 0; + for (; !refIt.IsAtEnd(); ++refIt, ++shiftedIt) + { + differing += refIt.Get() != shiftedIt.Get() ? 1 : 0; + } + if (differing != 0) + { + std::cerr << "Region starting at " << start << " differs from the origin-based result in " << differing + << " pixels." << std::endl; + status = EXIT_FAILURE; + } + } + return status; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6.cxx new file mode 100644 index 00000000..a2a3d6a0 --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6.cxx @@ -0,0 +1,171 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImage.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkTestingMacros.h" + +#include +#include + +namespace +{ +using ImageType = itk::Image; +using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + +ImageType::Pointer +CreateImage(itk::SizeValueType width, itk::SizeValueType height) +{ + ImageType::SizeType size{ { width, height } }; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType(size)); + image->Allocate(true); + + itk::ImageRegionIteratorWithIndex it(image, image->GetBufferedRegion()); + for (; !it.IsAtEnd(); ++it) + { + const auto index = it.GetIndex(); + const bool inside = index[0] >= 12 && index[0] < 40 && index[1] >= 12 && index[1] < 40; + it.Set(inside ? 10.0f + 2.0f * index[0] + index[1] : 0.0f); + } + return image; +} + +FilterType::Pointer +Run(const ImageType * image, float upper, float lower, const FilterType::SigmaArrayType & sigma) +{ + auto filter = FilterType::New(); + filter->SetInput(image); + filter->SetSigmaArray(sigma); + filter->SetUpperThreshold(upper); + filter->SetLowerThreshold(lower); + filter->Update(); + return filter; +} + +itk::SizeValueType +CountEdges(const ImageType * image, float upper, float lower, const FilterType::SigmaArrayType & sigma) +{ + const auto filter = Run(image, upper, lower, sigma); + + itk::SizeValueType count = 0; + itk::ImageRegionConstIterator it(filter->GetOutput(), filter->GetOutput()->GetBufferedRegion()); + for (; !it.IsAtEnd(); ++it) + { + count += it.Get() == 1.0f ? 1 : 0; + } + return count; +} + +// Sorted nonzero responses of the non-maximum suppression stage. +std::vector +Responses(const ImageType * image, const FilterType::SigmaArrayType & sigma) +{ + const auto filter = Run(image, 0.0f, 0.0f, sigma); + const ImageType * nms = filter->GetNonMaximumSuppressionImage(); + std::vector values; + itk::ImageRegionConstIterator it(nms, nms->GetBufferedRegion()); + for (; !it.IsAtEnd(); ++it) + { + if (it.Get() > 0.0f) + { + values.push_back(it.Get()); + } + } + std::sort(values.begin(), values.end()); + return values; +} +} // namespace + +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest6(int, char *[]) +{ + auto filter = FilterType::New(); + ITK_EXERCISE_BASIC_OBJECT_METHODS(filter, CannyEdgeDetectionRecursiveGaussianImageFilter, ImageToImageFilter); + + // Defaults. + ITK_TEST_SET_GET_VALUE(1.0, filter->GetSigma()); + ITK_TEST_SET_GET_VALUE(0.0f, filter->GetUpperThreshold()); + ITK_TEST_SET_GET_VALUE(0.0f, filter->GetLowerThreshold()); + + // Sigma accessors. + filter->SetSigma(2.0); + FilterType::SigmaArrayType expectedSigma(2.0); + ITK_TEST_SET_GET_VALUE(expectedSigma, filter->GetSigmaArray()); + + FilterType::SigmaArrayType anisotropic; + anisotropic[0] = 0.8; + anisotropic[1] = 1.6; + filter->SetSigmaArray(anisotropic); + ITK_TEST_SET_GET_VALUE(anisotropic, filter->GetSigmaArray()); + ITK_TEST_SET_GET_VALUE(0.8, filter->GetSigma()); + + const auto modifiedTime = filter->GetMTime(); + filter->SetSigmaArray(anisotropic); + ITK_TEST_SET_GET_VALUE(modifiedTime, filter->GetMTime()); + + + // No input. + ITK_TRY_EXPECT_EXCEPTION(filter->Update()); + + // The recursive Gaussian needs at least four pixels along each axis. + filter = FilterType::New(); + auto tiny = CreateImage(3, 40); + filter->SetInput(tiny); + ITK_TRY_EXPECT_EXCEPTION(filter->Update()); + + // Output matches the input geometry, and the non-maximum suppression image is available. + auto image = CreateImage(52, 44); + filter = FilterType::New(); + filter->SetInput(image); + filter->SetUpperThreshold(20.0f); + filter->SetLowerThreshold(10.0f); + ITK_TRY_EXPECT_NO_EXCEPTION(filter->Update()); + ITK_TEST_SET_GET_VALUE(image->GetLargestPossibleRegion(), filter->GetOutput()->GetLargestPossibleRegion()); + ITK_TEST_EXPECT_TRUE(filter->GetNonMaximumSuppressionImage() != nullptr); + ITK_TEST_SET_GET_VALUE(image->GetLargestPossibleRegion(), + filter->GetNonMaximumSuppressionImage()->GetLargestPossibleRegion()); + filter->Print(std::cout); + + // Strong seeds only, then the same seeds linked through weaker responses. + const FilterType::SigmaArrayType isotropic(1.0); + const auto responses = Responses(image, isotropic); + ITK_TEST_EXPECT_TRUE(responses.size() > 20); + const float weak = responses[responses.size() / 4]; + const float strong = responses[responses.size() * 9 / 10]; + + const auto all = CountEdges(image, weak, weak, isotropic); + const auto seedsOnly = CountEdges(image, strong, strong, isotropic); + const auto linked = CountEdges(image, strong, weak, isotropic); + const auto linkedMore = CountEdges(image, strong, 0.0f, isotropic); + const auto none = CountEdges(image, responses.back() + 1.0f, 0.0f, isotropic); + std::cout << "edges: all=" << all << " seedsOnly=" << seedsOnly << " linked=" << linked + << " linkedMore=" << linkedMore << " none=" << none << std::endl; + ITK_TEST_EXPECT_TRUE(seedsOnly > 0); + ITK_TEST_EXPECT_TRUE(linked > seedsOnly); + ITK_TEST_EXPECT_TRUE(linkedMore >= linked); + ITK_TEST_EXPECT_TRUE(all >= linked); + ITK_TEST_EXPECT_TRUE(none == 0); + + // Anisotropic smoothing changes the response. + const auto isotropicResponses = Responses(image, isotropic); + const auto anisotropicResponses = Responses(image, anisotropic); + ITK_TEST_EXPECT_TRUE(isotropicResponses != anisotropicResponses); + + std::cout << "Test finished." << std::endl; + return EXIT_SUCCESS; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7.cxx new file mode 100644 index 00000000..dab2ca42 --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7.cxx @@ -0,0 +1,55 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImage.h" +#include "itkTestingMacros.h" + +namespace +{ +using ImageType = itk::Image; + +ImageType::Pointer +CreateImage(itk::SizeValueType width, itk::SizeValueType height) +{ + ImageType::SizeType size{ { width, height } }; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType(size)); + image->Allocate(true); + return image; +} +} // namespace + +// One filter instance must handle successive inputs of different sizes. +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest7(int, char *[]) +{ + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + auto filter = FilterType::New(); + + for (const auto & size : { ImageType::SizeType{ { 52, 44 } }, + ImageType::SizeType{ { 70, 60 } }, + ImageType::SizeType{ { 30, 36 } }, + ImageType::SizeType{ { 52, 44 } } }) + { + auto image = CreateImage(size[0], size[1]); + filter->SetInput(image); + ITK_TRY_EXPECT_NO_EXCEPTION(filter->UpdateLargestPossibleRegion()); + ITK_TEST_SET_GET_VALUE(image->GetLargestPossibleRegion(), filter->GetOutput()->GetLargestPossibleRegion()); + } + + std::cout << "Test finished." << std::endl; + return EXIT_SUCCESS; +} diff --git a/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8.cxx b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8.cxx new file mode 100644 index 00000000..31609c2d --- /dev/null +++ b/test/itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8.cxx @@ -0,0 +1,90 @@ +/*========================================================================= + + Program: Lesion Sizing Toolkit + Module: itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8.cxx + + Copyright (c) Kitware Inc. + All rights reserved. + See Copyright.txt or https://www.kitware.com/Copyright.htm for details. + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the above copyright notice for more information. + +=========================================================================*/ + +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkAddImageFilter.h" +#include "itkCannyEdgeDetectionRecursiveGaussianImageFilter.h" +#include "itkImage.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkTestingMacros.h" + +namespace +{ +using ImageType = itk::Image; + +ImageType::Pointer +CreateImage() +{ + ImageType::SizeType size{ { 48, 40 } }; + auto image = ImageType::New(); + image->SetRegions(ImageType::RegionType(size)); + image->Allocate(true); + + itk::ImageRegionIteratorWithIndex it(image, image->GetBufferedRegion()); + for (; !it.IsAtEnd(); ++it) + { + const auto index = it.GetIndex(); + const bool inside = index[0] > 12 && index[0] < 36 && index[1] > 10 && index[1] < 30; + it.Set(inside ? 80.0f : 0.0f); + } + return image; +} + +double +Sum(const ImageType * image) +{ + double sum = 0.0; + itk::ImageRegionConstIterator it(image, image->GetBufferedRegion()); + for (; !it.IsAtEnd(); ++it) + { + sum += it.Get(); + } + return sum; +} +} // namespace + +// The internal sub-filters must not modify the output's pipeline state, which +// breaks a downstream consumer when the output may be released after use. +int +itkCannyEdgeDetectionRecursiveGaussianImageFilterTest8(int, char *[]) +{ + using FilterType = itk::CannyEdgeDetectionRecursiveGaussianImageFilter; + auto filter = FilterType::New(); + filter->SetInput(CreateImage()); + filter->SetSigma(1.0); + filter->SetUpperThreshold(20.0f); + filter->SetLowerThreshold(10.0f); + filter->ReleaseDataFlagOn(); + + auto add = itk::AddImageFilter::New(); + add->SetInput1(filter->GetOutput()); + add->SetInput2(filter->GetOutput()); + + double sums[3] = { 0.0, 0.0, 0.0 }; + for (double & sum : sums) + { + add->Modified(); + ITK_TRY_EXPECT_NO_EXCEPTION(add->Update()); + sum = Sum(add->GetOutput()); + } + + ITK_TEST_EXPECT_TRUE(sums[0] > 0.0); + ITK_TEST_EXPECT_TRUE(sums[0] == sums[1]); + ITK_TEST_EXPECT_TRUE(sums[1] == sums[2]); + + std::cout << "Test finished." << std::endl; + return EXIT_SUCCESS; +}