gevd_float_operators.cxx

Go to the documentation of this file.

// This is gel/gevd/gevd_float_operators.cxx
#include "gevd_float_operators.h"
//:
// \file

#include <vcl_iostream.h>
#include <vcl_fstream.h>
#include <vcl_algorithm.h>
#include <vcl_cassert.h>
#include <vnl/vnl_math.h> // for pi_over_2
#include <vnl/vnl_float_3.h>
#include <vnl/vnl_cross.h>

#include <vcl_cmath.h>
#include <vcl_cstring.h>

#include "gevd_pixel.h"
#include "gevd_xpixel.h"
#include "gevd_bufferxy.h"
#ifdef DEBUG
# include <vul/vul_timer.h>
#endif

#if defined(VCL_VC) || defined(VCL_BORLAND)
inline static double rint(double v)
{
  return  v - vcl_floor(v) < 0.5  ?  vcl_floor(v)  :  vcl_ceil(v);
}
#endif

const unsigned char DIR0 = 8, DIR1 = 9, DIR2 = 10, DIR3 = 11, DIR4 = 12;
// const int DIS[] = { 1, 1, 0,-1,-1,-1, 0, 1, // 8-connected neighbors
//                     1, 1, 0,-1,-1,-1, 0, 1, // wrapped by 2PI
//                     1, 1, 0,-1,-1,-1, 0, 1};
// const int DJS[] = { 0, 1, 1, 1, 0,-1,-1,-1,
//                     0, 1, 1, 1, 0,-1,-1,-1,
//                     0, 1, 1, 1, 0,-1,-1,-1};

//: Convolves \a from image with 2D filter and stores values in \a to image.
// O(m*n*k).
// The filter \a kernel has odd width and height,
// and is either even or odd along x- or y-axis.
// Assume image data is in row-major order.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.

float
gevd_float_operators::Convolve(const gevd_bufferxy& from,
                               const gevd_bufferxy& kernel,
                               gevd_bufferxy*& to)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Convolve image";
#endif
  bool overwrite = to == &from;
  gevd_bufferxy*& t = to;
  if (overwrite) to=0;
  to = gevd_float_operators::Allocate(to, from);
  const int wx = kernel.GetSizeX(), wy = kernel.GetSizeY();
  const int rx = wx/2, ry = wy/2;
  const int xhi = from.GetSizeX() - wx + 1;
  const int yhi = from.GetSizeY() - wy + 1;
  for (int y = 0; y < yhi; y++)
    for (int x = 0; x < xhi; x++) {
      float conv = 0;           // more efficient if use even/odd symmetry
      for (int j = 0; j < wy; j++)
        for (int i = 0; i < wx; i++)
          conv += floatPixel(from, x+i, y+j) * floatPixel(kernel, i, j);
      floatPixel(*to, x+rx, y+ry) = conv;
    }
  FillFrameX(*to, 0, rx);       // pad border with 0
  FillFrameY(*to, 0, ry);
#ifdef DEBUG
  vcl_cout << " in " << t.real() << " msecs.\n";
#endif
  if (overwrite) { gevd_float_operators::Update(*t,*to); delete to; to = t; }
  return 1;                     // extra scaling factor
}

//: Correlate \a from image with 2D filter and stores values in \a to image.
// O(m*n*k).
// The filter \a kernel has odd width and height,
// and is either even or odd along x- or y-axis.
// Assume image data is in row-major order.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.

float
gevd_float_operators::Correlation(const gevd_bufferxy& from,
                                  const gevd_bufferxy& kernel,
                                  gevd_bufferxy*& to)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Correlate image";
#endif
  bool overwrite = to == &from;
  gevd_bufferxy*& t = to;
  if (overwrite) to=0;
  to = gevd_float_operators::Allocate(to, from);
  const int wx = kernel.GetSizeX(), wy = kernel.GetSizeY();
  const int rx = wx/2, ry = wy/2;
  const int xhi = from.GetSizeX() - wx + 1;
  const int yhi = from.GetSizeY() - wy + 1;
  const int sum1 = wx*wy;
  for (int y = 0; y < yhi; y++)
    for (int x = 0; x < xhi; x++) {
      register double sumx=0, sumy=0, sumxx=0, sumyy=0, sumxy=0, xval, yval;
      for (int j = 0; j < wy; j++)
        for (int i = 0; i < wx; i++) {
          xval = floatPixel(from, x+i, y+j);
          yval = floatPixel(kernel, i, j);
        sumxy += xval * yval;           // accumulate correlation value
        sumx += xval;
        sumy += yval;
        sumxx += xval * xval;
        sumyy += yval * yval;
      }
      double varx = sum1 * sumxx - sumx * sumx; // all multiplied with sum1
      double vary = sum1 * sumyy - sumy * sumy;
      double cvar = sum1 * sumxy - sumx * sumy; // linear correlation coeft
      if (varx!=0 && vary!=0) cvar /= vcl_sqrt(varx * vary);
      floatPixel(*to, x+rx, y+ry) = (float)cvar;
    }
  FillFrameX(*to, 0, rx);       // pad border with 0
  FillFrameY(*to, 0, ry);
#ifdef DEBUG
  vcl_cout << " in " << t.real() << " msecs.\n";
#endif
  if (overwrite) { gevd_float_operators::Update(*t,*to); delete to; to = t; }
  return 1;                     // extra scaling factor
}


//: Correlate \a from image with 2D filter and stores values in \a to image.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.
float
gevd_float_operators::CorrelationAlongAxis(const gevd_bufferxy& from,
                                           const gevd_bufferxy& kernel,
                                           gevd_bufferxy*& to)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Correlate image";
