1 /*====================================================================*
2  -  Copyright (C) 2001 Leptonica.  All rights reserved.
3  -  This software is distributed in the hope that it will be
4  -  useful, but with NO WARRANTY OF ANY KIND.
5  -  No author or distributor accepts responsibility to anyone for the
6  -  consequences of using this software, or for whether it serves any
7  -  particular purpose or works at all, unless he or she says so in
8  -  writing.  Everyone is granted permission to copy, modify and
9  -  redistribute this source code, for commercial or non-commercial
10  -  purposes, with the following restrictions: (1) the origin of this
11  -  source code must not be misrepresented; (2) modified versions must
12  -  be plainly marked as such; and (3) this notice may not be removed
13  -  or altered from any source or modified source distribution.
14  *====================================================================*/
15 
16 /*
17  *  seedfilllow.c
18  *
19  *      Seedfill:
20  *      Gray seedfill (source: Luc Vincent:fast-hybrid-grayscale-reconstruction)
21  *               void   seedfillBinaryLow()
22  *               void   seedfillGrayLow()
23  *               void   seedfillGrayInvLow()
24  *               void   seedfillGrayLowSimple()
25  *               void   seedfillGrayInvLowSimple()
26  *
27  *      Distance function:
28  *               void   distanceFunctionLow()
29  *
30  *      Seed spread:
31  *               void   seedspreadLow()
32  *
33  */
34 
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <math.h>
38 #include "allheaders.h"
39 
40 struct L_Pixel
41 {
42     l_int32    x;
43     l_int32    y;
44 };
45 typedef struct L_Pixel  L_PIXEL;
46 
47 
48 /*-----------------------------------------------------------------------*
49  *                 Vincent's Iterative Binary Seedfill                   *
50  *-----------------------------------------------------------------------*/
51 /*!
52  *  seedfillBinaryLow()
53  *
54  *  Notes:
55  *      (1) This is an in-place fill, where the seed image is
56  *          filled, clipping to the filling mask, in one full
57  *          cycle of UL -> LR and LR -> UL raster scans.
58  *      (2) Assume the mask is a filling mask, not a blocking mask.
59  *      (3) Assume that the RHS pad bits of the mask
60  *          are properly set to 0.
61  *      (4) Clip to the smallest dimensions to avoid invalid reads.
62  */
63 void
seedfillBinaryLow(l_uint32 * datas,l_int32 hs,l_int32 wpls,l_uint32 * datam,l_int32 hm,l_int32 wplm,l_int32 connectivity)64 seedfillBinaryLow(l_uint32  *datas,
65                   l_int32    hs,
66                   l_int32    wpls,
67                   l_uint32  *datam,
68                   l_int32    hm,
69                   l_int32    wplm,
70                   l_int32    connectivity)
71 {
72 l_int32    i, j, h, wpl;
73 l_uint32   word, mask;
74 l_uint32   wordabove, wordleft, wordbelow, wordright;
75 l_uint32   wordprev;  /* test against this in previous iteration */
76 l_uint32  *lines, *linem;
77 
78     PROCNAME("seedfillBinaryLow");
79 
80     h = L_MIN(hs, hm);
81     wpl = L_MIN(wpls, wplm);
82 
83     switch (connectivity)
84     {
85     case 4:
86             /* UL --> LR scan */
87         for (i = 0; i < h; i++) {
88             lines = datas + i * wpls;
89             linem = datam + i * wplm;
90             for (j = 0; j < wpl; j++) {
91                 word = *(lines + j);
92                 mask = *(linem + j);
93 
94                     /* OR from word above and from word to left; mask */
95                 if (i > 0) {
96                     wordabove = *(lines - wpls + j);
97                     word |= wordabove;
98                 }
99                 if (j > 0) {
100                     wordleft = *(lines + j - 1);
101                     word |= wordleft << 31;
102                 }
103                 word &= mask;
104 
105                     /* No need to fill horizontally? */
106                 if (!word || !(~word)) {
107                     *(lines + j) = word;
108                     continue;
109                 }
110 
111                 while (1) {
112                     wordprev = word;
113                     word = (word | (word >> 1) | (word << 1)) & mask;
114                     if ((word ^ wordprev) == 0) {
115                         *(lines + j) = word;
116                         break;
117                     }
118                 }
119             }
120         }
121 
122             /* LR --> UL scan */
123         for (i = h - 1; i >= 0; i--) {
124             lines = datas + i * wpls;
125             linem = datam + i * wplm;
126             for (j = wpl - 1; j >= 0; j--) {
127                 word = *(lines + j);
128                 mask = *(linem + j);
129 
130                     /* OR from word below and from word to right; mask */
131                 if (i < h - 1) {
132                     wordbelow = *(lines + wpls + j);
133                     word |= wordbelow;
134                 }
135                 if (j < wpl - 1) {
136                     wordright = *(lines + j + 1);
137                     word |= wordright >> 31;
138                 }
139                 word &= mask;
140 
141                     /* No need to fill horizontally? */
142                 if (!word || !(~word)) {
143                     *(lines + j) = word;
144                     continue;
145                 }
146 
147                 while (1) {
148                     wordprev = word;
149                     word = (word | (word >> 1) | (word << 1)) & mask;
150                     if ((word ^ wordprev) == 0) {
151                         *(lines + j) = word;
152                         break;
153                     }
154                 }
155             }
156         }
157         break;
158 
159     case 8:
160             /* UL --> LR scan */
161         for (i = 0; i < h; i++) {
162             lines = datas + i * wpls;
163             linem = datam + i * wplm;
164             for (j = 0; j < wpl; j++) {
165                 word = *(lines + j);
166                 mask = *(linem + j);
167 
168                     /* OR from words above and from word to left; mask */
169                 if (i > 0) {
170                     wordabove = *(lines - wpls + j);
171                     word |= (wordabove | (wordabove << 1) | (wordabove >> 1));
172                     if (j > 0)
173                         word |= (*(lines - wpls + j - 1)) << 31;
174                     if (j < wpl - 1)
175                         word |= (*(lines - wpls + j + 1)) >> 31;
176                 }
177                 if (j > 0) {
178                     wordleft = *(lines + j - 1);
179                     word |= wordleft << 31;
180                 }
181                 word &= mask;
182 
183                     /* No need to fill horizontally? */
184                 if (!word || !(~word)) {
185                     *(lines + j) = word;
186                     continue;
187                 }
188 
189                 while (1) {
190                     wordprev = word;
191                     word = (word | (word >> 1) | (word << 1)) & mask;
192                     if ((word ^ wordprev) == 0) {
193                         *(lines + j) = word;
194                         break;
195                     }
196                 }
197             }
198         }
199 
200             /* LR --> UL scan */
201         for (i = h - 1; i >= 0; i--) {
202             lines = datas + i * wpls;
203             linem = datam + i * wplm;
204             for (j = wpl - 1; j >= 0; j--) {
205                 word = *(lines + j);
206                 mask = *(linem + j);
207 
208                     /* OR from words below and from word to right; mask */
209                 if (i < h - 1) {
210                     wordbelow = *(lines + wpls + j);
211                     word |= (wordbelow | (wordbelow << 1) | (wordbelow >> 1));
212                     if (j > 0)
213                         word |= (*(lines + wpls + j - 1)) << 31;
214                     if (j < wpl - 1)
215                         word |= (*(lines + wpls + j + 1)) >> 31;
216                 }
217                 if (j < wpl - 1) {
218                     wordright = *(lines + j + 1);
219                     word |= wordright >> 31;
220                 }
221                 word &= mask;
222 
223                     /* No need to fill horizontally? */
224                 if (!word || !(~word)) {
225                     *(lines + j) = word;
226                     continue;
227                 }
228 
229                 while (1) {
230                     wordprev = word;
231                     word = (word | (word >> 1) | (word << 1)) & mask;
232                     if ((word ^ wordprev) == 0) {
233                         *(lines + j) = word;
234                         break;
235                     }
236                 }
237             }
238         }
239         break;
240 
241     default:
242         L_ERROR("connectivity must be 4 or 8", procName);
243     }
244 
245     return;
246 }
247 
248 
249 
250 /*-----------------------------------------------------------------------*
251  *                 Vincent's Hybrid Grayscale Seedfill                *
252  *-----------------------------------------------------------------------*/
253 /*!
254  *  seedfillGrayLow()
255  *
256  *  Notes:
257  *      (1) The pixels are numbered as follows:
258  *              1  2  3
259  *              4  x  5
260  *              6  7  8
261  *          This low-level filling operation consists of two scans,
262  *          raster and anti-raster, covering the entire seed image.
263  *          This is followed by a breadth-first propagation operation to
264  *          complete the fill.
265  *          During the anti-raster scan, every pixel p whose current value
266  *          could still be propagated after the anti-raster scan is put into
267  *          the FIFO queue.
268  *          The propagation step is a breadth-first fill to completion.
269  *          Unlike the simple grayscale seedfill pixSeedfillGraySimple(),
270  *          where at least two full raster/anti-raster iterations are required
271  *          for completion and verification, the hybrid method uses only a
272  *          single raster/anti-raster set of scans.
273  *      (2) The filling action can be visualized from the following example.
274  *          Suppose the mask, which clips the fill, is a sombrero-shaped
275  *          surface, where the highest point is 200 and the low pixels
276  *          around the rim are 30.  Beyond the rim, the mask goes up a bit.
277  *          Suppose the seed, which is filled, consists of a single point
278  *          of height 150, located below the max of the mask, with
279  *          the rest 0.  Then in the raster scan, nothing happens until
280  *          the high seed point is encountered, and then this value is
281  *          propagated right and down, until it hits the side of the
282  *          sombrero.   The seed can never exceed the mask, so it fills
283  *          to the rim, going lower along the mask surface.  When it
284  *          passes the rim, the seed continues to fill at the rim
285  *          height to the edge of the seed image.  Then on the
286  *          anti-raster scan, the seed fills flat inside the
287  *          sombrero to the upper and left, and then out from the
288  *          rim as before.  The final result has a seed that is
289  *          flat outside the rim, and inside it fills the sombrero
290  *          but only up to 150.  If the rim height varies, the
291  *          filled seed outside the rim will be at the highest
292  *          point on the rim, which is a saddle point on the rim.
293  *      (3) Reference paper :
294  *            L. Vincent, Morphological grayscale reconstruction in image
295  *            analysis: applications and efficient algorithms, IEEE Transactions
296  *            on  Image Processing, vol. 2, no. 2, pp. 176-201, 1993.
297  */
298 void
seedfillGrayLow(l_uint32 * datas,l_int32 w,l_int32 h,l_int32 wpls,l_uint32 * datam,l_int32 wplm,l_int32 connectivity)299 seedfillGrayLow(l_uint32  *datas,
300                 l_int32    w,
301                 l_int32    h,
302                 l_int32    wpls,
303                 l_uint32  *datam,
304                 l_int32    wplm,
305                 l_int32    connectivity)
306 {
307 l_uint8    val1, val2, val3, val4, val5, val6, val7, val8;
308 l_uint8    val, maxval, maskval, boolval;
309 l_int32    i, j, imax, jmax, queue_size;
310 l_uint32  *lines, *linem;
311 L_PIXEL *pixel;
312 L_QUEUE  *lq_pixel;
313 
314     PROCNAME("seedfillGrayLow");
315 
316     imax = h - 1;
317     jmax = w - 1;
318 
319         /* In the worst case, most of the pixels could be pushed
320          * onto the FIFO queue during anti-raster scan.  However this
321          * will rarely happen, and we initialize the queue ptr size to
322          * the image perimeter. */
323     lq_pixel = lqueueCreate(2 * (w + h));
324 
325     switch (connectivity)
326     {
327     case 4:
328             /* UL --> LR scan  (Raster Order)
329              * If I : mask image
330              *    J : marker image
331              * Let p be the currect pixel;
332              * J(p) <- (max{J(p) union J(p) neighbors in raster order})
333              *          intersection I(p) */
334         for (i = 0; i < h; i++) {
335             lines = datas + i * wpls;
336             linem = datam + i * wplm;
337             for (j = 0; j < w; j++) {
338                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
339                     maxval = 0;
340                     if (i > 0)
341                         maxval = GET_DATA_BYTE(lines - wpls, j);
342                     if (j > 0) {
343                         val4 = GET_DATA_BYTE(lines, j - 1);
344                         maxval = L_MAX(maxval, val4);
345                     }
346                     val = GET_DATA_BYTE(lines, j);
347                     maxval = L_MAX(maxval, val);
348                     val = L_MIN(maxval, maskval);
349                     SET_DATA_BYTE(lines, j, val);
350                 }
351             }
352         }
353 
354             /* LR --> UL scan (anti-raster order)
355              * Let p be the currect pixel;
356              * J(p) <- (max{J(p) union J(p) neighbors in anti-raster order})
357              *          intersection I(p) */
358         for (i = imax; i >= 0; i--) {
359             lines = datas + i * wpls;
360             linem = datam + i * wplm;
361             for (j = jmax; j >= 0; j--) {
362                 boolval = FALSE;
363                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
364                     maxval = 0;
365                     if (i < imax)
366                         maxval = GET_DATA_BYTE(lines + wpls, j);
367                     if (j < jmax) {
368                         val5 = GET_DATA_BYTE(lines, j + 1);
369                         maxval = L_MAX(maxval, val5);
370                     }
371                     val = GET_DATA_BYTE(lines, j);
372                     maxval = L_MAX(maxval, val);
373                     val = L_MIN(maxval, maskval);
374                     SET_DATA_BYTE(lines, j, val);
375 
376                         /*
377                          * If there exists a point (q) which belongs to J(p)
378                          * neighbors in anti-raster order such that J(q) < J(p)
379                          * and J(q) < I(q) then
380                          * fifo_add(p) */
381                     if (i < imax) {
382                         val7 = GET_DATA_BYTE(lines + wpls, j);
383                         if ((val7 < val) &&
384                             (val7 < GET_DATA_BYTE(linem + wplm, j))) {
385                             boolval = TRUE;
386                         }
387                     }
388                     if (j < jmax) {
389                         val5 = GET_DATA_BYTE(lines, j + 1);
390                         if (!boolval && (val5 < val) &&
391                             (val5 < GET_DATA_BYTE(linem, j + 1))) {
392                             boolval = TRUE;
393                         }
394                     }
395                     if (boolval) {
396                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
397                         pixel->x = i;
398                         pixel->y = j;
399                         lqueueAdd(lq_pixel, pixel);
400                     }
401                 }
402             }
403         }
404 
405             /* Propagation step:
406              *        while fifo_empty = false
407              *          p <- fifo_first()
408              *          for every pixel (q) belong to neighbors of (p)
409              *            if J(q) < J(p) and I(q) != J(q)
410              *              J(q) <- min(J(p), I(q));
411              *              fifo_add(q);
412              *            end
413              *          end
414              *        end */
415         queue_size = lqueueGetCount(lq_pixel);
416         while (queue_size) {
417             pixel = (L_PIXEL *)lqueueRemove(lq_pixel);
418             i = pixel->x;
419             j = pixel->y;
420             FREE(pixel);
421             lines = datas + i * wpls;
422             linem = datam + i * wplm;
423 
424             if ((val = GET_DATA_BYTE(lines, j)) > 0) {
425                 if (i > 0) {
426                     val2 = GET_DATA_BYTE(lines - wpls, j);
427                     maskval = GET_DATA_BYTE(linem - wplm, j);
428                     if (val > val2 && val2 != maskval) {
429                         SET_DATA_BYTE(lines - wpls, j, L_MIN(val, maskval));
430                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
431                         pixel->x = i - 1;
432                         pixel->y = j;
433                         lqueueAdd(lq_pixel, pixel);
434                     }
435 
436                 }
437                 if (j > 0) {
438                     val4 = GET_DATA_BYTE(lines, j - 1);
439                     maskval = GET_DATA_BYTE(linem, j - 1);
440                     if (val > val4 && val4 != maskval) {
441                         SET_DATA_BYTE(lines, j - 1, L_MIN(val, maskval));
442                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
443                         pixel->x = i;
444                         pixel->y = j - 1;
445                         lqueueAdd(lq_pixel, pixel);
446                     }
447                 }
448                 if (i < imax) {
449                     val7 = GET_DATA_BYTE(lines + wpls, j);
450                     maskval = GET_DATA_BYTE(linem + wplm, j);
451                     if (val > val7 && val7 != maskval) {
452                         SET_DATA_BYTE(lines + wpls, j, L_MIN(val, maskval));
453                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
454                         pixel->x = i + 1;
455                         pixel->y = j;
456                         lqueueAdd(lq_pixel, pixel);
457                     }
458                 }
459                 if (j < jmax) {
460                     val5 = GET_DATA_BYTE(lines, j + 1);
461                     maskval = GET_DATA_BYTE(linem, j + 1);
462                     if (val > val5 && val5 != maskval) {
463                         SET_DATA_BYTE(lines, j + 1, L_MIN(val, maskval));
464                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
465                         pixel->x = i;
466                         pixel->y = j + 1;
467                         lqueueAdd(lq_pixel, pixel);
468                     }
469                 }
470             }
471 
472             queue_size = lqueueGetCount(lq_pixel);
473         }
474 
475         break;
476 
477     case 8:
478             /* UL --> LR scan  (Raster Order)
479              * If I : mask image
480              *    J : marker image
481              * Let p be the currect pixel;
482              * J(p) <- (max{J(p) union J(p) neighbors in raster order})
483              *          intersection I(p) */
484         for (i = 0; i < h; i++) {
485             lines = datas + i * wpls;
486             linem = datam + i * wplm;
487             for (j = 0; j < w; j++) {
488                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
489                     maxval = 0;
490                     if (i > 0) {
491                         if (j > 0)
492                             maxval = GET_DATA_BYTE(lines - wpls, j - 1);
493                         if (j < jmax) {
494                             val3 = GET_DATA_BYTE(lines - wpls, j + 1);
495                             maxval = L_MAX(maxval, val3);
496                         }
497                         val2 = GET_DATA_BYTE(lines - wpls, j);
498                         maxval = L_MAX(maxval, val2);
499                     }
500                     if (j > 0) {
501                         val4 = GET_DATA_BYTE(lines, j - 1);
502                         maxval = L_MAX(maxval, val4);
503                     }
504                     val = GET_DATA_BYTE(lines, j);
505                     maxval = L_MAX(maxval, val);
506                     val = L_MIN(maxval, maskval);
507                     SET_DATA_BYTE(lines, j, val);
508                 }
509             }
510         }
511 
512             /* LR --> UL scan (anti-raster order)
513              * Let p be the currect pixel;
514              * J(p) <- (max{J(p) union J(p) neighbors in anti-raster order})
515              *          intersection I(p) */
516         for (i = imax; i >= 0; i--) {
517             lines = datas + i * wpls;
518             linem = datam + i * wplm;
519             for (j = jmax; j >= 0; j--) {
520                 boolval = FALSE;
521                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
522                     maxval = 0;
523                     if (i < imax) {
524                         if (j > 0) {
525                             maxval = GET_DATA_BYTE(lines + wpls, j - 1);
526                         }
527                         if (j < jmax) {
528                             val8 = GET_DATA_BYTE(lines + wpls, j + 1);
529                             maxval = L_MAX(maxval, val8);
530                         }
531                         val7 = GET_DATA_BYTE(lines + wpls, j);
532                         maxval = L_MAX(maxval, val7);
533                     }
534                     if (j < jmax) {
535                         val5 = GET_DATA_BYTE(lines, j + 1);
536                         maxval = L_MAX(maxval, val5);
537                     }
538                     val = GET_DATA_BYTE(lines, j);
539                     maxval = L_MAX(maxval, val);
540                     val = L_MIN(maxval, maskval);
541                     SET_DATA_BYTE(lines, j, val);
542 
543                         /* If there exists a point (q) which belongs to J(p)
544                          * neighbors in anti-raster order such that J(q) < J(p)
545                          * and J(q) < I(q) then
546                          * fifo_add(p) */
547                     if (i < imax) {
548                         if (j > 0) {
549                             val6 = GET_DATA_BYTE(lines + wpls, j - 1);
550                             if ((val6 < val) &&
551                                 (val6 < GET_DATA_BYTE(linem + wplm, j - 1))) {
552                                 boolval = TRUE;
553                             }
554                         }
555                         if (j < jmax) {
556                             val8 = GET_DATA_BYTE(lines + wpls, j + 1);
557                             if (!boolval && (val8 < val) &&
558                                 (val8 < GET_DATA_BYTE(linem + wplm, j + 1))) {
559                                 boolval = TRUE;
560                             }
561                         }
562                         val7 = GET_DATA_BYTE(lines + wpls, j);
563                         if (!boolval && (val7 < val) &&
564                             (val7 < GET_DATA_BYTE(linem + wplm, j))) {
565                             boolval = TRUE;
566                         }
567                     }
568                     if (j < jmax) {
569                         val5 = GET_DATA_BYTE(lines, j + 1);
570                         if (!boolval && (val5 < val) &&
571                             (val5 < GET_DATA_BYTE(linem, j + 1))) {
572                             boolval = TRUE;
573                         }
574                     }
575                     if (boolval) {
576                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
577                         pixel->x = i;
578                         pixel->y = j;
579                         lqueueAdd(lq_pixel, pixel);
580                     }
581                 }
582             }
583         }
584 
585             /* Propagation step:
586              *        while fifo_empty = false
587              *          p <- fifo_first()
588              *          for every pixel (q) belong to neighbors of (p)
589              *            if J(q) < J(p) and I(q) != J(q)
590              *              J(q) <- min(J(p), I(q));
591              *              fifo_add(q);
592              *            end
593              *          end
594              *        end */
595         queue_size = lqueueGetCount(lq_pixel);
596         while (queue_size) {
597             pixel = (L_PIXEL *)lqueueRemove(lq_pixel);
598             i = pixel->x;
599             j = pixel->y;
600             FREE(pixel);
601             lines = datas + i * wpls;
602             linem = datam + i * wplm;
603 
604             if ((val = GET_DATA_BYTE(lines, j)) > 0) {
605                 if (i > 0) {
606                     if (j > 0) {
607                         val1 = GET_DATA_BYTE(lines - wpls, j - 1);
608                         maskval = GET_DATA_BYTE(linem - wplm, j - 1);
609                         if (val > val1 && val1 != maskval) {
610                             SET_DATA_BYTE(lines - wpls, j - 1, L_MIN(val, maskval));
611                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
612                             pixel->x = i - 1;
613                             pixel->y = j - 1;
614                             lqueueAdd(lq_pixel, pixel);
615                         }
616                     }
617                     if (j < jmax) {
618                         val3 = GET_DATA_BYTE(lines - wpls, j + 1);
619                         maskval = GET_DATA_BYTE(linem - wplm, j + 1);
620                         if (val > val3 && val3 != maskval) {
621                             SET_DATA_BYTE(lines - wpls, j + 1, L_MIN(val, maskval));
622                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
623                             pixel->x = i - 1;
624                             pixel->y = j + 1;
625                             lqueueAdd(lq_pixel, pixel);
626                         }
627                     }
628                     val2 = GET_DATA_BYTE(lines - wpls, j);
629                     maskval = GET_DATA_BYTE(linem - wplm, j);
630                     if (val > val2 && val2 != maskval) {
631                         SET_DATA_BYTE(lines - wpls, j, L_MIN(val, maskval));
632                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
633                         pixel->x = i - 1;
634                         pixel->y = j;
635                         lqueueAdd(lq_pixel, pixel);
636                     }
637 
638                 }
639                 if (j > 0) {
640                     val4 = GET_DATA_BYTE(lines, j - 1);
641                     maskval = GET_DATA_BYTE(linem, j - 1);
642                     if (val > val4 && val4 != maskval) {
643                         SET_DATA_BYTE(lines, j - 1, L_MIN(val, maskval));
644                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
645                         pixel->x = i;
646                         pixel->y = j - 1;
647                         lqueueAdd(lq_pixel, pixel);
648                     }
649                 }
650                 if (i < imax) {
651                     if (j > 0) {
652                         val6 = GET_DATA_BYTE(lines + wpls, j - 1);
653                         maskval = GET_DATA_BYTE(linem + wplm, j - 1);
654                         if (val > val6 && val6 != maskval) {
655                             SET_DATA_BYTE(lines + wpls, j - 1, L_MIN(val, maskval));
656                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
657                             pixel->x = i + 1;
658                             pixel->y = j - 1;
659                             lqueueAdd(lq_pixel, pixel);
660                         }
661                     }
662                     if (j < jmax) {
663                         val8 = GET_DATA_BYTE(lines + wpls, j + 1);
664                         maskval = GET_DATA_BYTE(linem + wplm, j + 1);
665                         if (val > val8 && val8 != maskval) {
666                             SET_DATA_BYTE(lines + wpls, j + 1, L_MIN(val, maskval));
667                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
668                             pixel->x = i + 1;
669                             pixel->y = j + 1;
670                             lqueueAdd(lq_pixel, pixel);
671                         }
672                     }
673                     val7 = GET_DATA_BYTE(lines + wpls, j);
674                     maskval = GET_DATA_BYTE(linem + wplm, j);
675                     if (val > val7 && val7 != maskval) {
676                         SET_DATA_BYTE(lines + wpls, j, L_MIN(val, maskval));
677                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
678                         pixel->x = i + 1;
679                         pixel->y = j;
680                         lqueueAdd(lq_pixel, pixel);
681                     }
682                 }
683                 if (j < jmax) {
684                     val5 = GET_DATA_BYTE(lines, j + 1);
685                     maskval = GET_DATA_BYTE(linem, j + 1);
686                     if (val > val5 && val5 != maskval) {
687                         SET_DATA_BYTE(lines, j + 1, L_MIN(val, maskval));
688                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
689                         pixel->x = i;
690                         pixel->y = j + 1;
691                         lqueueAdd(lq_pixel, pixel);
692                     }
693                 }
694             }
695 
696             queue_size = lqueueGetCount(lq_pixel);
697         }
698         break;
699 
700     default:
701         L_ERROR("connectivity must be 4 or 8", procName);
702         lqueueDestroy(&lq_pixel, TRUE);
703     }
704 
705     lqueueDestroy(&lq_pixel, TRUE);
706     return;
707 }
708 
709 
710 /*!
711  *  seedfillGrayInvLow()
712  *
713  *  Notes:
714  *      (1) The pixels are numbered as follows:
715  *              1  2  3
716  *              4  x  5
717  *              6  7  8
718  *          This low-level filling operation consists of two scans,
719  *          raster and anti-raster, covering the entire seed image.
720  *          During the anti-raster scan, every pixel p such that its
721  *          current value could still be propogated during the next
722  *          raster scanning is put into the FIFO-queue.
723  *          Next step is the propagation step where where we update
724  *          and propagate the values using FIFO structure created in
725  *          anti-raster scan.
726  *      (2) The "Inv" signifies the fact that in this case, filling
727  *          of the seed only takes place when the seed value is
728  *          greater than the mask value.  The mask will act to stop
729  *          the fill when it is higher than the seed level.  (This is
730  *          in contrast to conventional grayscale filling where the
731  *          seed always fills below the mask.)
732  *      (3) An example of use is a basin, described by the mask (pixm),
733  *          where within the basin, the seed pix (pixs) gets filled to the
734  *          height of the highest seed pixel that is above its
735  *          corresponding max pixel.  Filling occurs while the
736  *          propagating seed pixels in pixs are larger than the
737  *          corresponding mask values in pixm.
738  *      (4) Reference paper :
739  *            L. Vincent, Morphological grayscale reconstruction in image
740  *            analysis: applications and efficient algorithms, IEEE Transactions
741  *            on  Image Processing, vol. 2, no. 2, pp. 176-201, 1993.
742  */
743 void
seedfillGrayInvLow(l_uint32 * datas,l_int32 w,l_int32 h,l_int32 wpls,l_uint32 * datam,l_int32 wplm,l_int32 connectivity)744 seedfillGrayInvLow(l_uint32  *datas,
745                    l_int32    w,
746                    l_int32    h,
747                    l_int32    wpls,
748                    l_uint32  *datam,
749                    l_int32    wplm,
750                    l_int32    connectivity)
751 {
752 l_uint8    val1, val2, val3, val4, val5, val6, val7, val8;
753 l_uint8    val, maxval, maskval, boolval;
754 l_int32    i, j, imax, jmax, queue_size;
755 l_uint32  *lines, *linem;
756 L_PIXEL *pixel;
757 L_QUEUE  *lq_pixel;
758 
759     PROCNAME("seedfillGrayInvLow");
760 
761     imax = h - 1;
762     jmax = w - 1;
763 
764         /* In the worst case, most of the pixels could be pushed
765          * onto the FIFO queue during anti-raster scan.  However this
766          * will rarely happen, and we initialize the queue ptr size to
767          * the image perimeter. */
768     lq_pixel = lqueueCreate(2 * (w + h));
769 
770     switch (connectivity)
771     {
772     case 4:
773             /* UL --> LR scan  (Raster Order)
774              * If I : mask image
775              *    J : marker image
776              * Let p be the currect pixel;
777              * tmp <- max{J(p) union J(p) neighbors in raster order}
778              * if (tmp > I(p))
779              *   J(p) <- tmp
780              * end */
781         for (i = 0; i < h; i++) {
782             lines = datas + i * wpls;
783             linem = datam + i * wplm;
784             for (j = 0; j < w; j++) {
785                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
786                     maxval = GET_DATA_BYTE(lines, j);
787                     if (i > 0) {
788                         val2 = GET_DATA_BYTE(lines - wpls, j);
789                         maxval = L_MAX(maxval, val2);
790                     }
791                     if (j > 0) {
792                         val4 = GET_DATA_BYTE(lines, j - 1);
793                         maxval = L_MAX(maxval, val4);
794                     }
795                     if (maxval > maskval)
796                         SET_DATA_BYTE(lines, j, maxval);
797                 }
798             }
799         }
800 
801             /* LR --> UL scan (anti-raster order)
802              * If I : mask image
803              *    J : marker image
804              * Let p be the currect pixel;
805              * tmp <- max{J(p) union J(p) neighbors in anti-raster order}
806              * if (tmp > I(p))
807              *   J(p) <- tmp
808              * end */
809         for (i = imax; i >= 0; i--) {
810             lines = datas + i * wpls;
811             linem = datam + i * wplm;
812             for (j = jmax; j >= 0; j--) {
813                 boolval = FALSE;
814                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
815                     val = maxval = GET_DATA_BYTE(lines, j);
816                     if (i < imax) {
817                         val7 = GET_DATA_BYTE(lines + wpls, j);
818                         maxval = L_MAX(maxval, val7);
819                     }
820                     if (j < jmax) {
821                         val5 = GET_DATA_BYTE(lines, j + 1);
822                         maxval = L_MAX(maxval, val5);
823                     }
824                     if (maxval > maskval)
825                         SET_DATA_BYTE(lines, j, maxval);
826                     val = GET_DATA_BYTE(lines, j);
827 
828                         /*
829                          * If there exists a point (q) which belongs to J(p)
830                          * neighbors in anti-raster order such that J(q) < J(p)
831                          * and J(p) > I(q) then
832                          * fifo_add(p) */
833                     if (i < imax) {
834                         val7 = GET_DATA_BYTE(lines + wpls, j);
835                         if ((val7 < val) &&
836                             (val > GET_DATA_BYTE(linem + wplm, j))) {
837                             boolval = TRUE;
838                         }
839                     }
840                     if (j < jmax) {
841                         val5 = GET_DATA_BYTE(lines, j + 1);
842                         if (!boolval && (val5 < val) &&
843                             (val > GET_DATA_BYTE(linem, j + 1))) {
844                             boolval = TRUE;
845                         }
846                     }
847                     if (boolval) {
848                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
849                         pixel->x = i;
850                         pixel->y = j;
851                         lqueueAdd(lq_pixel, pixel);
852                     }
853                 }
854             }
855         }
856 
857             /* Propagation step:
858              *        while fifo_empty = false
859              *          p <- fifo_first()
860              *          for every pixel (q) belong to neighbors of (p)
861              *            if J(q) < J(p) and J(p) > I(q)
862              *              J(q) <- min(J(p), I(q));
863              *              fifo_add(q);
864              *            end
865              *          end
866              *        end */
867         queue_size = lqueueGetCount(lq_pixel);
868         while (queue_size) {
869             pixel = (L_PIXEL *)lqueueRemove(lq_pixel);
870             i = pixel->x;
871             j = pixel->y;
872             FREE(pixel);
873             lines = datas + i * wpls;
874             linem = datam + i * wplm;
875 
876             if ((val = GET_DATA_BYTE(lines, j)) > 0) {
877                 if (i > 0) {
878                     val2 = GET_DATA_BYTE(lines - wpls, j);
879                     maskval = GET_DATA_BYTE(linem - wplm, j);
880                     if (val > val2 && val > maskval) {
881                         SET_DATA_BYTE(lines - wpls, j, val);
882                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
883                         pixel->x = i - 1;
884                         pixel->y = j;
885                         lqueueAdd(lq_pixel, pixel);
886                     }
887 
888                 }
889                 if (j > 0) {
890                     val4 = GET_DATA_BYTE(lines, j - 1);
891                     maskval = GET_DATA_BYTE(linem, j - 1);
892                     if (val > val4 && val > maskval) {
893                         SET_DATA_BYTE(lines, j - 1, val);
894                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
895                         pixel->x = i;
896                         pixel->y = j - 1;
897                         lqueueAdd(lq_pixel, pixel);
898                     }
899                 }
900                 if (i < imax) {
901                     val7 = GET_DATA_BYTE(lines + wpls, j);
902                     maskval = GET_DATA_BYTE(linem + wplm, j);
903                     if (val > val7 && val > maskval) {
904                         SET_DATA_BYTE(lines + wpls, j, val);
905                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
906                         pixel->x = i + 1;
907                         pixel->y = j;
908                         lqueueAdd(lq_pixel, pixel);
909                     }
910                 }
911                 if (j < jmax) {
912                     val5 = GET_DATA_BYTE(lines, j + 1);
913                     maskval = GET_DATA_BYTE(linem, j + 1);
914                     if (val > val5 && val > maskval) {
915                         SET_DATA_BYTE(lines, j + 1, val);
916                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
917                         pixel->x = i;
918                         pixel->y = j + 1;
919                         lqueueAdd(lq_pixel, pixel);
920                     }
921                 }
922             }
923 
924             queue_size = lqueueGetCount(lq_pixel);
925         }
926 
927         break;
928 
929     case 8:
930             /* UL --> LR scan  (Raster Order)
931              * If I : mask image
932              *    J : marker image
933              * Let p be the currect pixel;
934              * tmp <- max{J(p) union J(p) neighbors in raster order}
935              * if (tmp > I(p))
936              *   J(p) <- tmp
937              * end */
938         for (i = 0; i < h; i++) {
939             lines = datas + i * wpls;
940             linem = datam + i * wplm;
941             for (j = 0; j < w; j++) {
942                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
943                     maxval = GET_DATA_BYTE(lines, j);
944                     if (i > 0) {
945                         if (j > 0) {
946                             val1 = GET_DATA_BYTE(lines - wpls, j - 1);
947                             maxval = L_MAX(maxval, val1);
948                         }
949                         if (j < jmax) {
950                             val3 = GET_DATA_BYTE(lines - wpls, j + 1);
951                             maxval = L_MAX(maxval, val3);
952                         }
953                         val2 = GET_DATA_BYTE(lines - wpls, j);
954                         maxval = L_MAX(maxval, val2);
955                     }
956                     if (j > 0) {
957                         val4 = GET_DATA_BYTE(lines, j - 1);
958                         maxval = L_MAX(maxval, val4);
959                     }
960                     if (maxval > maskval)
961                         SET_DATA_BYTE(lines, j, maxval);
962                 }
963             }
964         }
965 
966             /* LR --> UL scan (anti-raster order)
967              * If I : mask image
968              *    J : marker image
969              * Let p be the currect pixel;
970              * tmp <- max{J(p) union J(p) neighbors in anti-raster order}
971              * if (tmp > I(p))
972              *   J(p) <- tmp
973              * end */
974         for (i = imax; i >= 0; i--) {
975             lines = datas + i * wpls;
976             linem = datam + i * wplm;
977             for (j = jmax; j >= 0; j--) {
978                 boolval = FALSE;
979                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
980                     maxval = GET_DATA_BYTE(lines, j);
981                     if (i < imax) {
982                         if (j > 0) {
983                             val6 = GET_DATA_BYTE(lines + wpls, j - 1);
984                             maxval = L_MAX(maxval, val6);
985                         }
986                         if (j < jmax) {
987                             val8 = GET_DATA_BYTE(lines + wpls, j + 1);
988                             maxval = L_MAX(maxval, val8);
989                         }
990                         val7 = GET_DATA_BYTE(lines + wpls, j);
991                         maxval = L_MAX(maxval, val7);
992                     }
993                     if (j < jmax) {
994                         val5 = GET_DATA_BYTE(lines, j + 1);
995                         maxval = L_MAX(maxval, val5);
996                     }
997                     if (maxval > maskval)
998                         SET_DATA_BYTE(lines, j, maxval);
999                     val = GET_DATA_BYTE(lines, j);
1000 
1001                         /*
1002                          * If there exists a point (q) which belongs to J(p)
1003                          * neighbors in anti-raster order such that J(q) < J(p)
1004                          * and J(p) > I(q) then
1005                          * fifo_add(p) */
1006                     if (i < imax) {
1007                         if (j > 0) {
1008                             val6 = GET_DATA_BYTE(lines + wpls, j - 1);
1009                             if ((val6 < val) &&
1010                                 (val > GET_DATA_BYTE(linem + wplm, j - 1))) {
1011                                 boolval = TRUE;
1012                             }
1013                         }
1014                         if (j < jmax) {
1015                             val8 = GET_DATA_BYTE(lines + wpls, j + 1);
1016                             if (!boolval && (val8 < val) &&
1017                                 (val > GET_DATA_BYTE(linem + wplm, j + 1))) {
1018                                 boolval = TRUE;
1019                             }
1020                         }
1021                         val7 = GET_DATA_BYTE(lines + wpls, j);
1022                         if (!boolval && (val7 < val) &&
1023                             (val > GET_DATA_BYTE(linem + wplm, j))) {
1024                             boolval = TRUE;
1025                         }
1026                     }
1027                     if (j < jmax) {
1028                         val5 = GET_DATA_BYTE(lines, j + 1);
1029                         if (!boolval && (val5 < val) &&
1030                             (val > GET_DATA_BYTE(linem, j + 1))) {
1031                             boolval = TRUE;
1032                         }
1033                     }
1034                     if (boolval) {
1035                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1036                         pixel->x = i;
1037                         pixel->y = j;
1038                         lqueueAdd(lq_pixel, pixel);
1039                     }
1040                 }
1041             }
1042         }
1043 
1044             /* Propagation step:
1045              *        while fifo_empty = false
1046              *          p <- fifo_first()
1047              *          for every pixel (q) belong to neighbors of (p)
1048              *            if J(q) < J(p) and J(p) > I(q)
1049              *              J(q) <- min(J(p), I(q));
1050              *              fifo_add(q);
1051              *            end
1052              *          end
1053              *        end */
1054         queue_size = lqueueGetCount(lq_pixel);
1055         while (queue_size) {
1056             pixel = (L_PIXEL *)lqueueRemove(lq_pixel);
1057             i = pixel->x;
1058             j = pixel->y;
1059             FREE(pixel);
1060             lines = datas + i * wpls;
1061             linem = datam + i * wplm;
1062 
1063             if ((val = GET_DATA_BYTE(lines, j)) > 0) {
1064                 if (i > 0) {
1065                     if (j > 0) {
1066                         val1 = GET_DATA_BYTE(lines - wpls, j - 1);
1067                         maskval = GET_DATA_BYTE(linem - wplm, j - 1);
1068                         if (val > val1 && val > maskval) {
1069                             SET_DATA_BYTE(lines - wpls, j - 1, val);
1070                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1071                             pixel->x = i - 1;
1072                             pixel->y = j - 1;
1073                             lqueueAdd(lq_pixel, pixel);
1074                         }
1075                     }
1076                     if (j < jmax) {
1077                         val3 = GET_DATA_BYTE(lines - wpls, j + 1);
1078                         maskval = GET_DATA_BYTE(linem - wplm, j + 1);
1079                         if (val > val3 && val > maskval) {
1080                             SET_DATA_BYTE(lines - wpls, j + 1, val);
1081                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1082                             pixel->x = i - 1;
1083                             pixel->y = j + 1;
1084                             lqueueAdd(lq_pixel, pixel);
1085                         }
1086                     }
1087                     val2 = GET_DATA_BYTE(lines - wpls, j);
1088                     maskval = GET_DATA_BYTE(linem - wplm, j);
1089                     if (val > val2 && val > maskval) {
1090                         SET_DATA_BYTE(lines - wpls, j, val);
1091                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1092                         pixel->x = i - 1;
1093                         pixel->y = j;
1094                         lqueueAdd(lq_pixel, pixel);
1095                     }
1096 
1097                 }
1098                 if (j > 0) {
1099                     val4 = GET_DATA_BYTE(lines, j - 1);
1100                     maskval = GET_DATA_BYTE(linem, j - 1);
1101                     if (val > val4 && val > maskval) {
1102                         SET_DATA_BYTE(lines, j - 1, val);
1103                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1104                         pixel->x = i;
1105                         pixel->y = j - 1;
1106                         lqueueAdd(lq_pixel, pixel);
1107                     }
1108                 }
1109                 if (i < imax) {
1110                     if (j > 0) {
1111                         val6 = GET_DATA_BYTE(lines + wpls, j - 1);
1112                         maskval = GET_DATA_BYTE(linem + wplm, j - 1);
1113                         if (val > val6 && val > maskval) {
1114                             SET_DATA_BYTE(lines + wpls, j - 1, val);
1115                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1116                             pixel->x = i + 1;
1117                             pixel->y = j - 1;
1118                             lqueueAdd(lq_pixel, pixel);
1119                         }
1120                     }
1121                     if (j < jmax) {
1122                         val8 = GET_DATA_BYTE(lines + wpls, j + 1);
1123                         maskval = GET_DATA_BYTE(linem + wplm, j + 1);
1124                         if (val > val8 && val > maskval) {
1125                             SET_DATA_BYTE(lines + wpls, j + 1, val);
1126                             pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1127                             pixel->x = i + 1;
1128                             pixel->y = j + 1;
1129                             lqueueAdd(lq_pixel, pixel);
1130                         }
1131                     }
1132                     val7 = GET_DATA_BYTE(lines + wpls, j);
1133                     maskval = GET_DATA_BYTE(linem + wplm, j);
1134                     if (val > val7 && val > maskval) {
1135                         SET_DATA_BYTE(lines + wpls, j, val);
1136                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1137                         pixel->x = i + 1;
1138                         pixel->y = j;
1139                         lqueueAdd(lq_pixel, pixel);
1140                     }
1141                 }
1142                 if (j < jmax) {
1143                     val5 = GET_DATA_BYTE(lines, j + 1);
1144                     maskval = GET_DATA_BYTE(linem, j + 1);
1145                     if (val > val5 && val > maskval) {
1146                         SET_DATA_BYTE(lines, j + 1, val);
1147                         pixel = (L_PIXEL *)CALLOC(1, sizeof(L_PIXEL));
1148                         pixel->x = i;
1149                         pixel->y = j + 1;
1150                         lqueueAdd(lq_pixel, pixel);
1151                     }
1152                 }
1153             }
1154 
1155             queue_size = lqueueGetCount(lq_pixel);
1156         }
1157         break;
1158 
1159     default:
1160         lqueueDestroy(&lq_pixel, TRUE);
1161         L_ERROR("connectivity must be 4 or 8", procName);
1162     }
1163 
1164     lqueueDestroy(&lq_pixel, TRUE);
1165     return;
1166 }
1167 
1168 /*-----------------------------------------------------------------------*
1169  *                 Vincent's Iterative Grayscale Seedfill                *
1170  *-----------------------------------------------------------------------*/
1171 /*!
1172  *  seedfillGrayLowSimple()
1173  *
1174  *  Notes:
1175  *      (1) The pixels are numbered as follows:
1176  *              1  2  3
1177  *              4  x  5
1178  *              6  7  8
1179  *          This low-level filling operation consists of two scans,
1180  *          raster and anti-raster, covering the entire seed image.
1181  *          The caller typically iterates until the filling is
1182  *          complete.
1183  *      (2) The filling action can be visualized from the following example.
1184  *          Suppose the mask, which clips the fill, is a sombrero-shaped
1185  *          surface, where the highest point is 200 and the low pixels
1186  *          around the rim are 30.  Beyond the rim, the mask goes up a bit.
1187  *          Suppose the seed, which is filled, consists of a single point
1188  *          of height 150, located below the max of the mask, with
1189  *          the rest 0.  Then in the raster scan, nothing happens until
1190  *          the high seed point is encountered, and then this value is
1191  *          propagated right and down, until it hits the side of the
1192  *          sombrero.   The seed can never exceed the mask, so it fills
1193  *          to the rim, going lower along the mask surface.  When it
1194  *          passes the rim, the seed continues to fill at the rim
1195  *          height to the edge of the seed image.  Then on the
1196  *          anti-raster scan, the seed fills flat inside the
1197  *          sombrero to the upper and left, and then out from the
1198  *          rim as before.  The final result has a seed that is
1199  *          flat outside the rim, and inside it fills the sombrero
1200  *          but only up to 150.  If the rim height varies, the
1201  *          filled seed outside the rim will be at the highest
1202  *          point on the rim, which is a saddle point on the rim.
1203  */
1204 void
seedfillGrayLowSimple(l_uint32 * datas,l_int32 w,l_int32 h,l_int32 wpls,l_uint32 * datam,l_int32 wplm,l_int32 connectivity)1205 seedfillGrayLowSimple(l_uint32  *datas,
1206                       l_int32    w,
1207                       l_int32    h,
1208                       l_int32    wpls,
1209                       l_uint32  *datam,
1210                       l_int32    wplm,
1211                       l_int32    connectivity)
1212 {
1213 l_uint8    val2, val3, val4, val5, val7, val8;
1214 l_uint8    val, maxval, maskval;
1215 l_int32    i, j, imax, jmax;
1216 l_uint32  *lines, *linem;
1217 
1218     PROCNAME("seedfillGrayLowSimple");
1219 
1220     imax = h - 1;
1221     jmax = w - 1;
1222 
1223     switch (connectivity)
1224     {
1225     case 4:
1226             /* UL --> LR scan */
1227         for (i = 0; i < h; i++) {
1228             lines = datas + i * wpls;
1229             linem = datam + i * wplm;
1230             for (j = 0; j < w; j++) {
1231                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
1232                     maxval = 0;
1233                     if (i > 0)
1234                         maxval = GET_DATA_BYTE(lines - wpls, j);
1235                     if (j > 0) {
1236                         val4 = GET_DATA_BYTE(lines, j - 1);
1237                         maxval = L_MAX(maxval, val4);
1238                     }
1239                     val = GET_DATA_BYTE(lines, j);
1240                     maxval = L_MAX(maxval, val);
1241                     val = L_MIN(maxval, maskval);
1242                     SET_DATA_BYTE(lines, j, val);
1243                 }
1244             }
1245         }
1246 
1247             /* LR --> UL scan */
1248         for (i = imax; i >= 0; i--) {
1249             lines = datas + i * wpls;
1250             linem = datam + i * wplm;
1251             for (j = jmax; j >= 0; j--) {
1252                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
1253                     maxval = 0;
1254                     if (i < imax)
1255                         maxval = GET_DATA_BYTE(lines + wpls, j);
1256                     if (j < jmax) {
1257                         val5 = GET_DATA_BYTE(lines, j + 1);
1258                         maxval = L_MAX(maxval, val5);
1259                     }
1260                     val = GET_DATA_BYTE(lines, j);
1261                     maxval = L_MAX(maxval, val);
1262                     val = L_MIN(maxval, maskval);
1263                     SET_DATA_BYTE(lines, j, val);
1264                 }
1265             }
1266         }
1267         break;
1268 
1269     case 8:
1270             /* UL --> LR scan */
1271         for (i = 0; i < h; i++) {
1272             lines = datas + i * wpls;
1273             linem = datam + i * wplm;
1274             for (j = 0; j < w; j++) {
1275                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
1276                     maxval = 0;
1277                     if (i > 0) {
1278                         if (j > 0)
1279                             maxval = GET_DATA_BYTE(lines - wpls, j - 1);
1280                         if (j < jmax) {
1281                             val2 = GET_DATA_BYTE(lines - wpls, j + 1);
1282                             maxval = L_MAX(maxval, val2);
1283                         }
1284                         val3 = GET_DATA_BYTE(lines - wpls, j);
1285                         maxval = L_MAX(maxval, val3);
1286                     }
1287                     if (j > 0) {
1288                         val4 = GET_DATA_BYTE(lines, j - 1);
1289                         maxval = L_MAX(maxval, val4);
1290                     }
1291                     val = GET_DATA_BYTE(lines, j);
1292                     maxval = L_MAX(maxval, val);
1293                     val = L_MIN(maxval, maskval);
1294                     SET_DATA_BYTE(lines, j, val);
1295                 }
1296             }
1297         }
1298 
1299             /* LR --> UL scan */
1300         for (i = imax; i >= 0; i--) {
1301             lines = datas + i * wpls;
1302             linem = datam + i * wplm;
1303             for (j = jmax; j >= 0; j--) {
1304                 if ((maskval = GET_DATA_BYTE(linem, j)) > 0) {
1305                     maxval = 0;
1306                     if (i < imax) {
1307                         if (j > 0)
1308                             maxval = GET_DATA_BYTE(lines + wpls, j - 1);
1309                         if (j < jmax) {
1310                             val8 = GET_DATA_BYTE(lines + wpls, j + 1);
1311                             maxval = L_MAX(maxval, val8);
1312                         }
1313                         val7 = GET_DATA_BYTE(lines + wpls, j);
1314                         maxval = L_MAX(maxval, val7);
1315                     }
1316                     if (j < jmax) {
1317                         val5 = GET_DATA_BYTE(lines, j + 1);
1318                         maxval = L_MAX(maxval, val5);
1319                     }
1320                     val = GET_DATA_BYTE(lines, j);
1321                     maxval = L_MAX(maxval, val);
1322                     val = L_MIN(maxval, maskval);
1323                     SET_DATA_BYTE(lines, j, val);
1324                 }
1325             }
1326         }
1327         break;
1328 
1329     default:
1330         L_ERROR("connectivity must be 4 or 8", procName);
1331     }
1332 
1333     return;
1334 }
1335 
1336 
1337 /*!
1338  *  seedfillGrayInvLowSimple()
1339  *
1340  *  Notes:
1341  *      (1) The pixels are numbered as follows:
1342  *              1  2  3
1343  *              4  x  5
1344  *              6  7  8
1345  *          This low-level filling operation consists of two scans,
1346  *          raster and anti-raster, covering the entire seed image.
1347  *          The caller typically iterates until the filling is
1348  *          complete.
1349  *      (2) The "Inv" signifies the fact that in this case, filling
1350  *          of the seed only takes place when the seed value is
1351  *          greater than the mask value.  The mask will act to stop
1352  *          the fill when it is higher than the seed level.  (This is
1353  *          in contrast to conventional grayscale filling where the
1354  *          seed always fills below the mask.)
1355  *      (3) An example of use is a basin, described by the mask (pixm),
1356  *          where within the basin, the seed pix (pixs) gets filled to the
1357  *          height of the highest seed pixel that is above its
1358  *          corresponding max pixel.  Filling occurs while the
1359  *          propagating seed pixels in pixs are larger than the
1360  *          corresponding mask values in pixm.
1361  */
1362 void
seedfillGrayInvLowSimple(l_uint32 * datas,l_int32 w,l_int32 h,l_int32 wpls,l_uint32 * datam,l_int32 wplm,l_int32 connectivity)1363 seedfillGrayInvLowSimple(l_uint32  *datas,
1364                          l_int32    w,
1365                          l_int32    h,
1366                          l_int32    wpls,
1367                          l_uint32  *datam,
1368                          l_int32    wplm,
1369                          l_int32    connectivity)
1370 {
1371 l_uint8    val1, val2, val3, val4, val5, val6, val7, val8;
1372 l_uint8    maxval, maskval;
1373 l_int32    i, j, imax, jmax;
1374 l_uint32  *lines, *linem;
1375 
1376     PROCNAME("seedfillGrayInvLowSimple");
1377 
1378     imax = h - 1;
1379     jmax = w - 1;
1380 
1381     switch (connectivity)
1382     {
1383     case 4:
1384             /* UL --> LR scan */
1385         for (i = 0; i < h; i++) {
1386             lines = datas + i * wpls;
1387             linem = datam + i * wplm;
1388             for (j = 0; j < w; j++) {
1389                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
1390                     maxval = GET_DATA_BYTE(lines, j);
1391                     if (i > 0) {
1392                         val2 = GET_DATA_BYTE(lines - wpls, j);
1393                         maxval = L_MAX(maxval, val2);
1394                     }
1395                     if (j > 0) {
1396                         val4 = GET_DATA_BYTE(lines, j - 1);
1397                         maxval = L_MAX(maxval, val4);
1398                     }
1399                     if (maxval > maskval)
1400                         SET_DATA_BYTE(lines, j, maxval);
1401                 }
1402             }
1403         }
1404 
1405             /* LR --> UL scan */
1406         for (i = imax; i >= 0; i--) {
1407             lines = datas + i * wpls;
1408             linem = datam + i * wplm;
1409             for (j = jmax; j >= 0; j--) {
1410                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
1411                     maxval = GET_DATA_BYTE(lines, j);
1412                     if (i < imax) {
1413                         val7 = GET_DATA_BYTE(lines + wpls, j);
1414                         maxval = L_MAX(maxval, val7);
1415                     }
1416                     if (j < jmax) {
1417                         val5 = GET_DATA_BYTE(lines, j + 1);
1418                         maxval = L_MAX(maxval, val5);
1419                     }
1420                     if (maxval > maskval)
1421                         SET_DATA_BYTE(lines, j, maxval);
1422                 }
1423             }
1424         }
1425         break;
1426 
1427     case 8:
1428             /* UL --> LR scan */
1429         for (i = 0; i < h; i++) {
1430             lines = datas + i * wpls;
1431             linem = datam + i * wplm;
1432             for (j = 0; j < w; j++) {
1433                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
1434                     maxval = GET_DATA_BYTE(lines, j);
1435                     if (i > 0) {
1436                         if (j > 0) {
1437                             val1 = GET_DATA_BYTE(lines - wpls, j - 1);
1438                             maxval = L_MAX(maxval, val1);
1439                         }
1440                         if (j < jmax) {
1441                             val2 = GET_DATA_BYTE(lines - wpls, j + 1);
1442                             maxval = L_MAX(maxval, val2);
1443                         }
1444                         val3 = GET_DATA_BYTE(lines - wpls, j);
1445                         maxval = L_MAX(maxval, val3);
1446                     }
1447                     if (j > 0) {
1448                         val4 = GET_DATA_BYTE(lines, j - 1);
1449                         maxval = L_MAX(maxval, val4);
1450                     }
1451                     if (maxval > maskval)
1452                         SET_DATA_BYTE(lines, j, maxval);
1453                 }
1454             }
1455         }
1456 
1457             /* LR --> UL scan */
1458         for (i = imax; i >= 0; i--) {
1459             lines = datas + i * wpls;
1460             linem = datam + i * wplm;
1461             for (j = jmax; j >= 0; j--) {
1462                 if ((maskval = GET_DATA_BYTE(linem, j)) < 255) {
1463                     maxval = GET_DATA_BYTE(lines, j);
1464                     if (i < imax) {
1465                         if (j > 0) {
1466                             val6 = GET_DATA_BYTE(lines + wpls, j - 1);
1467                             maxval = L_MAX(maxval, val6);
1468                         }
1469                         if (j < jmax) {
1470                             val8 = GET_DATA_BYTE(lines + wpls, j + 1);
1471                             maxval = L_MAX(maxval, val8);
1472                         }
1473                         val7 = GET_DATA_BYTE(lines + wpls, j);
1474                         maxval = L_MAX(maxval, val7);
1475                     }
1476                     if (j < jmax) {
1477                         val5 = GET_DATA_BYTE(lines, j + 1);
1478                         maxval = L_MAX(maxval, val5);
1479                     }
1480                     if (maxval > maskval)
1481                         SET_DATA_BYTE(lines, j, maxval);
1482                 }
1483             }
1484         }
1485         break;
1486 
1487     default:
1488         L_ERROR("connectivity must be 4 or 8", procName);
1489     }
1490 
1491     return;
1492 }
1493 
1494 
1495 /*-----------------------------------------------------------------------*
1496  *                   Vincent's Distance Function method                  *
1497  *-----------------------------------------------------------------------*/
1498 /*!
1499  *  distanceFunctionLow()
1500  */
1501 void
distanceFunctionLow(l_uint32 * datad,l_int32 w,l_int32 h,l_int32 d,l_int32 wpld,l_int32 connectivity)1502 distanceFunctionLow(l_uint32  *datad,
1503                     l_int32    w,
1504                     l_int32    h,
1505                     l_int32    d,
1506                     l_int32    wpld,
1507                     l_int32    connectivity)
1508 {
1509 l_int32    val1, val2, val3, val4, val5, val6, val7, val8, minval, val;
1510 l_int32    i, j, imax, jmax;
1511 l_uint32  *lined;
1512 
1513     PROCNAME("distanceFunctionLow");
1514 
1515         /* One raster scan followed by one anti-raster scan.
1516          * This does not re-set the 1-boundary of pixels that
1517          * were initialized to either 0 or maxval. */
1518     imax = h - 1;
1519     jmax = w - 1;
1520     switch (connectivity)
1521     {
1522     case 4:
1523         if (d == 8) {
1524                 /* UL --> LR scan */
1525             for (i = 1; i < imax; i++) {
1526                 lined = datad + i * wpld;
1527                 for (j = 1; j < jmax; j++) {
1528                     if ((val = GET_DATA_BYTE(lined, j)) > 0) {
1529                         val2 = GET_DATA_BYTE(lined - wpld, j);
1530                         val4 = GET_DATA_BYTE(lined, j - 1);
1531                         minval = L_MIN(val2, val4);
1532                         minval = L_MIN(minval, 254);
1533                         SET_DATA_BYTE(lined, j, minval + 1);
1534                     }
1535                 }
1536             }
1537 
1538                 /* LR --> UL scan */
1539             for (i = imax - 1; i > 0; i--) {
1540                 lined = datad + i * wpld;
1541                 for (j = jmax - 1; j > 0; j--) {
1542                     if ((val = GET_DATA_BYTE(lined, j)) > 0) {
1543                         val7 = GET_DATA_BYTE(lined + wpld, j);
1544                         val5 = GET_DATA_BYTE(lined, j + 1);
1545                         minval = L_MIN(val5, val7);
1546                         minval = L_MIN(minval + 1, val);
1547                         SET_DATA_BYTE(lined, j, minval);
1548                     }
1549                 }
1550             }
1551         }
1552         else {  /* d == 16 */
1553                 /* UL --> LR scan */
1554             for (i = 1; i < imax; i++) {
1555                 lined = datad + i * wpld;
1556                 for (j = 1; j < jmax; j++) {
1557                     if ((val = GET_DATA_TWO_BYTES(lined, j)) > 0) {
1558                         val2 = GET_DATA_TWO_BYTES(lined - wpld, j);
1559                         val4 = GET_DATA_TWO_BYTES(lined, j - 1);
1560                         minval = L_MIN(val2, val4);
1561                         minval = L_MIN(minval, 0xfffe);
1562                         SET_DATA_TWO_BYTES(lined, j, minval + 1);
1563                     }
1564                 }
1565             }
1566 
1567                 /* LR --> UL scan */
1568             for (i = imax - 1; i > 0; i--) {
1569                 lined = datad + i * wpld;
1570                 for (j = jmax - 1; j > 0; j--) {
1571                     if ((val = GET_DATA_TWO_BYTES(lined, j)) > 0) {
1572                         val7 = GET_DATA_TWO_BYTES(lined + wpld, j);
1573                         val5 = GET_DATA_TWO_BYTES(lined, j + 1);
1574                         minval = L_MIN(val5, val7);
1575                         minval = L_MIN(minval + 1, val);
1576                         SET_DATA_TWO_BYTES(lined, j, minval);
1577                     }
1578                 }
1579             }
1580         }
1581         break;
1582 
1583     case 8:
1584         if (d == 8) {
1585                 /* UL --> LR scan */
1586             for (i = 1; i < imax; i++) {
1587                 lined = datad + i * wpld;
1588                 for (j = 1; j < jmax; j++) {
1589                     if ((val = GET_DATA_BYTE(lined, j)) > 0) {
1590                         val1 = GET_DATA_BYTE(lined - wpld, j - 1);
1591                         val2 = GET_DATA_BYTE(lined - wpld, j);
1592                         val3 = GET_DATA_BYTE(lined - wpld, j + 1);
1593                         val4 = GET_DATA_BYTE(lined, j - 1);
1594                         minval = L_MIN(val1, val2);
1595                         minval = L_MIN(minval, val3);
1596                         minval = L_MIN(minval, val4);
1597                         minval = L_MIN(minval, 254);
1598                         SET_DATA_BYTE(lined, j, minval + 1);
1599                     }
1600                 }
1601             }
1602 
1603                 /* LR --> UL scan */
1604             for (i = imax - 1; i > 0; i--) {
1605                 lined = datad + i * wpld;
1606                 for (j = jmax - 1; j > 0; j--) {
1607                     if ((val = GET_DATA_BYTE(lined, j)) > 0) {
1608                         val8 = GET_DATA_BYTE(lined + wpld, j + 1);
1609                         val7 = GET_DATA_BYTE(lined + wpld, j);
1610                         val6 = GET_DATA_BYTE(lined + wpld, j - 1);
1611                         val5 = GET_DATA_BYTE(lined, j + 1);
1612                         minval = L_MIN(val8, val7);
1613                         minval = L_MIN(minval, val6);
1614                         minval = L_MIN(minval, val5);
1615                         minval = L_MIN(minval + 1, val);
1616                         SET_DATA_BYTE(lined, j, minval);
1617                     }
1618                 }
1619             }
1620         }
1621         else {  /* d == 16 */
1622                 /* UL --> LR scan */
1623             for (i = 1; i < imax; i++) {
1624                 lined = datad + i * wpld;
1625                 for (j = 1; j < jmax; j++) {
1626                     if ((val = GET_DATA_TWO_BYTES(lined, j)) > 0) {
1627                         val1 = GET_DATA_TWO_BYTES(lined - wpld, j - 1);
1628                         val2 = GET_DATA_TWO_BYTES(lined - wpld, j);
1629                         val3 = GET_DATA_TWO_BYTES(lined - wpld, j + 1);
1630                         val4 = GET_DATA_TWO_BYTES(lined, j - 1);
1631                         minval = L_MIN(val1, val2);
1632                         minval = L_MIN(minval, val3);
1633                         minval = L_MIN(minval, val4);
1634                         minval = L_MIN(minval, 0xfffe);
1635                         SET_DATA_TWO_BYTES(lined, j, minval + 1);
1636                     }
1637                 }
1638             }
1639 
1640                 /* LR --> UL scan */
1641             for (i = imax - 1; i > 0; i--) {
1642                 lined = datad + i * wpld;
1643                 for (j = jmax - 1; j > 0; j--) {
1644                     if ((val = GET_DATA_TWO_BYTES(lined, j)) > 0) {
1645                         val8 = GET_DATA_TWO_BYTES(lined + wpld, j + 1);
1646                         val7 = GET_DATA_TWO_BYTES(lined + wpld, j);
1647                         val6 = GET_DATA_TWO_BYTES(lined + wpld, j - 1);
1648                         val5 = GET_DATA_TWO_BYTES(lined, j + 1);
1649                         minval = L_MIN(val8, val7);
1650                         minval = L_MIN(minval, val6);
1651                         minval = L_MIN(minval, val5);
1652                         minval = L_MIN(minval + 1, val);
1653                         SET_DATA_TWO_BYTES(lined, j, minval);
1654                     }
1655                 }
1656             }
1657         }
1658         break;
1659 
1660     default:
1661         L_ERROR("connectivity must be 4 or 8", procName);
1662         break;
1663     }
1664 
1665     return;
1666 }
1667 
1668 
1669 /*-----------------------------------------------------------------------*
1670  *                 Seed spread (based on distance function)              *
1671  *-----------------------------------------------------------------------*/
1672 /*!
1673  *  seedspreadLow()
1674  *
1675  *    See pixSeedspread() for a brief description of the algorithm here.
1676  */
1677 void
seedspreadLow(l_uint32 * datad,l_int32 w,l_int32 h,l_int32 wpld,l_uint32 * datat,l_int32 wplt,l_int32 connectivity)1678 seedspreadLow(l_uint32  *datad,
1679               l_int32    w,
1680               l_int32    h,
1681               l_int32    wpld,
1682               l_uint32  *datat,
1683               l_int32    wplt,
1684               l_int32    connectivity)
1685 {
1686 l_int32    val1t, val2t, val3t, val4t, val5t, val6t, val7t, val8t;
1687 l_int32    i, j, imax, jmax, minval, valt, vald;
1688 l_uint32  *linet, *lined;
1689 
1690     PROCNAME("seedspreadLow");
1691 
1692         /* One raster scan followed by one anti-raster scan.
1693          * pixt is initialized to have 0 on pixels where the
1694          * input is specified in pixd, and to have 1 on all
1695          * other pixels.  We only change pixels in pixt and pixd
1696          * that are non-zero in pixt. */
1697     imax = h - 1;
1698     jmax = w - 1;
1699     switch (connectivity)
1700     {
1701     case 4:
1702             /* UL --> LR scan */
1703         for (i = 1; i < h; i++) {
1704             linet = datat + i * wplt;
1705             lined = datad + i * wpld;
1706             for (j = 1; j < jmax; j++) {
1707                 if ((valt = GET_DATA_TWO_BYTES(linet, j)) > 0) {
1708                     val2t = GET_DATA_TWO_BYTES(linet - wplt, j);
1709                     val4t = GET_DATA_TWO_BYTES(linet, j - 1);
1710                     minval = L_MIN(val2t, val4t);
1711                     minval = L_MIN(minval, 0xfffe);
1712                     SET_DATA_TWO_BYTES(linet, j, minval + 1);
1713                     if (val2t < val4t)
1714                         vald = GET_DATA_BYTE(lined - wpld, j);
1715                     else
1716                         vald = GET_DATA_BYTE(lined, j - 1);
1717                     SET_DATA_BYTE(lined, j, vald);
1718                 }
1719             }
1720         }
1721 
1722             /* LR --> UL scan */
1723         for (i = imax - 1; i > 0; i--) {
1724             linet = datat + i * wplt;
1725             lined = datad + i * wpld;
1726             for (j = jmax - 1; j > 0; j--) {
1727                 if ((valt = GET_DATA_TWO_BYTES(linet, j)) > 0) {
1728                     val7t = GET_DATA_TWO_BYTES(linet + wplt, j);
1729                     val5t = GET_DATA_TWO_BYTES(linet, j + 1);
1730                     minval = L_MIN(val5t, val7t);
1731                     minval = L_MIN(minval + 1, valt);
1732                     if (valt > minval) {  /* replace */
1733                         SET_DATA_TWO_BYTES(linet, j, minval);
1734                         if (val5t < val7t)
1735                             vald = GET_DATA_BYTE(lined, j + 1);
1736                         else
1737                             vald = GET_DATA_BYTE(lined + wplt, j);
1738                         SET_DATA_BYTE(lined, j, vald);
1739                     }
1740                 }
1741             }
1742         }
1743         break;
1744     case 8:
1745             /* UL --> LR scan */
1746         for (i = 1; i < h; i++) {
1747             linet = datat + i * wplt;
1748             lined = datad + i * wpld;
1749             for (j = 1; j < jmax; j++) {
1750                 if ((valt = GET_DATA_TWO_BYTES(linet, j)) > 0) {
1751                     val1t = GET_DATA_TWO_BYTES(linet - wplt, j - 1);
1752                     val2t = GET_DATA_TWO_BYTES(linet - wplt, j);
1753                     val3t = GET_DATA_TWO_BYTES(linet - wplt, j + 1);
1754                     val4t = GET_DATA_TWO_BYTES(linet, j - 1);
1755                     minval = L_MIN(val1t, val2t);
1756                     minval = L_MIN(minval, val3t);
1757                     minval = L_MIN(minval, val4t);
1758                     minval = L_MIN(minval, 0xfffe);
1759                     SET_DATA_TWO_BYTES(linet, j, minval + 1);
1760                     if (minval == val1t)
1761                         vald = GET_DATA_BYTE(lined - wpld, j - 1);
1762                     else if (minval == val2t)
1763                         vald = GET_DATA_BYTE(lined - wpld, j);
1764                     else if (minval == val3t)
1765                         vald = GET_DATA_BYTE(lined - wpld, j + 1);
1766                     else  /* minval == val4t */
1767                         vald = GET_DATA_BYTE(lined, j - 1);
1768                     SET_DATA_BYTE(lined, j, vald);
1769                 }
1770             }
1771         }
1772 
1773             /* LR --> UL scan */
1774         for (i = imax - 1; i > 0; i--) {
1775             linet = datat + i * wplt;
1776             lined = datad + i * wpld;
1777             for (j = jmax - 1; j > 0; j--) {
1778                 if ((valt = GET_DATA_TWO_BYTES(linet, j)) > 0) {
1779                     val8t = GET_DATA_TWO_BYTES(linet + wplt, j + 1);
1780                     val7t = GET_DATA_TWO_BYTES(linet + wplt, j);
1781                     val6t = GET_DATA_TWO_BYTES(linet + wplt, j - 1);
1782                     val5t = GET_DATA_TWO_BYTES(linet, j + 1);
1783                     minval = L_MIN(val8t, val7t);
1784                     minval = L_MIN(minval, val6t);
1785                     minval = L_MIN(minval, val5t);
1786                     minval = L_MIN(minval + 1, valt);
1787                     if (valt > minval) {  /* replace */
1788                         SET_DATA_TWO_BYTES(linet, j, minval);
1789                         if (minval == val5t + 1)
1790                             vald = GET_DATA_BYTE(lined, j + 1);
1791                         else if (minval == val6t + 1)
1792                             vald = GET_DATA_BYTE(lined + wpld, j - 1);
1793                         else if (minval == val7t + 1)
1794                             vald = GET_DATA_BYTE(lined + wpld, j);
1795                         else  /* minval == val8t + 1 */
1796                             vald = GET_DATA_BYTE(lined + wpld, j + 1);
1797                         SET_DATA_BYTE(lined, j, vald);
1798                     }
1799                 }
1800             }
1801         }
1802         break;
1803     default:
1804         L_ERROR("connectivity must be 4 or 8", procName);
1805         break;
1806     }
1807 
1808     return;
1809 }
1810