// ==========================================================
// 4th order (cubic) b-spline smoothing
//
// Design and implementation by
// - Hervй Drolon (drolon@infonie.fr)
// - Detlev Vendt (detlev.vendt@brillit.de)
//
// This file is part of FreeImage 3
//
// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
// THIS DISCLAIMER.
//
// Use at your own risk!
// ==========================================================

// This algorithm was taken from the FreeImage sourcecodes
// and simplified for a non-rescaling.
//
//	Copyright (C) 2007-2008:
//	monday2000	monday2000@yandex.ru

// Main reference:
// Paul Heckbert, C code to zoom raster images up or down, with nice filtering. 
// UC Berkeley, August 1989. [online] http://www-2.cs.cmu.edu/afs/cs.cmu.edu/Web/People/ph/heckbert.html
// 
// Heckbert references:
// 
// Oppenheim A.V., Schafer R.W., Digital Signal Processing, Prentice-Hall, 1975
// Hamming R.W., Digital Filters, Prentice-Hall, Englewood Cliffs, NJ, 1983
// Pratt W.K., Digital Image Processing, John Wiley and Sons, 1978
// Hou H.S., Andrews H.C., "Cubic Splines for Image Interpolation and Digital Filtering", 
// IEEE Trans. Acoustics, Speech, and Signal Proc., vol. ASSP-26, no. 6, pp. 508-517, Dec. 1978.
// 

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

/////////////////////////////////////////////////////////////////////////////////////////////
//
// 4th order (cubic) b-spline
//
double CBSplineFilter(double dVal)
{ 
//  Default fixed width = 2

	dVal = fabs(dVal);
	if(dVal < 1) return (4 + dVal*dVal*(-6 + 3*dVal)) / 6;
	if(dVal < 2) {
		double t = 2 - dVal;
		return (t*t*t / 6);
	}
	return 0;
}

/////////////////////////////////////////////////////////////////////////////////////////////
//
//  Filter weights table.
//  This class stores contribution information for an entire line (row or column).
//
class CWeightsTable
{
// 
//  Sampled filter weight table.
//  Contribution information for a single pixel
//
typedef struct
 {
	/// Normalized weights of neighboring pixels
	double *Weights;
	/// Bounds of source pixels window
	int Left, Right;   
} Contribution;  

private:
	/// Row (or column) of contribution weights 
	Contribution *m_WeightTable;
	/// Filter window size (of affecting source pixels) 
	DWORD m_WindowSize;
	/// Length of line (no. of rows / cols) 
	DWORD m_LineLength;

public:
	/** 
	Constructor
	Allocate and compute the weights table
	@param uSize Length (in pixels) of the line buffer
	*/
	CWeightsTable(DWORD uSize);

	/**
	Destructor
	Destroy the weights table
	*/
	~CWeightsTable();

	/** Retrieve a filter weight, given source and destination positions
	@param dst_pos Pixel position in destination line buffer
	@param src_pos Pixel position in source line buffer
	@return Returns the filter weight
	*/
	double getWeight(int dst_pos, int src_pos) {
		return m_WeightTable[dst_pos].Weights[src_pos];
	}

	/** Retrieve left boundary of source line buffer
	@param dst_pos Pixel position in destination line buffer
	@return Returns the left boundary of source line buffer
	*/
	int getLeftBoundary(int dst_pos) {
		return m_WeightTable[dst_pos].Left;
	}

	/** Retrieve right boundary of source line buffer
	@param dst_pos Pixel position in destination line buffer
	@return Returns the right boundary of source line buffer
	*/
	int getRightBoundary(int dst_pos) {
		return m_WeightTable[dst_pos].Right;
	}
};

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

