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_h 19 #define itkNarrowBandLevelSetImageFilter_h 20 21 #include "itkNarrowBandImageFilterBase.h" 22 #include "itkSegmentationLevelSetFunction.h" 23 #include "itkFastChamferDistanceImageFilter.h" 24 #include "itkIsoContourDistanceImageFilter.h" 25 #include "itkMath.h" 26 27 namespace itk 28 { 29 /** 30 * 31 * \class NarrowBandLevelSetImageFilter 32 * 33 * \brief A base class which defines the API for implementing a special class of 34 * image segmentation filters using level set methods. 35 * 36 * \par OVERVIEW 37 * This object defines the framework for a class of segmentation filters which 38 * use level set methods. These filters work by constructing a "feature image" 39 * onto which the evolving level set locks as it moves. In the feature image, 40 * values that are close to zero are associated with object boundaries. An 41 * original (or preprocessed) image is given to the filter as the feature image 42 * and a seed for the level set is given as the input of the filter. The seed is 43 * converted into a level set embedding which propagates according to the features 44 * calculated from the original image. 45 * 46 * \par TEMPLATE PARAMETERS 47 * There are two required and two optional template parameter for these 48 * filters. Of the optional parameters, the last, TOutputImage, should not be 49 * changed from its default. It is only there to instantiate the parent class 50 * correctly. 51 * 52 * TInputImage is the image type of the initial model you will input to the 53 * filter using SetInput() or SetInitialImage(). 54 * 55 * TFeatureImage is the image type of the image from which the filter will 56 * calculate the speed term for segmentation (see INPUTS). 57 * 58 * TOutputPixelType is the data type used for the output image phi, the implicit 59 * level set image. This should really only ever be set as float (default) or 60 * double. 61 * 62 * \par INPUTS 63 * The input to any subclass of this filter is the seed image for the initial 64 * level set embedding. As with other subclasses of the 65 * SparseLevelSetImageFilter, the type of the input image is is not important. 66 * The (RequestedRegion) size of the seed image must, however, match the 67 * (RequestedRegion) size of the feature image. 68 * 69 * You must identify the initial front (surface) in the input image. You do 70 * this by specifying its isovalue through the method SetIsoSurfaceValue(float 71 * f). The algorithm will then initialize its solution using the front represented by 72 * value f. Note that the front is always represented by isosurface zero in 73 * the output and not the isosurface you specified for the input. This is 74 * because, for simplicity, the filter will shift your input image so that the 75 * active front has zero values. 76 * 77 * \par 78 * Depending on the particular application and filter that you are using, the 79 * feature image should be preprocessed with some type of noise reduction 80 * filtering. The feature image input can be of any type, but it will be cast 81 * to floating point before calculations are done. 82 * 83 * \par OUTPUTS 84 * The output of any subclass of this filter is a level set embedding as 85 * described in SparseFieldLevelSetImageFilter. The zero crossings of the output 86 * image give the pixels closest to the level set boundary. By ITK convention, 87 * NEGATIVE values are pixels INSIDE the segmented region and POSITIVE values are 88 * pixels OUTSIDE the segmented region. 89 * 90 * \par PARAMETERS 91 * The MaximumRMSChange parameter is used to determine when the solution has 92 * converged. A lower value will result in a tighter-fitting solution, but 93 * will require more computations. Too low a value could put the solver into 94 * an infinite loop unless a reasonable NumberOfIterations parameter is set. 95 * Values should always be greater than 0.0 and less than 1.0. 96 * 97 * \par 98 * The NumberOfIterations parameter can be used to halt the solution after a 99 * specified number of iterations, overriding the MaximumRMSChange halting 100 * criteria. 101 * 102 * \par 103 * The standard convention for ITK level-set segmentation filters is that 104 * POSITIVE propagation (speed) and advection terms cause the surface to EXPAND 105 * while negative terms cause the surface to CONTRACT. When the 106 * ReverseExpansionDirection parameter is set to TRUE (on), it tells the 107 * function object to reverse the standard ITK convention so that NEGATIVE 108 * terms cause EXPANSION and positive terms cause CONTRACTION. 109 * 110 * This parameter can be safely changed as appropriate for a particular 111 * application or data set to achieve the desired behavior. 112 * 113 * \par 114 * The FeatureScaling parameter controls the magnitude of the features calculated 115 * for use in the level set propagation and advection speeds. This value simply 116 * sets both parameters to equal values at once. Some filters may only use on of 117 * these two terms and this method is a generic way to set either or both without 118 * having to know which is in use. 119 * 120 * \par 121 * The CurvatureScaling parameter controls the magnitude of the curvature values 122 * which are calculated on the evolving isophote. This is important in 123 * controlling the relative effect of curvature in the calculation. Default 124 * value is 1.0. Higher values relative to the other level set equation terms 125 * (propagation and advection) will give a smoother result. 126 * 127 * \par 128 * The PropagationScaling parameter controls the scaling of the scalar 129 * propagation (speed) term relative to other terms in the level set 130 * equation. Setting this value will override any value already set by 131 * FeatureScaling. 132 * 133 * \par 134 * The AdvectionScaling parameter controls the scaling of the vector 135 * advection field term relative to other terms in the level set 136 * equation. Setting this value will override any value already set by 137 * FeatureScaling. 138 * 139 * \par 140 * See LevelSetFunction for more information. 141 * \ingroup ITKLevelSets 142 */ 143 template< typename TInputImage, 144 typename TFeatureImage, 145 typename TOutputPixelType = float, 146 typename TOutputImage = Image< TOutputPixelType, 147 TInputImage::ImageDimension > > 148 class ITK_TEMPLATE_EXPORT NarrowBandLevelSetImageFilter: 149 public NarrowBandImageFilterBase< TInputImage, TOutputImage > 150 { 151 public: 152 ITK_DISALLOW_COPY_AND_ASSIGN(NarrowBandLevelSetImageFilter); 153 154 /** Standard class type aliases */ 155 using Self = NarrowBandLevelSetImageFilter; 156 using Superclass = NarrowBandImageFilterBase< TInputImage, TOutputImage >; 157 using Pointer = SmartPointer< Self >; 158 using ConstPointer = SmartPointer< const Self >; 159 160 /** Inherited type alias from the superclass. */ 161 using ValueType = typename Superclass::ValueType; 162 using IndexType = typename Superclass::IndexType; 163 using TimeStepType = typename Superclass::TimeStepType; 164 using InputImageType = typename Superclass::InputImageType; 165 166 /** Local image type alias */ 167 using OutputImageType = TOutputImage; 168 using FeatureImageType = TFeatureImage; 169 170 /** The generic level set function type */ 171 using SegmentationFunctionType = 172 SegmentationLevelSetFunction< OutputImageType, FeatureImageType >; 173 174 /** The type used for the advection field */ 175 using VectorImageType = typename SegmentationFunctionType::VectorImageType; 176 177 /** Run-time type information (and related methods). */ 178 itkTypeMacro(NarrowBandLevelSetImageFilter, NarrowBandImageFilterBase); 179 180 /** Set/Get the feature image to be used for speed function of the level set 181 * equation. Equivalent to calling Set/GetInput(1, ..) */ SetFeatureImage(const FeatureImageType * f)182 virtual void SetFeatureImage(const FeatureImageType *f) 183 { 184 this->ProcessObject::SetNthInput( 1, const_cast< FeatureImageType * >( f ) ); 185 m_SegmentationFunction->SetFeatureImage(f); 186 } 187 GetFeatureImage()188 virtual FeatureImageType * GetFeatureImage() 189 { 190 return ( static_cast< FeatureImageType * >( this->ProcessObject::GetInput(1) ) ); 191 } 192 193 /** Set/Get the initial level set model. Equivalent to calling SetInput(..) 194 */ SetInitialImage(InputImageType * f)195 virtual void SetInitialImage(InputImageType *f) 196 { 197 this->SetInput(f); 198 } 199 GetSpeedImage()200 virtual const typename SegmentationFunctionType::ImageType * GetSpeedImage() const 201 { return m_SegmentationFunction->GetSpeedImage(); } 202 GetAdvectionImage()203 virtual const typename SegmentationFunctionType::VectorImageType * GetAdvectionImage() const 204 { return m_SegmentationFunction->GetAdvectionImage(); } 205 206 /** THIS METHOD IS DEPRECATED AND SHOULD NOT BE USED. This method reverses 207 * the speed function direction, effectively changing inside feature values to 208 * outside feature values and vice versa. */ SetUseNegativeFeaturesOn()209 void SetUseNegativeFeaturesOn() 210 { 211 itkWarningMacro( 212 << "SetUseNegativeFeaturesOn has been deprecated. Please use ReverseExpansionDirectionOn() instead"); 213 this->ReverseExpansionDirectionOn(); 214 } 215 SetUseNegativeFeaturesOff()216 void SetUseNegativeFeaturesOff() 217 { 218 itkWarningMacro( 219 << "SetUseNegativeFeaturesOff has been deprecated. Please use ReverseExpansionDirectionOff() instead"); 220 this->ReverseExpansionDirectionOff(); 221 } 222 223 /** Set/Get the value of the UseNegativeFeatures flag. This method is 224 * deprecated. Use Set/Get ReverseExpansionDirection instead. */ SetUseNegativeFeatures(bool u)225 void SetUseNegativeFeatures(bool u) 226 { 227 itkWarningMacro(<< "SetUseNegativeFeatures has been deprecated. Please use SetReverseExpansionDirection instead"); 228 if ( u == true ) 229 { 230 this->SetReverseExpansionDirection(false); 231 } 232 else 233 { 234 this->SetReverseExpansionDirection(true); 235 } 236 } 237 GetUseNegativeFeatures()238 bool GetUseNegativeFeatures() const 239 { 240 itkWarningMacro(<< "GetUseNegativeFeatures has been deprecated. Please use GetReverseExpansionDirection() instead"); 241 if ( this->GetReverseExpansionDirection() == false ) 242 { 243 return true; 244 } 245 else 246 { 247 return false; 248 } 249 } 250 251 /** Turn On/Off the flag which determines whether Positive or Negative speed 252 * terms will cause surface expansion. If set to TRUE then negative speed 253 * terms will cause the surface to expand and positive speed terms will cause 254 * the surface to contract. If set to FALSE (default) then positive speed terms will 255 * cause the surface to expand and negative speed terms will cause the 256 * surface to contract. This method can be safely used to reverse the 257 * expansion/contraction as appropriate to a particular application or data 258 * set. */ 259 itkSetMacro(ReverseExpansionDirection, bool); 260 itkGetConstMacro(ReverseExpansionDirection, bool); 261 itkBooleanMacro(ReverseExpansionDirection); 262 263 /** Combined scaling of the propagation and advection speed 264 terms. You should use either this -or- Get/SetPropagationScaling and 265 Get/SetAdvectionScaling (if appropriate). See subclasses for details 266 on when and whether to set these parameters. */ SetFeatureScaling(ValueType v)267 void SetFeatureScaling(ValueType v) 268 { 269 if ( v != m_SegmentationFunction->GetPropagationWeight() ) 270 { 271 this->SetPropagationScaling(v); 272 } 273 if ( v != m_SegmentationFunction->GetAdvectionWeight() ) 274 { 275 this->SetAdvectionScaling(v); 276 } 277 } 278 279 /** Set/Get the scaling of the propagation speed. Setting the FeatureScaling 280 parameter overrides any previous values set for PropagationScaling. */ SetPropagationScaling(ValueType v)281 void SetPropagationScaling(ValueType v) 282 { 283 if ( Math::NotExactlyEquals(v, m_SegmentationFunction->GetPropagationWeight()) ) 284 { 285 m_SegmentationFunction->SetPropagationWeight(v); 286 } 287 } 288 GetPropagationScaling()289 ValueType GetPropagationScaling() const 290 { 291 return m_SegmentationFunction->GetPropagationWeight(); 292 } 293 294 /** Set/Get the scaling of the advection field. Setting the FeatureScaling 295 parameter will override any existing value for AdvectionScaling. */ SetAdvectionScaling(ValueType v)296 void SetAdvectionScaling(ValueType v) 297 { 298 if ( Math::NotExactlyEquals(v, m_SegmentationFunction->GetAdvectionWeight()) ) 299 { 300 m_SegmentationFunction->SetAdvectionWeight(v); 301 } 302 } 303 GetAdvectionScaling()304 ValueType GetAdvectionScaling() const 305 { 306 return m_SegmentationFunction->GetAdvectionWeight(); 307 } 308 309 /** Set/Get the scaling of the curvature. Use this parameter to increase the 310 * influence of curvature on the movement of the surface. Higher 311 * values relative to Advection and Propagation values will give 312 * smoother surfaces. */ SetCurvatureScaling(ValueType v)313 void SetCurvatureScaling(ValueType v) 314 { 315 if ( Math::NotExactlyEquals(v, m_SegmentationFunction->GetCurvatureWeight()) ) 316 { 317 m_SegmentationFunction->SetCurvatureWeight(v); 318 } 319 } 320 GetCurvatureScaling()321 ValueType GetCurvatureScaling() const 322 { 323 return m_SegmentationFunction->GetCurvatureWeight(); 324 } 325 326 /** Set the segmentation function. In general, this should only be called by a subclass 327 * of this object. It is made public to allow itk::Command objects access. */ 328 virtual void SetSegmentationFunction(SegmentationFunctionType *s); 329 GetSegmentationFunction()330 virtual SegmentationFunctionType * GetSegmentationFunction() 331 { return m_SegmentationFunction; } 332 333 /** Set/Get the maximum number of iterations allowed for the solver. This 334 * prevents infinite loops if a solution "bounces". */ SetMaximumIterations(unsigned int i)335 void SetMaximumIterations(unsigned int i) 336 { 337 itkWarningMacro("SetMaximumIterations is deprecated. Please use SetNumberOfIterations instead."); 338 this->SetNumberOfIterations(i); 339 } 340 GetMaximumIterations()341 unsigned int GetMaximumIterations() 342 { 343 itkWarningMacro("GetMaximumIterations is deprecated. Please use GetNumberOfIterations instead."); 344 return this->GetNumberOfIterations(); 345 } 346 SetMaximumRMSError(const double)347 void SetMaximumRMSError(const double) override 348 { 349 itkWarningMacro( 350 "The current implmentation of this solver does not compute maximum RMS change. The maximum RMS error value will not be set or used."); 351 } 352 353 #ifdef ITK_USE_CONCEPT_CHECKING 354 // Begin concept checking 355 itkConceptMacro( OutputHasNumericTraitsCheck, 356 ( Concept::HasNumericTraits< typename TOutputImage::PixelType > ) ); 357 // End concept checking 358 #endif 359 360 protected: 361 ~NarrowBandLevelSetImageFilter() override = default; 362 NarrowBandLevelSetImageFilter(); 363 364 void PrintSelf(std::ostream & os, Indent indent) const override; 365 366 367 /** Overrides parent implementation */ InitializeIteration()368 void InitializeIteration() override 369 { 370 Superclass::InitializeIteration(); 371 // Estimate the progress of the filter 372 this->UpdateProgress( (float)( (float)this->GetElapsedIterations() 373 / (float)this->GetNumberOfIterations() ) ); 374 } 375 376 /** Tells the solver how to reinitialize the narrowband when the reinitialization 377 * criterion meets */ 378 379 void CreateNarrowBand() override; 380 381 /** Overridden from ProcessObject to set certain values before starting the 382 * finite difference solver and then create an appropriate output */ 383 void GenerateData() override; 384 385 /** Flag which sets the inward/outward direction of propagation speed. See 386 SetReverseExpansionDirection for more information. */ 387 bool m_ReverseExpansionDirection; 388 389 /** Reinitialization filters **/ 390 /** Internal filter types used for reinitialization */ 391 using IsoFilterType = 392 IsoContourDistanceImageFilter< OutputImageType, OutputImageType >; 393 using ChamferFilterType = 394 FastChamferDistanceImageFilter< OutputImageType, OutputImageType >; 395 396 typename IsoFilterType::Pointer m_IsoFilter; 397 398 typename ChamferFilterType::Pointer m_ChamferFilter; 399 400 private: 401 SegmentationFunctionType *m_SegmentationFunction; 402 }; 403 } // end namespace itk 404 405 #ifndef ITK_MANUAL_INSTANTIATION 406 #include "itkNarrowBandLevelSetImageFilter.hxx" 407 #endif 408 409 #endif 410