ptrdiff_t BinarizeSavuola2Inplace ( quint8*   pageBits,
                                    ptrdiff_t pageStride,
                                    ptrdiff_t pageWidth,
                                    ptrdiff_t pageHeight,
                                    ptrdiff_t radius,
                                    double    factor
                                  ) throw() {
  assert(pageBits);

  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 swStride  = ( pageWidth + 15 ) & (~15); // align by 16 bytes each row
  ptrdiff_t wmovesize = 2 * radius * swStride;

  quint8  *swindow_unaligned = static_cast<quint8*> ( malloc( swStride + wmovesize + 16 ) );
  if (!swindow_unaligned)
    return 0;

  quint8 *swindow     = (quint8*) ( (ptrdiff_t ( swindow_unaligned ) + 15) & ( ~15 ) );
  quint8 *currWRow    = swindow + radius * swStride;
  quint8 *bottomWRow  = swindow + wmovesize;

  // fill swindow
  quint8* cwRow = swindow;

  // top part - is a copy of first row radius times
  while ( cwRow != currWRow ) {
    memcpy(cwRow, pageBits, pageWidth );
    cwRow += swStride;
  }

  // bottom part - is a copy of input image
  quint8 *srcRow = pageBits;
  while ( cwRow != (bottomWRow + swStride) ) {
    memcpy(cwRow, srcRow, pageWidth );
    cwRow  += swStride;
    srcRow += pageStride;
  }

  ptrdiff_t y,x, ymin,ymax, xmin, xmax;
  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;
    }

  quint8*   lRow = pageBits + pageStride * qMin( pageHeight-1, radius );

  for ( y = 0, srcRow = pageBits; y < pageHeight; ++y, srcRow += pageStride ) {
    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
      srcRow [ x ] = currWRow [ 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 ];
      }
    }

    // Maintain sliding window
    memmove ( swindow,    swindow + swStride, wmovesize );
    memcpy  ( bottomWRow, lRow,               pageWidth );

    // Housekeeping for top and bottom boundary cases
    if ( ymin ) {
      // Subtract top row from row sums
      for ( x = 0; x < pageWidth; ++x ) {
        cpix = swindow [ x ];
        pISum [ x ] -= cpix;
        pISSum[ x ] -= cpix * cpix;
      }
    }

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

  free(swindow_unaligned);

  return 1;
}