Leptonica 1.84.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
106#ifdef HAVE_CONFIG_H
107#include <config_auto.h>
108#endif /* HAVE_CONFIG_H */
109
110#include <string.h>
111#include <math.h>
112#include "allheaders.h"
113
114 /* Small enough to consider equal to 0.0, for plot output */
115static const l_float32 TINY = 0.00001f;
116
117static l_ok findHistoGridDimensions(l_int32 n, l_int32 w, l_int32 h,
118 l_int32 *pnx, l_int32 *pny, l_int32 debug);
119static l_ok pixCompareTilesByHisto(PIX *pix1, PIX *pix2, l_int32 maxgray,
120 l_int32 factor, l_int32 n,
121 l_float32 *pscore, PIXA *pixadebug);
122
123/*------------------------------------------------------------------*
124 * Test for pix equality *
125 *------------------------------------------------------------------*/
156l_ok
157pixEqual(PIX *pix1,
158 PIX *pix2,
159 l_int32 *psame)
160{
161 return pixEqualWithAlpha(pix1, pix2, 0, psame);
162}
163
164
182l_ok
183pixEqualWithAlpha(PIX *pix1,
184 PIX *pix2,
185 l_int32 use_alpha,
186 l_int32 *psame)
187{
188l_int32 w1, h1, d1, w2, h2, d2, wpl1, wpl2;
189l_int32 spp1, spp2, i, j, color, mismatch, opaque;
190l_int32 fullwords, linebits, endbits;
191l_uint32 endmask, wordmask;
192l_uint32 *data1, *data2, *line1, *line2;
193PIX *pixs1, *pixs2, *pixt1, *pixt2, *pixalpha;
194PIXCMAP *cmap1, *cmap2;
195
196 if (!psame)
197 return ERROR_INT("psame not defined", __func__, 1);
198 *psame = 0; /* init to not equal */
199 if (!pix1 || !pix2)
200 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
201 pixGetDimensions(pix1, &w1, &h1, &d1);
202 pixGetDimensions(pix2, &w2, &h2, &d2);
203 if (w1 != w2 || h1 != h2) {
204 L_INFO("pix sizes differ\n", __func__);
205 return 0;
206 }
207
208 /* Suppose the use_alpha flag is true.
209 * If only one of two 32 bpp images has spp == 4, we call that
210 * a "mismatch" of the alpha component. In the case of a mismatch,
211 * if the 4 bpp pix does not have all alpha components opaque (255),
212 * the images are not-equal. However if they are all opaque,
213 * this image is equivalent to spp == 3, so we allow the
214 * comparison to go forward, testing only for the RGB equality. */
215 spp1 = pixGetSpp(pix1);
216 spp2 = pixGetSpp(pix2);
217 mismatch = 0;
218 if (use_alpha && d1 == 32 && d2 == 32) {
219 mismatch = ((spp1 == 4 && spp2 != 4) || (spp1 != 4 && spp2 == 4));
220 if (mismatch) {
221 pixalpha = (spp1 == 4) ? pix1 : pix2;
222 pixAlphaIsOpaque(pixalpha, &opaque);
223 if (!opaque) {
224 L_INFO("just one pix has a non-opaque alpha layer\n", __func__);
225 return 0;
226 }
227 }
228 }
229
230 cmap1 = pixGetColormap(pix1);
231 cmap2 = pixGetColormap(pix2);
232 if (!cmap1 && !cmap2 && (d1 != d2) && (d1 == 32 || d2 == 32)) {
233 L_INFO("no colormaps, pix depths unequal, and one of them is RGB\n",
234 __func__);
235 return 0;
236 }
237
238 if (cmap1 && cmap2 && (d1 == d2)) /* use special function */
239 return pixEqualWithCmap(pix1, pix2, psame);
240
241 /* Must remove colormaps if they exist, and in the process
242 * end up with the resulting images having the same depth. */
243 if (cmap1 && !cmap2) {
244 pixUsesCmapColor(pix1, &color);
245 if (color && d2 <= 8) /* can't be equal */
246 return 0;
247 if (d2 < 8)
248 pixs2 = pixConvertTo8(pix2, FALSE);
249 else
250 pixs2 = pixClone(pix2);
251 if (d2 <= 8)
252 pixs1 = pixRemoveColormap(pix1, REMOVE_CMAP_TO_GRAYSCALE);
253 else
254 pixs1 = pixRemoveColormap(pix1, REMOVE_CMAP_TO_FULL_COLOR);
255 } else if (!cmap1 && cmap2) {
256 pixUsesCmapColor(pix2, &color);
257 if (color && d1 <= 8) /* can't be equal */
258 return 0;
259 if (d1 < 8)
260 pixs1 = pixConvertTo8(pix1, FALSE);
261 else
262 pixs1 = pixClone(pix1);
263 if (d1 <= 8)
264 pixs2 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_GRAYSCALE);
265 else
266 pixs2 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_FULL_COLOR);
267 } else if (cmap1 && cmap2) { /* depths not equal; use rgb */
268 pixs1 = pixRemoveColormap(pix1, REMOVE_CMAP_TO_FULL_COLOR);
269 pixs2 = pixRemoveColormap(pix2, REMOVE_CMAP_TO_FULL_COLOR);
270 } else { /* no colormaps */
271 pixs1 = pixClone(pix1);
272 pixs2 = pixClone(pix2);
273 }
274
275 /* OK, we have no colormaps, but the depths may still be different */
276 d1 = pixGetDepth(pixs1);
277 d2 = pixGetDepth(pixs2);
278 if (d1 != d2) {
279 if (d1 == 16 || d2 == 16) {
280 L_INFO("one pix is 16 bpp\n", __func__);
281 pixDestroy(&pixs1);
282 pixDestroy(&pixs2);
283 return 0;
284 }
285 pixt1 = pixConvertLossless(pixs1, 8);
286 pixt2 = pixConvertLossless(pixs2, 8);
287 if (!pixt1 || !pixt2) {
288 L_INFO("failure to convert to 8 bpp\n", __func__);
289 pixDestroy(&pixs1);
290 pixDestroy(&pixs2);
291 pixDestroy(&pixt1);
292 pixDestroy(&pixt2);
293 return 0;
294 }
295 } else {
296 pixt1 = pixClone(pixs1);
297 pixt2 = pixClone(pixs2);
298 }
299 pixDestroy(&pixs1);
300 pixDestroy(&pixs2);
301
302 /* No colormaps, equal depths; do pixel comparisons */
303 d1 = pixGetDepth(pixt1);
304 d2 = pixGetDepth(pixt2);
305 wpl1 = pixGetWpl(pixt1);
306 wpl2 = pixGetWpl(pixt2);
307 data1 = pixGetData(pixt1);
308 data2 = pixGetData(pixt2);
309
310 if (d1 == 32) { /* test either RGB or RGBA pixels */
311 if (use_alpha && !mismatch)
312 wordmask = (spp1 == 3) ? 0xffffff00 : 0xffffffff;
313 else
314 wordmask = 0xffffff00;
315 for (i = 0; i < h1; i++) {
316 line1 = data1 + wpl1 * i;
317 line2 = data2 + wpl2 * i;
318 for (j = 0; j < wpl1; j++) {
319 if ((*line1 ^ *line2) & wordmask) {
320 pixDestroy(&pixt1);
321 pixDestroy(&pixt2);
322 return 0;
323 }
324 line1++;
325 line2++;
326 }
327 }
328 } else { /* all bits count */
329 linebits = d1 * w1;
330 fullwords = linebits / 32;
331 endbits = linebits & 31;
332 endmask = (endbits == 0) ? 0 : (0xffffffff << (32 - endbits));
333 for (i = 0; i < h1; i++) {
334 line1 = data1 + wpl1 * i;
335 line2 = data2 + wpl2 * i;
336 for (j = 0; j < fullwords; j++) {
337 if (*line1 ^ *line2) {
338 pixDestroy(&pixt1);
339 pixDestroy(&pixt2);
340 return 0;
341 }
342 line1++;
343 line2++;
344 }
345 if (endbits) {
346 if ((*line1 ^ *line2) & endmask) {
347 pixDestroy(&pixt1);
348 pixDestroy(&pixt2);
349 return 0;
350 }
351 }
352 }
353 }
354
355 pixDestroy(&pixt1);
356 pixDestroy(&pixt2);
357 *psame = 1;
358 return 0;
359}
360
361
382l_ok
383pixEqualWithCmap(PIX *pix1,
384 PIX *pix2,
385 l_int32 *psame)
386{
387l_int32 d, w, h, wpl1, wpl2, i, j, linebits, fullwords, endbits;
388l_int32 rval1, rval2, gval1, gval2, bval1, bval2, samecmaps;
389l_uint32 endmask, val1, val2;
390l_uint32 *data1, *data2, *line1, *line2;
391PIXCMAP *cmap1, *cmap2;
392
393 if (!psame)
394 return ERROR_INT("&same not defined", __func__, 1);
395 *psame = 0;
396 if (!pix1)
397 return ERROR_INT("pix1 not defined", __func__, 1);
398 if (!pix2)
399 return ERROR_INT("pix2 not defined", __func__, 1);
400
401 if (pixSizesEqual(pix1, pix2) == 0)
402 return 0;
403 cmap1 = pixGetColormap(pix1);
404 cmap2 = pixGetColormap(pix2);
405 if (!cmap1 || !cmap2) {
406 L_INFO("both images don't have colormap\n", __func__);
407 return 0;
408 }
409 pixGetDimensions(pix1, &w, &h, &d);
410 if (d != 1 && d != 2 && d != 4 && d != 8) {
411 L_INFO("pix depth not in {1, 2, 4, 8}\n", __func__);
412 return 0;
413 }
414
415 cmapEqual(cmap1, cmap2, 3, &samecmaps);
416 if (samecmaps == TRUE) { /* colormaps are identical; compare by words */
417 linebits = d * w;
418 wpl1 = pixGetWpl(pix1);
419 wpl2 = pixGetWpl(pix2);
420 data1 = pixGetData(pix1);
421 data2 = pixGetData(pix2);
422 fullwords = linebits / 32;
423 endbits = linebits & 31;
424 endmask = (endbits == 0) ? 0 : (0xffffffff << (32 - endbits));
425 for (i = 0; i < h; i++) {
426 line1 = data1 + wpl1 * i;
427 line2 = data2 + wpl2 * i;
428 for (j = 0; j < fullwords; j++) {
429 if (*line1 ^ *line2)
430 return 0;
431 line1++;
432 line2++;
433 }
434 if (endbits) {
435 if ((*line1 ^ *line2) & endmask)
436 return 0;
437 }
438 }
439 *psame = 1;
440 return 0;
441 }
442
443 /* Colormaps aren't identical; compare pixel by pixel */
444 for (i = 0; i < h; i++) {
445 for (j = 0; j < w; j++) {
446 pixGetPixel(pix1, j, i, &val1);
447 pixGetPixel(pix2, j, i, &val2);
448 pixcmapGetColor(cmap1, val1, &rval1, &gval1, &bval1);
449 pixcmapGetColor(cmap2, val2, &rval2, &gval2, &bval2);
450 if (rval1 != rval2 || gval1 != gval2 || bval1 != bval2)
451 return 0;
452 }
453 }
454
455 *psame = 1;
456 return 0;
457}
458
459
476l_ok
477cmapEqual(PIXCMAP *cmap1,
478 PIXCMAP *cmap2,
479 l_int32 ncomps,
480 l_int32 *psame)
481{
482l_int32 n1, n2, i, rval1, rval2, gval1, gval2, bval1, bval2, aval1, aval2;
483
484 if (!psame)
485 return ERROR_INT("&same not defined", __func__, 1);
486 *psame = FALSE;
487 if (!cmap1)
488 return ERROR_INT("cmap1 not defined", __func__, 1);
489 if (!cmap2)
490 return ERROR_INT("cmap2 not defined", __func__, 1);
491 if (ncomps != 3 && ncomps != 4)
492 return ERROR_INT("ncomps not 3 or 4", __func__, 1);
493
494 n1 = pixcmapGetCount(cmap1);
495 n2 = pixcmapGetCount(cmap2);
496 if (n1 != n2) {
497 L_INFO("colormap sizes are different\n", __func__);
498 return 0;
499 }
500
501 for (i = 0; i < n1; i++) {
502 pixcmapGetRGBA(cmap1, i, &rval1, &gval1, &bval1, &aval1);
503 pixcmapGetRGBA(cmap2, i, &rval2, &gval2, &bval2, &aval2);
504 if (rval1 != rval2 || gval1 != gval2 || bval1 != bval2)
505 return 0;
506 if (ncomps == 4 && aval1 != aval2)
507 return 0;
508 }
509 *psame = TRUE;
510 return 0;
511}
512
513
532l_ok
533pixUsesCmapColor(PIX *pixs,
534 l_int32 *pcolor)
535{
536l_int32 n, i, rval, gval, bval, numpix;
537NUMA *na;
538PIXCMAP *cmap;
539
540 if (!pcolor)
541 return ERROR_INT("&color not defined", __func__, 1);
542 *pcolor = 0;
543 if (!pixs)
544 return ERROR_INT("pixs not defined", __func__, 1);
545
546 if ((cmap = pixGetColormap(pixs)) == NULL)
547 return 0;
548
549 pixcmapHasColor(cmap, pcolor);
550 if (*pcolor == 0) /* no color */
551 return 0;
552
553 /* The cmap has color entries. Are they used? */
554 na = pixGetGrayHistogram(pixs, 1);
555 n = pixcmapGetCount(cmap);
556 for (i = 0; i < n; i++) {
557 pixcmapGetColor(cmap, i, &rval, &gval, &bval);
558 numaGetIValue(na, i, &numpix);
559 if ((rval != gval || rval != bval) && numpix) { /* color found! */
560 *pcolor = 1;
561 break;
562 }
563 }
564 numaDestroy(&na);
565
566 return 0;
567}
568
569
570/*------------------------------------------------------------------*
571 * Binary correlation *
572 *------------------------------------------------------------------*/
596l_ok
597pixCorrelationBinary(PIX *pix1,
598 PIX *pix2,
599 l_float32 *pval)
600{
601l_int32 count1, count2, countn;
602l_int32 *tab8;
603PIX *pixn;
604
605 if (!pval)
606 return ERROR_INT("&pval not defined", __func__, 1);
607 *pval = 0.0;
608 if (!pix1)
609 return ERROR_INT("pix1 not defined", __func__, 1);
610 if (!pix2)
611 return ERROR_INT("pix2 not defined", __func__, 1);
612
613 tab8 = makePixelSumTab8();
614 pixCountPixels(pix1, &count1, tab8);
615 pixCountPixels(pix2, &count2, tab8);
616 if (count1 == 0 || count2 == 0) {
617 LEPT_FREE(tab8);
618 return 0;
619 }
620 pixn = pixAnd(NULL, pix1, pix2);
621 pixCountPixels(pixn, &countn, tab8);
622 *pval = (l_float32)countn * (l_float32)countn /
623 ((l_float32)count1 * (l_float32)count2);
624 LEPT_FREE(tab8);
625 pixDestroy(&pixn);
626 return 0;
627}
628
629
630/*------------------------------------------------------------------*
631 * Difference of two images *
632 *------------------------------------------------------------------*/
656PIX *
657pixDisplayDiff(PIX *pix1,
658 PIX *pix2,
659 l_int32 showall,
660 l_int32 mindiff,
661 l_uint32 diffcolor)
662{
663l_int32 i, j, w1, h1, d1, w2, h2, d2, minw, minh, wpl1, wpl2, wpl3;
664l_int32 rval1, gval1, bval1, rval2, gval2, bval2;
665l_uint32 val1, val2;
666l_uint32 *data1, *data2, *data3, *line1, *line2, *line3;
667PIX *pix3 = NULL, *pix4 = NULL, *pixd;
668PIXA *pixa1;
669
670 if (!pix1 || !pix2)
671 return (PIX *)ERROR_PTR("pix1, pix2 not both defined", __func__, NULL);
672 pixGetDimensions(pix1, &w1, &h1, &d1);
673 pixGetDimensions(pix2, &w2, &h2, &d2);
674 if (d1 != d2)
675 return (PIX *)ERROR_PTR("unequal depths", __func__, NULL);
676 if (mindiff <= 0)
677 return (PIX *)ERROR_PTR("mindiff must be > 0", __func__, NULL);
678
679 if (d1 == 1) {
680 pix3 = pixDisplayDiffBinary(pix1, pix2);
681 pixd = pixConvertTo32(pix3);
682 pixDestroy(&pix3);
683 } else {
684 minw = L_MIN(w1, w2);
685 minh = L_MIN(h1, h2);
686 pix3 = pixConvertTo32(pix1);
687 pix4 = pixConvertTo32(pix2);
688 pixd = pixCreate(minw, minh, 32);
689 pixRasterop(pixd, 0, 0, minw, minh, PIX_SRC, pix3, 0, 0);
690 data1 = pixGetData(pix3);
691 wpl1 = pixGetWpl(pix3);
692 data2 = pixGetData(pix4);
693 wpl2 = pixGetWpl(pix4);
694 data3 = pixGetData(pixd);
695 wpl3 = pixGetWpl(pixd);
696 for (i = 0; i < minh; i++) {
697 line1 = data1 + i * wpl1;
698 line2 = data2 + i * wpl2;
699 line3 = data3 + i * wpl3;
700 for (j = 0; j < minw; j++) {
701 val1 = GET_DATA_FOUR_BYTES(line1, j);
702 val2 = GET_DATA_FOUR_BYTES(line2, j);
703 extractRGBValues(val1, &rval1, &gval1, &bval1);
704 extractRGBValues(val2, &rval2, &gval2, &bval2);
705 if (L_ABS(rval1 - rval2) >= mindiff ||
706 L_ABS(gval1 - gval2) >= mindiff ||
707 L_ABS(bval1 - bval2) >= mindiff)
708 SET_DATA_FOUR_BYTES(line3, j, diffcolor);
709 }
710 }
711 }
712
713 if (showall) {
714 pixa1 = pixaCreate(3);
715 if (d1 == 1) {
716 pixaAddPix(pixa1, pix1, L_COPY);
717 pixaAddPix(pixa1, pix2, L_COPY);
718 } else {
719 pixaAddPix(pixa1, pix3, L_INSERT);
720 pixaAddPix(pixa1, pix4, L_INSERT);
721 }
722 pixaAddPix(pixa1, pixd, L_INSERT); /* save diff image */
723 pixd = pixaDisplayTiledInColumns(pixa1, 2, 1.0, 30, 2); /* all 3 */
724 pixaDestroy(&pixa1);
725 }
726 return pixd;
727}
728
729
749PIX *
750pixDisplayDiffBinary(PIX *pix1,
751 PIX *pix2)
752{
753l_int32 w1, h1, d1, w2, h2, d2, minw, minh;
754PIX *pixt, *pixd;
755PIXCMAP *cmap;
756
757 if (!pix1 || !pix2)
758 return (PIX *)ERROR_PTR("pix1, pix2 not both defined", __func__, NULL);
759 pixGetDimensions(pix1, &w1, &h1, &d1);
760 pixGetDimensions(pix2, &w2, &h2, &d2);
761 if (d1 != 1 || d2 != 1)
762 return (PIX *)ERROR_PTR("pix1 and pix2 not 1 bpp", __func__, NULL);
763 minw = L_MIN(w1, w2);
764 minh = L_MIN(h1, h2);
765
766 pixd = pixCreate(minw, minh, 4);
767 cmap = pixcmapCreate(4);
768 pixcmapAddColor(cmap, 255, 255, 255); /* initialized to white */
769 pixcmapAddColor(cmap, 0, 0, 0);
770 pixcmapAddColor(cmap, 255, 0, 0);
771 pixcmapAddColor(cmap, 0, 255, 0);
772 pixSetColormap(pixd, cmap);
773
774 pixt = pixAnd(NULL, pix1, pix2);
775 pixPaintThroughMask(pixd, pixt, 0, 0, 0x0); /* black */
776 pixSubtract(pixt, pix1, pix2);
777 pixPaintThroughMask(pixd, pixt, 0, 0, 0xff000000); /* red */
778 pixSubtract(pixt, pix2, pix1);
779 pixPaintThroughMask(pixd, pixt, 0, 0, 0x00ff0000); /* green */
780 pixDestroy(&pixt);
781 return pixd;
782}
783
784
804l_ok
805pixCompareBinary(PIX *pix1,
806 PIX *pix2,
807 l_int32 comptype,
808 l_float32 *pfract,
809 PIX **ppixdiff)
810{
811l_int32 w, h, count;
812PIX *pixt;
813
814 if (ppixdiff) *ppixdiff = NULL;
815 if (!pfract)
816 return ERROR_INT("&pfract not defined", __func__, 1);
817 *pfract = 1.0; /* initialize to max difference */
818 if (!pix1 || pixGetDepth(pix1) != 1)
819 return ERROR_INT("pix1 not defined or not 1 bpp", __func__, 1);
820 if (!pix2 || pixGetDepth(pix2) != 1)
821 return ERROR_INT("pix2 not defined or not 1 bpp", __func__, 1);
822 if (comptype != L_COMPARE_XOR && comptype != L_COMPARE_SUBTRACT)
823 return ERROR_INT("invalid comptype", __func__, 1);
824
825 if (comptype == L_COMPARE_XOR)
826 pixt = pixXor(NULL, pix1, pix2);
827 else /* comptype == L_COMPARE_SUBTRACT) */
828 pixt = pixSubtract(NULL, pix1, pix2);
829 pixCountPixels(pixt, &count, NULL);
830 pixGetDimensions(pix1, &w, &h, NULL);
831 *pfract = (l_float32)(count) / (l_float32)(w * h);
832
833 if (ppixdiff)
834 *ppixdiff = pixt;
835 else
836 pixDestroy(&pixt);
837 return 0;
838}
839
840
882l_ok
883pixCompareGrayOrRGB(PIX *pix1,
884 PIX *pix2,
885 l_int32 comptype,
886 l_int32 plottype,
887 l_int32 *psame,
888 l_float32 *pdiff,
889 l_float32 *prmsdiff,
890 PIX **ppixdiff)
891{
892l_int32 retval, d1, d2;
893PIX *pixt1, *pixt2, *pixs1, *pixs2;
894
895 if (psame) *psame = 0;
896 if (pdiff) *pdiff = 255.0;
897 if (prmsdiff) *prmsdiff = 255.0;
898 if (ppixdiff) *ppixdiff = NULL;
899 if (!pix1 || pixGetDepth(pix1) == 1)
900 return ERROR_INT("pix1 not defined or 1 bpp", __func__, 1);
901 if (!pix2 || pixGetDepth(pix2) == 1)
902 return ERROR_INT("pix2 not defined or 1 bpp", __func__, 1);
903 if (comptype != L_COMPARE_SUBTRACT && comptype != L_COMPARE_ABS_DIFF)
904 return ERROR_INT("invalid comptype", __func__, 1);
905 if (plottype < 0 || plottype >= NUM_GPLOT_OUTPUTS)
906 return ERROR_INT("invalid plottype", __func__, 1);
907
908 pixt1 = pixRemoveColormap(pix1, REMOVE_CMAP_BASED_ON_SRC);
909 pixt2 = pixRemoveColormap(pix2, REMOVE_CMAP_BASED_ON_SRC);
910 d1 = pixGetDepth(pixt1);
911 d2 = pixGetDepth(pixt2);
912 if (d1 < 8)
913 pixs1 = pixConvertTo8(pixt1, FALSE);
914 else
915 pixs1 = pixClone(pixt1);
916 if (d2 < 8)
917 pixs2 = pixConvertTo8(pixt2, FALSE);
918 else
919 pixs2 = pixClone(pixt2);
920 pixDestroy(&pixt1);
921 pixDestroy(&pixt2);
922 d1 = pixGetDepth(pixs1);
923 d2 = pixGetDepth(pixs2);
924 if (d1 != d2) {
925 pixDestroy(&pixs1);
926 pixDestroy(&pixs2);
927 return ERROR_INT("intrinsic depths are not equal", __func__, 1);
928 }
929
930 if (d1 == 8 || d1 == 16)
931 retval = pixCompareGray(pixs1, pixs2, comptype, plottype, psame,
932 pdiff, prmsdiff, ppixdiff);
933 else /* d1 == 32 */
934 retval = pixCompareRGB(pixs1, pixs2, comptype, plottype, psame,
935 pdiff, prmsdiff, ppixdiff);
936 pixDestroy(&pixs1);
937 pixDestroy(&pixs2);
938 return retval;
939}
940
941
963l_ok
964pixCompareGray(PIX *pix1,
965 PIX *pix2,
966 l_int32 comptype,
967 l_int32 plottype,
968 l_int32 *psame,
969 l_float32 *pdiff,
970 l_float32 *prmsdiff,
971 PIX **ppixdiff)
972{
973char buf[64];
974static l_atomic index = 0;
975l_int32 d1, d2, same, first, last;
976GPLOT *gplot;
977NUMA *na, *nac;
978PIX *pixt;
979
980 if (psame) *psame = 0;
981 if (pdiff) *pdiff = 255.0;
982 if (prmsdiff) *prmsdiff = 255.0;
983 if (ppixdiff) *ppixdiff = NULL;
984 if (!pix1)
985 return ERROR_INT("pix1 not defined", __func__, 1);
986 if (!pix2)
987 return ERROR_INT("pix2 not defined", __func__, 1);
988 d1 = pixGetDepth(pix1);
989 d2 = pixGetDepth(pix2);
990 if ((d1 != d2) || (d1 != 8 && d1 != 16))
991 return ERROR_INT("depths unequal or not 8 or 16 bpp", __func__, 1);
992 if (pixGetColormap(pix1) || pixGetColormap(pix2))
993 return ERROR_INT("pix1 and/or pix2 are colormapped", __func__, 1);
994 if (comptype != L_COMPARE_SUBTRACT && comptype != L_COMPARE_ABS_DIFF)
995 return ERROR_INT("invalid comptype", __func__, 1);
996 if (plottype < 0 || plottype >= NUM_GPLOT_OUTPUTS)
997 return ERROR_INT("invalid plottype", __func__, 1);
998
999 lept_mkdir("lept/comp");
1000
1001 if (comptype == L_COMPARE_SUBTRACT)
1002 pixt = pixSubtractGray(NULL, pix1, pix2);
1003 else /* comptype == L_COMPARE_ABS_DIFF) */
1004 pixt = pixAbsDifference(pix1, pix2);
1005
1006 pixZero(pixt, &same);
1007 if (same)
1008 L_INFO("Images are pixel-wise identical\n", __func__);
1009 if (psame) *psame = same;
1010
1011 if (pdiff)
1012 pixGetAverageMasked(pixt, NULL, 0, 0, 1, L_MEAN_ABSVAL, pdiff);
1013
1014 /* Don't bother to plot if the images are the same */
1015 if (plottype && !same) {
1016 L_INFO("Images differ: output plots will be generated\n", __func__);
1017 na = pixGetGrayHistogram(pixt, 1);
1018 numaGetNonzeroRange(na, TINY, &first, &last);
1019 nac = numaClipToInterval(na, 0, last);
1020 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_gray%d", index);
1021 gplot = gplotCreate(buf, plottype,
1022 "Pixel Difference Histogram", "diff val",
1023 "number of pixels");
1024 gplotAddPlot(gplot, NULL, nac, GPLOT_LINES, "gray");
1025 gplotMakeOutput(gplot);
1026 gplotDestroy(&gplot);
1027 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_gray%d.png",
1028 index++);
1029 l_fileDisplay(buf, 100, 100, 1.0);
1030 numaDestroy(&na);
1031 numaDestroy(&nac);
1032 }
1033
1034 if (ppixdiff)
1035 *ppixdiff = pixCopy(NULL, pixt);
1036
1037 if (prmsdiff) {
1038 if (comptype == L_COMPARE_SUBTRACT) { /* wrong type for rms diff */
1039 pixDestroy(&pixt);
1040 pixt = pixAbsDifference(pix1, pix2);
1041 }
1042 pixGetAverageMasked(pixt, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, prmsdiff);
1043 }
1044
1045 pixDestroy(&pixt);
1046 return 0;
1047}
1048
1049
1070l_ok
1071pixCompareRGB(PIX *pix1,
1072 PIX *pix2,
1073 l_int32 comptype,
1074 l_int32 plottype,
1075 l_int32 *psame,
1076 l_float32 *pdiff,
1077 l_float32 *prmsdiff,
1078 PIX **ppixdiff)
1079{
1080char buf[64];
1081static l_atomic index = 0;
1082l_int32 rsame, gsame, bsame, same, first, rlast, glast, blast, last;
1083l_float32 rdiff, gdiff, bdiff;
1084GPLOT *gplot;
1085NUMA *nar, *nag, *nab, *narc, *nagc, *nabc;
1086PIX *pixr1, *pixr2, *pixg1, *pixg2, *pixb1, *pixb2;
1087PIX *pixr, *pixg, *pixb;
1088
1089 if (psame) *psame = 0;
1090 if (pdiff) *pdiff = 0.0;
1091 if (prmsdiff) *prmsdiff = 0.0;
1092 if (ppixdiff) *ppixdiff = NULL;
1093 if (!pix1 || pixGetDepth(pix1) != 32)
1094 return ERROR_INT("pix1 not defined or not 32 bpp", __func__, 1);
1095 if (!pix2 || pixGetDepth(pix2) != 32)
1096 return ERROR_INT("pix2 not defined or not ew bpp", __func__, 1);
1097 if (comptype != L_COMPARE_SUBTRACT && comptype != L_COMPARE_ABS_DIFF)
1098 return ERROR_INT("invalid comptype", __func__, 1);
1099 if (plottype < 0 || plottype >= NUM_GPLOT_OUTPUTS)
1100 return ERROR_INT("invalid plottype", __func__, 1);
1101
1102 lept_mkdir("lept/comp");
1103
1104 pixr1 = pixGetRGBComponent(pix1, COLOR_RED);
1105 pixr2 = pixGetRGBComponent(pix2, COLOR_RED);
1106 pixg1 = pixGetRGBComponent(pix1, COLOR_GREEN);
1107 pixg2 = pixGetRGBComponent(pix2, COLOR_GREEN);
1108 pixb1 = pixGetRGBComponent(pix1, COLOR_BLUE);
1109 pixb2 = pixGetRGBComponent(pix2, COLOR_BLUE);
1110 if (comptype == L_COMPARE_SUBTRACT) {
1111 pixr = pixSubtractGray(NULL, pixr1, pixr2);
1112 pixg = pixSubtractGray(NULL, pixg1, pixg2);
1113 pixb = pixSubtractGray(NULL, pixb1, pixb2);
1114 } else { /* comptype == L_COMPARE_ABS_DIFF) */
1115 pixr = pixAbsDifference(pixr1, pixr2);
1116 pixg = pixAbsDifference(pixg1, pixg2);
1117 pixb = pixAbsDifference(pixb1, pixb2);
1118 }
1119
1120 pixZero(pixr, &rsame);
1121 pixZero(pixg, &gsame);
1122 pixZero(pixb, &bsame);
1123 same = rsame && gsame && bsame;
1124 if (same)
1125 L_INFO("Images are pixel-wise identical\n", __func__);
1126 if (psame) *psame = same;
1127
1128 if (pdiff) {
1129 pixGetAverageMasked(pixr, NULL, 0, 0, 1, L_MEAN_ABSVAL, &rdiff);
1130 pixGetAverageMasked(pixg, NULL, 0, 0, 1, L_MEAN_ABSVAL, &gdiff);
1131 pixGetAverageMasked(pixb, NULL, 0, 0, 1, L_MEAN_ABSVAL, &bdiff);
1132 *pdiff = (rdiff + gdiff + bdiff) / 3.0;
1133 }
1134
1135 /* Don't bother to plot if the images are the same */
1136 if (plottype && !same) {
1137 L_INFO("Images differ: output plots will be generated\n", __func__);
1138 nar = pixGetGrayHistogram(pixr, 1);
1139 nag = pixGetGrayHistogram(pixg, 1);
1140 nab = pixGetGrayHistogram(pixb, 1);
1141 numaGetNonzeroRange(nar, TINY, &first, &rlast);
1142 numaGetNonzeroRange(nag, TINY, &first, &glast);
1143 numaGetNonzeroRange(nab, TINY, &first, &blast);
1144 last = L_MAX(rlast, glast);
1145 last = L_MAX(last, blast);
1146 narc = numaClipToInterval(nar, 0, last);
1147 nagc = numaClipToInterval(nag, 0, last);
1148 nabc = numaClipToInterval(nab, 0, last);
1149 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_rgb%d", index);
1150 gplot = gplotCreate(buf, plottype,
1151 "Pixel Difference Histogram", "diff val",
1152 "number of pixels");
1153 gplotAddPlot(gplot, NULL, narc, GPLOT_LINES, "red");
1154 gplotAddPlot(gplot, NULL, nagc, GPLOT_LINES, "green");
1155 gplotAddPlot(gplot, NULL, nabc, GPLOT_LINES, "blue");
1156 gplotMakeOutput(gplot);
1157 gplotDestroy(&gplot);
1158 snprintf(buf, sizeof(buf), "/tmp/lept/comp/compare_rgb%d.png",
1159 index++);
1160 l_fileDisplay(buf, 100, 100, 1.0);
1161 numaDestroy(&nar);
1162 numaDestroy(&nag);
1163 numaDestroy(&nab);
1164 numaDestroy(&narc);
1165 numaDestroy(&nagc);
1166 numaDestroy(&nabc);
1167 }
1168
1169 if (ppixdiff)
1170 *ppixdiff = pixCreateRGBImage(pixr, pixg, pixb);
1171
1172 if (prmsdiff) {
1173 if (comptype == L_COMPARE_SUBTRACT) {
1174 pixDestroy(&pixr);
1175 pixDestroy(&pixg);
1176 pixDestroy(&pixb);
1177 pixr = pixAbsDifference(pixr1, pixr2);
1178 pixg = pixAbsDifference(pixg1, pixg2);
1179 pixb = pixAbsDifference(pixb1, pixb2);
1180 }
1181 pixGetAverageMasked(pixr, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, &rdiff);
1182 pixGetAverageMasked(pixg, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, &gdiff);
1183 pixGetAverageMasked(pixb, NULL, 0, 0, 1, L_ROOT_MEAN_SQUARE, &bdiff);
1184 *prmsdiff = (rdiff + gdiff + bdiff) / 3.0;
1185 }
1186
1187 pixDestroy(&pixr1);
1188 pixDestroy(&pixr2);
1189 pixDestroy(&pixg1);
1190 pixDestroy(&pixg2);
1191 pixDestroy(&pixb1);
1192 pixDestroy(&pixb2);
1193 pixDestroy(&pixr);
1194 pixDestroy(&pixg);
1195 pixDestroy(&pixb);
1196 return 0;
1197}
1198
1199
1224l_ok
1225pixCompareTiled(PIX *pix1,
1226 PIX *pix2,
1227 l_int32 sx,
1228 l_int32 sy,
1229 l_int32 type,
1230 PIX **ppixdiff)
1231{
1232l_int32 d1, d2, w, h;
1233PIX *pixt, *pixr, *pixg, *pixb;
1234PIX *pixrdiff, *pixgdiff, *pixbdiff;
1235PIXACC *pixacc;
1236
1237 if (!ppixdiff)
1238 return ERROR_INT("&pixdiff not defined", __func__, 1);
1239 *ppixdiff = NULL;
1240 if (!pix1)
1241 return ERROR_INT("pix1 not defined", __func__, 1);
1242 if (!pix2)
1243 return ERROR_INT("pix2 not defined", __func__, 1);
1244 d1 = pixGetDepth(pix1);
1245 d2 = pixGetDepth(pix2);
1246 if (d1 != d2)
1247 return ERROR_INT("depths not equal", __func__, 1);
1248 if (d1 != 8 && d1 != 32)
1249 return ERROR_INT("pix1 not 8 or 32 bpp", __func__, 1);
1250 if (d2 != 8 && d2 != 32)
1251 return ERROR_INT("pix2 not 8 or 32 bpp", __func__, 1);
1252 if (sx < 2 || sy < 2)
1253 return ERROR_INT("sx and sy not both > 1", __func__, 1);
1254 if (type != L_MEAN_ABSVAL && type != L_ROOT_MEAN_SQUARE)
1255 return ERROR_INT("invalid type", __func__, 1);
1256
1257 pixt = pixAbsDifference(pix1, pix2);
1258 if (d1 == 8) {
1259 *ppixdiff = pixGetAverageTiled(pixt, sx, sy, type);
1260 } else { /* d1 == 32 */
1261 pixr = pixGetRGBComponent(pixt, COLOR_RED);
1262 pixg = pixGetRGBComponent(pixt, COLOR_GREEN);
1263 pixb = pixGetRGBComponent(pixt, COLOR_BLUE);
1264 pixrdiff = pixGetAverageTiled(pixr, sx, sy, type);
1265 pixgdiff = pixGetAverageTiled(pixg, sx, sy, type);
1266 pixbdiff = pixGetAverageTiled(pixb, sx, sy, type);
1267 pixGetDimensions(pixrdiff, &w, &h, NULL);
1268 pixacc = pixaccCreate(w, h, 0);
1269 pixaccAdd(pixacc, pixrdiff);
1270 pixaccAdd(pixacc, pixgdiff);
1271 pixaccAdd(pixacc, pixbdiff);
1272 pixaccMultConst(pixacc, 1.f / 3.f);
1273 *ppixdiff = pixaccFinal(pixacc, 8);
1274 pixDestroy(&pixr);
1275 pixDestroy(&pixg);
1276 pixDestroy(&pixb);
1277 pixDestroy(&pixrdiff);
1278 pixDestroy(&pixgdiff);
1279 pixDestroy(&pixbdiff);
1280 pixaccDestroy(&pixacc);
1281 }
1282 pixDestroy(&pixt);
1283 return 0;
1284}
1285
1286
1287/*------------------------------------------------------------------*
1288 * Other measures of the difference of two images *
1289 *------------------------------------------------------------------*/
1316NUMA *
1317pixCompareRankDifference(PIX *pix1,
1318 PIX *pix2,
1319 l_int32 factor)
1320{
1321l_int32 i;
1322l_float32 *array1, *array2;
1323NUMA *nah, *nan, *nad;
1324
1325 if (!pix1)
1326 return (NUMA *)ERROR_PTR("pix1 not defined", __func__, NULL);
1327 if (!pix2)
1328 return (NUMA *)ERROR_PTR("pix2 not defined", __func__, NULL);
1329
1330 if ((nah = pixGetDifferenceHistogram(pix1, pix2, factor)) == NULL)
1331 return (NUMA *)ERROR_PTR("na not made", __func__, NULL);
1332
1333 nan = numaNormalizeHistogram(nah, 1.0);
1334 array1 = numaGetFArray(nan, L_NOCOPY);
1335
1336 nad = numaCreate(256);
1337 numaSetCount(nad, 256); /* all initialized to 0.0 */
1338 array2 = numaGetFArray(nad, L_NOCOPY);
1339
1340 /* Do rank accumulation on normalized histo of diffs */
1341 array2[0] = 1.0;
1342 for (i = 1; i < 256; i++)
1343 array2[i] = array2[i - 1] - array1[i - 1];
1344
1345 numaDestroy(&nah);
1346 numaDestroy(&nan);
1347 return nad;
1348}
1349
1350
1399l_ok
1400pixTestForSimilarity(PIX *pix1,
1401 PIX *pix2,
1402 l_int32 factor,
1403 l_int32 mindiff,
1404 l_float32 maxfract,
1405 l_float32 maxave,
1406 l_int32 *psimilar,
1407 l_int32 details)
1408{
1409l_float32 fractdiff, avediff;
1410
1411 if (!psimilar)
1412 return ERROR_INT("&similar not defined", __func__, 1);
1413 *psimilar = 0;
1414 if (!pix1)
1415 return ERROR_INT("pix1 not defined", __func__, 1);
1416 if (!pix2)
1417 return ERROR_INT("pix2 not defined", __func__, 1);
1418 if (pixSizesEqual(pix1, pix2) == 0)
1419 return ERROR_INT("pix sizes not equal", __func__, 1);
1420 if (mindiff <= 0)
1421 return ERROR_INT("mindiff must be > 0", __func__, 1);
1422
1423 if (pixGetDifferenceStats(pix1, pix2, factor, mindiff,
1424 &fractdiff, &avediff, details))
1425 return ERROR_INT("diff stats not found", __func__, 1);
1426
1427 if (maxave <= 0.0) maxave = 256.0;
1428 if (fractdiff <= maxfract && avediff <= maxave)
1429 *psimilar = 1;
1430 return 0;
1431}
1432
1433
1476l_ok
1477pixGetDifferenceStats(PIX *pix1,
1478 PIX *pix2,
1479 l_int32 factor,
1480 l_int32 mindiff,
1481 l_float32 *pfractdiff,
1482 l_float32 *pavediff,
1483 l_int32 details)
1484{
1485l_int32 i, first, last, diff;
1486l_float32 fract, ave;
1487l_float32 *array;
1488NUMA *nah, *nan, *nac;
1489
1490 if (pfractdiff) *pfractdiff = 0.0;
1491 if (pavediff) *pavediff = 0.0;
1492 if (!pfractdiff)
1493 return ERROR_INT("&fractdiff not defined", __func__, 1);
1494 if (!pavediff)
1495 return ERROR_INT("&avediff not defined", __func__, 1);
1496 if (!pix1)
1497 return ERROR_INT("pix1 not defined", __func__, 1);
1498 if (!pix2)
1499 return ERROR_INT("pix2 not defined", __func__, 1);
1500 if (mindiff <= 0)
1501 return ERROR_INT("mindiff must be > 0", __func__, 1);
1502
1503 if ((nah = pixGetDifferenceHistogram(pix1, pix2, factor)) == NULL)
1504 return ERROR_INT("na not made", __func__, 1);
1505
1506 if ((nan = numaNormalizeHistogram(nah, 1.0)) == NULL) {
1507 numaDestroy(&nah);
1508 return ERROR_INT("nan not made", __func__, 1);
1509 }
1510 array = numaGetFArray(nan, L_NOCOPY);
1511
1512 if (details) {
1513 lept_mkdir("lept/comp");
1514 numaGetNonzeroRange(nan, 0.0, &first, &last);
1515 nac = numaClipToInterval(nan, first, last);
1516 gplotSimple1(nac, GPLOT_PNG, "/tmp/lept/comp/histo",
1517 "Difference histogram");
1518 l_fileDisplay("/tmp/lept/comp/histo.png", 500, 0, 1.0);
1519 lept_stderr("\nNonzero values in normalized histogram:");
1520 numaWriteStderr(nac);
1521 numaDestroy(&nac);
1522 lept_stderr(" Mindiff fractdiff avediff\n");
1523 lept_stderr(" -----------------------------------\n");
1524 for (diff = 1; diff < L_MIN(2 * mindiff, last); diff++) {
1525 fract = 0.0;
1526 ave = 0.0;
1527 for (i = diff; i <= last; i++) {
1528 fract += array[i];
1529 ave += (l_float32)i * array[i];
1530 }
1531 ave = (fract == 0.0) ? 0.0 : ave / fract;
1532 ave -= diff;
1533 lept_stderr("%5d %7.4f %7.4f\n",
1534 diff, fract, ave);
1535 }
1536 lept_stderr(" -----------------------------------\n");
1537 }
1538
1539 fract = 0.0;
1540 ave = 0.0;
1541 for (i = mindiff; i < 256; i++) {
1542 fract += array[i];
1543 ave += (l_float32)i * array[i];
1544 }
1545 ave = (fract == 0.0) ? 0.0 : ave / fract;
1546 ave -= mindiff;
1547
1548 *pfractdiff = fract;
1549 *pavediff = ave;
1550
1551 numaDestroy(&nah);
1552 numaDestroy(&nan);
1553 return 0;
1554}
1555
1556
1576NUMA *
1577pixGetDifferenceHistogram(PIX *pix1,
1578 PIX *pix2,
1579 l_int32 factor)
1580{
1581l_int32 w1, h1, d1, w2, h2, d2, w, h, wpl1, wpl2;
1582l_int32 i, j, val, val1, val2;
1583l_int32 rval1, rval2, gval1, gval2, bval1, bval2;
1584l_int32 rdiff, gdiff, bdiff, maxdiff;
1585l_uint32 *data1, *data2, *line1, *line2;
1586l_float32 *array;
1587NUMA *na;
1588PIX *pixt1, *pixt2;
1589
1590 if (!pix1)
1591 return (NUMA *)ERROR_PTR("pix1 not defined", __func__, NULL);
1592 if (!pix2)
1593 return (NUMA *)ERROR_PTR("pix2 not defined", __func__, NULL);
1594 d1 = pixGetDepth(pix1);
1595 d2 = pixGetDepth(pix2);
1596 if (d1 == 16 || d2 == 16)
1597 return (NUMA *)ERROR_PTR("d == 16 not supported", __func__, NULL);
1598 if (d1 < 8 && !pixGetColormap(pix1))
1599 return (NUMA *)ERROR_PTR("pix1 depth < 8 bpp and not cmapped",
1600 __func__, NULL);
1601 if (d2 < 8 && !pixGetColormap(pix2))
1602 return (NUMA *)ERROR_PTR("pix2 depth < 8 bpp and not cmapped",
1603 __func__, NULL);
1604 pixt1 = pixRemoveColormap(pix1, REMOVE_CMAP_BASED_ON_SRC);
1605 pixt2 = pixRemoveColormap(pix2, REMOVE_CMAP_BASED_ON_SRC);
1606 pixGetDimensions(pixt1, &w1, &h1, &d1);
1607 pixGetDimensions(pixt2, &w2, &h2, &d2);
1608 if (d1 != d2) {
1609 pixDestroy(&pixt1);
1610 pixDestroy(&pixt2);
1611 return (NUMA *)ERROR_PTR("pix depths not equal", __func__, NULL);
1612 }
1613 if (factor < 1) factor = 1;
1614
1615 na = numaCreate(256);
1616 numaSetCount(na, 256); /* all initialized to 0.0 */
1617 array = numaGetFArray(na, L_NOCOPY);
1618 w = L_MIN(w1, w2);
1619 h = L_MIN(h1, h2);
1620 data1 = pixGetData(pixt1);
1621 data2 = pixGetData(pixt2);
1622 wpl1 = pixGetWpl(pixt1);
1623 wpl2 = pixGetWpl(pixt2);
1624 if (d1 == 8) {
1625 for (i = 0; i < h; i += factor) {
1626 line1 = data1 + i * wpl1;
1627 line2 = data2 + i * wpl2;
1628 for (j = 0; j < w; j += factor) {
1629 val1 = GET_DATA_BYTE(line1, j);
1630 val2 = GET_DATA_BYTE(line2, j);
1631 val = L_ABS(val1 - val2);
1632 array[val]++;
1633 }
1634 }
1635 } else { /* d1 == 32 */
1636 for (i = 0; i < h; i += factor) {
1637 line1 = data1 + i * wpl1;
1638 line2 = data2 + i * wpl2;
1639 for (j = 0; j < w; j += factor) {
1640 extractRGBValues(line1[j], &rval1, &gval1, &bval1);
1641 extractRGBValues(line2[j], &rval2, &gval2, &bval2);
1642 rdiff = L_ABS(rval1 - rval2);
1643 gdiff = L_ABS(gval1 - gval2);
1644 bdiff = L_ABS(bval1 - bval2);
1645 maxdiff = L_MAX(rdiff, gdiff);
1646 maxdiff = L_MAX(maxdiff, bdiff);
1647 array[maxdiff]++;
1648 }
1649 }
1650 }
1651
1652 pixDestroy(&pixt1);
1653 pixDestroy(&pixt2);
1654 return na;
1655}
1656
1657
1705l_ok
1706pixGetPerceptualDiff(PIX *pixs1,
1707 PIX *pixs2,
1708 l_int32 sampling,
1709 l_int32 dilation,
1710 l_int32 mindiff,
1711 l_float32 *pfract,
1712 PIX **ppixdiff1,
1713 PIX **ppixdiff2)
1714{
1715l_int32 d1, d2, w, h, count;
1716PIX *pix1, *pix2, *pix3, *pix4, *pix5, *pix6, *pix7, *pix8, *pix9;
1717PIX *pix10, *pix11;
1718
1719 if (ppixdiff1) *ppixdiff1 = NULL;
1720 if (ppixdiff2) *ppixdiff2 = NULL;
1721 if (!pfract)
1722 return ERROR_INT("&fract not defined", __func__, 1);
1723 *pfract = 1.0; /* init to completely different */
1724 if ((dilation & 1) == 0)
1725 return ERROR_INT("dilation must be odd", __func__, 1);
1726 if (!pixs1)
1727 return ERROR_INT("pixs1 not defined", __func__, 1);
1728 if (!pixs2)
1729 return ERROR_INT("pixs2 not defined", __func__, 1);
1730 d1 = pixGetDepth(pixs1);
1731 d2 = pixGetDepth(pixs2);
1732 if (!pixGetColormap(pixs1) && d1 < 8)
1733 return ERROR_INT("pixs1 not cmapped and < 8 bpp", __func__, 1);
1734 if (!pixGetColormap(pixs2) && d2 < 8)
1735 return ERROR_INT("pixs2 not cmapped and < 8 bpp", __func__, 1);
1736
1737 /* Integer downsample if requested */
1738 if (sampling > 1) {
1739 pix1 = pixScaleByIntSampling(pixs1, sampling);
1740 pix2 = pixScaleByIntSampling(pixs2, sampling);
1741 } else {
1742 pix1 = pixClone(pixs1);
1743 pix2 = pixClone(pixs2);
1744 }
1745
1746 /* Remove colormaps */
1747 if (pixGetColormap(pix1)) {
1748 pix3 = pixRemoveColormap(pix1, REMOVE_CMAP_BASED_ON_SRC);
1749 d1 = pixGetDepth(pix3);
1750 } else {
1751 pix3 = pixClone(pix1);
1752 }
1753 if (pixGetColormap(pix2)) {
1754 pix4 = pixRemoveColormap(pix2, REMOVE_CMAP_BASED_ON_SRC);
1755 d2 = pixGetDepth(pix4);
1756 } else {
1757 pix4 = pixClone(pix2);
1758 }
1759 pixDestroy(&pix1);
1760 pixDestroy(&pix2);
1761 if (d1 != d2 || (d1 != 8 && d1 != 32)) {
1762 pixDestroy(&pix3);
1763 pixDestroy(&pix4);
1764 L_INFO("depths unequal or not in {8,32}: d1 = %d, d2 = %d\n",
1765 __func__, d1, d2);
1766 return 1;
1767 }
1768
1769 /* In each direction, do a small dilation and subtract the dilated
1770 * image from the other image to get a one-sided difference.
1771 * Then take the max of the differences for each direction
1772 * and clipping each component to 255 if necessary. Note that
1773 * for RGB images, the dilations and max selection are done
1774 * component-wise, and the conversion to grayscale also uses the
1775 * maximum component. The resulting grayscale images are
1776 * thresholded using %mindiff. */
1777 if (d1 == 8) {
1778 pix5 = pixDilateGray(pix3, dilation, dilation);
1779 pixCompareGray(pix4, pix5, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1780 &pix7);
1781 pix6 = pixDilateGray(pix4, dilation, dilation);
1782 pixCompareGray(pix3, pix6, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1783 &pix8);
1784 pix9 = pixMinOrMax(NULL, pix7, pix8, L_CHOOSE_MAX);
1785 pix10 = pixThresholdToBinary(pix9, mindiff);
1786 pixInvert(pix10, pix10);
1787 pixCountPixels(pix10, &count, NULL);
1788 pixGetDimensions(pix10, &w, &h, NULL);
1789 *pfract = (w <= 0 || h <= 0) ? 0.0 :
1790 (l_float32)count / (l_float32)(w * h);
1791 pixDestroy(&pix5);
1792 pixDestroy(&pix6);
1793 pixDestroy(&pix7);
1794 pixDestroy(&pix8);
1795 if (ppixdiff1)
1796 *ppixdiff1 = pix9;
1797 else
1798 pixDestroy(&pix9);
1799 if (ppixdiff2)
1800 *ppixdiff2 = pix10;
1801 else
1802 pixDestroy(&pix10);
1803 } else { /* d1 == 32 */
1804 pix5 = pixColorMorph(pix3, L_MORPH_DILATE, dilation, dilation);
1805 pixCompareRGB(pix4, pix5, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1806 &pix7);
1807 pix6 = pixColorMorph(pix4, L_MORPH_DILATE, dilation, dilation);
1808 pixCompareRGB(pix3, pix6, L_COMPARE_SUBTRACT, 0, NULL, NULL, NULL,
1809 &pix8);
1810 pix9 = pixMinOrMax(NULL, pix7, pix8, L_CHOOSE_MAX);
1811 pix10 = pixConvertRGBToGrayMinMax(pix9, L_CHOOSE_MAX);
1812 pix11 = pixThresholdToBinary(pix10, mindiff);
1813 pixInvert(pix11, pix11);
1814 pixCountPixels(pix11, &count, NULL);
1815 pixGetDimensions(pix11, &w, &h, NULL);
1816 *pfract = (w <= 0 || h <= 0) ? 0.0 :
1817 (l_float32)count / (l_float32)(w * h);
1818 pixDestroy(&pix5);
1819 pixDestroy(&pix6);
1820 pixDestroy(&pix7);
1821 pixDestroy(&pix8);
1822 pixDestroy(&pix10);
1823 if (ppixdiff1)
1824 *ppixdiff1 = pix9;
1825 else
1826 pixDestroy(&pix9);
1827 if (ppixdiff2)
1828 *ppixdiff2 = pix11;
1829 else
1830 pixDestroy(&pix11);
1831
1832 }
1833 pixDestroy(&pix3);
1834 pixDestroy(&pix4);
1835 return 0;
1836}
1837
1838
1870l_ok
1871pixGetPSNR(PIX *pix1,
1872 PIX *pix2,
1873 l_int32 factor,
1874 l_float32 *ppsnr)
1875{
1876l_int32 same, i, j, w, h, d, wpl1, wpl2, v1, v2, r1, g1, b1, r2, g2, b2;
1877l_uint32 *data1, *data2, *line1, *line2;
1878l_float32 mse; /* mean squared error */
1879
1880 if (!ppsnr)
1881 return ERROR_INT("&psnr not defined", __func__, 1);
1882 *ppsnr = 0.0;
1883 if (!pix1 || !pix2)
1884 return ERROR_INT("empty input pix", __func__, 1);
1885 if (!pixSizesEqual(pix1, pix2))
1886 return ERROR_INT("pix sizes unequal", __func__, 1);
1887 if (pixGetColormap(pix1))
1888 return ERROR_INT("pix1 has colormap", __func__, 1);
1889 if (pixGetColormap(pix2))
1890 return ERROR_INT("pix2 has colormap", __func__, 1);
1891 pixGetDimensions(pix1, &w, &h, &d);
1892 if (d != 8 && d != 32)
1893 return ERROR_INT("pix not 8 or 32 bpp", __func__, 1);
1894 if (factor < 1)
1895 return ERROR_INT("invalid sampling factor", __func__, 1);
1896
1897 pixEqual(pix1, pix2, &same);
1898 if (same) {
1899 *ppsnr = 1000.0; /* crazy big exponent */
1900 return 0;
1901 }
1902
1903 data1 = pixGetData(pix1);
1904 data2 = pixGetData(pix2);
1905 wpl1 = pixGetWpl(pix1);
1906 wpl2 = pixGetWpl(pix2);
1907 mse = 0.0;
1908 if (d == 8) {
1909 for (i = 0; i < h; i += factor) {
1910 line1 = data1 + i * wpl1;
1911 line2 = data2 + i * wpl2;
1912 for (j = 0; j < w; j += factor) {
1913 v1 = GET_DATA_BYTE(line1, j);
1914 v2 = GET_DATA_BYTE(line2, j);
1915 mse += (l_float32)(v1 - v2) * (v1 - v2);
1916 }
1917 }
1918 } else { /* d == 32 */
1919 for (i = 0; i < h; i += factor) {
1920 line1 = data1 + i * wpl1;
1921 line2 = data2 + i * wpl2;
1922 for (j = 0; j < w; j += factor) {
1923 extractRGBValues(line1[j], &r1, &g1, &b1);
1924 extractRGBValues(line2[j], &r2, &g2, &b2);
1925 mse += ((l_float32)(r1 - r2) * (r1 - r2) +
1926 (g1 - g2) * (g1 - g2) +
1927 (b1 - b2) * (b1 - b2)) / 3.0;
1928 }
1929 }
1930 }
1931 mse = mse / ((l_float32)(w) * h);
1932
1933 *ppsnr = -4.3429448 * log(mse / (255 * 255));
1934 return 0;
1935}
1936
1937
1938/*------------------------------------------------------------------*
1939 * Comparison of photo regions by histogram *
1940 *------------------------------------------------------------------*/
1993l_ok
1994pixaComparePhotoRegionsByHisto(PIXA *pixa,
1995 l_float32 minratio,
1996 l_float32 textthresh,
1997 l_int32 factor,
1998 l_int32 n,
1999 l_float32 simthresh,
2000 NUMA **pnai,
2001 l_float32 **pscores,
2002 PIX **ppixd,
2003 l_int32 debug)
2004{
2005char *text;
2006l_int32 i, j, nim, w, h, w1, h1, w2, h2, ival, index, classid;
2007l_float32 score;
2008l_float32 *scores;
2009NUMA *nai, *naw, *nah;
2010NUMAA *naa;
2011NUMAA **n3a; /* array of naa */
2012PIX *pix;
2013
2014 if (pscores) *pscores = NULL;
2015 if (ppixd) *ppixd = NULL;
2016 if (!pnai)
2017 return ERROR_INT("&na not defined", __func__, 1);
2018 *pnai = NULL;
2019 if (!pixa)
2020 return ERROR_INT("pixa not defined", __func__, 1);
2021 if (minratio < 0.0 || minratio > 1.0)
2022 return ERROR_INT("minratio not in [0.0 ... 1.0]", __func__, 1);
2023 if (textthresh <= 0.0) textthresh = 1.3f;
2024 if (factor < 1)
2025 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2026 if (n < 1 || n > 7) {
2027 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2028 n = 4;
2029 }
2030 if (simthresh <= 0.0) simthresh = 0.25;
2031 if (simthresh > 1.0)
2032 return ERROR_INT("simthresh invalid; should be near 0.25", __func__, 1);
2033
2034 /* Prepare the histograms */
2035 nim = pixaGetCount(pixa);
2036 if ((n3a = (NUMAA **)LEPT_CALLOC(nim, sizeof(NUMAA *))) == NULL)
2037 return ERROR_INT("calloc fail for n3a", __func__, 1);
2038 naw = numaCreate(0);
2039 nah = numaCreate(0);
2040 for (i = 0; i < nim; i++) {
2041 pix = pixaGetPix(pixa, i, L_CLONE);
2042 text = pixGetText(pix);
2043 pixSetResolution(pix, 150, 150);
2044 index = (debug) ? i : 0;
2045 pixGenPhotoHistos(pix, NULL, factor, textthresh, n,
2046 &naa, &w, &h, index);
2047 n3a[i] = naa;
2048 numaAddNumber(naw, w);
2049 numaAddNumber(nah, h);
2050 if (naa)
2051 lept_stderr("Image %s is photo\n", text);
2052 else
2053 lept_stderr("Image %s is NOT photo\n", text);
2054 pixDestroy(&pix);
2055 }
2056
2057 /* Do the comparisons. We are making a set of classes, where
2058 * all similar images are placed in the same class. There are
2059 * 'nim' input images. The classes are labeled by 'classid' (all
2060 * similar images get the same 'classid' value), and 'nai' maps
2061 * the classid of the image in the input array to the classid
2062 * of the similarity class. */
2063 if ((scores =
2064 (l_float32 *)LEPT_CALLOC((size_t)nim * nim, sizeof(l_float32)))
2065 == NULL) {
2066 L_ERROR("calloc fail for scores\n", __func__);
2067 goto cleanup;
2068 }
2069 nai = numaMakeConstant(-1, nim); /* classid array */
2070 for (i = 0, classid = 0; i < nim; i++) {
2071 scores[nim * i + i] = 1.0;
2072 numaGetIValue(nai, i, &ival);
2073 if (ival != -1) /* already set */
2074 continue;
2075 numaSetValue(nai, i, classid);
2076 if (n3a[i] == NULL) { /* not a photo */
2077 classid++;
2078 continue;
2079 }
2080 numaGetIValue(naw, i, &w1);
2081 numaGetIValue(nah, i, &h1);
2082 for (j = i + 1; j < nim; j++) {
2083 numaGetIValue(nai, j, &ival);
2084 if (ival != -1) /* already set */
2085 continue;
2086 if (n3a[j] == NULL) /* not a photo */
2087 continue;
2088 numaGetIValue(naw, j, &w2);
2089 numaGetIValue(nah, j, &h2);
2090 compareTilesByHisto(n3a[i], n3a[j], minratio, w1, h1, w2, h2,
2091 &score, NULL);
2092 scores[nim * i + j] = score;
2093 scores[nim * j + i] = score; /* the score array is symmetric */
2094/* lept_stderr("score = %5.3f\n", score); */
2095 if (score > simthresh) {
2096 numaSetValue(nai, j, classid);
2097 lept_stderr(
2098 "Setting %d similar to %d, in class %d; score %5.3f\n",
2099 j, i, classid, score);
2100 }
2101 }
2102 classid++;
2103 }
2104 *pnai = nai;
2105
2106 /* Debug: optionally save and display the score array.
2107 * All images that are photos are represented by a point on
2108 * the diagonal. Other images in the same similarity class
2109 * are on the same horizontal raster line to the right.
2110 * The array has been symmetrized, so images in the same
2111 * same similarity class also appear on the same column below. */
2112 if (pscores) {
2113 l_int32 wpl, fact;
2114 l_uint32 *line, *data;
2115 PIX *pix2, *pix3;
2116 pix2 = pixCreate(nim, nim, 8);
2117 data = pixGetData(pix2);
2118 wpl = pixGetWpl(pix2);
2119 for (i = 0; i < nim; i++) {
2120 line = data + i * wpl;
2121 for (j = 0; j < nim; j++) {
2122 SET_DATA_BYTE(line, j,
2123 L_MIN(255, 4.0 * 255 * scores[nim * i + j]));
2124 }
2125 }
2126 fact = L_MAX(2, 1000 / nim);
2127 pix3 = pixExpandReplicate(pix2, fact);
2128 lept_stderr("Writing to /tmp/lept/comp/scorearray.png\n");
2129 lept_mkdir("lept/comp");
2130 pixWrite("/tmp/lept/comp/scorearray.png", pix3, IFF_PNG);
2131 pixDestroy(&pix2);
2132 pixDestroy(&pix3);
2133 *pscores = scores;
2134 } else {
2135 LEPT_FREE(scores);
2136 }
2137
2138 /* Debug: optionally display and save the image comparisons.
2139 * Image similarity classes are displayed by column; similar
2140 * images are displayed in the same column. */
2141 if (ppixd)
2142 *ppixd = pixaDisplayTiledByIndex(pixa, nai, 200, 20, 2, 6, 0x0000ff00);
2143
2144cleanup:
2145 numaDestroy(&naw);
2146 numaDestroy(&nah);
2147 for (i = 0; i < nim; i++)
2148 numaaDestroy(&n3a[i]);
2149 LEPT_FREE(n3a);
2150 return 0;
2151}
2152
2153
2213l_ok
2214pixComparePhotoRegionsByHisto(PIX *pix1,
2215 PIX *pix2,
2216 BOX *box1,
2217 BOX *box2,
2218 l_float32 minratio,
2219 l_int32 factor,
2220 l_int32 n,
2221 l_float32 *pscore,
2222 l_int32 debugflag)
2223{
2224l_int32 w1, h1, w2, h2, w1c, h1c, w2c, h2c, debugindex;
2225l_float32 wratio, hratio;
2226NUMAA *naa1, *naa2;
2227PIX *pix3, *pix4;
2228PIXA *pixa;
2229
2230 if (!pscore)
2231 return ERROR_INT("&score not defined", __func__, 1);
2232 *pscore = 0.0;
2233 if (!pix1 || !pix2)
2234 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
2235 if (minratio < 0.5 || minratio > 1.0)
2236 return ERROR_INT("minratio not in [0.5 ... 1.0]", __func__, 1);
2237 if (factor < 1)
2238 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2239 if (n < 1 || n > 7) {
2240 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2241 n = 4;
2242 }
2243
2244 debugindex = 0;
2245 if (debugflag) {
2246 lept_mkdir("lept/comp");
2247 debugindex = 666; /* arbitrary number used for naming output */
2248 }
2249
2250 /* Initial filter by size */
2251 if (box1)
2252 boxGetGeometry(box1, NULL, NULL, &w1, &h1);
2253 else
2254 pixGetDimensions(pix1, &w1, &h1, NULL);
2255 if (box2)
2256 boxGetGeometry(box2, NULL, NULL, &w2, &h2);
2257 else
2258 pixGetDimensions(pix1, &w2, &h2, NULL);
2259 wratio = (w1 < w2) ? (l_float32)w1 / (l_float32)w2 :
2260 (l_float32)w2 / (l_float32)w1;
2261 hratio = (h1 < h2) ? (l_float32)h1 / (l_float32)h2 :
2262 (l_float32)h2 / (l_float32)h1;
2263 if (wratio < minratio || hratio < minratio)
2264 return 0;
2265
2266 /* Initial crop, if necessary, and make histos */
2267 if (box1)
2268 pix3 = pixClipRectangle(pix1, box1, NULL);
2269 else
2270 pix3 = pixClone(pix1);
2271 pixGenPhotoHistos(pix3, NULL, factor, 0, n, &naa1, &w1c, &h1c, debugindex);
2272 pixDestroy(&pix3);
2273 if (!naa1) return 0;
2274 if (box2)
2275 pix4 = pixClipRectangle(pix2, box2, NULL);
2276 else
2277 pix4 = pixClone(pix2);
2278 pixGenPhotoHistos(pix4, NULL, factor, 0, n, &naa2, &w2c, &h2c, debugindex);
2279 pixDestroy(&pix4);
2280 if (!naa2) return 0;
2281
2282 /* Compare histograms */
2283 pixa = (debugflag) ? pixaCreate(0) : NULL;
2284 compareTilesByHisto(naa1, naa2, minratio, w1c, h1c, w2c, h2c, pscore, pixa);
2285 pixaDestroy(&pixa);
2286 return 0;
2287}
2288
2289
2328l_ok
2329pixGenPhotoHistos(PIX *pixs,
2330 BOX *box,
2331 l_int32 factor,
2332 l_float32 thresh,
2333 l_int32 n,
2334 NUMAA **pnaa,
2335 l_int32 *pw,
2336 l_int32 *ph,
2337 l_int32 debugindex)
2338{
2339char buf[64];
2340NUMAA *naa;
2341PIX *pix1, *pix2, *pix3, *pixm;
2342PIXA *pixa;
2343
2344 if (pnaa) *pnaa = NULL;
2345 if (pw) *pw = 0;
2346 if (ph) *ph = 0;
2347 if (!pnaa)
2348 return ERROR_INT("&naa not defined", __func__, 1);
2349 if (!pw || !ph)
2350 return ERROR_INT("&w and &h not both defined", __func__, 1);
2351 if (!pixs || pixGetDepth(pixs) == 1)
2352 return ERROR_INT("pixs not defined or 1 bpp", __func__, 1);
2353 if (factor < 1)
2354 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2355 if (thresh <= 0.0) thresh = 1.3f; /* default */
2356 if (n < 1 || n > 7) {
2357 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2358 n = 4;
2359 }
2360
2361 pixa = NULL;
2362 if (debugindex > 0) {
2363 pixa = pixaCreate(0);
2364 lept_mkdir("lept/comp");
2365 }
2366
2367 /* Initial crop, if necessary */
2368 if (box)
2369 pix1 = pixClipRectangle(pixs, box, NULL);
2370 else
2371 pix1 = pixClone(pixs);
2372
2373 /* Convert to 8 bpp and pad to center the centroid */
2374 pix2 = pixConvertTo8(pix1, FALSE);
2375 pix3 = pixPadToCenterCentroid(pix2, factor);
2376
2377 /* Set to 255 all pixels above 230. Do this so that light gray
2378 * pixels do not enter into the comparison. */
2379 pixm = pixThresholdToBinary(pix3, 230);
2380 pixInvert(pixm, pixm);
2381 pixSetMaskedGeneral(pix3, pixm, 255, 0, 0);
2382 pixDestroy(&pixm);
2383
2384 if (debugindex > 0) {
2385 PIX *pix4, *pix5, *pix6, *pix7, *pix8;
2386 PIXA *pixa2;
2387 pix4 = pixConvertTo32(pix2);
2388 pix5 = pixConvertTo32(pix3);
2389 pix6 = pixScaleToSize(pix4, 400, 0);
2390 pix7 = pixScaleToSize(pix5, 400, 0);
2391 pixa2 = pixaCreate(2);
2392 pixaAddPix(pixa2, pix6, L_INSERT);
2393 pixaAddPix(pixa2, pix7, L_INSERT);
2394 pix8 = pixaDisplayTiledInRows(pixa2, 32, 1000, 1.0, 0, 50, 3);
2395 pixaAddPix(pixa, pix8, L_INSERT);
2396 pixDestroy(&pix4);
2397 pixDestroy(&pix5);
2398 pixaDestroy(&pixa2);
2399 }
2400 pixDestroy(&pix1);
2401 pixDestroy(&pix2);
2402
2403 /* Test if this is a photoimage */
2404 pixDecideIfPhotoImage(pix3, factor, thresh, n, &naa, pixa);
2405 if (naa) {
2406 *pnaa = naa;
2407 *pw = pixGetWidth(pix3);
2408 *ph = pixGetHeight(pix3);
2409 }
2410
2411 if (pixa) {
2412 snprintf(buf, sizeof(buf), "/tmp/lept/comp/tiledhistos.%d.pdf",
2413 debugindex);
2414 lept_stderr("Writing to %s\n", buf);
2415 pixaConvertToPdf(pixa, 300, 1.0, L_FLATE_ENCODE, 0, NULL, buf);
2416 pixaDestroy(&pixa);
2417 }
2418
2419 pixDestroy(&pix3);
2420 return 0;
2421}
2422
2423
2439PIX *
2440pixPadToCenterCentroid(PIX *pixs,
2441 l_int32 factor)
2442
2443{
2444l_float32 cx, cy;
2445l_int32 xs, ys, delx, dely, icx, icy, ws, hs, wd, hd;
2446PIX *pix1, *pixd;
2447
2448 if (!pixs)
2449 return (PIX *)ERROR_PTR("pixs not defined", __func__, NULL);
2450 if (factor < 1)
2451 return (PIX *)ERROR_PTR("invalid sampling factor", __func__, NULL);
2452
2453 pix1 = pixConvertTo8(pixs, FALSE);
2454 pixCentroid8(pix1, factor, &cx, &cy);
2455 icx = (l_int32)(cx + 0.5);
2456 icy = (l_int32)(cy + 0.5);
2457 pixGetDimensions(pix1, &ws, &hs, NULL);
2458 delx = ws - 2 * icx;
2459 dely = hs - 2 * icy;
2460 xs = L_MAX(0, delx);
2461 ys = L_MAX(0, dely);
2462 wd = 2 * L_MAX(icx, ws - icx);
2463 hd = 2 * L_MAX(icy, hs - icy);
2464 pixd = pixCreate(wd, hd, 8);
2465 pixSetAll(pixd); /* to white */
2466 pixCopyResolution(pixd, pixs);
2467 pixRasterop(pixd, xs, ys, ws, hs, PIX_SRC, pix1, 0, 0);
2468 pixDestroy(&pix1);
2469 return pixd;
2470}
2471
2472
2491l_ok
2492pixCentroid8(PIX *pixs,
2493 l_int32 factor,
2494 l_float32 *pcx,
2495 l_float32 *pcy)
2496{
2497l_int32 i, j, w, h, wpl, val;
2498l_float32 sumx, sumy, sumv;
2499l_uint32 *data, *line;
2500PIX *pix1;
2501
2502 if (pcx) *pcx = 0.0;
2503 if (pcy) *pcy = 0.0;
2504 if (!pixs || pixGetDepth(pixs) != 8)
2505 return ERROR_INT("pixs undefined or not 8 bpp", __func__, 1);
2506 if (factor < 1)
2507 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2508 if (!pcx || !pcy)
2509 return ERROR_INT("&cx and &cy not both defined", __func__, 1);
2510
2511 pix1 = pixInvert(NULL, pixs);
2512 pixGetDimensions(pix1, &w, &h, NULL);
2513 data = pixGetData(pix1);
2514 wpl = pixGetWpl(pix1);
2515 sumx = sumy = sumv = 0.0;
2516 for (i = 0; i < h; i++) {
2517 line = data + i * wpl;
2518 for (j = 0; j < w; j++) {
2519 val = GET_DATA_BYTE(line, j);
2520 sumx += val * j;
2521 sumy += val * i;
2522 sumv += val;
2523 }
2524 }
2525 pixDestroy(&pix1);
2526
2527 if (sumv == 0) {
2528 L_INFO("input image is white\n", __func__);
2529 *pcx = (l_float32)(w) / 2;
2530 *pcy = (l_float32)(h) / 2;
2531 } else {
2532 *pcx = sumx / sumv;
2533 *pcy = sumy / sumv;
2534 }
2535
2536 return 0;
2537}
2538
2539
2575l_ok
2576pixDecideIfPhotoImage(PIX *pix,
2577 l_int32 factor,
2578 l_float32 thresh,
2579 l_int32 n,
2580 NUMAA **pnaa,
2581 PIXA *pixadebug)
2582{
2583char buf[64];
2584l_int32 i, w, h, nx, ny, ngrids, istext, isphoto;
2585l_float32 maxval, sum1, sum2, ratio;
2586L_BMF *bmf;
2587NUMA *na1, *na2, *na3, *narv;
2588NUMAA *naa;
2589PIX *pix1;
2590PIXA *pixa1, *pixa2, *pixa3;
2591
2592 if (!pnaa)
2593 return ERROR_INT("&naa not defined", __func__, 1);
2594 *pnaa = NULL;
2595 if (!pix || pixGetDepth(pix) != 8 || pixGetColormap(pix))
2596 return ERROR_INT("pix undefined or invalid", __func__, 1);
2597 if (n < 1 || n > 7) {
2598 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2599 n = 4;
2600 }
2601 if (thresh <= 0.0) thresh = 1.3f; /* default */
2602
2603 /* Look for text lines */
2604 pixDecideIfText(pix, NULL, &istext, pixadebug);
2605 if (istext) {
2606 L_INFO("Image is text\n", __func__);
2607 return 0;
2608 }
2609
2610 /* Determine grid from n */
2611 pixGetDimensions(pix, &w, &h, NULL);
2612 if (w == 0 || h == 0)
2613 return ERROR_INT("invalid pix dimension", __func__, 1);
2614 findHistoGridDimensions(n, w, h, &nx, &ny, 1);
2615
2616 /* Evaluate histograms in each tile */
2617 pixa1 = pixaSplitPix(pix, nx, ny, 0, 0);
2618 ngrids = nx * ny;
2619 bmf = (pixadebug) ? bmfCreate(NULL, 6) : NULL;
2620 naa = numaaCreate(ngrids);
2621 if (pixadebug) {
2622 lept_rmdir("lept/compplot");
2623 lept_mkdir("lept/compplot");
2624 }
2625 for (i = 0; i < ngrids; i++) {
2626 pix1 = pixaGetPix(pixa1, i, L_CLONE);
2627
2628 /* Get histograms, set white count to 0, normalize max to 255 */
2629 na1 = pixGetGrayHistogram(pix1, factor);
2630 numaSetValue(na1, 255, 0);
2631 na2 = numaWindowedMean(na1, 5); /* do some smoothing */
2632 numaGetMax(na2, &maxval, NULL);
2633 na3 = numaTransform(na2, 0, 255.0 / maxval);
2634 if (pixadebug) {
2635 snprintf(buf, sizeof(buf), "/tmp/lept/compplot/plot.%d", i);
2636 gplotSimple1(na3, GPLOT_PNG, buf, "Histos");
2637 }
2638
2639 numaaAddNuma(naa, na3, L_INSERT);
2640 numaDestroy(&na1);
2641 numaDestroy(&na2);
2642 pixDestroy(&pix1);
2643 }
2644 if (pixadebug) {
2645 pix1 = pixaDisplayTiledInColumns(pixa1, nx, 1.0, 30, 2);
2646 pixaAddPix(pixadebug, pix1, L_INSERT);
2647 pixa2 = pixaReadFiles("/tmp/lept/compplot", ".png");
2648 pixa3 = pixaScale(pixa2, 0.4f, 0.4f);
2649 pix1 = pixaDisplayTiledInColumns(pixa3, nx, 1.0, 30, 2);
2650 pixaAddPix(pixadebug, pix1, L_INSERT);
2651 pixaDestroy(&pixa2);
2652 pixaDestroy(&pixa3);
2653 }
2654
2655 /* Compute the standard deviation between these histos to decide
2656 * if the image is photo or something more like line art,
2657 * which does not support good comparison by tiled histograms. */
2658 grayInterHistogramStats(naa, 5, NULL, NULL, NULL, &narv);
2659
2660 /* For photos, the root variance has a larger weight of
2661 * values in the range [50 ... 150] compared to [200 ... 230],
2662 * than text or line art. For the latter, most of the variance
2663 * between tiles is in the lightest parts of the image, well
2664 * above 150. */
2665 numaGetSumOnInterval(narv, 50, 150, &sum1);
2666 numaGetSumOnInterval(narv, 200, 230, &sum2);
2667 if (sum2 == 0.0) { /* shouldn't happen */
2668 ratio = 0.001f; /* anything very small for debug output */
2669 isphoto = 0; /* be conservative */
2670 } else {
2671 ratio = sum1 / sum2;
2672 isphoto = (ratio > thresh) ? 1 : 0;
2673 }
2674 if (pixadebug) {
2675 if (isphoto)
2676 L_INFO("ratio %f > %f; isphoto is true\n",
2677 __func__, ratio, thresh);
2678 else
2679 L_INFO("ratio %f < %f; isphoto is false\n",
2680 __func__, ratio, thresh);
2681 }
2682 if (isphoto)
2683 *pnaa = naa;
2684 else
2685 numaaDestroy(&naa);
2686 bmfDestroy(&bmf);
2687 numaDestroy(&narv);
2688 pixaDestroy(&pixa1);
2689 return 0;
2690}
2691
2692
2718static l_ok
2719findHistoGridDimensions(l_int32 n,
2720 l_int32 w,
2721 l_int32 h,
2722 l_int32 *pnx,
2723 l_int32 *pny,
2724 l_int32 debug)
2725{
2726l_int32 nx, ny, max;
2727l_float32 ratio;
2728
2729 ratio = (l_float32)w / (l_float32)h;
2730 max = n * n;
2731 nx = ny = n;
2732 while (nx > 1 && ny > 1) {
2733 if (ratio > 2.0) { /* reduce ny */
2734 ny--;
2735 nx = max / ny;
2736 if (debug)
2737 lept_stderr("nx = %d, ny = %d, ratio w/h = %4.2f\n",
2738 nx, ny, ratio);
2739 } else if (ratio < 0.5) { /* reduce nx */
2740 nx--;
2741 ny = max / nx;
2742 if (debug)
2743 lept_stderr("nx = %d, ny = %d, ratio w/h = %4.2f\n",
2744 nx, ny, ratio);
2745 } else { /* we're ok */
2746 if (debug)
2747 lept_stderr("nx = %d, ny = %d, ratio w/h = %4.2f\n",
2748 nx, ny, ratio);
2749 break;
2750 }
2751 ratio = (l_float32)(ny * w) / (l_float32)(nx * h);
2752 }
2753 *pnx = nx;
2754 *pny = ny;
2755 return 0;
2756}
2757
2758
2781l_ok
2782compareTilesByHisto(NUMAA *naa1,
2783 NUMAA *naa2,
2784 l_float32 minratio,
2785 l_int32 w1,
2786 l_int32 h1,
2787 l_int32 w2,
2788 l_int32 h2,
2789 l_float32 *pscore,
2790 PIXA *pixadebug)
2791{
2792char buf1[128], buf2[128];
2793l_int32 i, n;
2794l_float32 wratio, hratio, score, minscore, dist;
2795L_BMF *bmf;
2796NUMA *na1, *na2, *nadist, *nascore;
2797
2798 if (!pscore)
2799 return ERROR_INT("&score not defined", __func__, 1);
2800 *pscore = 0.0;
2801 if (!naa1 || !naa2)
2802 return ERROR_INT("naa1 and naa2 not both defined", __func__, 1);
2803
2804 /* Filter for different sizes */
2805 wratio = (w1 < w2) ? (l_float32)w1 / (l_float32)w2 :
2806 (l_float32)w2 / (l_float32)w1;
2807 hratio = (h1 < h2) ? (l_float32)h1 / (l_float32)h2 :
2808 (l_float32)h2 / (l_float32)h1;
2809 if (wratio < minratio || hratio < minratio) {
2810 if (pixadebug)
2811 L_INFO("Sizes differ: wratio = %f, hratio = %f\n",
2812 __func__, wratio, hratio);
2813 return 0;
2814 }
2815 n = numaaGetCount(naa1);
2816 if (n != numaaGetCount(naa2)) { /* due to differing w/h ratio */
2817 L_INFO("naa1 and naa2 sizes are different\n", __func__);
2818 return 0;
2819 }
2820
2821 if (pixadebug) {
2822 lept_rmdir("lept/comptile");
2823 lept_mkdir("lept/comptile");
2824 }
2825
2826
2827 /* Evaluate histograms in each tile. Remove white before
2828 * computing EMD, because there are may be a lot of white
2829 * pixels due to padding, and we don't want to include them.
2830 * This also makes the debug histo plots more informative. */
2831 minscore = 1.0;
2832 nadist = numaCreate(n);
2833 nascore = numaCreate(n);
2834 bmf = (pixadebug) ? bmfCreate(NULL, 6) : NULL;
2835 for (i = 0; i < n; i++) {
2836 na1 = numaaGetNuma(naa1, i, L_CLONE);
2837 na2 = numaaGetNuma(naa2, i, L_CLONE);
2838 numaSetValue(na1, 255, 0.0);
2839 numaSetValue(na2, 255, 0.0);
2840
2841 /* To compare histograms, use the normalized earthmover distance.
2842 * Further normalize to get the EM distance as a fraction of the
2843 * maximum distance in the histogram (255). Finally, scale this
2844 * up by 10.0, and subtract from 1.0 to get a similarity score. */
2845 numaEarthMoverDistance(na1, na2, &dist);
2846 score = L_MAX(0.0, 1.0 - 10.0 * (dist / 255.));
2847 numaAddNumber(nadist, dist);
2848 numaAddNumber(nascore, score);
2849 minscore = L_MIN(minscore, score);
2850 if (pixadebug) {
2851 snprintf(buf1, sizeof(buf1), "/tmp/lept/comptile/plot.%d", i);
2852 gplotSimple2(na1, na2, GPLOT_PNG, buf1, "Histos");
2853 }
2854 numaDestroy(&na1);
2855 numaDestroy(&na2);
2856 }
2857 *pscore = minscore;
2858
2859 if (pixadebug) {
2860 for (i = 0; i < n; i++) {
2861 PIX *pix1, *pix2;
2862 snprintf(buf1, sizeof(buf1), "/tmp/lept/comptile/plot.%d.png", i);
2863 pix1 = pixRead(buf1);
2864 numaGetFValue(nadist, i, &dist);
2865 numaGetFValue(nascore, i, &score);
2866 snprintf(buf2, sizeof(buf2),
2867 "Image %d\ndist = %5.3f, score = %5.3f", i, dist, score);
2868 pix2 = pixAddTextlines(pix1, bmf, buf2, 0x0000ff00, L_ADD_BELOW);
2869 pixaAddPix(pixadebug, pix2, L_INSERT);
2870 pixDestroy(&pix1);
2871 }
2872 lept_stderr("Writing to /tmp/lept/comptile/comparegray.pdf\n");
2873 pixaConvertToPdf(pixadebug, 300, 1.0, L_FLATE_ENCODE, 0, NULL,
2874 "/tmp/lept/comptile/comparegray.pdf");
2875 numaWriteDebug("/tmp/lept/comptile/scores.na", nascore);
2876 numaWriteDebug("/tmp/lept/comptile/dists.na", nadist);
2877 }
2878
2879 bmfDestroy(&bmf);
2880 numaDestroy(&nadist);
2881 numaDestroy(&nascore);
2882 return 0;
2883}
2884
2885
2957l_ok
2958pixCompareGrayByHisto(PIX *pix1,
2959 PIX *pix2,
2960 BOX *box1,
2961 BOX *box2,
2962 l_float32 minratio,
2963 l_int32 maxgray,
2964 l_int32 factor,
2965 l_int32 n,
2966 l_float32 *pscore,
2967 l_int32 debugflag)
2968{
2969l_int32 w1, h1, w2, h2;
2970l_float32 wratio, hratio;
2971BOX *box3, *box4;
2972PIX *pix3, *pix4, *pix5, *pix6, *pix7, *pix8;
2973PIXA *pixa;
2974
2975 if (!pscore)
2976 return ERROR_INT("&score not defined", __func__, 1);
2977 *pscore = 0.0;
2978 if (!pix1 || !pix2)
2979 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
2980 if (minratio < 0.5 || minratio > 1.0)
2981 return ERROR_INT("minratio not in [0.5 ... 1.0]", __func__, 1);
2982 if (maxgray < 200)
2983 return ERROR_INT("invalid maxgray; should be >= 200", __func__, 1);
2984 maxgray = L_MIN(255, maxgray);
2985 if (factor < 1)
2986 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
2987 if (n < 1 || n > 7) {
2988 L_WARNING("n = %d is invalid; setting to 4\n", __func__, n);
2989 n = 4;
2990 }
2991
2992 if (debugflag)
2993 lept_mkdir("lept/comp");
2994
2995 /* Initial filter by size */
2996 if (box1)
2997 boxGetGeometry(box1, NULL, NULL, &w1, &h1);
2998 else
2999 pixGetDimensions(pix1, &w1, &h1, NULL);
3000 if (box2)
3001 boxGetGeometry(box2, NULL, NULL, &w2, &h2);
3002 else
3003 pixGetDimensions(pix1, &w2, &h2, NULL);
3004 wratio = (w1 < w2) ? (l_float32)w1 / (l_float32)w2 :
3005 (l_float32)w2 / (l_float32)w1;
3006 hratio = (h1 < h2) ? (l_float32)h1 / (l_float32)h2 :
3007 (l_float32)h2 / (l_float32)h1;
3008 if (wratio < minratio || hratio < minratio)
3009 return 0;
3010
3011 /* Initial crop, if necessary */
3012 if (box1)
3013 pix3 = pixClipRectangle(pix1, box1, NULL);
3014 else
3015 pix3 = pixClone(pix1);
3016 if (box2)
3017 pix4 = pixClipRectangle(pix2, box2, NULL);
3018 else
3019 pix4 = pixClone(pix2);
3020
3021 /* Convert to 8 bpp, align centroids and do maximal crop */
3022 pix5 = pixConvertTo8(pix3, FALSE);
3023 pix6 = pixConvertTo8(pix4, FALSE);
3024 pixCropAlignedToCentroid(pix5, pix6, factor, &box3, &box4);
3025 pix7 = pixClipRectangle(pix5, box3, NULL);
3026 pix8 = pixClipRectangle(pix6, box4, NULL);
3027 pixa = (debugflag) ? pixaCreate(0) : NULL;
3028 if (debugflag) {
3029 PIX *pix9, *pix10, *pix11, *pix12, *pix13;
3030 PIXA *pixa2;
3031 pix9 = pixConvertTo32(pix5);
3032 pix10 = pixConvertTo32(pix6);
3033 pixRenderBoxArb(pix9, box3, 2, 255, 0, 0);
3034 pixRenderBoxArb(pix10, box4, 2, 255, 0, 0);
3035 pix11 = pixScaleToSize(pix9, 400, 0);
3036 pix12 = pixScaleToSize(pix10, 400, 0);
3037 pixa2 = pixaCreate(2);
3038 pixaAddPix(pixa2, pix11, L_INSERT);
3039 pixaAddPix(pixa2, pix12, L_INSERT);
3040 pix13 = pixaDisplayTiledInRows(pixa2, 32, 1000, 1.0, 0, 50, 0);
3041 pixaAddPix(pixa, pix13, L_INSERT);
3042 pixDestroy(&pix9);
3043 pixDestroy(&pix10);
3044 pixaDestroy(&pixa2);
3045 }
3046 pixDestroy(&pix3);
3047 pixDestroy(&pix4);
3048 pixDestroy(&pix5);
3049 pixDestroy(&pix6);
3050 boxDestroy(&box3);
3051 boxDestroy(&box4);
3052
3053 /* Tile and compare histograms */
3054 pixCompareTilesByHisto(pix7, pix8, maxgray, factor, n, pscore, pixa);
3055 pixaDestroy(&pixa);
3056 pixDestroy(&pix7);
3057 pixDestroy(&pix8);
3058 return 0;
3059}
3060
3061
3082static l_ok
3083pixCompareTilesByHisto(PIX *pix1,
3084 PIX *pix2,
3085 l_int32 maxgray,
3086 l_int32 factor,
3087 l_int32 n,
3088 l_float32 *pscore,
3089 PIXA *pixadebug)
3090{
3091char buf[64];
3092l_int32 w, h, i, j, nx, ny, ngr;
3093l_float32 score, minscore, maxval1, maxval2, dist;
3094L_BMF *bmf;
3095NUMA *na1, *na2, *na3, *na4, *na5, *na6, *na7;
3096PIX *pix3, *pix4;
3097PIXA *pixa1, *pixa2;
3098
3099 if (!pscore)
3100 return ERROR_INT("&score not defined", __func__, 1);
3101 *pscore = 0.0;
3102 if (!pix1 || !pix2)
3103 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
3104
3105 /* Determine grid from n */
3106 pixGetDimensions(pix1, &w, &h, NULL);
3107 findHistoGridDimensions(n, w, h, &nx, &ny, 1);
3108 ngr = nx * ny;
3109
3110 /* Evaluate histograms in each tile */
3111 pixa1 = pixaSplitPix(pix1, nx, ny, 0, 0);
3112 pixa2 = pixaSplitPix(pix2, nx, ny, 0, 0);
3113 na7 = (pixadebug) ? numaCreate(ngr) : NULL;
3114 bmf = (pixadebug) ? bmfCreate(NULL, 6) : NULL;
3115 minscore = 1.0;
3116 for (i = 0; i < ngr; i++) {
3117 pix3 = pixaGetPix(pixa1, i, L_CLONE);
3118 pix4 = pixaGetPix(pixa2, i, L_CLONE);
3119
3120 /* Get histograms, set white count to 0, normalize max to 255 */
3121 na1 = pixGetGrayHistogram(pix3, factor);
3122 na2 = pixGetGrayHistogram(pix4, factor);
3123 if (maxgray < 255) {
3124 for (j = maxgray + 1; j <= 255; j++) {
3125 numaSetValue(na1, j, 0);
3126 numaSetValue(na2, j, 0);
3127 }
3128 }
3129 na3 = numaWindowedMean(na1, 5);
3130 na4 = numaWindowedMean(na2, 5);
3131 numaGetMax(na3, &maxval1, NULL);
3132 numaGetMax(na4, &maxval2, NULL);
3133 na5 = numaTransform(na3, 0, 255.0 / maxval1);
3134 na6 = numaTransform(na4, 0, 255.0 / maxval2);
3135 if (pixadebug) {
3136 gplotSimple2(na5, na6, GPLOT_PNG, "/tmp/lept/comp/plot1", "Histos");
3137 }
3138
3139 /* To compare histograms, use the normalized earthmover distance.
3140 * Further normalize to get the EM distance as a fraction of the
3141 * maximum distance in the histogram (255). Finally, scale this
3142 * up by 10.0, and subtract from 1.0 to get a similarity score. */
3143 numaEarthMoverDistance(na5, na6, &dist);
3144 score = L_MAX(0.0, 1.0 - 8.0 * (dist / 255.));
3145 if (pixadebug) numaAddNumber(na7, score);
3146 minscore = L_MIN(minscore, score);
3147 if (pixadebug) {
3148 PIX *pix5, *pix6, *pix7, *pix8, *pix9, *pix10;
3149 PIXA *pixa3;
3150 l_int32 w, h, wscale;
3151 pixa3 = pixaCreate(3);
3152 pixGetDimensions(pix3, &w, &h, NULL);
3153 wscale = (w > h) ? 700 : 400;
3154 pix5 = pixScaleToSize(pix3, wscale, 0);
3155 pix6 = pixScaleToSize(pix4, wscale, 0);
3156 pixaAddPix(pixa3, pix5, L_INSERT);
3157 pixaAddPix(pixa3, pix6, L_INSERT);
3158 pix7 = pixRead("/tmp/lept/comp/plot1.png");
3159 pix8 = pixScaleToSize(pix7, 700, 0);
3160 snprintf(buf, sizeof(buf), "%5.3f", score);
3161 pix9 = pixAddTextlines(pix8, bmf, buf, 0x0000ff00, L_ADD_RIGHT);
3162 pixaAddPix(pixa3, pix9, L_INSERT);
3163 pix10 = pixaDisplayTiledInRows(pixa3, 32, 1000, 1.0, 0, 50, 0);
3164 pixaAddPix(pixadebug, pix10, L_INSERT);
3165 pixDestroy(&pix7);
3166 pixDestroy(&pix8);
3167 pixaDestroy(&pixa3);
3168 }
3169 numaDestroy(&na1);
3170 numaDestroy(&na2);
3171 numaDestroy(&na3);
3172 numaDestroy(&na4);
3173 numaDestroy(&na5);
3174 numaDestroy(&na6);
3175 pixDestroy(&pix3);
3176 pixDestroy(&pix4);
3177 }
3178 *pscore = minscore;
3179
3180 if (pixadebug) {
3181 pixaConvertToPdf(pixadebug, 300, 1.0, L_FLATE_ENCODE, 0, NULL,
3182 "/tmp/lept/comp/comparegray.pdf");
3183 numaWriteDebug("/tmp/lept/comp/tilescores.na", na7);
3184 }
3185
3186 bmfDestroy(&bmf);
3187 numaDestroy(&na7);
3188 pixaDestroy(&pixa1);
3189 pixaDestroy(&pixa2);
3190 return 0;
3191}
3192
3193
3210l_ok
3211pixCropAlignedToCentroid(PIX *pix1,
3212 PIX *pix2,
3213 l_int32 factor,
3214 BOX **pbox1,
3215 BOX **pbox2)
3216{
3217l_float32 cx1, cy1, cx2, cy2;
3218l_int32 w1, h1, w2, h2, icx1, icy1, icx2, icy2;
3219l_int32 xm, xm1, xm2, xp, xp1, xp2, ym, ym1, ym2, yp, yp1, yp2;
3220PIX *pix3, *pix4;
3221
3222 if (pbox1) *pbox1 = NULL;
3223 if (pbox2) *pbox2 = NULL;
3224 if (!pix1 || !pix2)
3225 return ERROR_INT("pix1 and pix2 not both defined", __func__, 1);
3226 if (factor < 1)
3227 return ERROR_INT("subsampling factor must be >= 1", __func__, 1);
3228 if (!pbox1 || !pbox2)
3229 return ERROR_INT("&box1 and &box2 not both defined", __func__, 1);
3230
3231 pix3 = pixConvertTo8(pix1, FALSE);
3232 pix4 = pixConvertTo8(pix2, FALSE);
3233 pixCentroid8(pix3, factor, &cx1, &cy1);
3234 pixCentroid8(pix4, factor, &cx2, &cy2);
3235 pixGetDimensions(pix3, &w1, &h1, NULL);
3236 pixGetDimensions(pix4, &w2, &h2, NULL);
3237 pixDestroy(&pix3);
3238 pixDestroy(&pix4);
3239
3240 icx1 = (l_int32)(cx1 + 0.5);
3241 icy1 = (l_int32)(cy1 + 0.5);
3242 icx2 = (l_int32)(cx2 + 0.5);
3243 icy2 = (l_int32)(cy2 + 0.5);
3244 xm = L_MIN(icx1, icx2);
3245 xm1 = icx1 - xm;
3246 xm2 = icx2 - xm;
3247 xp = L_MIN(w1 - icx1, w2 - icx2); /* one pixel beyond to the right */
3248 xp1 = icx1 + xp;
3249 xp2 = icx2 + xp;
3250 ym = L_MIN(icy1, icy2);
3251 ym1 = icy1 - ym;
3252 ym2 = icy2 - ym;
3253 yp = L_MIN(h1 - icy1, h2 - icy2); /* one pixel below the bottom */
3254 yp1 = icy1 + yp;
3255 yp2 = icy2 + yp;
3256 *pbox1 = boxCreate(xm1, ym1, xp1 - xm1, yp1 - ym1);
3257 *pbox2 = boxCreate(xm2, ym2, xp2 - xm2, yp2 - ym2);
3258 return 0;
3259}
3260
3261
3283l_uint8 *
3284l_compressGrayHistograms(NUMAA *naa,
3285 l_int32 w,
3286 l_int32 h,
3287 size_t *psize)
3288{
3289l_uint8 *bytea;
3290l_int32 i, j, n, nn, ival;
3291l_float32 maxval;
3292NUMA *na1, *na2;
3293
3294 if (!psize)
3295 return (l_uint8 *)ERROR_PTR("&size not defined", __func__, NULL);
3296 *psize = 0;
3297 if (!naa)
3298 return (l_uint8 *)ERROR_PTR("naa not defined", __func__, NULL);
3299 n = numaaGetCount(naa);
3300 for (i = 0; i < n; i++) {
3301 nn = numaaGetNumaCount(naa, i);
3302 if (nn != 256) {
3303 L_ERROR("%d numbers in numa[%d]\n", __func__, nn, i);
3304 return NULL;
3305 }
3306 }
3307
3308 if ((bytea = (l_uint8 *)LEPT_CALLOC(8 + 256 * n, sizeof(l_uint8))) == NULL)
3309 return (l_uint8 *)ERROR_PTR("bytea not made", __func__, NULL);
3310 *psize = 8 + 256 * n;
3311 l_setDataFourBytes(bytea, 0, w);
3312 l_setDataFourBytes(bytea, 1, h);
3313 for (i = 0; i < n; i++) {
3314 na1 = numaaGetNuma(naa, i, L_COPY);
3315 numaGetMax(na1, &maxval, NULL);
3316 na2 = numaTransform(na1, 0, 255.0 / maxval);
3317 for (j = 0; j < 256; j++) {
3318 numaGetIValue(na2, j, &ival);
3319 bytea[8 + 256 * i + j] = ival;
3320 }
3321 numaDestroy(&na1);
3322 numaDestroy(&na2);
3323 }
3324
3325 return bytea;
3326}
3327
3328
3349NUMAA *
3350l_uncompressGrayHistograms(l_uint8 *bytea,
3351 size_t size,
3352 l_int32 *pw,
3353 l_int32 *ph)
3354{
3355l_int32 i, j, n;
3356NUMA *na;
3357NUMAA *naa;
3358
3359 if (pw) *pw = 0;
3360 if (ph) *ph = 0;
3361 if (!pw || !ph)
3362 return (NUMAA *)ERROR_PTR("&w and &h not both defined", __func__, NULL);
3363 if (!bytea)
3364 return (NUMAA *)ERROR_PTR("bytea not defined", __func__, NULL);
3365 n = (size - 8) / 256;
3366 if ((size - 8) % 256 != 0)
3367 return (NUMAA *)ERROR_PTR("bytea size is invalid", __func__, NULL);
3368
3369 *pw = l_getDataFourBytes(bytea, 0);
3370 *ph = l_getDataFourBytes(bytea, 1);
3371 naa = numaaCreate(n);
3372 for (i = 0; i < n; i++) {
3373 na = numaCreate(256);
3374 for (j = 0; j < 256; j++)
3375 numaAddNumber(na, bytea[8 + 256 * i + j]);
3376 numaaAddNuma(naa, na, L_INSERT);
3377 }
3378
3379 return naa;
3380}
3381
3382
3383/*------------------------------------------------------------------*
3384 * Translated images at the same resolution *
3385 *------------------------------------------------------------------*/
3416l_ok
3417pixCompareWithTranslation(PIX *pix1,
3418 PIX *pix2,
3419 l_int32 thresh,
3420 l_int32 *pdelx,
3421 l_int32 *pdely,
3422 l_float32 *pscore,
3423 l_int32 debugflag)
3424{
3425l_uint8 *subtab;
3426l_int32 i, level, area1, area2, delx, dely;
3427l_int32 etransx, etransy, maxshift, dbint;
3428l_int32 *stab, *ctab;
3429l_float32 cx1, cx2, cy1, cy2, score;
3430PIX *pixb1, *pixb2, *pixt1, *pixt2, *pixt3, *pixt4;
3431PIXA *pixa1, *pixa2, *pixadb = NULL;
3432
3433 if (pdelx) *pdelx = 0;
3434 if (pdely) *pdely = 0;
3435 if (pscore) *pscore = 0.0;
3436 if (!pdelx || !pdely)
3437 return ERROR_INT("&delx and &dely not defined", __func__, 1);
3438 if (!pscore)
3439 return ERROR_INT("&score not defined", __func__, 1);
3440 if (!pix1)
3441 return ERROR_INT("pix1 not defined", __func__, 1);
3442 if (!pix2)
3443 return ERROR_INT("pix2 not defined", __func__, 1);
3444
3445 /* Make tables */
3446 subtab = makeSubsampleTab2x();
3447 stab = makePixelSumTab8();
3448 ctab = makePixelCentroidTab8();
3449
3450 /* Binarize each image */
3451 pixb1 = pixConvertTo1(pix1, thresh);
3452 pixb2 = pixConvertTo1(pix2, thresh);
3453
3454 /* Make a cascade of 2x reduced images for each, thresholding
3455 * with level 2 (neutral), down to 8x reduction */
3456 pixa1 = pixaCreate(4);
3457 pixa2 = pixaCreate(4);
3458 if (debugflag)
3459 pixadb = pixaCreate(4);
3460 pixaAddPix(pixa1, pixb1, L_INSERT);
3461 pixaAddPix(pixa2, pixb2, L_INSERT);
3462 for (i = 0; i < 3; i++) {
3463 pixt1 = pixReduceRankBinary2(pixb1, 2, subtab);
3464 pixt2 = pixReduceRankBinary2(pixb2, 2, subtab);
3465 pixaAddPix(pixa1, pixt1, L_INSERT);
3466 pixaAddPix(pixa2, pixt2, L_INSERT);
3467 pixb1 = pixt1;
3468 pixb2 = pixt2;
3469 }
3470
3471 /* At the lowest level, use the centroids with a maxshift of 6
3472 * to search for the best alignment. Then at higher levels,
3473 * use the result from the level below as the initial approximation
3474 * for the alignment, and search with a maxshift of 2. */
3475 for (level = 3; level >= 0; level--) {
3476 pixt1 = pixaGetPix(pixa1, level, L_CLONE);
3477 pixt2 = pixaGetPix(pixa2, level, L_CLONE);
3478 pixCountPixels(pixt1, &area1, stab);
3479 pixCountPixels(pixt2, &area2, stab);
3480 if (level == 3) {
3481 pixCentroid(pixt1, ctab, stab, &cx1, &cy1);
3482 pixCentroid(pixt2, ctab, stab, &cx2, &cy2);
3483 etransx = lept_roundftoi(cx1 - cx2);
3484 etransy = lept_roundftoi(cy1 - cy2);
3485 maxshift = 6;
3486 } else {
3487 etransx = 2 * delx;
3488 etransy = 2 * dely;
3489 maxshift = 2;
3490 }
3491 dbint = (debugflag) ? level + 1 : 0;
3492 pixBestCorrelation(pixt1, pixt2, area1, area2, etransx, etransy,
3493 maxshift, stab, &delx, &dely, &score, dbint);
3494 if (debugflag) {
3495 lept_stderr("Level %d: delx = %d, dely = %d, score = %7.4f\n",
3496 level, delx, dely, score);
3497 pixRasteropIP(pixt2, delx, dely, L_BRING_IN_WHITE);
3498 pixt3 = pixDisplayDiffBinary(pixt1, pixt2);
3499 pixt4 = pixExpandReplicate(pixt3, 8 / (1 << (3 - level)));
3500 pixaAddPix(pixadb, pixt4, L_INSERT);
3501 pixDestroy(&pixt3);
3502 }
3503 pixDestroy(&pixt1);
3504 pixDestroy(&pixt2);
3505 }
3506
3507 if (debugflag) {
3508 pixaConvertToPdf(pixadb, 300, 1.0, L_FLATE_ENCODE, 0, NULL,
3509 "/tmp/lept/comp/compare.pdf");
3510 convertFilesToPdf("/tmp/lept/comp", "correl_", 30, 1.0, L_FLATE_ENCODE,
3511 0, "Correlation scores at levels 1 through 5",
3512 "/tmp/lept/comp/correl.pdf");
3513 pixaDestroy(&pixadb);
3514 }
3515
3516 *pdelx = delx;
3517 *pdely = dely;
3518 *pscore = score;
3519 pixaDestroy(&pixa1);
3520 pixaDestroy(&pixa2);
3521 LEPT_FREE(subtab);
3522 LEPT_FREE(stab);
3523 LEPT_FREE(ctab);
3524 return 0;
3525}
3526
3527
3568l_ok
3569pixBestCorrelation(PIX *pix1,
3570 PIX *pix2,
3571 l_int32 area1,
3572 l_int32 area2,
3573 l_int32 etransx,
3574 l_int32 etransy,
3575 l_int32 maxshift,
3576 l_int32 *tab8,
3577 l_int32 *pdelx,
3578 l_int32 *pdely,
3579 l_float32 *pscore,
3580 l_int32 debugflag)
3581{
3582l_int32 shiftx, shifty, delx, dely;
3583l_int32 *tab;
3584l_float32 maxscore, score;
3585FPIX *fpix = NULL;
3586PIX *pix3, *pix4;
3587
3588 if (pdelx) *pdelx = 0;
3589 if (pdely) *pdely = 0;
3590 if (pscore) *pscore = 0.0;
3591 if (!pix1 || pixGetDepth(pix1) != 1)
3592 return ERROR_INT("pix1 not defined or not 1 bpp", __func__, 1);
3593 if (!pix2 || pixGetDepth(pix2) != 1)
3594 return ERROR_INT("pix2 not defined or not 1 bpp", __func__, 1);
3595 if (!area1 || !area2)
3596 return ERROR_INT("areas must be > 0", __func__, 1);
3597
3598 if (debugflag > 0)
3599 fpix = fpixCreate(2 * maxshift + 1, 2 * maxshift + 1);
3600
3601 if (!tab8)
3602 tab = makePixelSumTab8();
3603 else
3604 tab = tab8;
3605
3606 /* Search over a set of {shiftx, shifty} for the max */
3607 maxscore = 0;
3608 delx = etransx;
3609 dely = etransy;
3610 for (shifty = -maxshift; shifty <= maxshift; shifty++) {
3611 for (shiftx = -maxshift; shiftx <= maxshift; shiftx++) {
3612 pixCorrelationScoreShifted(pix1, pix2, area1, area2,
3613 etransx + shiftx,
3614 etransy + shifty, tab, &score);
3615 if (debugflag > 0) {
3616 fpixSetPixel(fpix, maxshift + shiftx, maxshift + shifty,
3617 1000.0 * score);
3618/* lept_stderr("(sx, sy) = (%d, %d): score = %6.4f\n",
3619 shiftx, shifty, score); */
3620 }
3621 if (score > maxscore) {
3622 maxscore = score;
3623 delx = etransx + shiftx;
3624 dely = etransy + shifty;
3625 }
3626 }
3627 }
3628
3629 if (debugflag > 0) {
3630 char buf[128];
3631 lept_mkdir("lept/comp");
3632 pix3 = fpixDisplayMaxDynamicRange(fpix);
3633 pix4 = pixExpandReplicate(pix3, 20);
3634 snprintf(buf, sizeof(buf), "/tmp/lept/comp/correl_%d.png",
3635 debugflag);
3636 pixWrite(buf, pix4, IFF_PNG);
3637 pixDestroy(&pix3);
3638 pixDestroy(&pix4);
3639 fpixDestroy(&fpix);
3640 }
3641
3642 if (pdelx) *pdelx = delx;
3643 if (pdely) *pdely = dely;
3644 if (pscore) *pscore = maxscore;
3645 if (!tab8) LEPT_FREE(tab);
3646 return 0;
3647}
SET_DATA_FOUR_BYTES
#define SET_DATA_FOUR_BYTES(pdata, n, val)
Definition arrayaccess.h:235
GET_DATA_BYTE
#define GET_DATA_BYTE(pdata, n)
Definition arrayaccess.h:188
GET_DATA_FOUR_BYTES
#define GET_DATA_FOUR_BYTES(pdata, n)
Definition arrayaccess.h:231
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:883
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:2782
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:1994
pixDisplayDiff
PIX * pixDisplayDiff(PIX *pix1, PIX *pix2, l_int32 showall, l_int32 mindiff, l_uint32 diffcolor)
pixDisplayDiff()
Definition compare.c:657
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:1477
pixCropAlignedToCentroid
l_ok pixCropAlignedToCentroid(PIX *pix1, PIX *pix2, l_int32 factor, BOX **pbox1, BOX **pbox2)
pixCropAlignedToCentroid()
Definition compare.c:3211
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:2958
pixDecideIfPhotoImage
l_ok pixDecideIfPhotoImage(PIX *pix, l_int32 factor, l_float32 thresh, l_int32 n, NUMAA **pnaa, PIXA *pixadebug)
pixDecideIfPhotoImage()
Definition compare.c:2576
pixDisplayDiffBinary
PIX * pixDisplayDiffBinary(PIX *pix1, PIX *pix2)
pixDisplayDiffBinary()
Definition compare.c:750
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:1071
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:2329
pixGetPSNR
l_ok pixGetPSNR(PIX *pix1, PIX *pix2, l_int32 factor, l_float32 *ppsnr)
pixGetPSNR()
Definition compare.c:1871
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:3569
cmapEqual
l_ok cmapEqual(PIXCMAP *cmap1, PIXCMAP *cmap2, l_int32 ncomps, l_int32 *psame)
cmapEqual()
Definition compare.c:477
pixGetDifferenceHistogram
NUMA * pixGetDifferenceHistogram(PIX *pix1, PIX *pix2, l_int32 factor)
pixGetDifferenceHistogram()
Definition compare.c:1577
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:3083
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:3417
l_compressGrayHistograms
l_uint8 * l_compressGrayHistograms(NUMAA *naa, l_int32 w, l_int32 h, size_t *psize)
l_compressGrayHistograms()
Definition compare.c:3284
pixEqual
l_ok pixEqual(PIX *pix1, PIX *pix2, l_int32 *psame)
pixEqual()
Definition compare.c:157
pixPadToCenterCentroid
PIX * pixPadToCenterCentroid(PIX *pixs, l_int32 factor)
pixPadToCenterCentroid()
Definition compare.c:2440
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:2214
pixEqualWithCmap
l_ok pixEqualWithCmap(PIX *pix1, PIX *pix2, l_int32 *psame)
pixEqualWithCmap()
Definition compare.c:383
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:964
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:1400
pixCompareRankDifference
NUMA * pixCompareRankDifference(PIX *pix1, PIX *pix2, l_int32 factor)
pixCompareRankDifference()
Definition compare.c:1317
pixCorrelationBinary
l_ok pixCorrelationBinary(PIX *pix1, PIX *pix2, l_float32 *pval)
pixCorrelationBinary()
Definition compare.c:597
pixCompareTiled
l_ok pixCompareTiled(PIX *pix1, PIX *pix2, l_int32 sx, l_int32 sy, l_int32 type, PIX **ppixdiff)
pixCompareTiled()
Definition compare.c:1225
pixCompareBinary
l_ok pixCompareBinary(PIX *pix1, PIX *pix2, l_int32 comptype, l_float32 *pfract, PIX **ppixdiff)
pixCompareBinary()
Definition compare.c:805
pixUsesCmapColor
l_ok pixUsesCmapColor(PIX *pixs, l_int32 *pcolor)
pixUsesCmapColor()
Definition compare.c:533
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:2719
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:1706
l_uncompressGrayHistograms
NUMAA * l_uncompressGrayHistograms(l_uint8 *bytea, size_t size, l_int32 *pw, l_int32 *ph)
l_uncompressGrayHistograms()
Definition compare.c:3350
pixCentroid8
l_ok pixCentroid8(PIX *pixs, l_int32 factor, l_float32 *pcx, l_float32 *pcy)
pixCentroid8()
Definition compare.c:2492
pixEqualWithAlpha
l_ok pixEqualWithAlpha(PIX *pix1, PIX *pix2, l_int32 use_alpha, l_int32 *psame)
pixEqualWithAlpha()
Definition compare.c:183
L_FLATE_ENCODE
@ L_FLATE_ENCODE
Definition imageio.h:161
COLOR_BLUE
@ COLOR_BLUE
Definition pix.h:330
COLOR_RED
@ COLOR_RED
Definition pix.h:328
COLOR_GREEN
@ COLOR_GREEN
Definition pix.h:329
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
L_ADD_BELOW
@ L_ADD_BELOW
Definition pix.h:1003
L_ADD_RIGHT
@ L_ADD_RIGHT
Definition pix.h:1005
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_SRC
#define PIX_SRC
Definition pix.h:444
L_ROOT_MEAN_SQUARE
@ L_ROOT_MEAN_SQUARE
Definition pix.h:765
L_MEAN_ABSVAL
@ L_MEAN_ABSVAL
Definition pix.h:761
L_BRING_IN_WHITE
@ L_BRING_IN_WHITE
Definition pix.h:662
GPlot
Definition gplot.h:77
L_Bmf
Definition bmf.h:47