CWeightsTable::CWeightsTable(DWORD uSize)
{
	DWORD u;
	double dWidth = 2; //  Default fixed width = 2 - for 4th order (cubic) b-spline
	
	// allocate a new line contributions structure
	//
	// window size is the number of sampled pixels
	m_WindowSize = 2 * (int)ceil(dWidth) + 1; 
	m_LineLength = uSize; 
	// allocate list of contributions 
	m_WeightTable = (Contribution*)malloc(m_LineLength * sizeof(Contribution));
	for(u = 0 ; u < m_LineLength ; u++) {
		// allocate contributions for every pixel
		m_WeightTable[u].Weights = (double*)malloc(m_WindowSize * sizeof(double));
	}
	
	
	for(u = 0; u < m_LineLength; u++) {
		// scan through line of contributions
		double dCenter = (double)u;   // reverse mapping
		// find the significant edge points that affect the pixel
		int iLeft = MAX (0, (int)floor (dCenter - dWidth)); 
		int iRight = MIN ((int)ceil (dCenter + dWidth), int(uSize) - 1); 
		
		// cut edge points to fit in filter window in case of spill-off
		if((iRight - iLeft + 1) > int(m_WindowSize)) {
			if(iLeft < (int(uSize) - 1 / 2)) {
				iLeft++; 
			} else {
				iRight--; 
			}
		}
		
		m_WeightTable[u].Left = iLeft; 
		m_WeightTable[u].Right = iRight;
		
		int iSrc = 0;
		double dTotalWeight = 0;  // zero sum of weights
		for(iSrc = iLeft; iSrc <= iRight; iSrc++) {
			// calculate weights
			double weight = CBSplineFilter(dCenter - (double)iSrc);
			m_WeightTable[u].Weights[iSrc-iLeft] = weight;
			dTotalWeight += weight;
		}
		if((dTotalWeight > 0) && (dTotalWeight != 1)) {
			// normalize weight of neighbouring points
			for(iSrc = iLeft; iSrc <= iRight; iSrc++) {
				// normalize point
				m_WeightTable[u].Weights[iSrc-iLeft] /= dTotalWeight; 
			}
			// simplify the filter, discarding null weights at the right
			iSrc = iRight - iLeft;
			while(m_WeightTable[u].Weights[iSrc] == 0){
				m_WeightTable[u].Right--;
				iSrc--;
				if(m_WeightTable[u].Right == m_WeightTable[u].Left)
					break;
			}			
		}
	} 
}

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

CWeightsTable::~CWeightsTable() {
	for(DWORD u = 0; u < m_LineLength; u++) {
		// free contributions for every pixel
		free(m_WeightTable[u].Weights);
	}
	// free list of pixels contributions
	free(m_WeightTable);
}

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

/// Performs horizontal image filtering
void horizontalFilter(FIBITMAP *src, unsigned width, unsigned height, FIBITMAP *dst)
{	
	// allocate and calculate the contributions
	CWeightsTable weightsTable(width);	
	
	BYTE* src_bits = (BYTE*)FreeImage_GetBits(src); // The image raster
	
	BYTE* dst_bits = (BYTE*)FreeImage_GetBits(dst); // The image raster
	
	BYTE* lines, *lined;

	unsigned pitch = FreeImage_GetPitch(src);
		
	// Calculate the number of bytes per pixel (1 for 8-bit, 3 for 24-bit)
	unsigned btpp = FreeImage_GetBPP(src) / 8;

	// step through rows
	for(unsigned y = 0; y < height; y++)
	{		
		lined = dst_bits + y * pitch;

		lines = src_bits + y * pitch;		

		// process each row
		for(unsigned x = 0; x < width; x++)
		{
			// loop through row
			double value[3] = {0, 0, 0};					// 3 = 24bpp max
			int iLeft = weightsTable.getLeftBoundary(x);    // retrieve left boundary
			int iRight = weightsTable.getRightBoundary(x);  // retrieve right boundary
			
			for(int i = iLeft; i <= iRight; i++)
			{
				// scan between boundaries
				// accumulate weighted effect of each neighboring pixel
				double weight = weightsTable.getWeight(x, i-iLeft);
								
				for (unsigned j = 0; j < btpp; j++)
				
					value[j] += (weight * (double)lines[i * btpp + j]);				
			} 
			
			// clamp and place result in destination pixel
			for (unsigned j = 0; j < btpp; j++)
			
				lined[x * btpp + j] = (BYTE)MIN(MAX((int)0, (int)(value[j] + 0.5)), (int)255);						
		}
	}
}

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

