%ffp

// Recursive gaussian demo code
// implemented by Tom Fiddaman, www.sd3.info/pf828
// from the article: 
//    Title: Recursive implementation of the Gaussian filter
//    Authors: I.T. Young, L.J. van Vliet
//    in: Signal Processing, vol. 44, no. 2, 1995, 139-151.
// http://www.ph.tn.tudelft.nl/~lucas/publications/1995/SP95TYLV/SP95TYLV.html

// changes by Alois Zingl:
// continues boundary condition added instead of reflected boundaries
//    (the boarder pixel is copied beyond the boarder)
// the image is processed in tiles

Category: "easy.Filter"
Title:    "Recursive Gaussian"
Version:  "1.0"
Filename: "recursivegauss.8bf"
Author:   "Tom Fiddaman, Alois Zingl"
URL:      "http://free.pages.at/easyfilter/gauss.html"
Description:"Recursive implementation of the Gaussian filter."

ctl(0): "Blur (pixel)"

OnFilterstart:
{
  isTileable = !doingProxy; // split large images
  needPadding = 4*(doingProxy?(ctl(0)+scaleFactor/2)/scaleFactor:ctl(0));
  bandWidth = 100+4*needPadding;
  return false;
}

ForEveryTile:
{
  int radius = needPadding;
  double b3 = max(radius-4, radius*2/3)/4.0, wn, wn1, wn2, wn3; 
  double b0 = 1.57825 + b3*(2.44413 + b3*(1.4281 + 0.422205*b3));
  double b1 = b3*(2.44413 + b3*(2.85619 + 1.26661*b3))/b0;
  double b2 = (-1.4281 - 1.26661*b3)*b3*b3/b0;
  b3 = 0.422205*b3*b3*b3/b0;
  b0 = 1.0 - (b1+b2+b3);

  for (z = 0; z < Z; z++)  // loop through all channels
  { 
    //Update progress bar and cancel if ESC key was pressed
    if (updateProgress((1-doingProxy)*y_start*Z*2+(y_end-y_start)*(2*z+1), 
                       (doingProxy?y_end-y_start:Y)*Z*2))
       abort();

    // DO COLUMNS
  
    for (x = x_start; x < x_end; x++)       // loop through all rows
    {
      y = y_start-radius;
      wn1 = wn2 = wn3 = src(x,y,z); // initialize value
  
      while (++y < y_start)  // initialize row
      {   
        wn = b0*src(x,y,z) + b1*wn1+b2*wn2+b3*wn3;
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 

      while (++y < y_end)   // real forward pass
      {
        wn = b0*src(x,y,z) + b1*wn1+b2*wn2+b3*wn3;
        pset(x,y,z, (int)wn);
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 

      while (++y < y_end+radius)  // post process forward 
      {
        wn = b0*src(x,y,z) + b1*wn1+b2*wn2+b3*wn3;
        put((int)wn, y);
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 
      wn1 = wn2 = wn3 = src(x,y-1,z); // initialize values
    
      while (--y > y_end)   // initialize backward
      {
        wn = b0*get(y) + b1*wn1+b2*wn2+b3*wn3;
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 

      while (y-- > y_start)   // real backward pass
      {
        wn = b0*pget(x,y,z) + b1*wn1+b2*wn2+b3*wn3;
        pset(x,y,z, (int)wn);
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 
    } // for x
  
    //Update progress bar and cancel if ESC key was pressed
    if (updateProgress((1-doingProxy)*y_start*Z*2+(y_end-y_start)*(2*z+2), 
                       (doingProxy?y_end-y_start:Y)*Z*2))
       abort();
  
    // DO ROWS (a bit simpler since rows aren't split in tiles)
    
    for (y = y_start; y < y_end; y++)   // loop through all columns
    {
      int p = pget(x_end-1,y,z);   // remember boarder pixel
      wn1 = wn2 = wn3 = pget(x_start,y,z);
  
      for (x = x_start; x < x_end; x++)  // real forward pass
      {   
        wn = b0*pget(x,y,z) + b1*wn1+b2*wn2+b3*wn3;
        pset(x,y,z, (int)wn);
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 

      for (; x < x_end+radius; x++)  // post process forward
      {   
        wn = b0*p + b1*wn1+b2*wn2+b3*wn3;
        put((int)wn, x);
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 
      wn1 = wn2 = wn3 = p; // initialize values
  
      while (--x > x_end)   // initialize backward
      {   
        wn = b0*get(x) + b1*wn1+b2*wn2+b3*wn3;
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 

      while (x-- > x_start)  // real backward pass
      {   
        wn = b0*pget(x,y,z) + b1*wn1+b2*wn2+b3*wn3;
        pset(x,y,z, (int)wn);
        wn3 = wn2; wn2 = wn1; wn1 = wn;
      } 
    } // for y
  } // for z

  //Stop processing and apply the effect
  return true;
}

OnFilterEnd:
{
  updateProgress(0, 1);
  return false;
}