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 itkScanlineFilterCommon_h 19 #define itkScanlineFilterCommon_h 20 21 #include "itkImageToImageFilter.h" 22 #include "itkConstShapedNeighborhoodIterator.h" 23 #include <atomic> 24 #include <deque> 25 #include <functional> 26 #include <mutex> 27 #include <vector> 28 29 namespace itk 30 { 31 /** 32 * \class ScanlineFilterCommon 33 * \brief Helper class for a group of filters which operate on scan-lines. 34 * 35 * This implementation was taken from the Insight Journal paper: 36 * https://hdl.handle.net/1926/584 or 37 * http://www.insight-journal.org/browse/publication/176 38 * 39 * \ingroup ITKImageLabel 40 */ 41 template< typename TInputImage, typename TOutputImage > 42 class ScanlineFilterCommon 43 { 44 public: 45 ITK_DISALLOW_COPY_AND_ASSIGN(ScanlineFilterCommon); 46 47 using Self = ScanlineFilterCommon; 48 using Pointer = SmartPointer< Self >; 49 using ConstPointer = SmartPointer< const Self >; Register()50 void Register() const 51 { 52 Object* obj = static_cast< Object* >( m_EnclosingFilter.GetPointer() ); 53 obj->Register(); 54 } UnRegister()55 void UnRegister() const noexcept 56 { 57 Object* obj = static_cast< Object* >( m_EnclosingFilter.GetPointer() ); 58 obj->UnRegister(); 59 } New()60 static Pointer New() 61 { 62 Pointer smartPtr = ObjectFactory< Self >::Create(); 63 if ( smartPtr == nullptr ) 64 { 65 smartPtr = new Self( nullptr ); 66 } 67 smartPtr->UnRegister(); 68 return smartPtr; 69 } 70 71 /** 72 * Extract some information from the image types. Dimensionality 73 * of the two images is assumed to be the same. 74 */ 75 using OutputPixelType = typename TOutputImage::PixelType; 76 using InputPixelType = typename TInputImage::PixelType; 77 static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension; 78 static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension; 79 static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension; 80 using InputImageType = TInputImage; 81 using InputImagePointer = typename InputImageType::Pointer; 82 using InputImageConstPointer = typename InputImageType::ConstPointer; 83 using IndexType = typename TInputImage::IndexType; 84 using SizeType = typename TInputImage::SizeType; 85 using OffsetType = typename TInputImage::OffsetType; 86 using OutputImageType = TOutputImage; 87 using OutputImagePointer = typename OutputImageType::Pointer; 88 using OutputRegionType = typename TOutputImage::RegionType; 89 using RegionType = OutputRegionType; 90 using OutputIndexType = typename TOutputImage::IndexType; 91 using OutputSizeType = typename TOutputImage::SizeType; 92 using OutputOffsetType = typename TOutputImage::OffsetType; 93 using OutputImagePixelType = typename TOutputImage::PixelType; 94 95 #ifdef ITK_USE_CONCEPT_CHECKING 96 // Concept checking -- input and output dimensions must be the same 97 itkConceptMacro( SameDimension, 98 ( Concept::SameDimension< InputImageDimension, 99 OutputImageDimension > ) ); 100 #endif 101 102 using EnclosingFilter = ImageToImageFilter< TInputImage, TOutputImage >; 103 ScanlineFilterCommon(EnclosingFilter * enclosingFilter)104 ScanlineFilterCommon(EnclosingFilter* enclosingFilter): 105 m_EnclosingFilter( enclosingFilter ), 106 m_FullyConnected( false ) 107 { 108 } 109 ~ScanlineFilterCommon() = default; 110 111 protected: 112 113 using InternalLabelType = SizeValueType; 114 using OutSizeType = typename TOutputImage::RegionType::SizeType; 115 116 struct RunLength 117 { 118 SizeValueType length; 119 typename InputImageType::IndexType where; 120 InternalLabelType label; 121 122 RunLength( SizeValueType iLength, const IndexType& iWhere, 123 InternalLabelType iLabel = 0 ) : lengthRunLength124 length( iLength ), where( iWhere ), label( iLabel ) {} 125 }; 126 127 using LineEncodingType = std::vector< RunLength >; 128 using LineEncodingIterator = typename LineEncodingType::iterator; 129 using LineEncodingConstIterator = typename LineEncodingType::const_iterator; 130 131 using OffsetVectorType = std::vector< OffsetValueType >; 132 using OffsetVectorConstIterator = typename OffsetVectorType::const_iterator; 133 134 using LineMapType = std::vector< LineEncodingType >; 135 136 using UnionFindType = std::vector< InternalLabelType >; 137 using ConsecutiveVectorType = std::vector< OutputPixelType >; 138 IndexToLinearIndex(const IndexType & index)139 SizeValueType IndexToLinearIndex( const IndexType& index ) const 140 { 141 SizeValueType linearIndex = 0; 142 SizeValueType stride = 1; 143 RegionType requestedRegion = m_EnclosingFilter->GetOutput()->GetRequestedRegion(); 144 // ignore x axis, which is always full size 145 for ( unsigned dim = 1; dim < ImageDimension; dim++ ) 146 { 147 itkAssertOrThrowMacro( requestedRegion.GetIndex( dim ) <= index[dim], 148 "Index must be within the requested region!" ); 149 linearIndex += ( index[dim] - requestedRegion.GetIndex( dim ) ) * stride; 150 stride *= requestedRegion.GetSize( dim ); 151 } 152 return linearIndex; 153 } 154 InitUnion(InternalLabelType numberOfLabels)155 void InitUnion(InternalLabelType numberOfLabels) 156 { 157 m_UnionFind = UnionFindType( numberOfLabels + 1 ); 158 159 typename LineMapType::iterator MapBegin, MapEnd, LineIt; 160 MapBegin = m_LineMap.begin(); 161 MapEnd = m_LineMap.end(); 162 LineIt = MapBegin; 163 InternalLabelType label = 1; 164 for ( LineIt = MapBegin; LineIt != MapEnd; ++LineIt ) 165 { 166 LineEncodingIterator cIt; 167 for ( cIt = LineIt->begin(); cIt != LineIt->end(); ++cIt ) 168 { 169 cIt->label = label; 170 m_UnionFind[label] = label; 171 label++; 172 } 173 } 174 } 175 LookupSet(const InternalLabelType label)176 InternalLabelType LookupSet(const InternalLabelType label) 177 { 178 InternalLabelType l = label; 179 while (l != m_UnionFind[l]) 180 { 181 l = m_UnionFind[l]; //transitively sets equivalence 182 } 183 return l; 184 } 185 LinkLabels(const InternalLabelType label1,const InternalLabelType label2)186 void LinkLabels(const InternalLabelType label1, const InternalLabelType label2) 187 { 188 std::lock_guard<std::mutex> mutexHolder(m_Mutex); 189 InternalLabelType E1 = this->LookupSet(label1); 190 InternalLabelType E2 = this->LookupSet(label2); 191 192 if ( E1 < E2 ) 193 { 194 m_UnionFind[E2] = E1; 195 } 196 else 197 { 198 m_UnionFind[E1] = E2; 199 } 200 } 201 CreateConsecutive(OutputPixelType backgroundValue)202 SizeValueType CreateConsecutive(OutputPixelType backgroundValue) 203 { 204 const size_t N = m_UnionFind.size(); 205 206 m_Consecutive = ConsecutiveVectorType( N ); 207 m_Consecutive[ 0 ] = backgroundValue; 208 209 OutputPixelType consecutiveLabel = 0; 210 SizeValueType count = 0; 211 212 for ( size_t i = 1; i < N; i++ ) 213 { 214 const auto label = static_cast< size_t >( m_UnionFind[i] ); 215 if ( label == i ) 216 { 217 if ( consecutiveLabel == backgroundValue ) 218 { 219 ++consecutiveLabel; 220 } 221 m_Consecutive[label] = consecutiveLabel; 222 ++consecutiveLabel; 223 ++count; 224 } 225 } 226 return count; 227 } 228 CheckNeighbors(const OutputIndexType & A,const OutputIndexType & B)229 bool CheckNeighbors(const OutputIndexType & A, const OutputIndexType & B) const 230 { 231 // This checks whether the line encodings are really neighbors. The first 232 // dimension gets ignored because the encodings are along that axis. 233 for ( unsigned i = 1; i < OutputImageDimension; i++ ) 234 { 235 if ( Math::abs(A[i] - B[i]) > 1 ) 236 { 237 return false; 238 } 239 } 240 return true; 241 } 242 243 using CompareLinesCallback = std::function<void( 244 const LineEncodingConstIterator& currentRun, 245 const LineEncodingConstIterator& neighborRun, 246 OffsetValueType oStart, 247 OffsetValueType oLast 248 )>; 249 CompareLines(const LineEncodingType & current,const LineEncodingType & Neighbour,bool sameLineOffset,bool labelCompare,OutputPixelType background,CompareLinesCallback callback)250 void CompareLines(const LineEncodingType & current, const LineEncodingType & Neighbour, 251 bool sameLineOffset, bool labelCompare, OutputPixelType background, 252 CompareLinesCallback callback) 253 { 254 bool sameLine = sameLineOffset; 255 if ( sameLineOffset ) 256 { 257 OutputOffsetType Off = current[0].where - Neighbour[0].where; 258 259 for ( unsigned int i = 1; i < ImageDimension; i++ ) 260 { 261 if ( Off[i] != 0 ) 262 { 263 sameLine = false; 264 break; 265 } 266 } 267 } 268 269 OffsetValueType offset = 0; 270 if ( m_FullyConnected || sameLine ) 271 { 272 offset = 1; 273 } 274 275 LineEncodingConstIterator nIt, mIt, cIt; 276 277 mIt = Neighbour.begin(); // out marker iterator 278 279 for ( cIt = current.begin(); cIt != current.end(); ++cIt ) 280 { 281 if ( !labelCompare || cIt->label != InternalLabelType(background) ) 282 { 283 OffsetValueType cStart = cIt->where[0]; // the start x position 284 OffsetValueType cLast = cStart + cIt->length - 1; 285 286 if (labelCompare) 287 { 288 mIt = Neighbour.begin(); 289 } 290 291 for ( nIt = mIt; nIt != Neighbour.end(); ++nIt ) 292 { 293 if ( !labelCompare || cIt->label != nIt->label ) 294 { 295 OffsetValueType nStart = nIt->where[0]; 296 OffsetValueType nLast = nStart + nIt->length - 1; 297 298 // there are a few ways that neighbouring lines might overlap 299 // neighbor S------------------E 300 // current S------------------------E 301 //------------- 302 // neighbor S------------------E 303 // current S----------------E 304 //------------- 305 // neighbor S------------------E 306 // current S------------------E 307 //------------- 308 // neighbor S------------------E 309 // current S-------E 310 //------------- 311 OffsetValueType ss1 = nStart - offset; 312 // OffsetValueType ss2 = nStart + offset; 313 OffsetValueType ee1 = nLast - offset; 314 OffsetValueType ee2 = nLast + offset; 315 316 bool eq = false; 317 OffsetValueType oStart = 0; 318 OffsetValueType oLast = 0; 319 320 // the logic here can probably be improved a lot 321 if ( ( ss1 >= cStart ) && ( ee2 <= cLast ) ) 322 { 323 // case 1 324 eq = true; 325 oStart = ss1; 326 oLast = ee2; 327 } 328 else if ( ( ss1 <= cStart ) && ( ee2 >= cLast ) ) 329 { 330 // case 4 331 eq = true; 332 oStart = cStart; 333 oLast = cLast; 334 } 335 else if ( ( ss1 <= cLast ) && ( ee2 >= cLast ) ) 336 { 337 // case 2 338 eq = true; 339 oStart = ss1; 340 oLast = cLast; 341 } 342 else if ( ( ss1 <= cStart ) && ( ee2 >= cStart ) ) 343 { 344 // case 3 345 eq = true; 346 oStart = cStart; 347 oLast = ee2; 348 } 349 350 if ( eq ) 351 { 352 callback(cIt, nIt, oStart, oLast); 353 if ( sameLineOffset && oStart == cStart && oLast == cLast ) 354 { 355 mIt = nIt; 356 break; 357 } 358 } 359 360 if ( !sameLineOffset && ee1 >= cLast ) 361 { 362 // No point looking for more overlaps with the current run 363 // because the neighbor run is either case 2 or 4 364 mIt = nIt; 365 break; 366 } 367 } 368 } 369 } 370 } 371 } 372 SetupLineOffsets(bool wholeNeighborhood)373 void SetupLineOffsets(bool wholeNeighborhood) 374 { 375 // Create a neighborhood so that we can generate a table of offsets 376 // to "previous" line indexes 377 // We are going to mis-use the neighborhood iterators to compute the 378 // offset for us. All this messing around produces an array of 379 // offsets that will be used to index the map 380 typename TOutputImage::Pointer output = m_EnclosingFilter->GetOutput(); 381 using PretendImageType = Image< OffsetValueType, TOutputImage::ImageDimension - 1 >; 382 using PretendSizeType = typename PretendImageType::RegionType::SizeType; 383 using PretendIndexType = typename PretendImageType::RegionType::IndexType; 384 using LineNeighborhoodType = ConstShapedNeighborhoodIterator< PretendImageType >; 385 386 typename PretendImageType::Pointer fakeImage; 387 fakeImage = PretendImageType::New(); 388 389 typename PretendImageType::RegionType LineRegion; 390 391 OutSizeType OutSize = output->GetRequestedRegion().GetSize(); 392 393 PretendSizeType PretendSize; 394 // The first dimension has been collapsed 395 for ( SizeValueType i = 0; i < PretendSize.GetSizeDimension(); i++ ) 396 { 397 PretendSize[i] = OutSize[i + 1]; 398 } 399 400 LineRegion.SetSize(PretendSize); 401 fakeImage->SetRegions(LineRegion); 402 PretendSizeType kernelRadius; 403 kernelRadius.Fill(1); 404 LineNeighborhoodType lnit(kernelRadius, fakeImage, LineRegion); 405 406 if (wholeNeighborhood) 407 { 408 setConnectivity(&lnit, m_FullyConnected); 409 } 410 else 411 { 412 setConnectivityPrevious(&lnit, m_FullyConnected); 413 } 414 415 typename LineNeighborhoodType::IndexListType ActiveIndexes; 416 ActiveIndexes = lnit.GetActiveIndexList(); 417 418 typename LineNeighborhoodType::IndexListType::const_iterator LI; 419 420 PretendIndexType idx = LineRegion.GetIndex(); 421 OffsetValueType offset = fakeImage->ComputeOffset(idx); 422 423 for ( LI = ActiveIndexes.begin(); LI != ActiveIndexes.end(); ++LI ) 424 { 425 m_LineOffsets.push_back(fakeImage->ComputeOffset( idx + lnit.GetOffset(*LI) ) - offset); 426 } 427 428 if (wholeNeighborhood) 429 { 430 m_LineOffsets.push_back(0); //center pixel 431 } 432 } 433 434 WeakPointer<EnclosingFilter> m_EnclosingFilter; 435 436 struct WorkUnitData 437 { 438 SizeValueType firstLine; 439 SizeValueType lastLine; 440 }; 441 CreateWorkUnitData(const RegionType & outputRegionForThread)442 WorkUnitData CreateWorkUnitData( const RegionType& outputRegionForThread ) 443 { 444 const SizeValueType xsizeForThread = outputRegionForThread.GetSize()[0]; 445 const SizeValueType numberOfLines = outputRegionForThread.GetNumberOfPixels() / xsizeForThread; 446 447 const SizeValueType firstLine = this->IndexToLinearIndex( outputRegionForThread.GetIndex() ); 448 const SizeValueType lastLine = firstLine + numberOfLines - 1; 449 450 return WorkUnitData{ firstLine, lastLine }; 451 } 452 453 /* Process the map and make appropriate entries in an equivalence table */ ComputeEquivalence(const SizeValueType workUnitResultsIndex,bool strictlyLess)454 void ComputeEquivalence(const SizeValueType workUnitResultsIndex, bool strictlyLess) 455 { 456 const OffsetValueType linecount = m_LineMap.size(); 457 WorkUnitData wud = m_WorkUnitResults[workUnitResultsIndex]; 458 SizeValueType lastLine = wud.lastLine; 459 if (!strictlyLess) 460 { 461 lastLine++; 462 //make sure we are not wrapping around 463 itkAssertInDebugAndIgnoreInReleaseMacro( lastLine >= wud.lastLine ); 464 } 465 for ( SizeValueType thisIdx = wud.firstLine; thisIdx < lastLine; ++thisIdx ) 466 { 467 if ( !m_LineMap[thisIdx].empty() ) 468 { 469 typename OffsetVectorType::const_iterator it = this->m_LineOffsets.begin(); 470 while ( it != this->m_LineOffsets.end() ) 471 { 472 OffsetValueType neighIdx = thisIdx + ( *it ); 473 // check if the neighbor is in the map 474 if ( neighIdx >= 0 && neighIdx < linecount && !m_LineMap[neighIdx].empty() ) 475 { 476 // Now check whether they are really neighbors 477 bool areNeighbors = this->CheckNeighbors(m_LineMap[thisIdx][0].where, m_LineMap[neighIdx][0].where); 478 if ( areNeighbors ) 479 { 480 this->CompareLines( 481 m_LineMap[thisIdx], 482 m_LineMap[neighIdx], 483 false, 484 false, 485 0, 486 [this]( 487 const LineEncodingConstIterator& currentRun, 488 const LineEncodingConstIterator& neighborRun, 489 OffsetValueType, 490 OffsetValueType) 491 { 492 this->LinkLabels(neighborRun->label, currentRun->label); 493 }); 494 } 495 } 496 ++it; 497 } 498 } 499 } 500 } 501 502 protected: 503 bool m_FullyConnected; 504 OffsetVectorType m_LineOffsets; 505 UnionFindType m_UnionFind; 506 ConsecutiveVectorType m_Consecutive; 507 std::mutex m_Mutex; 508 509 std::atomic< SizeValueType > m_NumberOfLabels; 510 std::deque< WorkUnitData > m_WorkUnitResults; 511 LineMapType m_LineMap; 512 }; 513 } // end namespace itk 514 515 #endif 516