/// Performs vertical image filtering
void verticalFilter(FIBITMAP *src, unsigned width, unsigned height, FIBITMAP *dst)
{	
	// allocate and calculate the contributions
	CWeightsTable weightsTable(height);	

	BYTE* src_bits = (BYTE*)FreeImage_GetBits(src); // The image raster
	
	BYTE* dst_bits = (BYTE*)FreeImage_GetBits(dst); // The image raster
	
	BYTE* lines, *lined;

	unsigned pitch = FreeImage_GetPitch(src);	
	
	// Calculate the number of bytes per pixel (1 for 8-bit, 3 for 24-bit or 4 for 32-bit)
	unsigned btpp = FreeImage_GetBPP(src) / 8;
	
	// step through columns	
	for(unsigned x = 0; x < width; x++)
	{
		// process each column
		for(unsigned y = 0; y < height; y++)
		{
			// loop through column
			double value[3] = {0, 0, 0};					// 3 = 24bpp max
			int iLeft = weightsTable.getLeftBoundary(y);    // retrieve left boundary
			int iRight = weightsTable.getRightBoundary(y);  // retrieve right boundary
			
			lined = dst_bits + y * pitch;			
			
			for(int i = iLeft; i <= iRight; i++)
			{
				// scan between boundaries
				// accumulate weighted effect of each neighboring pixel
				double weight = weightsTable.getWeight(y, i-iLeft);

				lines = src_bits + i * pitch;
				
				for (unsigned j = 0; j < btpp; j++)

					value[j] += (weight * (double)lines[x * btpp + j]);				
			}
			
			// clamp and place result in destination pixel
			for (unsigned j = 0; j < btpp; j++)

				lined[x * btpp + j] = (BYTE)MIN(MAX((int)0, (int)(value[j] + 0.5)), (int)255);
		}		
	} 
}	

/////////////////////////////////////////////////////////////////////////////////////////////
//
// This function performs B-Spline smoothing. It works with 8- and 24- buffers.
//
//  References : 
//  [1] Paul Heckbert, C code to zoom raster images up or down, with nice filtering. 
//  UC Berkeley, August 1989. [online] http://www-2.cs.cmu.edu/afs/cs.cmu.edu/Web/People/ph/heckbert.html
//  [2] Eran Yariv, Two Pass Scaling using Filters. The Code Project, December 1999. 
//  [online] http://www.codeproject.com/bitmap/2_pass_scaling.asp
//

FIBITMAP* ProcessFilter(FIBITMAP *src)
{ 
	DWORD width  = FreeImage_GetWidth(src); 
	DWORD height = FreeImage_GetHeight(src);
	
	unsigned bpp = FreeImage_GetBPP(src);
	
	FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(src);
	
	// allocate the dst image
	FIBITMAP *dst = FreeImage_Allocate(width, height, bpp);
	if(!dst) return NULL;
	
	if(bpp == 8)
	{
		if(FreeImage_GetColorType(src) == FIC_MINISWHITE) {
			// build an inverted greyscale palette
			RGBQUAD *dst_pal = FreeImage_GetPalette(dst);
			for(int i = 0; i < 256; i++) {
				dst_pal[i].rgbRed = dst_pal[i].rgbGreen =
					dst_pal[i].rgbBlue = (BYTE)(255 - i);
			}
		} else {
			// build a greyscale palette
			RGBQUAD *dst_pal = FreeImage_GetPalette(dst);
			for(int i = 0; i < 256; i++) {
				dst_pal[i].rgbRed = dst_pal[i].rgbGreen =
					dst_pal[i].rgbBlue = (BYTE)i;
			}
		}
	}
	
	// allocate a temporary image
	FIBITMAP *tmp = FreeImage_Allocate(width, height, bpp);
	if(!tmp) {
		FreeImage_Unload(dst);
		return NULL;
	}
	
	// smooth source image horizontally into temporary image
	horizontalFilter(src, width, height, tmp);
	
	// smooth temporary image vertically into result image    
	verticalFilter(tmp, width, height, dst);
	
	// free temporary image
	FreeImage_Unload(tmp);

	// Copying the DPI...
	
	FreeImage_SetDotsPerMeterX(dst, FreeImage_GetDotsPerMeterX(src));
	
	FreeImage_SetDotsPerMeterY(dst, FreeImage_GetDotsPerMeterY(src));
	
	return dst;
} 

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

/**
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 != 2) {
		printf("Usage : bspline <input_file>\n");
		return 0;
	}
	
	FIBITMAP *dib = GenericLoader(argv[1], 0);
	
	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);
				
				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;
}