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''
18  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS AND CONTRIBUTORS BE
21  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
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 {
50  template <class T1, class T2>
51  RES_T dist(const Image<T1> &imIn, Image<T2> &imOut,
52  const StrElt &se = DEFAULT_SE)
53  {
54  int st = se.getType();
55  int size = se.size;
56 
57  if (size > 1)
58  return distGeneric(imIn, imOut, se);
59  switch (st) {
60  case SE_Cross:
61  return distCross(imIn, imOut);
62  case SE_Cross3D:
63  return distCross3d(imIn, imOut);
64  case SE_Squ:
65  return distSquare(imIn, imOut);
66  default:
67  return distGeneric(imIn, imOut, se);
68  }
69  }
70 
74  template <class T1, class T2>
75  RES_T distGeneric(const Image<T1> &imIn, Image<T2> &imOut,
76  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  // Demi-Gradient to remove sources inside cluster of sources.
99  ASSERT(erode(tmp, tmp2, se) == RES_OK);
100  ASSERT(sub(tmp, tmp2, tmp) == RES_OK);
101 
102  ASSERT(copy(tmp, imOut) == RES_OK);
103 
104  queue<size_t> *level = new queue<size_t>();
105  queue<size_t> *next_level = new queue<size_t>();
106  queue<size_t> *swap;
107  T2 cur_level = T2(2);
108 
109  int size[3];
110  imIn.getSize(size);
111 
112  pixelsIn = imIn.getPixels();
113 
114  for (long i = 0, iMax = size[2] * size[1] * size[0]; i < iMax; ++i) {
115  if (pixelsOut[i] > T2(0)) {
116  level->push(i);
117  pixelsOut[i] = T2(1);
118  }
119  }
120 
121  size_t cur;
122  long int x, y, z, n_x, n_y, n_z;
123 
124  vector<IntPoint> sePoints = se.points;
125  vector<IntPoint>::iterator pt;
126 
127  bool oddLine;
128 
129  do {
130  while (!level->empty()) {
131  cur = level->front();
132  pt = sePoints.begin();
133 
134  z = cur / (size[1] * size[0]);
135  y = (cur - z * size[1] * size[0]) / size[0];
136  x = cur - y * size[0] - z * size[1] * size[0];
137 
138  oddLine = se.odd && (y % 2);
139 
140  while (pt != sePoints.end()) {
141  n_x = x + pt->x;
142  n_y = y + pt->y;
143  n_x += (oddLine && ((n_y + 1) % 2) != 0) ? 1 : 0;
144  n_z = z + pt->z;
145 
146  if (n_x >= 0 && n_x < (int) size[0] && n_y >= 0 &&
147  n_y < (int) size[1] && n_z >= 0 && n_z < (int) size[2] &&
148  pixelsOut[n_x + (n_y) *size[0] + (n_z) *size[1] * size[0]] ==
149  T2(0) &&
150  pixelsIn[n_x + (n_y) *size[0] + (n_z) *size[1] * size[0]] >
151  T1(0)) {
152  pixelsOut[n_x + (n_y) *size[0] + (n_z) *size[1] * size[0]] =
153  T2(cur_level);
154  next_level->push(n_x + (n_y) *size[0] + (n_z) *size[1] * size[0]);
155  }
156  ++pt;
157  }
158  level->pop();
159  }
160  ++cur_level;
161 
162  swap = level;
163  level = next_level;
164  next_level = swap;
165  } while (!level->empty());
166 
167  return RES_OK;
168  }
169 
173  template <class T1, class T2>
174  RES_T distCross3d(const Image<T1> &imIn, Image<T2> &imOut)
175  {
176  ASSERT_ALLOCATED(&imIn, &imOut);
177  ASSERT_SAME_SIZE(&imIn, &imOut);
178 
179  ImageFreezer freeze(imOut);
180  Image<T1> tmp(imIn);
181  ASSERT(inf(imIn, T1(1), tmp) == RES_OK);
182 
183  typedef Image<T1> imageInType;
184  typedef typename imageInType::lineType lineInType;
185  typedef Image<T2> imageOutType;
186  typedef typename imageOutType::lineType lineOutType;
187 
188  lineInType pixelsIn = tmp.getPixels();
189  lineOutType pixelsOut = imOut.getPixels();
190 
191  int size[3];
192  imIn.getSize(size);
193  size_t offset;
194  long int x, y, z;
195  T2 infinite = ImDtTypes<T2>::max();
196  long int min;
197 
198  for (z = 0; z < size[2]; ++z) {
199 #ifdef USE_OPEN_MP
200 #pragma omp for private(offset, x, y, min)
201 #endif // USE_OPEN_MP
202  for (x = 0; x < size[0]; ++x) {
203  offset = z * size[1] * size[0] + x;
204  if (pixelsIn[offset] == T1(0)) {
205  pixelsOut[offset] = T2(0);
206  } else {
207  pixelsOut[offset] = infinite;
208  }
209 
210  for (y = 1; y < size[1]; ++y) {
211  if (pixelsIn[offset + y * size[0]] == T1(0)) {
212  pixelsOut[offset + y * size[0]] = T2(0);
213  } else {
214  pixelsOut[offset + y * size[0]] =
215  (1 + pixelsOut[offset + (y - 1) * size[0]] > infinite)
216  ? infinite
217  : 1 + pixelsOut[offset + (y - 1) * size[0]];
218  }
219  }
220 
221  for (y = size[1] - 2; y >= 0; --y) {
222  min = (pixelsOut[offset + (y + 1) * size[0]] + 1 > infinite)
223  ? infinite
224  : pixelsOut[offset + (y + 1) * size[0]] + 1;
225  if (min < pixelsOut[offset + y * size[0]])
226  pixelsOut[offset + y * size[0]] =
227  (1 + pixelsOut[offset + (y + 1) * size[0]]);
228  }
229  }
230 
231 #ifdef USE_OPEN_MP
232 #pragma omp for private(x, y, offset)
233 #endif // USE_OPEN_MP
234  for (y = 0; y < size[1]; ++y) {
235  offset = z * size[1] * size[0] + y * size[0];
236  for (x = 1; x < size[0]; ++x) {
237  if (pixelsOut[offset + x] != 0 &&
238  pixelsOut[offset + x] > pixelsOut[offset + x - 1]) {
239  pixelsOut[offset + x] = pixelsOut[offset + x - 1] + 1;
240  }
241  }
242  for (x = size[0] - 2; x >= 0; --x) {
243  if (pixelsOut[offset + x] != 0 &&
244  pixelsOut[offset + x] > pixelsOut[offset + x + 1]) {
245  pixelsOut[offset + x] = pixelsOut[offset + x + 1] + 1;
246  }
247  }
248  }
249  }
250  for (y = 0; y < size[1]; ++y) {
251 #ifdef USE_OPEN_MP
252 #pragma omp for private(x, z, offset)
253 #endif // USE_OPEN_MP
254  for (x = 0; x < size[0]; ++x) {
255  offset = y * size[0] + x;
256  for (z = 1; z < size[2]; ++z) {
257  if (pixelsOut[offset + z * size[1] * size[0]] != 0 &&
258  pixelsOut[offset + z * size[1] * size[0]] >
259  pixelsOut[offset + (z - 1) * size[1] * size[0]]) {
260  pixelsOut[offset + z * size[1] * size[0]] =
261  pixelsOut[offset + (z - 1) * size[1] * size[0]] + 1;
262  }
263  }
264  for (z = size[2] - 2; z >= 0; --z) {
265  if (pixelsOut[offset + z * size[1] * size[0]] != 0 &&
266  pixelsOut[offset + z * size[1] * size[0]] >
267  pixelsOut[offset + (z + 1) * size[1] * size[0]]) {
268  pixelsOut[offset + z * size[1] * size[0]] =
269  pixelsOut[offset + (z + 1) * size[1] * size[0]] + 1;
270  }
271  }
272  }
273  }
274  return RES_OK;
275  }
276 
280  template <class T1, class T2>
281  RES_T distCross(const Image<T1> &imIn, Image<T2> &imOut)
282  {
283  ASSERT_ALLOCATED(&imIn, &imOut);
284  ASSERT_SAME_SIZE(&imIn, &imOut);
285 
286  ImageFreezer freeze(imOut);
287  Image<T1> tmp(imIn);
288  ASSERT(inf(imIn, T1(1), tmp) == RES_OK);
289 
290  typedef Image<T1> imageInType;
291  typedef typename imageInType::lineType lineInType;
292  typedef Image<T2> imageOutType;
293  typedef typename imageOutType::lineType lineOutType;
294 
295  lineInType pixelsIn = tmp.getPixels();
296  lineOutType pixelsOut = imOut.getPixels();
297 
298  int size[3];
299  imIn.getSize(size);
300  size_t offset;
301  long int x, y, z;
302  T2 infinite = ImDtTypes<T2>::max();
303  long int min;
304 
305  for (z = 0; z < size[2]; ++z) {
306 #ifdef USE_OPEN_MP
307 #pragma omp for private(offset, x, y, min)
308 #endif // USE_OPEN_MP
309  for (x = 0; x < size[0]; ++x) {
310  offset = z * size[1] * size[0] + x;
311  if (pixelsIn[offset] == T1(0)) {
312  pixelsOut[offset] = T2(0);
313  } else {
314  pixelsOut[offset] = infinite;
315  }
316 
317  for (y = 1; y < size[1]; ++y) {
318  if (pixelsIn[offset + y * size[0]] == T1(0)) {
319  pixelsOut[offset + y * size[0]] = T2(0);
320  } else {
321  pixelsOut[offset + y * size[0]] =
322  (1 + pixelsOut[offset + (y - 1) * size[0]] > infinite)
323  ? infinite
324  : 1 + pixelsOut[offset + (y - 1) * size[0]];
325  }
326  }
327 
328  for (y = size[1] - 2; y >= 0; --y) {
329  min = (pixelsOut[offset + (y + 1) * size[0]] + 1 > infinite)
330  ? infinite
331  : pixelsOut[offset + (y + 1) * size[0]] + 1;
332  if (min < pixelsOut[offset + y * size[0]])
333  pixelsOut[offset + y * size[0]] =
334  (1 + pixelsOut[offset + (y + 1) * size[0]]);
335  }
336  }
337 
338 #ifdef USE_OPEN_MP
339 #pragma omp for private(x, y, offset)
340 #endif // USE_OPEN_MP
341  for (y = 0; y < size[1]; ++y) {
342  offset = z * size[1] * size[0] + y * size[0];
343  for (x = 1; x < size[0]; ++x) {
344  if (pixelsOut[offset + x] != 0 &&
345  pixelsOut[offset + x] > pixelsOut[offset + x - 1]) {
346  pixelsOut[offset + x] = pixelsOut[offset + x - 1] + 1;
347  }
348  }
349  for (x = size[0] - 2; x >= 0; --x) {
350  if (pixelsOut[offset + x] != 0 &&
351  pixelsOut[offset + x] > pixelsOut[offset + x + 1]) {
352  pixelsOut[offset + x] = pixelsOut[offset + x + 1] + 1;
353  }
354  }
355  }
356  }
357  return RES_OK;
358  }
359 
363  template <class T1, class T2>
364  RES_T distSquare(const Image<T1> &imIn, Image<T2> &imOut)
365  {
366  ASSERT_ALLOCATED(&imIn, &imOut);
367  ASSERT_SAME_SIZE(&imIn, &imOut);
368 
369  ImageFreezer freeze(imOut);
370  Image<T2> tmp(imIn);
371 
372  typedef Image<T1> imageInType;
373  typedef typename imageInType::lineType lineInType;
374  typedef Image<T2> imageOutType;
375  typedef typename imageOutType::lineType lineOutType;
376 
377  lineInType pixelsIn = imIn.getPixels();
378  lineOutType pixelsOut = imOut.getPixels();
379  lineOutType pixelsTmp = tmp.getPixels();
380 
381  int size[3];
382  imIn.getSize(size);
383  size_t offset;
384  int x, y, z;
385  T2 infinite = ImDtTypes<T2>::max();
386  long int min;
387 
388  // H(x,u) is a minimizer, = MIN(h: 0 <= h < u & Any (i: 0 <= i < u : f(x,h)
389  // <= f(x,i)) : h )
390  long int size_array = MAX(size[0], size[1]);
391  vector<long int> s(size_array); // sets of the least minimizers that occurs
392  // during the scan from left to right.
393  vector<long int> t(
394  size_array); // sets of points with the same least minimizer
395  s[0] = 0;
396  t[0] = 0;
397  long int q = 0;
398  long int w;
399 
400  long int tmpdist, tmpdist2;
401 
402  for (z = 0; z < size[2]; ++z) {
403 #ifdef USE_OPEN_MP
404 #pragma omp for private(offset, x, y, min)
405 #endif // USE_OPEN_MP
406  for (x = 0; x < size[0]; ++x) {
407  offset = z * size[1] * size[0] + x;
408  if (pixelsIn[offset] == T1(0)) {
409  pixelsTmp[offset] = T2(0);
410  } else {
411  pixelsTmp[offset] = infinite;
412  }
413  // SCAN 1
414  for (y = 1; y < size[1]; ++y) {
415  if (pixelsIn[offset + y * size[0]] == T1(0)) {
416  pixelsTmp[offset + y * size[0]] = T2(0);
417  } else {
418  pixelsTmp[offset + y * size[0]] =
419  (1 + pixelsTmp[offset + (y - 1) * size[0]] > infinite)
420  ? infinite
421  : 1 + pixelsTmp[offset + (y - 1) * size[0]];
422  }
423  }
424  // SCAN 2
425  for (y = size[1] - 2; y >= 0; --y) {
426  min = (pixelsTmp[offset + (y + 1) * size[0]] + 1 > infinite)
427  ? infinite
428  : pixelsTmp[offset + (y + 1) * size[0]] + 1;
429  if (min < pixelsTmp[offset + y * size[0]])
430  pixelsTmp[offset + y * size[0]] =
431  (1 + pixelsTmp[offset + (y + 1) * size[0]]);
432  }
433  }
434 #ifdef USE_OPEN_MP
435 #pragma omp for private(offset, y, q, w, tmpdist, tmpdist2)
436 #endif // USE_OPEN_MP
437  for (y = 0; y < size[1]; ++y) {
438  offset = z * size[1] * size[0] + y * size[0];
439  q = 0;
440  t[0] = 0;
441  s[0] = 0;
442  // SCAN 3
443  for (int u = 1; u < size[0]; ++u) {
444  tmpdist = (t[q] > s[q]) ? t[q] - s[q] : s[q] - t[q];
445  tmpdist = (tmpdist >= pixelsTmp[offset + s[q]])
446  ? tmpdist
447  : pixelsTmp[offset + s[q]];
448  tmpdist2 = (t[q] > u) ? t[q] - u : u - t[q];
449  tmpdist2 = (tmpdist2 >= pixelsTmp[offset + u])
450  ? tmpdist2
451  : pixelsTmp[offset + u];
452 
453  while (q >= 0 && tmpdist > tmpdist2) {
454  q--;
455  if (q >= 0) {
456  tmpdist = (t[q] > s[q]) ? t[q] - s[q] : s[q] - t[q];
457  tmpdist = (tmpdist >= pixelsTmp[offset + s[q]])
458  ? tmpdist
459  : pixelsTmp[offset + s[q]];
460  tmpdist2 = (t[q] > u) ? t[q] - u : u - t[q];
461  tmpdist2 = (tmpdist2 >= pixelsTmp[offset + u])
462  ? tmpdist2
463  : pixelsTmp[offset + u];
464  }
465  }
466  if (q < 0) {
467  q = 0;
468  s[0] = u;
469  } else {
470  if (pixelsTmp[offset + s[q]] <= pixelsTmp[offset + u]) {
471  w = (s[q] + pixelsTmp[offset + u] >= (s[q] + u) / 2)
472  ? s[q] + pixelsTmp[offset + u]
473  : (s[q] + u) / 2;
474  } else {
475  w = (u - pixelsTmp[offset + s[q]] >= (s[q] + u) / 2)
476  ? (s[q] + u) / 2
477  : u - pixelsTmp[offset + s[q]];
478  }
479  w = 1 + w;
480  if (w < size[0]) {
481  q++;
482  s[q] = u;
483  t[q] = w;
484  }
485  }
486  }
487  // SCAN 4
488  for (int u = size[0] - 1; u >= 0; --u) {
489  pixelsOut[offset + u] = (u > s[q]) ? u - s[q] : s[q] - u;
490  pixelsOut[offset + u] =
491  (pixelsOut[offset + u] >= pixelsTmp[offset + s[q]])
492  ? pixelsOut[offset + u]
493  : pixelsTmp[offset + s[q]];
494  if (u == t[q])
495  q--;
496  }
497  }
498  }
499  return RES_OK;
500  }
501 
505  template <class T1, class T2>
506  RES_T distEuclidean(const Image<T1> &imIn, Image<T2> &imOut)
507  {
508  ASSERT_ALLOCATED(&imIn, &imOut);
509  ASSERT_SAME_SIZE(&imIn, &imOut);
510 
511  ImageFreezer freeze(imOut);
512  Image<T2> tmp(imIn);
513 
514  typedef Image<T1> imageInType;
515  typedef typename imageInType::lineType lineInType;
516  typedef Image<T2> imageOutType;
517  typedef typename imageOutType::lineType lineOutType;
518 
519  lineInType pixelsIn = imIn.getPixels();
520  lineOutType pixelsOut = imOut.getPixels();
521  lineOutType pixelsTmp = tmp.getPixels();
522 
523  int size[3];
524  imIn.getSize(size);
525  size_t nbrPixelsPerSlice = size[0] * size[1];
526  size_t offset;
527  int x, y, z;
528  T2 infinite = ImDtTypes<T2>::max();
529  // JOE long int min;
530 
531  // H(x,u) is a minimizer, = MIN(h: 0 <= h < u & Any (i: 0 <= i < u : f(x,h)
532  // <= f(x,i)) : h )
533  vector<long int> s(size[0]); // sets of the least minimizers that occurs
534  // during the scan from left to right.
535  vector<long int> t(size[0]); // sets of points with the same least minimizer
536  long int q = 0;
537  long int w;
538 
539  for (z = 0; z < size[2]; ++z) {
540  // #ifdef USE_OPEN_MP
541  // #pragma omp for private(offset,x,y,min)
542  // #endif // USE_OPEN_MP
543  for (x = 0; x < size[0]; ++x) {
544  offset = z * nbrPixelsPerSlice + x;
545  if (pixelsIn[offset] == T1(0)) {
546  pixelsOut[offset] = T2(0);
547  } else {
548  pixelsOut[offset] = infinite;
549  }
550  // SCAN 1
551  for (y = 1; y < size[1]; ++y) {
552  if (pixelsIn[offset + y * size[0]] == T1(0)) {
553  pixelsOut[offset + y * size[0]] = T2(0);
554  } else {
555  pixelsOut[offset + y * size[0]] =
556  (1 + pixelsOut[offset + (y - 1) * size[0]] > infinite)
557  ? infinite
558  : 1 + pixelsOut[offset + (y - 1) * size[0]];
559  }
560  }
561  // SCAN 2
562  for (y = size[1] - 2; y >= 0 && y < size[1]; --y) {
563  if (pixelsOut[offset + (y + 1) * size[0]] <
564  pixelsOut[offset + y * size[0]])
565  pixelsOut[offset + y * size[0]] =
566  (1 + pixelsOut[offset + (y + 1) * size[0]] < infinite)
567  ? 1 + pixelsOut[offset + (y + 1) * size[0]]
568  : infinite;
569  }
570  }
571  }
572 
573  copy(imOut, tmp);
574 
575 #define __f_euclidean(x, i) \
576  (x - i) * (x - i) + pixelsTmp[offset + i] * pixelsTmp[offset + i]
577 #define __sep(a, b) \
578  (b * b - a * a + pixelsTmp[offset + b] * pixelsTmp[offset + b] - \
579  pixelsTmp[offset + a] * pixelsTmp[offset + a]) / \
580  (2 * (b - a))
581 
582  for (z = 0; z < size[2]; ++z) {
583 #ifdef USE_OPEN_MP
584 #pragma omp for private(offset, y, x, q, w)
585 #endif // USE_OPEN_MP
586  for (y = 0; y < size[1]; ++y) {
587  offset = z * nbrPixelsPerSlice + y * size[0];
588  q = 0;
589  t[0] = 0;
590  s[0] = 0;
591  // SCAN 3
592  for (x = 1; x < size[0]; ++x) {
593  while (q >= 0 && __f_euclidean(t[q], s[q]) > __f_euclidean(t[q], x)) {
594  q--;
595  }
596  if (q < 0) {
597  q = 0;
598  s[0] = x;
599  } else {
600  w = 1 + __sep(s[q], x);
601  if (w < size[0]) {
602  q++;
603  s[q] = x;
604  t[q] = w;
605  }
606  }
607  }
608  // SCAN 4
609  for (x = size[0] - 1; x >= 0 && x < size[0]; --x) {
610  pixelsOut[offset + x] = __f_euclidean(x, s[q]);
611  if (x == t[q])
612  --q;
613  }
614  }
615  }
616 #undef __f_euclidean
617 #undef __sep
618 
619  copy(imOut, tmp);
620 // The previous pixels are already squarred ...
621 #define __f_euclidean(a, i) \
622  (a - i) * (a - i) + pixelsTmp[offset + i * nbrPixelsPerSlice]
623 #define __sep(a, b) \
624  (b * b - a * a + pixelsTmp[offset + b * nbrPixelsPerSlice] - \
625  pixelsTmp[offset + a * nbrPixelsPerSlice]) / \
626  (2 * (b - a))
627  for (y = 0; y < size[1]; ++y) {
628 #ifdef USE_OPENMP
629 #pragma omp for private(x, z, offset)
630 #endif
631  for (x = 0; x < size[0]; ++x) {
632  offset = y * size[0] + x;
633  q = 0;
634  t[0] = 0;
635  s[0] = 0;
636  for (z = 1; z < size[2]; ++z) {
637  while (q >= 0 && __f_euclidean(t[q], s[q]) > __f_euclidean(t[q], z)) {
638  q--;
639  }
640  if (q < 0) {
641  q = 0;
642  s[0] = z;
643  } else {
644  w = 1 + __sep(s[q], z);
645  if (w < size[2]) {
646  q++;
647  s[q] = z;
648  t[q] = w;
649  }
650  }
651  }
652  for (z = size[2] - 1; z >= 0 && z < size[2]; --z) {
653  pixelsOut[offset + z * nbrPixelsPerSlice] = __f_euclidean(z, s[q]);
654  if (z == t[q]) {
655  --q;
656  }
657  }
658  }
659  }
660 #undef __f_euclidean
661 #undef __sep
662 
663  return RES_OK;
664  }
665 
669  template <class T>
670  RES_T distV0(const Image<T> &imIn, Image<T> &imOut,
671  const StrElt &se = DEFAULT_SE)
672  {
673  ASSERT_ALLOCATED(&imIn, &imOut);
674  ASSERT_SAME_SIZE(&imIn, &imOut);
675 
676  ImageFreezer freeze(imOut);
677 
678  Image<T> tmpIm(imIn);
679 
680  // Set image to 1 when pixels are !=0
681  ASSERT((inf(imIn, T(1), tmpIm) == RES_OK));
682 
683  ASSERT((copy(tmpIm, imOut) == RES_OK));
684 
685  do {
686  ASSERT((erode(tmpIm, tmpIm, se) == RES_OK));
687  ASSERT((add(tmpIm, imOut, imOut) == RES_OK));
688 
689  } while (vol(tmpIm) != 0);
690 
691  return RES_OK;
692  }
693 
696 } // namespace smil
697 
698 #endif // _D_MORPHO_DISTANCE_HPP
Author Vincent Morard Date : 7 march 2011 See : Morard V. and Dokladal P. and Decenciere E...
Definition: DApplyThreshold.hpp:36
RES_T sub(const Image< T > &imIn1, const Image< T > &imIn2, Image< T > &imOut)
Subtraction.
Definition: DImageArith.hpp:446
RES_T distSquare(const Image< T1 > &imIn, Image< T2 > &imOut)
Distance Square function (???).
Definition: DMorphoDistance.hpp:364
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
RES_T distEuclidean(const Image< T1 > &imIn, Image< T2 > &imOut)
Euclidean Distance function.
Definition: DMorphoDistance.hpp:506
RES_T distGeneric(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Generic Distance function.
Definition: DMorphoDistance.hpp:75
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
RES_T distV0(const Image< T > &imIn, Image< T > &imOut, const StrElt &se=DEFAULT_SE)
Base distance function performed with successive erosions.
Definition: DMorphoDistance.hpp:670
Definition: DTypes.hpp:80
RES_T distCross(const Image< T1 > &imIn, Image< T2 > &imOut)
Distance Cross function (???).
Definition: DMorphoDistance.hpp:281
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
lineType getPixels() const
Get the pixels as a 1D array.
Definition: DImage.hpp:110
RES_T distCross3d(const Image< T1 > &imIn, Image< T2 > &imOut)
Distance Cross3D function (???).
Definition: DMorphoDistance.hpp:174
RES_T dist(const Image< T1 > &imIn, Image< T2 > &imOut, const StrElt &se=DEFAULT_SE)
Distance function.
Definition: DMorphoDistance.hpp:51