//==============================================
//  copyright            : (C) 2003-2005 by Will Stokes
//==============================================
//  This program is free software; you can redistribute it
//  and/or modify it under the terms of the GNU General
//  Public License as published by the Free Software
//  Foundation; either version 2 of the License, or
//  (at your option) any later version.
//==============================================

//Systemwide includes
#include <qimage.h>
#include <qstring.h>
#include <math.h>

//Projectwide includes
#include "blur.h"

//----------------------------------------------
// Inputs:
// -------
// QImage& - image to be blurred
// Float sigma - amount of blurring desired
//
// Description:
// ------------
// This method blurs a color image using the approach published
// by Ian T. Young and Lucas J. van Vliet, "Recursive implementation
// of the Gaussian filter", 1995 Signal Processing.
//
// The first solution one might take to blurring an image
// is to convolve a gaussian filter with an image. This is expensive since
// as the gaussian filter matrix grows, processing becomes quite expensive
// ( O(M^2 * N^2) where M is the width/height of the gaussian matrix and N is the 
// width/height of the image. )
//
// An alternative is to approximate a gaussian by repeatedly convolving a kernel 
// such as uniform filter. While faster, this is only an approximation of a 
// gaussian, and the weights one uses for each iteration must be hand picked in
// order to follow the Well's equation:
// (W_1^2 + W_2^2 + W_3^2 + W_4^2 = 12*sigma^2 + 4)
//
// This process is tedious and limits us in applying intermediate blur amounts.
//
// In theory Young and Vliet's technique (which is O(4N^2) = O(N^2)) is not only
// more accurate but faster as well. The only penalty we pay is the extra storage
// of a row or column while that row or column is being processed. For a 
// 2048x2048 image my implementation will require an extra 16kb for this extra 
// storage space, a drop in the bucket compared to the memory required to hold
// the entire image in float memory while being processed.
//
// Before being blurred color values are converted from integer [0-255] to float
// [0.0-1.0] space. This is necessary to avoid compound rounding errors. In order
// to minimize space required to store the float buffer, the image is processed 
// in three passes for each of the red, green, and blue color channels.
//
// Since the code is broken up nicely, well commented, and is a direct 
// implementation of Young and Vliet's techinque I point you to their paper for a
// better understanding of why/how this techinque actually works. Vliet has a 
// copy on his personal web site, last seen at:
// http://www.ph.tn.tudelft.nl/~lucas/
//----------------------------------------------

void computeCoeffs( float sigma );
void fillBuffer( QImage &image, int channel );
void blurBuffer();

void blurRow( int row );
void blurColumn( int column );

void blurRegionsInRow( int y );
void blurRegionsInCol( int x );

void resetImageData( QImage &image, int channel, bool blurEdges);

float edgeValue(int x, int y);

float q, b0, b1, b2, b3, B;
int width, height;
float* buffer;
float* rowBuffer;
float* colBuffer;

float* regionRowBuffer;
float* regionColBuffer;

QImage* edgeImage;
int* regionMap;
int regionCount;
QPoint displayOffset;
QSize fullRes;

