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  // Propagate label on plateaus
206  if (inPixels[nbOffset]==inPixels[curOffset] && curLbl!=0)
207  this->lblPixels[nbOffset] = curLbl;
208  else
209  this->lblPixels[nbOffset] = STAT_QUEUED;
210  }
211  else if (nbLbl<STAT_QUEUED)
212  {
213  if (curLbl==0)
214  this->lblPixels[curOffset] = this->lblPixels[nbOffset];
215  }
216 
217 
218  }
219 
220 
221 
222 
223  protected:
224  HQ_Type hq;
225 
226  };
227 
228 
229  template <class T, class labelT, class HQ_Type=HierarchicalQueue<T> >
231 #ifndef SWIG
232  : public BaseFlooding<T, labelT, HQ_Type>
233 #endif // SWIG
234  {
235  protected:
236  vector<size_t> tmpOffsets;
237  typename ImDtTypes<T>::lineType wsPixels;
238  const T STAT_LABELED, STAT_QUEUED, STAT_CANDIDATE, STAT_WS_LINE;
239 
240  public:
242  : STAT_LABELED(1),
243  STAT_QUEUED(2),
244  STAT_CANDIDATE(ImDtTypes<T>::max()-1),
245  STAT_WS_LINE(ImDtTypes<T>::max())
246  {
247  }
248  virtual ~WatershedFlooding()
249  {
250  }
251 
252  Image<T> *imgWS;
253 #ifdef SWIG
254  const Image<T> *imgIn;
255  Image<labelT> *imgLbl;
256 #endif // SWIG
257 
258  virtual RES_T flood(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<T> &imOut, Image<labelT> &imBasinsOut, const StrElt &se)
259  {
260  ASSERT_ALLOCATED(&imIn, &imMarkers, &imBasinsOut);
261  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imBasinsOut);
262 
263  ImageFreezer freeze(imBasinsOut);
264  ImageFreezer freeze2(imOut);
265 
266  copy(imMarkers, imBasinsOut);
267 
268  initialize(imIn, imBasinsOut, imOut, se);
269  this->processImage(imIn, imBasinsOut, se);
270 
271  // Finalize the image containing the ws lines
272  // Potential remaining candidate points (points surrounded by WS_LINE points) with status STAT_CANDIDATE are added to the WS_LINE
273  for (size_t i=0;i<imIn.getPixelCount();i++)
274  {
275  if (wsPixels[i]>=STAT_CANDIDATE)
276  wsPixels[i] = STAT_WS_LINE;
277  else
278  wsPixels[i] = 0;
279  }
280  return RES_OK;
281  }
282 
283  virtual RES_T initialize(const Image<T> &imIn, Image<labelT> &imLbl, Image<T> &imOut, const StrElt &se)
284  {
286 
287  imgWS = &imOut;
288  imgWS->setSize(this->imSize);
289  test(imLbl, STAT_LABELED, STAT_CANDIDATE, *imgWS);
290  wsPixels = imgWS->getPixels();
291 
292  tmpOffsets.clear();
293 
294  return RES_OK;
295  }
296 
297  inline virtual void processPixel(const size_t &curOffset)
298  {
299  wsPixels[curOffset] = STAT_LABELED;
300 
302 
303  if (!tmpOffsets.empty())
304  {
305  if (this->wsPixels[curOffset]!=STAT_WS_LINE)
306  {
307  size_t *offsets = tmpOffsets.data();
308  for (UINT i=0;i<tmpOffsets.size();i++)
309  {
310  this->hq.push(this->inPixels[*offsets], *offsets);
311  this->wsPixels[*offsets] = STAT_QUEUED;
312 
313  offsets++;
314  }
315  }
316  tmpOffsets.clear();
317  }
318 
319  }
320  inline virtual void processNeighbor(const size_t &curOffset, const size_t &nbOffset)
321  {
322  T nbStat = this->wsPixels[nbOffset];
323 
324  if (nbStat==STAT_CANDIDATE) // Add it to the tmp offsets queue
325  {
326  tmpOffsets.push_back(nbOffset);
327  }
328  else if (nbStat==STAT_LABELED)
329  {
330  if (this->lblPixels[curOffset]==0)
331  {
332  this->lblPixels[curOffset] = this->lblPixels[nbOffset];
333  }
334  else if (this->lblPixels[curOffset]!=this->lblPixels[nbOffset])
335  this->wsPixels[curOffset] = STAT_WS_LINE;
336  }
337  }
338 #ifndef SWIG
341 #endif // SWIG
342  };
343 
356  template <class T, class labelT>
357  RES_T basins(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<labelT> &imBasinsOut, const StrElt &se=DEFAULT_SE)
358  {
359  BaseFlooding<T, labelT> flooding;
360  return flooding.flood(imIn, imMarkers, imBasinsOut, se);
361  }
362 
363  template <class T, class labelT>
364  RES_T basins(const Image<T> &imIn, Image<labelT> &imBasinsInOut, const StrElt &se=DEFAULT_SE)
365  {
366  ASSERT_ALLOCATED(&imIn);
367  ASSERT_SAME_SIZE(&imIn, &imBasinsInOut);
368 
369  Image<labelT> imLbl(imIn);
370  minimaLabeled(imIn, imLbl, se);
371 
372  return basins(imIn, imLbl, imBasinsInOut);
373  }
374 
375 
390  template <class T, class labelT>
391  RES_T watershed(const Image<T> &imIn, const Image<labelT> &imMarkers, Image<T> &imOut, Image<labelT> &imBasinsOut, const StrElt &se=DEFAULT_SE)
392  {
393  ASSERT_ALLOCATED(&imIn, &imMarkers);
394  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imOut, &imBasinsOut);
395 
397  return flooding.flood(imIn, imMarkers, imOut, imBasinsOut, se);
398  }
399 
400  template <class T, class labelT>
401  RES_T watershed(const Image<T> &imIn, Image<labelT> &imMarkers, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
402  {
403  ASSERT_ALLOCATED(&imIn, &imMarkers);
404  ASSERT_SAME_SIZE(&imIn, &imMarkers, &imOut);
405 
406  Image<labelT> imBasinsOut(imMarkers);
407  return watershed(imIn, imMarkers, imOut, imBasinsOut, se);
408  }
409 
410  template <class T>
411  RES_T watershed(const Image<T> &imIn, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
412  {
413  ASSERT_ALLOCATED(&imIn);
414  ASSERT_SAME_SIZE(&imIn, &imOut);
415 
416  Image<UINT> imLbl(imIn);
417  minimaLabeled(imIn, imLbl, se);
418  return watershed(imIn, imLbl, imOut, se);
419  }
420 
428  template <class T>
429  RES_T lblSkiz(Image<T> &labelIm1, Image<T> &labelIm2, const Image<T> &maskIm, const StrElt &se=DEFAULT_SE)
430  {
431  ASSERT_ALLOCATED(&labelIm1, &labelIm2, &maskIm);
432  ASSERT_SAME_SIZE(&labelIm1, &labelIm2, &maskIm);
433 
434  Image<T> tmpIm1(labelIm1);
435  Image<T> tmpIm2(labelIm1);
436  Image<T> sumIm(labelIm1);
437 
438  double vol1 = -1, vol2 = 0;
439 
440  ImageFreezer freeze1(labelIm1);
441  ImageFreezer freeze2(labelIm2);
442 
443  T threshMin = ImDtTypes<T>::min(), threshMax = ImDtTypes<T>::max()-T(1);
444 
445  while(vol1<vol2)
446  {
447  dilate(labelIm1, tmpIm1, se);
448  dilate(labelIm2, tmpIm2, se);
449  addNoSat(labelIm1, labelIm2, sumIm);
450  threshold(sumIm, threshMin, threshMax, sumIm);
451  // if (&maskIm) // a reference cannot be NULL
452  inf(maskIm, sumIm, sumIm);
453  mask(tmpIm1, sumIm, tmpIm1);
454  mask(tmpIm2, sumIm, tmpIm2);
455  sup(labelIm1, tmpIm1, labelIm1);
456  sup(labelIm2, tmpIm2, labelIm2);
457 
458  vol1 = vol2;
459  vol2 = vol(labelIm1);
460  }
461  return RES_OK;
462  }
463 
471  template <class T1, class T2>
472  RES_T inflBasins(const Image<T1> &imIn, Image<T2> &basinsOut, const StrElt &se=DEFAULT_SE)
473  {
474  ASSERT_ALLOCATED(&imIn, &basinsOut);
475  ASSERT_SAME_SIZE(&imIn, &basinsOut);
476 
477  Image<T2> imLbl2(basinsOut);
478  Image<T2> nullIm(basinsOut);
479  fill(nullIm, ImDtTypes<T2>::max()); // invariant to liblSkiz inf
480 
481  // Create the label images
482  label(imIn, basinsOut, se);
483  inv(basinsOut, imLbl2);
484  mask(imLbl2, basinsOut, imLbl2);
485 
486  ASSERT(lblSkiz(basinsOut, imLbl2, nullIm, se)==RES_OK);
487 
488  // Clean result image
489  open(basinsOut, basinsOut);
490 
491  return RES_OK;
492  }
493 
501  template <class T>
502  RES_T inflZones(const Image<T> &imIn, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
503  {
504  ASSERT_ALLOCATED(&imIn, &imOut);
505  ASSERT_SAME_SIZE(&imIn, &imOut);
506 
507  Image<UINT16> basinsIm(imIn);
508 
509  ASSERT(inflBasins(imIn, basinsIm, se)==RES_OK);
510  gradient(basinsIm, basinsIm, se, StrElt());
511  threshold(basinsIm, UINT16(ImDtTypes<T>::min()+T(1)), UINT16(ImDtTypes<T>::max()), basinsIm);
512  copy(basinsIm, imOut);
513 
514  return RES_OK;
515  }
516 
521  template <class T>
522  RES_T waterfall(const Image<T> &gradIn, const Image<T> &wsIn, Image<T> &imGradOut, Image<T> &imWsOut, const StrElt &se=DEFAULT_SE)
523  {
524  ASSERT_ALLOCATED(&gradIn, &wsIn, &imGradOut, &imWsOut);
525  ASSERT_SAME_SIZE(&gradIn, &wsIn, &imGradOut, &imWsOut);
526 
527  ImageFreezer freeze(imWsOut);
528 
529  test(wsIn, gradIn, ImDtTypes<T>::max(), imWsOut);
530  dualBuild(imWsOut, gradIn, imGradOut, se);
531  watershed(imGradOut, imWsOut, se);
532 
533  return RES_OK;
534  }
535 
540  template <class T>
541  RES_T waterfall(const Image<T> &gradIn, UINT nLevel, Image<T> &imWsOut, const StrElt &se=DEFAULT_SE)
542  {
543  ASSERT_ALLOCATED(&gradIn, &imWsOut);
544  ASSERT_SAME_SIZE(&gradIn, &imWsOut);
545 
546  ImageFreezer freeze(imWsOut);
547 
548  Image<T> tmpGradIm(gradIn, 1); //clone
549  Image<T> tmpGradIm2(gradIn); //clone
550 
551  watershed(gradIn, imWsOut, se);
552 
553  for (UINT i=0;i<nLevel;i++)
554  {
555  waterfall(tmpGradIm, imWsOut, tmpGradIm2, imWsOut, se);
556  copy(tmpGradIm2, tmpGradIm);
557  }
558 
559  return RES_OK;
560  }
561 
564 } // namespace smil
565 
566 
567 #endif // _D_MORPHO_WATERSHED_HPP
568 
569 
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 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
Definition: DColorConvert.h:38
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: DMeasures.hpp:96
vector< IntPoint > points
List of neighbor points.
Definition: DStructuringElement.h:125
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
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:357
RES_T inflBasins(const Image< T1 > &imIn, Image< T2 > &basinsOut, const StrElt &se=DEFAULT_SE)
Influences basins.
Definition: DMorphoWatershed.hpp:472
Definition: DMorphoWatershed.hpp:230
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:502
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:429
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:78
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:44
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
T maxVal(const Image< T > &imIn, bool onlyNonZero=false)
Max value of an image.
Definition: 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 waterfall(const Image< T > &gradIn, UINT nLevel, Image< T > &imWsOut, const StrElt &se=DEFAULT_SE)
Waterfall.
Definition: DMorphoWatershed.hpp:541
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:391
Definition: DCommon.h:104