SMIL  0.9.1
DMorphoWatershed.hpp
1 /*
2  * Copyright (c) 2011-2016, Matthieu FAESSEL and ARMINES
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright
9  * notice, this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  * * Neither the name of Matthieu FAESSEL, or ARMINES nor the
14  * names of its contributors may be used to endorse or promote products
15  * derived from this software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS'' AND ANY
18  * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20  * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE LIABLE FOR ANY
21  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27  */
28 
29 
30 #ifndef _D_MORPHO_WATERSHED_HPP
31 #define _D_MORPHO_WATERSHED_HPP
32 
33 
34 #include "DMorphoHierarQ.hpp"
35 #include "DMorphoExtrema.hpp"
36 #include "DMorphoLabel.hpp"
37 #include "DMorphoResidues.hpp"
38 #include "Core/include/DTypes.h"
39 
40 
41 namespace smil
42 {
49  template <class T, class labelT, class HQ_Type=HierarchicalQueue<T> >
51  {
52  public:
53  BaseFlooding()
54  : STAT_QUEUED(ImDtTypes<labelT>::max())
55  {
56  }
57  virtual ~BaseFlooding()
58  {
59  }
60 
61  protected:
62 
63  const labelT STAT_QUEUED;
64 
65  typename ImDtTypes<T>::lineType inPixels;
66  typename ImDtTypes<labelT>::lineType lblPixels;
67 
68  size_t imSize[3], pixPerSlice;
69  inline void getCoordsFromOffset(size_t off, size_t &x, size_t &y, size_t &z) const
70  {
71  z = off / pixPerSlice;
72  y = (off % pixPerSlice) / imSize[0];
73  x = off % imSize[0];
74  }
75 
76  vector<IntPoint> sePts;
77  UINT sePtsNbr;
78  bool oddSE;
79  vector<int> dOffsets;
80 
81  T currentLevel;
82 
83  public:
84 
85  const Image<T> *imgIn;
86  Image<labelT> *imgLbl;
87 
88  virtual RES_T flood(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<labelT> &imBasinsOut, const StrElt &se=DEFAULT_SE)
89  {
90  ASSERT_ALLOCATED(&imIn, &imMarkers, &imBasinsOut);
91  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imBasinsOut);
92 
93  ASSERT(maxVal(imMarkers)<STAT_QUEUED);
94 
95  ImageFreezer freeze(imBasinsOut);
96 
97  copy(imMarkers, imBasinsOut);
98 
99  initialize(imIn, imBasinsOut, se);
100  processImage(imIn, imBasinsOut, se);
101 
102  return RES_OK;
103  }
104 
105  virtual RES_T initialize(const Image<T> &imIn, Image<labelT> &imLbl, const StrElt &se)
106  {
107  imgIn = &imIn;
108  imgLbl = &imLbl;
109 
110  // Empty the priority queue
111  hq.initialize(imIn);
112 
113  inPixels = imIn.getPixels();
114  lblPixels = imLbl.getPixels();
115 
116  imIn.getSize(imSize);
117  pixPerSlice = imSize[0]*imSize[1];
118 
119  dOffsets.clear();
120  sePts.clear();
121  oddSE = se.odd;
122 
123  // set an offset distance for each se point (!=0,0,0)
124  for(vector<IntPoint>::const_iterator it = se.points.begin() ; it!=se.points.end() ; it++)
125  if (it->x!=0 || it->y!=0 || it->z!=0)
126  {
127  sePts.push_back(*it);
128  dOffsets.push_back(it->x + it->y*imSize[0] + it->z*imSize[0]*imSize[1]);
129  }
130 
131  sePtsNbr = sePts.size();
132 
133  return RES_OK;
134  }
135 
136  virtual RES_T processImage(const Image<T> &/*imIn*/, Image<labelT> &/*imLbl*/, const StrElt &/*se*/)
137  {
138  // Put the marker pixels in the HQ
139  size_t offset = 0;
140  for (size_t k=0;k<imSize[2];k++)
141  for (size_t j=0;j<imSize[1];j++)
142  for (size_t i=0;i<imSize[0];i++)
143  {
144  if (lblPixels[offset]!=0)
145  {
146  hq.push(inPixels[offset], offset);
147  }
148  offset++;
149  }
150 
151 
152 
153  currentLevel = ImDtTypes<T>::min();
154 
155  while(!hq.isEmpty())
156  {
157  this->processPixel(hq.pop());
158  }
159 
160  return RES_OK;
161  }
162 
163  inline virtual void processPixel(const size_t &curOffset)
164  {
165  size_t x0, y0, z0;
166 
167  getCoordsFromOffset(curOffset, x0, y0, z0);
168 
169  bool oddLine = oddSE && ((y0)%2);
170 
171  int x, y, z;
172  size_t nbOffset;
173 
174  for(UINT i=0;i<sePtsNbr;i++)
175  {
176  IntPoint &pt = sePts[i];
177  x = x0 + pt.x;
178  y = y0 + pt.y;
179  z = z0 + pt.z;
180 
181  if (oddLine)
182  x += (((y+1)%2)!=0);
183 
184  if (x>=0 && x<(int)imSize[0] && y>=0 && y<(int)imSize[1] && z>=0 && z<(int)imSize[2])
185  {
186  nbOffset = curOffset + dOffsets[i];
187 
188  if (oddLine)
189  nbOffset += (((y+1)%2)!=0);
190 
191  processNeighbor(curOffset, nbOffset);
192 
193  }
194  }
195  }
196 
197  inline virtual void processNeighbor(const size_t &curOffset, const size_t &nbOffset)
198  {
199  labelT nbLbl = this->lblPixels[nbOffset];
200  labelT curLbl = lblPixels[curOffset];//==STAT_QUEUED ? 0 : lblPixels[curOffset];
201 
202  if (nbLbl==0) // Add it to the tmp offsets queue
203  {
204  hq.push(inPixels[nbOffset], nbOffset);
205  this->lblPixels[nbOffset] = curLbl;
206  }
207  }
208 
209  protected:
210  HQ_Type hq;
211 
212  };
213 
214 
215  template <class T, class labelT, class HQ_Type=HierarchicalQueue<T> >
217 #ifndef SWIG
218  : public BaseFlooding<T, labelT, HQ_Type>
219 #endif // SWIG
220  {
221  protected:
222  vector<size_t> tmpOffsets;
223  typename ImDtTypes<T>::lineType wsPixels;
224  const T STAT_LABELED, STAT_QUEUED, STAT_CANDIDATE, STAT_WS_LINE;
225 
226  public:
228  : STAT_LABELED(1),
229  STAT_QUEUED(2),
230  STAT_CANDIDATE(ImDtTypes<T>::max()-1),
231  STAT_WS_LINE(ImDtTypes<T>::max())
232  {
233  }
234  virtual ~WatershedFlooding()
235  {
236  }
237 
238  Image<T> *imgWS;
239 #ifdef SWIG
240  const Image<T> *imgIn;
241  Image<labelT> *imgLbl;
242 #endif // SWIG
243 
244  virtual RES_T flood(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<T> &imOut, Image<labelT> &imBasinsOut, const StrElt &se)
245  {
246  ASSERT_ALLOCATED(&imIn, &imMarkers, &imBasinsOut);
247  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imBasinsOut);
248 
249  ImageFreezer freeze(imBasinsOut);
250  ImageFreezer freeze2(imOut);
251 
252  copy(imMarkers, imBasinsOut);
253 
254  initialize(imIn, imBasinsOut, imOut, se);
255  this->processImage(imIn, imBasinsOut, se);
256 
257  // Finalize the image containing the ws lines
258  // Potential remaining candidate points (points surrounded by WS_LINE points) with status STAT_CANDIDATE are added to the WS_LINE
259  for (size_t i=0;i<imIn.getPixelCount();i++)
260  {
261  if (wsPixels[i]>=STAT_CANDIDATE)
262  wsPixels[i] = STAT_WS_LINE;
263  else
264  wsPixels[i] = 0;
265  }
266  return RES_OK;
267  }
268 
269  virtual RES_T initialize(const Image<T> &imIn, Image<labelT> &imLbl, Image<T> &imOut, const StrElt &se)
270  {
272 
273  imgWS = &imOut;
274  imgWS->setSize(this->imSize);
275  test(imLbl, STAT_LABELED, STAT_CANDIDATE, *imgWS);
276  wsPixels = imgWS->getPixels();
277 
278  tmpOffsets.clear();
279 
280  return RES_OK;
281  }
282 
283  inline virtual void processPixel(const size_t &curOffset)
284  {
285  wsPixels[curOffset] = STAT_LABELED;
286 
288 
289  if (!tmpOffsets.empty())
290  {
291  if (this->wsPixels[curOffset]!=STAT_WS_LINE)
292  {
293  size_t *offsets = tmpOffsets.data();
294  for (UINT i=0;i<tmpOffsets.size();i++)
295  {
296  this->hq.push(this->inPixels[*offsets], *offsets);
297  this->wsPixels[*offsets] = STAT_QUEUED;
298 
299  offsets++;
300  }
301  }
302  tmpOffsets.clear();
303  }
304 
305  }
306  inline virtual void processNeighbor(const size_t &curOffset, const size_t &nbOffset)
307  {
308  T nbStat = this->wsPixels[nbOffset];
309 
310  if (nbStat==STAT_CANDIDATE) // Add it to the tmp offsets queue
311  {
312  tmpOffsets.push_back(nbOffset);
313  }
314  else if (nbStat==STAT_LABELED)
315  {
316  if (this->lblPixels[curOffset]==0)
317  {
318  this->lblPixels[curOffset] = this->lblPixels[nbOffset];
319  }
320  else if (this->lblPixels[curOffset]!=this->lblPixels[nbOffset])
321  this->wsPixels[curOffset] = STAT_WS_LINE;
322  }
323  }
324 #ifndef SWIG
327 #endif // SWIG
328  };
329 
342  template <class T, class labelT>
343  RES_T basins(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<labelT> &imBasinsOut, const StrElt &se=DEFAULT_SE)
344  {
345  BaseFlooding<T, labelT> flooding;
346  return flooding.flood(imIn, imMarkers, imBasinsOut, se);
347  }
348 
349  template <class T, class labelT>
350  RES_T basins(const Image<T> &imIn, Image<labelT> &imBasinsInOut, const StrElt &se=DEFAULT_SE)
351  {
352  ASSERT_ALLOCATED(&imIn);
353  ASSERT_SAME_SIZE(&imIn, &imBasinsInOut);
354 
355  Image<labelT> imLbl(imIn);
356  minimaLabeled(imIn, imLbl, se);
357 
358  return basins(imIn, imLbl, imBasinsInOut);
359  }
360 
361 
376  template <class T, class labelT>
377  RES_T watershed(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<T> &imOut, Image<labelT> &imBasinsOut, const StrElt &se=DEFAULT_SE)
378  {
379  ASSERT_ALLOCATED(&imIn, &imMarkers);
380  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imOut, &imBasinsOut);
381 
383  return flooding.flood(imIn, imMarkers, imOut, imBasinsOut, se);
384  }
385 
386  template <class T, class labelT>
387  RES_T watershed(const Image<T> &imIn, Image<labelT> &imMarkers, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
388  {
389  ASSERT_ALLOCATED(&imIn, &imMarkers);
390  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imOut);
391 
392  Image<labelT> imBasinsOut(imMarkers);
393  return watershed(imIn, imMarkers, imOut, imBasinsOut, se);
394  }
395 
396  template <class T>
397  RES_T watershed(const Image<T> &imIn, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
398  {
399  ASSERT_ALLOCATED(&imIn);
400  ASSERT_SAME_SIZE(&imIn, &imOut);
401 
402  Image<UINT> imLbl(imIn);
403  minimaLabeled(imIn, imLbl, se);
404  return watershed(imIn, imLbl, imOut, se);
405  }
406 
414  template <class T>
415  RES_T lblSkiz(Image<T> &labelIm1, Image<T> &labelIm2, const Image<T> &maskIm, const StrElt &se=DEFAULT_SE)
416  {
417  ASSERT_ALLOCATED(&labelIm1, &labelIm2, &maskIm);
418  ASSERT_SAME_SIZE(&labelIm1, &labelIm2, &maskIm);
419 
420  Image<T> tmpIm1(labelIm1);
421  Image<T> tmpIm2(labelIm1);
422  Image<T> sumIm(labelIm1);
423 
424  double vol1 = -1, vol2 = 0;
425 
426  ImageFreezer freeze1(labelIm1);
427  ImageFreezer freeze2(labelIm2);
428 
429  T threshMin = ImDtTypes<T>::min(), threshMax = ImDtTypes<T>::max()-T(1);
430 
431  while(vol1<vol2)
432  {
433  dilate(labelIm1, tmpIm1, se);
434  dilate(labelIm2, tmpIm2, se);
435  addNoSat(labelIm1, labelIm2, sumIm);
436  threshold(sumIm, threshMin, threshMax, sumIm);
437  // if (&maskIm) // a reference cannot be NULL
438  inf(maskIm, sumIm, sumIm);
439  mask(tmpIm1, sumIm, tmpIm1);
440  mask(tmpIm2, sumIm, tmpIm2);
441  sup(labelIm1, tmpIm1, labelIm1);
442  sup(labelIm2, tmpIm2, labelIm2);
443 
444  vol1 = vol2;
445  vol2 = vol(labelIm1);
446  }
447  return RES_OK;
448  }
449 
457  template <class T1, class T2>
458  RES_T inflBasins(const Image<T1> &imIn, Image<T2> &basinsOut, const StrElt &se=DEFAULT_SE)
459  {
460  ASSERT_ALLOCATED(&imIn, &basinsOut);
461  ASSERT_SAME_SIZE(&imIn, &basinsOut);
462 
463  Image<T2> imLbl2(basinsOut);
464  Image<T2> nullIm(basinsOut);
465  fill(nullIm, ImDtTypes<T2>::max()); // invariant to liblSkiz inf
466 
467  // Create the label images
468  label(imIn, basinsOut, se);
469  inv(basinsOut, imLbl2);
470  mask(imLbl2, basinsOut, imLbl2);
471 
472  ASSERT(lblSkiz(basinsOut, imLbl2, nullIm, se)==RES_OK);
473 
474  // Clean result image
475  open(basinsOut, basinsOut);
476 
477  return RES_OK;
478  }
479 
487  template <class T>
488  RES_T inflZones(const Image<T> &imIn, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
489  {
490  ASSERT_ALLOCATED(&imIn, &imOut);
491  ASSERT_SAME_SIZE(&imIn, &imOut);
492 
493  Image<UINT16> basinsIm(imIn);
494 
495  ASSERT(inflBasins(imIn, basinsIm, se)==RES_OK);
496  gradient(basinsIm, basinsIm, se, StrElt());
497  threshold(basinsIm, UINT16(ImDtTypes<T>::min()+T(1)), UINT16(ImDtTypes<T>::max()), basinsIm);
498  copy(basinsIm, imOut);
499 
500  return RES_OK;
501  }
502 
507  template <class T>
508  RES_T waterfall(const Image<T> &gradIn, const Image<T> &wsIn, Image<T> &imGradOut, Image<T> &imWsOut, const StrElt &se=DEFAULT_SE)
509  {
510  ASSERT_ALLOCATED(&gradIn, &wsIn, &imGradOut, &imWsOut);
511  ASSERT_SAME_SIZE(&gradIn, &wsIn, &imGradOut, &imWsOut);
512 
513  ImageFreezer freeze(imWsOut);
514 
515  test(wsIn, gradIn, ImDtTypes<T>::max(), imWsOut);
516  dualBuild(imWsOut, gradIn, imGradOut, se);
517  watershed(imGradOut, imWsOut, se);
518 
519  return RES_OK;
520  }
521 
526  template <class T>
527  RES_T waterfall(const Image<T> &gradIn, UINT nLevel, Image<T> &imWsOut, const StrElt &se=DEFAULT_SE)
528  {
529  ASSERT_ALLOCATED(&gradIn, &imWsOut);
530  ASSERT_SAME_SIZE(&gradIn, &imWsOut);
531 
532  ImageFreezer freeze(imWsOut);
533 
534  Image<T> tmpGradIm(gradIn, 1); //clone
535  Image<T> tmpGradIm2(gradIn); //clone
536 
537  watershed(gradIn, imWsOut, se);
538 
539  for (UINT i=0;i<nLevel;i++)
540  {
541  waterfall(tmpGradIm, imWsOut, tmpGradIm2, imWsOut, se);
542  copy(tmpGradIm2, tmpGradIm);
543  }
544 
545  return RES_OK;
546  }
547 
550 } // namespace smil
551 
552 
553 #endif // _D_MORPHO_WATERSHED_HPP
554 
555 
RES_T minimaLabeled(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Calculate the minima and labelize them.
Definition: DMorphoExtrema.hpp:99
RES_T waterfall(const Image< T > &gradIn, const Image< T > &wsIn, Image< T > &imGradOut, Image< T > &imWsOut, const StrElt &se=DEFAULT_SE)
Waterfall.
Definition: DMorphoWatershed.hpp:508
RES_T threshold(const Image< T > &imIn, T minVal, T maxVal, T_out trueVal, T_out falseVal, Image< T_out > &imOut)
Image threshold.
Definition: DImageHistogram.hpp:170
RES_T dilate(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE, T borderVal=ImDtTypes< T >::min())
Morphological grayscale dilation.
Definition: DMorphoBase.hpp:65
Author Vincent Morard Date : 7 march 2011 See : Morard V. and Dokladal P. and Decenciere E...
Definition: DApplyThreshold.hpp:36
RES_T copy(const Image< T1 > &imIn, size_t startX, size_t startY, size_t startZ, size_t sizeX, size_t sizeY, size_t sizeZ, Image< T2 > &imOut, size_t outStartX=0, size_t outStartY=0, size_t outStartZ=0)
Copy image.
Definition: DImageArith.hpp:114
Definition: DBaseImage.h:235
double vol(const Image< T > &imIn)
Volume of an image.
Definition: Base/include/private/DMeasures.hpp:96
vector< IntPoint > points
List of neighbor points.
Definition: DStructuringElement.h:125
RES_T inv(const Image< T > &imIn, Image< T > &imOut)
Invert an image.
Definition: DImageArith.hpp:375
virtual RES_T setSize(size_t w, size_t h, size_t d=1, bool doAllocate=true)
Set the size of image.
Definition: DImage.hxx:319
RES_T basins(const Image< T > &imIn, const Image< labelT > &imMarkers, Image< labelT > &imBasinsOut, const StrElt &se=DEFAULT_SE)
Constrained basins.
Definition: DMorphoWatershed.hpp:343
RES_T inflBasins(const Image< T1 > &imIn, Image< T2 > &basinsOut, const StrElt &se=DEFAULT_SE)
Influences basins.
Definition: DMorphoWatershed.hpp:458
Definition: DMorphoWatershed.hpp:216
RES_T fill(Image< T > &imOut, const T &value)
Fill an image with a given value.
Definition: DImageArith.hpp:62
RES_T sup(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Sup of two images.
Definition: DImageArith.hpp:505
RES_T inflZones(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Influences zones.
Definition: DMorphoWatershed.hpp:488
RES_T lblSkiz(Image< T > &labelIm1, Image< T > &labelIm2, const Image< T > &maskIm, const StrElt &se=DEFAULT_SE)
Skiz on label image.
Definition: DMorphoWatershed.hpp:415
Base structuring element.
Definition: DStructuringElement.h:51
size_t label(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Image labelization.
Definition: DMorphoLabel.hpp:503
RES_T addNoSat(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Addition (without saturation check)
Definition: DImageArith.hpp:419
size_t getPixelCount() const
Get the number of pixels.
Definition: DBaseImage.h:150
Definition: DTypes.hpp:80
RES_T mask(const Image< T > &imIn, const Image< T > &imMask, Image< T > &imOut)
Image mask.
Definition: DImageArith.hpp:1060
RES_T dualBuild(const Image< T > &imIn, const Image< T > &imMask, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Dual reconstruction (using hierarchical queues).
Definition: DMorphoGeodesic.hpp:296
Main Image class.
Definition: DQVtkViewer.hpp:42
RES_T open(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Morphological grayscale opening.
Definition: DMorphoBase.hpp:159
RES_T inf(const Image< T > &imIn1, const T &value, Image< T > &imOut)
Inf of two images.
Definition: DImageArith.hpp:538
Definition: DMorphoWatershed.hpp:50
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
Max value of an image.
Definition: Base/include/private/DMeasures.hpp:260
RES_T gradient(const Image< T > &imIn, Image< T > &imOut, const StrElt &dilSe, const StrElt &eroSe)
Morphological gradient.
Definition: DMorphoResidues.hpp:49
RES_T test(const Image< T1 > &imIn1, const Image< T2 > &imIn2, const Image< T2 > &imIn3, Image< T2 > &imOut)
Test.
Definition: DImageArith.hpp:920
RES_T watershed(const Image< T > &imIn, const Image< labelT > &imMarkers, Image< T > &imOut, Image< labelT > &imBasinsOut, const StrElt &se=DEFAULT_SE)
Constrained watershed.
Definition: DMorphoWatershed.hpp:377
Definition: DCommon.h:113