/* 
This file is part of the TRI2 source code.

Copyright (c) 2010, 2011, Gray Institute University of Oxford

    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 3 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, see <http://www.gnu.org/licenses/>.
 */


/*********************************************************************
 * Aortic Ring Processing
 *
 * Processes an image to remove a large bright central aortic ring
 * and the binarise the remaining vessels that surround it.
 * This can then be used to mask the input image for further processing.
 *
 * P Barber
 * 2011
 *********************************************************************/


#include "main.h"
#include "windowMenu.h"

#include <utility.h>
#include <userint.h>
#include "imageProcessing_ui.h"

#include "freeimagealgorithms.h"
#include "freeimagealgorithms_utilities.h"
#include "freeimagealgorithms_filters.h"
#include "freeimagealgorithms_arithmetic.h"
#include "FreeImageAlgorithms_Convolution.h"
#include "FreeImageAlgorithms_Morphology.h"
#include "FreeImageAlgorithms_Particle.h"

/*********************************************************************
 * RemoveAorticRing
 * Takes in an image and returns the processed binary image of the
 * vessels around the aortic ring.
 * FIBITMAP *imin - the input image
 * returns FIBITMAP* - the output image, created in this function
 *********************************************************************/
FIBITMAP* RemoveAorticRing(FIBITMAP *imin)
{
	FIBITMAP *temp1=NULL;
	FIBITMAP *temp2=NULL;
	FIBITMAP *temp3=NULL;
	int i, n=10;
	const double vals[] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
	FilterKernel kernel;
	FIABITMAP *fiab=NULL;

	temp1 = FreeImage_Clone(imin);
	temp2 = FIA_Blur(imin, FIA_KERNEL_GAUSSIAN, 7);

	if (temp1==NULL || temp2==NULL)	return NULL;

	FIA_InPlaceConvertToGreyscaleFloatType(&temp1, FIT_FLOAT);
	if (temp1==NULL || temp2==NULL)	return NULL;
	FIA_SubtractGreyLevelImages(temp1, temp2);							FreeImage_Unload(temp2);
	if (temp1==NULL || temp2==NULL)	return NULL;
	FIA_InPlaceConvertToStandardType(&temp1, 0);
	if (temp1==NULL || temp2==NULL)	return NULL;
	FIA_InPlaceThreshold(temp1, 8.0, 255.0, 255.0);
	if (temp1==NULL || temp2==NULL)	return NULL;

//	DisplayImageIPI(1, temp1, "Thresholded", 0, 0);

	// Closing
	kernel = FIA_NewKernel(2, 2, vals, 25.0);
	fiab = FIA_SetBorder (temp1, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); // FreeImage_Unload(temp1);  // keep temp1
	temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryErosion(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryErosion(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryErosion(fiab, kernel); FIA_Unload(fiab);

	if (temp1==NULL || temp2==NULL)	return NULL;

	temp3 = FIA_Fillholes(temp2, 1);
	FIA_InPlaceConvertToGreyscaleFloatType(&temp3, FIT_FLOAT);
	FIA_SubtractGreyLevelImages(temp3, temp2);							FreeImage_Unload(temp2);
	FIA_InPlaceConvertToStandardType(&temp3, 0);

	temp2 = FIA_Fillholes(temp3, 1);	FreeImage_Unload(temp3);

	for (i=0, n=10; i<n; i++){
		fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
		temp2 = FIA_BinaryErosion(fiab, kernel); FIA_Unload(fiab);
	}
	for (i=0; i<n; i++){
		fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
		temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);
	}

	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);
	fiab = FIA_SetBorder (temp2, kernel.x_radius, kernel.y_radius, BorderType_Constant, 0.0); FreeImage_Unload(temp2);
	temp2 = FIA_BinaryDilation(fiab, kernel); FIA_Unload(fiab);

	FIA_InPlaceConvertToGreyscaleFloatType(&temp2, FIT_FLOAT);
	FIA_MultiplyGreyLevelImageConstant(temp2, 255.0);
	FIA_InPlaceConvertToStandardType(&temp2, 0);
//	DisplayImageIPI(1, temp2, "ring", 0, 0);

	// Opening
	temp3 = FIA_Binary3x3Opening(temp1);

	FIA_InPlaceConvertToGreyscaleFloatType(&temp3, FIT_FLOAT);
	FIA_MultiplyGreyLevelImageConstant(temp3, 255.0);
	FIA_InPlaceConvertToStandardType(&temp3, 0);
//	DisplayImageIPI(1, temp3, "opened - vessels", 0, 0);

	FIA_InPlaceConvertToGreyscaleFloatType(&temp3, FIT_FLOAT);
	FIA_SubtractGreyLevelImages(temp3, temp2);			FreeImage_Unload(temp2);
	FIA_InPlaceConvertToStandardType(&temp3, 0);

	FIA_InPlaceConvertToGreyscaleFloatType(&temp3, FIT_FLOAT);
	FIA_MultiplyGreyLevelImageConstant(temp3, 255.0);
	FIA_InPlaceConvertToStandardType(&temp3, 0);

	return temp3;
}

/*********************************************************************
 * doAorticRing
 * Called from the main program, performs the processing and creates 
 * a new image space to display the result in.
 * int wSpace - the TRI2 workspace of the image to process
 * int iSpace - the TRI2 imagespace of the image to process
 * returns int - 0 = no error, -ve = error occurred
 *********************************************************************/
int doAorticRing (int wSpace, int iSpace)
{	// Called from the main program, performs the processing and creates a new image space to display the result in
	int input, ispace, newISpace;
	FIBITMAP *image=NULL, *dest_image=NULL;
	int progress = CreateProgressDialog("Aortic Ring Processing", "Percent complete", 1, VAL_NO_INNER_MARKERS, "");

	input = iSpace;

	dest_image = RemoveAorticRing(gWS[wSpace].gIS[input].fib_image);

	if (dest_image==NULL) 
	{
		DiscardProgressDialog(progress);
		return -1;
	}

	UpdateProgressDialog(progress, 80, 0);
	
	ispace = createISpaceInWSpaceWithColTable (&wSpace, &newISpace, gWS[wSpace].gIS[input].colTable, gWS[wSpace].gIS[input].nColours);
	if (ispace<0) return(-1);

	gWS[wSpace].gIS[newISpace].fib_image = dest_image;
	sprintf(gWS[wSpace].gIS[newISpace].title, "%s - AorticAssay", gWS[wSpace].gIS[input].title);
	
	initISpaceDisplay(wSpace, dest_image, newISpace, NULL, -1, -1);
	displayISpace(wSpace, newISpace);

	DiscardProgressDialog(progress);

	return 0;
}