/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2001 Topi Mäenpää and Jaakko Viertola
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision: 1.7 $
 *********************************************************************/

#include "../ConvolutionMask.h"
#include "../Image.h"
#include <math.h>
#include <Math.h>
#include <String.h>
#include <AsciiCodec.h>

#include "SmallFilters.h"
#include "JainAndHealeyFeatureExtractor.h"

using namespace util;

namespace prapi { namespace extras {									
	double JHExtractor::calculateFeature(Matrix<double>& mat)
	{
		int r,c;
		float tmp;
		double sum = 0.0;
	
		AllItems(r,c,mat)
			{
				tmp = mat(r,c);
				sum += tmp*tmp;
			}

		return sqrt(sum);
	}

	//Make sure mat1 is always at least as small as mat2
	double JHExtractor::calculateOpponentFeature(Matrix<double>& mat1, double val1, Matrix<double>& mat2, double val2)
	{
		int r,c,offset = (mat2.getRows()-mat1.getRows())>>1;
		double sum = 0.0;
		double divisor = val1*val2;

		if(divisor != 0.0)
			{
				if(offset >= 0)
					AllItems(r,c,mat1)
						sum += mat1(r,c)*mat2(r+offset,c+offset)/divisor;
				else
					AllItems(r,c,mat2)
						sum += mat1(r-offset,c-offset)*mat2(r,c)/divisor;

				return sqrt(2-2*sum);
			}
		else return 0.0;
	}

	Matrix<double> JHExtractor::getFilteredImage(const Matrix<double>& mat, int orientation, int scale) throw (InvalidArgumentException&)
	{
		//cerr << "getFilteredImage\n";
		if (scale > 2 || scale < 0 || orientation > 3 || orientation < 0)
			throw InvalidArgumentException("JHExtractor::getFilteredImage(const Matrix<double>&, int, int): scale or orientation is incorrect");

		//cerr << "Creating convolution mask.\n";
		ConvolutionMask<double> mask(getGaborMask(orientation, scale));
		//cerr << "Mask created.\n";
		return Matrix<double>(Image::convolve(mat,mask));
	}


	//Generates a feature vector containing 120 unichrome and opponent
	//color features computed from the outputs of a Gabor filter bank.
	List<double> JHExtractor::getFeatureVector(const RGBColorImage<>& img) throw (FeatureExtractionException&)
	{
		List<double> result(120);
		result.setLength(120);
		int i,j,k,l;
	
		Matrix<double> channels[3] = {
			Matrix<double>(img.getRed()),
			Matrix<double>(img.getGreen()),
			Matrix<double>(img.getBlue())
		};

		Matrix<double> filteredImages[36]; //Unichrome feature images
		double unichromeFeatures[36];
		double opponentFeatures[84];

		
		ConvolutionMask<double>* masks[12];
		for(int angle=0;angle<4;angle++)
			for(int scale=0;scale<3;scale++)
				masks[3*angle+scale] = new ConvolutionMask<double>(getGaborMask(angle,scale));
				
		for (i=0;i<3;i++)
			for (j=0;j<12;j++)
				filteredImages[i*12+j] = Image::convolve(channels[i],*masks[j]);
		
		
		for (i=0;i<36;i++)
			unichromeFeatures[i] = calculateFeature(filteredImages[i]);

		//The allowed combinations of different scales for the opponent features
		int cc[7][2] = { {0,0}, {1,1}, {2,2}, {0,1}, {1,0}, {1,2}, {2,1} };
	
		for (l=0, i=0;i<3;i++) //Channel combinations
			{
				int c1=0,c2=1;
				switch(i)
					{
					case 1: c1=0; c2=2; break;
					case 2: c1=1; c2=2; break;
					}
				for (j=0;j<4;j++) //Directions
					{
						int index1 = c1*12 + j*3, index2 = c2*12 + j*3;
						for (k=0;k<7;k++, l++)
							{
								int i1 = index1+cc[k][0], i2 = index2+cc[k][1];
								opponentFeatures[l] = calculateOpponentFeature(filteredImages[i1],unichromeFeatures[i1],
																															 filteredImages[i2],unichromeFeatures[i2]);
							}
					}
			}

		for (i=0;i<36;i++)result[i] = unichromeFeatures[i];
		for (i=0;i<84;i++)result[i+36] = opponentFeatures[i];
		
		return result;
	}

	Matrix<double> JHExtractor::getGaborMask(int orientation, int scale)
		throw (InvalidArgumentException&)
	{
		if(orientation < 0 | orientation > 3 | scale < 0 | scale > 2)
			throw InvalidArgumentException("JHExtractor::getGaborMask(int,int): Invalid parameters");
		
		int size=4;
		if(scale==0)size=17;
		else if(scale==1)size=8;
	
		double bt = M_PI/9; // 40 degrees / 2
		double tanBt = tan(bt);
		double um = pow(0.5,double(3-scale))/(1.0+tanBt);

		double sm = sqrt(2.0*log(2.0)) / (2.0 * M_PI * um * tanBt);
		double tn = double(orientation) * M_PI_4; // orientation in radians
		
		double Pi2Um = 2.0*M_PI*um;
		double Sm2 = 2.0*sm*sm;
		double Sm2Pi = 1/(Sm2*M_PI);
		
		// give x and y values at the range -size -> size
		int length=2*size+1;
		List<double> x,y;
		for(int i=0,j=-size;i<length;i++,x+=double(j),y+=double(j),j++);
		
		Matrix<double> filter(length,length);
		double cosTn = cos(tn); // calculate the values ready because they will be used quite
		double sinTn = sin(tn); // many times.
		double x2=0,xCosTn=0;

		for(int i=0;i<length;i++)
			{
				x2 = x[i]*x[i];
				xCosTn = x[i]*cosTn;
				for(int j=0;j<length;j++)
					filter(j,i) = Sm2Pi * exp(-(x2 + y[j]*y[j])/Sm2)*cos(Pi2Um* (xCosTn + y[j]*sinTn));
			}

		return filter;
		         
	}
	
}}