/*
*
* Copyright (C) 2001-2005 Ichiro Fujinaga, Michael Droettboom, and Karl MacMillan
*
* 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.
*/

/*
________________________________________________________________

 bin_ab.c
 $Id: threshold.hpp 962 2007-05-21 13:29:29Z cdalitz $
 Copyright 1990, Blab, UiO
 Image processing lab, Department of Informatics
 University of Oslo
 E-mail: blab@ifi.uio.no
 ________________________________________________________________
 
  Permission to use, copy, modify and distribute this software and its
  documentation for any purpose and without fee is hereby granted, 
  provided that this copyright notice appear in all copies and that 
  both that copyright notice and this permission notice appear in supporting
  documentation and that the name of B-lab, Department of Informatics or
  University of Oslo not be used in advertising or publicity pertaining 
  to distribution of the software without specific, written prior permission.
  
  B-LAB DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL
  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL B-LAB
  BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION
  OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN 
  CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
   
  References:
  &Ahmed S. Abutaleb
  "Automatic thresholding of gray-level pictures using
  two-dimensional entropy",
  Computer Vision, Graphics, and Image Processing,
  vol 47, pp 22-32, 1989.
	
  Original author:
  Øivind Due Trier
*/

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

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

#include <limits>

using namespace std;

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

inline void SetPixel(BYTE *bits, unsigned x, BYTE* value)
{   // this function is simplified from FreeImage_SetPixelIndex
	
	*value ? bits[x >> 3] |= (0x80 >> (x & 0x7)) : bits[x >> 3] &= (0xFF7F >> (x & 0x7));
}

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

void mean_filter(FIBITMAP* src_dib, int region_size, double* p_mean)
{
	unsigned width = FreeImage_GetWidth(src_dib);
	
	unsigned height = FreeImage_GetHeight(src_dib);
	
	unsigned pitch = FreeImage_GetPitch(src_dib);
	
	unsigned bpp = FreeImage_GetBPP(src_dib);
	
	unsigned x, y, i, j;
	
	int k_index;
	
	BYTE* src_bits = (BYTE*)FreeImage_GetBits(src_dib); // The image raster
	
	BYTE* src_end_row = src_bits + (height-1) * pitch;
	
	int end_col = width - 1;
	
	BYTE* lines, *linek;
	
	int sum;
	
	unsigned area = region_size * region_size;
	
	
	for (y = 0; y < height; y++)
	{
		lines = src_bits + y * pitch;
		
		for (x = 0; x < width; x++)
		{
			sum = 0;
			
			// kernel processing
			for (i = 0; i < region_size; i++)
				for (j = 0; j < region_size; j++)
				{					
					linek = lines + (i-1) * pitch;
					
					if (linek < src_bits) linek = src_bits;
					if (linek > src_end_row) linek = src_end_row;
					
					k_index = x+j-1;
					
					if (k_index < 0) k_index = 0;
					if (k_index > end_col) k_index = end_col;
					
					sum += linek[k_index];
				}
				
				p_mean[y*width+x] = (double)sum / area;
		}
	}
}

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

