Leptonica 1.68
C Image Processing Library

convolvelow.c

Go to the documentation of this file.
00001 /*====================================================================*
00002  -  Copyright (C) 2001 Leptonica.  All rights reserved.
00003  -  This software is distributed in the hope that it will be
00004  -  useful, but with NO WARRANTY OF ANY KIND.
00005  -  No author or distributor accepts responsibility to anyone for the
00006  -  consequences of using this software, or for whether it serves any
00007  -  particular purpose or works at all, unless he or she says so in
00008  -  writing.  Everyone is granted permission to copy, modify and
00009  -  redistribute this source code, for commercial or non-commercial
00010  -  purposes, with the following restrictions: (1) the origin of this
00011  -  source code must not be misrepresented; (2) modified versions must
00012  -  be plainly marked as such; and (3) this notice may not be removed
00013  -  or altered from any source or modified source distribution.
00014  *====================================================================*/
00015 
00016 
00017 /*
00018  *  convolvelow.c
00019  *
00020  *      Grayscale block convolution
00021  *          void      blockconvLow()
00022  *          void      blockconvAccumLow()
00023  *
00024  *      Binary block sum and rank filter
00025  *          void      blocksumLow()
00026  */
00027 
00028 #include <stdio.h>
00029 #include "allheaders.h"
00030 
00031 
00032 /*----------------------------------------------------------------------*
00033  *                     Grayscale Block Convolution                      *
00034  *----------------------------------------------------------------------*/
00035 /*!
00036  *  blockconvLow()
00037  *
00038  *      Input:  data   (data of input image, to be convolved)
00039  *              w, h, wpl
00040  *              dataa    (data of 32 bpp accumulator)
00041  *              wpla     (accumulator)
00042  *              wc      (convolution "half-width")
00043  *              hc      (convolution "half-height")
00044  *      Return: void
00045  *
00046  *  Notes:
00047  *      (1) The full width and height of the convolution kernel
00048  *          are (2 * wc + 1) and (2 * hc + 1).
00049  *      (2) The lack of symmetry between the handling of the
00050  *          first (hc + 1) lines and the last (hc) lines,
00051  *          and similarly with the columns, is due to fact that
00052  *          for the pixel at (x,y), the accumulator values are
00053  *          taken at (x + wc, y + hc), (x - wc - 1, y + hc),
00054  *          (x + wc, y - hc - 1) and (x - wc - 1, y - hc - 1).
00055  *      (3) We compute sums, normalized as if there were no reduced
00056  *          area at the boundary.  This under-estimates the value
00057  *          of the boundary pixels, so we multiply them by another
00058  *          normalization factor that is greater than 1.
00059  *      (4) This second normalization is done first for the first
00060  *          hc + 1 lines; then for the last hc lines; and finally
00061  *          for the first wc + 1 and last wc columns in the intermediate
00062  *          lines.
00063  *      (5) The caller should verify that wc < w and hc < h.
00064  *          Under those conditions, illegal reads and writes can occur.
00065  *      (6) Implementation note: to get the same results in the interior
00066  *          between this function and pixConvolve(), it is necessary to
00067  *          add 0.5 for roundoff in the main loop that runs over all pixels.
00068  *          However, if we do that and have white (255) pixels near the
00069  *          image boundary, some overflow occurs for pixels very close
00070  *          to the boundary.  We can't fix this by subtracting from the
00071  *          normalized values for the boundary pixels, because this results
00072  *          in underflow if the boundary pixels are black (0).  Empirically,
00073  *          adding 0.25 (instead of 0.5) before truncating in the main
00074  *          loop will not cause overflow, but this gives some
00075  *          off-by-1-level errors in interior pixel values.  So we add
00076  *          0.5 for roundoff in the main loop, and for pixels within a
00077  *          half filter width of the boundary, use a L_MIN of the
00078  *          computed value and 255 to avoid overflow during normalization.
00079  */
00080 void
00081 blockconvLow(l_uint32  *data,
00082              l_int32    w,
00083              l_int32    h,
00084              l_int32    wpl,
00085              l_uint32  *dataa,
00086              l_int32    wpla,
00087              l_int32    wc,
00088              l_int32    hc)
00089 {
00090 l_int32    i, j, imax, imin, jmax, jmin;
00091 l_int32    wn, hn, fwc, fhc, wmwc, hmhc;
00092 l_float32  norm, normh, normw;
00093 l_uint32   val;
00094 l_uint32  *linemina, *linemaxa, *line;
00095 
00096     PROCNAME("blockconvLow");
00097 
00098     wmwc = w - wc;
00099     hmhc = h - hc;
00100     if (wmwc <= 0 || hmhc <= 0) {
00101         L_ERROR("wc >= w || hc >=h", procName);
00102         return;
00103     }
00104     fwc = 2 * wc + 1;
00105     fhc = 2 * hc + 1;
00106     norm = 1. / (fwc * fhc);
00107 
00108         /*------------------------------------------------------------*
00109          *  compute, using b.c. only to set limits on the accum image *
00110          *------------------------------------------------------------*/
00111     for (i = 0; i < h; i++) {
00112         imin = L_MAX(i - 1 - hc, 0);
00113         imax = L_MIN(i + hc, h - 1);
00114         line = data + wpl * i;
00115         linemina = dataa + wpla * imin;
00116         linemaxa = dataa + wpla * imax;
00117         for (j = 0; j < w; j++) {
00118             jmin = L_MAX(j - 1 - wc, 0);
00119             jmax = L_MIN(j + wc, w - 1);
00120             val = linemaxa[jmax] - linemaxa[jmin]
00121                   + linemina[jmin] - linemina[jmax];
00122             val = (l_uint8)(norm * val + 0.5);  /* see comment above */
00123             SET_DATA_BYTE(line, j, val);
00124         }
00125     }
00126 
00127         /*------------------------------------------------------------*
00128          *          now fix normalization for boundary pixels         *
00129          *------------------------------------------------------------*/
00130     for (i = 0; i <= hc; i++) {    /* first hc + 1 lines */
00131         hn = hc + i;
00132         normh = (l_float32)fhc / (l_float32)hn;   /* > 1 */
00133         line = data + wpl * i;
00134         for (j = 0; j <= wc; j++) {
00135             wn = wc + j;
00136             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00137             val = GET_DATA_BYTE(line, j);
00138             val = (l_uint8)L_MIN(val * normh * normw, 255);
00139             SET_DATA_BYTE(line, j, val);
00140         }
00141         for (j = wc + 1; j < wmwc; j++) {
00142             val = GET_DATA_BYTE(line, j);
00143             val = (l_uint8)L_MIN(val * normh, 255);
00144             SET_DATA_BYTE(line, j, val);
00145         }
00146         for (j = wmwc; j < w; j++) {
00147             wn = wc + w - j;
00148             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00149             val = GET_DATA_BYTE(line, j);
00150             val = (l_uint8)L_MIN(val * normh * normw, 255);
00151             SET_DATA_BYTE(line, j, val);
00152         }
00153     }
00154 
00155     for (i = hmhc; i < h; i++) {  /* last hc lines */
00156         hn = hc + h - i;
00157         normh = (l_float32)fhc / (l_float32)hn;   /* > 1 */
00158         line = data + wpl * i;
00159         for (j = 0; j <= wc; j++) {
00160             wn = wc + j;
00161             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00162             val = GET_DATA_BYTE(line, j);
00163             val = (l_uint8)L_MIN(val * normh * normw, 255);
00164             SET_DATA_BYTE(line, j, val);
00165         }
00166         for (j = wc + 1; j < wmwc; j++) {
00167             val = GET_DATA_BYTE(line, j);
00168             val = (l_uint8)L_MIN(val * normh, 255);
00169             SET_DATA_BYTE(line, j, val);
00170         }
00171         for (j = wmwc; j < w; j++) {
00172             wn = wc + w - j;
00173             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00174             val = GET_DATA_BYTE(line, j);
00175             val = (l_uint8)L_MIN(val * normh * normw, 255);
00176             SET_DATA_BYTE(line, j, val);
00177         }
00178     }
00179 
00180     for (i = hc + 1; i < hmhc; i++) {    /* intermediate lines */
00181         line = data + wpl * i;
00182         for (j = 0; j <= wc; j++) {   /* first wc + 1 columns */
00183             wn = wc + j;
00184             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00185             val = GET_DATA_BYTE(line, j);
00186             val = (l_uint8)L_MIN(val * normw, 255);
00187             SET_DATA_BYTE(line, j, val);
00188         }
00189         for (j = wmwc; j < w; j++) {   /* last wc columns */
00190             wn = wc + w - j;
00191             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00192             val = GET_DATA_BYTE(line, j);
00193             val = (l_uint8)L_MIN(val * normw, 255);
00194             SET_DATA_BYTE(line, j, val);
00195         }
00196     }
00197 
00198     return;
00199 }
00200 
00201 
00202 
00203 /*
00204  *  blockconvAccumLow()
00205  *
00206  *      Input:  datad  (32 bpp dest)
00207  *              w, h, wpld (of 32 bpp dest)
00208  *              datas (1, 8 or 32 bpp src)
00209  *              d (bpp of src)
00210  *              wpls (of src)
00211  *      Return: void
00212  *
00213  *  Notes:
00214  *      (1) The general recursion relation is
00215  *             a(i,j) = v(i,j) + a(i-1, j) + a(i, j-1) - a(i-1, j-1)
00216  *          For the first line, this reduces to the special case
00217  *             a(i,j) = v(i,j) + a(i, j-1)
00218  *          For the first column, the special case is
00219  *             a(i,j) = v(i,j) + a(i-1, j)
00220  */
00221 void
00222 blockconvAccumLow(l_uint32  *datad,
00223                   l_int32    w,
00224                   l_int32    h,
00225                   l_int32    wpld,
00226                   l_uint32  *datas,
00227                   l_int32    d,
00228                   l_int32    wpls)
00229 {
00230 l_uint8    val;
00231 l_int32    i, j;
00232 l_uint32   val32;
00233 l_uint32  *lines, *lined, *linedp;
00234 
00235     PROCNAME("blockconvAccumLow");
00236 
00237     lines = datas;
00238     lined = datad;
00239 
00240     if (d == 1) {
00241             /* Do the first line */
00242         for (j = 0; j < w; j++) {
00243             val = GET_DATA_BIT(lines, j);
00244             if (j == 0)
00245                 lined[0] = val;
00246             else
00247                 lined[j] = lined[j - 1] + val;
00248         }
00249 
00250             /* Do the other lines */
00251         for (i = 1; i < h; i++) {
00252             lines = datas + i * wpls;
00253             lined = datad + i * wpld;  /* curr dest line */
00254             linedp = lined - wpld;   /* prev dest line */
00255             for (j = 0; j < w; j++) {
00256                 val = GET_DATA_BIT(lines, j);
00257                 if (j == 0)
00258                     lined[0] = val + linedp[0];
00259                 else 
00260                     lined[j] = val + lined[j - 1] + linedp[j] - linedp[j - 1];
00261             }
00262         }
00263     }
00264     else if (d == 8) {
00265             /* Do the first line */
00266         for (j = 0; j < w; j++) {
00267             val = GET_DATA_BYTE(lines, j);
00268             if (j == 0)
00269                 lined[0] = val;
00270             else
00271                 lined[j] = lined[j - 1] + val;
00272         }
00273 
00274             /* Do the other lines */
00275         for (i = 1; i < h; i++) {
00276             lines = datas + i * wpls;
00277             lined = datad + i * wpld;  /* curr dest line */
00278             linedp = lined - wpld;   /* prev dest line */
00279             for (j = 0; j < w; j++) {
00280                 val = GET_DATA_BYTE(lines, j);
00281                 if (j == 0)
00282                     lined[0] = val + linedp[0];
00283                 else 
00284                     lined[j] = val + lined[j - 1] + linedp[j] - linedp[j - 1];
00285             }
00286         }
00287     }
00288     else if (d == 32) {
00289             /* Do the first line */
00290         for (j = 0; j < w; j++) {
00291             val32 = lines[j];
00292             if (j == 0)
00293                 lined[0] = val32;
00294             else
00295                 lined[j] = lined[j - 1] + val32;
00296         }
00297 
00298             /* Do the other lines */
00299         for (i = 1; i < h; i++) {
00300             lines = datas + i * wpls;
00301             lined = datad + i * wpld;  /* curr dest line */
00302             linedp = lined - wpld;   /* prev dest line */
00303             for (j = 0; j < w; j++) {
00304                 val32 = lines[j];
00305                 if (j == 0)
00306                     lined[0] = val32 + linedp[0];
00307                 else 
00308                     lined[j] = val32 + lined[j - 1] + linedp[j] - linedp[j - 1];
00309             }
00310         }
00311     }
00312     else
00313         L_ERROR("depth not 1, 8 or 32 bpp", procName);
00314 
00315     return;
00316 }
00317 
00318 
00319 /*----------------------------------------------------------------------*
00320  *                        Binary Block Sum/Rank                         *
00321  *----------------------------------------------------------------------*/
00322 /*!
00323  *  blocksumLow()
00324  *
00325  *      Input:  datad  (of 8 bpp dest)
00326  *              w, h, wpl  (of 8 bpp dest)
00327  *              dataa (of 32 bpp accum)
00328  *              wpla  (of 32 bpp accum)
00329  *              wc, hc  (convolution "half-width" and "half-height")
00330  *      Return: void
00331  *
00332  *  Notes:
00333  *      (1) The full width and height of the convolution kernel
00334  *          are (2 * wc + 1) and (2 * hc + 1).
00335  *      (2) The lack of symmetry between the handling of the
00336  *          first (hc + 1) lines and the last (hc) lines,
00337  *          and similarly with the columns, is due to fact that
00338  *          for the pixel at (x,y), the accumulator values are
00339  *          taken at (x + wc, y + hc), (x - wc - 1, y + hc),
00340  *          (x + wc, y - hc - 1) and (x - wc - 1, y - hc - 1).
00341  *      (3) Compute sums of ON pixels within the block filter size,
00342  *          normalized between 0 and 255, as if there were no reduced
00343  *          area at the boundary.  This under-estimates the value
00344  *          of the boundary pixels, so we multiply them by another
00345  *          normalization factor that is greater than 1.
00346  *      (4) This second normalization is done first for the first
00347  *          hc + 1 lines; then for the last hc lines; and finally
00348  *          for the first wc + 1 and last wc columns in the intermediate
00349  *          lines.
00350  *      (5) The caller should verify that wc < w and hc < h.
00351  *          Under those conditions, illegal reads and writes can occur.
00352  */
00353 void
00354 blocksumLow(l_uint32  *datad,
00355             l_int32    w,
00356             l_int32    h,
00357             l_int32    wpl,
00358             l_uint32  *dataa,
00359             l_int32    wpla,
00360             l_int32    wc,
00361             l_int32    hc)
00362 {
00363 l_int32    i, j, imax, imin, jmax, jmin;
00364 l_int32    wn, hn, fwc, fhc, wmwc, hmhc;
00365 l_float32  norm, normh, normw;
00366 l_uint32   val;
00367 l_uint32  *linemina, *linemaxa, *lined;
00368 
00369     PROCNAME("blocksumLow");
00370 
00371     wmwc = w - wc;
00372     hmhc = h - hc;
00373     if (wmwc <= 0 || hmhc <= 0) {
00374         L_ERROR("wc >= w || hc >=h", procName);
00375         return;
00376     }
00377     fwc = 2 * wc + 1;
00378     fhc = 2 * hc + 1;
00379     norm = 255. / (fwc * fhc);
00380 
00381         /*------------------------------------------------------------*
00382          *  compute, using b.c. only to set limits on the accum image *
00383          *------------------------------------------------------------*/
00384     for (i = 0; i < h; i++) {
00385         imin = L_MAX(i - 1 - hc, 0);
00386         imax = L_MIN(i + hc, h - 1);
00387         lined = datad + wpl * i;
00388         linemina = dataa + wpla * imin;
00389         linemaxa = dataa + wpla * imax;
00390         for (j = 0; j < w; j++) {
00391             jmin = L_MAX(j - 1 - wc, 0);
00392             jmax = L_MIN(j + wc, w - 1);
00393             val = linemaxa[jmax] - linemaxa[jmin]
00394                   - linemina[jmax] + linemina[jmin];
00395             val = (l_uint8)(norm * val);
00396             SET_DATA_BYTE(lined, j, val);
00397         }
00398     }
00399 
00400         /*------------------------------------------------------------*
00401          *          now fix normalization for boundary pixels         *
00402          *------------------------------------------------------------*/
00403     for (i = 0; i <= hc; i++) {    /* first hc + 1 lines */
00404         hn = hc + i;
00405         normh = (l_float32)fhc / (l_float32)hn;   /* > 1 */
00406         lined = datad + wpl * i;
00407         for (j = 0; j <= wc; j++) {
00408             wn = wc + j;
00409             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00410             val = GET_DATA_BYTE(lined, j);
00411             val = (l_uint8)(val * normh * normw);
00412             SET_DATA_BYTE(lined, j, val);
00413         }
00414         for (j = wc + 1; j < wmwc; j++) {
00415             val = GET_DATA_BYTE(lined, j);
00416             val = (l_uint8)(val * normh);
00417             SET_DATA_BYTE(lined, j, val);
00418         }
00419         for (j = wmwc; j < w; j++) {
00420             wn = wc + w - j;
00421             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00422             val = GET_DATA_BYTE(lined, j);
00423             val = (l_uint8)(val * normh * normw);
00424             SET_DATA_BYTE(lined, j, val);
00425         }
00426     }
00427 
00428     for (i = hmhc; i < h; i++) {  /* last hc lines */
00429         hn = hc + h - i;
00430         normh = (l_float32)fhc / (l_float32)hn;   /* > 1 */
00431         lined = datad + wpl * i;
00432         for (j = 0; j <= wc; j++) {
00433             wn = wc + j;
00434             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00435             val = GET_DATA_BYTE(lined, j);
00436             val = (l_uint8)(val * normh * normw);
00437             SET_DATA_BYTE(lined, j, val);
00438         }
00439         for (j = wc + 1; j < wmwc; j++) {
00440             val = GET_DATA_BYTE(lined, j);
00441             val = (l_uint8)(val * normh);
00442             SET_DATA_BYTE(lined, j, val);
00443         }
00444         for (j = wmwc; j < w; j++) {
00445             wn = wc + w - j;
00446             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00447             val = GET_DATA_BYTE(lined, j);
00448             val = (l_uint8)(val * normh * normw);
00449             SET_DATA_BYTE(lined, j, val);
00450         }
00451     }
00452 
00453     for (i = hc + 1; i < hmhc; i++) {    /* intermediate lines */
00454         lined = datad + wpl * i;
00455         for (j = 0; j <= wc; j++) {   /* first wc + 1 columns */
00456             wn = wc + j;
00457             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00458             val = GET_DATA_BYTE(lined, j);
00459             val = (l_uint8)(val * normw);
00460             SET_DATA_BYTE(lined, j, val);
00461         }
00462         for (j = wmwc; j < w; j++) {   /* last wc columns */
00463             wn = wc + w - j;
00464             normw = (l_float32)fwc / (l_float32)wn;   /* > 1 */
00465             val = GET_DATA_BYTE(lined, j);
00466             val = (l_uint8)(val * normw);
00467             SET_DATA_BYTE(lined, j, val);
00468         }
00469     }
00470 
00471     return;
00472 }
00473 
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Defines