//==============================================
void blurImage( QImage &image, float sigma )
{
  //supply dummy data for edges, notably NULL for the edge image pointer.
  //other values have no effect
  blurImage( image, sigma, QPoint(0,0), image.size(), NULL, NULL, 0, false );
}
//==============================================
void blurImage( QImage &image, float sigma,
                QPoint offset, QSize fullImageRes,
                QImage* edges, int* regions, int numRegions,
                bool targetEdges)
{
  edgeImage = edges;
  regionMap = regions;
  regionCount = numRegions;
  displayOffset = offset;
  fullRes = fullImageRes;
  
  //compute blurring coeffecients
  computeCoeffs(sigma);
  
  //store image dimensions
  width = image.width();
  height = image.height();
  
  //Construct float buffer that is the size of the image/
  //In order to conserve memory process image three times, once for
  //each color channel.
  buffer = new float[ width * height ];

  rowBuffer = new float[width];
  colBuffer = new float[height];
  
  regionRowBuffer = new float[width * numRegions];
  regionColBuffer = new float[height * numRegions];
  
  //iterate over each color channel
  int channel;
  for( channel = 0; channel <=2; channel++)
  {
    //copy color data into float buffer
    fillBuffer( image, channel );
    
    //blur buffer data
    blurBuffer();
    
    //reset image data used blurred buffer
    resetImageData(image, channel, targetEdges);
  }
  
  //delete buffer
  delete[] buffer;
  delete[] rowBuffer;
  delete[] colBuffer;
}
//==============================================
void computeCoeffs( float sigma )
{
  //compute q as a function of sigma
  if( sigma >= 2.5f )
  {
    q = 0.98711f*sigma - 0.96330f;  
  }
  else
  {
    q = 3.97156f - 4.14554f * sqrt( 1.0f - 0.26891f*sigma );
  }

  //compute b0, b1, b2, and b3
  b0 = 1.57825f + (2.44413f*q) + (1.4281f * q*q ) + (0.422205f * q*q*q );
  b1 = (2.44413f * q) + (2.85619f * q*q) + (1.26661 * q*q*q );
  b2 = -((1.4281 * q*q) + (1.26661 * q*q*q));
  b3 = 0.422205 * q*q*q;
  
  //compute B
  B = 1.0f - ((b1 + b2 + b3) / b0); 
}
//==============================================
void fillBuffer( QImage &image, int channel )
{
  //precompute 1/255
  float multiplier = 1.0f / 255.0f;
  
  //iterate over each selected scanline 
  int x, y;
  QRgb* rgb;
  uchar* scanLine;
  for( y=0; y<image.height(); y++)
  {   
    //iterate over each pixel in scanline
    scanLine = image.scanLine(y);
    for( x=0; x<image.width(); x++)
    {
      //get handle on rgb value in image
      rgb = ((QRgb*)scanLine+x);

      //compute index where float value is stored in buffer
      int index = x + y*image.width();
      
      //convert and store correct channel in buffer
      if( channel == 0 )
        buffer[index] = multiplier * qRed( *rgb );
      else if( channel == 1 )
        buffer[index] = multiplier * qGreen( *rgb );
      else
        buffer[index] = multiplier * qBlue( *rgb );
    } //x
  } //y
}
//==============================================
void blurBuffer()
{
  //blur rows, then columns
  int index;
  
  if(regionMap == NULL || edgeImage == NULL )
  {
    for(index=0; index < height; index++)
    { blurRow( index ); }
    
    for(index=0; index< width; index++)
    { blurColumn( index ); }
  }
  else
  {
    for(index=0; index < height; index++)
    { blurRegionsInRow( index ); }
    
    for(index=0; index< width; index++)
    { blurRegionsInCol( index ); }
  }
}
//==============================================
int regionIndex(int x, int y)
{
  int edgeX = ((edgeImage->width()-1) * (x+displayOffset.x())) / (fullRes.width()-1);
  int edgeY = ((edgeImage->height()-1) * (y+displayOffset.y())) / (fullRes.height()-1);
  return edgeY*edgeImage->width() + edgeX;
}
//==============================================
float edgeValue(int x, int y)
{
  //compute floating point x and y coordinates for edge image
  float edgeX = ((edgeImage->width()-1.0f) * (x+displayOffset.x())) / (fullRes.width()-1);
  float edgeY = ((edgeImage->height()-1.0f) * (y+displayOffset.y())) / (fullRes.height()-1);
  
  //compute 4 int values of coordinates
  int x1 = (int)edgeX;
  int y1 = (int)edgeY;
  int x2, y2;
  if( edgeX > x1 )
    x2 = x1+1;
  else
    x2 = x1;
  if( edgeY > y1 )
    y2 = y1+1;
  else
    y2 = y1;
  
  //compute the four indices
  int index1, index2, index3, index4;
  index1 = x1 + y1*edgeImage->width();
  index2 = x2 + y1*edgeImage->width();
  index3 = x1 + y2*edgeImage->width();
  index4 = x2 + y2*edgeImage->width();
  
  //find edge quantity for each corner
  float v1, v2, v3, v4; 
  uchar* scanline = edgeImage->scanLine( y1 );
  QRgb* rgb = ((QRgb*)scanline+x1);
  v1 = ((float) qRed( *rgb )) / 255.0f;
  rgb = ((QRgb*)scanline+x2);
  v2 = ((float) qRed( *rgb )) / 255.0f;

  scanline = edgeImage->scanLine( y2 );
  rgb = ((QRgb*)scanline+x1);
  v3 = ((float) qRed( *rgb )) / 255.0f;
  rgb = ((QRgb*)scanline+x2);
  v4 = ((float) qRed( *rgb )) / 255.0f;
  
  //blur combine left-right
  v1 = (edgeX-x1)*v2 + (1 - edgeX + x1)*v1;
  v3 = (edgeX-x1)*v4 + (1 - edgeX + x1)*v3;
  
  //combine top-bottom
  v1 = (edgeY-y1)*v3 + (1 - edgeY + y1)*v1;
  
  //return result
  return v1;
}
//==============================================
void blurRow( int row )
{
  int i;
  int rtw = row*width;
  
  //forward
  rowBuffer[0] = buffer[ 0 + rtw ];
  for(i=1; i<width; i++)
  {
     rowBuffer[i] = B*buffer[ i + rtw ] +
     ( b1*rowBuffer[ QMAX(i-1, 0) ] +
       b2 * rowBuffer[ QMAX(i-2, 0) ] +
       b3 * rowBuffer[ QMAX(i-3, 0) ]) / b0;    
    
  }
  
  //reverse
  for(i=width-1; i>=0; i--)
  {
     buffer[ i + rtw ] = B*rowBuffer[ i ] +
     ( b1 * buffer[ QMIN(i+1, width-1) + rtw ] +
       b2 * buffer[ QMIN(i+2, width-1) + rtw ] +
       b3 * buffer[ QMIN(i+3, width-1) + rtw ]) / b0;    
  }
}
//==============================================
void blurRegionsInRow( int y )
{
  //---------------------------------
  //populate region row buffer. a row has been allocated for 
  //each region. Pixels between regions
  //take the closest pixel value in that row from that region
  int yTimesWidth = y*width;
  int regionTimesWidth;
  int region,x,x2;
  
  //for each region
  for(region=0; region<regionCount; region++)
  {
    regionTimesWidth = region*width;
    int lastX = -1;
    for(x=0; x<width; x++)
    {
      //if pixel belongs to this region then update lastX index and copy value over
      //if lastX is mroe than one pixel away than fill inbetween region
      if( region == regionMap[regionIndex(x, y)] )
      {
        //fill empty region preceeding this region blob
        if( lastX < x-1)
        {
          //no preceeding region, spread left!
          if(lastX == -1)
          {
            for(x2=0; x2<x; x2++) { regionRowBuffer[x2 + regionTimesWidth] = buffer[x + yTimesWidth]; }
          }
          //else spread from both left and right of empty stretch
          else
          {
            int xMid = lastX + ((x-1) - lastX)/2;
            
            for(x2=lastX+1; x2<=xMid; x2++)
            { regionRowBuffer[x2 + regionTimesWidth] = buffer[lastX + yTimesWidth]; }

            for(x2=xMid+1; x2<x; x2++)
            { regionRowBuffer[x2 + regionTimesWidth] = buffer[x + yTimesWidth]; }
          }
        }
        
        regionRowBuffer[x + regionTimesWidth] = buffer[x + yTimesWidth];
        lastX = x;
      } 
    } //x
    
    //if last stretch is empty, fill right
    if( region != regionMap[regionIndex(width-1, y)] )
    {
      for(x2=lastX+1; x2<width; x2++)
      { regionRowBuffer[x2 + regionTimesWidth] = buffer[lastX + yTimesWidth]; }
    }
    
  } //region
  //---------------------------------
  //blur the region row buffers

  //for each region
  for(region=0; region<regionCount; region++)
  {
    regionTimesWidth = region*width;
    
    //forward
    rowBuffer[0] = regionRowBuffer[ 0 + regionTimesWidth ];
    for(x=1; x<width; x++)
    {
      rowBuffer[x] = B*regionRowBuffer[ x + regionTimesWidth ] +
      ( b1*rowBuffer[ QMAX(x-1, 0) ] +
        b2 * rowBuffer[ QMAX(x-2, 0) ] +
        b3 * rowBuffer[ QMAX(x-3, 0) ]) / b0;    
    }
    
    //reverse
    for(x=width-1; x>=0; x--)
    {
      regionRowBuffer[ x + regionTimesWidth ] = B*rowBuffer[ x ] +
      ( b1 * regionRowBuffer[ QMIN(x+1, width-1) + regionTimesWidth ] +
        b2 * regionRowBuffer[ QMIN(x+2, width-1) + regionTimesWidth ] +
        b3 * regionRowBuffer[ QMIN(x+3, width-1) + regionTimesWidth ]) / b0;    
    }
  }
  //---------------------------------
  //copy data from the region row buffers back to the
  //buffer. for each pixel we choose the correct region
  //row buffer basedon the original regionidentity of hte pixel
  for(x=0; x<width; x++)
  {
    int ri = regionIndex(x,y);
    int region = regionMap[ri];
    float bufferVal = regionRowBuffer[ x + region*width ];
    buffer[x + yTimesWidth] = bufferVal;
    
//    buffer[x + yTimesWidth] = regionRowBuffer[x + regionMap[regionIndex(x,y)]*width]; 
  }
  //---------------------------------
}
//==============================================
void blurColumn( int column )
{
  int i;

  //forward
  colBuffer[0] = buffer[ column + 0*width ];
  for(i=1; i<height; i++)
  {
     colBuffer[i] = B*buffer[ column + i*width ] +
     ( b1 * colBuffer[ QMAX(i-1, 0) ] +
       b2 * colBuffer[ QMAX(i-2, 0) ] +
       b3 * colBuffer[ QMAX(i-3, 0) ]) / b0;    
  }
  
  //reverse
  for(i=height-1; i>=0; i--)
  {
     buffer[ column + i*width ] = B*colBuffer[ i ] +
     ( b1 * buffer[ column + QMIN(i+1, height-1)*width ] +
       b2 * buffer[ column + QMIN(i+2, height-1)*width ] +
       b3 * buffer[ column + QMIN(i+3, height-1)*width ]) / b0;        
  }
  
}
//==============================================
void blurRegionsInCol( int x )
{
  //---------------------------------
  //populate region col buffer. a col has been allocated for 
  //each region. Pixels between regions
  //take the closest pixel value in that col from that region
//  int yTimesWidth = y*width;
  int regionTimesHeight;
  int region,y,y2;
  
  //for each region
  for(region=0; region<regionCount; region++)
  {
    regionTimesHeight = region*height;
    int lastY = -1;
    for(y=0; y<height; y++)
    {
      //if pixel belongs to this region then update lastY index and copy value over
      //if lastY is more than one pixel away than fill inbetween region
      if( region == regionMap[regionIndex(x, y)] )
      {
        //fill empty region preceeding this region blob
        if( lastY < y-1)
        {
          //no preceeding region, spread left!
          if(lastY == -1)
          {
            for(y2=0; y2<y; y2++) { regionColBuffer[y2 + regionTimesHeight] = buffer[x + y*width]; }
          }
          //else spread from both left and right of empty stretch
          else
          {
            int yMid = lastY + ((y-1) - lastY)/2;
            
            for(y2=lastY+1; y2<=yMid; y2++)
            { regionColBuffer[y2 + regionTimesHeight] = buffer[x + lastY*width]; }
            
            for(y2=yMid+1; y2<y; y2++)
            { regionColBuffer[y2 + regionTimesHeight] = buffer[x + y*width]; }
          }
        }
        
        regionColBuffer[y + regionTimesHeight] = buffer[x + y*width];
        lastY = y;
      } 
    } //y
    
    //if last stretch is empty, fill right
    if( region != regionMap[regionIndex(x, height-1)] )
    {
      for(y2=lastY+1; y2<height; y2++)
      { regionColBuffer[y2 + regionTimesHeight] = buffer[x + lastY*width]; }
    }
    
  } //region
  //---------------------------------
  //blur the region col buffers
  
  //for each region
  for(region=0; region<regionCount; region++)
  {
    regionTimesHeight = region*height;
    
    //forward
    colBuffer[0] = regionColBuffer[ 0 + regionTimesHeight ];
    for(y=1; y<height; y++)
    {
      colBuffer[y] = B*regionColBuffer[ y + regionTimesHeight ] +
      ( b1 * colBuffer[ QMAX(y-1, 0) ] +
        b2 * colBuffer[ QMAX(y-2, 0) ] +
        b3 * colBuffer[ QMAX(y-3, 0) ]) / b0;    
    }
    
    //reverse
    for(y=height-1; y>=0; y--)
    {
      regionColBuffer[ y + regionTimesHeight ] = B*colBuffer[ y ] +
      ( b1 * regionColBuffer[ QMIN(y+1, height-1) + regionTimesHeight ] +
        b2 * regionColBuffer[ QMIN(y+2, height-1) + regionTimesHeight ] +
        b3 * regionColBuffer[ QMIN(y+3, height-1) + regionTimesHeight ]) / b0;    
    }
  }
  //---------------------------------
  //copy data from the region row buffers back to the
  //buffer. for each pixel we choose the correct region
  //row buffer basedon the original regionidentity of hte pixel
  for(y=0; y<height; y++)
  {
    buffer[x + y*width] = regionColBuffer[y + regionMap[regionIndex(x,y)]*height]; 
  }
  //---------------------------------
}
//==============================================
void resetImageData( QImage &image, int channel, bool blurEdges)
{
  //iterate over each selected scanline 
  int x, y;
  QRgb *rgb;
  uchar* imageScanline = NULL;
  for( y=0; y<image.height(); y++)
  {   
    imageScanline = image.scanLine(y);
    for( x=0; x<image.width(); x++)
    {
      //get handle on rgb value in image
      rgb = ((QRgb*)imageScanline+x);
      
      //compute index where float value is stored in buffer
      int index = x + y*image.width();

      //convert blured value to 0-255 range
      int blurredColor = QMAX( QMIN( ((int) (255*buffer[index])), 255 ), 0 );
      
      //blur the entire thing!
      float alpha;
      if( edgeImage == NULL)
        alpha = 1.0f;
      else
      {
        alpha = edgeValue( x, y ); 
        if(!blurEdges)
          alpha = 1.0f - alpha;
      }

      //convert and store correct channel in buffer
      if( channel == 0 )
        *rgb = qRgb( (int) (alpha*blurredColor + (1-alpha)*qRed(*rgb)),
                     qGreen(*rgb), qBlue(*rgb) );
      else if( channel == 1 )
        *rgb = qRgb( qRed(*rgb), 
                     (int) (alpha*blurredColor + (1-alpha)*qGreen(*rgb)),
                     qBlue(*rgb) );
      else
        *rgb = qRgb( qRed(*rgb), qGreen(*rgb), 
                     (int) (alpha*blurredColor + (1-alpha)*qBlue(*rgb)) );
    } //x
  } //y
}
//==============================================