#endif
  bool overwrite = to == &from;
  gevd_bufferxy*& t = to;
  if (overwrite) to=0;
  to = gevd_float_operators::Allocate(to, from);
  to->Clear();
  const int wx = kernel.GetSizeX(), wy = kernel.GetSizeY();
  const int rx = wx/2, ry = wy/2;
  const int xhi = from.GetSizeX() - wx + 1;
  const int yhi = from.GetSizeY() - wy + 1;
  const int sum1 = wx*wy;
  for (int x = 0, y = yhi/2; x < xhi; ++x) {
    register double sumx=0, sumy=0, sumxx=0, sumyy=0, sumxy=0, xval, yval;
    for (int j = 0; j < wy; j++)
      for (int i = 0; i < wx; i++) {
        xval = floatPixel(from, x+i, y+j);
        yval = floatPixel(kernel, i, j);
        sumxy += xval * yval;           // accumulate correlation value
        sumx += xval;
        sumy += yval;
        sumxx += xval * xval;
        sumyy += yval * yval;
      }
    double varx = sum1 * sumxx - sumx * sumx; // all multiplied with sum1
    double vary = sum1 * sumyy - sumy * sumy;
    double cvar = sum1 * sumxy - sumx * sumy;   // linear correlation coeft
    if (varx!=0 && vary!=0) cvar /= vcl_sqrt(varx * vary);
    floatPixel(*to, x+rx, y+ry) = (float)cvar;
  }
  for (int x = xhi/2, y = 0; y < yhi; ++y) {
    register double sumx=0, sumy=0, sumxx=0, sumyy=0, sumxy=0, xval, yval;
    for (int j = 0; j < wy; j++)
      for (int i = 0; i < wx; i++) {
        xval = floatPixel(from, x+i, y+j);
        yval = floatPixel(kernel, i, j);
        sumxy += xval * yval;           // accumulate correlation value
        sumx += xval;
        sumy += yval;
        sumxx += xval * xval;
        sumyy += yval * yval;
      }
    double varx = sum1 * sumxx - sumx * sumx; // all multiplied with sum1
    double vary = sum1 * sumyy - sumy * sumy;
    double cvar = sum1 * sumxy - sumx * sumy;   // linear correlation coeft
    if (varx!=0 && vary!=0) cvar /= vcl_sqrt(varx * vary);
    floatPixel(*to, x+rx, y+ry) = (float)cvar;
  }
#ifdef DEBUG
  vcl_cout << " in " << t.real() << " msecs.\n";
#endif
  if (overwrite) { gevd_float_operators::Update(*t,*to); delete to; to = t; }
  return 1;                     // extra scaling factor
}


//: Read 2D kernel from a file: width height k_x_y .....

gevd_bufferxy*
gevd_float_operators::Read2dKernel(const char* filename)
{
  vcl_ifstream infile(filename, vcl_ios_in); // open the file
  if (!infile)
    return NULL;
  int width, height;
  infile >> width;
  infile >> height;
  if (width < 1 || height < 1) return NULL;
  gevd_bufferxy* kernel = new gevd_bufferxy(width, height, bits_per_float);
  for (int y = 0; y < height; y++)
    for (int x = 0; x < width; x++)
      infile >> floatPixel(*kernel, x, y);
  return kernel;
}

//: Convolves \a from image with two separable filters, along directions x and y, and stores values in \a to image.
// O(m*n*k).
// The filter \a kernel has length 2*radius + 1, and is either even or odd.
// Assume image data is in row-major order.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.

float
gevd_float_operators::Convolve(const gevd_bufferxy& from, gevd_bufferxy*& to,
                               const float* xkernel, const int xradius,
                               const bool xevenp,
                               const float* ykernel, const int yradius,
                               const bool yevenp,
                               const bool xwrap, const bool ywrap)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Convolve image";
