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 itkContourExtractor2DImageFilter_h 19 #define itkContourExtractor2DImageFilter_h 20 21 #include "itkImageToPathFilter.h" 22 #include "itkPolyLineParametricPath.h" 23 #include "itkConceptChecking.h" 24 #include <unordered_map> 25 #include <deque> 26 #include <list> 27 28 namespace itk 29 { 30 /** \class ContourExtractor2DImageFilter 31 * \brief Computes a list of PolyLineParametricPath objects from the contours in 32 * a 2D image. 33 * 34 * Uses the "marching squares" method to compute a the iso-valued contours of 35 * the input 2D image for a given intensity value. Multiple outputs may be 36 * produced because an image can have multiple contours at a given level, so it 37 * is advised to call GetNumberOfIndexedOutputs() and GetOutput(n) to retrieve all of 38 * the contours. The contour value to be extracted can be set with 39 * SetContourValue(). Image intensities will be linearly interpolated to provide 40 * sub-pixel resolution for the output contours. 41 * 42 * The marching squares algorithm is a special case of the marching cubes 43 * algorithm (Lorensen, William and Harvey E. Cline. Marching Cubes: A High 44 * Resolution 3D Surface Construction Algorithm. Computer Graphics (SIGGRAPH 87 45 * Proceedings) 21(4) July 1987, p. 163-170). A simple explanation is available 46 * here: http://users.polytech.unice.fr/~lingrand/MarchingCubes/algo.html 47 * 48 * There is a single ambiguous case in the marching squares algorithm: if a 49 * given 2x2-pixel square has two high-valued and two low-valued pixels, each 50 * pair diagonally adjacent. (Where high- and low-valued is with respect to the 51 * contour value sought.) In this case, either the high-valued pixels can be 52 * connected into the same "object" (where groups of pixels encircled by a given 53 * contour are considered an object), or the low-valued pixels can be connected. 54 * This is the "face connected" versus "face + vertex connected" (or 4- versus 55 * 4-connected) distinction: high-valued pixels most be treated as one, and 56 * low-valued as the other. By default, high-valued pixels are treated as 57 * "face-connected" and low-valued pixels are treated as "face + vertex" 58 * connected. To reverse this, call VertexConnectHighPixelsOn(); 59 * 60 * Outputs are not guaranteed to be closed paths: contours which intersect the 61 * image edge will be left open. All other paths will be closed. (The 62 * closedness of a path can be tested by checking whether the beginning point 63 * is the same as the end point.) 64 * 65 * Produced paths are oriented. Following the path from beginning to end, image 66 * intensity values lower than the contour value are to the left of the path and 67 * intensity values greater than the contour value are to the right. In other 68 * words, the image gradient at a path segment is (approximately) in the direct 69 * of that segment rotated right by 90 degrees, because the image intensity 70 * values increase from left-to-right across the segment. This means that the 71 * generated contours will circle clockwise around "hills" of 72 * above-contour-value intensity, and counter-clockwise around "depressions" of 73 * below-contour-value intensity. This convention can be reversed by calling 74 * ReverseContourOrientationOn(). 75 * 76 * By default the input image's largest possible region will be processed; call 77 * SetRequestedRegion() to process a different region, or ClearRequestedRegion() 78 * to revert to the default value. Note that the requested regions are usually 79 * set on the output; however since paths have no notion of a "region", this 80 * must be set at the filter level. 81 * 82 * This class was contributed to the Insight Journal by Zachary Pincus. 83 * https://hdl.handle.net/1926/165 84 * 85 * \sa Image 86 * \sa Path 87 * \sa PolyLineParametricPath 88 * 89 * \ingroup ITKPath 90 * 91 * \wiki 92 * \wikiexample{Segmentation/ContourExtractor2DImageFilter,Extract contours from an image} 93 * \endwiki 94 */ 95 96 template< typename TInputImage > 97 class ITK_TEMPLATE_EXPORT ContourExtractor2DImageFilter: 98 public ImageToPathFilter< TInputImage, PolyLineParametricPath< 2 > > 99 { 100 public: 101 ITK_DISALLOW_COPY_AND_ASSIGN(ContourExtractor2DImageFilter); 102 103 /** Extract dimension from input and output image. */ 104 static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; 105 106 /** Convenient type alias for simplifying declarations. */ 107 using InputImageType = TInputImage; 108 using OutputPathType = PolyLineParametricPath< 2 >; 109 110 /** Standard class type aliases. */ 111 using Self = ContourExtractor2DImageFilter; 112 using Superclass = ImageToPathFilter< InputImageType, OutputPathType >; 113 using Pointer = SmartPointer< Self >; 114 using ConstPointer = SmartPointer< const Self >; 115 116 /** Method for creation through the object factory. */ 117 itkNewMacro(Self); 118 119 /** Run-time type information (and related methods). */ 120 itkTypeMacro(ContourExtractor2DImageFilter, ImageToPathFilter); 121 122 /** Image and path type alias support */ 123 using InputImagePointer = typename InputImageType::Pointer; 124 using InputPixelType = typename InputImageType::PixelType; 125 using InputIndexType = typename InputImageType::IndexType; 126 using InputOffsetType = typename InputImageType::OffsetType; 127 using InputRegionType = typename InputImageType::RegionType; 128 using OutputPathPointer = typename OutputPathType::Pointer; 129 using VertexType = typename OutputPathType::VertexType; 130 using VertexListType = typename OutputPathType::VertexListType; 131 132 /** Real type associated to the input pixel type. */ 133 using InputRealType = typename NumericTraits< InputPixelType >::RealType; 134 135 using VertexListConstPointer = typename VertexListType::ConstPointer; 136 137 /** Control the orientation of the contours with reference to the image 138 * gradient. (See class documentation.) */ 139 itkSetMacro(ReverseContourOrientation, bool); 140 itkGetConstReferenceMacro(ReverseContourOrientation, bool); 141 itkBooleanMacro(ReverseContourOrientation); 142 143 /** Control whether high- or low-valued pixels are vertex-connected. 144 * Default is for low-valued pixels to be vertex-connected. 145 * (See class documentation.) */ 146 itkSetMacro(VertexConnectHighPixels, bool); 147 itkGetConstReferenceMacro(VertexConnectHighPixels, bool); 148 itkBooleanMacro(VertexConnectHighPixels); 149 150 /** Control whether the largest possible input region is used, or if a 151 * custom requested region is to be used. */ 152 void SetRequestedRegion(const InputRegionType region); 153 154 itkGetConstReferenceMacro(RequestedRegion, InputRegionType); 155 void ClearRequestedRegion(); 156 157 /** Set/Get the image intensity value that the contours should follow. 158 * This is the equivalent of an iso-value in Marching Squares. */ 159 itkSetMacro(ContourValue, InputRealType); 160 itkGetConstReferenceMacro(ContourValue, InputRealType); 161 162 #ifdef ITK_USE_CONCEPT_CHECKING 163 // Begin concept checking 164 itkConceptMacro( DimensionShouldBe2, 165 ( Concept::SameDimension< Self::InputImageDimension, 2 > ) ); 166 itkConceptMacro( InputPixelTypeComparable, 167 ( Concept::Comparable< InputPixelType > ) ); 168 itkConceptMacro( InputHasPixelTraitsCheck, 169 ( Concept::HasPixelTraits< InputPixelType > ) ); 170 itkConceptMacro( InputHasNumericTraitsCheck, 171 ( Concept::HasNumericTraits< InputPixelType > ) ); 172 // End concept checking 173 #endif 174 175 protected: 176 177 ContourExtractor2DImageFilter(); 178 ~ContourExtractor2DImageFilter() override = default; 179 void PrintSelf(std::ostream & os, Indent indent) const override; 180 181 void GenerateData() override; 182 183 /** ContourExtractor2DImageFilter manually controls the input requested 184 * region via SetRequestedRegion and ClearRequestedRegion, so it must 185 * override the superclass method. */ 186 void GenerateInputRequestedRegion() override; 187 188 private: 189 VertexType InterpolateContourPosition(InputPixelType fromValue, 190 InputPixelType toValue, 191 InputIndexType fromIndex, 192 InputOffsetType toOffset); 193 194 void AddSegment(const VertexType from, const VertexType to); 195 196 void FillOutputs(); 197 198 InputRealType m_ContourValue; 199 bool m_ReverseContourOrientation; 200 bool m_VertexConnectHighPixels; 201 bool m_UseCustomRegion; 202 InputRegionType m_RequestedRegion; 203 unsigned int m_NumberOfContoursCreated; 204 205 // Represent each contour as deque of vertices to facilitate addition of 206 // nodes at beginning or end. At the end of the processing, we will copy 207 // the contour into a PolyLineParametricPath. 208 // We subclass the deque to store an additional bit of information: an 209 // identification number for each growing contour. We use this number so 210 // that when it becomes necessary to merge two growing contours, we can 211 // merge the newer one into the older one. This helps because then we can 212 // guarantee that the output contour list is ordered from left to right, 213 // top to bottom (in terms of the first pixel of the contour encountered 214 // by the marching squares). Currently we make no guarantees that this 215 // pixel is the first pixel in the contour list, just that the contours 216 // are so ordered in the output. Ensuring this latter condition (first 217 // pixel traversed = first pixel in contour) would be possible by either 218 // changing the merging rules, which would make the contouring operation 219 //slower, or by storing additional data as to which pixel was first. 220 class ContourType:public std::deque< VertexType > 221 { 222 public: 223 unsigned int m_ContourNumber; 224 }; 225 226 // Store all the growing contours in a list. We may need to delete contours 227 // from anywhere in the sequence (when we merge them together), so we need to 228 // use a list instead of a vector or similar. 229 using ContourContainer = std::list< ContourType >; 230 using ContourRef = typename ContourContainer::iterator; 231 232 // declare the hash function we are using for the hash_map. 233 struct VertexHash { 234 using CoordinateType = typename VertexType::CoordRepType; operatorVertexHash235 inline SizeValueType operator()(const VertexType & k) const 236 { 237 // Xor the hashes of the vertices together, after multiplying the 238 // first by some number, so that identical (x,y) vertex indices 239 // don't all hash to the same bucket. This is a decent if not 240 // optimal hash. 241 const SizeValueType hashVertex1 = this->float_hash(k[0] * 0xbeef); 242 const SizeValueType hashVertex2 = this->float_hash(k[1]); 243 const SizeValueType hashValue = hashVertex1 ^ hashVertex2; 244 245 return hashValue; 246 } 247 248 // Define hash function for floats. Based on method from 249 // http://www.brpreiss.com/books/opus4/html/page217.html float_hashVertexHash250 inline SizeValueType float_hash(const CoordinateType & k) const 251 { 252 if ( k == 0 ) 253 { 254 return 0; 255 } 256 int exponent; 257 CoordinateType mantissa = std::frexp(k, &exponent); 258 auto value = static_cast< SizeValueType >( std::fabs(mantissa) ); 259 value = ( 2 * value - 1 ) * ~0U; 260 return value; 261 } 262 }; 263 264 // We use a hash to associate the endpoints of each contour with the 265 // contour itself. This makes it easy to look up which contour we should add 266 // a new arc to. 267 // We can't store the contours themselves in the hashtable because we 268 // need to have two tables (one to hash from beginpoint -> contour and one 269 // for endpoint -> contour), and sometimes will remove a contour from the 270 // tables (if it has been closed or merged with another contour). So in the 271 // hash table we store a reference to the contour. Because sometimes we will 272 // need to merge contours, we need to be able to quickly remove contours 273 // from our list when they have been merged into another. Thus, we store 274 // an iterator pointing to the contour in the list. 275 276 using VertexToContourMap = std::unordered_map< VertexType, ContourRef, VertexHash >; 277 using VertexMapIterator = typename VertexToContourMap::iterator; 278 using VertexContourRefPair = typename VertexToContourMap::value_type; 279 280 // The contours we find in the image are stored here 281 ContourContainer m_Contours; 282 283 // And indexed by their beginning and ending points here 284 VertexToContourMap m_ContourStarts; 285 VertexToContourMap m_ContourEnds; 286 }; 287 } // end namespace itk 288 289 #ifndef ITK_MANUAL_INSTANTIATION 290 #include "itkContourExtractor2DImageFilter.hxx" 291 #endif 292 293 #endif 294