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 
176  template <class T1, class T2>
177  RES_T dist_cross_3d (const Image<T1> &imIn, Image<T2> &imOut) {
178  ASSERT_ALLOCATED (&imIn, &imOut);
179  ASSERT_SAME_SIZE (&imIn, &imOut);
180 
181  ImageFreezer freeze (imOut);
182  Image<T1> tmp(imIn);
183  ASSERT (inf(imIn, T1(1), tmp) == RES_OK);
184 
185  typedef Image<T1> imageInType;
186  typedef typename imageInType::lineType lineInType;
187  typedef Image<T2> imageOutType;
188  typedef typename imageOutType::lineType lineOutType;
189 
190  lineInType pixelsIn = tmp.getPixels () ;
191  lineOutType pixelsOut = imOut.getPixels () ;
192 
193  int size[3];
194  imIn.getSize (size) ;
195  size_t offset ;
196  int x,y,z;
197  T2 infinite=ImDtTypes<T2>::max();
198  long int min;
199 
200  for (z=0; z<size[2]; ++z) {
201  #ifdef USE_OPEN_MP
202  #pragma omp for private(offset,x,y,min)
203  #endif // USE_OPEN_MP
204  for (x=0; x<size[0];++x) {
205  offset = z*size[1]*size[0]+x;
206  if (pixelsIn[offset] == T1(0)) {
207  pixelsOut[offset] = T2(0);
208  } else {
209  pixelsOut[offset] = infinite;
210  }
211 
212  for (y=1; y<size[1]; ++y) {
213  if (pixelsIn[offset+y*size[0]] == T1(0)) {
214  pixelsOut[offset+y*size[0]] = T2(0);
215  } else {
216  pixelsOut[offset+y*size[0]] = (1 + pixelsOut[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsOut[offset+(y-1)*size[0]];
217  }
218  }
219 
220  for (y=size[1]-2; y>=0; --y) {
221  min = (pixelsOut[offset+(y+1)*size[0]]+1 > infinite) ? infinite : pixelsOut[offset+(y+1)*size[0]]+1;
222  if (min < pixelsOut[offset+y*size[0]])
223  pixelsOut[offset+y*size[0]] = (1+pixelsOut[offset+(y+1)*size[0]]);
224  }
225  }
226 
227  #ifdef USE_OPEN_MP
228  #pragma omp for private(x,y,offset)
229  #endif // USE_OPEN_MP
230  for (y=0; y<size[1]; ++y) {
231  offset = z*size[1]*size[0]+y*size[0];
232  for (x=1; x<size[0]; ++x) {
233  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x-1]) {
234  pixelsOut[offset+x] = pixelsOut[offset+x-1]+1;
235  }
236  }
237  for (x=size[0]-2; x>=0; --x) {
238  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x+1]) {
239  pixelsOut[offset+x] = pixelsOut[offset+x+1]+1;
240  }
241  }
242  }
243  }
244  for (y=0; y<size[1]; ++y) {
245  #ifdef USE_OPEN_MP
246  #pragma omp for private(x,z,offset)
247  #endif // USE_OPEN_MP
248  for (x=0; x<size[0]; ++x) {
249  offset = y*size[0]+x;
250  for (z=1; z<size[2]; ++z) {
251  if (pixelsOut[offset+z*size[1]*size[0]] != 0 && pixelsOut[offset+z*size[1]*size[0]] > pixelsOut[offset+(z-1)*size[1]*size[0]]) {
252  pixelsOut[offset+z*size[1]*size[0]] = pixelsOut[offset+(z-1)*size[1]*size[0]]+1;
253  }
254  }
255  for (z=size[2]-2; z>=0; --z) {
256  if (pixelsOut[offset+z*size[1]*size[0]] != 0 && pixelsOut[offset+z*size[1]*size[0]] > pixelsOut[offset+(z+1)*size[1]*size[0]]) {
257  pixelsOut[offset+z*size[1]*size[0]] = pixelsOut[offset+(z+1)*size[1]*size[0]]+1;
258  }
259  }
260  }
261  }
262  return RES_OK;
263  }
264 
268  template <class T1, class T2>
269  RES_T dist_cross (const Image<T1> &imIn, Image<T2> &imOut) {
270  ASSERT_ALLOCATED (&imIn, &imOut);
271  ASSERT_SAME_SIZE (&imIn, &imOut);
272 
273  ImageFreezer freeze (imOut);
274  Image<T1> tmp(imIn);
275  ASSERT (inf(imIn, T1(1), tmp) == RES_OK);
276 
277  typedef Image<T1> imageInType;
278  typedef typename imageInType::lineType lineInType;
279  typedef Image<T2> imageOutType;
280  typedef typename imageOutType::lineType lineOutType;
281 
282  lineInType pixelsIn = tmp.getPixels () ;
283  lineOutType pixelsOut = imOut.getPixels () ;
284 
285  int size[3];
286  imIn.getSize (size) ;
287  size_t offset ;
288  int x,y,z;
289  T2 infinite=ImDtTypes<T2>::max();
290  long int min;
291 
292  for (z=0; z<size[2]; ++z) {
293  #ifdef USE_OPEN_MP
294  #pragma omp for private(offset,x,y,min)
295  #endif // USE_OPEN_MP
296  for (x=0; x<size[0];++x) {
297  offset = z*size[1]*size[0]+x;
298  if (pixelsIn[offset] == T1(0)) {
299  pixelsOut[offset] = T2(0);
300  } else {
301  pixelsOut[offset] = infinite;
302  }
303 
304  for (y=1; y<size[1]; ++y) {
305  if (pixelsIn[offset+y*size[0]] == T1(0)) {
306  pixelsOut[offset+y*size[0]] = T2(0);
307  } else {
308  pixelsOut[offset+y*size[0]] = (1 + pixelsOut[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsOut[offset+(y-1)*size[0]];
309  }
310  }
311 
312  for (y=size[1]-2; y>=0; --y) {
313  min = (pixelsOut[offset+(y+1)*size[0]]+1 > infinite) ? infinite : pixelsOut[offset+(y+1)*size[0]]+1;
314  if (min < pixelsOut[offset+y*size[0]])
315  pixelsOut[offset+y*size[0]] = (1+pixelsOut[offset+(y+1)*size[0]]);
316  }
317  }
318 
319  #ifdef USE_OPEN_MP
320  #pragma omp for private(x,y,offset)
321  #endif // USE_OPEN_MP
322  for (y=0; y<size[1]; ++y) {
323  offset = z*size[1]*size[0]+y*size[0];
324  for (x=1; x<size[0]; ++x) {
325  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x-1]) {
326  pixelsOut[offset+x] = pixelsOut[offset+x-1]+1;
327  }
328  }
329  for (x=size[0]-2; x>=0; --x) {
330  if (pixelsOut[offset+x] != 0 && pixelsOut[offset+x] > pixelsOut[offset+x+1]) {
331  pixelsOut[offset+x] = pixelsOut[offset+x+1]+1;
332  }
333  }
334  }
335  }
336  return RES_OK;
337  }
338 
342  template <class T1, class T2>
343  RES_T dist_square (const Image<T1> &imIn, Image<T2> &imOut) {
344  ASSERT_ALLOCATED (&imIn, &imOut);
345  ASSERT_SAME_SIZE (&imIn, &imOut);
346 
347  ImageFreezer freeze (imOut);
348  Image<T2> tmp(imIn);
349 
350  typedef Image<T1> imageInType;
351  typedef typename imageInType::lineType lineInType;
352  typedef Image<T2> imageOutType;
353  typedef typename imageOutType::lineType lineOutType;
354 
355  lineInType pixelsIn = imIn.getPixels () ;
356  lineOutType pixelsOut = imOut.getPixels () ;
357  lineOutType pixelsTmp = tmp.getPixels () ;
358 
359  int size[3];
360  imIn.getSize (size) ;
361  size_t offset ;
362  int x,y,z;
363  T2 infinite=ImDtTypes<T2>::max();
364  long int min;
365 
366  // H(x,u) is a minimizer, = MIN(h: 0 <= h < u & Any (i: 0 <= i < u : f(x,h) <= f(x,i)) : h )
367  long int size_array = MAX(size[0],size[1]);
368  vector<long int> s(size_array); // sets of the least minimizers that occurs during the scan from left to right.
369  vector<long int> t(size_array); // sets of points with the same least minimizer
370  s[0] = 0;
371  t[0] = 0;
372  long int q = 0;
373  long int w;
374 
375  long int tmpdist, tmpdist2;
376 
377  for (z=0; z<size[2]; ++z) {
378  #ifdef USE_OPEN_MP
379  #pragma omp for private(offset,x,y,min)
380  #endif // USE_OPEN_MP
381  for (x=0; x<size[0];++x) {
382  offset = z*size[1]*size[0]+x;
383  if (pixelsIn[offset] == T1(0)) {
384  pixelsTmp[offset] = T2(0);
385  } else {
386  pixelsTmp[offset] = infinite;
387  }
388  // SCAN 1
389  for (y=1; y<size[1]; ++y) {
390  if (pixelsIn[offset+y*size[0]] == T1(0)) {
391  pixelsTmp[offset+y*size[0]] = T2(0);
392  } else {
393  pixelsTmp[offset+y*size[0]] = (1 + pixelsTmp[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsTmp[offset+(y-1)*size[0]];
394  }
395  }
396  // SCAN 2
397  for (y=size[1]-2; y>=0; --y) {
398  min = (pixelsTmp[offset+(y+1)*size[0]]+1 > infinite) ? infinite : pixelsTmp[offset+(y+1)*size[0]]+1;
399  if (min < pixelsTmp[offset+y*size[0]])
400  pixelsTmp[offset+y*size[0]] = (1+pixelsTmp[offset+(y+1)*size[0]]);
401  }
402  }
403  #ifdef USE_OPEN_MP
404  #pragma omp for private(offset,y,q,w,tmpdist,tmpdist2)
405  #endif // USE_OPEN_MP
406  for (y=0; y<size[1]; ++y) {
407  offset = z*size[1]*size[0]+y*size[0];
408  q=0; t[0]=0; s[0]=0;
409  // SCAN 3
410  for (int u=1; u<size[0];++u) {
411  tmpdist = (t[q] > s[q]) ? t[q]-s[q] : s[q]-t[q];
412  tmpdist = (tmpdist >= pixelsTmp[offset+s[q]]) ? tmpdist : pixelsTmp[offset+s[q]];
413  tmpdist2 = (t[q] > u) ? t[q]-u : u-t[q];
414  tmpdist2 = (tmpdist2 >= pixelsTmp[offset+u]) ? tmpdist2 : pixelsTmp[offset+u];
415 
416  while (q>=0 && tmpdist > tmpdist2) {
417  q--;
418  if (q>=0) {
419  tmpdist = (t[q] > s[q]) ? t[q]-s[q] : s[q]-t[q];
420  tmpdist = (tmpdist >= pixelsTmp[offset+s[q]]) ? tmpdist : pixelsTmp[offset+s[q]];
421  tmpdist2 = (t[q] > u) ? t[q]-u : u-t[q];
422  tmpdist2 = (tmpdist2 >= pixelsTmp[offset+u]) ? tmpdist2 : pixelsTmp[offset+u];
423  }
424 
425  }
426  if (q<0) {
427  q=0; s[0]=u;}
428  else {
429  if (pixelsTmp[offset+s[q]] <= pixelsTmp[offset+u]) {
430  w = (s[q]+pixelsTmp[offset+u] >= (s[q]+u)/2) ? s[q]+pixelsTmp[offset+u] : (s[q]+u)/2;
431  } else {
432  w = (u-pixelsTmp[offset+s[q]] >= (s[q]+u)/2) ? (s[q]+u)/2 : u-pixelsTmp[offset+s[q]];
433  }
434  w=1+w;
435  if (w<size[0]) {q++; s[q]=u; t[q]=w;}
436  }
437  }
438  // SCAN 4
439  for (int u=size[0]-1; u>=0; --u) {
440  pixelsOut[offset+u] = (u > s[q]) ? u-s[q] : s[q]-u ;
441  pixelsOut[offset+u] = (pixelsOut[offset+u] >= pixelsTmp[offset+s[q]]) ? pixelsOut[offset+u] : pixelsTmp[offset+s[q]];
442  if (u == t[q])
443  q--;
444  }
445 
446  }
447  }
448  return RES_OK;
449  }
450 
454  template <class T1, class T2>
455  RES_T dist_euclidean (const Image<T1> &imIn, Image<T2> &imOut) {
456  ASSERT_ALLOCATED (&imIn, &imOut);
457  ASSERT_SAME_SIZE (&imIn, &imOut);
458 
459  ImageFreezer freeze (imOut);
460  Image<T2> tmp(imIn);
461 
462  typedef Image<T1> imageInType;
463  typedef typename imageInType::lineType lineInType;
464  typedef Image<T2> imageOutType;
465  typedef typename imageOutType::lineType lineOutType;
466 
467  lineInType pixelsIn = imIn.getPixels () ;
468  lineOutType pixelsOut = imOut.getPixels () ;
469  lineOutType pixelsTmp = tmp.getPixels () ;
470 
471  size_t size[3];
472  imIn.getSize (size) ;
473  size_t nbrPixelsPerSlice = size[0]*size[1];
474  size_t offset ;
475  size_t x,y,z;
476  T2 infinite= ImDtTypes<T2>::max();
477  // JOE long int min;
478 
479  // H(x,u) is a minimizer, = MIN(h: 0 <= h < u & Any (i: 0 <= i < u : f(x,h) <= f(x,i)) : h )
480  vector<size_t> s(size[0]); // sets of the least minimizers that occurs during the scan from left to right.
481  vector<size_t> t(size[0]); // sets of points with the same least minimizer
482  size_t q = 0;
483  size_t w;
484 
485  for (z=0; z<size[2]; ++z) {
486 // #ifdef USE_OPEN_MP
487 // #pragma omp for private(offset,x,y,min)
488 // #endif // USE_OPEN_MP
489  for (x=0; x<size[0];++x) {
490  offset = z*nbrPixelsPerSlice+x;
491  if (pixelsIn[offset] == T1(0)) {
492  pixelsOut[offset] = T2(0);
493  } else {
494  pixelsOut[offset] = infinite;
495  }
496  // SCAN 1
497  for (y=1; y<size[1]; ++y) {
498  if (pixelsIn[offset+y*size[0]] == T1(0)) {
499  pixelsOut[offset+y*size[0]] = T2(0);
500  } else {
501  pixelsOut[offset+y*size[0]] = (1 + pixelsOut[offset+(y-1)*size[0]] > infinite) ? infinite : 1 + pixelsOut[offset+(y-1)*size[0]];
502  }
503  }
504  // SCAN 2
505  for (y=size[1]-2; y>=0; --y) {
506  if (pixelsOut[offset+(y+1)*size[0]] < pixelsOut[offset+y*size[0]])
507  pixelsOut[offset+y*size[0]] = (1+pixelsOut[offset+(y+1)*size[0]] < infinite) ? 1+pixelsOut[offset+(y+1)*size[0]] : infinite;
508  }
509  }
510  }
511 
512  copy (imOut, tmp);
513 
514  #define __f_euclidean(x,i) (x-i)*(x-i)+pixelsTmp[offset+i]*pixelsTmp[offset+i]
515  #define __sep(a,b) (b*b - a*a + pixelsTmp[offset+b]*pixelsTmp[offset+b] - pixelsTmp[offset+a] * pixelsTmp[offset+a]) / (2*(b-a))
516 
517  for (z=0; z<size[2]; ++z) {
518 
519  #ifdef USE_OPEN_MP
520  #pragma omp for private(offset,y,x,q,w)
521  #endif // USE_OPEN_MP
522  for (y=0; y<size[1]; ++y) {
523  offset = z*nbrPixelsPerSlice+y*size[0];
524  q=0; t[0]=0; s[0]=0;
525  // SCAN 3
526  for (x=1; x<size[0];++x) {
527  while (q>=0 && __f_euclidean (t[q], s[q]) > __f_euclidean (t[q], x)) {
528  q--;
529  }
530  if (q<0) {q=0; s[0]=x;}
531  else {
532  w= 1 + __sep(s[q], x);
533  if (w<size[0]) {q++; s[q]=x; t[q]=w;}
534  }
535  }
536  // SCAN 4
537  for (x=size[0]-1; x>=0; --x) {
538  pixelsOut[offset+x] = __f_euclidean (x, s[q]) ;
539  if (x == t[q])
540  --q;
541  }
542  }
543 
544  }
545  #undef __f_euclidean
546  #undef __sep
547 
548  copy (imOut, tmp);
549  // The previous pixels are already squarred ...
550  #define __f_euclidean(a,i) (a-i)*(a-i)+pixelsTmp[offset+i*nbrPixelsPerSlice]
551  #define __sep(a,b) (b*b - a*a + pixelsTmp[offset+b*nbrPixelsPerSlice] - pixelsTmp[offset+a*nbrPixelsPerSlice]) / (2*(b-a))
552  for (y=0; y<size[1]; ++y) {
553 
554  #ifdef USE_OPENMP
555  #pragma omp for private(x,z, offset)
556  #endif
557  for (x=0; x<size[0]; ++x) {
558  offset = y*size[0]+x;
559  q=0; t[0]=0; s[0]=0;
560  for (z=1; z<size[2]; ++z) {
561  while (q>=0 && __f_euclidean (t[q], s[q]) > __f_euclidean (t[q], z)) {
562  q--;
563  }
564  if (q<0) {q=0; s[0]=z;}
565  else {
566 
567  w= 1 + __sep(s[q], z);
568  if (w<size[2]) {q++; s[q]=z; t[q]=w;}
569  }
570  }
571  for (z=size[2]-1; z>=0; --z) {
572  pixelsOut[offset+z*nbrPixelsPerSlice] = __f_euclidean (z, s[q]);
573  if (z == t[q]) {
574  --q;
575  }
576  }
577  }
578  }
579  #undef __f_euclidean
580  #undef __sep
581 
582  return RES_OK;
583  }
584 
588  template <class T>
589  RES_T dist_v0(const Image<T> &imIn, Image<T> &imOut, const StrElt &se=DEFAULT_SE)
590  {
591  ASSERT_ALLOCATED(&imIn, &imOut);
592  ASSERT_SAME_SIZE(&imIn, &imOut);
593 
594  ImageFreezer freeze(imOut);
595 
596  Image<T> tmpIm(imIn);
597 
598  // Set image to 1 when pixels are !=0
599  ASSERT((inf(imIn, T(1), tmpIm)==RES_OK));
600 
601  ASSERT((copy(tmpIm, imOut)==RES_OK));
602 
603  do
604  {
605  ASSERT((erode(tmpIm, tmpIm, se)==RES_OK));
606  ASSERT((add(tmpIm, imOut, imOut)==RES_OK));
607 
608  } while (vol(tmpIm)!=0);
609 
610  return RES_OK;
611  }
612 
615 } // namespace smil
616 
617 
618 
619 #endif // _D_MORPHO_DISTANCE_HPP
620 
RES_T dist_cross_3d(const Image< T1 > &imIn, Image< T2 > &imOut)
Distance Cross3D function (???).
Definition: DMorphoDistance.hpp:177
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 dist_cross(const Image< T1 > &imIn, Image< T2 > &imOut)
Distance Cross function (???).
Definition: DMorphoDistance.hpp:269
RES_T add(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Addition.
Definition: DImageArith.hpp:393
RES_T dist_euclidean(const Image< T1 > &imIn, Image< T2 > &imOut)
Euclidean Distance function.
Definition: DMorphoDistance.hpp:455
Base structuring element.
Definition: DStructuringElement.h:51
RES_T dist_square(const Image< T1 > &imIn, Image< T2 > &imOut)
Distance Square function (???).
Definition: DMorphoDistance.hpp:343
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:589
RES_T dist(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Distance function.
Definition: DMorphoDistance.hpp:52