#endif
  bool overwrite = to == &from;
  gevd_bufferxy*& t = to;
  if (overwrite) to=0;
  to = gevd_float_operators::Allocate(to, from);
  const int sizeX = to->GetSizeX(), sizeY = to->GetSizeY();
  const int ylo = yradius, yhi = sizeY - yradius;
  const int kborder = 2*yradius;

  // 1. Setup the pipeline of 4*yradius+1 lines, convolved along x-axis
  float** cache = new float*[4*yradius+1];
  float** pipeline = cache+yradius;
  float* row;
  for (int p = 0; p <= kborder; ++p) {
    pipeline[p] = row = new float[sizeX];
    gevd_float_operators::Convolve(&floatPixel(from, 0, p), // row-major order
                                   row, sizeX,
                                   xkernel, xradius, xevenp, xwrap);
  }
  if (ywrap)                    // circular wrap at 2 end rows
    for (int r = 1; r <= yradius; r++) {
      pipeline[-r] = row = new float[sizeX];
      gevd_float_operators::Convolve(&floatPixel(from, 0, sizeY-r), // row-major order
                                     row, sizeX,
                                     xkernel, xradius, xevenp, xwrap);
      pipeline[kborder+r] = row = new float[sizeX];
      for (int i = 0; i < sizeX; i++) // copy previous convolved result
        row[i] = pipeline[r-1][i];
    }
  else                          // reflection at 2 end rows
    for (int r = 1; r <= yradius; r++) {
      pipeline[-r] = row = new float[sizeX];
      for (int i = 0; i < sizeX; i++) // copy previous convolved result
        row[i] = pipeline[r][i];
      pipeline[kborder+r] = row = new float[sizeX];
      gevd_float_operators::Convolve(&floatPixel(from, 0, sizeY-1-r), // row-major order
                                     row, sizeX,
                                     xkernel, xradius, xevenp, xwrap);
    }

  // 2. Convolve along y-axis, shifting pipeline by 1 each time.
  if (yevenp)
  {
    int y=0;
    for (int yy = 0; y < ylo; ++y, ++yy) { // reflect/wrap at ymin
      row = &floatPixel(*to, 0, y); // row-major order
      for (int x = 0; x < sizeX; x++) {
        float sum = ykernel[yradius] * pipeline[yy][x];
        for (int k = 1; k <= yradius; k++)
          sum += ykernel[yradius-k] * (pipeline[yy-k][x] + pipeline[yy+k][x]);
        row[x] = sum;
      }
    }
    int p = kborder+1;
    for ( ; y < yhi; ++y) {     // convolution along y-axis
      row = &floatPixel(*to, 0, y);
      for (int x = 0; x < sizeX; x++) {
        float sum = ykernel[yradius] * pipeline[yradius][x];
        for (int k = 1; k <= yradius; k++)
          sum += ykernel[yradius-k] * (pipeline[yradius-k][x] + // even kernel
                                       pipeline[yradius+k][x]);
        row[x] = sum;
      }
      if (p < sizeY) {
        row = pipeline[0];      // next line
        for (int k = 0; k < kborder; k++) // shift the lines of
          pipeline[k] = pipeline[k+1]; // the pipeline by 1
        gevd_float_operators::Convolve(&floatPixel(from, 0, p++), // row-major order
                                       row, sizeX,
                                       xkernel, xradius, xevenp, xwrap); // new line x-conv
        pipeline[kborder] = row; // update pipeline
      }
    }
    for (int yy = yradius+1; y < sizeY; y++, yy++) {  // reflect/wrap at ymax
      row = &floatPixel(*to, 0, y); // row-major order
      for (int x = 0; x < sizeX; x++) {
        float sum = ykernel[yradius] * pipeline[yy][x];
        for (int k = 1; k <= yradius; k++)  // convolution along y-axis
          sum += ykernel[yradius-k] * (pipeline[yy-k][x] + pipeline[yy+k][x]);
        row[x] = sum;
      }
    }
  }
  else
  {
    int y=0;
    for (int yy = 0; y < ylo; y++, yy++) { // reflect at ymin
      row = &floatPixel(*to, 0, y); // row-major order
      for (int x = 0; x < sizeX; x++) {
        float sum = ykernel[yradius] * pipeline[yy][x];
        for (int k = 1; k <= yradius; k++)
          sum += ykernel[yradius-k] * (pipeline[yy-k][x] - pipeline[yy+k][x]);
        row[x] = sum;
      }
    }
    int p = kborder+1;
    for ( ; y < yhi; y++) {
      row = &floatPixel(*to, 0, y);
      for (int x = 0; x < sizeX; x++) {
        float sum = ykernel[yradius] * pipeline[yradius][x];
        for (int k = 1; k <= yradius; k++)  // convolution along y-axis
          sum += ykernel[yradius-k] * (pipeline[yradius-k][x] - // odd kernel
                                       pipeline[yradius+k][x]);
        row[x] = sum;
      }
      if (p < sizeY) {
        row = pipeline[0];      // next line
        for (int k = 0; k < kborder; k++)   // shift the lines of
          pipeline[k] = pipeline[k+1];  // the pipeline by 1
        gevd_float_operators::Convolve(&floatPixel(from, 0, p++), // row-major order
                                       row, sizeX,
                                       xkernel, xradius, xevenp, xwrap); // new line x-conv
        pipeline[kborder] = row; // update pipeline
      }
    }
    for (int yy = yradius+1; y < sizeY; y++, yy++) { // reflect/wrap at ymax
      row = &floatPixel(*to, 0, y); // row-major order
      for (int x = 0; x < sizeX; x++) {
        float sum = ykernel[yradius] * pipeline[yy][x];
        for (int k = 1; k <= yradius; k++)  // convolution along y-axis
          sum += ykernel[yradius-k] * (pipeline[yy-k][x] - pipeline[yy+k][x]);
        row[x] = sum;
      }
    }
  }
  for (int p = 0; p <= 4*yradius; p++)
    delete[] cache[p];         // Free lines in pipeline cache
  delete[] cache;
#ifdef DEBUG
  vcl_cout << " in " << t.real() << " msecs.\n";
#endif
  if (overwrite) { gevd_float_operators::Update(*t,*to); delete to; to = t; }
  return 1;                     // assume normalized kernel
}


//: Convolves \a from image with a filter along x, and a running sum along y-axis, then stores values in \a to image.
// O(m*n*k).

float
gevd_float_operators::Convolve(const gevd_bufferxy& from, gevd_bufferxy*& to,
                               const float* xkernel, const int xradius,
                               const bool xevenp,
                               const int yradius,
                               const bool xwrap, const bool ywrap)
{
  float* ykernel = new float[2*yradius+1];
  const bool yevenp = true;
  for (int i = 0; i < 2*yradius+1; i++)
    ykernel[i] = 1;
  float fact = gevd_float_operators::Convolve(from, to,
                                              xkernel, xradius, xevenp,
                                              ykernel, yradius, yevenp,
                                              xwrap, ywrap);
  delete[] ykernel;
  return fact;
}

