Leptonica 1.83.1
Image processing and image analysis suite
Loading...
Searching...
No Matches
compare.c
Go to the documentation of this file.
1/*====================================================================*
2 - Copyright (C) 2001 Leptonica. All rights reserved.
3 -
4 - Redistribution and use in source and binary forms, with or without
5 - modification, are permitted provided that the following conditions
6 - are met:
7 - 1. Redistributions of source code must retain the above copyright
8 - notice, this list of conditions and the following disclaimer.
9 - 2. Redistributions in binary form must reproduce the above
10 - copyright notice, this list of conditions and the following
11 - disclaimer in the documentation and/or other materials
12 - provided with the distribution.
13 -
14 - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
15 - ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
16 - LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
17 - A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL ANY
18 - CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19 - EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20 - PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21 - PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
22 - OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23 - NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 *====================================================================*/
26
104
105#ifdef HAVE_CONFIG_H
106#include <config_auto.h>
107#endif /* HAVE_CONFIG_H */
108
109#include <string.h>
110#include <math.h>
111#include "allheaders.h"
112
113 /* Small enough to consider equal to 0.0, for plot output */
114static const l_float32 TINY = 0.00001;
115
116static l_ok findHistoGridDimensions(l_int32 n, l_int32 w, l_int32 h,
117 l_int32 *pnx, l_int32 *pny, l_int32 debug);
118static l_ok pixCompareTilesByHisto(PIX *pix1, PIX *pix2, l_int32 maxgray,
119 l_int32 factor, l_int32 n,
120 l_float32 *pscore, PIXA *pixadebug);
121
122/*------------------------------------------------------------------*
123 * Test for pix equality *
124 *------------------------------------------------------------------*/
155l_ok
156pixEqual(PIX *pix1,
157 PIX *pix2,
158 l_int32 *psame)
159{
160 return pixEqualWithAlpha(pix1, pix2, 0, psame);
161}
162
163
181l_ok
182pixEqualWithAlpha(PIX *pix1,
183 PIX *pix2,
184 l_int32 use_alpha,
185 l_int32 *psame)
186{
187l_int32 w1, h1, d1, w2, h2, d2, wpl1, wpl2;
188l_int32 spp1, spp2, i, j, color, mismatch, opaque;
189l_int32 fullwords, linebits, endbits;
190l_uint32 endmask, wordmask;
191l_uint32 *data1, *data2, *line1, *line2;
192PIX *pixs1, *pixs2, *pixt1, *pixt2, *pixalpha;
193PIXCMAP *cmap1, *cmap2;
194
195 if (!psame)
196 return ERROR_INT("psame not defined", __func__, 1);
197 *psame = 0; /* init to not equal */
198 if (!pix1 || !pix2)
199 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
200 pixGetDimensions(pix1, &w1, &h1, &d1);
201 pixGetDimensions(pix2, &w2, &h2, &d2);
202 if (w1 != w2 || h1 != h2) {
203 L_INFO("pix sizes differ\n", __func__);
204 return 0;
205 }
206
207 /* Suppose the use_alpha flag is true.
208 * If only one of two 32 bpp images has spp == 4, we call that
209 * a "mismatch" of the alpha component. In the case of a mismatch,
210 * if the 4 bpp pix does not have all alpha components opaque (255),
211 * the images are not-equal. However if they are all opaque,
212 * this image is equivalent to spp == 3, so we allow the
213 * comparison to go forward, testing only for the RGB equality. */
214 spp1 = pixGetSpp(pix1);
215 spp2 = pixGetSpp(pix2);
216 mismatch = 0;
217 if (use_alpha && d1 == 32 && d2 == 32) {
218 mismatch = ((spp1 == 4 && spp2 != 4) || (spp1 != 4 && spp2 == 4));
219 if (mismatch) {
220 pixalpha = (spp1 == 4) ? pix1 : pix2;
221 pixAlphaIsOpaque(pixalpha, &opaque);
222 if (!opaque) {
223 L_INFO("just one pix has a non-opaque alpha layer\n", __func__);
224 return 0;
225 }
226 }
227 }
228
229 cmap1 = pixGetColormap(pix1);
230 cmap2 = pixGetColormap(pix2);
231 if (!cmap1 && !cmap2 && (d1 != d2) && (d1 == 32 || d2 == 32)) {
232 L_INFO("no colormaps, pix depths unequal, and one of them is RGB\n",
233 __func__);
234 return 0;
235 }
236
237 if (cmap1 && cmap2 && (d1 == d2)) /* use special function */
238 return pixEqualWithCmap(pix1, pix2, psame);
239
240 /* Must remove colormaps if they exist, and in the process
241 * end up with the resulting images having the same depth. */
242 if (cmap1 && !cmap2) {
243 pixUsesCmapColor(pix1, &color);
244 if (color && d2 <= 8) /* can't be equal */
245 return 0;
246 if (d2 < 8)
247 pixs2 = pixConvertTo8(pix2, FALSE);
248 else
249 pixs2 = pixClone(pix2);
250 if (d2 <= 8)
251 pixs1 = pixRemoveColormap(pix1, REMOVE_CMAP_TO_GRAYSCALE);
252 else
253 pixs1 = pixRemoveColormap(pix1, REMOVE_CMAP_TO_FULL_COLOR);
254 } else if (!cmap1 && cmap2) {
255 pixUsesCmapColor(pix2, &color);
256 if (color && d1 <= 8) /* can't be equal */
257 return 0;
258 if (d1 < 8)
259 pixs1 = pixConvertTo8(pix1, FALSE);
260 else
261 pixs1 = pixClone(pix1);
262 if (d1 <= 8)
263 pixs2 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_GRAYSCALE);
264 else
265 pixs2 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_FULL_COLOR);
266 } else if (cmap1 && cmap2) { /* depths not equal; use rgb */
267 pixs1 = pixRemoveColormap(pix1, REMOVE_CMAP_TO_FULL_COLOR);
268 pixs2 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_FULL_COLOR);
269 } else { /* no colormaps */
270 pixs1 = pixClone(pix1);
271 pixs2 = pixClone(pix2);
272 }
273
274 /* OK, we have no colormaps, but the depths may still be different */
275 d1 = pixGetDepth(pixs1);
276 d2 = pixGetDepth(pixs2);
277 if (d1 != d2) {
278 if (d1 == 16 || d2 == 16) {
279 L_INFO("one pix is 16 bpp\n", __func__);
280 pixDestroy(&pixs1);
281 pixDestroy(&pixs2);
282 return 0;
283 }
284 pixt1 = pixConvertLossless(pixs1, 8);
285 pixt2 = pixConvertLossless(pixs2, 8);
286 if (!pixt1 || !pixt2) {
287 L_INFO("failure to convert to 8 bpp\n", __func__);
288 pixDestroy(&pixs1);
289 pixDestroy(&pixs2);
290 pixDestroy(&pixt1);
291 pixDestroy(&pixt2);
292 return 0;
293 }
294 } else {
295 pixt1 = pixClone(pixs1);
296 pixt2 = pixClone(pixs2);
297 }
298 pixDestroy(&pixs1);
299 pixDestroy(&pixs2);
300
301 /* No colormaps, equal depths; do pixel comparisons */
302 d1 = pixGetDepth(pixt1);
303 d2 = pixGetDepth(pixt2);
304 wpl1 = pixGetWpl(pixt1);
305 wpl2 = pixGetWpl(pixt2);
306 data1 = pixGetData(pixt1);
307 data2 = pixGetData(pixt2);
308
309 if (d1 == 32) { /* test either RGB or RGBA pixels */
310 if (use_alpha && !mismatch)
311 wordmask = (spp1 == 3) ? 0xffffff00 : 0xffffffff;
312 else
313 wordmask = 0xffffff00;
314 for (i = 0; i < h1; i++) {
315 line1 = data1 + wpl1 * i;
316 line2 = data2 + wpl2 * i;
317 for (j = 0; j < wpl1; j++) {
318 if ((*line1 ^ *line2) & wordmask) {
319 pixDestroy(&pixt1);
320 pixDestroy(&pixt2);
321 return 0;
322 }
323 line1++;
324 line2++;
325 }
326 }
327 } else { /* all bits count */
328 linebits = d1 * w1;
329 fullwords = linebits / 32;
330 endbits = linebits & 31;
331 endmask = (endbits == 0) ? 0 : (0xffffffff << (32 - endbits));
332 for (i = 0; i < h1; i++) {
333 line1 = data1 + wpl1 * i;
334 line2 = data2 + wpl2 * i;
335 for (j = 0; j < fullwords; j++) {
336 if (*line1 ^ *line2) {
337 pixDestroy(&pixt1);
338 pixDestroy(&pixt2);
339 return 0;
340 }
341 line1++;
342 line2++;
343 }
344 if (endbits) {
345 if ((*line1 ^ *line2) & endmask) {
346 pixDestroy(&pixt1);
347 pixDestroy(&pixt2);
348 return 0;
349 }
350 }
351 }
352 }
353
354 pixDestroy(&pixt1);
355 pixDestroy(&pixt2);
356 *psame = 1;
357 return 0;
358}
359
360
381l_ok
382pixEqualWithCmap(PIX *pix1,
383 PIX *pix2,
384 l_int32 *psame)
385{
386l_int32 d, w, h, wpl1, wpl2, i, j, linebits, fullwords, endbits;
387l_int32 rval1, rval2, gval1, gval2, bval1, bval2, samecmaps;
388l_uint32 endmask, val1, val2;
389l_uint32 *data1, *data2, *line1, *line2;
390PIXCMAP *cmap1, *cmap2;
391
392 if (!psame)
393 return ERROR_INT("&same not defined", __func__, 1);
394 *psame = 0;
395 if (!pix1)
396 return ERROR_INT("pix1 not defined", __func__, 1);
397 if (!pix2)
398 return ERROR_INT("pix2 not defined", __func__, 1);
399
400 if (pixSizesEqual(pix1, pix2) == 0)
401 return 0;
402 cmap1 = pixGetColormap(pix1);
403 cmap2 = pixGetColormap(pix2);
404 if (!cmap1 || !cmap2) {
405 L_INFO("both images don't have colormap\n", __func__);
406 return 0;
407 }
408 pixGetDimensions(pix1, &w, &h, &d);
409 if (d != 1 && d != 2 && d != 4 && d != 8) {
410 L_INFO("pix depth not in {1, 2, 4, 8}\n", __func__);
411 return 0;
412 }
413
414 cmapEqual(cmap1, cmap2, 3, &samecmaps);
415 if (samecmaps == TRUE) { /* colormaps are identical; compare by words */
416 linebits = d * w;
417 wpl1 = pixGetWpl(pix1);
418 wpl2 = pixGetWpl(pix2);
419 data1 = pixGetData(pix1);
420 data2 = pixGetData(pix2);
421 fullwords = linebits / 32;
422 endbits = linebits & 31;
423 endmask = (endbits == 0) ? 0 : (0xffffffff << (32 - endbits));
424 for (i = 0; i < h; i++) {
425 line1 = data1 + wpl1 * i;
426 line2 = data2 + wpl2 * i;
427 for (j = 0; j < fullwords; j++) {
428 if (*line1 ^ *line2)
429 return 0;
430 line1++;
431 line2++;
432 }
433 if (endbits) {
434 if ((*line1 ^ *line2) & endmask)
435 return 0;
436 }
437 }
438 *psame = 1;
439 return 0;
440 }
441
442 /* Colormaps aren't identical; compare pixel by pixel */
443 for (i = 0; i < h; i++) {
444 for (j = 0; j < w; j++) {
445 pixGetPixel(pix1, j, i, &val1);
446 pixGetPixel(pix2, j, i, &val2);
447 pixcmapGetColor(cmap1, val1, &rval1, &gval1, &bval1);
448 pixcmapGetColor(cmap2, val2, &rval2, &gval2, &bval2);
449 if (rval1 != rval2 || gval1 != gval2 || bval1 != bval2)
450 return 0;
451 }
452 }
453
454 *psame = 1;
455 return 0;
456}
457
458
475l_ok
476cmapEqual(PIXCMAP *cmap1,
477 PIXCMAP *cmap2,
478 l_int32 ncomps,
479 l_int32 *psame)
480{
481l_int32 n1, n2, i, rval1, rval2, gval1, gval2, bval1, bval2, aval1, aval2;
482
483 if (!psame)
484 return ERROR_INT("&same not defined", __func__, 1);
485 *psame = FALSE;
486 if (!cmap1)
487 return ERROR_INT("cmap1 not defined", __func__, 1);
488 if (!cmap2)
489 return ERROR_INT("cmap2 not defined", __func__, 1);
490 if (ncomps != 3 && ncomps != 4)
491 return ERROR_INT("ncomps not 3 or 4", __func__, 1);
492
493 n1 = pixcmapGetCount(cmap1);
494 n2 = pixcmapGetCount(cmap2);
495 if (n1 != n2) {
496 L_INFO("colormap sizes are different\n", __func__);
497 return 0;
498 }
499
500 for (i = 0; i < n1; i++) {
501 pixcmapGetRGBA(cmap1, i, &rval1, &gval1, &bval1, &aval1);
502 pixcmapGetRGBA(cmap2, i, &rval2, &gval2, &bval2, &aval2);
503 if (rval1 != rval2 || gval1 != gval2 || bval1 != bval2)
504 return 0;
505 if (ncomps == 4 && aval1 != aval2)
506 return 0;
507 }
508 *psame = TRUE;
509 return 0;
510}
511
512
531l_ok
532pixUsesCmapColor(PIX *pixs,
533 l_int32 *pcolor)
534{
535l_int32 n, i, rval, gval, bval, numpix;
536NUMA *na;
537PIXCMAP *cmap;
538
539 if (!pcolor)
540 return ERROR_INT("&color not defined", __func__, 1);
541 *pcolor = 0;
542 if (!pixs)
543 return ERROR_INT("pixs not defined", __func__, 1);
544
545 if ((cmap = pixGetColormap(pixs)) == NULL)
546 return 0;
547
548 pixcmapHasColor(cmap, pcolor);
549 if (*pcolor == 0) /* no color */
550 return 0;
551
552 /* The cmap has color entries. Are they used? */
553 na = pixGetGrayHistogram(pixs, 1);
554 n = pixcmapGetCount(cmap);
555 for (i = 0; i < n; i++) {
556 pixcmapGetColor(cmap, i, &rval, &gval, &bval);
557 numaGetIValue(na, i, &numpix);
558 if ((rval != gval || rval != bval) && numpix) { /* color found! */
559 *pcolor = 1;
560 break;
561 }
562 }
563 numaDestroy(&na);
564
565 return 0;
566}
567
568
569/*------------------------------------------------------------------*
570 * Binary correlation *
571 *------------------------------------------------------------------*/
595l_ok
596pixCorrelationBinary(PIX *pix1,
597 PIX *pix2,
598 l_float32 *pval)
599{
600l_int32 count1, count2, countn;
601l_int32 *tab8;
602PIX *pixn;
603
604 if (!pval)
605 return ERROR_INT("&pval not defined", __func__, 1);
606 *pval = 0.0;
607 if (!pix1)
608 return ERROR_INT("pix1 not defined", __func__, 1);
609 if (!pix2)
610 return ERROR_INT("pix2 not defined", __func__, 1);
611
612 tab8 = makePixelSumTab8();
613 pixCountPixels(pix1, &count1, tab8);
614 pixCountPixels(pix2, &count2, tab8);
615 if (count1 == 0 || count2 == 0) {
616 LEPT_FREE(tab8);
617 return 0;
618 }
619 pixn = pixAnd(NULL, pix1, pix2);
620 pixCountPixels(pixn, &countn, tab8);
621 *pval = (l_float32)countn * (l_float32)countn /
622 ((l_float32)count1 * (l_float32)count2);
623 LEPT_FREE(tab8);
624 pixDestroy(&pixn);
625 return 0;
626}
627
628
629/*------------------------------------------------------------------*
630 * Difference of two images *
631 *------------------------------------------------------------------*/
651PIX *
652pixDisplayDiffBinary(PIX *pix1,
653 PIX *pix2)
654{
655l_int32 w1, h1, d1, w2, h2, d2, minw, minh;
656PIX *pixt, *pixd;
657PIXCMAP *cmap;
658
659 if (!pix1 || !pix2)
660 return (PIX *)ERROR_PTR("pix1, pix2 not both defined", __func__, NULL);
661 pixGetDimensions(pix1, &w1, &h1, &d1);
662 pixGetDimensions(pix2, &w2, &h2, &d2);
663 if (d1 != 1 || d2 != 1)
664 return (PIX *)ERROR_PTR("pix1 and pix2 not 1 bpp", __func__, NULL);
665 minw = L_MIN(w1, w2);
666 minh = L_MIN(h1, h2);
667
668 pixd = pixCreate(minw, minh, 4);
669 cmap = pixcmapCreate(4);
670 pixcmapAddColor(cmap, 255, 255, 255); /* initialized to white */
671 pixcmapAddColor(cmap, 0, 0, 0);
672 pixcmapAddColor(cmap, 255, 0, 0);
673 pixcmapAddColor(cmap, 0, 255, 0);
674 pixSetColormap(pixd, cmap);
675
676 pixt = pixAnd(NULL, pix1, pix2);
677 pixPaintThroughMask(pixd, pixt, 0, 0, 0x0); /* black */
678 pixSubtract(pixt, pix1, pix2);
679 pixPaintThroughMask(pixd, pixt, 0, 0, 0xff000000); /* red */
680 pixSubtract(pixt, pix2, pix1);
681 pixPaintThroughMask(pixd, pixt, 0, 0, 0x00ff0000); /* green */
682 pixDestroy(&pixt);
683 return pixd;
684}
685
686
706l_ok
707pixCompareBinary(PIX *pix1,
708 PIX *pix2,
709 l_int32 comptype,
710 l_float32 *pfract,
711 PIX **ppixdiff)
712{
713l_int32 w, h, count;
714PIX *pixt;
715
716 if (ppixdiff) *ppixdiff = NULL;
717 if (!pfract)
718 return ERROR_INT("&pfract not defined", __func__, 1);
719 *pfract = 1.0; /* initialize to max difference */
720 if (!pix1 || pixGetDepth(pix1) != 1)
721 return ERROR_INT("pix1 not defined or not 1 bpp", __func__, 1);
722 if (!pix2 || pixGetDepth(pix2) != 1)
723 return ERROR_INT("pix2 not defined or not 1 bpp", __func__, 1);
724 if (comptype != L_COMPARE_XOR && comptype != L_COMPARE_SUBTRACT)
725 return ERROR_INT("invalid comptype", __func__, 1);
726
727 if (comptype == L_COMPARE_XOR)
728 pixt = pixXor(NULL, pix1, pix2);
729 else /* comptype == L_COMPARE_SUBTRACT) */
730 pixt = pixSubtract(NULL, pix1, pix2);
731 pixCountPixels(pixt, &count, NULL);
732 pixGetDimensions(pix1, &w, &h, NULL);
733 *pfract = (l_float32)(count) / (l_float32)(w * h);
734
735 if (ppixdiff)
736 *ppixdiff = pixt;
737 else
738 pixDestroy(&pixt);
739 return 0;
740}
741
742
784l_ok
785pixCompareGrayOrRGB(PIX *pix1,
786 PIX *pix2,
787 l_int32 comptype,
788 l_int32 plottype,
789 l_int32 *psame,
790 l_float32 *pdiff,
791 l_float32 *prmsdiff,
792 PIX **ppixdiff)
793{
794l_int32 retval, d1, d2;
795PIX *pixt1, *pixt2, *pixs1, *pixs2;
796
797 if (psame) *psame = 0;
798 if (pdiff) *pdiff = 255.0;
799 if (prmsdiff) *prmsdiff = 255.0;
800 if (ppixdiff) *ppixdiff = NULL;
801 if (!pix1 || pixGetDepth(pix1) == 1)
802 return ERROR_INT("pix1 not defined or 1 bpp", __func__, 1);
803 if (!pix2 || pixGetDepth(pix2) == 1)
804 return ERROR_INT("pix2 not defined or 1 bpp", __func__, 1);
805 if (comptype != L_COMPARE_SUBTRACT && comptype != L_COMPARE_ABS_DIFF)
806 return ERROR_INT("invalid comptype", __func__, 1);
807 if (plottype < 0 || plottype >= NUM_GPLOT_OUTPUTS)
808 return ERROR_INT("invalid plottype", __func__, 1);
809
810 pixt1 = pixRemoveColormap(pix1, REMOVE_CMAP_BASED_ON_SRC);
811 pixt2 = pixRemoveColormap(pix2, REMOVE_CMAP_BASED_ON_SRC);
812 d1 = pixGetDepth(pixt1);
813 d2 = pixGetDepth(pixt2);
814 if (d1 < 8)
815 pixs1 = pixConvertTo8(pixt1, FALSE);
816 else
817 pixs1 = pixClone(pixt1);
818 if (d2 < 8)
819 pixs2 = pixConvertTo8(pixt2, FALSE);
820 else
821 pixs2 = pixClone(pixt2);
822 pixDestroy(&pixt1);
823 pixDestroy(&pixt2);
824 d1 = pixGetDepth(pixs1);
825 d2 = pixGetDepth(pixs2);
826 if (d1 != d2) {
827 pixDestroy(&pixs1);
828 pixDestroy(&pixs2);
829 return ERROR_INT("intrinsic depths are not equal", __func__, 1);
830 }
831
832 if (d1 == 8 || d1 == 16)
833 retval = pixCompareGray(pixs1, pixs2, comptype, plottype, psame,
834 pdiff, prmsdiff, ppixdiff);
835 else /* d1 == 32 */
836 retval = pixCompareRGB(pixs1, pixs2, comptype, plottype, psame,
837 pdiff, prmsdiff, ppixdiff);
838 pixDestroy(&pixs1);
839 pixDestroy(&pixs2);
840 return retval;
841}
842
843
865l_ok
866pixCompareGray(PIX *pix1,
867 PIX *pix2,
868 l_int32 comptype,
869 l_int32 plottype,
870 l_int32 *psame,
871 l_float32 *pdiff,
872 l_float32 *prmsdiff,
873 PIX **ppixdiff)
874{
875char buf[64];
876static l_int32 index = 0;
877l_int32 d1, d2, same, first, last;
878GPLOT *gplot;
879NUMA *na, *nac;
880PIX *pixt;
881
882 if (psame) *psame = 0;
883 if (pdiff) *pdiff = 255.0;
884 if (prmsdiff) *prmsdiff = 255.0;
885 if (ppixdiff) *ppixdiff = NULL;
886 if (!pix1)
887 return ERROR_INT("pix1 not defined", __func__, 1);
888 if (!pix2)
889 return ERROR_INT("pix2 not defined", __func__, 1);
890 d1 = pixGetDepth(pix1);
891 d2 = pixGetDepth(pix2);
892 if ((d1 != d2) || (d1 != 8 && d1 != 16))
893 return ERROR_INT("depths unequal or not 8 or 16 bpp", __func__, 1);
894 if (pixGetColormap(pix1) || pixGetColormap(pix2))
895 return ERROR_INT("pix1 and/or pix2 are colormapped", __func__, 1);
896 if (comptype != L_COMPARE_SUBTRACT && comptype != L_COMPARE_ABS_DIFF)
897 return ERROR_INT("invalid comptype", __func__, 1);
898 if (plottype < 0 || plottype >= NUM_GPLOT_OUTPUTS)
899 return ERROR_INT("invalid plottype", __func__, 1);
900
901 lept_mkdir("lept/comp");
902
903 if (comptype == L_COMPARE_SUBTRACT)
904 pixt = pixSubtractGray(NULL, pix1, pix2);
905 else /* comptype == L_COMPARE_ABS_DIFF) */
906 pixt = pixAbsDifference(pix1, pix2);
907
908 pixZero(pixt, &same);
909 if (same)
910 L_INFO("Images are pixel-wise identical\n", __func__);
911 if (psame) *psame = same;
912
913 if (pdiff)
914 pixGetAverageMasked(pixt, NULL, 0, 0, 1, L_MEAN_ABSVAL, pdiff);
915
916 /* Don't bother to plot if the images are the same */
917 if (plottype && !same) {
918 L_INFO("Images differ: output plots will be generated\n", __func__);
919 na = pixGetGrayHistogram(pixt, 1);
920 numaGetNonzeroRange(na, TINY, &first, &last);
921 nac = numaClipToInterval(na, 0, last);
922 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_gray%d", index);
923 gplot = gplotCreate(buf, plottype,
924 "Pixel Difference Histogram", "diff val",
925 "number of pixels");
926 gplotAddPlot(gplot, NULL, nac, GPLOT_LINES, "gray");
927 gplotMakeOutput(gplot);
928 gplotDestroy(&gplot);
929 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_gray%d.png",
930 index++);
931 l_fileDisplay(buf, 100, 100, 1.0);
932 numaDestroy(&na);
933 numaDestroy(&nac);
934 }
935
936 if (ppixdiff)
937 *ppixdiff = pixCopy(NULL, pixt);
938
939 if (prmsdiff) {
940 if (comptype == L_COMPARE_SUBTRACT) { /* wrong type for rms diff */
941 pixDestroy(&pixt);
942 pixt = pixAbsDifference(pix1, pix2);
943 }
944 pixGetAverageMasked(pixt, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, prmsdiff);
945 }
946
947 pixDestroy(&pixt);
948 return 0;
949}
950
951
972l_ok
973pixCompareRGB(PIX *pix1,
974 PIX *pix2,
975 l_int32 comptype,
976 l_int32 plottype,
977 l_int32 *psame,
978 l_float32 *pdiff,
979 l_float32 *prmsdiff,
980 PIX **ppixdiff)
981{
982char buf[64];
983static l_int32 index = 0;
984l_int32 rsame, gsame, bsame, same, first, rlast, glast, blast, last;
985l_float32 rdiff, gdiff, bdiff;
986GPLOT *gplot;
987NUMA *nar, *nag, *nab, *narc, *nagc, *nabc;
988PIX *pixr1, *pixr2, *pixg1, *pixg2, *pixb1, *pixb2;
989PIX *pixr, *pixg, *pixb;
990
991 if (psame) *psame = 0;
992 if (pdiff) *pdiff = 0.0;
993 if (prmsdiff) *prmsdiff = 0.0;
994 if (ppixdiff) *ppixdiff = NULL;
995 if (!pix1 || pixGetDepth(pix1) != 32)
996 return ERROR_INT("pix1 not defined or not 32 bpp", __func__, 1);
997 if (!pix2 || pixGetDepth(pix2) != 32)
998 return ERROR_INT("pix2 not defined or not ew bpp", __func__, 1);
999 if (comptype != L_COMPARE_SUBTRACT && comptype != L_COMPARE_ABS_DIFF)
1000 return ERROR_INT("invalid comptype", __func__, 1);
1001 if (plottype < 0 || plottype >= NUM_GPLOT_OUTPUTS)
1002 return ERROR_INT("invalid plottype", __func__, 1);
1003
1004 lept_mkdir("lept/comp");
1005
1006 pixr1 = pixGetRGBComponent(pix1, COLOR_RED);
1007 pixr2 = pixGetRGBComponent(pix2, COLOR_RED);
1008 pixg1 = pixGetRGBComponent(pix1, COLOR_GREEN);
1009 pixg2 = pixGetRGBComponent(pix2, COLOR_GREEN);
1010 pixb1 = pixGetRGBComponent(pix1, COLOR_BLUE);
1011 pixb2 = pixGetRGBComponent(pix2, COLOR_BLUE);
1012 if (comptype == L_COMPARE_SUBTRACT) {
1013 pixr = pixSubtractGray(NULL, pixr1, pixr2);
1014 pixg = pixSubtractGray(NULL, pixg1, pixg2);
1015 pixb = pixSubtractGray(NULL, pixb1, pixb2);
1016 } else { /* comptype == L_COMPARE_ABS_DIFF) */
1017 pixr = pixAbsDifference(pixr1, pixr2);
1018 pixg = pixAbsDifference(pixg1, pixg2);
1019 pixb = pixAbsDifference(pixb1, pixb2);
1020 }
1021
1022 pixZero(pixr, &rsame);
1023 pixZero(pixg, &gsame);
1024 pixZero(pixb, &bsame);
1025 same = rsame && gsame && bsame;
1026 if (same)
1027 L_INFO("Images are pixel-wise identical\n", __func__);
1028 if (psame) *psame = same;
1029
1030 if (pdiff) {
1031 pixGetAverageMasked(pixr, NULL, 0, 0, 1, L_MEAN_ABSVAL, &rdiff);
1032 pixGetAverageMasked(pixg, NULL, 0, 0, 1, L_MEAN_ABSVAL, &gdiff);
1033 pixGetAverageMasked(pixb, NULL, 0, 0, 1, L_MEAN_ABSVAL, &bdiff);
1034 *pdiff = (rdiff + gdiff + bdiff) / 3.0;
1035 }
1036
1037 /* Don't bother to plot if the images are the same */
1038 if (plottype && !same) {
1039 L_INFO("Images differ: output plots will be generated\n", __func__);
1040 nar = pixGetGrayHistogram(pixr, 1);
1041 nag = pixGetGrayHistogram(pixg, 1);
1042 nab = pixGetGrayHistogram(pixb, 1);
1043 numaGetNonzeroRange(nar, TINY, &first, &rlast);
1044 numaGetNonzeroRange(nag, TINY, &first, &glast);
1045 numaGetNonzeroRange(nab, TINY, &first, &blast);
1046 last = L_MAX(rlast, glast);
1047 last = L_MAX(last, blast);
1048 narc = numaClipToInterval(nar, 0, last);
1049 nagc = numaClipToInterval(nag, 0, last);
1050 nabc = numaClipToInterval(nab, 0, last);
1051 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_rgb%d", index);
1052 gplot = gplotCreate(buf, plottype,
1053 "Pixel Difference Histogram", "diff val",
1054 "number of pixels");
1055 gplotAddPlot(gplot, NULL, narc, GPLOT_LINES, "red");
1056 gplotAddPlot(gplot, NULL, nagc, GPLOT_LINES, "green");
1057 gplotAddPlot(gplot, NULL, nabc, GPLOT_LINES, "blue");
1058 gplotMakeOutput(gplot);
1059 gplotDestroy(&gplot);
1060 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_rgb%d.png",
1061 index++);
1062 l_fileDisplay(buf, 100, 100, 1.0);
1063 numaDestroy(&nar);
1064 numaDestroy(&nag);
1065 numaDestroy(&nab);
1066 numaDestroy(&narc);
1067 numaDestroy(&nagc);
1068 numaDestroy(&nabc);
1069 }
1070
1071 if (ppixdiff)
1072 *ppixdiff = pixCreateRGBImage(pixr, pixg, pixb);
1073
1074 if (prmsdiff) {
1075 if (comptype == L_COMPARE_SUBTRACT) {
1076 pixDestroy(&pixr);
1077 pixDestroy(&pixg);
1078 pixDestroy(&pixb);
1079 pixr = pixAbsDifference(pixr1, pixr2);
1080 pixg = pixAbsDifference(pixg1, pixg2);
1081 pixb = pixAbsDifference(pixb1, pixb2);
1082 }
1083 pixGetAverageMasked(pixr, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, &rdiff);
1084 pixGetAverageMasked(pixg, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, &gdiff);
1085 pixGetAverageMasked(pixb, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, &bdiff);
1086 *prmsdiff = (rdiff + gdiff + bdiff) / 3.0;
1087 }
1088
1089 pixDestroy(&pixr1);
1090 pixDestroy(&pixr2);
1091 pixDestroy(&pixg1);
1092 pixDestroy(&pixg2);
1093 pixDestroy(&pixb1);
1094 pixDestroy(&pixb2);
1095 pixDestroy(&pixr);
1096 pixDestroy(&pixg);
1097 pixDestroy(&pixb);
1098 return 0;
1099}
1100
1101
1126l_ok
1127pixCompareTiled(PIX *pix1,
1128 PIX *pix2,
1129 l_int32 sx,
1130 l_int32 sy,
1131 l_int32 type,
1132 PIX **ppixdiff)
1133{
1134l_int32 d1, d2, w, h;
1135PIX *pixt, *pixr, *pixg, *pixb;
1136PIX *pixrdiff, *pixgdiff, *pixbdiff;
1137PIXACC *pixacc;
1138
1139 if (!ppixdiff)
1140 return ERROR_INT("&pixdiff not defined", __func__, 1);
1141 *ppixdiff = NULL;
1142 if (!pix1)
1143 return ERROR_INT("pix1 not defined", __func__, 1);
1144 if (!pix2)
1145 return ERROR_INT("pix2 not defined", __func__, 1);
1146 d1 = pixGetDepth(pix1);
1147 d2 = pixGetDepth(pix2);
1148 if (d1 != d2)
1149 return ERROR_INT("depths not equal", __func__, 1);
1150 if (d1 != 8 && d1 != 32)
1151 return ERROR_INT("pix1 not 8 or 32 bpp", __func__, 1);
1152 if (d2 != 8 && d2 != 32)
1153 return ERROR_INT("pix2 not 8 or 32 bpp", __func__, 1);
1154 if (sx < 2 || sy < 2)
1155 return ERROR_INT("sx and sy not both > 1", __func__, 1);
1156 if (type != L_MEAN_ABSVAL && type != L_ROOT_MEAN_SQUARE)
1157 return ERROR_INT("invalid type", __func__, 1);
1158
1159 pixt = pixAbsDifference(pix1, pix2);
1160 if (d1 == 8) {
1161 *ppixdiff = pixGetAverageTiled(pixt, sx, sy, type);
1162 } else { /* d1 == 32 */
1163 pixr = pixGetRGBComponent(pixt, COLOR_RED);
1164 pixg = pixGetRGBComponent(pixt, COLOR_GREEN);
1165 pixb = pixGetRGBComponent(pixt, COLOR_BLUE);
1166 pixrdiff = pixGetAverageTiled(pixr, sx, sy, type);
1167 pixgdiff = pixGetAverageTiled(pixg, sx, sy, type);
1168 pixbdiff = pixGetAverageTiled(pixb, sx, sy, type);
1169 pixGetDimensions(pixrdiff, &w, &h, NULL);
1170 pixacc = pixaccCreate(w, h, 0);
1171 pixaccAdd(pixacc, pixrdiff);
1172 pixaccAdd(pixacc, pixgdiff);
1173 pixaccAdd(pixacc, pixbdiff);
1174 pixaccMultConst(pixacc, 1. / 3.);
1175 *ppixdiff = pixaccFinal(pixacc, 8);
1176 pixDestroy(&pixr);
1177 pixDestroy(&pixg);
1178 pixDestroy(&pixb);
1179 pixDestroy(&pixrdiff);
1180 pixDestroy(&pixgdiff);
1181 pixDestroy(&pixbdiff);
1182 pixaccDestroy(&pixacc);
1183 }
1184 pixDestroy(&pixt);
1185 return 0;
1186}
1187
1188
1189/*------------------------------------------------------------------*
1190 * Other measures of the difference of two images *
1191 *------------------------------------------------------------------*/
1218NUMA *
1219pixCompareRankDifference(PIX *pix1,
1220 PIX *pix2,
1221 l_int32 factor)
1222{
1223l_int32 i;
1224l_float32 *array1, *array2;
1225NUMA *nah, *nan, *nad;
1226
1227 if (!pix1)
1228 return (NUMA *)ERROR_PTR("pix1 not defined", __func__, NULL);
1229 if (!pix2)
1230 return (NUMA *)ERROR_PTR("pix2 not defined", __func__, NULL);
1231
1232 if ((nah = pixGetDifferenceHistogram(pix1, pix2, factor)) == NULL)
1233 return (NUMA *)ERROR_PTR("na not made", __func__, NULL);
1234
1235 nan = numaNormalizeHistogram(nah, 1.0);
1236 array1 = numaGetFArray(nan, L_NOCOPY);
1237
1238 nad = numaCreate(256);
1239 numaSetCount(nad, 256); /* all initialized to 0.0 */
1240 array2 = numaGetFArray(nad, L_NOCOPY);
1241
1242 /* Do rank accumulation on normalized histo of diffs */
1243 array2[0] = 1.0;
1244 for (i = 1; i < 256; i++)
1245 array2[i] = array2[i - 1] - array1[i - 1];
1246
1247 numaDestroy(&nah);
1248 numaDestroy(&nan);
1249 return nad;
1250}
1251
1252
1301l_ok
1302pixTestForSimilarity(PIX *pix1,
1303 PIX *pix2,
1304 l_int32 factor,
1305 l_int32 mindiff,
1306 l_float32 maxfract,
1307 l_float32 maxave,
1308 l_int32 *psimilar,
1309 l_int32 details)
1310{
1311l_float32 fractdiff, avediff;
1312
1313 if (!psimilar)
1314 return ERROR_INT("&similar not defined", __func__, 1);
1315 *psimilar = 0;
1316 if (!pix1)
1317 return ERROR_INT("pix1 not defined", __func__, 1);
1318 if (!pix2)
1319 return ERROR_INT("pix2 not defined", __func__, 1);
1320 if (pixSizesEqual(pix1, pix2) == 0)
1321 return ERROR_INT("pix sizes not equal", __func__, 1);
1322 if (mindiff <= 0)
1323 return ERROR_INT("mindiff must be > 0", __func__, 1);
1324
1325 if (pixGetDifferenceStats(pix1, pix2, factor, mindiff,
1326 &fractdiff, &avediff, details))
1327 return ERROR_INT("diff stats not found", __func__, 1);
1328
1329 if (maxave <= 0.0) maxave = 256.0;
1330 if (fractdiff <= maxfract && avediff <= maxave)
1331 *psimilar = 1;
1332 return 0;
1333}
1334
1335
1378l_ok
1379pixGetDifferenceStats(PIX *pix1,
1380 PIX *pix2,
1381 l_int32 factor,
1382 l_int32 mindiff,
1383 l_float32 *pfractdiff,
1384 l_float32 *pavediff,
1385 l_int32 details)
1386{
1387l_int32 i, first, last, diff;
1388l_float32 fract, ave;
1389l_float32 *array;
1390NUMA *nah, *nan, *nac;
1391
1392 if (pfractdiff) *pfractdiff = 0.0;
1393 if (pavediff) *pavediff = 0.0;
1394 if (!pfractdiff)
1395 return ERROR_INT("&fractdiff not defined", __func__, 1);
1396 if (!pavediff)
1397 return ERROR_INT("&avediff not defined", __func__, 1);
1398 if (!pix1)
1399 return ERROR_INT("pix1 not defined", __func__, 1);
1400 if (!pix2)
1401 return ERROR_INT("pix2 not defined", __func__, 1);
1402 if (mindiff <= 0)
1403 return ERROR_INT("mindiff must be > 0", __func__, 1);
1404
1405 if ((nah = pixGetDifferenceHistogram(pix1, pix2, factor)) == NULL)
1406 return ERROR_INT("na not made", __func__, 1);
1407
1408 if ((nan = numaNormalizeHistogram(nah, 1.0)) == NULL) {
1409 numaDestroy(&nah);
1410 return ERROR_INT("nan not made", __func__, 1);
1411 }
1412 array = numaGetFArray(nan, L_NOCOPY);
1413
1414 if (details) {
1415 lept_mkdir("lept/comp");
1416 numaGetNonzeroRange(nan, 0.0, &first, &last);
1417 nac = numaClipToInterval(nan, first, last);
1418 gplotSimple1(nac, GPLOT_PNG, "/tmp/lept/comp/histo",
1419 "Difference histogram");
1420 l_fileDisplay("/tmp/lept/comp/histo.png", 500, 0, 1.0);
1421 lept_stderr("\nNonzero values in normalized histogram:");
1422 numaWriteStderr(nac);
1423 numaDestroy(&nac);
1424 lept_stderr(" Mindiff fractdiff avediff\n");
1425 lept_stderr(" -----------------------------------\n");
1426 for (diff = 1; diff < L_MIN(2 * mindiff, last); diff++) {
1427 fract = 0.0;
1428 ave = 0.0;
1429 for (i = diff; i <= last; i++) {
1430 fract += array[i];
1431 ave += (l_float32)i * array[i];
1432 }
1433 ave = (fract == 0.0) ? 0.0 : ave / fract;
1434 ave -= diff;
1435 lept_stderr("%5d %7.4f %7.4f\n",
1436 diff, fract, ave);
1437 }
1438 lept_stderr(" -----------------------------------\n");
1439 }
1440
1441 fract = 0.0;
1442 ave = 0.0;
1443 for (i = mindiff; i < 256; i++) {
1444 fract += array[i];
1445 ave += (l_float32)i * array[i];
1446 }
1447 ave = (fract == 0.0) ? 0.0 : ave / fract;
1448 ave -= mindiff;
1449
1450 *pfractdiff = fract;
1451 *pavediff = ave;
1452
1453 numaDestroy(&nah);
1454 numaDestroy(&nan);
1455 return 0;
1456}
1457
1458
1478NUMA *
1479pixGetDifferenceHistogram(PIX *pix1,
1480 PIX *pix2,
1481 l_int32 factor)
1482{
1483l_int32 w1, h1, d1, w2, h2, d2, w, h, wpl1, wpl2;
1484l_int32 i, j, val, val1, val2;
1485l_int32 rval1, rval2, gval1, gval2, bval1, bval2;
1486l_int32 rdiff, gdiff, bdiff, maxdiff;
1487l_uint32 *data1, *data2, *line1, *line2;
1488l_float32 *array;
1489NUMA *na;
1490PIX *pixt1, *pixt2;
1491
1492 if (!pix1)
1493 return (NUMA *)ERROR_PTR("pix1 not defined", __func__, NULL);
1494 if (!pix2)
1495 return (NUMA *)ERROR_PTR("pix2 not defined", __func__, NULL);
1496 d1 = pixGetDepth(pix1);
1497 d2 = pixGetDepth(pix2);
1498 if (d1 == 16 || d2 == 16)
1499 return (NUMA *)ERROR_PTR("d == 16 not supported", __func__, NULL);
1500 if (d1 < 8 && !pixGetColormap(pix1))
1501 return (NUMA *)ERROR_PTR("pix1 depth < 8 bpp and not cmapped",
1502 __func__, NULL);
1503 if (d2 < 8 && !pixGetColormap(pix2))
1504 return (NUMA *)ERROR_PTR("pix2 depth < 8 bpp and not cmapped",
1505 __func__, NULL);
1506 pixt1 = pixRemoveColormap(pix1, REMOVE_CMAP_BASED_ON_SRC);
1507 pixt2 = pixRemoveColormap(pix2, REMOVE_CMAP_BASED_ON_SRC);
1508 pixGetDimensions(pixt1, &w1, &h1, &d1);
1509 pixGetDimensions(pixt2, &w2, &h2, &d2);
1510 if (d1 != d2) {
1511 pixDestroy(&pixt1);
1512 pixDestroy(&pixt2);
1513 return (NUMA *)ERROR_PTR("pix depths not equal", __func__, NULL);
1514 }
1515 if (factor < 1) factor = 1;
1516
1517 na = numaCreate(256);
1518 numaSetCount(na, 256); /* all initialized to 0.0 */
1519 array = numaGetFArray(na, L_NOCOPY);
1520 w = L_MIN(w1, w2);
1521 h = L_MIN(h1, h2);
1522 data1 = pixGetData(pixt1);
1523 data2 = pixGetData(pixt2);
1524 wpl1 = pixGetWpl(pixt1);
1525 wpl2 = pixGetWpl(pixt2);
1526 if (d1 == 8) {
1527 for (i = 0; i < h; i += factor) {
1528 line1 = data1 + i * wpl1;
1529 line2 = data2 + i * wpl2;
1530 for (j = 0; j < w; j += factor) {
1531 val1 = GET_DATA_BYTE(line1, j);
1532 val2 = GET_DATA_BYTE(line2, j);
1533 val = L_ABS(val1 - val2);
1534 array[val]++;
1535 }
1536 }
1537 } else { /* d1 == 32 */
1538 for (i = 0; i < h; i += factor) {
1539 line1 = data1 + i * wpl1;
1540 line2 = data2 + i * wpl2;
1541 for (j = 0; j < w; j += factor) {
1542 extractRGBValues(line1[j], &rval1, &gval1, &bval1);
1543 extractRGBValues(line2[j], &rval2, &gval2, &bval2);
1544 rdiff = L_ABS(rval1 - rval2);
1545 gdiff = L_ABS(gval1 - gval2);
1546 bdiff = L_ABS(bval1 - bval2);
1547 maxdiff = L_MAX(rdiff, gdiff);
1548 maxdiff = L_MAX(maxdiff, bdiff);
1549 array[maxdiff]++;
1550 }
1551 }
1552 }
1553
1554 pixDestroy(&pixt1);
1555 pixDestroy(&pixt2);
1556 return na;
1557}
1558
1559
1607l_ok
1608pixGetPerceptualDiff(PIX *pixs1,
1609 PIX *pixs2,
1610 l_int32 sampling,
1611 l_int32 dilation,
1612 l_int32 mindiff,
1613 l_float32 *pfract,
1614 PIX **ppixdiff1,
1615 PIX **ppixdiff2)
1616{
1617l_int32 d1, d2, w, h, count;
1618PIX *pix1, *pix2, *pix3, *pix4, *pix5, *pix6, *pix7, *pix8, *pix9;
1619PIX *pix10, *pix11;
1620
1621 if (ppixdiff1) *ppixdiff1 = NULL;
1622 if (ppixdiff2) *ppixdiff2 = NULL;
1623 if (!pfract)
1624 return ERROR_INT("&fract not defined", __func__, 1);
1625 *pfract = 1.0; /* init to completely different */
1626 if ((dilation & 1) == 0)
1627 return ERROR_INT("dilation must be odd", __func__, 1);
1628 if (!pixs1)
1629 return ERROR_INT("pixs1 not defined", __func__, 1);
1630 if (!pixs2)
1631 return ERROR_INT("pixs2 not defined", __func__, 1);
1632 d1 = pixGetDepth(pixs1);
1633 d2 = pixGetDepth(pixs2);
1634 if (!pixGetColormap(pixs1) && d1 < 8)
1635 return ERROR_INT("pixs1 not cmapped and < 8 bpp", __func__, 1);
1636 if (!pixGetColormap(pixs2) && d2 < 8)
1637 return ERROR_INT("pixs2 not cmapped and < 8 bpp", __func__, 1);
1638
1639 /* Integer downsample if requested */
1640 if (sampling > 1) {
1641 pix1 = pixScaleByIntSampling(pixs1, sampling);
1642 pix2 = pixScaleByIntSampling(pixs2, sampling);
1643 } else {
1644 pix1 = pixClone(pixs1);
1645 pix2 = pixClone(pixs2);
1646 }
1647
1648 /* Remove colormaps */
1649 if (pixGetColormap(pix1)) {
1650 pix3 = pixRemoveColormap(pix1, REMOVE_CMAP_BASED_ON_SRC);
1651 d1 = pixGetDepth(pix3);
1652 } else {
1653 pix3 = pixClone(pix1);
1654 }
1655 if (pixGetColormap(pix2)) {
1656 pix4 = pixRemoveColormap(pix2, REMOVE_CMAP_BASED_ON_SRC);
1657 d2 = pixGetDepth(pix4);
1658 } else {
1659 pix4 = pixClone(pix2);
1660 }
1661 pixDestroy(&pix1);
1662 pixDestroy(&pix2);
1663 if (d1 != d2 || (d1 != 8 && d1 != 32)) {
1664 pixDestroy(&pix3);
1665 pixDestroy(&pix4);
1666 L_INFO("depths unequal or not in {8,32}: d1 = %d, d2 = %d\n",
1667 __func__, d1, d2);
1668 return 1;
1669 }
1670
1671 /* In each direction, do a small dilation and subtract the dilated
1672 * image from the other image to get a one-sided difference.
1673 * Then take the max of the differences for each direction
1674 * and clipping each component to 255 if necessary. Note that
1675 * for RGB images, the dilations and max selection are done
1676 * component-wise, and the conversion to grayscale also uses the
1677 * maximum component. The resulting grayscale images are
1678 * thresholded using %mindiff. */
1679 if (d1 == 8) {
1680 pix5 = pixDilateGray(pix3, dilation, dilation);
1681 pixCompareGray(pix4, pix5, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1682 &pix7);
1683 pix6 = pixDilateGray(pix4, dilation, dilation);
1684 pixCompareGray(pix3, pix6, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1685 &pix8);
1686 pix9 = pixMinOrMax(NULL, pix7, pix8, L_CHOOSE_MAX);
1687 pix10 = pixThresholdToBinary(pix9, mindiff);
1688 pixInvert(pix10, pix10);
1689 pixCountPixels(pix10, &count, NULL);
1690 pixGetDimensions(pix10, &w, &h, NULL);
1691 *pfract = (w <= 0 || h <= 0) ? 0.0 :
1692 (l_float32)count / (l_float32)(w * h);
1693 pixDestroy(&pix5);
1694 pixDestroy(&pix6);
1695 pixDestroy(&pix7);
1696 pixDestroy(&pix8);
1697 if (ppixdiff1)
1698 *ppixdiff1 = pix9;
1699 else
1700 pixDestroy(&pix9);
1701 if (ppixdiff2)
1702 *ppixdiff2 = pix10;
1703 else
1704 pixDestroy(&pix10);
1705 } else { /* d1 == 32 */
1706 pix5 = pixColorMorph(pix3, L_MORPH_DILATE, dilation, dilation);
1707 pixCompareRGB(pix4, pix5, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1708 &pix7);
1709 pix6 = pixColorMorph(pix4, L_MORPH_DILATE, dilation, dilation);
1710 pixCompareRGB(pix3, pix6, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1711 &pix8);
1712 pix9 = pixMinOrMax(NULL, pix7, pix8, L_CHOOSE_MAX);
1713 pix10 = pixConvertRGBToGrayMinMax(pix9, L_CHOOSE_MAX);
1714 pix11 = pixThresholdToBinary(pix10, mindiff);
1715 pixInvert(pix11, pix11);
1716 pixCountPixels(pix11, &count, NULL);
1717 pixGetDimensions(pix11, &w, &h, NULL);
1718 *pfract = (w <= 0 || h <= 0) ? 0.0 :
1719 (l_float32)count / (l_float32)(w * h);
1720 pixDestroy(&pix5);
1721 pixDestroy(&pix6);
1722 pixDestroy(&pix7);
1723 pixDestroy(&pix8);
1724 pixDestroy(&pix10);
1725 if (ppixdiff1)
1726 *ppixdiff1 = pix9;
1727 else
1728 pixDestroy(&pix9);
1729 if (ppixdiff2)
1730 *ppixdiff2 = pix11;
1731 else
1732 pixDestroy(&pix11);
1733
1734 }
1735 pixDestroy(&pix3);
1736 pixDestroy(&pix4);
1737 return 0;
1738}
1739
1740
1772l_ok
1773pixGetPSNR(PIX *pix1,
1774 PIX *pix2,
1775 l_int32 factor,
1776 l_float32 *ppsnr)
1777{
1778l_int32 same, i, j, w, h, d, wpl1, wpl2, v1, v2, r1, g1, b1, r2, g2, b2;
1779l_uint32 *data1, *data2, *line1, *line2;
1780l_float32 mse; /* mean squared error */
1781
1782 if (!ppsnr)
1783 return ERROR_INT("&psnr not defined", __func__, 1);
1784 *ppsnr = 0.0;
1785 if (!pix1 || !pix2)
1786 return ERROR_INT("empty input pix", __func__, 1);
1787 if (!pixSizesEqual(pix1, pix2))
1788 return ERROR_INT("pix sizes unequal", __func__, 1);
1789 if (pixGetColormap(pix1))
1790 return ERROR_INT("pix1 has colormap", __func__, 1);
1791 if (pixGetColormap(pix2))
1792 return ERROR_INT("pix2 has colormap", __func__, 1);
1793 pixGetDimensions(pix1, &w, &h, &d);
1794 if (d != 8 && d != 32)
1795 return ERROR_INT("pix not 8 or 32 bpp", __func__, 1);
1796 if (factor < 1)
1797 return ERROR_INT("invalid sampling factor", __func__, 1);
1798
1799 pixEqual(pix1, pix2, &same);
1800 if (same) {
1801 *ppsnr = 1000.0; /* crazy big exponent */
1802 return 0;
1803 }
1804
1805 data1 = pixGetData(pix1);
1806 data2 = pixGetData(pix2);
1807 wpl1 = pixGetWpl(pix1);
1808 wpl2 = pixGetWpl(pix2);
1809 mse = 0.0;
1810 if (d == 8) {
1811 for (i = 0; i < h; i += factor) {
1812 line1 = data1 + i * wpl1;
1813 line2 = data2 + i * wpl2;
1814 for (j = 0; j < w; j += factor) {
1815 v1 = GET_DATA_BYTE(line1, j);
1816 v2 = GET_DATA_BYTE(line2, j);
1817 mse += (l_float32)(v1 - v2) * (v1 - v2);
1818 }
1819 }
1820 } else { /* d == 32 */
1821 for (i = 0; i < h; i += factor) {
1822 line1 = data1 + i * wpl1;
1823 line2 = data2 + i * wpl2;
1824 for (j = 0; j < w; j += factor) {
1825 extractRGBValues(line1[j], &r1, &g1, &b1);
1826 extractRGBValues(line2[j], &r2, &g2, &b2);
1827 mse += ((l_float32)(r1 - r2) * (r1 - r2) +
1828 (g1 - g2) * (g1 - g2) +
1829 (b1 - b2) * (b1 - b2)) / 3.0;
1830 }
1831 }
1832 }
1833 mse = mse / ((l_float32)(w) * h);
1834
1835 *ppsnr = -4.3429448 * log(mse / (255 * 255));
1836 return 0;
1837}
1838
1839
1840/*------------------------------------------------------------------*
1841 * Comparison of photo regions by histogram *
1842 *------------------------------------------------------------------*/
1895l_ok
1896pixaComparePhotoRegionsByHisto(PIXA *pixa,
1897 l_float32 minratio,
1898 l_float32 textthresh,
1899 l_int32 factor,
1900 l_int32 n,
1901 l_float32 simthresh,
1902 NUMA **pnai,
1903 l_float32 **pscores,
1904 PIX **ppixd,
1905 l_int32 debug)
1906{
1907char *text;
1908l_int32 i, j, nim, w, h, w1, h1, w2, h2, ival, index, classid;
1909l_float32 score;
1910l_float32 *scores;
1911NUMA *nai, *naw, *nah;
1912NUMAA *naa;
1913NUMAA **n3a; /* array of naa */
1914PIX *pix;
1915
1916 if (pscores) *pscores = NULL;
1917 if (ppixd) *ppixd = NULL;
1918 if (!pnai)
1919 return ERROR_INT("&na not defined", __func__, 1);
1920 *pnai = NULL;
1921 if (!pixa)
1922 return ERROR_INT("pixa not defined", __func__, 1);
1923 if (minratio < 0.0 || minratio > 1.0)
1924 return ERROR_INT("minratio not in [0.0 ... 1.0]", __func__, 1);
1925 if (textthresh <= 0.0) textthresh = 1.3;
1926 if (factor < 1)
1927 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
1928 if (n < 1 || n > 7) {
1929 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
1930 n = 4;
1931 }
1932 if (simthresh <= 0.0) simthresh = 0.25;
1933 if (simthresh > 1.0)
1934 return ERROR_INT("simthresh invalid; should be near 0.25", __func__, 1);
1935
1936 /* Prepare the histograms */
1937 nim = pixaGetCount(pixa);
1938 if ((n3a = (NUMAA **)LEPT_CALLOC(nim, sizeof(NUMAA *))) == NULL)
1939 return ERROR_INT("calloc fail for n3a", __func__, 1);
1940 naw = numaCreate(0);
1941 nah = numaCreate(0);
1942 for (i = 0; i < nim; i++) {
1943 pix = pixaGetPix(pixa, i, L_CLONE);
1944 text = pixGetText(pix);
1945 pixSetResolution(pix, 150, 150);
1946 index = (debug) ? i : 0;
1947 pixGenPhotoHistos(pix, NULL, factor, textthresh, n,
1948 &naa, &w, &h, index);
1949 n3a[i] = naa;
1950 numaAddNumber(naw, w);
1951 numaAddNumber(nah, h);
1952 if (naa)
1953 lept_stderr("Image %s is photo\n", text);
1954 else
1955 lept_stderr("Image %s is NOT photo\n", text);
1956 pixDestroy(&pix);
1957 }
1958
1959 /* Do the comparisons. We are making a set of classes, where
1960 * all similar images are placed in the same class. There are
1961 * 'nim' input images. The classes are labeled by 'classid' (all
1962 * similar images get the same 'classid' value), and 'nai' maps
1963 * the classid of the image in the input array to the classid
1964 * of the similarity class. */
1965 if ((scores =
1966 (l_float32 *)LEPT_CALLOC((size_t)nim * nim, sizeof(l_float32)))
1967 == NULL) {
1968 L_ERROR("calloc fail for scores\n", __func__);
1969 goto cleanup;
1970 }
1971 nai = numaMakeConstant(-1, nim); /* classid array */
1972 for (i = 0, classid = 0; i < nim; i++) {
1973 scores[nim * i + i] = 1.0;
1974 numaGetIValue(nai, i, &ival);
1975 if (ival != -1) /* already set */
1976 continue;
1977 numaSetValue(nai, i, classid);
1978 if (n3a[i] == NULL) { /* not a photo */
1979 classid++;
1980 continue;
1981 }
1982 numaGetIValue(naw, i, &w1);
1983 numaGetIValue(nah, i, &h1);
1984 for (j = i + 1; j < nim; j++) {
1985 numaGetIValue(nai, j, &ival);
1986 if (ival != -1) /* already set */
1987 continue;
1988 if (n3a[j] == NULL) /* not a photo */
1989 continue;
1990 numaGetIValue(naw, j, &w2);
1991 numaGetIValue(nah, j, &h2);
1992 compareTilesByHisto(n3a[i], n3a[j], minratio, w1, h1, w2, h2,
1993 &score, NULL);
1994 scores[nim * i + j] = score;
1995 scores[nim * j + i] = score; /* the score array is symmetric */
1996/* lept_stderr("score = %5.3f\n", score); */
1997 if (score > simthresh) {
1998 numaSetValue(nai, j, classid);
1999 lept_stderr(
2000 "Setting %d similar to %d, in class %d; score %5.3f\n",
2001 j, i, classid, score);
2002 }
2003 }
2004 classid++;
2005 }
2006 *pnai = nai;
2007
2008 /* Debug: optionally save and display the score array.
2009 * All images that are photos are represented by a point on
2010 * the diagonal. Other images in the same similarity class
2011 * are on the same horizontal raster line to the right.
2012 * The array has been symmetrized, so images in the same
2013 * same similarity class also appear on the same column below. */
2014 if (pscores) {
2015 l_int32 wpl, fact;
2016 l_uint32 *line, *data;
2017 PIX *pix2, *pix3;
2018 pix2 = pixCreate(nim, nim, 8);
2019 data = pixGetData(pix2);
2020 wpl = pixGetWpl(pix2);
2021 for (i = 0; i < nim; i++) {
2022 line = data + i * wpl;
2023 for (j = 0; j < nim; j++) {
2024 SET_DATA_BYTE(line, j,
2025 L_MIN(255, 4.0 * 255 * scores[nim * i + j]));
2026 }
2027 }
2028 fact = L_MAX(2, 1000 / nim);
2029 pix3 = pixExpandReplicate(pix2, fact);
2030 lept_stderr("Writing to /tmp/lept/comp/scorearray.png\n");
2031 lept_mkdir("lept/comp");
2032 pixWrite("/tmp/lept/comp/scorearray.png", pix3, IFF_PNG);
2033 pixDestroy(&pix2);
2034 pixDestroy(&pix3);
2035 *pscores = scores;
2036 } else {
2037 LEPT_FREE(scores);
2038 }
2039
2040 /* Debug: optionally display and save the image comparisons.
2041 * Image similarity classes are displayed by column; similar
2042 * images are displayed in the same column. */
2043 if (ppixd)
2044 *ppixd = pixaDisplayTiledByIndex(pixa, nai, 200, 20, 2, 6, 0x0000ff00);
2045
2046cleanup:
2047 numaDestroy(&naw);
2048 numaDestroy(&nah);
2049 for (i = 0; i < nim; i++)
2050 numaaDestroy(&n3a[i]);
2051 LEPT_FREE(n3a);
2052 return 0;
2053}
2054
2055
2115l_ok
2116pixComparePhotoRegionsByHisto(PIX *pix1,
2117 PIX *pix2,
2118 BOX *box1,
2119 BOX *box2,
2120 l_float32 minratio,
2121 l_int32 factor,
2122 l_int32 n,
2123 l_float32 *pscore,
2124 l_int32 debugflag)
2125{
2126l_int32 w1, h1, w2, h2, w1c, h1c, w2c, h2c, debugindex;
2127l_float32 wratio, hratio;
2128NUMAA *naa1, *naa2;
2129PIX *pix3, *pix4;
2130PIXA *pixa;
2131
2132 if (!pscore)
2133 return ERROR_INT("&score not defined", __func__, 1);
2134 *pscore = 0.0;
2135 if (!pix1 || !pix2)
2136 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
2137 if (minratio < 0.5 || minratio > 1.0)
2138 return ERROR_INT("minratio not in [0.5 ... 1.0]", __func__, 1);
2139 if (factor < 1)
2140 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2141 if (n < 1 || n > 7) {
2142 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2143 n = 4;
2144 }
2145
2146 debugindex = 0;
2147 if (debugflag) {
2148 lept_mkdir("lept/comp");
2149 debugindex = 666; /* arbitrary number used for naming output */
2150 }
2151
2152 /* Initial filter by size */
2153 if (box1)
2154 boxGetGeometry(box1, NULL, NULL, &w1, &h1);
2155 else
2156 pixGetDimensions(pix1, &w1, &h1, NULL);
2157 if (box2)
2158 boxGetGeometry(box2, NULL, NULL, &w2, &h2);
2159 else
2160 pixGetDimensions(pix1, &w2, &h2, NULL);
2161 wratio = (w1 < w2) ? (l_float32)w1 / (l_float32)w2 :
2162 (l_float32)w2 / (l_float32)w1;
2163 hratio = (h1 < h2) ? (l_float32)h1 / (l_float32)h2 :
2164 (l_float32)h2 / (l_float32)h1;
2165 if (wratio < minratio || hratio < minratio)
2166 return 0;
2167
2168 /* Initial crop, if necessary, and make histos */
2169 if (box1)
2170 pix3 = pixClipRectangle(pix1, box1, NULL);
2171 else
2172 pix3 = pixClone(pix1);
2173 pixGenPhotoHistos(pix3, NULL, factor, 0, n, &naa1, &w1c, &h1c, debugindex);
2174 pixDestroy(&pix3);
2175 if (!naa1) return 0;
2176 if (box2)
2177 pix4 = pixClipRectangle(pix2, box2, NULL);
2178 else
2179 pix4 = pixClone(pix2);
2180 pixGenPhotoHistos(pix4, NULL, factor, 0, n, &naa2, &w2c, &h2c, debugindex);
2181 pixDestroy(&pix4);
2182 if (!naa2) return 0;
2183
2184 /* Compare histograms */
2185 pixa = (debugflag) ? pixaCreate(0) : NULL;
2186 compareTilesByHisto(naa1, naa2, minratio, w1c, h1c, w2c, h2c, pscore, pixa);
2187 pixaDestroy(&pixa);
2188 return 0;
2189}
2190
2191
2230l_ok
2231pixGenPhotoHistos(PIX *pixs,
2232 BOX *box,
2233 l_int32 factor,
2234 l_float32 thresh,
2235 l_int32 n,
2236 NUMAA **pnaa,
2237 l_int32 *pw,
2238 l_int32 *ph,
2239 l_int32 debugindex)
2240{
2241char buf[64];
2242NUMAA *naa;
2243PIX *pix1, *pix2, *pix3, *pixm;
2244PIXA *pixa;
2245
2246 if (pnaa) *pnaa = NULL;
2247 if (pw) *pw = 0;
2248 if (ph) *ph = 0;
2249 if (!pnaa)
2250 return ERROR_INT("&naa not defined", __func__, 1);
2251 if (!pw || !ph)
2252 return ERROR_INT("&w and &h not both defined", __func__, 1);
2253 if (!pixs || pixGetDepth(pixs) == 1)
2254 return ERROR_INT("pixs not defined or 1 bpp", __func__, 1);
2255 if (factor < 1)
2256 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2257 if (thresh <= 0.0) thresh = 1.3; /* default */
2258 if (n < 1 || n > 7) {
2259 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2260 n = 4;
2261 }
2262
2263 pixa = NULL;
2264 if (debugindex > 0) {
2265 pixa = pixaCreate(0);
2266 lept_mkdir("lept/comp");
2267 }
2268
2269 /* Initial crop, if necessary */
2270 if (box)
2271 pix1 = pixClipRectangle(pixs, box, NULL);
2272 else
2273 pix1 = pixClone(pixs);
2274
2275 /* Convert to 8 bpp and pad to center the centroid */
2276 pix2 = pixConvertTo8(pix1, FALSE);
2277 pix3 = pixPadToCenterCentroid(pix2, factor);
2278
2279 /* Set to 255 all pixels above 230. Do this so that light gray
2280 * pixels do not enter into the comparison. */
2281 pixm = pixThresholdToBinary(pix3, 230);
2282 pixInvert(pixm, pixm);
2283 pixSetMaskedGeneral(pix3, pixm, 255, 0, 0);
2284 pixDestroy(&pixm);
2285
2286 if (debugindex > 0) {
2287 PIX *pix4, *pix5, *pix6, *pix7, *pix8;
2288 PIXA *pixa2;
2289 pix4 = pixConvertTo32(pix2);
2290 pix5 = pixConvertTo32(pix3);
2291 pix6 = pixScaleToSize(pix4, 400, 0);
2292 pix7 = pixScaleToSize(pix5, 400, 0);
2293 pixa2 = pixaCreate(2);
2294 pixaAddPix(pixa2, pix6, L_INSERT);
2295 pixaAddPix(pixa2, pix7, L_INSERT);
2296 pix8 = pixaDisplayTiledInRows(pixa2, 32, 1000, 1.0, 0, 50, 3);
2297 pixaAddPix(pixa, pix8, L_INSERT);
2298 pixDestroy(&pix4);
2299 pixDestroy(&pix5);
2300 pixaDestroy(&pixa2);
2301 }
2302 pixDestroy(&pix1);
2303 pixDestroy(&pix2);
2304
2305 /* Test if this is a photoimage */
2306 pixDecideIfPhotoImage(pix3, factor, thresh, n, &naa, pixa);
2307 if (naa) {
2308 *pnaa = naa;
2309 *pw = pixGetWidth(pix3);
2310 *ph = pixGetHeight(pix3);
2311 }
2312
2313 if (pixa) {
2314 snprintf(buf, sizeof(buf), "/tmp/lept/comp/tiledhistos.%d.pdf",
2315 debugindex);
2316 lept_stderr("Writing to %s\n", buf);
2317 pixaConvertToPdf(pixa, 300, 1.0, L_FLATE_ENCODE, 0, NULL, buf);
2318 pixaDestroy(&pixa);
2319 }
2320
2321 pixDestroy(&pix3);
2322 return 0;
2323}
2324
2325
2341PIX *
2342pixPadToCenterCentroid(PIX *pixs,
2343 l_int32 factor)
2344
2345{
2346l_float32 cx, cy;
2347l_int32 xs, ys, delx, dely, icx, icy, ws, hs, wd, hd;
2348PIX *pix1, *pixd;
2349
2350 if (!pixs)
2351 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2352 if (factor < 1)
2353 return (PIX *)ERROR_PTR("invalid sampling factor", __func__, NULL);
2354
2355 pix1 = pixConvertTo8(pixs, FALSE);
2356 pixCentroid8(pix1, factor, &cx, &cy);
2357 icx = (l_int32)(cx + 0.5);
2358 icy = (l_int32)(cy + 0.5);
2359 pixGetDimensions(pix1, &ws, &hs, NULL);
2360 delx = ws - 2 * icx;
2361 dely = hs - 2 * icy;
2362 xs = L_MAX(0, delx);
2363 ys = L_MAX(0, dely);
2364 wd = 2 * L_MAX(icx, ws - icx);
2365 hd = 2 * L_MAX(icy, hs - icy);
2366 pixd = pixCreate(wd, hd, 8);
2367 pixSetAll(pixd); /* to white */
2368 pixCopyResolution(pixd, pixs);
2369 pixRasterop(pixd, xs, ys, ws, hs, PIX_SRC, pix1, 0, 0);
2370 pixDestroy(&pix1);
2371 return pixd;
2372}
2373
2374
2393l_ok
2394pixCentroid8(PIX *pixs,
2395 l_int32 factor,
2396 l_float32 *pcx,
2397 l_float32 *pcy)
2398{
2399l_int32 i, j, w, h, wpl, val;
2400l_float32 sumx, sumy, sumv;
2401l_uint32 *data, *line;
2402PIX *pix1;
2403
2404 if (pcx) *pcx = 0.0;
2405 if (pcy) *pcy = 0.0;
2406 if (!pixs || pixGetDepth(pixs) != 8)
2407 return ERROR_INT("pixs undefined or not 8 bpp", __func__, 1);
2408 if (factor < 1)
2409 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2410 if (!pcx || !pcy)
2411 return ERROR_INT("&cx and &cy not both defined", __func__, 1);
2412
2413 pix1 = pixInvert(NULL, pixs);
2414 pixGetDimensions(pix1, &w, &h, NULL);
2415 data = pixGetData(pix1);
2416 wpl = pixGetWpl(pix1);
2417 sumx = sumy = sumv = 0.0;
2418 for (i = 0; i < h; i++) {
2419 line = data + i * wpl;
2420 for (j = 0; j < w; j++) {
2421 val = GET_DATA_BYTE(line, j);
2422 sumx += val * j;
2423 sumy += val * i;
2424 sumv += val;
2425 }
2426 }
2427 pixDestroy(&pix1);
2428
2429 if (sumv == 0) {
2430 L_INFO("input image is white\n", __func__);
2431 *pcx = (l_float32)(w) / 2;
2432 *pcy = (l_float32)(h) / 2;
2433 } else {
2434 *pcx = sumx / sumv;
2435 *pcy = sumy / sumv;
2436 }
2437
2438 return 0;
2439}
2440
2441
2477l_ok
2478pixDecideIfPhotoImage(PIX *pix,
2479 l_int32 factor,
2480 l_float32 thresh,
2481 l_int32 n,
2482 NUMAA **pnaa,
2483 PIXA *pixadebug)
2484{
2485char buf[64];
2486l_int32 i, w, h, nx, ny, ngrids, istext, isphoto;
2487l_float32 maxval, sum1, sum2, ratio;
2488L_BMF *bmf;
2489NUMA *na1, *na2, *na3, *narv;
2490NUMAA *naa;
2491PIX *pix1;
2492PIXA *pixa1, *pixa2, *pixa3;
2493
2494 if (!pnaa)
2495 return ERROR_INT("&naa not defined", __func__, 1);
2496 *pnaa = NULL;
2497 if (!pix || pixGetDepth(pix) != 8 || pixGetColormap(pix))
2498 return ERROR_INT("pix undefined or invalid", __func__, 1);
2499 if (n < 1 || n > 7) {
2500 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2501 n = 4;
2502 }
2503 if (thresh <= 0.0) thresh = 1.3; /* default */
2504
2505 /* Look for text lines */
2506 pixDecideIfText(pix, NULL, &istext, pixadebug);
2507 if (istext) {
2508 L_INFO("Image is text\n", __func__);
2509 return 0;
2510 }
2511
2512 /* Determine grid from n */
2513 pixGetDimensions(pix, &w, &h, NULL);
2514 if (w == 0 || h == 0)
2515 return ERROR_INT("invalid pix dimension", __func__, 1);
2516 findHistoGridDimensions(n, w, h, &nx, &ny, 1);
2517
2518 /* Evaluate histograms in each tile */
2519 pixa1 = pixaSplitPix(pix, nx, ny, 0, 0);
2520 ngrids = nx * ny;
2521 bmf = (pixadebug) ? bmfCreate(NULL, 6) : NULL;
2522 naa = numaaCreate(ngrids);
2523 if (pixadebug) {
2524 lept_rmdir("lept/compplot");
2525 lept_mkdir("lept/compplot");
2526 }
2527 for (i = 0; i < ngrids; i++) {
2528 pix1 = pixaGetPix(pixa1, i, L_CLONE);
2529
2530 /* Get histograms, set white count to 0, normalize max to 255 */
2531 na1 = pixGetGrayHistogram(pix1, factor);
2532 numaSetValue(na1, 255, 0);
2533 na2 = numaWindowedMean(na1, 5); /* do some smoothing */
2534 numaGetMax(na2, &maxval, NULL);
2535 na3 = numaTransform(na2, 0, 255.0 / maxval);
2536 if (pixadebug) {
2537 snprintf(buf, sizeof(buf), "/tmp/lept/compplot/plot.%d", i);
2538 gplotSimple1(na3, GPLOT_PNG, buf, "Histos");
2539 }
2540
2541 numaaAddNuma(naa, na3, L_INSERT);
2542 numaDestroy(&na1);
2543 numaDestroy(&na2);
2544 pixDestroy(&pix1);
2545 }
2546 if (pixadebug) {
2547 pix1 = pixaDisplayTiledInColumns(pixa1, nx, 1.0, 30, 2);
2548 pixaAddPix(pixadebug, pix1, L_INSERT);
2549 pixa2 = pixaReadFiles("/tmp/lept/compplot", ".png");
2550 pixa3 = pixaScale(pixa2, 0.4, 0.4);
2551 pix1 = pixaDisplayTiledInColumns(pixa3, nx, 1.0, 30, 2);
2552 pixaAddPix(pixadebug, pix1, L_INSERT);
2553 pixaDestroy(&pixa2);
2554 pixaDestroy(&pixa3);
2555 }
2556
2557 /* Compute the standard deviation between these histos to decide
2558 * if the image is photo or something more like line art,
2559 * which does not support good comparison by tiled histograms. */
2560 grayInterHistogramStats(naa, 5, NULL, NULL, NULL, &narv);
2561
2562 /* For photos, the root variance has a larger weight of
2563 * values in the range [50 ... 150] compared to [200 ... 230],
2564 * than text or line art. For the latter, most of the variance
2565 * between tiles is in the lightest parts of the image, well
2566 * above 150. */
2567 numaGetSumOnInterval(narv, 50, 150, &sum1);
2568 numaGetSumOnInterval(narv, 200, 230, &sum2);
2569 if (sum2 == 0.0) { /* shouldn't happen */
2570 ratio = 0.001; /* anything very small for debug output */
2571 isphoto = 0; /* be conservative */
2572 } else {
2573 ratio = sum1 / sum2;
2574 isphoto = (ratio > thresh) ? 1 : 0;
2575 }
2576 if (pixadebug) {
2577 if (isphoto)
2578 L_INFO("ratio %f > %f; isphoto is true\n",
2579 __func__, ratio, thresh);
2580 else
2581 L_INFO("ratio %f < %f; isphoto is false\n",
2582 __func__, ratio, thresh);
2583 }
2584 if (isphoto)
2585 *pnaa = naa;
2586 else
2587 numaaDestroy(&naa);
2588 bmfDestroy(&bmf);
2589 numaDestroy(&narv);
2590 pixaDestroy(&pixa1);
2591 return 0;
2592}
2593
2594
2620static l_ok
2621findHistoGridDimensions(l_int32 n,
2622 l_int32 w,
2623 l_int32 h,
2624 l_int32 *pnx,
2625 l_int32 *pny,
2626 l_int32 debug)
2627{
2628l_int32 nx, ny, max;
2629l_float32 ratio;
2630
2631 ratio = (l_float32)w / (l_float32)h;
2632 max = n * n;
2633 nx = ny = n;
2634 while (nx > 1 && ny > 1) {
2635 if (ratio > 2.0) { /* reduce ny */
2636 ny--;
2637 nx = max / ny;
2638 if (debug)
2639 lept_stderr("nx = %d, ny = %d, ratio w/h = %4.2f\n",
2640 nx, ny, ratio);
2641 } else if (ratio < 0.5) { /* reduce nx */
2642 nx--;
2643 ny = max / nx;
2644 if (debug)
2645 lept_stderr("nx = %d, ny = %d, ratio w/h = %4.2f\n",
2646 nx, ny, ratio);
2647 } else { /* we're ok */
2648 if (debug)
2649 lept_stderr("nx = %d, ny = %d, ratio w/h = %4.2f\n",
2650 nx, ny, ratio);
2651 break;
2652 }
2653 ratio = (l_float32)(ny * w) / (l_float32)(nx * h);
2654 }
2655 *pnx = nx;
2656 *pny = ny;
2657 return 0;
2658}
2659
2660
2683l_ok
2684compareTilesByHisto(NUMAA *naa1,
2685 NUMAA *naa2,
2686 l_float32 minratio,
2687 l_int32 w1,
2688 l_int32 h1,
2689 l_int32 w2,
2690 l_int32 h2,
2691 l_float32 *pscore,
2692 PIXA *pixadebug)
2693{
2694char buf1[128], buf2[128];
2695l_int32 i, n;
2696l_float32 wratio, hratio, score, minscore, dist;
2697L_BMF *bmf;
2698NUMA *na1, *na2, *nadist, *nascore;
2699
2700 if (!pscore)
2701 return ERROR_INT("&score not defined", __func__, 1);
2702 *pscore = 0.0;
2703 if (!naa1 || !naa2)
2704 return ERROR_INT("naa1 and naa2 not both defined", __func__, 1);
2705
2706 /* Filter for different sizes */
2707 wratio = (w1 < w2) ? (l_float32)w1 / (l_float32)w2 :
2708 (l_float32)w2 / (l_float32)w1;
2709 hratio = (h1 < h2) ? (l_float32)h1 / (l_float32)h2 :
2710 (l_float32)h2 / (l_float32)h1;
2711 if (wratio < minratio || hratio < minratio) {
2712 if (pixadebug)
2713 L_INFO("Sizes differ: wratio = %f, hratio = %f\n",
2714 __func__, wratio, hratio);
2715 return 0;
2716 }
2717 n = numaaGetCount(naa1);
2718 if (n != numaaGetCount(naa2)) { /* due to differing w/h ratio */
2719 L_INFO("naa1 and naa2 sizes are different\n", __func__);
2720 return 0;
2721 }
2722
2723 if (pixadebug) {
2724 lept_rmdir("lept/comptile");
2725 lept_mkdir("lept/comptile");
2726 }
2727
2728
2729 /* Evaluate histograms in each tile. Remove white before
2730 * computing EMD, because there are may be a lot of white
2731 * pixels due to padding, and we don't want to include them.
2732 * This also makes the debug histo plots more informative. */
2733 minscore = 1.0;
2734 nadist = numaCreate(n);
2735 nascore = numaCreate(n);
2736 bmf = (pixadebug) ? bmfCreate(NULL, 6) : NULL;
2737 for (i = 0; i < n; i++) {
2738 na1 = numaaGetNuma(naa1, i, L_CLONE);
2739 na2 = numaaGetNuma(naa2, i, L_CLONE);
2740 numaSetValue(na1, 255, 0.0);
2741 numaSetValue(na2, 255, 0.0);
2742
2743 /* To compare histograms, use the normalized earthmover distance.
2744 * Further normalize to get the EM distance as a fraction of the
2745 * maximum distance in the histogram (255). Finally, scale this
2746 * up by 10.0, and subtract from 1.0 to get a similarity score. */
2747 numaEarthMoverDistance(na1, na2, &dist);
2748 score = L_MAX(0.0, 1.0 - 10.0 * (dist / 255.));
2749 numaAddNumber(nadist, dist);
2750 numaAddNumber(nascore, score);
2751 minscore = L_MIN(minscore, score);
2752 if (pixadebug) {
2753 snprintf(buf1, sizeof(buf1), "/tmp/lept/comptile/plot.%d", i);
2754 gplotSimple2(na1, na2, GPLOT_PNG, buf1, "Histos");
2755 }
2756 numaDestroy(&na1);
2757 numaDestroy(&na2);
2758 }
2759 *pscore = minscore;
2760
2761 if (pixadebug) {
2762 for (i = 0; i < n; i++) {
2763 PIX *pix1, *pix2;
2764 snprintf(buf1, sizeof(buf1), "/tmp/lept/comptile/plot.%d.png", i);
2765 pix1 = pixRead(buf1);
2766 numaGetFValue(nadist, i, &dist);
2767 numaGetFValue(nascore, i, &score);
2768 snprintf(buf2, sizeof(buf2),
2769 "Image %d\ndist = %5.3f, score = %5.3f", i, dist, score);
2770 pix2 = pixAddTextlines(pix1, bmf, buf2, 0x0000ff00, L_ADD_BELOW);
2771 pixaAddPix(pixadebug, pix2, L_INSERT);
2772 pixDestroy(&pix1);
2773 }
2774 lept_stderr("Writing to /tmp/lept/comptile/comparegray.pdf\n");
2775 pixaConvertToPdf(pixadebug, 300, 1.0, L_FLATE_ENCODE, 0, NULL,
2776 "/tmp/lept/comptile/comparegray.pdf");
2777 numaWriteDebug("/tmp/lept/comptile/scores.na", nascore);
2778 numaWriteDebug("/tmp/lept/comptile/dists.na", nadist);
2779 }
2780
2781 bmfDestroy(&bmf);
2782 numaDestroy(&nadist);
2783 numaDestroy(&nascore);
2784 return 0;
2785}
2786
2787
2859l_ok
2860pixCompareGrayByHisto(PIX *pix1,
2861 PIX *pix2,
2862 BOX *box1,
2863 BOX *box2,
2864 l_float32 minratio,
2865 l_int32 maxgray,
2866 l_int32 factor,
2867 l_int32 n,
2868 l_float32 *pscore,
2869 l_int32 debugflag)
2870{
2871l_int32 w1, h1, w2, h2;
2872l_float32 wratio, hratio;
2873BOX *box3, *box4;
2874PIX *pix3, *pix4, *pix5, *pix6, *pix7, *pix8;
2875PIXA *pixa;
2876
2877 if (!pscore)
2878 return ERROR_INT("&score not defined", __func__, 1);
2879 *pscore = 0.0;
2880 if (!pix1 || !pix2)
2881 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
2882 if (minratio < 0.5 || minratio > 1.0)
2883 return ERROR_INT("minratio not in [0.5 ... 1.0]", __func__, 1);
2884 if (maxgray < 200)
2885 return ERROR_INT("invalid maxgray; should be >= 200", __func__, 1);
2886 maxgray = L_MIN(255, maxgray);
2887 if (factor < 1)
2888 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2889 if (n < 1 || n > 7) {
2890 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2891 n = 4;
2892 }
2893
2894 if (debugflag)
2895 lept_mkdir("lept/comp");
2896
2897 /* Initial filter by size */
2898 if (box1)
2899 boxGetGeometry(box1, NULL, NULL, &w1, &h1);
2900 else
2901 pixGetDimensions(pix1, &w1, &h1, NULL);
2902 if (box2)
2903 boxGetGeometry(box2, NULL, NULL, &w2, &h2);
2904 else
2905 pixGetDimensions(pix1, &w2, &h2, NULL);
2906 wratio = (w1 < w2) ? (l_float32)w1 / (l_float32)w2 :
2907 (l_float32)w2 / (l_float32)w1;
2908 hratio = (h1 < h2) ? (l_float32)h1 / (l_float32)h2 :
2909 (l_float32)h2 / (l_float32)h1;
2910 if (wratio < minratio || hratio < minratio)
2911 return 0;
2912
2913 /* Initial crop, if necessary */
2914 if (box1)
2915 pix3 = pixClipRectangle(pix1, box1, NULL);
2916 else
2917 pix3 = pixClone(pix1);
2918 if (box2)
2919 pix4 = pixClipRectangle(pix2, box2, NULL);
2920 else
2921 pix4 = pixClone(pix2);
2922
2923 /* Convert to 8 bpp, align centroids and do maximal crop */
2924 pix5 = pixConvertTo8(pix3, FALSE);
2925 pix6 = pixConvertTo8(pix4, FALSE);
2926 pixCropAlignedToCentroid(pix5, pix6, factor, &box3, &box4);
2927 pix7 = pixClipRectangle(pix5, box3, NULL);
2928 pix8 = pixClipRectangle(pix6, box4, NULL);
2929 pixa = (debugflag) ? pixaCreate(0) : NULL;
2930 if (debugflag) {
2931 PIX *pix9, *pix10, *pix11, *pix12, *pix13;
2932 PIXA *pixa2;
2933 pix9 = pixConvertTo32(pix5);
2934 pix10 = pixConvertTo32(pix6);
2935 pixRenderBoxArb(pix9, box3, 2, 255, 0, 0);
2936 pixRenderBoxArb(pix10, box4, 2, 255, 0, 0);
2937 pix11 = pixScaleToSize(pix9, 400, 0);
2938 pix12 = pixScaleToSize(pix10, 400, 0);
2939 pixa2 = pixaCreate(2);
2940 pixaAddPix(pixa2, pix11, L_INSERT);
2941 pixaAddPix(pixa2, pix12, L_INSERT);
2942 pix13 = pixaDisplayTiledInRows(pixa2, 32, 1000, 1.0, 0, 50, 0);
2943 pixaAddPix(pixa, pix13, L_INSERT);
2944 pixDestroy(&pix9);
2945 pixDestroy(&pix10);
2946 pixaDestroy(&pixa2);
2947 }
2948 pixDestroy(&pix3);
2949 pixDestroy(&pix4);
2950 pixDestroy(&pix5);
2951 pixDestroy(&pix6);
2952 boxDestroy(&box3);
2953 boxDestroy(&box4);
2954
2955 /* Tile and compare histograms */
2956 pixCompareTilesByHisto(pix7, pix8, maxgray, factor, n, pscore, pixa);
2957 pixaDestroy(&pixa);
2958 pixDestroy(&pix7);
2959 pixDestroy(&pix8);
2960 return 0;
2961}
2962
2963
2984static l_ok
2985pixCompareTilesByHisto(PIX *pix1,
2986 PIX *pix2,
2987 l_int32 maxgray,
2988 l_int32 factor,
2989 l_int32 n,
2990 l_float32 *pscore,
2991 PIXA *pixadebug)
2992{
2993char buf[64];
2994l_int32 w, h, i, j, nx, ny, ngr;
2995l_float32 score, minscore, maxval1, maxval2, dist;
2996L_BMF *bmf;
2997NUMA *na1, *na2, *na3, *na4, *na5, *na6, *na7;
2998PIX *pix3, *pix4;
2999PIXA *pixa1, *pixa2;
3000
3001 if (!pscore)
3002 return ERROR_INT("&score not defined", __func__, 1);
3003 *pscore = 0.0;
3004 if (!pix1 || !pix2)
3005 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
3006
3007 /* Determine grid from n */
3008 pixGetDimensions(pix1, &w, &h, NULL);
3009 findHistoGridDimensions(n, w, h, &nx, &ny, 1);
3010 ngr = nx * ny;
3011
3012 /* Evaluate histograms in each tile */
3013 pixa1 = pixaSplitPix(pix1, nx, ny, 0, 0);
3014 pixa2 = pixaSplitPix(pix2, nx, ny, 0, 0);
3015 na7 = (pixadebug) ? numaCreate(ngr) : NULL;
3016 bmf = (pixadebug) ? bmfCreate(NULL, 6) : NULL;
3017 minscore = 1.0;
3018 for (i = 0; i < ngr; i++) {
3019 pix3 = pixaGetPix(pixa1, i, L_CLONE);
3020 pix4 = pixaGetPix(pixa2, i, L_CLONE);
3021
3022 /* Get histograms, set white count to 0, normalize max to 255 */
3023 na1 = pixGetGrayHistogram(pix3, factor);
3024 na2 = pixGetGrayHistogram(pix4, factor);
3025 if (maxgray < 255) {
3026 for (j = maxgray + 1; j <= 255; j++) {
3027 numaSetValue(na1, j, 0);
3028 numaSetValue(na2, j, 0);
3029 }
3030 }
3031 na3 = numaWindowedMean(na1, 5);
3032 na4 = numaWindowedMean(na2, 5);
3033 numaGetMax(na3, &maxval1, NULL);
3034 numaGetMax(na4, &maxval2, NULL);
3035 na5 = numaTransform(na3, 0, 255.0 / maxval1);
3036 na6 = numaTransform(na4, 0, 255.0 / maxval2);
3037 if (pixadebug) {
3038 gplotSimple2(na5, na6, GPLOT_PNG, "/tmp/lept/comp/plot1", "Histos");
3039 }
3040
3041 /* To compare histograms, use the normalized earthmover distance.
3042 * Further normalize to get the EM distance as a fraction of the
3043 * maximum distance in the histogram (255). Finally, scale this
3044 * up by 10.0, and subtract from 1.0 to get a similarity score. */
3045 numaEarthMoverDistance(na5, na6, &dist);
3046 score = L_MAX(0.0, 1.0 - 8.0 * (dist / 255.));
3047 if (pixadebug) numaAddNumber(na7, score);
3048 minscore = L_MIN(minscore, score);
3049 if (pixadebug) {
3050 PIX *pix5, *pix6, *pix7, *pix8, *pix9, *pix10;
3051 PIXA *pixa3;
3052 l_int32 w, h, wscale;
3053 pixa3 = pixaCreate(3);
3054 pixGetDimensions(pix3, &w, &h, NULL);
3055 wscale = (w > h) ? 700 : 400;
3056 pix5 = pixScaleToSize(pix3, wscale, 0);
3057 pix6 = pixScaleToSize(pix4, wscale, 0);
3058 pixaAddPix(pixa3, pix5, L_INSERT);
3059 pixaAddPix(pixa3, pix6, L_INSERT);
3060 pix7 = pixRead("/tmp/lept/comp/plot1.png");
3061 pix8 = pixScaleToSize(pix7, 700, 0);
3062 snprintf(buf, sizeof(buf), "%5.3f", score);
3063 pix9 = pixAddTextlines(pix8, bmf, buf, 0x0000ff00, L_ADD_RIGHT);
3064 pixaAddPix(pixa3, pix9, L_INSERT);
3065 pix10 = pixaDisplayTiledInRows(pixa3, 32, 1000, 1.0, 0, 50, 0);
3066 pixaAddPix(pixadebug, pix10, L_INSERT);
3067 pixDestroy(&pix7);
3068 pixDestroy(&pix8);
3069 pixaDestroy(&pixa3);
3070 }
3071 numaDestroy(&na1);
3072 numaDestroy(&na2);
3073 numaDestroy(&na3);
3074 numaDestroy(&na4);
3075 numaDestroy(&na5);
3076 numaDestroy(&na6);
3077 pixDestroy(&pix3);
3078 pixDestroy(&pix4);
3079 }
3080 *pscore = minscore;
3081
3082 if (pixadebug) {
3083 pixaConvertToPdf(pixadebug, 300, 1.0, L_FLATE_ENCODE, 0, NULL,
3084 "/tmp/lept/comp/comparegray.pdf");
3085 numaWriteDebug("/tmp/lept/comp/tilescores.na", na7);
3086 }
3087
3088 bmfDestroy(&bmf);
3089 numaDestroy(&na7);
3090 pixaDestroy(&pixa1);
3091 pixaDestroy(&pixa2);
3092 return 0;
3093}
3094
3095
3112l_ok
3113pixCropAlignedToCentroid(PIX *pix1,
3114 PIX *pix2,
3115 l_int32 factor,
3116 BOX **pbox1,
3117 BOX **pbox2)
3118{
3119l_float32 cx1, cy1, cx2, cy2;
3120l_int32 w1, h1, w2, h2, icx1, icy1, icx2, icy2;
3121l_int32 xm, xm1, xm2, xp, xp1, xp2, ym, ym1, ym2, yp, yp1, yp2;
3122PIX *pix3, *pix4;
3123
3124 if (pbox1) *pbox1 = NULL;
3125 if (pbox2) *pbox2 = NULL;
3126 if (!pix1 || !pix2)
3127 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
3128 if (factor < 1)
3129 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
3130 if (!pbox1 || !pbox2)
3131 return ERROR_INT("&box1 and &box2 not both defined", __func__, 1);
3132
3133 pix3 = pixConvertTo8(pix1, FALSE);
3134 pix4 = pixConvertTo8(pix2, FALSE);
3135 pixCentroid8(pix3, factor, &cx1, &cy1);
3136 pixCentroid8(pix4, factor, &cx2, &cy2);
3137 pixGetDimensions(pix3, &w1, &h1, NULL);
3138 pixGetDimensions(pix4, &w2, &h2, NULL);
3139 pixDestroy(&pix3);
3140 pixDestroy(&pix4);
3141
3142 icx1 = (l_int32)(cx1 + 0.5);
3143 icy1 = (l_int32)(cy1 + 0.5);
3144 icx2 = (l_int32)(cx2 + 0.5);
3145 icy2 = (l_int32)(cy2 + 0.5);
3146 xm = L_MIN(icx1, icx2);
3147 xm1 = icx1 - xm;
3148 xm2 = icx2 - xm;
3149 xp = L_MIN(w1 - icx1, w2 - icx2); /* one pixel beyond to the right */
3150 xp1 = icx1 + xp;
3151 xp2 = icx2 + xp;
3152 ym = L_MIN(icy1, icy2);
3153 ym1 = icy1 - ym;
3154 ym2 = icy2 - ym;
3155 yp = L_MIN(h1 - icy1, h2 - icy2); /* one pixel below the bottom */
3156 yp1 = icy1 + yp;
3157 yp2 = icy2 + yp;
3158 *pbox1 = boxCreate(xm1, ym1, xp1 - xm1, yp1 - ym1);
3159 *pbox2 = boxCreate(xm2, ym2, xp2 - xm2, yp2 - ym2);
3160 return 0;
3161}
3162
3163
3185l_uint8 *
3186l_compressGrayHistograms(NUMAA *naa,
3187 l_int32 w,
3188 l_int32 h,
3189 size_t *psize)
3190{
3191l_uint8 *bytea;
3192l_int32 i, j, n, nn, ival;
3193l_float32 maxval;
3194NUMA *na1, *na2;
3195
3196 if (!psize)
3197 return (l_uint8 *)ERROR_PTR("&size not defined", __func__, NULL);
3198 *psize = 0;
3199 if (!naa)
3200 return (l_uint8 *)ERROR_PTR("naa not defined", __func__, NULL);
3201 n = numaaGetCount(naa);
3202 for (i = 0; i < n; i++) {
3203 nn = numaaGetNumaCount(naa, i);
3204 if (nn != 256) {
3205 L_ERROR("%d numbers in numa[%d]\n", __func__, nn, i);
3206 return NULL;
3207 }
3208 }
3209
3210 if ((bytea = (l_uint8 *)LEPT_CALLOC(8 + 256 * n, sizeof(l_uint8))) == NULL)
3211 return (l_uint8 *)ERROR_PTR("bytea not made", __func__, NULL);
3212 *psize = 8 + 256 * n;
3213 l_setDataFourBytes(bytea, 0, w);
3214 l_setDataFourBytes(bytea, 1, h);
3215 for (i = 0; i < n; i++) {
3216 na1 = numaaGetNuma(naa, i, L_COPY);
3217 numaGetMax(na1, &maxval, NULL);
3218 na2 = numaTransform(na1, 0, 255.0 / maxval);
3219 for (j = 0; j < 256; j++) {
3220 numaGetIValue(na2, j, &ival);
3221 bytea[8 + 256 * i + j] = ival;
3222 }
3223 numaDestroy(&na1);
3224 numaDestroy(&na2);
3225 }
3226
3227 return bytea;
3228}
3229
3230
3251NUMAA *
3252l_uncompressGrayHistograms(l_uint8 *bytea,
3253 size_t size,
3254 l_int32 *pw,
3255 l_int32 *ph)
3256{
3257l_int32 i, j, n;
3258NUMA *na;
3259NUMAA *naa;
3260
3261 if (pw) *pw = 0;
3262 if (ph) *ph = 0;
3263 if (!pw || !ph)
3264 return (NUMAA *)ERROR_PTR("&w and &h not both defined", __func__, NULL);
3265 if (!bytea)
3266 return (NUMAA *)ERROR_PTR("bytea not defined", __func__, NULL);
3267 n = (size - 8) / 256;
3268 if ((size - 8) % 256 != 0)
3269 return (NUMAA *)ERROR_PTR("bytea size is invalid", __func__, NULL);
3270
3271 *pw = l_getDataFourBytes(bytea, 0);
3272 *ph = l_getDataFourBytes(bytea, 1);
3273 naa = numaaCreate(n);
3274 for (i = 0; i < n; i++) {
3275 na = numaCreate(256);
3276 for (j = 0; j < 256; j++)
3277 numaAddNumber(na, bytea[8 + 256 * i + j]);
3278 numaaAddNuma(naa, na, L_INSERT);
3279 }
3280
3281 return naa;
3282}
3283
3284
3285/*------------------------------------------------------------------*
3286 * Translated images at the same resolution *
3287 *------------------------------------------------------------------*/
3318l_ok
3319pixCompareWithTranslation(PIX *pix1,
3320 PIX *pix2,
3321 l_int32 thresh,
3322 l_int32 *pdelx,
3323 l_int32 *pdely,
3324 l_float32 *pscore,
3325 l_int32 debugflag)
3326{
3327l_uint8 *subtab;
3328l_int32 i, level, area1, area2, delx, dely;
3329l_int32 etransx, etransy, maxshift, dbint;
3330l_int32 *stab, *ctab;
3331l_float32 cx1, cx2, cy1, cy2, score;
3332PIX *pixb1, *pixb2, *pixt1, *pixt2, *pixt3, *pixt4;
3333PIXA *pixa1, *pixa2, *pixadb;
3334
3335 if (pdelx) *pdelx = 0;
3336 if (pdely) *pdely = 0;
3337 if (pscore) *pscore = 0.0;
3338 if (!pdelx || !pdely)
3339 return ERROR_INT("&delx and &dely not defined", __func__, 1);
3340 if (!pscore)
3341 return ERROR_INT("&score not defined", __func__, 1);
3342 if (!pix1)
3343 return ERROR_INT("pix1 not defined", __func__, 1);
3344 if (!pix2)
3345 return ERROR_INT("pix2 not defined", __func__, 1);
3346
3347 /* Make tables */
3348 subtab = makeSubsampleTab2x();
3349 stab = makePixelSumTab8();
3350 ctab = makePixelCentroidTab8();
3351
3352 /* Binarize each image */
3353 pixb1 = pixConvertTo1(pix1, thresh);
3354 pixb2 = pixConvertTo1(pix2, thresh);
3355
3356 /* Make a cascade of 2x reduced images for each, thresholding
3357 * with level 2 (neutral), down to 8x reduction */
3358 pixa1 = pixaCreate(4);
3359 pixa2 = pixaCreate(4);
3360 if (debugflag)
3361 pixadb = pixaCreate(4);
3362 pixaAddPix(pixa1, pixb1, L_INSERT);
3363 pixaAddPix(pixa2, pixb2, L_INSERT);
3364 for (i = 0; i < 3; i++) {
3365 pixt1 = pixReduceRankBinary2(pixb1, 2, subtab);
3366 pixt2 = pixReduceRankBinary2(pixb2, 2, subtab);
3367 pixaAddPix(pixa1, pixt1, L_INSERT);
3368 pixaAddPix(pixa2, pixt2, L_INSERT);
3369 pixb1 = pixt1;
3370 pixb2 = pixt2;
3371 }
3372
3373 /* At the lowest level, use the centroids with a maxshift of 6
3374 * to search for the best alignment. Then at higher levels,
3375 * use the result from the level below as the initial approximation
3376 * for the alignment, and search with a maxshift of 2. */
3377 for (level = 3; level >= 0; level--) {
3378 pixt1 = pixaGetPix(pixa1, level, L_CLONE);
3379 pixt2 = pixaGetPix(pixa2, level, L_CLONE);
3380 pixCountPixels(pixt1, &area1, stab);
3381 pixCountPixels(pixt2, &area2, stab);
3382 if (level == 3) {
3383 pixCentroid(pixt1, ctab, stab, &cx1, &cy1);
3384 pixCentroid(pixt2, ctab, stab, &cx2, &cy2);
3385 etransx = lept_roundftoi(cx1 - cx2);
3386 etransy = lept_roundftoi(cy1 - cy2);
3387 maxshift = 6;
3388 } else {
3389 etransx = 2 * delx;
3390 etransy = 2 * dely;
3391 maxshift = 2;
3392 }
3393 dbint = (debugflag) ? level + 1 : 0;
3394 pixBestCorrelation(pixt1, pixt2, area1, area2, etransx, etransy,
3395 maxshift, stab, &delx, &dely, &score, dbint);
3396 if (debugflag) {
3397 lept_stderr("Level %d: delx = %d, dely = %d, score = %7.4f\n",
3398 level, delx, dely, score);
3399 pixRasteropIP(pixt2, delx, dely, L_BRING_IN_WHITE);
3400 pixt3 = pixDisplayDiffBinary(pixt1, pixt2);
3401 pixt4 = pixExpandReplicate(pixt3, 8 / (1 << (3 - level)));
3402 pixaAddPix(pixadb, pixt4, L_INSERT);
3403 pixDestroy(&pixt3);
3404 }
3405 pixDestroy(&pixt1);
3406 pixDestroy(&pixt2);
3407 }
3408
3409 if (debugflag) {
3410 pixaConvertToPdf(pixadb, 300, 1.0, L_FLATE_ENCODE, 0, NULL,
3411 "/tmp/lept/comp/compare.pdf");
3412 convertFilesToPdf("/tmp/lept/comp", "correl_", 30, 1.0, L_FLATE_ENCODE,
3413 0, "Correlation scores at levels 1 through 5",
3414 "/tmp/lept/comp/correl.pdf");
3415 pixaDestroy(&pixadb);
3416 }
3417
3418 *pdelx = delx;
3419 *pdely = dely;
3420 *pscore = score;
3421 pixaDestroy(&pixa1);
3422 pixaDestroy(&pixa2);
3423 LEPT_FREE(subtab);
3424 LEPT_FREE(stab);
3425 LEPT_FREE(ctab);
3426 return 0;
3427}
3428
3429
3470l_ok
3471pixBestCorrelation(PIX *pix1,
3472 PIX *pix2,
3473 l_int32 area1,
3474 l_int32 area2,
3475 l_int32 etransx,
3476 l_int32 etransy,
3477 l_int32 maxshift,
3478 l_int32 *tab8,
3479 l_int32 *pdelx,
3480 l_int32 *pdely,
3481 l_float32 *pscore,
3482 l_int32 debugflag)
3483{
3484l_int32 shiftx, shifty, delx, dely;
3485l_int32 *tab;
3486l_float32 maxscore, score;
3487FPIX *fpix;
3488PIX *pix3, *pix4;
3489
3490 if (pdelx) *pdelx = 0;
3491 if (pdely) *pdely = 0;
3492 if (pscore) *pscore = 0.0;
3493 if (!pix1 || pixGetDepth(pix1) != 1)
3494 return ERROR_INT("pix1 not defined or not 1 bpp", __func__, 1);
3495 if (!pix2 || pixGetDepth(pix2) != 1)
3496 return ERROR_INT("pix2 not defined or not 1 bpp", __func__, 1);
3497 if (!area1 || !area2)
3498 return ERROR_INT("areas must be > 0", __func__, 1);
3499
3500 if (debugflag > 0)
3501 fpix = fpixCreate(2 * maxshift + 1, 2 * maxshift + 1);
3502
3503 if (!tab8)
3504 tab = makePixelSumTab8();
3505 else
3506 tab = tab8;
3507
3508 /* Search over a set of {shiftx, shifty} for the max */
3509 maxscore = 0;
3510 delx = etransx;
3511 dely = etransy;
3512 for (shifty = -maxshift; shifty <= maxshift; shifty++) {
3513 for (shiftx = -maxshift; shiftx <= maxshift; shiftx++) {
3514 pixCorrelationScoreShifted(pix1, pix2, area1, area2,
3515 etransx + shiftx,
3516 etransy + shifty, tab, &score);
3517 if (debugflag > 0) {
3518 fpixSetPixel(fpix, maxshift + shiftx, maxshift + shifty,
3519 1000.0 * score);
3520/* lept_stderr("(sx, sy) = (%d, %d): score = %6.4f\n",
3521 shiftx, shifty, score); */
3522 }
3523 if (score > maxscore) {
3524 maxscore = score;
3525 delx = etransx + shiftx;
3526 dely = etransy + shifty;
3527 }
3528 }
3529 }
3530
3531 if (debugflag > 0) {
3532 char buf[128];
3533 lept_mkdir("lept/comp");
3534 pix3 = fpixDisplayMaxDynamicRange(fpix);
3535 pix4 = pixExpandReplicate(pix3, 20);
3536 snprintf(buf, sizeof(buf), "/tmp/lept/comp/correl_%d.png",
3537 debugflag);
3538 pixWrite(buf, pix4, IFF_PNG);
3539 pixDestroy(&pix3);
3540 pixDestroy(&pix4);
3541 fpixDestroy(&fpix);
3542 }
3543
3544 if (pdelx) *pdelx = delx;
3545 if (pdely) *pdely = dely;
3546 if (pscore) *pscore = maxscore;
3547 if (!tab8) LEPT_FREE(tab);
3548 return 0;
3549}
NUMAA
struct Numaa NUMAA
Definition array.h:69
NUMA
struct Numa NUMA
Definition array.h:66
GET_DATA_BYTE
#define GET_DATA_BYTE(pdata, n)
Definition arrayaccess.h:188
SET_DATA_BYTE
#define SET_DATA_BYTE(pdata, n, val)
Definition arrayaccess.h:198
pixCompareGrayOrRGB
l_ok pixCompareGrayOrRGB(PIX *pix1, PIX *pix2, l_int32 comptype, l_int32 plottype, l_int32 *psame, l_float32 *pdiff, l_float32 *prmsdiff, PIX **ppixdiff)
pixCompareGrayOrRGB()
Definition compare.c:785
compareTilesByHisto
l_ok compareTilesByHisto(NUMAA *naa1, NUMAA *naa2, l_float32 minratio, l_int32 w1, l_int32 h1, l_int32 w2, l_int32 h2, l_float32 *pscore, PIXA *pixadebug)
compareTilesByHisto()
Definition compare.c:2684
pixaComparePhotoRegionsByHisto
l_ok pixaComparePhotoRegionsByHisto(PIXA *pixa, l_float32 minratio, l_float32 textthresh, l_int32 factor, l_int32 n, l_float32 simthresh, NUMA **pnai, l_float32 **pscores, PIX **ppixd, l_int32 debug)
pixaComparePhotoRegionsByHisto()
Definition compare.c:1896
pixGetDifferenceStats
l_ok pixGetDifferenceStats(PIX *pix1, PIX *pix2, l_int32 factor, l_int32 mindiff, l_float32 *pfractdiff, l_float32 *pavediff, l_int32 details)
pixGetDifferenceStats()
Definition compare.c:1379
pixCropAlignedToCentroid
l_ok pixCropAlignedToCentroid(PIX *pix1, PIX *pix2, l_int32 factor, BOX **pbox1, BOX **pbox2)
pixCropAlignedToCentroid()
Definition compare.c:3113
pixCompareGrayByHisto
l_ok pixCompareGrayByHisto(PIX *pix1, PIX *pix2, BOX *box1, BOX *box2, l_float32 minratio, l_int32 maxgray, l_int32 factor, l_int32 n, l_float32 *pscore, l_int32 debugflag)
pixCompareGrayByHisto()
Definition compare.c:2860
pixDecideIfPhotoImage
l_ok pixDecideIfPhotoImage(PIX *pix, l_int32 factor, l_float32 thresh, l_int32 n, NUMAA **pnaa, PIXA *pixadebug)
pixDecideIfPhotoImage()
Definition compare.c:2478
pixDisplayDiffBinary
PIX * pixDisplayDiffBinary(PIX *pix1, PIX *pix2)
pixDisplayDiffBinary()
Definition compare.c:652
pixCompareRGB
l_ok pixCompareRGB(PIX *pix1, PIX *pix2, l_int32 comptype, l_int32 plottype, l_int32 *psame, l_float32 *pdiff, l_float32 *prmsdiff, PIX **ppixdiff)
pixCompareRGB()
Definition compare.c:973
pixGenPhotoHistos
l_ok pixGenPhotoHistos(PIX *pixs, BOX *box, l_int32 factor, l_float32 thresh, l_int32 n, NUMAA **pnaa, l_int32 *pw, l_int32 *ph, l_int32 debugindex)
pixGenPhotoHistos()
Definition compare.c:2231
pixGetPSNR
l_ok pixGetPSNR(PIX *pix1, PIX *pix2, l_int32 factor, l_float32 *ppsnr)
pixGetPSNR()
Definition compare.c:1773
pixBestCorrelation
l_ok pixBestCorrelation(PIX *pix1, PIX *pix2, l_int32 area1, l_int32 area2, l_int32 etransx, l_int32 etransy, l_int32 maxshift, l_int32 *tab8, l_int32 *pdelx, l_int32 *pdely, l_float32 *pscore, l_int32 debugflag)
pixBestCorrelation()
Definition compare.c:3471
cmapEqual
l_ok cmapEqual(PIXCMAP *cmap1, PIXCMAP *cmap2, l_int32 ncomps, l_int32 *psame)
cmapEqual()
Definition compare.c:476
pixGetDifferenceHistogram
NUMA * pixGetDifferenceHistogram(PIX *pix1, PIX *pix2, l_int32 factor)
pixGetDifferenceHistogram()
Definition compare.c:1479
pixCompareTilesByHisto
static l_ok pixCompareTilesByHisto(PIX *pix1, PIX *pix2, l_int32 maxgray, l_int32 factor, l_int32 n, l_float32 *pscore, PIXA *pixadebug)
pixCompareTilesByHisto()
Definition compare.c:2985
pixCompareWithTranslation
l_ok pixCompareWithTranslation(PIX *pix1, PIX *pix2, l_int32 thresh, l_int32 *pdelx, l_int32 *pdely, l_float32 *pscore, l_int32 debugflag)
pixCompareWithTranslation()
Definition compare.c:3319
l_compressGrayHistograms
l_uint8 * l_compressGrayHistograms(NUMAA *naa, l_int32 w, l_int32 h, size_t *psize)
l_compressGrayHistograms()
Definition compare.c:3186
pixEqual
l_ok pixEqual(PIX *pix1, PIX *pix2, l_int32 *psame)
pixEqual()
Definition compare.c:156
pixPadToCenterCentroid
PIX * pixPadToCenterCentroid(PIX *pixs, l_int32 factor)
pixPadToCenterCentroid()
Definition compare.c:2342
pixComparePhotoRegionsByHisto
l_ok pixComparePhotoRegionsByHisto(PIX *pix1, PIX *pix2, BOX *box1, BOX *box2, l_float32 minratio, l_int32 factor, l_int32 n, l_float32 *pscore, l_int32 debugflag)
pixComparePhotoRegionsByHisto()
Definition compare.c:2116
pixEqualWithCmap
l_ok pixEqualWithCmap(PIX *pix1, PIX *pix2, l_int32 *psame)
pixEqualWithCmap()
Definition compare.c:382
pixCompareGray
l_ok pixCompareGray(PIX *pix1, PIX *pix2, l_int32 comptype, l_int32 plottype, l_int32 *psame, l_float32 *pdiff, l_float32 *prmsdiff, PIX **ppixdiff)
pixCompareGray()
Definition compare.c:866
pixTestForSimilarity
l_ok pixTestForSimilarity(PIX *pix1, PIX *pix2, l_int32 factor, l_int32 mindiff, l_float32 maxfract, l_float32 maxave, l_int32 *psimilar, l_int32 details)
pixTestForSimilarity()
Definition compare.c:1302
pixCompareRankDifference
NUMA * pixCompareRankDifference(PIX *pix1, PIX *pix2, l_int32 factor)
pixCompareRankDifference()
Definition compare.c:1219
pixCorrelationBinary
l_ok pixCorrelationBinary(PIX *pix1, PIX *pix2, l_float32 *pval)
pixCorrelationBinary()
Definition compare.c:596
pixCompareTiled
l_ok pixCompareTiled(PIX *pix1, PIX *pix2, l_int32 sx, l_int32 sy, l_int32 type, PIX **ppixdiff)
pixCompareTiled()
Definition compare.c:1127
pixCompareBinary
l_ok pixCompareBinary(PIX *pix1, PIX *pix2, l_int32 comptype, l_float32 *pfract, PIX **ppixdiff)
pixCompareBinary()
Definition compare.c:707
pixUsesCmapColor
l_ok pixUsesCmapColor(PIX *pixs, l_int32 *pcolor)
pixUsesCmapColor()
Definition compare.c:532
findHistoGridDimensions
static l_ok findHistoGridDimensions(l_int32 n, l_int32 w, l_int32 h, l_int32 *pnx, l_int32 *pny, l_int32 debug)
findHistoGridDimensions()
Definition compare.c:2621
pixGetPerceptualDiff
l_ok pixGetPerceptualDiff(PIX *pixs1, PIX *pixs2, l_int32 sampling, l_int32 dilation, l_int32 mindiff, l_float32 *pfract, PIX **ppixdiff1, PIX **ppixdiff2)
pixGetPerceptualDiff()
Definition compare.c:1608
l_uncompressGrayHistograms
NUMAA * l_uncompressGrayHistograms(l_uint8 *bytea, size_t size, l_int32 *pw, l_int32 *ph)
l_uncompressGrayHistograms()
Definition compare.c:3252
pixCentroid8
l_ok pixCentroid8(PIX *pixs, l_int32 factor, l_float32 *pcx, l_float32 *pcy)
pixCentroid8()
Definition compare.c:2394
pixEqualWithAlpha
l_ok pixEqualWithAlpha(PIX *pix1, PIX *pix2, l_int32 use_alpha, l_int32 *psame)
pixEqualWithAlpha()
Definition compare.c:182
L_FLATE_ENCODE
@ L_FLATE_ENCODE
Definition imageio.h:161
L_BRING_IN_WHITE
@ L_BRING_IN_WHITE
Definition pix.h:662
L_ROOT_MEAN_SQUARE
@ L_ROOT_MEAN_SQUARE
Definition pix.h:765
L_MEAN_ABSVAL
@ L_MEAN_ABSVAL
Definition pix.h:761
REMOVE_CMAP_TO_FULL_COLOR
@ REMOVE_CMAP_TO_FULL_COLOR
Definition pix.h:382
REMOVE_CMAP_TO_GRAYSCALE
@ REMOVE_CMAP_TO_GRAYSCALE
Definition pix.h:381
REMOVE_CMAP_BASED_ON_SRC
@ REMOVE_CMAP_BASED_ON_SRC
Definition pix.h:384
PIXACC
struct Pixacc PIXACC
Definition pix.h:273
FPIX
struct FPix FPIX
Definition pix.h:285
L_COPY
@ L_COPY
Definition pix.h:505
L_CLONE
@ L_CLONE
Definition pix.h:506
L_NOCOPY
@ L_NOCOPY
Definition pix.h:503
L_INSERT
@ L_INSERT
Definition pix.h:504
PIX
struct Pix PIX
Definition pix.h:228
BOX
struct Box BOX
Definition pix.h:252
PIX_SRC
#define PIX_SRC
Definition pix.h:444
PIXCMAP
struct PixColormap PIXCMAP
Definition pix.h:231
L_ADD_BELOW
@ L_ADD_BELOW
Definition pix.h:1003
L_ADD_RIGHT
@ L_ADD_RIGHT
Definition pix.h:1005
COLOR_BLUE
@ COLOR_BLUE
Definition pix.h:330
COLOR_RED
@ COLOR_RED
Definition pix.h:328
COLOR_GREEN
@ COLOR_GREEN
Definition pix.h:329
PIXA
struct Pixa PIXA
Definition pix.h:243