// 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
// http://www.gnu.org/copyleft/gpl.html

// AForge Image Processing Library
// AForge.NET framework
//
// Copyright © Andrew Kirillov, 2005-2007
// andrew.kirillov@gmail.com
//
// Found description in
// "An Edge Detection Technique Using the Facet
// Model and Parameterized Relaxation Labeling"
// by Ioannis Matalas, Student Member, IEEE, Ralph Benjamin, and Richard Kitney

// Adaptive Smoothing - noise removal with edges preserving.

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

// double radius = 3.0; // default value

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

FIBITMAP* ProcessFilter(FIBITMAP* src_dib, double radius)
{
	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);
	
	
	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, *lined, *linek;

	int size_k = 3;

	int i, j, kin, d;

	BYTE* src_end_row = src_bits + (height-1) * pitch;

	int end_col = width - 1;
		
	// precalculated factor value
	double f = -8 * radius * radius;
	
	// gradient and weights
	double* p_gx = new double[btpp];
	double* p_gy = new double[btpp];
	double* p_weight = new double[btpp];
	double* p_weightTotal = new double[btpp];
	double* p_total = new double[btpp];
		
	for ( int y = 0; y < height; y++ )
	{
		lines = src_bits + y * pitch;
		
		lined = dst_bits + y * pitch;

		for ( int x = 0; x < width; x++ )
		{		
			// original formulas for weight calculation:
			// w(x, y) = exp( -1 * (Gx^2 + Gy^2) / (2 * factor^2) )
			// Gx(x, y) = (I(x + 1, y) - I(x - 1, y)) / 2
			// Gy(x, y) = (I(x, y + 1) - I(x, y - 1)) / 2
			//
			// here is a little bit optimized version

			memset(p_weightTotal, 0, sizeof(double)*btpp);
			memset(p_total, 0, sizeof(double)*btpp);

			// kernel processing
			for (i = 0; i < size_k; i++)
			{
				linek = lines + (i-1) * pitch;

				// 1 aperture y-boundary checkup
				if (linek < src_bits) linek = src_bits;
				if (linek > src_end_row) linek = src_end_row;

				for (j = 0; j < size_k; j++)
				{
					// 2 aperture y-boundary checkup
					if ((linek-pitch) < src_bits) linek = src_bits+pitch;
					if ((linek+pitch) > src_end_row) linek = src_end_row-pitch;

					kin = x+j-1;
					
					// 1 aperture x-boundary checkup
					if (kin < 0) kin = 0;
					if (kin > end_col) kin = end_col;

					// 2 aperture x-boundary checkup
					if ((kin-1) < 0) kin = 1;
					if ((kin+1) > end_col) kin = end_col-1;


					for (d=0; d<btpp; d++)
					{

					p_gx[d] = linek[(kin+1) * btpp + d] - linek[(kin-1) * btpp + d];

					p_gy[d] = (linek + pitch)[kin * btpp + d] - (linek - pitch)[kin * btpp + d];

					p_weight[d] = exp( ( p_gx[d] * p_gx[d] + p_gy[d] * p_gy[d] ) / f );

					p_total[d] += p_weight[d] * linek[kin * btpp + d];

					p_weightTotal[d] += p_weight[d];
					}
				}
			}

			for (d=0; d<btpp; d++)
			{
			if (p_weightTotal[d] == 0.0)

				lined[x * btpp + d] = lines[x * btpp + d];

			else

			// clamp and place result in destination pixel
			lined[x * btpp + d] = (BYTE)MIN(MAX((int)0, (int)(p_total[d] / p_weightTotal[d] + 0.5)), (int)255);
			}
		}
	}

	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));

	delete [] p_gx;
	delete [] p_gy;
	delete [] p_weight;
	delete [] p_weightTotal;
	delete [] p_total;

	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;
}

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

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 != 3) {
		printf("Usage : ad_smooth <input_file> <radius (double)>\n");
		return 0;
	}
	
	
	FIBITMAP *dib = GenericLoader(argv[1], 0);
	
	double radius = atof(argv[2]);
	
	if (dib)
	{		
		// bitmap is successfully loaded!
		
		if (FreeImage_GetImageType(dib) == FIT_BITMAP)
		{
			if (FreeImage_GetBPP(dib) == 8 || FreeImage_GetBPP(dib) == 24)
			{
				FIBITMAP* dst_dib = ProcessFilter(dib, radius);
				
				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;