#if 0 // commented out
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Convolve image";
#endif
  to = gevd_float_operators::Allocate(to, from);
  const int sizeX = to->GetSizeX(), sizeY = to->GetSizeY();
  const int ylo = yradius, yhi = sizeY - yradius;
  const int kborder = 2*yradius;

  // 1. Setup the pipeline of 4*yradius+1 lines, convolved along x-axis
  float** cache = new float*[4*yradius+1];
  float** pipeline = cache+yradius;
  double* rsum = new double[sizeX]; // avoid addition/subtraction errors
  float* row = NULL;            // current row
  for (int p = 0; p <= kborder; p++) {
    pipeline[p] = row = new float[sizeX];
    gevd_float_operators::Convolve(&floatPixel(from, 0, p), // row-major order
                                   row, sizeX,
                                   xkernel, xradius, xevenp, xwrap);
  }
  // **** wrapping of the from buffer at the end ****
  if (ywrap)                    // circular wrap at 2 end rows
    for (int r = 1; r <= yradius; r++) {
      pipeline[-r] = row = new float[sizeX];
      gevd_float_operators::Convolve(&floatPixel(from, 0, sizeY-r), // row-major order
                                     row, sizeX,
                                     xkernel, xradius, xevenp, xwrap);
      pipeline[kborder+r] = row = new float[sizeX];
      for (int i = 0; i < sizeX; i++) // copy previous convolved result
        row[i] = pipeline[r-1][i];
    }
  else                          // reflection at 2 end rows
    for (int r = 1; r <= yradius; r++) {
      pipeline[-r] = row = new float[sizeX];
      for (int i = 0; i < sizeX; i++) // copy previous convolved result
        row[i] = pipeline[r][i];
      pipeline[kborder+r] = row = new float[sizeX];
      gevd_float_operators::Convolve(&floatPixel(from, 0, sizeY-1-r), // row-major order
                                     row, sizeX,
                                     xkernel, xradius, xevenp, xwrap);
    }

  // 2. Running sum along y-axis, shifting pipeline by 1 each time.
  int y = 0, yy = 0;
  row = &floatPixel(*to, 0, y); // row-major order
  for (int x = 0; x < sizeX; x++) { // running sum at first row
    double sum = pipeline[yy][x];
    for (int k = 1; k <= yradius; k++)
      sum += (pipeline[yy-k][x] + pipeline[yy+k][x]);
    row[x] = float(rsum[x] = sum);
  }
  for (y++, yy++; y < ylo; y++, yy++) { // reflect/wrap at ymin
    row = &floatPixel(*to, 0, y); // row-major order
    for (int x = 0; x < sizeX; x++)
      row[x] = float(rsum[x] = rsum[x] -
                     pipeline[yy-yradius-1][x] + pipeline[yy+yradius][x]);
  }
  for ( ; y < yhi; y++) {       // convolution along y-axis
    row = &floatPixel(*to, 0, y);
    for (int x = 0; x < sizeX; x++)
      row[x] = float(rsum[x] = rsum[x] -
                     pipeline[-1][x] + pipeline[kborder][x]);
    if (p < sizeY) {
      row = pipeline[0];        // next line
      for (int k = 0; k < kborder; k++) // shift the lines of
        pipeline[k] = pipeline[k+1]; // the pipeline by 1
      gevd_float_operators::Convolve(&floatPixel(from, 0, p++), // row-major order
                                     row, sizeX,
                                     xkernel, xradius, xevenp, xwrap); // new x-convolution
      pipeline[kborder] = row; // update pipeline
    }
  }
  for (yy = yradius+1; y < sizeY; y++, yy++) {  // reflect/wrap at ymax
    row = &floatPixel(*to, 0, y); // row-major order
    for (int x = 0; x < sizeX; x++)
      row[x] = float(rsum[x] = rsum[x] -
                     pipeline[yy-yradius-1][x] + pipeline[yy+yradius][x]);
  }
  for (int p = 0; p <= 4*yradius; p++)
    delete[] cache[p];         // Free lines in pipeline cache
  delete[] cache;
  delete[] rsum;
#ifdef DEBUG
  vcl_cout << " in " << t.real() << " msecs.\n";
#endif
  return 2*yradius+1;           // return multiplication factor
}
#endif

//: Convolve a linear array of data, wrapping at the 2 ends.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.

float
gevd_float_operators::Convolve(const float* from, float*& to, const int len,
                               const float* kernel, const int kradius,
                               const bool evenp, const bool wrap)
{
  bool overwrite = to == from;
  float*& t = to;
  if (!to || overwrite) to = new float[len];
  const int xlo = kradius, xhi = len - kradius;
  const int kborder = 2*kradius;
  float *cache, *pipeline;
  int p = gevd_float_operators::SetupPipeline(from, len, kradius, wrap,
                                              cache, pipeline);

  // Convolution along x with above pipeline
  float sum;
  if (evenp)
  {
    int x=0;
    for (int xx = 0; x < xlo; ++x, ++xx) {
      sum = kernel[kradius] * pipeline[xx];
      for (int k = 1; k <= kradius; k++)
        sum += kernel[kradius-k] * (pipeline[xx-k] + pipeline[xx+k]);
      to[x] = sum;
    }
    for ( ; x < xhi; x++) {
      sum = kernel[kradius] * pipeline[kradius];
      for (int k = 1; k <= kradius; k++)
        sum += kernel[kradius-k] * (pipeline[kradius-k] + // plus for
                                    pipeline[kradius+k]); // symmetric
      to[x] = sum;
      if (p < len) {
        for (int k = 0; k < kborder; k++)   // shift the floats of
          pipeline[k] = pipeline[k+1];  // the pipeline by 1
        pipeline[kborder] = from[p++];  // update pipeline
      }
    }
    for (int xx = kradius+1; x < len; x++, xx++) {
      sum = kernel[kradius] * pipeline[xx];
      for (int k = 1; k <= kradius; k++)
        sum += kernel[kradius-k] * (pipeline[xx-k] + pipeline[xx+k]);
      to[x] = sum;
    }
  }
  else
  {
    int x=0;
    for (int xx = 0; x < xlo; ++x, ++xx) {
      sum = kernel[kradius] * pipeline[xx];
      for (int k = 1; k <= kradius; k++)
        sum += kernel[kradius-k] * (pipeline[xx-k] - pipeline[xx+k]);
      to[x] = sum;
    }
    for ( ; x < xhi; x++) {
      sum = kernel[kradius] * from[x];
      for (int k = 1; k <= kradius; k++)
        sum += kernel[kradius-k] * (pipeline[kradius-k] - // minus for
                                    pipeline[kradius+k]); // antisymmetric
      to[x] = sum;
      if (p < len) {
        for (int k = 0; k < kborder; k++)   // shift the floats of
          pipeline[k] = pipeline[k+1];  // the pipeline by 1
        pipeline[kborder] = from[p++];  // update pipeline
      }
    }
    for (int xx = kradius+1; x < len; x++, xx++) {
      sum = kernel[kradius] * pipeline[xx];
      for (int k = 1; k <= kradius; k++)
        sum += kernel[kradius-k] * (pipeline[xx-k] - pipeline[xx+k]);
      to[x] = sum;
    }
  }
  delete[] cache;
  if (overwrite) {
    for (int x=0 ; x < len; ++x) t[x] = to[x];
    delete[] to;
    to = t;
  }
  return 1;                     // assume normalized kernel
}

