Leptonica  1.82.0
Image processing and image analysis suite
skew.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 
98 #ifdef HAVE_CONFIG_H
99 #include <config_auto.h>
100 #endif /* HAVE_CONFIG_H */
101 
102 #include <math.h>
103 #include "allheaders.h"
104 
105  /* Default sweep angle parameters for pixFindSkew() */
106 static const l_float32 DefaultSweepRange = 7.0; /* degrees */
107 static const l_float32 DefaultSweepDelta = 1.0; /* degrees */
108 
109  /* Default final angle difference parameter for binary
110  * search in pixFindSkew(). The expected accuracy is
111  * not better than the inverse image width in pixels,
112  * say, 1/2000 radians, or about 0.03 degrees. */
113 static const l_float32 DefaultMinbsDelta = 0.01; /* degrees */
114 
115  /* Default scale factors for pixFindSkew() */
116 static const l_int32 DefaultSweepReduction = 4; /* sweep part; 4 is good */
117 static const l_int32 DefaultBsReduction = 2; /* binary search part */
118 
119  /* Minimum angle for deskewing in pixDeskew() */
120 static const l_float32 MinDeskewAngle = 0.1; /* degree */
121 
122  /* Minimum allowed confidence (ratio) for deskewing in pixDeskew() */
123 static const l_float32 MinAllowedConfidence = 3.0;
124 
125  /* Minimum allowed maxscore to give nonzero confidence */
126 static const l_int32 MinValidMaxscore = 10000;
127 
128  /* Constant setting threshold for minimum allowed minscore
129  * to give nonzero confidence; multiply this constant by
130  * (height * width^2) */
131 static const l_float32 MinscoreThreshFactor = 0.000002;
132 
133  /* Default binarization threshold value */
134 static const l_int32 DefaultBinaryThreshold = 130;
135 
136 #ifndef NO_CONSOLE_IO
137 #define DEBUG_PRINT_SCORES 0
138 #define DEBUG_PRINT_SWEEP 0
139 #define DEBUG_PRINT_BINARY 0
140 #define DEBUG_PRINT_ORTH 0
141 #define DEBUG_THRESHOLD 0
142 #define DEBUG_PLOT_SCORES 0 /* requires the gnuplot executable */
143 #endif /* ~NO_CONSOLE_IO */
144 
145 
146 
147 /*-----------------------------------------------------------------------*
148  * Top-level deskew interfaces *
149  *-----------------------------------------------------------------------*/
166 PIX *
168  l_int32 redsearch)
169 {
170 PIX *pix1, *pix2, *pix3, *pix4;
171 
172  PROCNAME("pixDeskewBoth");
173 
174  if (!pixs)
175  return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
176  if (redsearch == 0)
177  redsearch = DefaultBsReduction;
178  else if (redsearch != 1 && redsearch != 2 && redsearch != 4)
179  return (PIX *)ERROR_PTR("redsearch not in {1,2,4}", procName, NULL);
180 
181  pix1 = pixDeskew(pixs, redsearch);
182  pix2 = pixRotate90(pix1, 1);
183  pix3 = pixDeskew(pix2, redsearch);
184  pix4 = pixRotate90(pix3, -1);
185  pixDestroy(&pix1);
186  pixDestroy(&pix2);
187  pixDestroy(&pix3);
188  return pix4;
189 }
190 
191 
209 PIX *
211  l_int32 redsearch)
212 {
213  PROCNAME("pixDeskew");
214 
215  if (!pixs)
216  return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
217  if (redsearch == 0)
218  redsearch = DefaultBsReduction;
219  else if (redsearch != 1 && redsearch != 2 && redsearch != 4)
220  return (PIX *)ERROR_PTR("redsearch not in {1,2,4}", procName, NULL);
221 
222  return pixDeskewGeneral(pixs, 0, 0.0, 0.0, redsearch, 0, NULL, NULL);
223 }
224 
225 
245 PIX *
247  l_int32 redsearch,
248  l_float32 *pangle,
249  l_float32 *pconf)
250 {
251  PROCNAME("pixFindSkewAndDeskew");
252 
253  if (!pixs)
254  return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
255  if (redsearch == 0)
256  redsearch = DefaultBsReduction;
257  else if (redsearch != 1 && redsearch != 2 && redsearch != 4)
258  return (PIX *)ERROR_PTR("redsearch not in {1,2,4}", procName, NULL);
259 
260  return pixDeskewGeneral(pixs, 0, 0.0, 0.0, redsearch, 0, pangle, pconf);
261 }
262 
263 
289 PIX *
291  l_int32 redsweep,
292  l_float32 sweeprange,
293  l_float32 sweepdelta,
294  l_int32 redsearch,
295  l_int32 thresh,
296  l_float32 *pangle,
297  l_float32 *pconf)
298 {
299 l_int32 ret, depth;
300 l_float32 angle, conf, deg2rad;
301 PIX *pixb, *pixd;
302 
303  PROCNAME("pixDeskewGeneral");
304 
305  if (pangle) *pangle = 0.0;
306  if (pconf) *pconf = 0.0;
307  if (!pixs)
308  return (PIX *)ERROR_PTR("pixs not defined", procName, NULL);
309  if (redsweep == 0)
310  redsweep = DefaultSweepReduction;
311  else if (redsweep != 1 && redsweep != 2 && redsweep != 4)
312  return (PIX *)ERROR_PTR("redsweep not in {1,2,4}", procName, NULL);
313  if (sweeprange == 0.0)
314  sweeprange = DefaultSweepRange;
315  if (sweepdelta == 0.0)
316  sweepdelta = DefaultSweepDelta;
317  if (redsearch == 0)
318  redsearch = DefaultBsReduction;
319  else if (redsearch != 1 && redsearch != 2 && redsearch != 4)
320  return (PIX *)ERROR_PTR("redsearch not in {1,2,4}", procName, NULL);
321  if (thresh == 0)
322  thresh = DefaultBinaryThreshold;
323 
324  deg2rad = 3.1415926535 / 180.;
325 
326  /* Binarize if necessary */
327  depth = pixGetDepth(pixs);
328  if (depth == 1)
329  pixb = pixClone(pixs);
330  else
331  pixb = pixConvertTo1(pixs, thresh);
332 
333  /* Use the 1 bpp image to find the skew */
334  ret = pixFindSkewSweepAndSearch(pixb, &angle, &conf, redsweep, redsearch,
335  sweeprange, sweepdelta,
336  DefaultMinbsDelta);
337  pixDestroy(&pixb);
338  if (pangle) *pangle = angle;
339  if (pconf) *pconf = conf;
340  if (ret)
341  return pixClone(pixs);
342 
343  if (L_ABS(angle) < MinDeskewAngle || conf < MinAllowedConfidence)
344  return pixClone(pixs);
345 
346  if ((pixd = pixRotate(pixs, deg2rad * angle, L_ROTATE_AREA_MAP,
347  L_BRING_IN_WHITE, 0, 0)) == NULL)
348  return pixClone(pixs);
349  else
350  return pixd;
351 }
352 
353 
354 /*-----------------------------------------------------------------------*
355  * Simple top-level angle-finding interface *
356  *-----------------------------------------------------------------------*/
374 l_ok
376  l_float32 *pangle,
377  l_float32 *pconf)
378 {
379  PROCNAME("pixFindSkew");
380 
381  if (pangle) *pangle = 0.0;
382  if (pconf) *pconf = 0.0;
383  if (!pangle || !pconf)
384  return ERROR_INT("&angle and/or &conf not defined", procName, 1);
385  if (!pixs)
386  return ERROR_INT("pixs not defined", procName, 1);
387  if (pixGetDepth(pixs) != 1)
388  return ERROR_INT("pixs not 1 bpp", procName, 1);
389 
390  return pixFindSkewSweepAndSearch(pixs, pangle, pconf,
391  DefaultSweepReduction,
392  DefaultBsReduction,
393  DefaultSweepRange,
394  DefaultSweepDelta,
395  DefaultMinbsDelta);
396 }
397 
398 
399 /*-----------------------------------------------------------------------*
400  * Basic angle-finding functions *
401  *-----------------------------------------------------------------------*/
418 l_ok
420  l_float32 *pangle,
421  l_int32 reduction,
422  l_float32 sweeprange,
423  l_float32 sweepdelta)
424 {
425 l_int32 ret, bzero, i, nangles;
426 l_float32 deg2rad, theta;
427 l_float32 sum, maxscore, maxangle;
428 NUMA *natheta, *nascore;
429 PIX *pix, *pixt;
430 
431  PROCNAME("pixFindSkewSweep");
432 
433  if (!pangle)
434  return ERROR_INT("&angle not defined", procName, 1);
435  *pangle = 0.0;
436  if (!pixs)
437  return ERROR_INT("pixs not defined", procName, 1);
438  if (pixGetDepth(pixs) != 1)
439  return ERROR_INT("pixs not 1 bpp", procName, 1);
440  if (reduction != 1 && reduction != 2 && reduction != 4 && reduction != 8)
441  return ERROR_INT("reduction must be in {1,2,4,8}", procName, 1);
442 
443  deg2rad = 3.1415926535 / 180.;
444  ret = 0;
445 
446  /* Generate reduced image, if requested */
447  if (reduction == 1)
448  pix = pixClone(pixs);
449  else if (reduction == 2)
450  pix = pixReduceRankBinaryCascade(pixs, 1, 0, 0, 0);
451  else if (reduction == 4)
452  pix = pixReduceRankBinaryCascade(pixs, 1, 1, 0, 0);
453  else /* reduction == 8 */
454  pix = pixReduceRankBinaryCascade(pixs, 1, 1, 2, 0);
455 
456  pixZero(pix, &bzero);
457  if (bzero) {
458  pixDestroy(&pix);
459  return 1;
460  }
461 
462  nangles = (l_int32)((2. * sweeprange) / sweepdelta + 1);
463  natheta = numaCreate(nangles);
464  nascore = numaCreate(nangles);
465  pixt = pixCreateTemplate(pix);
466 
467  if (!pix || !pixt) {
468  ret = ERROR_INT("pix and pixt not both made", procName, 1);
469  goto cleanup;
470  }
471  if (!natheta || !nascore) {
472  ret = ERROR_INT("natheta and nascore not both made", procName, 1);
473  goto cleanup;
474  }
475 
476  for (i = 0; i < nangles; i++) {
477  theta = -sweeprange + i * sweepdelta; /* degrees */
478 
479  /* Shear pix about the UL corner and put the result in pixt */
480  pixVShearCorner(pixt, pix, deg2rad * theta, L_BRING_IN_WHITE);
481 
482  /* Get score */
483  pixFindDifferentialSquareSum(pixt, &sum);
484 
485 #if DEBUG_PRINT_SCORES
486  L_INFO("sum(%7.2f) = %7.0f\n", procName, theta, sum);
487 #endif /* DEBUG_PRINT_SCORES */
488 
489  /* Save the result in the output arrays */
490  numaAddNumber(nascore, sum);
491  numaAddNumber(natheta, theta);
492  }
493 
494  /* Find the location of the maximum (i.e., the skew angle)
495  * by fitting the largest data point and its two neighbors
496  * to a quadratic, using lagrangian interpolation. */
497  numaFitMax(nascore, &maxscore, natheta, &maxangle);
498  *pangle = maxangle;
499 
500 #if DEBUG_PRINT_SWEEP
501  L_INFO(" From sweep: angle = %7.3f, score = %7.3f\n", procName,
502  maxangle, maxscore);
503 #endif /* DEBUG_PRINT_SWEEP */
504 
505 #if DEBUG_PLOT_SCORES
506  /* Plot the result -- the scores versus rotation angle --
507  * using gnuplot with GPLOT_LINES (lines connecting data points).
508  * The GPLOT data structure is first created, with the
509  * appropriate data incorporated from the two input NUMAs,
510  * and then the function gplotMakeOutput() uses gnuplot to
511  * generate the output plot. This can be either a .png file
512  * or a .ps file, depending on whether you use GPLOT_PNG
513  * or GPLOT_PS. */
514  {GPLOT *gplot;
515  gplot = gplotCreate("sweep_output", GPLOT_PNG,
516  "Sweep. Variance of difference of ON pixels vs. angle",
517  "angle (deg)", "score");
518  gplotAddPlot(gplot, natheta, nascore, GPLOT_LINES, "plot1");
519  gplotAddPlot(gplot, natheta, nascore, GPLOT_POINTS, "plot2");
520  gplotMakeOutput(gplot);
521  gplotDestroy(&gplot);
522  }
523 #endif /* DEBUG_PLOT_SCORES */
524 
525 cleanup:
526  pixDestroy(&pix);
527  pixDestroy(&pixt);
528  numaDestroy(&nascore);
529  numaDestroy(&natheta);
530  return ret;
531 }
532 
533 
562 l_ok
564  l_float32 *pangle,
565  l_float32 *pconf,
566  l_int32 redsweep,
567  l_int32 redsearch,
568  l_float32 sweeprange,
569  l_float32 sweepdelta,
570  l_float32 minbsdelta)
571 {
572  return pixFindSkewSweepAndSearchScore(pixs, pangle, pconf, NULL,
573  redsweep, redsearch, 0.0, sweeprange,
574  sweepdelta, minbsdelta);
575 }
576 
577 
616 l_ok
618  l_float32 *pangle,
619  l_float32 *pconf,
620  l_float32 *pendscore,
621  l_int32 redsweep,
622  l_int32 redsearch,
623  l_float32 sweepcenter,
624  l_float32 sweeprange,
625  l_float32 sweepdelta,
626  l_float32 minbsdelta)
627 {
628  return pixFindSkewSweepAndSearchScorePivot(pixs, pangle, pconf, pendscore,
629  redsweep, redsearch, 0.0,
630  sweeprange, sweepdelta,
631  minbsdelta,
633 }
634 
635 
665 l_ok
667  l_float32 *pangle,
668  l_float32 *pconf,
669  l_float32 *pendscore,
670  l_int32 redsweep,
671  l_int32 redsearch,
672  l_float32 sweepcenter,
673  l_float32 sweeprange,
674  l_float32 sweepdelta,
675  l_float32 minbsdelta,
676  l_int32 pivot)
677 {
678 l_int32 ret, bzero, i, nangles, n, ratio, maxindex, minloc;
679 l_int32 width, height;
680 l_float32 deg2rad, theta, delta;
681 l_float32 sum, maxscore, maxangle;
682 l_float32 centerangle, leftcenterangle, rightcenterangle;
683 l_float32 lefttemp, righttemp;
684 l_float32 bsearchscore[5];
685 l_float32 minscore, minthresh;
686 l_float32 rangeleft;
687 NUMA *natheta, *nascore;
688 PIX *pixsw, *pixsch, *pixt1, *pixt2;
689 
690  PROCNAME("pixFindSkewSweepAndSearchScorePivot");
691 
692  if (pendscore) *pendscore = 0.0;
693  if (pangle) *pangle = 0.0;
694  if (pconf) *pconf = 0.0;
695  if (!pangle || !pconf)
696  return ERROR_INT("&angle and/or &conf not defined", procName, 1);
697  if (!pixs || pixGetDepth(pixs) != 1)
698  return ERROR_INT("pixs not defined or not 1 bpp", procName, 1);
699  if (redsweep != 1 && redsweep != 2 && redsweep != 4 && redsweep != 8)
700  return ERROR_INT("redsweep must be in {1,2,4,8}", procName, 1);
701  if (redsearch != 1 && redsearch != 2 && redsearch != 4 && redsearch != 8)
702  return ERROR_INT("redsearch must be in {1,2,4,8}", procName, 1);
703  if (redsearch > redsweep)
704  return ERROR_INT("redsearch must not exceed redsweep", procName, 1);
705  if (pivot != L_SHEAR_ABOUT_CORNER && pivot != L_SHEAR_ABOUT_CENTER)
706  return ERROR_INT("invalid pivot", procName, 1);
707 
708  deg2rad = 3.1415926535 / 180.;
709  ret = 0;
710 
711  /* Generate reduced image for binary search, if requested */
712  if (redsearch == 1)
713  pixsch = pixClone(pixs);
714  else if (redsearch == 2)
715  pixsch = pixReduceRankBinaryCascade(pixs, 1, 0, 0, 0);
716  else if (redsearch == 4)
717  pixsch = pixReduceRankBinaryCascade(pixs, 1, 1, 0, 0);
718  else /* redsearch == 8 */
719  pixsch = pixReduceRankBinaryCascade(pixs, 1, 1, 2, 0);
720 
721  pixZero(pixsch, &bzero);
722  if (bzero) {
723  pixDestroy(&pixsch);
724  return 1;
725  }
726 
727  /* Generate reduced image for sweep, if requested */
728  ratio = redsweep / redsearch;
729  if (ratio == 1) {
730  pixsw = pixClone(pixsch);
731  } else { /* ratio > 1 */
732  if (ratio == 2)
733  pixsw = pixReduceRankBinaryCascade(pixsch, 1, 0, 0, 0);
734  else if (ratio == 4)
735  pixsw = pixReduceRankBinaryCascade(pixsch, 1, 2, 0, 0);
736  else /* ratio == 8 */
737  pixsw = pixReduceRankBinaryCascade(pixsch, 1, 2, 2, 0);
738  }
739 
740  pixt1 = pixCreateTemplate(pixsw);
741  if (ratio == 1)
742  pixt2 = pixClone(pixt1);
743  else
744  pixt2 = pixCreateTemplate(pixsch);
745 
746  nangles = (l_int32)((2. * sweeprange) / sweepdelta + 1);
747  natheta = numaCreate(nangles);
748  nascore = numaCreate(nangles);
749 
750  if (!pixsch || !pixsw) {
751  ret = ERROR_INT("pixsch and pixsw not both made", procName, 1);
752  goto cleanup;
753  }
754  if (!pixt1 || !pixt2) {
755  ret = ERROR_INT("pixt1 and pixt2 not both made", procName, 1);
756  goto cleanup;
757  }
758  if (!natheta || !nascore) {
759  ret = ERROR_INT("natheta and nascore not both made", procName, 1);
760  goto cleanup;
761  }
762 
763  /* Do sweep */
764  rangeleft = sweepcenter - sweeprange;
765  for (i = 0; i < nangles; i++) {
766  theta = rangeleft + i * sweepdelta; /* degrees */
767 
768  /* Shear pix and put the result in pixt1 */
769  if (pivot == L_SHEAR_ABOUT_CORNER)
770  pixVShearCorner(pixt1, pixsw, deg2rad * theta, L_BRING_IN_WHITE);
771  else
772  pixVShearCenter(pixt1, pixsw, deg2rad * theta, L_BRING_IN_WHITE);
773 
774  /* Get score */
775  pixFindDifferentialSquareSum(pixt1, &sum);
776 
777 #if DEBUG_PRINT_SCORES
778  L_INFO("sum(%7.2f) = %7.0f\n", procName, theta, sum);
779 #endif /* DEBUG_PRINT_SCORES */
780 
781  /* Save the result in the output arrays */
782  numaAddNumber(nascore, sum);
783  numaAddNumber(natheta, theta);
784  }
785 
786  /* Find the largest of the set (maxscore at maxangle) */
787  numaGetMax(nascore, &maxscore, &maxindex);
788  numaGetFValue(natheta, maxindex, &maxangle);
789 
790 #if DEBUG_PRINT_SWEEP
791  L_INFO(" From sweep: angle = %7.3f, score = %7.3f\n", procName,
792  maxangle, maxscore);
793 #endif /* DEBUG_PRINT_SWEEP */
794 
795 #if DEBUG_PLOT_SCORES
796  /* Plot the sweep result -- the scores versus rotation angle --
797  * using gnuplot with GPLOT_LINES (lines connecting data points). */
798  {GPLOT *gplot;
799  gplot = gplotCreate("sweep_output", GPLOT_PNG,
800  "Sweep. Variance of difference of ON pixels vs. angle",
801  "angle (deg)", "score");
802  gplotAddPlot(gplot, natheta, nascore, GPLOT_LINES, "plot1");
803  gplotAddPlot(gplot, natheta, nascore, GPLOT_POINTS, "plot2");
804  gplotMakeOutput(gplot);
805  gplotDestroy(&gplot);
806  }
807 #endif /* DEBUG_PLOT_SCORES */
808 
809  /* Check if the max is at the end of the sweep. */
810  n = numaGetCount(natheta);
811  if (maxindex == 0 || maxindex == n - 1) {
812  L_WARNING("max found at sweep edge\n", procName);
813  goto cleanup;
814  }
815 
816  /* Empty the numas for re-use */
817  numaEmpty(nascore);
818  numaEmpty(natheta);
819 
820  /* Do binary search to find skew angle.
821  * First, set up initial three points. */
822  centerangle = maxangle;
823  if (pivot == L_SHEAR_ABOUT_CORNER) {
824  pixVShearCorner(pixt2, pixsch, deg2rad * centerangle, L_BRING_IN_WHITE);
825  pixFindDifferentialSquareSum(pixt2, &bsearchscore[2]);
826  pixVShearCorner(pixt2, pixsch, deg2rad * (centerangle - sweepdelta),
828  pixFindDifferentialSquareSum(pixt2, &bsearchscore[0]);
829  pixVShearCorner(pixt2, pixsch, deg2rad * (centerangle + sweepdelta),
831  pixFindDifferentialSquareSum(pixt2, &bsearchscore[4]);
832  } else {
833  pixVShearCenter(pixt2, pixsch, deg2rad * centerangle, L_BRING_IN_WHITE);
834  pixFindDifferentialSquareSum(pixt2, &bsearchscore[2]);
835  pixVShearCenter(pixt2, pixsch, deg2rad * (centerangle - sweepdelta),
837  pixFindDifferentialSquareSum(pixt2, &bsearchscore[0]);
838  pixVShearCenter(pixt2, pixsch, deg2rad * (centerangle + sweepdelta),
840  pixFindDifferentialSquareSum(pixt2, &bsearchscore[4]);
841  }
842 
843  numaAddNumber(nascore, bsearchscore[2]);
844  numaAddNumber(natheta, centerangle);
845  numaAddNumber(nascore, bsearchscore[0]);
846  numaAddNumber(natheta, centerangle - sweepdelta);
847  numaAddNumber(nascore, bsearchscore[4]);
848  numaAddNumber(natheta, centerangle + sweepdelta);
849 
850  /* Start the search */
851  delta = 0.5 * sweepdelta;
852  while (delta >= minbsdelta)
853  {
854  /* Get the left intermediate score */
855  leftcenterangle = centerangle - delta;
856  if (pivot == L_SHEAR_ABOUT_CORNER)
857  pixVShearCorner(pixt2, pixsch, deg2rad * leftcenterangle,
859  else
860  pixVShearCenter(pixt2, pixsch, deg2rad * leftcenterangle,
862  pixFindDifferentialSquareSum(pixt2, &bsearchscore[1]);
863  numaAddNumber(nascore, bsearchscore[1]);
864  numaAddNumber(natheta, leftcenterangle);
865 
866  /* Get the right intermediate score */
867  rightcenterangle = centerangle + delta;
868  if (pivot == L_SHEAR_ABOUT_CORNER)
869  pixVShearCorner(pixt2, pixsch, deg2rad * rightcenterangle,
871  else
872  pixVShearCenter(pixt2, pixsch, deg2rad * rightcenterangle,
874  pixFindDifferentialSquareSum(pixt2, &bsearchscore[3]);
875  numaAddNumber(nascore, bsearchscore[3]);
876  numaAddNumber(natheta, rightcenterangle);
877 
878  /* Find the maximum of the five scores and its location.
879  * Note that the maximum must be in the center
880  * three values, not in the end two. */
881  maxscore = bsearchscore[1];
882  maxindex = 1;
883  for (i = 2; i < 4; i++) {
884  if (bsearchscore[i] > maxscore) {
885  maxscore = bsearchscore[i];
886  maxindex = i;
887  }
888  }
889 
890  /* Set up score array to interpolate for the next iteration */
891  lefttemp = bsearchscore[maxindex - 1];
892  righttemp = bsearchscore[maxindex + 1];
893  bsearchscore[2] = maxscore;
894  bsearchscore[0] = lefttemp;
895  bsearchscore[4] = righttemp;
896 
897  /* Get new center angle and delta for next iteration */
898  centerangle = centerangle + delta * (maxindex - 2);
899  delta = 0.5 * delta;
900  }
901  *pangle = centerangle;
902 
903 #if DEBUG_PRINT_SCORES
904  L_INFO(" Binary search score = %7.3f\n", procName, bsearchscore[2]);
905 #endif /* DEBUG_PRINT_SCORES */
906 
907  if (pendscore) /* save if requested */
908  *pendscore = bsearchscore[2];
909 
910  /* Return the ratio of Max score over Min score
911  * as a confidence value. Don't trust if the Min score
912  * is too small, which can happen if the image is all black
913  * with only a few white pixels interspersed. In that case,
914  * we get a contribution from the top and bottom edges when
915  * vertically sheared, but this contribution becomes zero when
916  * the shear angle is zero. For zero shear angle, the only
917  * contribution will be from the white pixels. We expect that
918  * the signal goes as the product of the (height * width^2),
919  * so we compute a (hopefully) normalized minimum threshold as
920  * a function of these dimensions. */
921  numaGetMin(nascore, &minscore, &minloc);
922  width = pixGetWidth(pixsch);
923  height = pixGetHeight(pixsch);
924  minthresh = MinscoreThreshFactor * width * width * height;
925 
926 #if DEBUG_THRESHOLD
927  L_INFO(" minthresh = %10.2f, minscore = %10.2f\n", procName,
928  minthresh, minscore);
929  L_INFO(" maxscore = %10.2f\n", procName, maxscore);
930 #endif /* DEBUG_THRESHOLD */
931 
932  if (minscore > minthresh)
933  *pconf = maxscore / minscore;
934  else
935  *pconf = 0.0;
936 
937  /* Don't trust it if too close to the edge of the sweep
938  * range or if maxscore is small */
939  if ((centerangle > rangeleft + 2 * sweeprange - sweepdelta) ||
940  (centerangle < rangeleft + sweepdelta) ||
941  (maxscore < MinValidMaxscore))
942  *pconf = 0.0;
943 
944 #if DEBUG_PRINT_BINARY
945  lept_stderr("Binary search: angle = %7.3f, score ratio = %6.2f\n",
946  *pangle, *pconf);
947  lept_stderr(" max score = %8.0f\n", maxscore);
948 #endif /* DEBUG_PRINT_BINARY */
949 
950 #if DEBUG_PLOT_SCORES
951  /* Plot the result -- the scores versus rotation angle --
952  * using gnuplot with GPLOT_POINTS. Because the data
953  * points are not ordered by theta (increasing or decreasing),
954  * using GPLOT_LINES would be confusing! */
955  {GPLOT *gplot;
956  gplot = gplotCreate("search_output", GPLOT_PNG,
957  "Binary search. Variance of difference of ON pixels vs. angle",
958  "angle (deg)", "score");
959  gplotAddPlot(gplot, natheta, nascore, GPLOT_POINTS, "plot1");
960  gplotMakeOutput(gplot);
961  gplotDestroy(&gplot);
962  }
963 #endif /* DEBUG_PLOT_SCORES */
964 
965 cleanup:
966  pixDestroy(&pixsw);
967  pixDestroy(&pixsch);
968  pixDestroy(&pixt1);
969  pixDestroy(&pixt2);
970  numaDestroy(&nascore);
971  numaDestroy(&natheta);
972  return ret;
973 }
974 
975 
976 /*---------------------------------------------------------------------*
977  * Search over arbitrary range of angles in orthogonal directions *
978  *---------------------------------------------------------------------*/
979 /*
980  * \brief pixFindSkewOrthogonalRange()
981  *
982  * \param[in] pixs 1 bpp
983  * \param[out] pangle angle required to deskew; in degrees cw
984  * \param[out] pconf confidence given by ratio of max/min score
985  * \param[in] redsweep sweep reduction factor = 1, 2, 4 or 8
986  * \param[in] redsearch binary search reduction factor = 1, 2, 4 or 8;
987  * and must not exceed redsweep
988  * \param[in] sweeprange half the full range in each orthogonal
989  * direction, taken about 0, in degrees
990  * \param[in] sweepdelta angle increment of sweep; in degrees
991  * \param[in] minbsdelta min binary search increment angle; in degrees
992  * \param[in] confprior amount by which confidence of 90 degree rotated
993  * result is reduced when comparing with unrotated
994  * confidence value
995  * \return 0 if OK, 1 on error or if angle measurement not valid
996  *
997  * <pre>
998  * Notes:
999  * (1) This searches for the skew angle, first in the range
1000  * [-sweeprange, sweeprange], and then in
1001  * [90 - sweeprange, 90 + sweeprange], with angles measured
1002  * clockwise. For exploring the full range of possibilities,
1003  * suggest using sweeprange = 47.0 degrees, giving some overlap
1004  * at 45 and 135 degrees. From these results, and discounting
1005  * the the second confidence by %confprior, it selects the
1006  * angle for maximal differential variance. If the angle
1007  * is larger than pi/4, the angle found after 90 degree rotation
1008  * is selected.
1009  * (2) The larger the confidence value, the greater the probability
1010  * that the proper alignment is given by the angle that maximizes
1011  * variance. It should be compared to a threshold, which depends
1012  * on the application. Values between 3.0 and 6.0 are common.
1013  * (3) Allowing for both portrait and landscape searches is more
1014  * difficult, because if the signal from the text lines is weak,
1015  * a signal from vertical rules can be larger!
1016  * The most difficult documents to deskew have some or all of:
1017  * (a) Multiple columns, not aligned
1018  * (b) Black lines along the vertical edges
1019  * (c) Text from two pages, and at different angles
1020  * Rule of thumb for resolution:
1021  * (a) If the margins are clean, you can work at 75 ppi,
1022  * although 100 ppi is safer.
1023  * (b) If there are vertical lines in the margins, do not
1024  * work below 150 ppi. The signal from the text lines must
1025  * exceed that from the margin lines.
1026  * (4) Choosing the %confprior parameter depends on knowing something
1027  * about the source of image. However, we're not using
1028  * real probabilities here, so its use is qualitative.
1029  * If landscape and portrait are equally likely, use
1030  * %confprior = 0.0. If the likelihood of portrait (non-rotated)
1031  * is 100 times higher than that of landscape, we want to reduce
1032  * the chance that we rotate to landscape in a situation where
1033  * the landscape signal is accidentally larger than the
1034  * portrait signal. To do this use a positive value of
1035  * %confprior; say 1.5.
1036  * </pre>
1037  */
1038 l_int32
1039 pixFindSkewOrthogonalRange(PIX *pixs,
1040  l_float32 *pangle,
1041  l_float32 *pconf,
1042  l_int32 redsweep,
1043  l_int32 redsearch,
1044  l_float32 sweeprange,
1045  l_float32 sweepdelta,
1046  l_float32 minbsdelta,
1047  l_float32 confprior)
1048 {
1049 l_float32 angle1, conf1, score1, angle2, conf2, score2;
1050 PIX *pixr;
1051 
1052  PROCNAME("pixFindSkewOrthogonalRange");
1053 
1054  if (pangle) *pangle = 0.0;
1055  if (pconf) *pconf = 0.0;
1056  if (!pangle || !pconf)
1057  return ERROR_INT("&angle and/or &conf not defined", procName, 1);
1058  if (!pixs || pixGetDepth(pixs) != 1)
1059  return ERROR_INT("pixs not defined or not 1 bpp", procName, 1);
1060 
1061  pixFindSkewSweepAndSearchScorePivot(pixs, &angle1, &conf1, &score1,
1062  redsweep, redsearch, 0.0,
1063  sweeprange, sweepdelta, minbsdelta,
1065  pixr = pixRotateOrth(pixs, 1);
1066  pixFindSkewSweepAndSearchScorePivot(pixr, &angle2, &conf2, &score2,
1067  redsweep, redsearch, 0.0,
1068  sweeprange, sweepdelta, minbsdelta,
1070  pixDestroy(&pixr);
1071 
1072  if (conf1 > conf2 - confprior) {
1073  *pangle = angle1;
1074  *pconf = conf1;
1075  } else {
1076  *pangle = -90.0 + angle2;
1077  *pconf = conf2;
1078  }
1079 
1080 #if DEBUG_PRINT_ORTH
1081  lept_stderr(" About 0: angle1 = %7.3f, conf1 = %7.3f, score1 = %f\n",
1082  angle1, conf1, score1);
1083  lept_stderr(" About 90: angle2 = %7.3f, conf2 = %7.3f, score2 = %f\n",
1084  angle2, conf2, score2);
1085  lept_stderr(" Final: angle = %7.3f, conf = %7.3f\n", *pangle, *pconf);
1086 #endif /* DEBUG_PRINT_ORTH */
1087 
1088  return 0;
1089 }
1090 
1091 
1092 
1093 /*----------------------------------------------------------------*
1094  * Differential square sum function *
1095  *----------------------------------------------------------------*/
1111 l_ok
1113  l_float32 *psum)
1114 {
1115 l_int32 i, n;
1116 l_int32 w, h, skiph, skip, nskip;
1117 l_float32 val1, val2, diff, sum;
1118 NUMA *na;
1119 
1120  PROCNAME("pixFindDifferentialSquareSum");
1121 
1122  if (!psum)
1123  return ERROR_INT("&sum not defined", procName, 1);
1124  *psum = 0.0;
1125  if (!pixs)
1126  return ERROR_INT("pixs not defined", procName, 1);
1127 
1128  /* Generate a number array consisting of the sum
1129  * of pixels in each row of pixs */
1130  if ((na = pixCountPixelsByRow(pixs, NULL)) == NULL)
1131  return ERROR_INT("na not made", procName, 1);
1132 
1133  /* Compute the number of rows at top and bottom to omit.
1134  * We omit these to avoid getting a spurious signal from
1135  * the top and bottom of a (nearly) all black image. */
1136  w = pixGetWidth(pixs);
1137  h = pixGetHeight(pixs);
1138  skiph = (l_int32)(0.05 * w); /* skip for max shear of 0.025 radians */
1139  skip = L_MIN(h / 10, skiph); /* don't remove more than 10% of image */
1140  nskip = L_MAX(skip / 2, 1); /* at top & bot; skip at least one line */
1141 
1142  /* Sum the squares of differential row sums, on the
1143  * allowed rows. Note that nskip must be >= 1. */
1144  n = numaGetCount(na);
1145  sum = 0.0;
1146  for (i = nskip; i < n - nskip; i++) {
1147  numaGetFValue(na, i - 1, &val1);
1148  numaGetFValue(na, i, &val2);
1149  diff = val2 - val1;
1150  sum += diff * diff;
1151  }
1152  numaDestroy(&na);
1153  *psum = sum;
1154  return 0;
1155 }
1156 
1157 
1158 /*----------------------------------------------------------------*
1159  * Normalized square sum *
1160  *----------------------------------------------------------------*/
1184 l_ok
1186  l_float32 *phratio,
1187  l_float32 *pvratio,
1188  l_float32 *pfract)
1189 {
1190 l_int32 i, w, h, empty;
1191 l_float32 sum, sumsq, uniform, val;
1192 NUMA *na;
1193 PIX *pixt;
1194 
1195  PROCNAME("pixFindNormalizedSquareSum");
1196 
1197  if (phratio) *phratio = 0.0;
1198  if (pvratio) *pvratio = 0.0;
1199  if (pfract) *pfract = 0.0;
1200  if (!phratio && !pvratio)
1201  return ERROR_INT("nothing to do", procName, 1);
1202  if (!pixs || pixGetDepth(pixs) != 1)
1203  return ERROR_INT("pixs not defined or not 1 bpp", procName, 1);
1204  pixGetDimensions(pixs, &w, &h, NULL);
1205 
1206  empty = 0;
1207  if (phratio) {
1208  na = pixCountPixelsByRow(pixs, NULL);
1209  numaGetSum(na, &sum); /* fg pixels */
1210  if (pfract) *pfract = sum / (l_float32)(w * h);
1211  if (sum != 0.0) {
1212  uniform = sum * sum / h; /* h*(sum / h)^2 */
1213  sumsq = 0.0;
1214  for (i = 0; i < h; i++) {
1215  numaGetFValue(na, i, &val);
1216  sumsq += val * val;
1217  }
1218  *phratio = sumsq / uniform;
1219  } else {
1220  empty = 1;
1221  }
1222  numaDestroy(&na);
1223  }
1224 
1225  if (pvratio) {
1226  if (empty == 1) return 1;
1227  pixt = pixRotateOrth(pixs, 1);
1228  na = pixCountPixelsByRow(pixt, NULL);
1229  numaGetSum(na, &sum);
1230  if (pfract) *pfract = sum / (l_float32)(w * h);
1231  if (sum != 0.0) {
1232  uniform = sum * sum / w;
1233  sumsq = 0.0;
1234  for (i = 0; i < w; i++) {
1235  numaGetFValue(na, i, &val);
1236  sumsq += val * val;
1237  }
1238  *pvratio = sumsq / uniform;
1239  } else {
1240  empty = 1;
1241  }
1242  pixDestroy(&pixt);
1243  numaDestroy(&na);
1244  }
1245 
1246  return empty;
1247 }
l_ok numaGetSum(NUMA *na, l_float32 *psum)
numaGetSum()
Definition: numafunc1.c:538
void gplotDestroy(GPLOT **pgplot)
gplotDestroy()
Definition: gplot.c:255
NUMA * pixCountPixelsByRow(PIX *pix, l_int32 *tab8)
pixCountPixelsByRow()
Definition: pix3.c:2143
PIX * pixConvertTo1(PIX *pixs, l_int32 threshold)
pixConvertTo1()
Definition: pixconv.c:3026
l_ok numaGetFValue(NUMA *na, l_int32 index, l_float32 *pval)
numaGetFValue()
Definition: numabasic.c:719
l_ok numaFitMax(NUMA *na, l_float32 *pmaxval, NUMA *naloc, l_float32 *pmaxloc)
numaFitMax()
Definition: numafunc1.c:2162
PIX * pixDeskew(PIX *pixs, l_int32 redsearch)
pixDeskew()
Definition: skew.c:210
l_ok pixFindDifferentialSquareSum(PIX *pixs, l_float32 *psum)
pixFindDifferentialSquareSum()
Definition: skew.c:1112
PIX * pixCreateTemplate(const PIX *pixs)
pixCreateTemplate()
Definition: pix1.c:383
l_ok numaAddNumber(NUMA *na, l_float32 val)
numaAddNumber()
Definition: numabasic.c:478
l_ok gplotMakeOutput(GPLOT *gplot)
gplotMakeOutput()
Definition: gplot.c:466
GPLOT * gplotCreate(const char *rootname, l_int32 outformat, const char *title, const char *xlabel, const char *ylabel)
gplotCreate()
Definition: gplot.c:187
PIX * pixDeskewBoth(PIX *pixs, l_int32 redsearch)
pixDeskewBoth()
Definition: skew.c:167
void lept_stderr(const char *fmt,...)
lept_stderr()
Definition: utils1.c:306
l_ok pixFindSkewSweep(PIX *pixs, l_float32 *pangle, l_int32 reduction, l_float32 sweeprange, l_float32 sweepdelta)
pixFindSkewSweep()
Definition: skew.c:419
l_ok pixFindSkewSweepAndSearchScore(PIX *pixs, l_float32 *pangle, l_float32 *pconf, l_float32 *pendscore, l_int32 redsweep, l_int32 redsearch, l_float32 sweepcenter, l_float32 sweeprange, l_float32 sweepdelta, l_float32 minbsdelta)
pixFindSkewSweepAndSearchScore()
Definition: skew.c:617
PIX * pixFindSkewAndDeskew(PIX *pixs, l_int32 redsearch, l_float32 *pangle, l_float32 *pconf)
pixFindSkewAndDeskew()
Definition: skew.c:246
l_ok pixFindSkew(PIX *pixs, l_float32 *pangle, l_float32 *pconf)
pixFindSkew()
Definition: skew.c:375
NUMA * numaCreate(l_int32 n)
numaCreate()
Definition: numabasic.c:194
PIX * pixVShearCenter(PIX *pixd, PIX *pixs, l_float32 radang, l_int32 incolor)
pixVShearCenter()
Definition: shear.c:433
PIX * pixVShearCorner(PIX *pixd, PIX *pixs, l_float32 radang, l_int32 incolor)
pixVShearCorner()
Definition: shear.c:371
Definition: array.h:70
l_ok pixFindSkewSweepAndSearchScorePivot(PIX *pixs, l_float32 *pangle, l_float32 *pconf, l_float32 *pendscore, l_int32 redsweep, l_int32 redsearch, l_float32 sweepcenter, l_float32 sweeprange, l_float32 sweepdelta, l_float32 minbsdelta, l_int32 pivot)
pixFindSkewSweepAndSearchScorePivot()
Definition: skew.c:666
l_int32 numaGetCount(NUMA *na)
numaGetCount()
Definition: numabasic.c:658
l_ok pixFindNormalizedSquareSum(PIX *pixs, l_float32 *phratio, l_float32 *pvratio, l_float32 *pfract)
pixFindNormalizedSquareSum()
Definition: skew.c:1185
PIX * pixDeskewGeneral(PIX *pixs, l_int32 redsweep, l_float32 sweeprange, l_float32 sweepdelta, l_int32 redsearch, l_int32 thresh, l_float32 *pangle, l_float32 *pconf)
pixDeskewGeneral()
Definition: skew.c:290
Definition: gplot.h:76
l_ok gplotAddPlot(GPLOT *gplot, NUMA *nax, NUMA *nay, l_int32 plotstyle, const char *plotlabel)
gplotAddPlot()
Definition: gplot.c:320
PIX * pixClone(PIX *pixs)
pixClone()
Definition: pix1.c:593
l_ok numaEmpty(NUMA *na)
numaEmpty()
Definition: numabasic.c:454
void pixDestroy(PIX **ppix)
pixDestroy()
Definition: pix1.c:621
PIX * pixRotateOrth(PIX *pixs, l_int32 quads)
pixRotateOrth()
Definition: rotateorth.c:75
void numaDestroy(NUMA **pna)
numaDestroy()
Definition: numabasic.c:366
l_ok pixGetDimensions(const PIX *pix, l_int32 *pw, l_int32 *ph, l_int32 *pd)
pixGetDimensions()
Definition: pix1.c:1113
l_ok numaGetMin(NUMA *na, l_float32 *pminval, l_int32 *piminloc)
numaGetMin()
Definition: numafunc1.c:453
Definition: pix.h:138
l_ok pixZero(PIX *pix, l_int32 *pempty)
pixZero()
Definition: pix3.c:1815
l_ok pixFindSkewSweepAndSearch(PIX *pixs, l_float32 *pangle, l_float32 *pconf, l_int32 redsweep, l_int32 redsearch, l_float32 sweeprange, l_float32 sweepdelta, l_float32 minbsdelta)
pixFindSkewSweepAndSearch()
Definition: skew.c:563
PIX * pixRotate90(PIX *pixs, l_int32 direction)
pixRotate90()
Definition: rotateorth.c:166
l_ok numaGetMax(NUMA *na, l_float32 *pmaxval, l_int32 *pimaxloc)
numaGetMax()
Definition: numafunc1.c:496
PIX * pixReduceRankBinaryCascade(PIX *pixs, l_int32 level1, l_int32 level2, l_int32 level3, l_int32 level4)
pixReduceRankBinaryCascade()
Definition: binreduce.c:152
PIX * pixRotate(PIX *pixs, l_float32 angle, l_int32 type, l_int32 incolor, l_int32 width, l_int32 height)
pixRotate()
Definition: rotate.c:101