/* Selective gaussian blur filter for GIMP, version 0.1
* Adapted from the original gaussian blur filter by Spencer Kimball and
* Peter Mattis.
*
* Copyright (C) 1995 Spencer Kimball and Peter Mattis
* Copyright (C) 1999 Thom van Os <thom@vanos.com>
* Copyright (C) 2006 Loren Merritt
*
* 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.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*
* To do:
*      - support for horizontal or vertical only blur
*      - use memory more efficiently, smaller regions at a time
*      - integrating with other convolution matrix based filters ?
*      - create more selective and adaptive filters
*      - threading
*      - optimization
*/

// This algorithm was taken from the GIMP v2.4.5 sourcecodes
// and adopted for the FreeImage library
//
//  Copyright (C) 2007-2008:
//  monday2000  monday2000@yandex.ru

#include "FreeImage.h"
#include "Utilities.h"

#define PI  ((double)3.14159265358979323846264338327950288419716939937510)

/////////////////////////////////////////////////////////////////////////////////////////////

static void init_matrix (double radius, double *mat, int num)
{
    int    dx;
    double sd, c1, c2;
    
    // This formula isn't really correct, but it'll do
    sd = radius / 3.329042969;
    c1 = 1.0 / sqrt (2.0 * PI * sd);
    c2 = -2.0 * (sd * sd);
    
    for (dx = 0; dx < num; dx++)
        mat[dx] = c1 * exp ((dx * dx)/ c2);
}

/////////////////////////////////////////////////////////////////////////////////////////////

static FIBITMAP* sel_gauss (FIBITMAP* src_dib, double radius, int maxdelta)
{
	double	   *mat;
	int 		numrad;
	
	unsigned bpp = FreeImage_GetBPP(src_dib);
	
	unsigned btpp = bpp/8;
	
	unsigned width = FreeImage_GetWidth(src_dib);
	
	unsigned height = FreeImage_GetHeight(src_dib);
	
	unsigned pitch = FreeImage_GetPitch(src_dib);
	
	unsigned d;
	
	FIBITMAP* dst_dib = FreeImage_Allocate(width, height, bpp);	
	
	BYTE* src_bits = (BYTE*)FreeImage_GetBits(src_dib); // The image raster
	
	BYTE* dst_bits = (BYTE*)FreeImage_GetBits(dst_dib); // The image raster
	
	BYTE* lines, *linek, *lined;
	
	
	numrad = (int)(radius + 1.0);
	
	mat = (double*) calloc (numrad, sizeof(double));
	
	init_matrix(radius, mat, numrad);
	
    unsigned short    *imat;
    double     fsum, fscale;
    int        i, j, x, y;
	
    
    imat = (unsigned short*)calloc (2 * numrad, sizeof(unsigned short));
    
    fsum = 0.0;
    for (y = 1 - numrad; y < numrad; y++)
        fsum += mat[abs(y)];
    
    // Ensure that the sum fits in 32bits.
    fscale = 0x1000 / fsum;
    for (y = 0; y < numrad; y++)
        imat[numrad - y] = imat[numrad + y] = mat[y] * fscale;
	
	/////////////////////////////////////////
	
    unsigned* p_sum = new unsigned[btpp];
    unsigned* p_fact = new unsigned[btpp];
	
    unsigned* p_rowsum = new unsigned[btpp];
    unsigned* p_rowfact = new unsigned[btpp];
	
	int di;
	
	int tmp;
	
	int flag = 0;
	
	/////////////////////////////////////////
    
    for (y = 0; y < height; y++)
    {
		lines = src_bits + y * pitch;
		
		lined = dst_bits + y * pitch;
		
        for (x = 0; x < width; x++)
        {	
			memset(p_sum, 0, sizeof(unsigned)*btpp);
			memset(p_fact, 0, sizeof(unsigned)*btpp);
			
			linek = src_bits + pitch * (y - numrad);		
			
			for (j = 1 - numrad; j < numrad; j++)
			{
				memset(p_rowsum, 0, sizeof(unsigned)*btpp);
				memset(p_rowfact, 0, sizeof(unsigned)*btpp);
				
				linek += pitch;
				
				if (y + j < 0 || y + j >= height)
					continue;				
				
				for (i = 1 - numrad; i < numrad; i++)
				{
					if ((x + i) < 0 || (x + i) >= width)					
						continue;
					
					for (d=0; d<btpp; d++)				
					{					
						
						tmp = lines[x * btpp + d] - linek[(x + i) * btpp + d];
						
						if (tmp > maxdelta || tmp < -maxdelta)
						{
							continue;
						}
						
						di = imat[numrad+i];
						
						p_rowsum[d] += di * linek[(x + i) * btpp + d];
						p_rowfact[d] += di;			
						
					}
				}
				
				di = imat[numrad+j];
				
				for (d=0; d<btpp; d++)			
				{				
					p_sum[d] += di * p_rowsum[d];
					p_fact[d] += di * p_rowfact[d];
				}
			}
			
			for (d=0; d<btpp; d++)			
			{
				if (p_fact[d] == 0)
					lined[x * btpp + d] = lines[x * btpp + d];
				else
					lined[x * btpp + d] = p_sum[d] / p_fact[d];
			}
		}
    }
    
	// free up buffers
    free (imat);
	free (mat);
	
	delete [] p_sum;
	
	delete [] p_fact;
	
	delete [] p_rowsum;
	
	delete [] p_rowfact;
	
	if(bpp == 8)
	{
		// copy the original palette to the destination bitmap
		
		RGBQUAD *src_pal = FreeImage_GetPalette(src_dib);
		RGBQUAD *dst_pal = FreeImage_GetPalette(dst_dib);
		memcpy(&dst_pal[0], &src_pal[0], 256 * sizeof(RGBQUAD));		
	}
	
	// Copying the DPI...
	
	FreeImage_SetDotsPerMeterX(dst_dib, FreeImage_GetDotsPerMeterX(src_dib));
	
	FreeImage_SetDotsPerMeterY(dst_dib, FreeImage_GetDotsPerMeterY(src_dib));
	
	return dst_dib;
}

