SMIL  0.9.1
DMorphoDistance.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_DISTANCE_HPP
31 #define _D_MORPHO_DISTANCE_HPP
32 
33 #include "DMorphImageOperations.hpp"
34 #include "DMorphoHierarQ.hpp"
35 #include "Base/include/private/DImageDraw.hpp"
36 #include "Base/include/private/DImageHistogram.hpp"
37 #include "Morpho/include/private/DMorphoBase.hpp"
38 
39 namespace smil
40 {
51  template <class T1, class T2>
52  RES_T dist(const Image<T1> &imIn, Image<T2> &imOut, const StrElt &se=DEFAULT_SE)
53  {
54  int st = se.getType () ;
55  int size = se.size ;
56 
57  if (size > 1)
58  return dist_generic(imIn, imOut, se);
59  switch (st)
60  {
61  case SE_Cross:
62  return dist_cross (imIn, imOut);
63  case SE_Cross3D:
64  return dist_cross_3d (imIn, imOut);
65  case SE_Squ:
66  return dist_square (imIn, imOut);
67  default:
68  return dist_generic (imIn, imOut, se);
69  }
70  }
71 
75  template <class T1, class T2>
76  RES_T dist_generic(const Image<T1> &imIn, Image<T2> &imOut, const StrElt &se=DEFAULT_SE)
77  {
78  ASSERT_ALLOCATED(&imIn, &imOut);
79  ASSERT_SAME_SIZE(&imIn, &imOut);
80 
81  ImageFreezer freeze(imOut);
82 
83  typedef Image<T1> imageInType;
84  typedef typename imageInType::lineType lineInType;
85  typedef Image<T2> imageOutType;
86  typedef typename imageOutType::lineType lineOutType;
87 
88  lineInType pixelsIn ;
89  lineOutType pixelsOut = imOut.getPixels () ;
90 
91  Image<T1> tmp(imIn);
92  Image<T1> tmp2(imIn);
93 
94  // Set image to 1 when pixels are !=0
95  ASSERT(inf(imIn, T1(1), tmp)==RES_OK);
96  ASSERT(mul(tmp, T1(255), tmp)==RES_OK);
97 
98 
99  // Demi-Gradient to remove sources inside cluster of sources.
100  ASSERT(erode (tmp, tmp2, se)==RES_OK);
101  ASSERT(sub(tmp, tmp2, tmp)==RES_OK);
102 
103  ASSERT(copy(tmp, imOut)==RES_OK);
104 
105  queue<size_t> *level = new queue<size_t>();
106  queue<size_t> *next_level = new queue<size_t>();
107  queue<size_t> *swap ;
108  T2 cur_level=T2(2);
109 
110  int size[3];
111  imIn.getSize (size) ;
112 
113  pixelsIn = imIn.getPixels ();
114 
115  int i=0;
116 
117  for (i=0; i<size[2]*size[1]*size[0]; ++i)
118  {
119  if (pixelsOut[i] > T2(0)) {
120  level->push (i);
121  pixelsOut[i] = T2(1);
122  }
123  }
124 
125  size_t cur;
126  int x,y,z, n_x, n_y, n_z;
127 
128  vector<IntPoint> sePoints = se.points;
129  vector<IntPoint>::iterator pt ;
130 
131  bool oddLine;
132 
133  do {
134  while (!level->empty()) {
135  cur = level->front();
136  pt = sePoints.begin();
137 
138  z = cur / (size[1] * size[0]);
139  y = (cur - z * size[1] * size[0]) / size[0];
140  x = cur - y *size[0] - z * size[1] * size[0];
141 
142  oddLine = se.odd && (y%2);
143 
144  while (pt!=sePoints.end()) {
145  n_x = x+pt->x;
146  n_y = y+pt->y;
147  n_x += (oddLine && ((n_y+1)%2) != 0) ? 1 : 0 ;
148  n_z = z+pt->z;
149 
150  if ( n_x >= 0 && n_x < (int)size[0] &&
151  n_y >= 0 && n_y < (int)size[1] &&
152  n_z >= 0 && n_z < (int)size[2] &&
153  pixelsOut [n_x + (n_y)*size[0] + (n_z)*size[1]*size[0]] == T2(0) &&
154  pixelsIn [n_x + (n_y)*size[0] + (n_z)*size[1]*size[0]] > T1(0))
155  {
156  pixelsOut [n_x + (n_y)*size[0] + (n_z)*size[1]*size[0]] = T2(cur_level);
157  next_level->push (n_x + (n_y)*size[0] + (n_z)*size[1]*size[0]);
158  }
159  ++pt;
160  }
161  level->pop();
162  }
163  ++cur_level;
164 
165  swap = level;
166  level = next_level;
167  next_level = swap;
168  } while (!level->empty());
169 
170  return RES_OK;
171  }
172 
173  template <class T1, class T2>
174  RES_T dist_cross_3d (const Image<T1> &imIn, Image<T2> &imOut) {
175  ASSERT_ALLOCATED (&imIn, &imOut);
176  ASSERT_SAME_SIZE (&imIn, &imOut);
177 
178  ImageFreezer freeze (imOut);
179  Image<T1> tmp(imIn);
180  ASSERT (inf(imIn, T1(1), tmp) == RES_OK);
181 
182  typedef Image<T1> imageInType;
183  typedef typename imageInType::lineType lineInType;
184  typedef Image<T2> imageOutType;
185  typedef typename imageOutType::lineType lineOutType;
186 
187  lineInType pixelsIn = tmp.getPixels () ;
188  lineOutType pixelsOut = imOut.getPixels () ;
189 
190  int size[3];
191  imIn.getSize (size) ;
192  size_t offset ;
193  int x,y,z;
194  T2 infinite=ImDtTypes<T2>::max();
195  long int min;
196 
197  for (z=0; z<size[2]; ++z) {
198  #ifdef USE_OPEN_MP
199  #pragma omp for private(offset,x,y,min)
200  #endif // USE_OPEN_MP
201  for (x=0; x<size[0];++x) {
202  offset = z*size[1]*size[0]+x;
203  if (pixelsIn[offset] == T1(0)) {
204  pixelsOut[offset] = T2(0);
205  } else {
206  pixelsOut[offset] = infinite;
207  }
208 
209  for (y=1; y<size[1]; ++y) {
210  if (pixelsIn[offset+y*size[0]] == T1(0)) {
211  pixelsOut[offset+y*size[0]] = T2(0);
212  } else {
213  pixelsOut[offset+y*size[0]] = (1 + pixelsOut[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsOut[offset+(y-1)*size[0]];
214  }
215  }
216 
217  for (y=size[1]-2; y>=0; --y) {
218  min = (pixelsOut[offset+(y+1)*size[0]]+1 > infinite) ? infinite : pixelsOut[offset+(y+1)*size[0]]+1;
219  if (min < pixelsOut[offset+y*size[0]])
220  pixelsOut[offset+y*size[0]] = (1+pixelsOut[offset+(y+1)*size[0]]);
221  }
222  }
223 
224  #ifdef USE_OPEN_MP
225  #pragma omp for private(x,y,offset)
226  #endif // USE_OPEN_MP
227  for (y=0; y<size[1]; ++y) {
228  offset = z*size[1]*size[0]+y*size[0];
229  for (x=1; x<size[0]; ++x) {
230  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x-1]) {
231  pixelsOut[offset+x] = pixelsOut[offset+x-1]+1;
232  }
233  }
234  for (x=size[0]-2; x>=0; --x) {
235  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x+1]) {
236  pixelsOut[offset+x] = pixelsOut[offset+x+1]+1;
237  }
238  }
239  }
240  }
241  for (y=0; y<size[1]; ++y) {
242  #ifdef USE_OPEN_MP
243  #pragma omp for private(x,z,offset)
244  #endif // USE_OPEN_MP
245  for (x=0; x<size[0]; ++x) {
246  offset = y*size[0]+x;
247  for (z=1; z<size[2]; ++z) {
248  if (pixelsOut[offset+z*size[1]*size[0]] != 0 && pixelsOut[offset+z*size[1]*size[0]] > pixelsOut[offset+(z-1)*size[1]*size[0]]) {
249  pixelsOut[offset+z*size[1]*size[0]] = pixelsOut[offset+(z-1)*size[1]*size[0]]+1;
250  }
251  }
252  for (z=size[2]-2; z>=0; --z) {
253  if (pixelsOut[offset+z*size[1]*size[0]] != 0 && pixelsOut[offset+z*size[1]*size[0]] > pixelsOut[offset+(z+1)*size[1]*size[0]]) {
254  pixelsOut[offset+z*size[1]*size[0]] = pixelsOut[offset+(z+1)*size[1]*size[0]]+1;
255  }
256  }
257  }
258  }
259  return RES_OK;
260  }
261 
262  template <class T1, class T2>
263  RES_T dist_cross (const Image<T1> &imIn, Image<T2> &imOut) {
264  ASSERT_ALLOCATED (&imIn, &imOut);
265  ASSERT_SAME_SIZE (&imIn, &imOut);
266 
267  ImageFreezer freeze (imOut);
268  Image<T1> tmp(imIn);
269  ASSERT (inf(imIn, T1(1), tmp) == RES_OK);
270 
271  typedef Image<T1> imageInType;
272  typedef typename imageInType::lineType lineInType;
273  typedef Image<T2> imageOutType;
274  typedef typename imageOutType::lineType lineOutType;
275 
276  lineInType pixelsIn = tmp.getPixels () ;
277  lineOutType pixelsOut = imOut.getPixels () ;
278 
279  int size[3];
280  imIn.getSize (size) ;
281  size_t offset ;
282  int x,y,z;
283  T2 infinite=ImDtTypes<T2>::max();
284  long int min;
285 
286  for (z=0; z<size[2]; ++z) {
287  #ifdef USE_OPEN_MP
288  #pragma omp for private(offset,x,y,min)
289  #endif // USE_OPEN_MP
290  for (x=0; x<size[0];++x) {
291  offset = z*size[1]*size[0]+x;
292  if (pixelsIn[offset] == T1(0)) {
293  pixelsOut[offset] = T2(0);
294  } else {
295  pixelsOut[offset] = infinite;
296  }
297 
298  for (y=1; y<size[1]; ++y) {
299  if (pixelsIn[offset+y*size[0]] == T1(0)) {
300  pixelsOut[offset+y*size[0]] = T2(0);
301  } else {
302  pixelsOut[offset+y*size[0]] = (1 + pixelsOut[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsOut[offset+(y-1)*size[0]];
303  }
304  }
305 
306  for (y=size[1]-2; y>=0; --y) {
307  min = (pixelsOut[offset+(y+1)*size[0]]+1 > infinite) ? infinite : pixelsOut[offset+(y+1)*size[0]]+1;
308  if (min < pixelsOut[offset+y*size[0]])
309  pixelsOut[offset+y*size[0]] = (1+pixelsOut[offset+(y+1)*size[0]]);
310  }
311  }
312 
313  #ifdef USE_OPEN_MP
314  #pragma omp for private(x,y,offset)
315  #endif // USE_OPEN_MP
316  for (y=0; y<size[1]; ++y) {
317  offset = z*size[1]*size[0]+y*size[0];
318  for (x=1; x<size[0]; ++x) {
319  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x-1]) {
320  pixelsOut[offset+x] = pixelsOut[offset+x-1]+1;
321  }
322  }
323  for (x=size[0]-2; x>=0; --x) {
324  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x+1]) {
325  pixelsOut[offset+x] = pixelsOut[offset+x+1]+1;
326  }
327  }
328  }
329  }
330  return RES_OK;
331  }
332 
333  template <class T1, class T2>
334  RES_T dist_square (const Image<T1> &imIn, Image<T2> &imOut) {
335  ASSERT_ALLOCATED (&imIn, &imOut);
336  ASSERT_SAME_SIZE (&imIn, &imOut);
337 
338  ImageFreezer freeze (imOut);
339  Image<T2> tmp(imIn);
340 
341  typedef Image<T1> imageInType;
342  typedef typename imageInType::lineType lineInType;
343  typedef Image<T2> imageOutType;
344  typedef typename imageOutType::lineType lineOutType;
345 
346  lineInType pixelsIn = imIn.getPixels () ;
347  lineOutType pixelsOut = imOut.getPixels () ;
348  lineOutType pixelsTmp = tmp.getPixels () ;
349 
350  int size[3];
351  imIn.getSize (size) ;
352  size_t offset ;
353  int x,y,z;
354  T2 infinite=ImDtTypes<T2>::max();
355  long int min;
356 
357  // H(x,u) is a minimizer, = MIN(h: 0 <= h < u & Any (i: 0 <= i < u : f(x,h) <= f(x,i)) : h )
358  long int size_array = MAX(size[0],size[1]);
359  vector<long int> s(size_array); // sets of the least minimizers that occurs during the scan from left to right.
360  vector<long int> t(size_array); // sets of points with the same least minimizer
361  s[0] = 0;
362  t[0] = 0;
363  long int q = 0;
364  long int w;
365 
366  long int tmpdist, tmpdist2;
367 
368  for (z=0; z<size[2]; ++z) {
369  #ifdef USE_OPEN_MP
370  #pragma omp for private(offset,x,y,min)
371  #endif // USE_OPEN_MP
372  for (x=0; x<size[0];++x) {
373  offset = z*size[1]*size[0]+x;
374  if (pixelsIn[offset] == T1(0)) {
375  pixelsTmp[offset] = T2(0);
376  } else {
377  pixelsTmp[offset] = infinite;
378  }
379  // SCAN 1
380  for (y=1; y<size[1]; ++y) {
381  if (pixelsIn[offset+y*size[0]] == T1(0)) {
382  pixelsTmp[offset+y*size[0]] = T2(0);
383  } else {
384  pixelsTmp[offset+y*size[0]] = (1 + pixelsTmp[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsTmp[offset+(y-1)*size[0]];
385  }
386  }
387  // SCAN 2
388  for (y=size[1]-2; y>=0; --y) {
389  min = (pixelsTmp[offset+(y+1)*size[0]]+1 > infinite) ? infinite : pixelsTmp[offset+(y+1)*size[0]]+1;
390  if (min < pixelsTmp[offset+y*size[0]])
391  pixelsTmp[offset+y*size[0]] = (1+pixelsTmp[offset+(y+1)*size[0]]);
392  }
393  }
394  #ifdef USE_OPEN_MP
395  #pragma omp for private(offset,y,q,w,tmpdist,tmpdist2)
396  #endif // USE_OPEN_MP
397  for (y=0; y<size[1]; ++y) {
398  offset = z*size[1]*size[0]+y*size[0];
399  q=0; t[0]=0; s[0]=0;
400  // SCAN 3
401  for (int u=1; u<size[0];++u) {
402  tmpdist = (t[q] > s[q]) ? t[q]-s[q] : s[q]-t[q];
403  tmpdist = (tmpdist >= pixelsTmp[offset+s[q]]) ? tmpdist : pixelsTmp[offset+s[q]];
404  tmpdist2 = (t[q] > u) ? t[q]-u : u-t[q];
405  tmpdist2 = (tmpdist2 >= pixelsTmp[offset+u]) ? tmpdist2 : pixelsTmp[offset+u];
406 
407  while (q>=0 && tmpdist > tmpdist2) {
408  q--;
409  if (q>=0) {
410  tmpdist = (t[q] > s[q]) ? t[q]-s[q] : s[q]-t[q];
411  tmpdist = (tmpdist >= pixelsTmp[offset+s[q]]) ? tmpdist : pixelsTmp[offset+s[q]];
412  tmpdist2 = (t[q] > u) ? t[q]-u : u-t[q];
413  tmpdist2 = (tmpdist2 >= pixelsTmp[offset+u]) ? tmpdist2 : pixelsTmp[offset+u];
414  }
415 
416  }
417  if (q<0) {
418  q=0; s[0]=u;}
419  else {
420  if (pixelsTmp[offset+s[q]] <= pixelsTmp[offset+u]) {
421  w = (s[q]+pixelsTmp[offset+u] >= (s[q]+u)/2) ? s[q]+pixelsTmp[offset+u] : (s[q]+u)/2;
422  } else {
423  w = (u-pixelsTmp[offset+s[q]] >= (s[q]+u)/2) ? (s[q]+u)/2 : u-pixelsTmp[offset+s[q]];
424  }
425  w=1+w;
426  if (w<size[0]) {q++; s[q]=u; t[q]=w;}
427  }
428  }
429  // SCAN 4
430  for (int u=size[0]-1; u>=0; --u) {
431  pixelsOut[offset+u] = (u > s[q]) ? u-s[q] : s[q]-u ;
432  pixelsOut[offset+u] = (pixelsOut[offset+u] >= pixelsTmp[offset+s[q]]) ? pixelsOut[offset+u] : pixelsTmp[offset+s[q]];
433  if (u == t[q])
434  q--;
435  }
436 
437  }
438  }
439  return RES_OK;
440  }
441 
442  template <class T1, class T2>
443  RES_T dist_euclidean (const Image<T1> &imIn, Image<T2> &imOut) {
444  ASSERT_ALLOCATED (&imIn, &imOut);
445  ASSERT_SAME_SIZE (&imIn, &imOut);
446 
447  ImageFreezer freeze (imOut);
448  Image<T2> tmp(imIn);
449 
450  typedef Image<T1> imageInType;
451  typedef typename imageInType::lineType lineInType;
452  typedef Image<T2> imageOutType;
453  typedef typename imageOutType::lineType lineOutType;
454 
455  lineInType pixelsIn = imIn.getPixels () ;
456  lineOutType pixelsOut = imOut.getPixels () ;
457  lineOutType pixelsTmp = tmp.getPixels () ;
458 
459  size_t size[3];
460  imIn.getSize (size) ;
461  size_t nbrPixelsPerSlice = size[0]*size[1];
462  size_t offset ;
463  size_t x,y,z;
464  T2 infinite= ImDtTypes<T2>::max();
465  // JOE long int min;
466 
467  // H(x,u) is a minimizer, = MIN(h: 0 <= h < u & Any (i: 0 <= i < u : f(x,h) <= f(x,i)) : h )
468  vector<size_t> s(size[0]); // sets of the least minimizers that occurs during the scan from left to right.
469  vector<size_t> t(size[0]); // sets of points with the same least minimizer
470  size_t q = 0;
471  size_t w;
472 
473  for (z=0; z<size[2]; ++z) {
474 // #ifdef USE_OPEN_MP
475 // #pragma omp for private(offset,x,y,min)
476 // #endif // USE_OPEN_MP
477  for (x=0; x<size[0];++x) {
478  offset = z*nbrPixelsPerSlice+x;
479  if (pixelsIn[offset] == T1(0)) {
480  pixelsOut[offset] = T2(0);
481  } else {
482  pixelsOut[offset] = infinite;
483  }
484  // SCAN 1
485  for (y=1; y<size[1]; ++y) {
486  if (pixelsIn[offset+y*size[0]] == T1(0)) {
487  pixelsOut[offset+y*size[0]] = T2(0);
488  } else {
489  pixelsOut[offset+y*size[0]] = (1 + pixelsOut[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsOut[offset+(y-1)*size[0]];
490  }
491  }
492  // SCAN 2
493  for (y=size[1]-2; y>=0; --y) {
494  if (pixelsOut[offset+(y+1)*size[0]] < pixelsOut[offset+y*size[0]])
495  pixelsOut[offset+y*size[0]] = (1+pixelsOut[offset+(y+1)*size[0]] < infinite) ? 1+pixelsOut[offset+(y+1)*size[0]] : infinite;
496  }
497  }
498  }
499 
500  copy (imOut, tmp);
501 
502  #define __f_euclidean(x,i) (x-i)*(x-i)+pixelsTmp[offset+i]*pixelsTmp[offset+i]
503  #define __sep(a,b) (b*b - a*a + pixelsTmp[offset+b]*pixelsTmp[offset+b] - pixelsTmp[offset+a] * pixelsTmp[offset+a]) / (2*(b-a))
504 
505  for (z=0; z<size[2]; ++z) {
506 
507  #ifdef USE_OPEN_MP
508  #pragma omp for private(offset,y,x,q,w)
509  #endif // USE_OPEN_MP
510  for (y=0; y<size[1]; ++y) {
511  offset = z*nbrPixelsPerSlice+y*size[0];
512  q=0; t[0]=0; s[0]=0;
513  // SCAN 3
514  for (x=1; x<size[0];++x) {
515  while (q>=0 && __f_euclidean (t[q], s[q]) > __f_euclidean (t[q], x)) {
516  q--;
517  }
518  if (q<0) {q=0; s[0]=x;}
519  else {
520  w= 1 + __sep(s[q], x);
521  if (w<size[0]) {q++; s[q]=x; t[q]=w;}
522  }
523  }
524  // SCAN 4
525  for (x=size[0]-1; x>=0; --x) {
526  pixelsOut[offset+x] = __f_euclidean (x, s[q]) ;
527  if (x == t[q])
528  --q;
529  }
530  }
531 
532  }
533  #undef __f_euclidean
534  #undef __sep
535 
536  copy (imOut, tmp);
537  // The previous pixels are already squarred ...
538  #define __f_euclidean(a,i) (a-i)*(a-i)+pixelsTmp[offset+i*nbrPixelsPerSlice]
539  #define __sep(a,b) (b*b - a*a + pixelsTmp[offset+b*nbrPixelsPerSlice] - pixelsTmp[offset+a*nbrPixelsPerSlice]) / (2*(b-a))
540  for (y=0; y<size[1]; ++y) {
541 
542  #ifdef USE_OPENMP
543  #pragma omp for private(x,z, offset)
544  #endif
545  for (x=0; x<size[0]; ++x) {
546  offset = y*size[0]+x;
547  q=0; t[0]=0; s[0]=0;
548  for (z=1; z<size[2]; ++z) {
549  while (q>=0 && __f_euclidean (t[q], s[q]) > __f_euclidean (t[q], z)) {
550  q--;
551  }
552  if (q<0) {q=0; s[0]=z;}
553  else {
554 
555  w= 1 + __sep(s[q], z);
556  if (w<size[2]) {q++; s[q]=z; t[q]=w;}
557  }
558  }
559  for (z=size[2]-1; z>=0; --z) {
560  pixelsOut[offset+z*nbrPixelsPerSlice] = __f_euclidean (z, s[q]);
561  if (z == t[q]) {
562  --q;
563  }
564  }
565  }
566  }
567  #undef __f_euclidean
568  #undef __sep
569 
570  return RES_OK;
571  }
572 
576  template <class T>
577  RES_T dist_v0(const Image<T> &imIn, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
578  {
579  ASSERT_ALLOCATED(&imIn, &imOut);
580  ASSERT_SAME_SIZE(&imIn, &imOut);
581 
582  ImageFreezer freeze(imOut);
583 
584  Image<T> tmpIm(imIn);
585 
586  // Set image to 1 when pixels are !=0
587  ASSERT((inf(imIn, T(1), tmpIm)==RES_OK));
588 
589  ASSERT((copy(tmpIm, imOut)==RES_OK));
590 
591  do
592  {
593  ASSERT((erode(tmpIm, tmpIm, se)==RES_OK));
594  ASSERT((add(tmpIm, imOut, imOut)==RES_OK));
595 
596  } while (vol(tmpIm)!=0);
597 
598  return RES_OK;
599  }
600 
603 } // namespace smil
604 
605 
606 
607 #endif // _D_MORPHO_DISTANCE_HPP
608 
Definition: DColorConvert.h:38
RES_T sub(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Subtraction.
Definition: DImageArith.hpp:446
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
RES_T dist_generic(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Generic Distance function.
Definition: DMorphoDistance.hpp:76
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
RES_T add(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Addition.
Definition: DImageArith.hpp:393
Base structuring element.
Definition: DStructuringElement.h:51
Definition: DTypes.hpp:78
RES_T erode(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE, T borderVal=ImDtTypes< T >::max())
Morphological grayscale erosion.
Definition: DMorphoBase.hpp:101
RES_T inf(const Image< T > &imIn1, const T &value, Image< T > &imOut)
Inf of two images.
Definition: DImageArith.hpp:538
RES_T mul(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Multiply.
Definition: DImageArith.hpp:770
RES_T dist_v0(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Base distance function performed with successive erosions.
Definition: DMorphoDistance.hpp:577
RES_T dist(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Distance function.
Definition: DMorphoDistance.hpp:52