1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkNarrowBandLevelSetImageFilter_hxx
19 #define itkNarrowBandLevelSetImageFilter_hxx
20
21 #include "itkNarrowBandLevelSetImageFilter.h"
22 #include <cstdio>
23 #include "itkMath.h"
24
25 namespace itk
26 {
27 template< typename TInputImage, typename TFeatureImage, typename TOutputPixelType, typename TOutputImage >
28 void
29 NarrowBandLevelSetImageFilter< TInputImage, TFeatureImage, TOutputPixelType, TOutputImage >
PrintSelf(std::ostream & os,Indent indent) const30 ::PrintSelf(std::ostream & os, Indent indent) const
31 {
32 Superclass::PrintSelf(os, indent);
33 os << indent << "m_ReverseExpansionDirection = " << m_ReverseExpansionDirection << std::endl;
34 os << indent << "m_SegmentationFunction = " << m_SegmentationFunction << std::endl;
35 }
36
37 template< typename TInputImage, typename TFeatureImage, typename TOutputPixelType, typename TOutputImage >
38 NarrowBandLevelSetImageFilter< TInputImage, TFeatureImage, TOutputPixelType, TOutputImage >
NarrowBandLevelSetImageFilter()39 ::NarrowBandLevelSetImageFilter()
40 {
41 this->SetNumberOfRequiredInputs(2);
42 //this->SetNarrowBandInnerRadius();
43 //this->SetNarrowBandTotalRadius();
44 m_SegmentationFunction = nullptr;
45
46 m_IsoFilter = IsoFilterType::New();
47 m_ChamferFilter = ChamferFilterType::New();
48
49 // Provide some reasonable defaults which will at least prevent infinite
50 // looping.
51 this->SetMaximumRMSError(0.02);
52 this->SetNumberOfIterations(1000);
53 m_ReverseExpansionDirection = false;
54 }
55
56 template< typename TInputImage, typename TFeatureImage, typename TOutputPixelType, typename TOutputImage >
57 void
58 NarrowBandLevelSetImageFilter< TInputImage, TFeatureImage, TOutputPixelType, TOutputImage >
SetSegmentationFunction(SegmentationFunctionType * s)59 ::SetSegmentationFunction(SegmentationFunctionType *s)
60 {
61 unsigned int i;
62
63 m_SegmentationFunction = s;
64
65 typename SegmentationFunctionType::RadiusType r;
66 for ( i = 0; i < Self::ImageDimension; ++i )
67 {
68 r[i] = 1;
69 }
70
71 m_SegmentationFunction->Initialize(r);
72 this->SetDifferenceFunction(m_SegmentationFunction);
73 this->Modified();
74 }
75
76 template< typename TInputImage, typename TFeatureImage, typename TOutputPixelType, typename TOutputImage >
77 void
78 NarrowBandLevelSetImageFilter< TInputImage, TFeatureImage, TOutputPixelType, TOutputImage >
GenerateData()79 ::GenerateData()
80 {
81 if ( m_SegmentationFunction == nullptr )
82 { itkExceptionMacro("No finite difference function was specified."); }
83
84 // A positive speed value causes surface expansion, the opposite of the
85 // default. Flip the sign of the propagation and advection weights.
86 if ( m_ReverseExpansionDirection == true )
87 {
88 this->GetSegmentationFunction()->ReverseExpansionDirection();
89 }
90
91 // Allocate the images from which speeds will be sampled.
92 if ( Math::NotExactlyEquals(this->GetSegmentationFunction()->GetPropagationWeight(), 0) )
93 {
94 m_SegmentationFunction->AllocateSpeedImage();
95 m_SegmentationFunction->CalculateSpeedImage();
96 }
97 if ( Math::NotExactlyEquals(this->GetSegmentationFunction()->GetAdvectionWeight(), 0) )
98 {
99 m_SegmentationFunction->AllocateAdvectionImage();
100 m_SegmentationFunction->CalculateAdvectionImage();
101 }
102
103 // Start the solver
104 Superclass::GenerateData();
105
106 // Reset all the signs of the weights.
107 if ( m_ReverseExpansionDirection == true )
108 {
109 this->GetSegmentationFunction()->ReverseExpansionDirection();
110 }
111 }
112
113 template< typename TInputImage, typename TFeatureImage, typename TOutputPixelType, typename TOutputImage >
114 void
115 NarrowBandLevelSetImageFilter< TInputImage, TFeatureImage, TOutputPixelType, TOutputImage >
CreateNarrowBand()116 ::CreateNarrowBand()
117 {
118 typename OutputImageType::Pointer levelset = this->GetOutput();
119
120 if ( !this->m_NarrowBand->Empty() )
121 {
122 m_IsoFilter->SetNarrowBand( this->m_NarrowBand.GetPointer() );
123 m_IsoFilter->NarrowBandingOn(); //Maybe we should check that the NarrowBand
124 // exits.
125 }
126 else
127 {
128 m_IsoFilter->NarrowBandingOff();
129 }
130
131 m_IsoFilter->SetFarValue(this->m_NarrowBand->GetTotalRadius() + 1);
132 m_IsoFilter->SetInput(levelset);
133 m_IsoFilter->Update();
134
135 m_ChamferFilter->SetInput( m_IsoFilter->GetOutput() );
136 m_ChamferFilter->SetMaximumDistance(this->m_NarrowBand->GetTotalRadius() + 1);
137 m_ChamferFilter->SetNarrowBand( this->m_NarrowBand.GetPointer() );
138 m_ChamferFilter->Update();
139
140 this->GraftOutput( m_ChamferFilter->GetOutput() );
141 m_IsoFilter->SetInput(nullptr);
142 m_ChamferFilter->SetInput(nullptr);
143 }
144 } // end namespace itk
145
146 #endif
147