Mercurial > hgrepos > Python2 > PyMuPDF

comparison mupdf-source/source/fitz/warp.c @ 2:b50eed0cc0ef upstream

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ADD: MuPDF v1.26.7: the MuPDF source as downloaded by a default build of PyMuPDF 1.26.4. The directory name has changed: no version number in the expanded directory now.
author Franz Glasner <fzglas.hg@dom66.de>
date Mon, 15 Sep 2025 11:43:07 +0200
parents
children
comparison
equal deleted inserted replaced
1:1d09e1dec1d9 2:b50eed0cc0ef
1 // Copyright (C) 2004-2025 Artifex Software, Inc.
2 //
3 // This file is part of MuPDF.
4 //
5 // MuPDF is free software: you can redistribute it and/or modify it under the
6 // terms of the GNU Affero General Public License as published by the Free
7 // Software Foundation, either version 3 of the License, or (at your option)
8 // any later version.
9 //
10 // MuPDF is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU Affero General Public License
16 // along with MuPDF. If not, see <https://www.gnu.org/licenses/agpl-3.0.en.html>
17 //
18 // Alternative licensing terms are available from the licensor.
19 // For commercial licensing, see <https://www.artifex.com/> or contact
20 // Artifex Software, Inc., 39 Mesa Street, Suite 108A, San Francisco,
21 // CA 94129, USA, for further information.
22
23 #include "mupdf/fitz.h"
24
25 #include "pixmap-imp.h"
26
27 #include <string.h>
28
29 /* Define TIMINGS to get timing information dumped to stdout. */
30 #undef TIMINGS
31
32 /* Define WARP_DEBUG to get debugging output (and PNGs saved). Note
33 * that this will affect timings! */
34 #undef WARP_DEBUG
35
36 /* Define WARP_SPEW_DEBUG to get even more debug output (and PNGs). */
37 #undef WARP_SPEW_DEBUG
38
39 /* One reference suggested doing histogram equalisation. */
40 #define DO_HISTEQ
41
42 /* Define DETECT_DOCUMENT_RGB, and edge detection on RGB documents will
43 * look for edges in just the R,G,B planes as well as the grey plane. */
44 #undef DETECT_DOCUMENT_RGB
45
46 #undef SLOW_INTERPOLATION
47 #undef SLOW_WARPING
48
49 #ifdef WARP_DEBUG
50 static void
51 debug_printf(fz_context *ctx, const char *fmt, ...)
52 {
53 char text[1024];
54 va_list list;
55 va_start(list, fmt);
56 vsnprintf(text, sizeof(text), fmt, list);
57 va_end(list);
58
59 #ifdef _WIN32
60 fz_write_string(ctx, fz_stdods(ctx), text);
61 #endif
62 fputs(text, stderr);
63 }
64 #else
65 #define debug_printf(CTX, FMT, ...) do {} while (0)
66 #endif
67
68 #if defined(TIMINGS) && defined(_WIN32)
69
70 #include "windows.h"
71
72 #define MAX_TIMERS 256
73 static struct {
74 int started;
75 int stackptr;
76 int stack[MAX_TIMERS];
77 const char *name[MAX_TIMERS];
78 LARGE_INTEGER time[MAX_TIMERS];
79 } timer;
80 #define START_TIME() \
81 do {\
82 int i = timer.stack[timer.stackptr++] = timer.started++;\
83 QueryPerformanceCounter(&timer.time[i]);\
84 } while (0)
85 #define END_TIME(NAME) \
86 do {\
87 LARGE_INTEGER end;\
88 int i;\
89 QueryPerformanceCounter(&end);\
90 i = timer.stack[--timer.stackptr];\
91 timer.time[i].QuadPart = end.QuadPart - timer.time[i].QuadPart;\
92 timer.name[i] = NAME;\
93 } while (0)
94 #define DUMP_TIMES() \
95 do {\
96 int i;\
97 LARGE_INTEGER freq;\
98 QueryPerformanceFrequency(&freq);\
99 float f = freq.QuadPart;\
100 for (i = 0; i < timer.started; i++)\
101 debug_printf(ctx, "%s: %g\n", timer.name[i], (int)1000*timer.time[i].QuadPart/f);\
102 } while (0)
103 #else
104 #define START_TIME() do {} while(0)
105 #define END_TIME(NAME) do {} while(0)
106 #define DUMP_TIMES() do {} while(0)
107 #endif
108
109 typedef struct
110 {
111 int x;
112 int y;
113 } fz_ipoint;
114
115 typedef struct
116 {
117 int i;
118 int f;
119 int di;
120 int df;
121 } fz_bresenham_core;
122
123 typedef struct
124 {
125 fz_bresenham_core c;
126 int n;
127 } fz_bresenham;
128
129 typedef struct
130 {
131 fz_bresenham_core x;
132 fz_bresenham_core y;
133 int n;
134 } fz_ipoint_bresenham;
135
136 typedef struct
137 {
138 fz_bresenham_core sx;
139 fz_bresenham_core sy;
140 fz_bresenham_core ex;
141 fz_bresenham_core ey;
142 int n;
143 } fz_ipoint2_bresenham;
144
145 static inline fz_bresenham_core
146 init_bresenham_core(int start, int end, int n)
147 {
148 fz_bresenham_core b;
149 int delta = end-start;
150
151 b.di = n == 0 ? 0 : delta/n;
152 b.df = delta - n*b.di; /* 0 <= b.df < n */
153 if (b.df < 0)
154 b.di--, b.df += n;
155 /* Starts with bi.i = start, bi.f = n, and then does a half
156 * step. */
157 b.i = start + (b.di>>1);
158 b.f = n - (((b.di & 1) * n + b.df)>>1);
159
160 return b;
161 }
162
163 #ifdef CURRENTLY_UNUSED
164 static inline fz_bresenham
165 init_bresenham(int start, int end, int n)
166 {
167 fz_bresenham b;
168
169 b.c = init_bresenham_core(start, end, n);
170 b.n = n;
171
172 return b;
173 }
174
175 static inline void
176 step(fz_bresenham *b)
177 {
178 step_core(&b->c, b->n);
179 }
180 #endif
181
182 static inline void
183 step_core(fz_bresenham_core *b, int n)
184 {
185 b->i += b->di;
186 b->f -= b->df;
187 if (b->f <= 0)
188 {
189 b->f += n;
190 b->i++;
191 }
192 }
193
194 static inline fz_ipoint_bresenham
195 init_ip_bresenham(fz_ipoint start, fz_ipoint end, int n)
196 {
197 fz_ipoint_bresenham b;
198
199 b.x = init_bresenham_core(start.x, end.x, n);
200 b.y = init_bresenham_core(start.y, end.y, n);
201 b.n = n;
202
203 return b;
204 }
205
206 static inline void
207 step_ip(fz_ipoint_bresenham *b)
208 {
209 step_core(&b->x, b->n);
210 step_core(&b->y, b->n);
211 }
212
213 static inline fz_ipoint
214 current_ip(const fz_ipoint_bresenham *b)
215 {
216 fz_ipoint ip;
217
218 ip.x = b->x.i;
219 ip.y = b->y.i;
220
221 return ip;
222 }
223
224 static inline fz_ipoint2_bresenham
225 init_ip2_bresenham(fz_ipoint ss, fz_ipoint se, fz_ipoint es, fz_ipoint ee, int n)
226 {
227 fz_ipoint2_bresenham b;
228
229 b.sx = init_bresenham_core(ss.x, se.x, n);
230 b.sy = init_bresenham_core(ss.y, se.y, n);
231 b.ex = init_bresenham_core(es.x, ee.x, n);
232 b.ey = init_bresenham_core(es.y, ee.y, n);
233 b.n = n;
234
235 return b;
236 }
237
238 static inline void
239 step_ip2(fz_ipoint2_bresenham *b)
240 {
241 step_core(&b->sx, b->n);
242 step_core(&b->sy, b->n);
243 step_core(&b->ex, b->n);
244 step_core(&b->ey, b->n);
245 }
246
247 static inline fz_ipoint
248 start_ip(const fz_ipoint2_bresenham *b)
249 {
250 fz_ipoint ip;
251
252 ip.x = b->sx.i;
253 ip.y = b->sy.i;
254
255 return ip;
256 }
257
258 static fz_forceinline fz_ipoint
259 end_ip(const fz_ipoint2_bresenham *b)
260 {
261 fz_ipoint ip;
262
263 ip.x = b->ex.i;
264 ip.y = b->ey.i;
265
266 return ip;
267 }
268
269 static void
270 interp_n(unsigned char *d, const unsigned char *s0,
271 const unsigned char *s1, int f, int n)
272 {
273 do
274 {
275 int a = *s0++;
276 int b = *s1++ - a;
277 *d++ = ((a<<8) + b*f + 128)>>8;
278 }
279 while (--n);
280 }
281
282 static void
283 interp2_n(unsigned char *d, const unsigned char *s0,
284 const unsigned char *s1, const unsigned char *s2,
285 int f0, int f1, int n)
286 {
287 do
288 {
289 int a = *s0++;
290 int b = *s1++ - a;
291 int c;
292 a = (a<<8) + b*f0;
293 c = (*s2++<<8) - a;
294 *d++ = ((a<<8) + c*f1 + (1<<15))>>16;
295 }
296 while (--n);
297 }
298
299 static inline void
300 copy_pixel(unsigned char *d, const fz_pixmap *src, fz_ipoint p)
301 {
302 int u = p.x>>8;
303 int v = p.y>>8;
304 int fu = p.x & 255;
305 int fv = p.y & 255;
306 int n = src->n;
307 const unsigned char *s;
308 ptrdiff_t stride = src->stride;
309
310 if (u < 0)
311 u = 0, fu = 0;
312 else if (u >= src->w-1)
313 u = src->w-1, fu = 0;
314
315 if (v < 0)
316 v = 0, fv = 0;
317 else if (v >= src->h-1)
318 v = src->h-1, fv = 0;
319
320 s = &src->samples[u * n + v * stride];
321
322 #ifdef SLOW_INTERPOLATION
323 {
324 int i;
325
326 for (i = 0; i < n; i++)
327 {
328 int v0 = s[0];
329 int v1 = s[n];
330 int v2 = s[stride];
331 int v3 = s[stride+n];
332 int v01, v23, v;
333 v01 = (v0<<8) + (v1-v0)*fu;
334 v23 = (v2<<8) + (v3-v2)*fu;
335 v = (v01<<8) + (v23-v01)*fv;
336 assert(v >= 0 && v < (1<<24)-32768);
337 *d++ = (v + 32768)>>16;
338 s++;
339 }
340 return;
341 }
342 #else
343 if (fu == 0)
344 {
345 if (fv == 0)
346 {
347 /* Copy single pixel */
348 memcpy(d, s, n);
349 return;
350 }
351 /* interpolate y pixels */
352 interp_n(d, s, s + stride, fv, n);
353 return;
354 }
355 if (fv == 0)
356 {
357 /* interpolate x pixels */
358 interp_n(d, s, s+n, fu, n);
359 return;
360 }
361
362 if (fu <= fv)
363 {
364 /* Top half of the trapezoid. */
365 interp2_n(d, s, s+stride, s+stride+n, fv, fu, n);
366 }
367 else
368 {
369 /* Bottom half of the trapezoid. */
370 interp2_n(d, s, s+n, s+stride+n, fu, fv, n);
371 }
372 #endif
373 }
374
375 static void
376 warp_core(unsigned char *d, int n, int width, int height, int stride,
377 const fz_ipoint corner[4], const fz_pixmap *src)
378 {
379 fz_ipoint2_bresenham row_bres;
380 int x;
381
382 /* We have a bresenham pair for how to move the start
383 * and end of the row each y step. */
384 row_bres = init_ip2_bresenham(corner[0], corner[3],
385 corner[1], corner[2], height);
386 stride -= width * n;
387
388 #ifdef SLOW_WARPING
389 {
390 int h;
391 for (h = 0 ; h < height ; h++)
392 {
393 int sx = corner[0].x + (corner[3].x - corner[0].x)*h/height;
394 int sy = corner[0].y + (corner[3].y - corner[0].y)*h/height;
395 int ex = corner[1].x + (corner[2].x - corner[1].x)*h/height;
396 int ey = corner[1].y + (corner[2].y - corner[1].y)*h/height;
397 for (x = 0; x < width; x++)
398 {
399 fz_ipoint p;
400 p.x = sx + (ex-sx)*x/width;
401 p.y = sy + (ey-sy)*x/width;
402 copy_pixel(d, src, p);
403 d += n;
404 }
405 d += stride;
406 }
407 }
408 #else
409 for (; height > 0; height--)
410 {
411 /* We have a bresenham for how to move the
412 * current pixel across the row. */
413 fz_ipoint_bresenham pix_bres;
414 pix_bres = init_ip_bresenham(start_ip(&row_bres),
415 end_ip(&row_bres),
416 width);
417 for (x = width; x > 0; x--)
418 {
419 /* Copy pixel */
420 copy_pixel(d, src, current_ip(&pix_bres));
421 d += n;
422 step_ip(&pix_bres);
423 }
424
425 /* step to the next line. */
426 step_ip2(&row_bres);
427 d += stride;
428 }
429 #endif
430 }
431
432 /*
433 points are clockwise from NW.
434
435 This performs simple affine warping.
436 */
437 fz_pixmap *
438 fz_warp_pixmap(fz_context *ctx, fz_pixmap *src, fz_quad points, int width, int height)
439 {
440 fz_pixmap *dst;
441
442 if (src == NULL)
443 return NULL;
444
445 if (width >= (1<<24) || width < 0 || height >= (1<<24) || height < 0)
446 fz_throw(ctx, FZ_ERROR_LIMIT, "Bad width/height");
447
448 dst = fz_new_pixmap(ctx, src->colorspace, width, height,
449 src->seps, src->alpha);
450 dst->xres = src->xres;
451 dst->yres = src->yres;
452
453 fz_try(ctx)
454 {
455 unsigned char *d = dst->samples;
456 int n = dst->n;
457 fz_ipoint corner[4];
458
459 /* Find the corner texture positions as fixed point */
460 corner[0].x = (int)(points.ul.x * 256 + 128);
461 corner[0].y = (int)(points.ul.y * 256 + 128);
462 corner[1].x = (int)(points.ur.x * 256 + 128);
463 corner[1].y = (int)(points.ur.y * 256 + 128);
464 corner[2].x = (int)(points.ll.x * 256 + 128);
465 corner[2].y = (int)(points.ll.y * 256 + 128);
466 corner[3].x = (int)(points.lr.x * 256 + 128);
467 corner[3].y = (int)(points.lr.y * 256 + 128);
468
469 warp_core(d, n, width, height, width * n, corner, src);
470 }
471 fz_catch(ctx)
472 {
473 fz_drop_pixmap(ctx, dst);
474 fz_rethrow(ctx);
475 }
476
477 return dst;
478 }
479
480 static float
481 dist(fz_point a, fz_point b)
482 {
483 float x = a.x-b.x;
484 float y = a.y-b.y;
485
486 return sqrtf(x*x+y*y);
487 }
488
489 /* Again, affine warping, but this time where the destination width/height
490 * are chosen automatically. */
491 fz_pixmap *
492 fz_autowarp_pixmap(fz_context *ctx, fz_pixmap *src, fz_quad points)
493 {
494 float w0 = dist(points.ur, points.ul);
495 float w1 = dist(points.ll, points.lr);
496 float h0 = dist(points.lr, points.ul);
497 float h1 = dist(points.ll, points.ur);
498 int w = (w0+w1+0.5)/2;
499 int h = (h0+h1+0.5)/2;
500
501 return fz_warp_pixmap(ctx, src, points, w, h);
502 }
503
504 /*
505 Corner detection: We shall steal the algorithm from the Dropbox
506 Document Scanner, as described here:
507
508 https://dropbox.tech/machine-learning/fast-and-accurate-document-detection-for-scanning
509
510 A good reference to the steps involved in the canny edge
511 detection process can be found here:
512
513 https://towardsdatascience.com/canny-edge-detection-step-by-step-in-python-computer-vision-b49c3a2d8123
514
515 This involves:
516 * Canny Edge Detection
517 * * Greyscale conversion
518 * * Noise reduction
519 * * Gradient Calculation
520 * * Non-maximum suppression
521 * * Double threshold
522 * * Edge tracking by Hysteresis
523 * Hough Transform to fix possible edges
524 * Computing intersections and scoring quads
525
526 We modify the gradient calculation with a simple scale to ensure we fill the range.
527 */
528
529 #ifdef DO_HISTEQ
530 static void
531 histeq(fz_pixmap *im)
532 {
533 uint32_t count[256];
534 unsigned char tbl[256];
535 int i;
536 unsigned char *s = im->samples;
537 int n = im->w*im->h;
538 int sigma;
539 int den;
540
541 memset(count, 0, sizeof(count));
542 for (i = n; i > 0; i--)
543 count[*s++]++;
544
545 /* Rather than doing a pure histogram equalisation, we bend
546 * the table so that 0 always stays as 0, and 255 always stays
547 * as 255. */
548 sigma = (count[0]>>1) - count[0];
549 den = sigma + n - (count[255]>>1);
550 for (i = 0; i < 256; i++)
551 {
552 int v = count[i];
553 sigma += v - (v>>1);
554 tbl[i] = (int)(255.0f * sigma / den + 0.5f);
555 sigma += (v>>1);
556 }
557
558 s = im->samples;
559 for (i = n; i > 0; i--)
560 *s = tbl[*s], s++;
561 }
562 #endif
563
564 /* The first functions apply a 5x5 gauss filter to blur the greyscale
565 * image and remove noise. The gauss filter is a convolution with
566 * weights:
567 *
568 * 2 4 5 4 2
569 * 4 9 12 9 4
570 * 5 12 15 12 5
571 * 4 9 12 9 4
572 * 2 4 5 4 2
573 *
574 * As you can see, there are 3 distinct lines of weights within that
575 * matrix. We walk each row of source pixels once, calculating each of
576 * those convolutions, storing the result in a rolling buffer of 5x3xw
577 * entries. Then we sum columns from the buffer to give us the results.
578 */
579
580 /* Read across a row of pixels of width w, from s, performing
581 * the horizontal portion of the convolutions, storing the results
582 * in 3 lines of a buffer, starting at d, &d[w], &d[2*w].
583 *
584 * d[0*w] uses weights (5 12 15 12 5)
585 * d[1*w] uses weights (4 9 12 9 4)
586 * d[2*w] uses weights (2 4 5 4 2)
587 */
588 static void
589 gauss5row(uint16_t *d, const unsigned char *s, int w)
590 {
591 int i;
592 int s0 = s[0];
593 int s1 = s[1];
594 int s2 = s[2];
595 int s3 = s[3];
596
597 s += 4;
598
599 d[2*w] = 11*s0 + 4*s1 + 2*s2;
600 d[1*w] = 25*s0 + 9*s1 + 4*s2;
601 *d++ = 32*s0 + 12*s1 + 5*s2;
602
603 d[2*w] = 6*s0 + 5*s1 + 4*s2 + 2*s3;
604 d[1*w] = 13*s0 + 12*s1 + 9*s2 + 4*s3;
605 *d++ = 17*s0 + 15*s1 + 12*s2 + 5*s3;
606
607 for (i = w - 4; i > 0; i--)
608 {
609 int d2 = 2*s0 + 4*s1 + 5*s2 + 4*s3;
610 int d1 = 4*s0 + 9*s1 + 12*s2 + 9*s3;
611 int d0 = 5*s0 + 12*s1 + 15*s2 + 12*s3;
612 s0 = s1;
613 s1 = s2;
614 s2 = s3;
615 s3 = *s++;
616 d[2*w] = d2 + 2*s3;
617 d[1*w] = d1 + 4*s3;
618 *d++ = d0 + 5*s3;
619 }
620
621 d[2*w] = 2*s0 + 4*s1 + 5*s2 + 6*s3;
622 d[1*w] = 4*s0 + 9*s1 + 12*s2 + 13*s3;
623 *d++ = 5*s0 + 12*s1 + 15*s2 + 17*s3;
624
625 d[2*w] = 2*s1 + 4*s2 + 11*s3;
626 d[1*w] = 4*s1 + 9*s2 + 25*s3;
627 *d = 5*s1 + 12*s2 + 32*s3;
628 }
629
630 /* Calculate the results for row y of the image, of width w, by
631 * summing results from the temporary buffer s, and writing into the
632 * original pixmap at d. */
633 static void
634 gauss5col(unsigned char *d, const uint16_t *s, int y, int w)
635 {
636 const uint16_t *s0, *s1, *s2, *s3, *s4;
637 y *= 3;
638 s0 = &s[((y+ 9+2)%15)*w];
639 s1 = &s[((y+12+1)%15)*w];
640 s2 = &s[( y %15)*w];
641 s3 = &s[((y+ 3+1)%15)*w];
642 s4 = &s[((y+ 6+2)%15)*w];
643
644 for (; w > 0; w--)
645 *d++ = (*s0++ + *s1++ + *s2++ + *s3++ + *s4++ + 79)/159;
646 }
647
648 static void
649 gauss5x5(fz_context *ctx, fz_pixmap *src)
650 {
651 int w = src->w;
652 int h = src->h;
653 uint16_t *buf;
654 unsigned char *s = src->samples;
655 int y;
656
657 if (w < 5 || h < 5)
658 fz_throw(ctx, FZ_ERROR_ARGUMENT, "Pixmap too small");
659
660 buf = fz_malloc(ctx, w*3*5 * sizeof(uint16_t));
661
662 gauss5row(&buf[0*3*w], &s[0*w], w);
663 gauss5row(&buf[1*3*w], &s[1*w], w);
664 memcpy(&buf[3*3*w], buf, w*3 * sizeof(uint16_t)); /* row -2 */
665 memcpy(&buf[4*3*w], buf, w*3 * sizeof(uint16_t)); /* row -1 */
666 for (y = 2; y < h; y++)
667 {
668 gauss5row(&buf[(y%5)*3*w], &s[2*w], w);
669
670 gauss5col(s, buf, y-2, w);
671 s += w;
672 }
673 for (; y < h+2; y++)
674 {
675 memcpy(&buf[(y%5)*3*w], &buf[((y+4)%5)*3*w], 3*w*sizeof(uint16_t));
676 gauss5col(s, buf, y-2, w);
677 s += w;
678 }
679
680 fz_free(ctx, buf);
681 }
682
683 #ifdef DETECT_DOCUMENT_RGB
684 /* Variant of the above that works on a single plane from rgb data */
685 static void
686 gauss5row3(uint16_t *d, const unsigned char *s, int w)
687 {
688 int i;
689 int s0 = s[0];
690 int s1 = s[3];
691 int s2 = s[6];
692 int s3 = s[9];
693
694 s += 4*3;
695
696 d[2*w] = 11*s0 + 4*s1 + 2*s2;
697 d[1*w] = 25*s0 + 9*s1 + 4*s2;
698 *d++ = 32*s0 + 12*s1 + 5*s2;
699
700 d[2*w] = 6*s0 + 5*s1 + 4*s2 + 2*s3;
701 d[1*w] = 13*s0 + 12*s1 + 9*s2 + 4*s3;
702 *d++ = 17*s0 + 15*s1 + 12*s2 + 5*s3;
703
704 for (i = w - 4; i > 0; i--)
705 {
706 int d2 = 2*s0 + 4*s1 + 5*s2 + 4*s3;
707 int d1 = 4*s0 + 9*s1 + 12*s2 + 9*s3;
708 int d0 = 5*s0 + 12*s1 + 15*s2 + 12*s3;
709 s0 = s1;
710 s1 = s2;
711 s2 = s3;
712 s3 = *s; s += 3;
713 d[2*w] = d2 + 2*s3;
714 d[1*w] = d1 + 4*s3;
715 *d++ = d0 + 5*s3;
716 }
717
718 d[2*w] = 2*s0 + 4*s1 + 5*s2 + 6*s3;
719 d[1*w] = 4*s0 + 9*s1 + 12*s2 + 13*s3;
720 *d++ = 5*s0 + 12*s1 + 15*s2 + 17*s3;
721
722 d[2*w] = 2*s1 + 4*s2 + 11*s3;
723 d[1*w] = 4*s1 + 9*s2 + 25*s3;
724 *d = 5*s1 + 12*s2 + 32*s3;
725 }
726
727 static void
728 gauss5x5_3(fz_context *ctx, fz_pixmap *dst, const fz_pixmap *src, int comp)
729 {
730 int w = src->w;
731 int h = src->h;
732 uint16_t *buf;
733 unsigned char *s = src->samples + comp;
734 unsigned char *d = dst->samples;
735 int y;
736
737 if (w < 5 || h < 5)
738 fz_throw(ctx, FZ_ERROR_ARGUMENT, "Pixmap too small");
739
740 buf = fz_malloc(ctx, w*3*5 * sizeof(uint16_t));
741
742 gauss5row3(&buf[0*3*w], s, w);
743 s += w*3;
744 gauss5row3(&buf[1*3*w], s, w);
745 s += w*3;
746 memcpy(&buf[3*3*w], buf, w*3 * sizeof(uint16_t)); /* row -2 */
747 memcpy(&buf[4*3*w], buf, w*3 * sizeof(uint16_t)); /* row -1 */
748 for (y = 2; y < h; y++)
749 {
750 gauss5row3(&buf[((y*3)%15)*w], s, w);
751
752 gauss5col(d, buf, y-2, w);
753 s += w*3;
754 d += w;
755 }
756 for (; y < h+2; y++)
757 {
758 memcpy(&buf[(y%5)*3*w], &buf[((y+4)%5)*3*w], 3*w*sizeof(uint16_t));
759 gauss5col(d, buf, y-2, w);
760 d += w;
761 }
762
763 fz_free(ctx, buf);
764 }
765 #endif
766
767 /* The next set of functions perform the gradient calculation.
768 * We convolve with Sobel kernels Kx and Ky respectively:
769 *
770 * Kx = -1 0 1 Ky = 1 2 1
771 * -2 0 2 0 0 0
772 * -1 0 1 -1 -2 -1
773 *
774 * We do this by using a rolling temporary buffer of int16_t's to hold
775 * 3 pairs of lines of weights scaled by (1 0 1) and (1 2 1).
776 *
777 * We can then sum entries from those lines to calculate kx and ky for
778 * each pixel in the image.
779 *
780 * Then by examining the values of x and y, we can figure out the
781 * "direction" of the edge (horizontal, vertical, or either diagonal),
782 * and the magnitude of the difference across that edge. These get
783 * encoded back into the original image storage using the 2 bottom bits
784 * for direction, and the top 6 bits for magnitude.
785 */
786
787 static void
788 pregradrow(int16_t *d, const unsigned char *s, int w)
789 {
790 int i;
791 unsigned char s0 = *s++;
792 unsigned char s1 = *s++;
793
794 d[w] = 3*s0 + s1;
795 *d++ = s1 - s0;
796 for (i = w-2; i > 0; i--)
797 {
798 int s2 = *s++;
799 d[w] = s0 + 2*s1 + s2;
800 *d++ = s2 - s0;
801 s0 = s1;
802 s1 = s2;
803 }
804 d[w] = s0 + 3*s1;
805 *d = s1 - s0;
806 }
807
808 static void
809 pregradcol(unsigned char *d, const int16_t *buf, int y, int w, uint32_t *max)
810 {
811 const int16_t *s0 = &buf[((y+2)%3)*w*2];
812 const int16_t *s1 = &buf[((y )%3)*w*2];
813 const int16_t *s2 = &buf[((y+1)%3)*w*2];
814 int i;
815
816 for (i = w; i > 0; i--)
817 {
818 uint32_t ax, ay, mag;
819 int x;
820 y = s0[w] - s2[w];
821 x = *s0++ + 2 * *s1++ + *s2++;
822 ax = x >= 0 ? x : -x;
823 ay = y >= 0 ? y : -y;
824 /* x and y are now both in the range -1020..1020 */
825 /* Now we calculate slope and gradient.
826 * angle = atan2(y, x);
827 * intensity = hypot(x, y);
828 * But wait, we don't need that accuracy. We only need
829 * to distinguish 4 directions...
830 *
831 * -22.5 < angle <= 22.5 = 0
832 * 22.5 < angle <= 67.5 = 1
833 * 67.5 < angle <= 112.5 = 2
834 * 115.5 < angle <= 157.5 = 3
835 * (and the reflections)
836 *
837 * 7 0 1 (x positive right, y positive downwards)
838 * 6 * 2
839 * 5 4 3
840 *
841 * tan(22.5)*65536 = 27146.
842 * And for magnitude, we consider the magnitude just
843 * along the 8 directions we've picked. So
844 * 65536/SQR(2) = 46341.
845 * We want magnitude in the 0...63 range.
846 */
847 if (ax<<16 < 27146*ay)
848 {
849 /* angle = 0 */
850 mag = ay<<16; /* 0 to 1020<<16 */
851 }
852 else if (ay<<16 < ax*27146)
853 {
854 /* angle = 2 */
855 mag = ax<<16; /* 0 to 1020<<16 */
856 }
857 else
858 {
859 /* angle = 1 or 3 */
860 mag = (46341*(ax+ay));
861 }
862 if (mag > *max)
863 *max = mag;
864 }
865 }
866
867 static uint32_t
868 pregrad(fz_context *ctx, fz_pixmap *src)
869 {
870 int w = src->w;
871 int h = src->h;
872 unsigned char *s = src->samples;
873 int16_t *buf = fz_malloc(ctx, w*2*3*sizeof(int16_t));
874 int y;
875 uint32_t max = 0;
876
877 pregradrow(buf, s, w); /* Line 0 */
878 memcpy(&buf[w*2*2], buf, w*2*sizeof(int16_t)); /* Line 1 */
879 s += w;
880 for (y = 1; y < h-1; y++)
881 {
882 pregradrow(&buf[(y%3)*w*2], s, w);
883 pregradcol(s-w, buf, y-1, w, &max);
884 s += w;
885 }
886 memcpy(&buf[((y+1)%3)*w*2], &buf[(y%3)*w*2], w*2*sizeof(int16_t)); /* Line h */
887 pregradcol(s-w, buf, h-2, w, &max);
888 pregradcol(s, buf, h-1, w, &max);
889
890 fz_free(ctx, buf);
891
892 if (max == 0)
893 return 1;
894 else
895 return 0x7FFFFFFFU/max;
896 }
897
898
899 static void
900 gradrow(int16_t *d, const unsigned char *s, int w)
901 {
902 int i;
903 unsigned char s0 = *s++;
904 unsigned char s1 = *s++;
905
906 d[w] = 3*s0 + s1;
907 *d++ = s1 - s0;
908 for (i = w-2; i > 0; i--)
909 {
910 int s2 = *s++;
911 d[w] = s0 + 2*s1 + s2;
912 *d++ = s2 - s0;
913 s0 = s1;
914 s1 = s2;
915 }
916 d[w] = s0 + 3*s1;
917 *d = s1 - s0;
918 }
919
920 static void
921 gradcol(unsigned char *d, const int16_t *buf, int y, int w, int scale)
922 {
923 const int16_t *s0 = &buf[((y+2)%3)*w*2];
924 const int16_t *s1 = &buf[((y )%3)*w*2];
925 const int16_t *s2 = &buf[((y+1)%3)*w*2];
926 int i;
927
928 for (i = w; i > 0; i--)
929 {
930 uint32_t ax, ay, mag, scaled;
931 int angle, x;
932 y = s0[w] - s2[w];
933 x = *s0++ + 2 * *s1++ + *s2++;
934 ax = x >= 0 ? x : -x;
935 ay = y >= 0 ? y : -y;
936 /* x and y are now both in the range -1020..1020 */
937 /* Now we calculate slope and gradient.
938 * angle = atan2(y, x);
939 * intensity = hypot(x, y);
940 * But wait, we don't need that accuracy. We only need
941 * to distinguish 4 directions...
942 *
943 * -22.5 < angle <= 22.5 = 0
944 * 22.5 < angle <= 67.5 = 1
945 * 67.5 < angle <= 112.5 = 2
946 * 115.5 < angle <= 157.5 = 3
947 * (and the reflections)
948 *
949 * 7 0 1 (x positive right, y positive downwards)
950 * 6 * 2
951 * 5 4 3
952 *
953 * tan(22.5)*65536 = 27146.
954 * And for magnitude, we consider the magnitude just
955 * along the 8 directions we've picked. So
956 * 65536/SQR(2) = 46341.
957 * We want magnitude in the 0...63 range.
958 */
959 if (ax<<16 < 27146*ay)
960 {
961 angle = 0;
962 mag = ay<<16; /* 0 to 1020<<16 */
963 }
964 else if (ay<<16 < ax*27146)
965 {
966 angle = 2;
967 mag = ax<<16; /* 0 to 1020<<16 */
968 }
969 else
970 {
971 /* 1 or 3 */
972 angle = (x^y) >= 0 ? 3 : 1;
973 mag = (46341*(ax+ay));
974 }
975 scaled = (mag * scale)>>25;
976 assert(scaled >= 0 && scaled <= 63);
977 *d++ = (scaled<<2) | angle;
978 }
979 }
980
981 static void
982 grad(fz_context *ctx, fz_pixmap *src, uint32_t scale)
983 {
984 int w = src->w;
985 int h = src->h;
986 unsigned char *s = src->samples;
987 int16_t *buf = fz_malloc(ctx, w*2*3*sizeof(int16_t));
988 int y;
989
990 gradrow(buf, s, w); /* Line 0 */
991 memcpy(&buf[w*2*2], buf, w*2*sizeof(int16_t)); /* Line 1 */
992 s += w;
993 for (y = 1; y < h-1; y++)
994 {
995 gradrow(&buf[(y%3)*w*2], s, w);
996 gradcol(s-w, buf, y-1, w, scale);
997 s += w;
998 }
999 memcpy(&buf[((y+1)%3)*w*2], &buf[(y%3)*w*2], w*2*sizeof(int16_t)); /* Line h */
1000 gradcol(s-w, buf, h-2, w, scale);
1001 gradcol(s, buf, h-1, w, scale);
1002
1003 fz_free(ctx, buf);
1004 }
1005
1006 #ifdef DETECT_DOCUMENT_RGB
1007 static void
1008 combine_grad(fz_pixmap *grey, const fz_pixmap *r, const fz_pixmap *g, const fz_pixmap *b)
1009 {
1010 int n;
1011 unsigned char *sd = grey->samples;
1012 const unsigned char *sr = r->samples;
1013 const unsigned char *sg = g->samples;
1014 const unsigned char *sb = b->samples;
1015
1016 for (n = g->w * g->h; n > 0; n--)
1017 {
1018 unsigned char vg = *sg++;
1019 unsigned char vr = *sr++;
1020 unsigned char vb = *sb++;
1021 unsigned char vd = *sd++;
1022 if (vr > vg)
1023 vg = vr;
1024 if (vb > vg)
1025 vg = vb;
1026 if (vg > vd)
1027 sd[-1] = vg;
1028 }
1029 }
1030 #endif
1031
1032 /* Next, we perform Non-Maximum Suppression and Double Thresholding,
1033 * both in the same phase.
1034 *
1035 * We walk the image, looking at the magnitude of the edges. Edges below
1036 * the 'weak' threshold are discarded. Otherwise, neighbouring pixels in
1037 * the direction of the edge are considered; if other pixels are stronger
1038 * then this pixel is discarded. If not, we classify ourself as either
1039 * 'strong' or 'weak'.
1040 */
1041 #define WEAK_EDGE 64
1042 #define STRONG_EDGE 128
1043 static void
1044 nonmax(fz_context *ctx, fz_pixmap *dst, const fz_pixmap *src, int pass)
1045 {
1046 int w = src->w;
1047 int h = src->h;
1048 const unsigned char *s0 = src->samples;
1049 const unsigned char *s1 = s0;
1050 const unsigned char *s2 = s0+w;
1051 unsigned char *d = dst->samples;
1052 int x, y;
1053 /* thresholds are in the 0 to 63 range.
1054 * WEAK is typically 0.1ish, STRONG 0.3ish
1055 */
1056 int weak = 6 - pass;
1057 int strong = 12 - pass*2;
1058
1059 /* On entry, pixels have the angle in the bottom 2 bits and the magnitude in the rest. */
1060 /* On exit, strong pixels have bit 7 set, weak pixels have bit 6, others are 0.
1061 * strong and weak pixels have the angle in bits 4 and 5. */
1062
1063 for (y = h-1; y >= 0;)
1064 {
1065 int lastmag;
1066 int ang = *s1++;
1067 int mag = ang>>2;
1068 int q, r;
1069
1070 /* Pixel 0 */
1071 if (mag <= weak)
1072 {
1073 /* Not even a weak edge. We'll never keep it. */
1074 *d++ = 0;
1075 s0++;
1076 s2++;
1077 }
1078 else
1079 {
1080 ang &= 3;
1081 switch (ang)
1082 {
1083 default:
1084 case 0:
1085 q = (*s0++)>>2;
1086 r = (*s2++)>>2;
1087 break;
1088 case 1:
1089 q = (*++s0)>>2;
1090 r = 0;
1091 s2++;
1092 break;
1093 case 2:
1094 s0++;
1095 s2++;
1096 q = 0;
1097 r = *s1>>2;
1098 break;
1099 case 3:
1100 q = 0;
1101 s0++;
1102 r = (*++s2)>>2;
1103 break;
1104 }
1105 if (mag < q || mag < r)
1106 {
1107 /* Neighbouring edges are stronger.
1108 * Lose this one. */
1109 *d++ = 0;
1110 }
1111 else if (mag < strong)
1112 {
1113 /* Weak edge. */
1114 *d++ = WEAK_EDGE | (ang<<4);
1115 }
1116 else
1117 {
1118 /* Strong edge */
1119 *d++ = STRONG_EDGE | (ang<<4);
1120 }
1121 }
1122 lastmag = mag;
1123 for (x = w-2; x > 0; x--)
1124 {
1125 ang = *s1++;
1126 mag = ang>>2;
1127 if (mag <= weak)
1128 {
1129 /* Not even a weak edge. We'll never keep it. */
1130 *d++ = 0;
1131 s0++;
1132 s2++;
1133 }
1134 else
1135 {
1136 ang &= 3;
1137 switch (ang)
1138 {
1139 default:
1140 case 0:
1141 q = (*s0++)>>2;
1142 r = (*s2++)>>2;
1143 break;
1144 case 1:
1145 q = (*++s0)>>2;
1146 r = ((s2++)[-1])>>2;
1147 break;
1148 case 2:
1149 s0++;
1150 s2++;
1151 q = lastmag;
1152 r = *s1>>2;
1153 break;
1154 case 3:
1155 q = ((s0++)[-1])>>2;
1156 r = (*++s2)>>2;
1157 break;
1158 }
1159 if (mag < q || mag < r)
1160 {
1161 /* Neighbouring edges are stronger.
1162 * Lose this one. */
1163 *d++ = 0;
1164 }
1165 else if (mag < strong)
1166 {
1167 /* Weak edge. */
1168 *d++ = WEAK_EDGE | (ang<<4);
1169 }
1170 else
1171 {
1172 /* Strong edge */
1173 *d++ = STRONG_EDGE | (ang<<4);
1174 }
1175 }
1176 lastmag = mag;
1177 }
1178 /* Pixel w-1 */
1179 ang = *s1++;
1180 mag = ang>>2;
1181 if (mag <= weak)
1182 {
1183 /* Not even a weak edge. We'll never keep it. */
1184 *d++ = 0;
1185 s0++;
1186 s2++;
1187 lastmag = 0;
1188 }
1189 else
1190 {
1191 ang &= 3;
1192 switch (ang)
1193 {
1194 default:
1195 case 0:
1196 q = (*s0++)>>2;
1197 r = (*s2++)>>2;
1198 break;
1199 case 1:
1200 q = 0;
1201 s0++;
1202 r = ((s2++)[-1])>>2;
1203 break;
1204 case 2:
1205 s0++;
1206 s2++;
1207 q = 0;
1208 r = *s1>>2;
1209 break;
1210 case 3:
1211 q = ((s0++)[-1])>>2;
1212 r = 0;
1213 s2++;
1214 break;
1215 }
1216 if (mag < q || mag < r)
1217 {
1218 /* Neighbouring edges are stronger.
1219 * Lose this one. */
1220 *d++ = 0;
1221 }
1222 else if (mag < strong)
1223 {
1224 /* Weak edge. */
1225 *d++ = WEAK_EDGE | (ang<<4);
1226 }
1227 else
1228 {
1229 /* Strong edge */
1230 *d++ = STRONG_EDGE | (ang<<4);
1231 }
1232 }
1233 s0 = s1-w;
1234 if (--y == 0)
1235 s2 = s1;
1236 }
1237 }
1238
1239 /* Next, we have the hysteresis phase. Here we bump any 'weak' pixel
1240 * that has at least one strong pixel around it up to being a 'strong'
1241 * pixel.
1242 *
1243 * On entry, strong pixels have bit 7 set, weak pixels have bit 6 set,
1244 * the angle is in bits 4 and 5, everything else is 0.
1245 *
1246 * We walk the rows of the image, and for each pixel we set bit 0 to be
1247 * the logical OR of all bit 7's of itself, and its horizontally
1248 * neighbouring pixels.
1249 *
1250 * Once we have done the first 2 rows like that, we can combine the
1251 * operation of generating the bottom bits for row i, with the
1252 * calculation of row i-1. Any given pixel on row i-1 should be promoted
1253 * to 'strong' if it was 'weak', and if the logical OR of the matching
1254 * pixels in row i, i-1 or i-2 has bit 0 set.
1255 *
1256 * At the end of this process any pixel with bit 7 set is 'strong'.
1257 * Bits 4 and 5 still have the angle.
1258 */
1259
1260 static void
1261 hysteresis(fz_context *ctx, fz_pixmap *src)
1262 {
1263 int w = src->w;
1264 int h = src->h;
1265 unsigned char *s0 = src->samples;
1266 unsigned char *s1 = s0;
1267 unsigned char *s2 = s0;
1268 unsigned char v0, v1, v2, r0, r1, r2;
1269 int x, y;
1270
1271 /* On entry, strong pixels have bit 7 set, weak pixels have bit 6, others are 0.
1272 * strong and weak pixels have the angle in bits 4 and 5. */
1273
1274 /* We make the bottom bit in every pixel be 1 iff the pixel
1275 * or the ones to either side of it are 'strong'. */
1276
1277 /* First row - just do the bottom bit. */
1278 /* Pixel 0 */
1279 v0 = *s0++;
1280 v1 = *s0++;
1281 s0[-2] = v0 | ((v0 | v1)>>7);
1282 /* Middle pixels */
1283 for (x = w-2; x > 0; x--)
1284 {
1285 v2 = *s0++;
1286 s0[-2] = v1 | ((v0 | v1 | v2)>>7);
1287 v0 = v1;
1288 v1 = v2;
1289 }
1290 /* Pixel w-1 */
1291 s0[-1] = v1 | ((v0 | v1)>>7);
1292 assert(s0 == src->samples + w);
1293
1294 /* Second row - do the "bottom bit" for the second row, and
1295 * perform hysteresis on the top row. */
1296 /* Pixel 0 */
1297 v0 = *s0++;
1298 v1 = *s0++;
1299 r0 = v0 | ((v0 | v1)>>7);
1300 s0[-2] = r0;
1301 r1 = *s1++;
1302 if ((r1>>6) & (r0 | r1))
1303 s1[-1] |= 128;
1304 /* Middle pixels */
1305 for (x = w-2; x > 0; x--)
1306 {
1307 v2 = *s0++;
1308 r0 = v1 | ((v0 | v1 | v2)>>7);
1309 s0[-2] = r0;
1310 r1 = *s1++;
1311 if ((r1>>6) & (r0 | r1))
1312 s1[-1] |= 128;
1313 v0 = v1;
1314 v1 = v2;
1315 }
1316 /* Pixel w-1 */
1317 r0 = v1 | ((v0 | v1)>>7);
1318 s0[-1] = r0;
1319 r1 = *s1++;
1320 if ((r1>>6) & (r0 | r1))
1321 s1[-1] |= 128;
1322 assert(s0 == s1 + w);
1323
1324 /* Now we get into the swing of things. We do the "bottom bit"
1325 * for row n+1, and do the actual processing for row n. */
1326 for (y = h-4; y > 0; y--)
1327 {
1328 /* Pixel 0 */
1329 v0 = *s0++;
1330 v1 = *s0++;
1331 r0 = v0 | ((v0 | v1)>>7);
1332 s0[-2] = r0;
1333 r1 = *s1++;
1334 r2 = *s2++;
1335 if ((r1>>6) & (r0 | r1 | r2))
1336 s1[-1] |= 128;
1337 /* Middle pixels */
1338 for (x = w-2; x > 0; x--)
1339 {
1340 v2 = *s0++;
1341 r0 = v1 | ((v0 | v1 | v2)>>7);
1342 s0[-2] = r0;
1343 r1 = *s1++;
1344 r2 = *s2++;
1345 if ((r1>>6) & (r0 | r1 | r2))
1346 s1[-1] |= 128;
1347 v0 = v1;
1348 v1 = v2;
1349 }
1350 /* Pixel w-1 */
1351 r0 = v1 | ((v0 | v1)>>7);
1352 s0[-1] = r0;
1353 r1 = *s1++;
1354 r2 = *s2++;
1355 if ((r1>>6) & (r0 | r1 | r2))
1356 s1[-1] |= 128;
1357 assert(s0 == s1 + w);
1358 assert(s1 == s2 + w);
1359 }
1360
1361 /* Final 2 rows together */
1362 /* Pixel 0 */
1363 v0 = *s0++;
1364 v1 = *s0++;
1365 r0 = v0 | ((v0 | v1)>>7);
1366 r1 = *s1++;
1367 r2 = *s2++;
1368 if ((r1>>6) & (r0 | r1 | r2))
1369 s1[-1] |= 128;
1370 if ((r0>>6) & (r0 | r1))
1371 s0[-2] |= 128;
1372 /* Middle pixels */
1373 for (x = w-2; x > 0; x--)
1374 {
1375 v2 = *s0++;
1376 r0 = v1 | ((v0 | v1 | v2)>>7);
1377 r1 = *s1++;
1378 r2 = *s2++;
1379 if ((r1>>6) & (r0 | r1 | r2))
1380 s1[-1] |= 128;
1381 if ((r0>>6) & (r0 | r1))
1382 s0[-2] |= 128;
1383 v0 = v1;
1384 v1 = v2;
1385 }
1386 /* Pixel w-1 */
1387 r0 = v1 | ((v0 | v1)>>7);
1388 r1 = *s1;
1389 r2 = *s2;
1390 if ((r1>>6) & (r0 | r1 | r2))
1391 s1[0] |= 128;
1392 if ((r0>>6) & (r0 | r1))
1393 s0[-1] |= 128;
1394 }
1395
1396 #ifdef WARP_DEBUG
1397 /* A simple function to keep just bit 7 of the image. This 'cleans' the
1398 * pixmap so that a grey version can be saved out for visual checking. */
1399 static void
1400 clean(fz_context *ctx, fz_pixmap *src)
1401 {
1402 int w = src->w;
1403 int h = src->h;
1404 unsigned char *s = src->samples;
1405
1406 for (w = w*h; w > 0; w--)
1407 *s = *s & 128, s++;
1408 }
1409 #endif
1410
1411 #define SINTABLE_SHIFT 14
1412 static int16_t sintable[270];
1413 #define costable (&sintable[90])
1414
1415 /* We have collected an array of edge data.
1416 * For each pixel, we know whether there is a 'strong' edge
1417 * there, and if so, in which of 4 directions it runs.
1418 *
1419 * We want to convert this into hough space.
1420 *
1421 * The literature describes points in Hough space as having the
1422 * form (r, theta). The value of any given point point in hough
1423 * space is the "strength" of the line described by:
1424 * x.cos(theta) + y.sin(theta) = r
1425 *
1426 * | \
1427 * | /\
1428 * | /\/\
1429 * | r/ \
1430 * | / \
1431 * | / \
1432 * |/theta \
1433 * -+--------------
1434 *
1435 * i.e. r is the shortest distance from the origin to the line,
1436 * and theta gives us the angle of that shortest line.
1437 *
1438 * But, we are using angles of theta from the vertical, so we need
1439 * a different formulation:
1440 *
1441 * | \
1442 * | /\
1443 * | /\/\
1444 * | r/ \ t = theta
1445 * | / \
1446 * |t/ \
1447 * |/ \
1448 * -+--------------
1449 *
1450 * So we're using 90-theta. cos(90-theta) = sin(theta),
1451 * and sin(90-theta) = cos(theta).
1452 *
1453 * So: x.sin(theta) + y.cos(theta) = r (for theta measured
1454 * clockwise from the y axis).
1455 *
1456 * We've been collecting angles according to their position in one
1457 * of 4 octants:
1458 *
1459 * Ang 0 = close to a horizontal edge (-22.5 to 22.5 degrees)
1460 * Ang 1 = close to diagonal edge (top left to bottom right) (22.5 to 67.5 degrees)
1461 * Ang 2 = close to a vertical edge (67.5 to 112.5 degrees)
1462 * Ang 3 = close to diagonal edge (bottom left to top right) (112.5 to 157.5 degrees)
1463 *
1464 * The other 4 octants mirror onto these.
1465 *
1466 * So, for each point in our (x,y) pixmap we whether we have a strong
1467 * pixel or not. If we have such a pixel, then we know that an edge
1468 * passes through that point, and which of those 4 octants it is in.
1469 *
1470 * We therefore consider all the possible angles within that octant,
1471 * and add a 'vote' for each of those lines into our hough transformed
1472 * space.
1473 */
1474
1475
1476 static void
1477 mark_hough(uint32_t *d, int x, int y, int maxlen, int reduce, int ang)
1478 {
1479 int theta;
1480 int stride = (maxlen*2)>>reduce;
1481 int minang, maxang;
1482
1483 switch (ang)
1484 {
1485 /* The angles are really 22.5, 67.5 etc, but we are working in ints
1486 * and specifying maxang as the first one greater than the one we
1487 * want to mark. */
1488 default:
1489 case 0:
1490 /* Vertical boundary. Lines through this boundary
1491 * go horizontally. So the perpendicular to them
1492 * is vertical. */
1493 minang = 0; maxang = 23;
1494 break;
1495 case 1:
1496 /* NE boundary */
1497 minang = 23; maxang = 68;
1498 break;
1499 case 2:
1500 /* Horizontal boundary. */
1501 minang = 68; maxang = 113;
1502 break;
1503 case 3:
1504 /* SE boundary */
1505 minang = 113; maxang = 158;
1506 break;
1507 case 4:
1508 /* For debugging: */
1509 minang = 0; maxang = 180;
1510 break;
1511 }
1512
1513 d += minang * stride;
1514 while (1)
1515 {
1516 for (theta = minang; theta < maxang; theta++)
1517 {
1518 int p = (x*sintable[theta] + y*costable[theta])>>SINTABLE_SHIFT;
1519 int v = (maxlen + p)>>reduce;
1520
1521 d[v]++;
1522 d += stride;
1523 }
1524 if (ang != 0)
1525 break;
1526 ang = 4;
1527 minang = 158;
1528 d += (minang - maxang) * stride;
1529 maxang = 180;
1530 }
1531 }
1532
1533 #ifdef WARP_DEBUG
1534 static void
1535 save_hough_debug(fz_context *ctx, uint32_t *hough, int stride)
1536 {
1537 uint32_t scale;
1538 uint32_t maxval;
1539 uint32_t *p;
1540 int y;
1541 fz_pixmap *dst = fz_new_pixmap(ctx, NULL, stride, 180, NULL, 0);
1542 unsigned char *d = dst->samples;
1543 /* Make the image of the hough space (for debugging) */
1544 maxval = 1; /* Avoid possible division by zero */
1545 p = hough;
1546 for (y = 180*stride; y > 0; y--)
1547 {
1548 uint32_t v = *p++;
1549 if (v > maxval)
1550 maxval = v;
1551 }
1552
1553 scale = 0xFFFFFFFFU/maxval;
1554 p = hough;
1555 for (y = 180*stride; y > 0; y--)
1556 {
1557 *d++ = (scale * *p++)>>24;
1558 }
1559 fz_save_pixmap_as_png(ctx, dst, "hough.png");
1560 fz_drop_pixmap(ctx, dst);
1561 }
1562 #endif
1563
1564 static uint32_t *do_hough(fz_context *ctx, const fz_pixmap *src, int stride, int maxlen, int reduce)
1565 {
1566 int w = src->w;
1567 int h = src->h;
1568 int x, y;
1569 const unsigned char *s = src->samples;
1570 uint32_t *hough = fz_calloc(ctx, sizeof(uint32_t), 180*stride);
1571
1572 START_TIME();
1573 /* Construct the hough space representation. */
1574 for (y = 0; y < h; y++)
1575 {
1576 for (x = 0; x < w; x++)
1577 {
1578 unsigned char v = *s++;
1579 if (v & 128)
1580 mark_hough(hough, x, y, maxlen, reduce, (v>>4)&3);
1581 }
1582 }
1583 END_TIME("Building hough");
1584
1585 #ifdef WARP_DEBUG
1586 save_hough_debug(ctx, hough, stride);
1587 #endif
1588
1589 return hough;
1590 }
1591
1592 typedef struct
1593 {
1594 int ang;
1595 int dis;
1596 int strength;
1597 int x0;
1598 int y0;
1599 int x1;
1600 int y1;
1601 } hough_edge_t;
1602
1603 typedef struct
1604 {
1605 int e0;
1606 int e1;
1607 int score;
1608 float x;
1609 float y;
1610 } hough_point_t;
1611
1612 static int
1613 intersect(hough_point_t *pt, hough_edge_t *edges, int e0, int e1)
1614 {
1615 int x1 = edges[e0].x0;
1616 int y1 = edges[e0].y0;
1617 int x2 = edges[e0].x1;
1618 int y2 = edges[e0].y1;
1619 int x3 = edges[e1].x0;
1620 int y3 = edges[e1].y0;
1621 int x4 = edges[e1].x1;
1622 int y4 = edges[e1].y1;
1623 int d = (x1-x2)*(y3-y4)-(y1-y2)*(x3-x4);
1624 int t = (x1-x3)*(y3-y4)-(y1-y3)*(x3-x4);
1625 float ft;
1626
1627 if (d < 0)
1628 {
1629 if (t < d || t > 0)
1630 return 0;
1631 }
1632 else
1633 {
1634 if (t < 0 || t > d)
1635 return 0;
1636 }
1637
1638 pt->e0 = e0;
1639 pt->e1 = e1;
1640 pt->score = edges[e0].strength + edges[e1].strength;
1641 ft = t / (float)d;
1642 pt->x = x1 + ft*(x2-x1);
1643 pt->y = y1 + ft*(y2-y1);
1644
1645 return 1;
1646 }
1647
1648 typedef struct
1649 {
1650 int w;
1651 int h;
1652 int edge[4];
1653 int point[4];
1654 int best_score;
1655 int best_edge[4];
1656 int best_point[4];
1657 } hough_route_t;
1658
1659 /* To test convexity, we use the dot product:
1660 * B
1661 * / \
1662 * / C
1663 * A
1664 *
1665 * The dot product of vectors u and v is:
1666 * dot(u,v) = |u|.|v|.cos(theta)
1667 * Helpfully, cos(theta) > 0 for -90 < theta < 90.
1668 * So if we can form perp(AB) = the vector given by rotating
1669 * AB 90 degrees clockwise, we can then look at the sign of:
1670 * dot(perp(AB),BC)
1671 * If we do that for each corner of the polygon, and the signs
1672 * of the results for all of them are the same, we have convexity.
1673 *
1674 * If any of the dot products are zero, the edges are parallel
1675 * (i.e. the same edge). This should never happen.
1676 *
1677 * If AB = (x,y), then perp(AB) = (y,-x)
1678 */
1679 static void
1680 score_corner(hough_point_t *points, hough_route_t *route, int p0, int *score, int *sign)
1681 {
1682 float x0, x1, x2, y0, y1, y2, ux, uy, vx, vy, dot;
1683 int s, len;
1684 float costheta;
1685
1686 /* If we've already decided this is duff, then bale. */
1687 if (*score < 0)
1688 return;
1689
1690 x0 = points[route->point[p0]].x;
1691 y0 = points[route->point[p0]].y;
1692 x1 = points[route->point[(p0+1)&3]].x;
1693 y1 = points[route->point[(p0+1)&3]].y;
1694 x2 = points[route->point[(p0+2)&3]].x;
1695 y2 = points[route->point[(p0+2)&3]].y;
1696
1697 ux = y1-y0; /* u = Perp(p0 to p1)*/
1698 uy = x0-x1;
1699 vx = x2-x1; /* v = p1 to p2 */
1700 vy = y2-y1;
1701
1702 dot = ux*vx + uy*vy;
1703 if (dot == 0)
1704 {
1705 *score = -1;
1706 return;
1707 }
1708
1709 s = (dot > 0) ? 1 : -1;
1710 if (*sign == 0)
1711 *sign = s;
1712 else if (*sign != s)
1713 {
1714 *score = -1;
1715 return;
1716 }
1717
1718 len = sqrt(ux*ux + uy*uy) * sqrt(vx*vx + vy*vy);
1719 costheta = dot / (float)len;
1720 if (costheta < 0)
1721 costheta = -costheta;
1722 if (costheta < 0.7)
1723 {
1724 *score = -1;
1725 return;
1726 }
1727 costheta *= costheta;
1728
1729 *score += points[route->point[(p0+1)%3]].score * costheta;
1730 }
1731
1732 /* route points to 8 ints:
1733 * 2*i = edge number
1734 * 2*i+1 = point number
1735 */
1736 static int
1737 score_route(hough_point_t *points, hough_route_t *route)
1738 {
1739 int score = 0;
1740 int sign = 0;
1741
1742 score_corner(points, route, 0, &score, &sign);
1743 score_corner(points, route, 1, &score, &sign);
1744 score_corner(points, route, 2, &score, &sign);
1745 score_corner(points, route, 3, &score, &sign);
1746
1747 return score;
1748 }
1749
1750 static float
1751 score_by_area(const hough_point_t *points, const hough_route_t *route)
1752 {
1753 float double_area_of_quad =
1754 points[route->point[0]].x * points[route->point[1]].y +
1755 points[route->point[1]].x * points[route->point[2]].y +
1756 points[route->point[2]].x * points[route->point[3]].y +
1757 points[route->point[3]].x * points[route->point[0]].y -
1758 points[route->point[1]].x * points[route->point[0]].y -
1759 points[route->point[2]].x * points[route->point[1]].y -
1760 points[route->point[3]].x * points[route->point[2]].y -
1761 points[route->point[0]].x * points[route->point[3]].y;
1762 float double_area = route->w * (float)route->h * 2;
1763 if (double_area_of_quad < 0)
1764 double_area_of_quad = -double_area_of_quad;
1765 /* Anything larger than a quarter of the screen is acceptable. */
1766 if (double_area_of_quad*4 > double_area)
1767 return 1;
1768 /* Anything smaller than a 16th of the screen is unacceptable in all circumstances. */
1769 if (double_area_of_quad*16 < double_area)
1770 return 0;
1771
1772 /* Otherwise, scale the score down by how much it's less than a quarter of the screen. */
1773 return (double_area_of_quad*4)/double_area;
1774 }
1775
1776 /* The first n+1 edges of the route are filled in, as are the first n
1777 * points.
1778 * 2*i = edge number
1779 * 2*i+1 = point number
1780 */
1781 static void
1782 find_route(fz_context *ctx, hough_point_t *points, int num_points, hough_route_t *route, int n)
1783 {
1784 int i;
1785
1786 for (i = 0; i < num_points; i++)
1787 {
1788 /* If this point continues our route (e0 == route->edge[n])
1789 * then the next point in the route might be e1. */
1790 int e0 = points[i].e0;
1791 int e1 = points[i].e1;
1792 if (e0 == route->edge[n])
1793 {
1794 }
1795 else if (e1 == route->edge[n])
1796 {
1797 int t = e0; e0 = e1; e1 = t;
1798 }
1799 else
1800 continue; /* Doesn't fit. Keep searching. */
1801
1802 /* If we've found 3 points already, then this is the
1803 * fourth. */
1804 if (n == 3)
1805 {
1806 int score = 0;
1807
1808 /* If we haven't looped back to our first edge,
1809 * no good. Keep searching. */
1810 if (route->edge[0] != e1)
1811 continue;
1812 route->point[3] = i;
1813
1814 score = score_route(points, route);
1815
1816 if (score > 0)
1817 score *= score_by_area(points, route);
1818
1819 #ifdef WARP_DEBUG
1820 debug_printf(ctx, "Found route: (point=%d %d %d %d) (edge %d %d %d %d) score=%d\n",
1821 route->point[0], route->point[1], route->point[2], route->point[3],
1822 route->edge[0], route->edge[1], route->edge[2], route->edge[3],
1823 score);
1824 #endif
1825
1826 /* We want our route to be convex */
1827 if (score < 0)
1828 continue;
1829
1830 /* Score the route */
1831 if (route->best_score < score)
1832 {
1833 route->best_score = score;
1834 memcpy(route->best_edge, route->edge, sizeof(*route->edge)*4);
1835 memcpy(route->best_point, route->point, sizeof(*route->point)*4);
1836 }
1837 }
1838 else
1839 {
1840 int j;
1841 for (j = 0; j < n && route->edge[j] != e1; j++);
1842 /* If we're about to loop back to any of the
1843 * previous edges, keep searching. */
1844 if (j != n)
1845 continue;
1846
1847 /* Possible. Extend the route, and recurse. */
1848 route->point[n] = i;
1849 route->edge[n+1] = e1;
1850 find_route(ctx, points, num_points, route, n+1);
1851 }
1852 }
1853 }
1854
1855 #define MAX_EDGES 32
1856 #define BLOT_ANG 10
1857 #define BLOT_DIS 10
1858 static int
1859 make_hough(fz_context *ctx, const fz_pixmap *src, fz_quad *corners)
1860 {
1861 int w = src->w;
1862 int h = src->h;
1863 int maxlen = (int)(sqrtf(w*w + h*h) + 0.5);
1864 uint32_t *hough;
1865 int x, y;
1866 int reduce;
1867 int stride;
1868 hough_edge_t edge[MAX_EDGES];
1869 hough_point_t points[MAX_EDGES * MAX_EDGES];
1870 hough_route_t route;
1871 int num_edges, num_points;
1872
1873 /* costable could (should) be statically inited. */
1874 {
1875 int i;
1876 for (i = 0; i < 270; i++)
1877 {
1878 float theta = i*M_PI/180;
1879 sintable[i] = (int16_t)((1<<SINTABLE_SHIFT)*sinf(theta) + 0.5f);
1880 }
1881 }
1882
1883 /* Figure out a suitable scale for the data. */
1884 reduce = 0;
1885 while ((maxlen*2>>reduce) > 720 && reduce < 16)
1886 reduce++;
1887
1888 stride = (maxlen*2)>>reduce;
1889 hough = do_hough(ctx, src, stride, maxlen, reduce);
1890
1891 /* We want to find the top n edges that aren't too close to
1892 * one another. */
1893 for (x = 0; x < MAX_EDGES; x++)
1894 {
1895 int ang, dis;
1896 int minang, maxang, mindis, maxdis;
1897 int where = 0;
1898 uint32_t *p = hough;
1899 uint32_t maxval = 0;
1900 for (y = 180*stride; y > 0; y--)
1901 {
1902 uint32_t v = *p++;
1903 if (v > maxval)
1904 maxval = v, where = y;
1905 }
1906 if (where == 0)
1907 break;
1908 where = 180*stride - where;
1909 ang = edge[x].ang = where/stride;
1910 dis = edge[x].dis = where - edge[x].ang*stride;
1911 edge[x].strength = hough[where];
1912 /* We don't want to find any other maxima that are too
1913 * close to this one, so we 'blot out' stuff around this
1914 * maxima. */
1915 #ifdef WARP_DEBUG
1916 debug_printf(ctx, "Maxima %d: dist=%d ang=%d strength=%d\n",
1917 x, (dis<<reduce)-maxlen, ang-90, edge[x].strength);
1918 #endif
1919 minang = ang - BLOT_ANG;
1920 if (minang < 0)
1921 minang = 0;
1922 maxang = ang + BLOT_ANG;
1923 if (maxang > 180)
1924 maxang = 180;
1925 mindis = dis - BLOT_DIS;
1926 if (mindis < 0)
1927 mindis = 0;
1928 maxdis = dis + BLOT_DIS;
1929 if (maxdis > stride)
1930 maxdis = stride;
1931 p = hough + minang*stride + mindis;
1932 maxdis = (maxdis-mindis)*sizeof(uint32_t);
1933 for (y = maxang-minang; y > 0; y--)
1934 {
1935 memset(p, 0, maxdis);
1936 p += stride;
1937 }
1938 #ifdef WARP_DEBUG
1939 //save_hough_debug(ctx, hough, stride);
1940 #endif
1941 }
1942 num_edges = x;
1943 if (num_edges == 0)
1944 return 0;
1945
1946 /* Find edges in terms of lines. */
1947 for (x = 0; x < num_edges; x++)
1948 {
1949 int ang = edge[x].ang;
1950 int dis = edge[x].dis;
1951 int p = (dis<<reduce) - maxlen;
1952 if (ang < 45 || ang > 135)
1953 {
1954 /* Mostly horizontal line */
1955 edge[x].x0 = 0;
1956 edge[x].x1 = w;
1957 edge[x].y0 = ((p<<SINTABLE_SHIFT) - edge[x].x0*sintable[ang])/costable[ang];
1958 edge[x].y1 = ((p<<SINTABLE_SHIFT) - edge[x].x1*sintable[ang])/costable[ang];
1959 }
1960 else
1961 {
1962 /* Mostly vertical line */
1963 edge[x].y0 = 0;
1964 edge[x].y1 = h;
1965 edge[x].x0 = ((p<<SINTABLE_SHIFT) - edge[x].y0*costable[ang])/sintable[ang];
1966 edge[x].x1 = ((p<<SINTABLE_SHIFT) - edge[x].y1*costable[ang])/sintable[ang];
1967 }
1968 }
1969
1970 /* Find the points of intersection */
1971 num_points = 0;
1972 for (x = 0; x < num_edges-1; x++)
1973 for (y = x+1; y < num_edges; y++)
1974 num_points += intersect(&points[num_points],
1975 edge, x, y);
1976
1977 #ifdef WARP_DEBUG
1978 {
1979 debug_printf(ctx, "%d edges, %d points\n", num_edges, num_points);
1980 for (x = 0; x < num_points; x++)
1981 {
1982 debug_printf(ctx, "p%d: %d %d (score %d, %d+%d)\n", x,
1983 (int)points[x].x, (int)points[x].y, points[x].score,
1984 points[x].e0, points[x].e1);
1985 }
1986 }
1987 #endif
1988
1989 /* Now, go looking for 'routes' A->B->C->D->A */
1990 {
1991 int i;
1992 route.w = src->w;
1993 route.h = src->h;
1994 route.best_score = -1;
1995 for (i = 0; i < num_points; i++)
1996 {
1997 route.edge[0] = points[i].e0;
1998 route.point[0] = i;
1999 route.edge[1] = points[i].e1;
2000 find_route(ctx, points, num_points, &route, 1);
2001 }
2002
2003 #ifdef WARP_DEBUG
2004 if (route.best_score >= 0)
2005 {
2006 debug_printf(ctx, "Score: %d, Edges=%d->%d->%d->%d, Points=%d->%d->%d->%d\n",
2007 route.best_score,
2008 route.best_edge[0],
2009 route.best_edge[1],
2010 route.best_edge[2],
2011 route.best_edge[3],
2012 route.best_point[0],
2013 route.best_point[1],
2014 route.best_point[2],
2015 route.best_point[3]);
2016 debug_printf(ctx, "(%d,%d)->(%d,%d)->(%d,%d)->(%d,%d)\n",
2017 (int)points[route.best_point[0]].x,
2018 (int)points[route.best_point[0]].y,
2019 (int)points[route.best_point[1]].x,
2020 (int)points[route.best_point[1]].y,
2021 (int)points[route.best_point[2]].x,
2022 (int)points[route.best_point[2]].y,
2023 (int)points[route.best_point[3]].x,
2024 (int)points[route.best_point[3]].y);
2025 }
2026 #endif
2027 }
2028
2029 #ifdef WARP_DEBUG
2030 /* Mark up the src (again, for debugging) */
2031 {
2032 fz_device *dev = fz_new_draw_device(ctx, fz_identity, src);
2033 fz_stroke_state *stroke = fz_new_stroke_state(ctx);
2034 float col = 1;
2035 fz_color_params params = { FZ_RI_PERCEPTUAL };
2036 #ifdef WARP_SPEW_DEBUG
2037 for (x = 0; x < num_edges; x++)
2038 {
2039 char text[64];
2040 fz_path *path = fz_new_path(ctx);
2041 fz_moveto(ctx, path, edge[x].x0, edge[x].y0);
2042 fz_lineto(ctx, path, edge[x].x1, edge[x].y1);
2043 fz_stroke_path(ctx, dev, path, stroke, fz_identity, fz_device_gray(ctx), &col, 1, params);
2044 fz_drop_path(ctx, path);
2045 debug_printf(ctx, "%d %d -> %d %d\n", edge[x].x0, edge[x].y0, edge[x].x1, edge[x].y1);
2046 sprintf(text, "line%d.png", x);
2047 fz_save_pixmap_as_png(ctx, src, text);
2048 }
2049 #endif
2050
2051 stroke->linewidth *= 4;
2052 if (route.best_score >= 0)
2053 {
2054 fz_path *path = fz_new_path(ctx);
2055 fz_moveto(ctx, path, points[route.best_point[0]].x, points[route.best_point[0]].y);
2056 fz_lineto(ctx, path, points[route.best_point[1]].x, points[route.best_point[1]].y);
2057 fz_lineto(ctx, path, points[route.best_point[2]].x, points[route.best_point[2]].y);
2058 fz_lineto(ctx, path, points[route.best_point[3]].x, points[route.best_point[3]].y);
2059 fz_closepath(ctx, path);
2060 fz_stroke_path(ctx, dev, path, stroke, fz_identity, fz_device_gray(ctx), &col, 1, params);
2061 fz_drop_path(ctx, path);
2062 }
2063
2064 fz_drop_stroke_state(ctx, stroke);
2065 fz_close_device(ctx, dev);
2066 fz_drop_device(ctx, dev);
2067 }
2068 #endif
2069
2070 fz_free(ctx, hough);
2071
2072 if (route.best_score == -1)
2073 return 0;
2074
2075 if (corners)
2076 {
2077 corners->ul.x = points[route.best_point[0]].x;
2078 corners->ul.y = points[route.best_point[0]].y;
2079 corners->ur.x = points[route.best_point[1]].x;
2080 corners->ur.y = points[route.best_point[1]].y;
2081 corners->ll.x = points[route.best_point[2]].x;
2082 corners->ll.y = points[route.best_point[2]].y;
2083 corners->lr.x = points[route.best_point[3]].x;
2084 corners->lr.y = points[route.best_point[3]].y;
2085 }
2086
2087 /* Discard any possible matches that aren't at least 1/8 of the pixmap. */
2088 {
2089 fz_rect r = fz_empty_rect;
2090 if (corners)
2091 {
2092 r = fz_include_point_in_rect(r, corners->ul);
2093 r = fz_include_point_in_rect(r, corners->ur);
2094 r = fz_include_point_in_rect(r, corners->ll);
2095 r = fz_include_point_in_rect(r, corners->lr);
2096 }
2097 if ((r.x1 - r.x0) * (r.y1 - r.y0) * 8 < (src->w * src->h))
2098 return 0;
2099 }
2100
2101 return 1;
2102 }
2103
2104 #define DOC_DETECT_MAXDIM 500
2105 int
2106 fz_detect_document(fz_context *ctx, fz_quad *points, fz_pixmap *orig_src)
2107 {
2108 fz_color_params p = {FZ_RI_PERCEPTUAL };
2109 fz_pixmap *grey = NULL;
2110 fz_pixmap *src = NULL;
2111 #ifdef DETECT_DOCUMENT_RGB
2112 fz_pixmap *r = NULL;
2113 fz_pixmap *g = NULL;
2114 fz_pixmap *b = NULL;
2115 #endif
2116 fz_pixmap *processed = NULL;
2117 int i;
2118 int found = 0;
2119
2120 /* Gauss function has std deviation of ~sqr(2), so a variance of
2121 * 2. Apply it twice and we get twice that. So:
2122 * n stddev variance
2123 * 1 1.4 2
2124 * 2 2 4
2125 * 3 2.8 8
2126 * 4 4 16
2127 * 5 5.6 32
2128 * 6 8 64
2129 * 7 11.3 128
2130 * 8 16 256
2131 * 9 22.6 512
2132 * 10 32 1024
2133 * 11 45 2048
2134 *
2135 * The stddev is also known as the "width" by some sources.
2136 *
2137 * Our testing indicates that a width of 30ish works well for a
2138 * image with it's major axis being ~3k pixels. So figure on a
2139 * width of about 1% of the major axis dimension.
2140 *
2141 * By subsampling the incoming image, we can reduce the work we
2142 * do in the gauss phase, as we get to use a smaller width.
2143 */
2144 int n = 10; /* Based on DOC_DETECT_MAXDIM */
2145 int l2factor = 0;
2146
2147 fz_var(src);
2148 fz_var(grey);
2149 #ifdef DETECT_DOCUMENT_RGB
2150 fz_var(r);
2151 fz_var(g);
2152 fz_var(b);
2153 #endif
2154 fz_var(processed);
2155
2156 fz_try(ctx)
2157 {
2158 START_TIME();
2159 {
2160 int maxdim = orig_src->w > orig_src->h ? orig_src->w : orig_src->h;
2161 while (maxdim > DOC_DETECT_MAXDIM)
2162 maxdim >>= 1, l2factor++;
2163 START_TIME();
2164 if (l2factor == 0)
2165 src = fz_keep_pixmap(ctx, orig_src);
2166 else
2167 {
2168 src = fz_clone_pixmap(ctx, orig_src);
2169 fz_subsample_pixmap(ctx, src, l2factor);
2170 }
2171 END_TIME("subsample");
2172 }
2173 START_TIME();
2174 if (src->n == 1 && src->alpha == 0)
2175 {
2176 if (src == orig_src)
2177 grey = fz_clone_pixmap(ctx, src);
2178 else
2179 grey = fz_keep_pixmap(ctx, src);
2180 }
2181 else
2182 {
2183 grey = fz_convert_pixmap(ctx, src,
2184 fz_default_gray(ctx, NULL), NULL, NULL, p, 0);
2185 }
2186 END_TIME("clone");
2187 START_TIME();
2188 for (i = 0; i < n; i++)
2189 gauss5x5(ctx, grey);
2190 END_TIME("gauss grey");
2191 #ifdef WARP_DEBUG
2192 fz_save_pixmap_as_png(ctx, grey, "gauss.png");
2193 #endif
2194 #ifdef DO_HISTEQ
2195 START_TIME();
2196 histeq(grey);
2197 END_TIME("histeq grey");
2198 #ifdef WARP_DEBUG
2199 fz_save_pixmap_as_png(ctx, grey, "hist.png");
2200 #endif
2201 #endif
2202 START_TIME();
2203 grad(ctx, grey, pregrad(ctx, grey));
2204 END_TIME("grad grey");
2205 #ifdef WARP_DEBUG
2206 fz_save_pixmap_as_png(ctx, grey, "grad.png");
2207 #endif
2208
2209 #ifdef DETECT_DOCUMENT_RGB
2210 if (src->n == 3 && src->alpha == 0)
2211 {
2212 START_TIME();
2213 r = fz_new_pixmap(ctx, NULL, src->w, src->h, NULL, 0);
2214 g = fz_new_pixmap(ctx, NULL, src->w, src->h, NULL, 0);
2215 b = fz_new_pixmap(ctx, NULL, src->w, src->h, NULL, 0);
2216 gauss5x5_3(ctx, r, src, 0);
2217 for (i = 1; i < n; i++)
2218 gauss5x5(ctx, r);
2219 gauss5x5_3(ctx, g, src, 1);
2220 for (i = 1; i < n; i++)
2221 gauss5x5(ctx, g);
2222 gauss5x5_3(ctx, b, src, 2);
2223 for (i = 1; i < n; i++)
2224 gauss5x5(ctx, b);
2225 #ifdef DO_HISTEQ
2226 histeq(r);
2227 histeq(g);
2228 histeq(b);
2229 #endif
2230 grad(ctx, r);
2231 grad(ctx, g);
2232 grad(ctx, b);
2233 #ifdef WARP_DEBUG
2234 fz_save_pixmap_as_png(ctx, r, "r.png");
2235 fz_save_pixmap_as_png(ctx, g, "g.png");
2236 fz_save_pixmap_as_png(ctx, b, "b.png");
2237 #endif
2238 combine_grad(grey, r, g, b);
2239 END_TIME("rgb");
2240 fz_drop_pixmap(ctx, r);
2241 fz_drop_pixmap(ctx, g);
2242 fz_drop_pixmap(ctx, b);
2243 r = NULL;
2244 g = NULL;
2245 b = NULL;
2246 }
2247 #ifdef WARP_DEBUG
2248 fz_save_pixmap_as_png(ctx, grey, "combined.png");
2249 #endif
2250 #endif
2251 processed = fz_new_pixmap(ctx, fz_device_gray(ctx), grey->w, grey->h, NULL, 0);
2252
2253 /* Do multiple passes if required, dropping the thresholds for
2254 * strong/weak pixels each time until we find a suitable result. */
2255 for (i = 0; i < 6; i++)
2256 {
2257 START_TIME();
2258 nonmax(ctx, processed, grey, i);
2259 END_TIME("nonmax");
2260 #ifdef WARP_DEBUG
2261 fz_save_pixmap_as_png(ctx, processed, "nonmax.png");
2262 #endif
2263 START_TIME();
2264 hysteresis(ctx, processed);
2265 END_TIME("hysteresis");
2266 #ifdef WARP_DEBUG
2267 fz_save_pixmap_as_png(ctx, processed, "hysteresis.png");
2268 #endif
2269 START_TIME();
2270 found = make_hough(ctx, processed, points);
2271 END_TIME("total hough");
2272 END_TIME("Total time");
2273 DUMP_TIMES();
2274 #ifdef WARP_DEBUG
2275 clean(ctx, processed);
2276 fz_save_pixmap_as_png(ctx, processed, "out.png");
2277 #endif
2278 if (found)
2279 break;
2280 }
2281 }
2282 fz_always(ctx)
2283 {
2284 fz_drop_pixmap(ctx, src);
2285 #ifdef DETECT_DOCUMENT_RGB
2286 fz_drop_pixmap(ctx, r);
2287 fz_drop_pixmap(ctx, g);
2288 fz_drop_pixmap(ctx, b);
2289 #endif
2290 fz_drop_pixmap(ctx, grey);
2291 fz_drop_pixmap(ctx, processed);
2292 }
2293 fz_catch(ctx)
2294 {
2295 fz_rethrow(ctx);
2296 }
2297
2298 if (found && points)
2299 {
2300 float f = (1<<l2factor);
2301 float r = f/2;
2302 points->ul.x = points->ul.x *f + r;
2303 points->ul.y = points->ul.y *f + r;
2304 points->ur.x = points->ur.x *f + r;
2305 points->ur.y = points->ur.y *f + r;
2306 points->ll.x = points->ll.x *f + r;
2307 points->ll.y = points->ll.y *f + r;
2308 points->lr.x = points->lr.x *f + r;
2309 points->lr.y = points->lr.y *f + r;
2310 }
2311
2312 return found;
2313 }