//: For large smoothing sigma, use running sum to avoid floating multiplications.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.
float
gevd_float_operators::RunningSum(float* from, float*& to, const int len,
                                 const int kradius, const bool wrap)
{
  bool overwrite = to == from;
  if (!to || overwrite) to = new float[len];
  const int xlo = kradius, xhi = len - kradius;
  const int kborder = 2*kradius;
  float *cache, *pipeline;
  int p = gevd_float_operators::SetupPipeline(from, len, kradius, wrap,
                                              cache, pipeline);
  // Running sum along x with above pipeline
  int x = 0, xx = 0;
  float sum = pipeline[x];      // pre-compute running sum
  for (int k = 1; k <= kradius; k++)
    sum += (pipeline[xx-k] + pipeline[xx+k]);
  to[x] = sum;
  for (x++, xx++; x < xlo; x++, xx++) {
    sum += pipeline[xx+kradius] - pipeline[xx-kradius-1];
    to[x] = sum;
  }
  for ( ; x < xhi; x++) {
    sum += pipeline[kborder] - pipeline[-1];
    to[x] = sum;
    if (p < len) {
      for (int k = -1; k < kborder; k++)    // shift the floats of
        pipeline[k] = pipeline[k+1];    // the pipeline by 1
      pipeline[kborder] = from[p++];    // update pipeline
    }
  }
  for ( ; x < len; x++, xx++) {
    sum += pipeline[xx+kradius] - pipeline[xx-kradius-1];
    to[x] = sum;
  }
  delete[] cache;
  if (overwrite) {
    for (int x=0 ; x < len; ++x) from[x] = to[x];
    delete[] to;
    to = from;
  }
  return float(2*kradius+1);            // magnification factor
}

//: Read 1d odd/even kernel from a file: width k_i .....

bool
gevd_float_operators::Read1dKernel(const char* filename,
                                   float*& kernel, int& radius, bool& evenp)
{
  vcl_ifstream infile(filename, vcl_ios_in); // open the file
  if (!infile)
    return false;
  int width;
  infile >> width;
  if (width < 1) return false;
  radius = (width - 1)/2;
  width = 2*radius + 1;
  delete[] kernel;
  kernel = new float[width];
  for (int i = 0; i < width; i++)
    infile >> kernel[i];
  evenp = true;
  for (int i = 1; i <= radius; i++)
    if (kernel[radius-i] != kernel[radius+i])
      evenp = false;
  if (evenp)                    // double check that kernel is even
    return true;
  for (int i = 1; i <= radius; i++)
    if (kernel[radius-i] != -kernel[radius+i])
      evenp = true;
  if (!evenp)                   // double check that kernel is odd
    return true;
  delete[] kernel; kernel = NULL;
  return false;                 // invalid kernel
}


//: Convolves the \a from image with a 2D Gaussian filter with standard deviation \a sigma.
// O(m*n*k).
// The 2D convolution is decomposed into 2 1D convolutions:
// Gxy * I = Gy * (Gx * I).
// The frame with width equal to the radius of the Gaussian kernel is
// filled with zero.

float
gevd_float_operators::Gaussian(gevd_bufferxy& from, gevd_bufferxy*& to, const float sigma,
                               const bool xwrap, const bool ywrap)
{
  float* kernel = NULL;
  int radius = 0;
  if (!gevd_float_operators::Find1dGaussianKernel(sigma, kernel, radius)) {
    to = gevd_float_operators::Allocate(to, from);
    if (to != &from)
      gevd_float_operators::Update(*to, from); // just a copy, no smoothing needed
  } else {
    bool evenp = true;
    // 2 1D convolutions O(k), instead of a 2D convolution O(k^2).
    gevd_float_operators::Convolve(from, to,
                                   kernel, radius, evenp,
                                   kernel, radius, evenp,
                                   xwrap, ywrap);
  }
  delete[] kernel;
  return 1;                     // return multiplication factor
}


//:
// Returns the kernel array for Gaussian with given standard deviation
// \a sigma [pixel], and cutoff at min/max = \a fuzz.
// The area under the discrete Gaussian is normalized to 1.

bool
gevd_float_operators::Find1dGaussianKernel(const float sigma,
                                           float*& kernel, int& radius,
                                           const float fuzz)
{
  for (radius = 0; Gaussian(float(radius), sigma) > fuzz; radius++)
    {;}                                         // find radius
  if (radius == 0)
    return false;

  kernel = new float[2*radius + 1];             // create kernel
  for (int i=0; i<=radius; ++i)
    kernel[radius+i] = kernel[radius-i] = Gaussian(float(i), sigma);
  float sum = 0;
  for (int i= 0; i <= 2*radius; ++i)
    sum += kernel[i];                           // find integral of weights
  for (int i= 0; i <= 2*radius; ++i)
    kernel[i] /= sum;                           // normalize by integral
#ifdef DEBUG
  vcl_cout << "Gaussian kernel = ";
  for (int i= 0; i <= 2*radius; ++i)
    vcl_cout << kernel[i] << ' ';
  vcl_cout << vcl_endl;
#endif
  return true;
}


float
gevd_float_operators::Gaussian(const float x, const float sigma)
{
  double x_on_sigma = x / sigma;
  return (float)vcl_exp(- x_on_sigma * x_on_sigma / 2);
}