/////////////////////////////////////////////////////////////////////////////////////////////
// ----------------------------------------------------------
/**
FreeImage error handler
@param fif Format / Plugin responsible for the error 
@param message Error message
*/
void FreeImageErrorHandler(FREE_IMAGE_FORMAT fif, const char *message) {
	printf("\n*** "); 
	printf("%s Format\n", FreeImage_GetFormatFromFIF(fif));
	printf(message);
	printf(" ***\n");
}

/////////////////////////////////////////////////////////////////////////////////////////////
// ----------------------------------------------------------

/** Generic image loader

  @param lpszPathName Pointer to the full file name
  @param flag Optional load flag constant
  @return Returns the loaded dib if successful, returns NULL otherwise
*/

FIBITMAP* GenericLoader(const char* lpszPathName, int flag)
{	
	FREE_IMAGE_FORMAT fif = FIF_UNKNOWN;
	// check the file signature and deduce its format
	// (the second argument is currently not used by FreeImage)
	
	fif = FreeImage_GetFileType(lpszPathName, 0);
	
	FIBITMAP* dib;
	
	if(fif == FIF_UNKNOWN)
	{
		// no signature ?
		// try to guess the file format from the file extension
		fif = FreeImage_GetFIFFromFilename(lpszPathName);
	}
	
	// check that the plugin has reading capabilities ...
	if((fif != FIF_UNKNOWN) && FreeImage_FIFSupportsReading(fif))
	{
		// ok, let's load the file
		dib = FreeImage_Load(fif, lpszPathName, flag);
		
		// unless a bad file format, we are done !
		if (!dib)
		{
			printf("%s%s%s\n","File \"", lpszPathName, "\" not found.");
			return NULL;
		}
	}	
	
	return dib;
}

/////////////////////////////////////////////////////////////////////////////////////////////

// Selective Gaussian Blur

// Blur neighboring pixels, but only in low-contrast areas,
// This filter functions similar to the regular
// gaussian blur filter except that neighbouring
// pixels that differ more than the given maxdelta
// parameter will not be blended with. This way with
// the correct parameters, an image can be smoothed
// out without losing details. However, this filter
// can be rather slow.
// Thom van Os,
// 1999

// double radius    // Radius of gaussian blur (in pixels, > 0)
// int max-delta  // Maximum delta

// 5.0  // default radius
// 50   // default max-delta

int main(int argc, char *argv[])
{
	
	// call this ONLY when linking with FreeImage as a static library
#ifdef FREEIMAGE_LIB
	FreeImage_Initialise();
#endif // FREEIMAGE_LIB
	
	// initialize your own FreeImage error handler
	
	FreeImage_SetOutputMessage(FreeImageErrorHandler);
	
	if(argc !=4) {
		printf("Usage : sel_gauss <input_file> <radius (double)> <max_delta (int)>\n");
		return 0;
	}
	
	FIBITMAP *dib = GenericLoader(argv[1], 0);
	
	double radius = atof(argv[2]);
	
	int maxdelta = atof(argv[3]);
	
	if (dib)
	{		
		// bitmap is successfully loaded!
		
		if (FreeImage_GetImageType(dib) == FIT_BITMAP)
		{
			if (FreeImage_GetBPP(dib) == 8 || FreeImage_GetBPP(dib) == 24)
			{
				FIBITMAP* dst_dib = sel_gauss (dib, radius, maxdelta);				
				
				if (dst_dib)
				{					
					// save the filtered bitmap
					const char *output_filename = "filtered.tif";
					
					// first, check the output format from the file name or file extension
					FREE_IMAGE_FORMAT out_fif = FreeImage_GetFIFFromFilename(output_filename);
					
					if(out_fif != FIF_UNKNOWN)
					{
						// then save the file
						FreeImage_Save(out_fif, dst_dib, output_filename, 0);
					}
					
					// free the loaded FIBITMAP
					FreeImage_Unload(dst_dib);					
				}
			}
			
			else
				
				printf("%s\n", "Unsupported color mode.");
		}
		
		else // non-FIT_BITMAP images are not supported.
			
			printf("%s\n", "Unsupported color mode.");
		
		FreeImage_Unload(dib);
	}	 
	
	// call this ONLY when linking with FreeImage as a static library
#ifdef FREEIMAGE_LIB
	FreeImage_DeInitialise();
#endif // FREEIMAGE_LIB
	
	return 0;