FIBITMAP* ProcessFilter(FIBITMAP* src_dib, int region_size)
{
	unsigned width = FreeImage_GetWidth(src_dib);
	
	unsigned height = FreeImage_GetHeight(src_dib);
	
	unsigned src_pitch = FreeImage_GetPitch(src_dib);
	
	unsigned bpp = FreeImage_GetBPP(src_dib);
	
	FIBITMAP* dst_dib = FreeImage_Allocate(width, height, 1);
	
	// Build a monochrome palette
	RGBQUAD *pal = FreeImage_GetPalette(dst_dib);
	pal[0].rgbRed = pal[0].rgbGreen = pal[0].rgbBlue = 0;
	pal[1].rgbRed = pal[1].rgbGreen = pal[1].rgbBlue = 255;
	
	unsigned dst_pitch = FreeImage_GetPitch(dst_dib);
	
	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;
	
	unsigned x, y;
	
	BYTE val;
	
	int int_size = height*width*sizeof(double);
	
	double* p_average = (double*)malloc(int_size);	
	
	memset(p_average,0,int_size);
	
	mean_filter(src_dib, region_size, p_average);

	int hw = 256; // square histogram width
	
	int hist_size = hw*hw*sizeof(double);	
	
	double* histogram = (double*)malloc(hist_size);
	double* P_histogram = (double*)malloc(hist_size);
	double* H_histogram = (double*)malloc(hist_size);
	
	memset(histogram,0,hist_size);
	memset(P_histogram,0,hist_size);
	memset(H_histogram,0,hist_size);
	
	unsigned a, b, s, t;
	
	for (y = 0; y < height; y++)
	{
		lines = src_bits + y * src_pitch;
		
		for (x = 0; x < width; x++)
		{
			a = lines[x];
			b = p_average[x*width+y];
			histogram[a*hw+b] = histogram[a*hw+b] + 1.0;
		}
	}	
	
	double one_over_area = 1.0 / (width * height);
	
	double P_sum, p, H_sum;
	
	for (b = 0; b < 256; b++)			
		for (a = 0; a < 256; a++)
			
			histogram[a*hw+b] = histogram[a*hw+b] * one_over_area;
		
		P_sum = 0.0;
		
		for (s = 0; s < 256; s++)
		{
			P_sum += histogram[s*hw];
			
			P_histogram[s*hw] = P_sum;
		}
		
		for (t = 1; t < 256; t++)
		{
			P_sum = 0.0;
			
			for (s = 0; s < 256; s++)
			{
				P_sum += histogram[s*256+t];
				
				P_histogram[s*hw+t] = P_histogram[s*hw+(t - 1)] + P_sum;
			}
		}
		
		H_sum = 0.0;
		
		for (s = 0; s < 256; s++)
		{
			p = histogram[s*hw];
			
			if (p != 0)
				H_sum -= p * log(p);
			
			H_histogram[s*hw] = H_sum;
		}
		
		for (t = 1; t < 256; t++)
		{
			H_sum = 0.0;
			
			for (s = 0; s < 256; ++s)
			{
				p = histogram[s*hw+t];
				
				if (p != 0)
					H_sum -= p * log(p);
				
				H_histogram[s*hw+t] = H_histogram[s*hw+(t - 1)] + H_sum;
			}
		}
		
		double Phi_max = std::numeric_limits<double>::min();
		double tiny = 1e-6;
		
		double H_end = H_histogram[255*hw+255];
		unsigned threshold = 0, avg_threshold = 0;
		double Phi, P, H;
		
		for (s = 0; s < 256; s++)
			for (t = 0; t < 256; t++)
			{
				P = P_histogram[s*hw+t];
				H = H_histogram[s*hw+t];
				
				if ((P > tiny) && ((1.0 - P) > tiny))
				{	
					Phi = log(P * (1.0 - P)) + H / P + (H_end - H) / (1.0 - P);
					
					if (Phi > Phi_max)
					{
						Phi_max = Phi;
						threshold = s;
						avg_threshold = t;
					}
				}
			}
			
			for (y = 0; y < height; y++)
			{
				lines = src_bits + y * src_pitch;
				
				lined = dst_bits + y * dst_pitch;
				
				for (x = 0; x < width; x++)
				{
					if (lines[x] <= threshold && p_average[x*width+y] <= avg_threshold)			
						
						val = 0;
					
					else
						
						val = 255;
					
					SetPixel(lined, x, &val);
				}
			}
			
			// Copying the DPI...
			
			FreeImage_SetDotsPerMeterX(dst_dib, FreeImage_GetDotsPerMeterX(src_dib));
			
			FreeImage_SetDotsPerMeterY(dst_dib, FreeImage_GetDotsPerMeterY(src_dib));
			
			free(p_average);	
			
			free(histogram);
			free(P_histogram);
			free(H_histogram);
			
			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 : abutaleb_thres <input_file> <region_size (int)>\n");
		return 0;
	}
	
	FIBITMAP *dib = GenericLoader(argv[1], 0);
	
	int region_size = atoi(argv[2]);
	
	if (dib)
	{		
		// bitmap is successfully loaded!
		
		if (FreeImage_GetImageType(dib) == FIT_BITMAP)
		{
			if (FreeImage_GetBPP(dib) == 8)
			{
				FIBITMAP* dst_dib = ProcessFilter(dib, region_size);
				
				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;
}