//:
// Setup the cache for reflection/wrapping at the borders,
// and the center pipeline. Only cache should be deleted after
// you're done, not the pipeline, since pipeline shares the same space.

int
gevd_float_operators::SetupPipeline(const float* data, const int len,
                                    const int kradius, const bool wrap,
                                    float*& cache, float*& pipeline)
{
  cache = new float[4*kradius+1]; // Setup the cache of 4*kradius+1 floats
  pipeline = cache+kradius;     // center pipeline is 2*kradius+1 floats
  int p = 0;
  for (; p <= 2*kradius; p++) pipeline[p] = data[p];
  if (wrap)                     // circular wrap at 2 end points
    for (int r = 1; r <= kradius; r++) {
      pipeline[-r] = data[len-r];
      pipeline[2*kradius+r] = data[r-1];
    }
  else                          // reflection at 2 end points
    for (int r = 1; r <= kradius; r++) {
      pipeline[-r] = data[r];
      pipeline[2*kradius+r] = data[len-1-r];
    }
  return p;                     // current index in pipeline
}

//: Create a buffer that wraps rows/columns.

gevd_bufferxy*
gevd_float_operators::WrapAlongX(const gevd_bufferxy& img)
{
  const int lx = img.GetSizeX()-1, sy = img.GetSizeY();
  gevd_bufferxy* wrap = gevd_float_operators::SimilarBuffer(img, bits_per_float,  4, sy);
  for (int i = 0; i < 2; i++)
    for (int j = 0; j < sy; j++) {
      floatPixel(*wrap, 2+i, j) = floatPixel(img, i, j);
      floatPixel(*wrap, 1-i, j) = floatPixel(img, lx-i, j);
    }
  return wrap;
}


gevd_bufferxy*
gevd_float_operators::WrapAlongY(const gevd_bufferxy& img)
{
  const int sx = img.GetSizeX(), ly = img.GetSizeY()-1;
  gevd_bufferxy* wrap = gevd_float_operators::SimilarBuffer(img, bits_per_float, sx, 4);
  for (int j = 0; j < 2; j++)
    for (int i = 0; i < sx; i++) {
      floatPixel(*wrap, i, 2+j) = floatPixel(img, i, j);
      floatPixel(*wrap, i, 1-j) = floatPixel(img, i, ly-j);
    }
  return wrap;
}

// Local gradient magnitude and direction in 3x3 window

inline void
LocalGradient(const gevd_bufferxy& smooth, const int i, const int j,
              float& mag, float& gx, float& gy)
{
  gx = floatPixel(smooth, i+1, j) - floatPixel(smooth, i-1, j);
  gy = floatPixel(smooth, i, j+1) - floatPixel(smooth, i, j-1);
  mag = vcl_sqrt(gx*gx + gy*gy);
}

//: Compute the gradient of the intensity surface.
// O(m*n).
// The intensity surface is assumed smoothed by a Gaussian filter,
// or convolved with a low-pass filter.
// Return at each pixel, the magnitude, the vector (dG * I),
// and the multiplication factor to normalize the magnitude.
// The product of detection |dG * I| and localization |dddG * I| of
// the step relative to noise is invariant to filter-size for the
// ideal step edge, and depends only on the shape of the step edge.
// As the filter size (sigma) increases by k, detection or
// signal-to-noise ratio increases by sqrt(k), but localization
// decreases by 1/sqrt(k). This invariance is consistent with the
// uncertainty principle.

float
gevd_float_operators::Gradient(const gevd_bufferxy& smooth,
                               gevd_bufferxy*& magnitude,
                               gevd_bufferxy*& gradx, gevd_bufferxy*& grady,
                               const bool xwrap, const bool ywrap)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Compute local gradient magnitude/direction";
#endif
  magnitude = gevd_float_operators::Allocate(magnitude, smooth);
  gradx = gevd_float_operators::Allocate(gradx, smooth);
  grady = gevd_float_operators::Allocate(grady, smooth);

  // 1. Inside image frame, compute values based on 3x3 window.
  const int frame = 1;
  const int highx = smooth.GetSizeX() - frame; // exclusive bounds
  const int highy = smooth.GetSizeY() - frame;
  for (int j = frame; j < highy; ++j)
    for (int i = frame; i < highx; ++i)
      LocalGradient(smooth, i, j,
                    floatPixel(*magnitude, i, j),
                    floatPixel(*gradx, i, j),
                    floatPixel(*grady, i, j));

  // 2. Along horizontal/vertical borders
  if (xwrap) {                  //  each row wrap at border
    const int lo = 2, hi = 1;
    gevd_bufferxy* pad = gevd_float_operators::WrapAlongX(smooth);
    for (int j = 1; j < highy; ++j) {
      LocalGradient(*pad, lo, j,
                    floatPixel(*magnitude, 0, j),
                    floatPixel(*gradx, 0, j),
                    floatPixel(*grady, 0, j));
      LocalGradient(*pad, hi, j,
                    floatPixel(*magnitude, highx, j),
                    floatPixel(*gradx, highx, j),
                    floatPixel(*grady, highx, j));
    }
    delete pad;
  } else {                      // zero by default
    gevd_float_operators::FillFrameX(*magnitude, 0, frame);
    gevd_float_operators::FillFrameX(*gradx, 0, frame);
    gevd_float_operators::FillFrameX(*grady, 0, frame);
  }
  if (ywrap) {                  // each column wrap at border
    const int lo = 2, hi = 1;
    gevd_bufferxy* pad = gevd_float_operators::WrapAlongY(smooth);
    for (int i = 1; i < highx; ++i) {
      LocalGradient(*pad, i, lo,
                    floatPixel(*magnitude, i, 0),
                    floatPixel(*gradx, i, 0),
                    floatPixel(*grady, i, 0));
      LocalGradient(*pad, i, hi,
                    floatPixel(*magnitude, i, highy),
                    floatPixel(*gradx, i, highy),
                    floatPixel(*grady, i, highy));
    }
    delete pad;
  } else {                      // zero by default
    gevd_float_operators::FillFrameY(*magnitude, 0, frame);
    gevd_float_operators::FillFrameY(*gradx, 0, frame);
    gevd_float_operators::FillFrameY(*grady, 0, frame);
  }
