ngry

ptrdiff_t BinarizeSavuola2 ( quint8 const* pageBits,
                             ptrdiff_t pageStride,
                             ptrdiff_t pageWidth,
                             ptrdiff_t pageHeight,
                             ptrdiff_t radius,
                             double    factor,
                             QImage&   dst
                            ) throw() {
  assert(pageBits);

  dst = QImage ( pageWidth, pageHeight, QImage::Format_Indexed8 );

  QVector<QRgb> colorTable (256);
  for (size_t i = 0; i < 256; ++i) {
    QColor col( i, i, i);
    colorTable[i] = col.rgb();
  }
  dst.setColorTable ( colorTable );

  QVarLengthArray <double, 2048> rowSum ( 2 * pageWidth );

  double *pISum  = rowSum.data(),     // sums
         *pISSum = pISum + pageWidth; // sums of squares

  memset( rowSum.data(), 0, 2 * pageWidth * sizeof(double) );

  ptrdiff_t y,x, ymin,ymax, xmin, xmax, dstStride = dst.bytesPerLine();
  quint8 const *srcRow;
  quint8 *dstRow;
  size_t    cpix;

  double mean, std, area, wsum, wssum, factor1 = 1. - factor;

  // initialize rowSum array

  for ( y = 0, srcRow = pageBits; y < radius; ++y, srcRow += pageStride )
    for ( x = 0; x < pageWidth; ++x ) {
      cpix = srcRow [ x ];
      pISum [ x ] += cpix;
      pISSum[ x ] += cpix * cpix;
    }

  for ( y = 0, srcRow = pageBits, dstRow = dst.bits(); y < pageHeight; ++y, srcRow += pageStride, dstRow += dstStride ) {

    ymin = qMax(ptrdiff_t(0),y-radius);
    ymax = qMin(pageHeight-1,y+radius);

    // initialize window sums
    for ( wsum = wssum = 0., x = 0; x < radius; ++x ) {
      wsum  += pISum [ x ];
      wssum += pISSum[ x ];
    }

    // loop on row
    for ( x = 0; x < pageWidth; ++x ) {

      xmin = qMax(ptrdiff_t(0),x-radius);
      xmax = qMin(pageWidth-1, x+radius);

      area = (xmax-xmin+1.)*(ymax-ymin+1.);

      mean = wsum / area;
      std  = sqrt((wssum - wsum*wsum/area)/area);

      // threshold the current pixel

      dstRow [ x ] = srcRow [ x ] < mean * ( factor1 + factor *  std * .0078125 )  ? 0 : 255;

      // Housekeeping for left and right boundary cases

      if ( xmin ) {
        // Subtract leftmost column from window sums
        wsum  -= pISum [ xmin ];
        wssum -= pISSum[ xmin ];
      }

      if ( xmax < pageWidth - 1 ) {
        // Add rightmost collumn to window histogram
        wsum  += pISum [ xmax ];
        wssum += pISSum[ xmax ];
      }
    }

    // Housekeeping for top and bottom boundary cases

    if ( ymin ) {
      // Subtract top row from row sums
      quint8 const* row = pageBits + ymin * pageStride;
      for ( x = 0; x < pageWidth; ++x ) {
        cpix = row [ x ];
        pISum [ x ] -= cpix;
        pISSum[ x ] -= cpix*cpix;
      }
    }

    if ( ymax < pageHeight - 1 ) {
      // Add bottom collumn to row sums
      quint8 const* row = pageBits + ymax * pageStride;
      for ( x = 0; x < pageWidth; ++x ) {
        cpix = row [ x ];
        pISum [ x ] += cpix;
        pISSum[ x ] += cpix*cpix;
      }
    }

  }

  return 1;
}