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 19 #ifndef itkFastMarchingBase_h 20 #define itkFastMarchingBase_h 21 22 #include "itkIntTypes.h" 23 #include "itkFastMarchingStoppingCriterionBase.h" 24 #include "itkFastMarchingTraits.h" 25 26 #include <queue> 27 #include <functional> 28 29 namespace itk 30 { 31 /** 32 * \class FastMarchingBase 33 * \brief Abstract class to solve an Eikonal based-equation using Fast Marching 34 * Method. 35 * 36 * Fast marching solves an Eikonal equation where the speed is always 37 * non-negative and depends on the position only. Starting from an 38 * initial position on the front, fast marching systematically moves the 39 * front forward one node at a time. 40 * 41 * Updates are preformed using an entropy satisfy scheme where only 42 * "upwind" neighborhoods are used. This implementation of Fast Marching 43 * uses a std::priority_queue to locate the next proper node to 44 * update. 45 * 46 * Fast Marching sweeps through N points in (N log N) steps to obtain 47 * the arrival time value as the front propagates through the domain. 48 * 49 * The initial front is specified by two containers: 50 * \li one containing the known nodes (Alive Nodes: nodes that are already 51 * part of the object), 52 * \li one containing the trial nodes (Trial Nodes: nodes that are 53 * considered for inclusion). 54 * 55 * In order for the filter to evolve, at least some trial nodes must be 56 * specified. These can for instance be specified as the layer of 57 * nodes around the alive ones. 58 * 59 * The algorithm is terminated early by setting an appropriate stopping 60 * criterion, or if there are no more nodes to process. 61 * 62 * \tparam TTraits traits which includes definition such as: 63 * \li InputDomainType (itk::Image or itk::QuadEdgeMesh) 64 * \li OutputDomainType (similar to InputDomainType) 65 * \li NodeType (itk::Index if itk::Image and PointIdentifier if 66 * itk::QuadEdgeMesh) 67 * \li NodePairType std::pair< NodeType, OutputPixelType > 68 * \li Superclass (itk::ImageToImageFilter or 69 * itk::QuadEdgeMeshToQuadEdgeMeshFilter ) 70 * 71 * \todo In the current implementation, std::priority_queue only allows 72 * taking nodes out from the front and putting nodes in from the back. 73 * Use itk::PriorityQueueContainer instead. 74 * 75 * \par Topology constraints: 76 * Additional flexibiility in this class includes the implementation of 77 * topology constraints for image-based fast marching. Further details 78 * can be found in the paper 79 * 80 * NJ Tustison, BA Avants, MF Siqueira, JC Gee. "Topological Well- 81 * Composedness and Glamorous Glue: A Digital Gluing Algorithm for 82 * Topologically Constrained Front Propagation, IEEE Transactions on 83 * Image Processing, 20(6):1756-1761, June 2011. 84 * 85 * Essentially, one can constrain the propagating front(s) such that 86 * they either: 87 * 1. don't merge (using the "Strict" option) 88 * 2. don't create handles (using the "NoHandles" option) 89 * 90 * Whereas the majority of related work uses the digital topological 91 * concept of "simple points" to constrain the evolving front, this 92 * filter uses the concept of "well-composedness". Advantages of 93 * the latter over the former includes being able to use the standard 94 * marching cubes algorithm to produce a mesh whose genus is identical 95 * to that of the evolved front(s). 96 * 97 * \sa FastMarchingStoppingCriterionBase 98 * 99 * \ingroup ITKFastMarching 100 */ 101 template< typename TInput, typename TOutput > 102 class ITK_TEMPLATE_EXPORT FastMarchingBase : public FastMarchingTraits<TInput, TOutput>::SuperclassType 103 { 104 public: 105 ITK_DISALLOW_COPY_AND_ASSIGN(FastMarchingBase); 106 107 using Traits = FastMarchingTraits<TInput, TOutput>; 108 using SuperclassType = typename Traits::SuperclassType; 109 110 using Self = FastMarchingBase; 111 using Superclass = typename FastMarchingTraits<TInput, TOutput>::SuperclassType; 112 using Pointer = SmartPointer< Self >; 113 using ConstPointer = SmartPointer< const Self >; 114 115 /** Run-time type information (and related methods). */ 116 itkTypeMacro(FastMarchingBase, FastMarchingTraits); 117 118 /** Input Domain related definitions */ 119 using InputDomainType = typename Traits::InputDomainType; 120 using InputDomainPointer = typename Traits::InputDomainPointer; 121 using InputPixelType = typename Traits::InputPixelType; 122 123 /** Output Domain related definitions */ 124 using OutputDomainType = typename Traits::OutputDomainType; 125 using OutputDomainPointer = typename Traits::OutputDomainPointer; 126 using OutputPixelType = typename Traits::OutputPixelType; 127 128 /** NodeType type of node */ 129 using NodeType = typename Traits::NodeType; 130 131 /** NodePairType pair of node and corresponding value */ 132 using NodePairType = typename Traits::NodePairType; 133 using NodePairContainerType = typename Traits::NodePairContainerType; 134 using NodePairContainerPointer = typename Traits::NodePairContainerPointer; 135 using NodePairContainerConstIterator = typename Traits::NodePairContainerConstIterator; 136 137 using LabelType = typename Traits::LabelType; 138 139 /** StoppingCriterionType stopping criterion */ 140 using StoppingCriterionType = FastMarchingStoppingCriterionBase< TInput, TOutput >; 141 using StoppingCriterionPointer = typename StoppingCriterionType::Pointer; 142 143 /* 144 using ElementIdentifier = long; 145 146 using PriorityQueueElementType = MinPriorityQueueElementWrapper< NodeType, 147 OutputPixelType, 148 ElementIdentifier >; 149 150 using PriorityQueueType = PriorityQueueContainer< PriorityQueueElementType, 151 PriorityQueueElementType, OutputPixelType, ElementIdentifier >; 152 using PriorityQueuePointer = typename PriorityQueueType::Pointer; 153 */ 154 155 /** \enum TopologyCheckType */ 156 enum TopologyCheckType { 157 /** \c Nothing */ 158 Nothing = 0, 159 /** \c NoHandles */ 160 NoHandles, 161 /** \c Strict */ 162 Strict }; 163 164 /** Set/Get the TopologyCheckType macro indicating whether the user 165 wants to check topology (and which one). */ 166 itkSetMacro( TopologyCheck, TopologyCheckType ); 167 itkGetConstReferenceMacro( TopologyCheck, TopologyCheckType ); 168 169 /** Set/Get TrialPoints */ 170 itkSetObjectMacro( TrialPoints, NodePairContainerType ); 171 itkGetModifiableObjectMacro(TrialPoints, NodePairContainerType ); 172 173 /** Set/Get AlivePoints */ 174 itkSetObjectMacro( AlivePoints, NodePairContainerType ); 175 itkGetModifiableObjectMacro(AlivePoints, NodePairContainerType ); 176 177 /** Set/Get ProcessedPoints */ 178 itkSetObjectMacro( ProcessedPoints, NodePairContainerType ); 179 itkGetModifiableObjectMacro(ProcessedPoints, NodePairContainerType ); 180 181 /** Set/Get ForbiddenPoints */ 182 itkSetObjectMacro( ForbiddenPoints, NodePairContainerType ); 183 itkGetModifiableObjectMacro(ForbiddenPoints, NodePairContainerType ); 184 185 /** \brief Set/Get the Stopping Criterion */ 186 itkSetObjectMacro( StoppingCriterion, StoppingCriterionType ); 187 itkGetModifiableObjectMacro(StoppingCriterion, StoppingCriterionType ); 188 189 /** \brief Set/Get SpeedConstant */ 190 itkGetMacro( SpeedConstant, double ); 191 itkSetMacro( SpeedConstant, double ); 192 193 /** \brief Set/Get NormalizationFactor */ 194 itkGetMacro( NormalizationFactor, double ); 195 itkSetMacro( NormalizationFactor, double ); 196 197 /** \brief Get the value reached by the front when it stops propagating */ 198 itkGetMacro( TargetReachedValue, OutputPixelType ); 199 200 /** Set the Collect Points flag. Instrument the algorithm to collect 201 * a container of all nodes which it has visited. Useful for 202 * creating Narrowbands for level set algorithms that supports 203 * narrow banding. */ 204 itkSetMacro(CollectPoints, bool); 205 206 /** Get the Collect Points flag. */ 207 itkGetConstReferenceMacro(CollectPoints, bool); 208 itkBooleanMacro(CollectPoints); 209 210 protected: 211 212 /** \brief Constructor */ 213 FastMarchingBase(); 214 215 /** \brief Destructor */ 216 ~FastMarchingBase() override = default; 217 218 StoppingCriterionPointer m_StoppingCriterion; 219 220 double m_SpeedConstant; 221 double m_InverseSpeed; 222 double m_NormalizationFactor; 223 224 OutputPixelType m_TargetReachedValue; 225 OutputPixelType m_LargeValue; 226 OutputPixelType m_TopologyValue; 227 228 NodePairContainerPointer m_TrialPoints; 229 NodePairContainerPointer m_AlivePoints; 230 NodePairContainerPointer m_ProcessedPoints; 231 NodePairContainerPointer m_ForbiddenPoints; 232 233 bool m_CollectPoints; 234 235 //PriorityQueuePointer m_Heap; 236 using HeapContainerType = std::vector< NodePairType >; 237 using NodeComparerType = std::greater< NodePairType >; 238 239 using PriorityQueueType = std::priority_queue< 240 NodePairType, HeapContainerType, NodeComparerType >; 241 242 PriorityQueueType m_Heap; 243 244 TopologyCheckType m_TopologyCheck; 245 246 /** \brief Get the total number of nodes in the domain */ 247 virtual IdentifierType GetTotalNumberOfNodes() const = 0; 248 249 /** \brief Get the output value (front value) for a given node */ 250 virtual const OutputPixelType GetOutputValue( OutputDomainType* oDomain, 251 const NodeType& iNode ) const = 0; 252 253 /** \brief Set the output value (front value) for a given node */ 254 virtual void SetOutputValue( OutputDomainType* oDomain, 255 const NodeType& iNode, 256 const OutputPixelType& iValue ) = 0; 257 258 /** \brief Get the LabelType Value for a given node 259 \param[in] iNode 260 \return its label value */ 261 virtual unsigned char 262 GetLabelValueForGivenNode( const NodeType& iNode ) const = 0; 263 264 /** \brief Set the Label Value for a given node 265 \param[in] iNode 266 \param[in] iLabel */ 267 virtual void SetLabelValueForGivenNode( const NodeType& iNode, 268 const LabelType& iLabel ) = 0; 269 270 /** \brief Update neighbors to a given node 271 \param[in] oDomain 272 \param[in] iNode 273 */ 274 virtual void UpdateNeighbors( OutputDomainType* oDomain, 275 const NodeType& iNode ) = 0; 276 277 /** \brief Update value for a given node 278 \param[in] oDomain 279 \param[in] iNode 280 */ 281 virtual void UpdateValue( OutputDomainType* oDomain, 282 const NodeType& iNode ) = 0; 283 284 /** \brief Check if the current node violate topological criterion. 285 \param[in] oDomain 286 \param[in] iNode 287 */ 288 virtual bool CheckTopology( OutputDomainType* oDomain, 289 const NodeType& iNode ) = 0; 290 291 /** \brief */ 292 void Initialize( OutputDomainType* oDomain ); 293 294 /** */ 295 virtual void InitializeOutput( OutputDomainType* oDomain ) = 0; 296 297 /** */ 298 void GenerateData() override; 299 300 /** \brief PrintSelf method */ 301 void PrintSelf(std::ostream & os, Indent indent) const override; 302 }; 303 } 304 305 #ifndef ITK_MANUAL_INSTANTIATION 306 #include "itkFastMarchingBase.hxx" 307 #endif 308 309 #endif 310