#ifdef DEBUG
  vcl_cout << ", in " << t.real() << " msecs.\n";
#endif
  return 2;                     // return multiplication factor
}

//: Compute slope or first-difference, in linear/circular array.
//
// Special care is taken when \a to and \a from are identical - Peter Vanroose.

float
gevd_float_operators::Slope(float* from, float*& to, const int len,
                            const bool wrap)
{
  bool overwrite = to == from;
  if (!to || overwrite) to = new float[len];
  const int xlo = 1, xhi = len - 1;
  float *cache = NULL, *pipeline;
  int p = gevd_float_operators::SetupPipeline(from, len, 1, wrap, // current pipeline and index
                                              cache, pipeline);
  to[0] = pipeline[+1] - pipeline[-1];
  for (int x = xlo; x < xhi; x++) {
    to[x] = pipeline[2] - pipeline[0];
    if (p < len) {
      for (int k = 0; k < 2; k++) // shift the floats of
        pipeline[k] = pipeline[k+1]; // the pipeline by 1
      pipeline[2] = from[p++];  // update pipeline
    }
  }
  to[xhi] = pipeline[3] - pipeline[1];
  delete[] cache;
  if (overwrite) {
    for (int x=0 ; x < len; ++x) from[x] = to[x];
    delete[] to;
    to = from;
  }
  return 2;                     // magnification factor
}

// Local Hessian magnitude and direction in 3x3 window
// corresponding to largest eigenvalue, in absolute magnitude.
// Encode the sign of the curvature in the angle (dirx, diry).

void
LocalHessian(const gevd_bufferxy& smooth, const int i, const int j,
             float& mag, float& dirx, float& diry)
{
  float two_pij = 2 * floatPixel(smooth, i, j);
  float ddx = floatPixel(smooth, i+1, j) + floatPixel(smooth, i-1, j) - two_pij;
  float ddy = floatPixel(smooth, i, j+1) + floatPixel(smooth, i, j-1) - two_pij;
  float two_dxdy = (floatPixel(smooth, i+1, j+1) +
                    floatPixel(smooth, i-1, j-1) -
                    floatPixel(smooth, i-1, j+1) -
                    floatPixel(smooth, i+1, j-1)) / 2;
  float ddx_plus_ddy = ddx + ddy;
  float ddx_minus_ddy = ddx - ddy;
  float theta = (two_dxdy==0 && ddx_minus_ddy==0) ? 0 : // DOMAIN cond. on atan2
                (float)vcl_atan2(two_dxdy, ddx_minus_ddy) / 2; // modulo PI
  if (ddx_plus_ddy < 0) {
    mag = - ddx_plus_ddy        // most negative eigenvalue
          + vcl_sqrt((ddx_minus_ddy * ddx_minus_ddy) + (two_dxdy * two_dxdy));
    theta += (float)vnl_math::pi_over_2;// angle in range [0 pi]
  } else {
    mag = + ddx_plus_ddy        // most positive eigenvalue
          + vcl_sqrt((ddx_minus_ddy * ddx_minus_ddy) + (two_dxdy * two_dxdy));
    if (theta > 0)
      theta -= (float)vnl_math::pi;    // angle in range [-pi 0]
  }
  dirx = (float)vcl_cos(theta); // eigenvector corresponding to
  diry = (float)vcl_sin(theta); // largest eigenvalue/curvature
}

//: Compute the Hessian of the intensity surface.
// O(m*n).
// The Hessian is directional, and described by the largest absolute
// eigenvalue, and its corresponding eigenvector.
// The intensity surface is assumed smoothed by a Gaussian filter,
// or convolved with a low-pass filter.
// Return at each pixel, the largest absolute value of the two
// eigenvalues, the corresponding eigenvector, and the
// multiplication factor to normalize the magnitude.

float
gevd_float_operators::Hessian(const gevd_bufferxy& smooth,
                              gevd_bufferxy*& magnitude,
                              gevd_bufferxy*& dirx, gevd_bufferxy*& diry,
                              const bool xwrap, const bool ywrap)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Compute local Hessian magnitude/direction";
#endif
  magnitude = gevd_float_operators::Allocate(magnitude, smooth);
  dirx = gevd_float_operators::Allocate(dirx, smooth);
  diry = gevd_float_operators::Allocate(diry, smooth);

  // 1. Inside image frame, compute values based on 3x3 window.
  const int frame = 1;
  const int highx = smooth.GetSizeX() - frame;  // exclusive bounds
  const int highy = smooth.GetSizeY() - frame;
  for (int j = frame; j < highy; ++j)
    for (int i = frame; i < highx; ++i)
      LocalHessian(smooth, i, j,
                   floatPixel(*magnitude, i, j),
                   floatPixel(*dirx, i, j),
                   floatPixel(*diry, i, j));

  // 2. Along horizontal/vertical borders
  if (xwrap) {                  //  each row wrap at border
    const int lo = 2, hi = 1;
    gevd_bufferxy* pad = gevd_float_operators::WrapAlongX(smooth);
    for (int j = 1; j < highy; ++j) {
      LocalHessian(*pad, lo, j,
                   floatPixel(*magnitude, 0, j),
                   floatPixel(*dirx, 0, j),
                   floatPixel(*diry, 0, j));
      LocalHessian(*pad, hi, j,
                   floatPixel(*magnitude, highx, j),
                   floatPixel(*dirx, highx, j),
                   floatPixel(*diry, highx, j));
    }
    delete pad;
  } else {                      // zero by default
    gevd_float_operators::FillFrameX(*magnitude, 0, frame);
    gevd_float_operators::FillFrameX(*dirx, 0, frame);
    gevd_float_operators::FillFrameX(*diry, 0, frame);
  }
  if (ywrap) {                  // each column wrap at border
    const int lo = 2, hi = 1;
    gevd_bufferxy* pad = gevd_float_operators::WrapAlongY(smooth);
    for (int i = 1; i < highx; ++i) {
      LocalHessian(*pad, i, lo,
                   floatPixel(*magnitude, i, 0),
                   floatPixel(*dirx, i, 0),
                   floatPixel(*diry, i, 0));
      LocalHessian(*pad, i, hi,
                   floatPixel(*magnitude, i, highy),
                   floatPixel(*dirx, i, highy),
                   floatPixel(*diry, i, highy));
    }
    delete pad;
  } else {                      // zero by default
    gevd_float_operators::FillFrameY(*magnitude, 0, frame);
    gevd_float_operators::FillFrameY(*dirx, 0, frame);
    gevd_float_operators::FillFrameY(*diry, 0, frame);
  }
