ITK  6.0.0
Insight Toolkit
itkScanlineFilterCommon.h
Go to the documentation of this file.
1/*=========================================================================
2 *
3 * Copyright NumFOCUS
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 * https://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
23#include <atomic>
24#include <deque>
25#include <functional>
26#include <mutex>
27#include <vector>
28
29namespace itk
30{
40template <typename TInputImage, typename TOutputImage>
42{
43public:
44 ITK_DISALLOW_COPY_AND_MOVE(ScanlineFilterCommon);
45
49 void
50 Register() const
51 {
52 auto * obj = static_cast<Object *>(m_EnclosingFilter.GetPointer());
53 obj->Register();
54 }
55 void
56 UnRegister() const noexcept
57 {
58 auto * obj = static_cast<Object *>(m_EnclosingFilter.GetPointer());
59 obj->UnRegister();
60 }
61 static Pointer
63 {
65 if (smartPtr == nullptr)
66 {
67 smartPtr = new Self(nullptr);
68 }
69 smartPtr->UnRegister();
70 return smartPtr;
71 }
72
77 using OutputPixelType = typename TOutputImage::PixelType;
78 using InputPixelType = typename TInputImage::PixelType;
79 static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
80 static constexpr unsigned int OutputImageDimension = TOutputImage::ImageDimension;
81 static constexpr unsigned int InputImageDimension = TInputImage::ImageDimension;
82 using InputImageType = TInputImage;
87 using OffsetType = typename TInputImage::OffsetType;
88 using OutputImageType = TOutputImage;
94 using OutputOffsetType = typename TOutputImage::OffsetType;
95 using OutputImagePixelType = typename TOutputImage::PixelType;
96
97#ifdef ITK_USE_CONCEPT_CHECKING
98 // Concept checking -- input and output dimensions must be the same
100#endif
101
103
105 : m_EnclosingFilter(enclosingFilter)
106 , m_FullyConnected(false)
107 {}
109
110protected:
113
115 {
119
120 RunLength(SizeValueType iLength, const IndexType & iWhere, InternalLabelType iLabel = 0)
121 : length(iLength)
122 , where(iWhere)
123 , label(iLabel)
124 {}
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
140 IndexToLinearIndex(const IndexType & index) const
141 {
142 SizeValueType linearIndex = 0;
143 SizeValueType stride = 1;
144 RegionType requestedRegion = m_EnclosingFilter->GetOutput()->GetRequestedRegion();
145 // ignore x axis, which is always full size
146 for (unsigned int dim = 1; dim < ImageDimension; ++dim)
147 {
148 itkAssertOrThrowMacro(requestedRegion.GetIndex(dim) <= index[dim], "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
155 void
157 {
158 m_UnionFind = UnionFindType(numberOfLabels + 1);
159
160 typename LineMapType::iterator MapBegin, MapEnd, LineIt;
161 MapBegin = m_LineMap.begin();
162 MapEnd = m_LineMap.end();
163 LineIt = MapBegin;
164 InternalLabelType label = 1;
165 for (LineIt = MapBegin; LineIt != MapEnd; ++LineIt)
166 {
168 for (cIt = LineIt->begin(); cIt != LineIt->end(); ++cIt)
169 {
170 cIt->label = label;
171 m_UnionFind[label] = label;
172 ++label;
173 }
174 }
175 }
176
179 {
180 InternalLabelType l = label;
181 while (l != m_UnionFind[l])
182 {
183 l = m_UnionFind[l]; // transitively sets equivalence
184 }
185 return l;
186 }
187
188 void
190 {
191 const std::lock_guard<std::mutex> lockGuard(m_Mutex);
192 InternalLabelType E1 = this->LookupSet(label1);
193 InternalLabelType E2 = this->LookupSet(label2);
194
195 if (E1 < E2)
196 {
197 m_UnionFind[E2] = E1;
198 }
199 else
200 {
201 m_UnionFind[E1] = E2;
202 }
203 }
204
207 {
208 const size_t N = m_UnionFind.size();
209
211 m_Consecutive[0] = backgroundValue;
212
213 OutputPixelType consecutiveLabel = 0;
214 SizeValueType count = 0;
215
216 for (size_t i = 1; i < N; ++i)
217 {
218 const auto label = static_cast<size_t>(m_UnionFind[i]);
219 if (label == i)
220 {
221 if (consecutiveLabel == backgroundValue)
222 {
223 ++consecutiveLabel;
224 }
225 m_Consecutive[label] = consecutiveLabel;
226 ++consecutiveLabel;
227 ++count;
228 }
229 }
230 return count;
231 }
232
233 bool
235 {
236 // This checks whether the line encodings are really neighbors. The first
237 // dimension gets ignored because the encodings are along that axis.
238 SizeValueType diffSum = 0;
239 for (unsigned int i = 1; i < OutputImageDimension; ++i)
240 {
241 SizeValueType diff = itk::Math::abs(A[i] - B[i]);
242 if (diff > 1)
243 {
244 return false;
245 }
246 diffSum += diff;
247 }
248
249 if (!this->m_FullyConnected)
250 {
251 return (diffSum <= 1); // indices can differ only along one dimension
252 }
253 return true;
254 }
255
256 using CompareLinesCallback = std::function<void(const LineEncodingConstIterator & currentRun,
257 const LineEncodingConstIterator & neighborRun,
258 OffsetValueType oStart,
259 OffsetValueType oLast)>;
260
261 void
263 const LineEncodingType & Neighbour,
264 bool sameLineOffset,
265 bool labelCompare,
266 OutputPixelType background,
267 CompareLinesCallback callback)
268 {
269 bool sameLine = sameLineOffset;
270 if (sameLineOffset)
271 {
272 OutputOffsetType Off = current[0].where - Neighbour[0].where;
273
274 for (unsigned int i = 1; i < ImageDimension; ++i)
275 {
276 if (Off[i] != 0)
277 {
278 sameLine = false;
279 break;
280 }
281 }
282 }
283
284 OffsetValueType offset = 0;
285 if (m_FullyConnected || sameLine)
286 {
287 offset = 1;
288 }
289
290 LineEncodingConstIterator nIt, mIt, cIt;
291
292 mIt = Neighbour.begin(); // out marker iterator
293
294 for (cIt = current.begin(); cIt != current.end(); ++cIt)
295 {
296 if (!labelCompare || cIt->label != InternalLabelType(background))
297 {
298 OffsetValueType cStart = cIt->where[0]; // the start x position
299 OffsetValueType cLast = cStart + cIt->length - 1;
300
301 if (labelCompare)
302 {
303 mIt = Neighbour.begin();
304 }
305
306 for (nIt = mIt; nIt != Neighbour.end(); ++nIt)
307 {
308 if (!labelCompare || cIt->label != nIt->label)
309 {
310 OffsetValueType nStart = nIt->where[0];
311 OffsetValueType nLast = nStart + nIt->length - 1;
312
313 // there are a few ways that neighbouring lines might overlap
314 // neighbor S------------------E
315 // current S------------------------E
316 //-------------
317 // neighbor S------------------E
318 // current S----------------E
319 //-------------
320 // neighbor S------------------E
321 // current S------------------E
322 //-------------
323 // neighbor S------------------E
324 // current S-------E
325 //-------------
326 OffsetValueType ss1 = nStart - offset;
327 // OffsetValueType ss2 = nStart + offset;
328 OffsetValueType ee1 = nLast - offset;
329 OffsetValueType ee2 = nLast + offset;
330
331 bool eq = false;
332 OffsetValueType oStart = 0;
333 OffsetValueType oLast = 0;
334
335 // the logic here can probably be improved a lot
336 if ((ss1 >= cStart) && (ee2 <= cLast))
337 {
338 // case 1
339 eq = true;
340 oStart = ss1;
341 oLast = ee2;
342 }
343 else if ((ss1 <= cStart) && (ee2 >= cLast))
344 {
345 // case 4
346 eq = true;
347 oStart = cStart;
348 oLast = cLast;
349 }
350 else if ((ss1 <= cLast) && (ee2 >= cLast))
351 {
352 // case 2
353 eq = true;
354 oStart = ss1;
355 oLast = cLast;
356 }
357 else if ((ss1 <= cStart) && (ee2 >= cStart))
358 {
359 // case 3
360 eq = true;
361 oStart = cStart;
362 oLast = ee2;
363 }
364
365 if (eq)
366 {
367 callback(cIt, nIt, oStart, oLast);
368 if (sameLineOffset && oStart == cStart && oLast == cLast)
369 {
370 mIt = nIt;
371 break;
372 }
373 }
374
375 if (!sameLineOffset && ee1 >= cLast)
376 {
377 // No point looking for more overlaps with the current run
378 // because the neighbor run is either case 2 or 4
379 mIt = nIt;
380 break;
381 }
382 }
383 }
384 }
385 }
386 }
387
388 void
389 SetupLineOffsets(bool wholeNeighborhood)
390 {
391 // Create a neighborhood so that we can generate a table of offsets
392 // to "previous" line indexes
393 // We are going to misuse the neighborhood iterators to compute the
394 // offset for us. All this messing around produces an array of
395 // offsets that will be used to index the map
396 typename TOutputImage::Pointer output = m_EnclosingFilter->GetOutput();
397 using PretendImageType = Image<OffsetValueType, TOutputImage::ImageDimension - 1>;
398 using PretendSizeType = typename PretendImageType::RegionType::SizeType;
399 using PretendIndexType = typename PretendImageType::RegionType::IndexType;
400 using LineNeighborhoodType = ConstShapedNeighborhoodIterator<PretendImageType>;
401
402 auto fakeImage = PretendImageType::New();
403
404 typename PretendImageType::RegionType LineRegion;
405
406 OutSizeType OutSize = output->GetRequestedRegion().GetSize();
407
408 PretendSizeType PretendSize;
409 // The first dimension has been collapsed
410 for (SizeValueType i = 0; i < PretendSize.GetSizeDimension(); ++i)
411 {
412 PretendSize[i] = OutSize[i + 1];
413 }
414
415 LineRegion.SetSize(PretendSize);
416 fakeImage->SetRegions(LineRegion);
417 auto kernelRadius = PretendSizeType::Filled(1);
418 LineNeighborhoodType lnit(kernelRadius, fakeImage, LineRegion);
419
420 if (wholeNeighborhood)
421 {
423 }
424 else
425 {
427 }
428
429 typename LineNeighborhoodType::IndexListType ActiveIndexes = lnit.GetActiveIndexList();
430
431 typename LineNeighborhoodType::IndexListType::const_iterator LI;
432
433 PretendIndexType idx = LineRegion.GetIndex();
434 OffsetValueType offset = fakeImage->ComputeOffset(idx);
435
436 for (LI = ActiveIndexes.begin(); LI != ActiveIndexes.end(); ++LI)
437 {
438 m_LineOffsets.push_back(fakeImage->ComputeOffset(idx + lnit.GetOffset(*LI)) - offset);
439 }
440
441 if (wholeNeighborhood)
442 {
443 m_LineOffsets.push_back(0); // center pixel
444 }
445 }
446
448
450 {
453 };
454
456 CreateWorkUnitData(const RegionType & outputRegionForThread)
457 {
458 const SizeValueType xsizeForThread = outputRegionForThread.GetSize()[0];
459 const SizeValueType numberOfLines = outputRegionForThread.GetNumberOfPixels() / xsizeForThread;
460
461 const SizeValueType firstLine = this->IndexToLinearIndex(outputRegionForThread.GetIndex());
462 const SizeValueType lastLine = firstLine + numberOfLines - 1;
463
464 return WorkUnitData{ firstLine, lastLine };
465 }
466
467 /* Process the map and make appropriate entries in an equivalence table */
468 void
469 ComputeEquivalence(const SizeValueType workUnitResultsIndex, bool strictlyLess)
470 {
471 const OffsetValueType linecount = m_LineMap.size();
472 WorkUnitData wud = m_WorkUnitResults[workUnitResultsIndex];
473 SizeValueType lastLine = wud.lastLine;
474 if (!strictlyLess)
475 {
476 ++lastLine;
477 // make sure we are not wrapping around
478 itkAssertInDebugAndIgnoreInReleaseMacro(lastLine >= wud.lastLine);
479 }
480 for (SizeValueType thisIdx = wud.firstLine; thisIdx < lastLine; ++thisIdx)
481 {
482 if (!m_LineMap[thisIdx].empty())
483 {
484 auto it = this->m_LineOffsets.begin();
485 while (it != this->m_LineOffsets.end())
486 {
487 OffsetValueType neighIdx = thisIdx + (*it);
488 // check if the neighbor is in the map
489 if (neighIdx >= 0 && neighIdx < linecount && !m_LineMap[neighIdx].empty())
490 {
491 // Now check whether they are really neighbors
492 bool areNeighbors = this->CheckNeighbors(m_LineMap[thisIdx][0].where, m_LineMap[neighIdx][0].where);
493 if (areNeighbors)
494 {
495 this->CompareLines(m_LineMap[thisIdx],
496 m_LineMap[neighIdx],
497 false,
498 false,
499 0,
500 [this](const LineEncodingConstIterator & currentRun,
501 const LineEncodingConstIterator & neighborRun,
503 OffsetValueType) { this->LinkLabels(neighborRun->label, currentRun->label); });
504 }
505 }
506 ++it;
507 }
508 }
509 }
510 }
511
512protected:
517 std::mutex m_Mutex;
518
519 std::atomic<SizeValueType> m_NumberOfLabels;
520 std::deque<WorkUnitData> m_WorkUnitResults;
522};
523} // end namespace itk
524
525#endif
Const version of ShapedNeighborhoodIterator, defining iteration of a local N-dimensional neighborhood...
const SizeType & GetSize() const
Base class for filters that take an image as input and produce an image as output.
Templated n-dimensional image class.
Definition: itkImage.h:89
static T::Pointer Create()
Base class for most ITK classes.
Definition: itkObject.h:62
void UnRegister() const noexcept override
void Register() const override
Helper class for a group of filters which operate on scan-lines.
typename TInputImage::IndexType IndexType
typename OffsetVectorType::const_iterator OffsetVectorConstIterator
typename TOutputImage::SizeType OutputSizeType
typename TOutputImage::RegionType::SizeType OutSizeType
typename LineEncodingType::const_iterator LineEncodingConstIterator
InternalLabelType LookupSet(const InternalLabelType label)
void InitUnion(InternalLabelType numberOfLabels)
ConsecutiveVectorType m_Consecutive
std::deque< WorkUnitData > m_WorkUnitResults
typename InputImageType::Pointer InputImagePointer
SizeValueType IndexToLinearIndex(const IndexType &index) const
typename TInputImage::OffsetType OffsetType
void ComputeEquivalence(const SizeValueType workUnitResultsIndex, bool strictlyLess)
void SetupLineOffsets(bool wholeNeighborhood)
bool CheckNeighbors(const OutputIndexType &A, const OutputIndexType &B) const
static constexpr unsigned int InputImageDimension
typename TOutputImage::PixelType OutputImagePixelType
typename LineEncodingType::iterator LineEncodingIterator
WorkUnitData CreateWorkUnitData(const RegionType &outputRegionForThread)
static constexpr unsigned int ImageDimension
void LinkLabels(const InternalLabelType label1, const InternalLabelType label2)
std::vector< LineEncodingType > LineMapType
SizeValueType CreateConsecutive(OutputPixelType backgroundValue)
typename TOutputImage::OffsetType OutputOffsetType
typename InputImageType::ConstPointer InputImageConstPointer
std::vector< RunLength > LineEncodingType
typename TInputImage::SizeType SizeType
std::vector< OffsetValueType > OffsetVectorType
std::vector< InternalLabelType > UnionFindType
std::vector< OutputPixelType > ConsecutiveVectorType
typename TOutputImage::RegionType OutputRegionType
void UnRegister() const noexcept
std::function< void(const LineEncodingConstIterator &currentRun, const LineEncodingConstIterator &neighborRun, OffsetValueType oStart, OffsetValueType oLast)> CompareLinesCallback
typename TOutputImage::IndexType OutputIndexType
std::atomic< SizeValueType > m_NumberOfLabels
typename TInputImage::PixelType InputPixelType
void CompareLines(const LineEncodingType &current, const LineEncodingType &Neighbour, bool sameLineOffset, bool labelCompare, OutputPixelType background, CompareLinesCallback callback)
WeakPointer< EnclosingFilter > m_EnclosingFilter
static constexpr unsigned int OutputImageDimension
typename OutputImageType::Pointer OutputImagePointer
ScanlineFilterCommon(EnclosingFilter *enclosingFilter)
typename TOutputImage::PixelType OutputPixelType
Implements a weak reference to an object.
SmartPointer< const Self > ConstPointer
static Pointer New()
SmartPointer< Self > Pointer
#define itkConceptMacro(name, concept)
bool abs(bool x)
Definition: itkMath.h:840
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
TIterator * setConnectivity(TIterator *it, bool fullyConnected=false)
unsigned long SizeValueType
Definition: itkIntTypes.h:86
TIterator * setConnectivityPrevious(TIterator *it, bool fullyConnected=false)
long OffsetValueType
Definition: itkIntTypes.h:97
RunLength(SizeValueType iLength, const IndexType &iWhere, InternalLabelType iLabel=0)
void SetSize(const SizeValueType val[VDimension])
Definition: itkSize.h:179