#ifdef DEBUG
  vcl_cout << ", in " << t.real() << " msecs.\n";
#endif
  return 2;                     // multiplication factor for magnitude
}

// Local Laplacian in 3x3 neighborhood, sum of curvature,
// and eigenvector corresponding to largest absolute eigenvalue.

void
LocalLaplacian(const gevd_bufferxy& smooth, const int i, const int j,
               float& mag, float& dirx, float& diry)
{
  float two_pij = 2 * floatPixel(smooth, i, j);
  float ddx = floatPixel(smooth, i+1, j) + floatPixel(smooth, i-1, j) - two_pij;
  float ddy = floatPixel(smooth, i, j+1) + floatPixel(smooth, i, j-1) - two_pij;
  float diag1 = floatPixel(smooth, i+1, j+1) + floatPixel(smooth, i-1, j-1);
  float diag2 = floatPixel(smooth, i-1, j+1) + floatPixel(smooth, i+1, j-1);
  mag = 4*(ddx + ddy) - 2*two_pij + diag1 + diag2; // save division by 6
#if 0 // commented out
  mag = (floatPixel(smooth, i+1, j+1) + floatPixel(smooth, i-1, j-1) +
         floatPixel(smooth, i-1, j+1) + floatPixel(smooth, i+1, j-1) +
         4 * (floatPixel(smooth, i+1, j) + floatPixel(smooth, i-1, j) +
              floatPixel(smooth, i, j+1) + floatPixel(smooth, i, j-1)) +
         -20 * floatPixel(smooth, i, j)) / 6;
#endif
  float theta = (diag1==diag2 && ddx==ddy) ? 0 : // DOMAIN condition on atan2
                (float)vcl_atan2((diag1 - diag2) / 2, ddx - ddy) / 2; // modulo PI
  if (mag < 0) {
    mag = - mag;                  // absolute magnitude
    theta += (float)vnl_math::pi_over_2; // other eigenvector
  }
  dirx = (float)vcl_cos(theta);   // eigenvector corresponding to
  diry = (float)vcl_sin(theta);   // largest eigenvalue/curvature
}

//: Compute the Laplacian of the intensity surface.
// O(m*n).
// The Laplacian is rotational symmetric, and is the sum of the 2
// eigenvalues/curvatures, with positive/negative sign for concave/convex.
// The intensity surface is assumed smoothed by a Gaussian filter,
// or convolved with a low-pass filter.
// Return at each pixel, the absolute value of the laplacian,
// the eigenvector corresponding to the largest absolute
// eigenvalue/curvature, and the multiplication factor to
// normalize the magnitude.

float
gevd_float_operators::Laplacian(const gevd_bufferxy& smooth,
                                gevd_bufferxy*& magnitude,
                                gevd_bufferxy*& dirx, gevd_bufferxy*& diry,
                                const bool xwrap, const bool ywrap)
{
#ifdef DEBUG
  vul_timer t;
  vcl_cout << "Compute local Laplacian magnitude/direction";
#endif
  magnitude = gevd_float_operators::Allocate(magnitude, smooth);
  dirx = gevd_float_operators::Allocate(dirx, smooth);
  diry = gevd_float_operators::Allocate(diry, smooth);

  // 1. Inside image frame, compute values based on 3x3 window.
  const int frame = 1;
  const int highx = smooth.GetSizeX() - frame;// exclusive bounds
  const int highy = smooth.GetSizeY() - frame;
  for (int j = frame; j < highy; ++j)
    for (int i = frame; i < highx; ++i)
      LocalLaplacian(smooth, i, j,
                     floatPixel(*magnitude, i, j),
                     floatPixel(*dirx, i, j),
                     floatPixel(*diry, i, j));

  // 2. Along horizontal/vertical borders
  if (xwrap) {                  //  each row wrap at border
    const int lo = 2, hi = 1;
    gevd_bufferxy* pad = gevd_float_operators::WrapAlongX(smooth);
    for (int j = 1; j < highy; ++j) {
      LocalLaplacian(*pad, lo, j,
                     floatPixel(*magnitude, 0, j),
                     floatPixel(*dirx, 0, j),
                     floatPixel(*diry, 0, j));
      LocalLaplacian(*pad, hi, j,
                     floatPixel(*magnitude, highx, j),
                     floatPixel(*dirx, highx, j),
                     floatPixel(*diry, highx, j));
    }
    delete pad;
  } else {                      // zero by default
    gevd_float_operators::FillFrameX(*magnitude, 0, frame);
    gevd_float_operators::FillFrameX(*dirx, 0, frame);
    gevd_float_operators::FillFrameX(*diry, 0, frame);
  }
  if (ywrap) {                  // each column wrap at border
    const int lo = 2, hi = 1;
    gevd_bufferxy* pad = gevd_float_operators::WrapAlongY(smooth);