/*********************************************************************
 * 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.8 $
 *********************************************************************/

#ifndef _EDGEFEATURES_H
#define _EDGEFEATURES_H

#include <math.h>
#include <SortedList.h>
#include <List.h>
#include <Pair.h>
#include <Math.h>
#include <Matrix.h>
#include "../FeatureExtractor.h"
#include "../binary/Labeling.h"
#include "../texture/EdgeDetector.h"

namespace prapi { namespace extras {

	template <class T> class EdgeFeatures : public FeatureExtractor<double, util::Matrix<T> >
	{
	public:
		/**
		 * This are all of the edge features which can be calculated
		 * here. (Descriptin of every one below)
		 *
		 * C0  Amount of edge segments in Matrix.
		 * C1  Average area of edge segment in Matrix.
		 * C2  Average width ration of edge segment to it's height.
		 * C3  Average ration of 
		 * C4  Median of angle Matrix.
		 * C5  Variance of angle Matrix.
		 * C6  Average amount of edge segments in Matrix.
		 * C7  The ration of average amount of pixels in edge segment to length of segment.
		 * C8  The curvature of Matrix calculated with 3*3 mask from gradient matrix.
		 * C9  The curvature of matrix weighted with the value of gradient
		 *     calculated with 3*3 mask from gradient Matrix.
		 * C10 The curvature of Matrix calculated with 3*3 mask from angle matrix.
		 * C11 The curvature of Matrix calculated with 3*3 mask from angle matrix
		 *     and weighted with the value of gradient.
		 * C12 C2 weighted with the mean of gradient Matrix.
		 * C13 C2 weighted with the square of mean in gradient Matrix.
		 * C14 The ration of pixels to area.
		 * C36 The average horizontal distance bethween two different edge chain.
		 * C37 The average vertical distance bethween two different edge chain.
		 **/
		enum Feature { C0,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14};

		/**
		 * The constructor of EdgeFeatures. You have to give gradient matrices to
		 * constructor that it can calculate features. There is no need to threshold
		 * gradient matrices because when you give the thresholded matrix to
		 * getFeature or getFeatureVector they will use Matrix's mask function
		 * for the gradient matrices (if gradient matrices is needeed.
		 *
		 * @param gradients The gradient images.
		 **/
		EdgeFeatures(const util::Matrix<Point<T> >& gradients):
			_matGradients(gradients){}

		/**
		 * Function removes all crossings from the matrix.
		 *
		 * @param mat util::Matrix where the crossings will be removed.
		 **/
		static util::Matrix<T> removeCrossings(const util::Matrix<T>& mat);
		
		/**
		 * Calculates the value of given Feature type in util::Matrix given
		 * in parameter.
		 * NOTE! Edge detection and thresholding must be done before
		 *       using this function. 
		 *
		 *@param mat util::Matrix where wanted to calculate feature.
		 *@param feature Feature which  wanted to calculate.
		 **/
		double getFeature(const util::Matrix<T>& matFeature, Feature feature)
			throw (FeatureExtractionException&);

		/**
		 * Calculates all the features described above and adds them to a
		 * list which will be returned. The order of the features is the
		 * same as in the <i>Feature</i> enumeration.
		 *
		 *@param mat util::Matrix where wanted to calculate features.
		 **/
		util::List<double> getFeatureVector(const util::Matrix<T>& pattern)
			throw (FeatureExtractionException&);

	private:
		/**
		 * The gradient matrices.
		 **/
		util::Matrix<Point<T> > _matGradients;
		util::Matrix<T> _matY;
		/**
		 * Calculates the Width of edge segment which is labeled as a parameter.
		 **/
		int calculateWidth(const util::Matrix<int>& mat, int parameter);
		/**
		 * Calculates the Height of edge segment which is labeled as a parameter.
		 **/
		int calculateHeight(const util::Matrix<int>& mat, int parameter);
		/**
		 * Calculate distance bethween the ends of edge labeled in parameter.
		 **/
		double lengthBethweenEnds(const util::Matrix<int>& mat, int parameter);
	};

	template <class T> util::List<double> EdgeFeatures<T>::getFeatureVector(const util::Matrix<T>& pattern)
		throw (FeatureExtractionException&)
	{
		util::List<double> result(15);
//  		result += getFeature(pattern,C0);
//  		cerr << "C0 made ";
//  		result += getFeature(pattern,C1);
//  		cerr << "C1 made ";
//  		result += getFeature(pattern,C2);
//  		cerr << "C2 made ";
//  		result += getFeature(pattern,C3);
//  		cerr << "C3 made ";
//  		result += getFeature(pattern,C4);
//  		cerr << "C4 made ";
//  		result += getFeature(pattern,C5);
//  		cerr << "C5 made ";
//  		result += getFeature(pattern,C6);
//  		cerr << "C6 made ";
//  		result += getFeature(pattern,C7);
//  		cerr << "C7 made ";
//  		result += getFeature(pattern,C8);
//  		cerr << "C8 made ";
//  		result += getFeature(pattern,C9);
//  		cerr << "C9 made ";
//  		result += getFeature(pattern,C10);
//  		cerr << "C10 made ";
//  		result += getFeature(pattern,C11);
//  		cerr << "C11 made ";
//  		result += getFeature(pattern,C12);
//  		cerr << "C12 made ";
//  		result += getFeature(pattern,C13);
//  		cerr << "C13 made ";
//  		result += getFeature(pattern,C14);
//  		cerr << "C14 made ";
		result.addElements(15,getFeature(pattern,C0),getFeature(pattern,C1),getFeature(pattern,C2),getFeature(pattern,C3),getFeature(pattern,C4),
											 getFeature(pattern,C5),getFeature(pattern,C6),getFeature(pattern,C7),getFeature(pattern,C8),getFeature(pattern,C9),
											 getFeature(pattern,C10),getFeature(pattern,C11),getFeature(pattern,C12),getFeature(pattern,C13),getFeature(pattern,C14));
		return result;
	}

	template <class T> double EdgeFeatures<T>::getFeature(const util::Matrix<T>& matFeature, Feature feature)
		throw (FeatureExtractionException&)
	{
		int matRows = matFeature.getRows();
		int matCols = matFeature.getColumns();
		int matGRows = _matGradients.getRows();
		int matGCols = _matGradients.getColumns();
		double result = 0.0;

		if((matRows != matGRows) && (matCols != matGCols))
			throw FeatureExtractionException("EdgeFeatures<T>::getFeature(const util::Matrix<T>&, Feature): util::Matrix sizes differs.");
		
		switch(feature)
			{
			case C0:
				{
					prapi::binary::Conectivity8<T,int> conectivity;
					result=util::Math::max(conectivity.getTransformedImage(matFeature));
				}
				break;

			case C1:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					for(int i=1;i<max+1;i++)
						result += calculateWidth(matResult,i)*calculateHeight(matResult,i);
					if(max != 0)result /= double(max);
					else result=0;
				}
				break;

			case C2:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					for(int i=1;i<max+1;i++)
						{
							int cHeight=calculateHeight(matResult,i);
							if(cHeight != 0)result += double(calculateWidth(matResult,i))/double(calculateHeight(matResult,i));
						}
					if(max != 0)result /= double(max);
					else result=0;
				}
				break;
					
			case C3:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					for(int i=1;i<max+1;i++)
						{
							int width = calculateWidth(matResult,i);
							int height = calculateHeight(matResult,i);
							if(width>height && height!=0)result=width/height;
							else if(width != 0)result=height/width;
							else result=0;
						}
					if(max != 0)result /= double(max);
					else result=0;
				}
				break;
				
			case C4:
				{
					util::Matrix<T> temp(prapi::texture::EdgeDetector<T>::getAngle(_matGradients));
					Image::mask(temp,matFeature);
					result = util::Math::mean(temp);
				}
				break;

			case C5:
				{
					util::Matrix<T> temp(prapi::texture::EdgeDetector<T>::getAngle(_matGradients));
					Image::mask(temp,matFeature);
					result = util::Math::variance(temp);
				}
				break;

			case C6:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					int sum = 0;
					int* data = matResult.getData();
					for(int i=0;i<matRows*matCols;i++,data++)
						if(*data != 0)sum++;
					if(max != 0)result = double(sum)/double(max);
					else result=0;
					
				}
				break;
					
			case C7:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					for(int i=1;i<max+1;i++)
						{
							double number = lengthBethweenEnds(matResult,i);
							if(number != 0)result += double(util::Math::numberOf(matResult,i))/number;
						}
					if(max != 0)result /= double(max);
					else result=0;
				}
				break;

			case C8:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					util::Matrix<T> matResult(matRows,matCols);
					int pixelSum = 0;
					// mask the matrix
					for(int i=1;i<matRows-1;i++)
						for(int j=1;j<matCols-1;j++)
							{ // check that the pixel is not zero and check that we havent visited there yet.
								if(matRemoved(i,j) != 0)
									{
										if(((matRemoved(i-1,j-1) != 0)&&(matRemoved(i+1,j-1) != 0)) || ((matRemoved(i-1,j-1) != 0)&&(matRemoved(i-1,j+1) != 0)) ||
											 ((matRemoved(i+1,j-1) != 0)&&(matRemoved(i+1,j+1) != 0)) || ((matRemoved(i-1,j+1) != 0)&&(matRemoved(i+1,j+1) != 0)))
											matResult(i,j)=90; // set the right value for the pixel.
										else if(((matRemoved(i-1,j-1) != 0)&&(matRemoved(i+1,j) != 0)) || ((matRemoved(i-1,j+1) != 0)&&(matRemoved(i+1,j) != 0)) ||
														((matRemoved(i-1,j) != 0)&&(matRemoved(i+1,j+1) != 0)) || ((matRemoved(i-1,j) != 0)&&(matRemoved(i+1,j+1) != 0)) ||
														((matRemoved(i-1,j-1) != 0)&&(matRemoved(i,j+1) != 0)) || ((matRemoved(i+1,j+1) != 0)&&(matRemoved(i,j+1) != 0)) ||
														((matRemoved(i-1,j+1) != 0)&&(matRemoved(i,j-1) != 0)) || ((matRemoved(i+1,j+1) != 0)&&(matRemoved(i,j-1) != 0)))
											matResult(i,j)=45; // set the right value for the pixel.
										pixelSum++;//calculate at the same time the amount of edge pixels. 
									}//if
							}//for
					//then calculate the mean of matrix.
					result=util::Math::sum(matResult);
					if(pixelSum == 0) result = 0;
					else result /= pixelSum;
				}
				break;
				
			case C9:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					util::Matrix<T> matResult(matRows,matCols);
					int pixelSum = 0;
					// mask the matrix
					for(int i=1;i<matRows-1;i++)
						for(int j=1;j<matCols-1;j++)
							{ // check that the pixel is not zero and check that we havent visited there yet.
								if(matRemoved(i,j) != 0)
									{
										if(((matRemoved(i-1,j-1) != 0)&&(matRemoved(i+1,j-1) != 0)) || ((matRemoved(i-1,j-1) != 0)&&(matRemoved(i-1,j+1) != 0)) ||
											 ((matRemoved(i+1,j-1) != 0)&&(matRemoved(i+1,j+1) != 0)) || ((matRemoved(i-1,j+1) != 0)&&(matRemoved(i+1,j+1) != 0)))
											matResult(i,j)=90*matRemoved(i,j); // set the right value for the pixel.
										else if(((matRemoved(i-1,j-1) != 0)&&(matRemoved(i+1,j) != 0)) || ((matRemoved(i-1,j+1) != 0)&&(matRemoved(i+1,j) != 0)) ||
														((matRemoved(i-1,j) != 0)&&(matRemoved(i+1,j+1) != 0)) || ((matRemoved(i-1,j) != 0)&&(matRemoved(i+1,j+1) != 0)) ||
														((matRemoved(i-1,j-1) != 0)&&(matRemoved(i,j+1) != 0)) || ((matRemoved(i+1,j+1) != 0)&&(matRemoved(i,j+1) != 0)) ||
														((matRemoved(i-1,j+1) != 0)&&(matRemoved(i,j-1) != 0)) || ((matRemoved(i+1,j+1) != 0)&&(matRemoved(i,j-1) != 0)))
											matResult(i,j)=45*matRemoved(i,j); // set the right value for the pixel.
										pixelSum++;//calculate at the same time the amount of edge pixels. 
									}//if
							}//for
					//then calculate the mean of matrix.
					result=util::Math::sum(matResult);
					if(pixelSum == 0)result=0;
					else result /= pixelSum;
				}
				break;

			case C10:
				{ // calculate angle matrix mask matrix and remove crossings
					util::Matrix<T> matRemoved(prapi::texture::EdgeDetector<T>::getAngle(_matGradients));
					Image::mask(matRemoved,matFeature);// mask the matrix
					// delete the crossings
					matRemoved=EdgeFeatures<T>::removeCrossings(matRemoved);
					util::Matrix<T> matResult(matRows,matCols);
					int pixelSum = 0;
					// mask the matrix

					for(int i=1;i<matRows-1;i++)
						for(int j=1;j<matCols-1;j++)
							{ // check that the pixel is not zero and check that we havent visited there yet.
								if(matRemoved(i,j) != 0)
									{ // calculate the rifht value for pixel.
										if((matRemoved(i-1,j-1) != 0)&&(matRemoved(i+1,j-1) != 0))
											matResult(i,j)=matRemoved(i-1,j-1)-2*matRemoved(i,j)+matRemoved(i+1,j-1);
										else if((matRemoved(i-1,j-1) != 0)&&(matRemoved(i-1,j+1) != 0))
											matResult(i,j)=matRemoved(i-1,j-1)-2*matRemoved(i,j)+matRemoved(i-1,j+1);
										else if((matRemoved(i+1,j-1) != 0)&&(matRemoved(i+1,j+1) != 0))
											matResult(i,j)=matRemoved(i+1,j-1)-2*matRemoved(i,j)+matRemoved(i+1,j+1);
										else if((matRemoved(i-1,j+1) != 0)&&(matRemoved(i+1,j+1) != 0))
											matResult(i,j)=matRemoved(i-1,j+1)-2*matRemoved(i,j)+matRemoved(i+1,j+1);
										//else cerr << "This is the reason why"<< endl;
										pixelSum++;//calculate at the same time the amount of edge pixels. 
									}//if
							}//for
					//then calculate the mean of matrix.
					result=util::Math::sum(matResult);

					if(pixelSum == 0)result=0;
					else result /= pixelSum;
					
				}
				break;
					
			case C11:
				{ // calculate angle matrix mask matrix and remove crossings
					util::Matrix<T> matRemoved(prapi::texture::EdgeDetector<T>::getAngle(_matGradients));
					Image::mask(matRemoved,matFeature);// mask the matrix
					// delete the crossings
					matRemoved=EdgeFeatures<T>::removeCrossings(matRemoved);										
					util::Matrix<T> matResult(matRows,matCols);
					int pixelSum = 0;
					// mask the matrix
					for(int i=1;i<matRows-1;i++)
						for(int j=1;j<matCols-1;j++)
							{ // check that the pixel is not zero and check that we havent visited there yet.
								if(matRemoved(i,j) != 0)
									{ // calculate the rifht value for pixel.
										if((matRemoved(i-1,j-1) != 0)&&(matRemoved(i+1,j-1) != 0))
											matResult(i,j)=(matRemoved(i-1,j-1)-2*matRemoved(i,j)+matRemoved(i+1,j-1))*matRemoved(i,j);
										else if((matRemoved(i-1,j-1) != 0)&&(matRemoved(i-1,j+1) != 0))
											matResult(i,j)=(matRemoved(i-1,j-1)-2*matRemoved(i,j)+matRemoved(i-1,j+1))*matRemoved(i,j);
										else if((matRemoved(i+1,j-1) != 0)&&(matRemoved(i+1,j+1) != 0))
											matResult(i,j)=(matRemoved(i+1,j-1)-2*matRemoved(i,j)+matRemoved(i+1,j+1))*matRemoved(i,j);
										else if((matRemoved(i-1,j+1) != 0)&&(matRemoved(i+1,j+1) != 0))
											matResult(i,j)=(matRemoved(i-1,j+1)-2*matRemoved(i,j)+matRemoved(i+1,j+1))*matRemoved(i,j);
										pixelSum++;//calculate at the same time the amount of edge pixels. 
									}//if
						}//for
					//then calculate the mean of matrix.
					result=util::Math::sum(matResult);
					if(pixelSum == 0)result=0;
					else result /= pixelSum;
				}
				break;
				
			case C12:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					for(int i=1;i<max+1;i++)
						{
							int cHeight=calculateHeight(matResult,i);
							if(cHeight != 0)result = double(calculateWidth(matResult,i))/double(calculateHeight(matResult,i));
						}
					if(max != 0)result /= double(max);
					else result=0;
					result *= util::Math::mean(matFeature);
				}
				break;
					
			case C13:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					for(int i=1;i<max+1;i++)
						{
							int cHeight=calculateHeight(matResult,i);
							if(cHeight != 0)result += double(calculateWidth(matResult,i))/double(calculateHeight(matResult,i));
						}
					if(max != 0)result /= double(max);
					else result=0;
					result *= pow(util::Math::mean(matFeature),2);
				}
				break;

			case C14:
				{
					util::Matrix<T> matRemoved(EdgeFeatures<T>::removeCrossings(matFeature));// remove crossings
					prapi::binary::Conectivity8<T,int> conectivity; // make labeling to matrix
					util::Matrix<int> matResult(conectivity.getTransformedImage(matRemoved));
					int max = util::Math::max(matResult);
					int sum = 0;
					int* data = matResult.getData();
					for(int i=0;i<matRows*matCols;i++,data++)// calculate the pixelsum
						if(*data != 0)sum++;
					if(max != 0)result = double(sum)/double(max);
					else result=0; 
					double temp = 0;
					for(int i=1;i<max+1;i++) // calculate the width*height
						temp += calculateWidth(matResult,i)*calculateHeight(matResult,i);
					if(max != 0)temp/= double(max);
					else result=0;
					if(temp !=0)result /= temp;
					else result=0;
				}
				break;
			}//switch
		
		return result;
	}
	
	template <class T> int EdgeFeatures<T>::calculateWidth(const util::Matrix<int>& mat, int parameter)
	{
		int startPoint = -1,endPoint=0; 
		// go throw the whole matrix
		for(int j=0;j<mat.getColumns();j++)
			for(int i=0;i<mat.getRows();i++)
				{
					if(mat(i,j)==parameter && startPoint == -1)startPoint=j;
					else if(mat(i,j)== parameter)endPoint=j;
				}
		return (endPoint-startPoint+1);
	}

	template <class T> int EdgeFeatures<T>::calculateHeight(const util::Matrix<int>& mat, int parameter)
	{
		int startPoint = -1,endPoint=0; 
		// go throw the whole matrix
		for(int i=0;i<mat.getRows();i++)
			for(int j=0;j<mat.getColumns();j++)
				{
					if(mat(i,j)==parameter && startPoint == -1)startPoint=i;
					else if(mat(i,j)== parameter)endPoint=i;
				}
		return (endPoint-startPoint+1);
	}

	template <class T> double EdgeFeatures<T>::lengthBethweenEnds(const util::Matrix<int>& mat, int parameter)
	{
		util::List<int> points(4);
		// the function just finds two pixels which have only one "friend".
		for(int i=0;i<mat.getRows();i++)
			for(int j=0;j<mat.getColumns();j++)
				{
					if(mat(i,j)==parameter)
						{
							util::List<int> lst(8);
							lst.addElements(8,mat(i-1,j-1),mat(i-1,j),mat(i-1,j+1),mat(i,j-1),mat(i,j+1),mat(i+1,j-1),mat(i+1,j),mat(i+1,j+1));
							if(util::Math::numberOf(lst,parameter) == 1)points.addElements(2,i,j);
						}		
				}
		// calculate the distance bethween points.
		return sqrt(pow((points[0]-points[2]),2)+pow((points[1]-points[3]),2));
	}

	template <class T> util::Matrix<T> EdgeFeatures<T>::removeCrossings(const util::Matrix<T>& mat)
	{
		int matRows = mat.getRows();
		int matCols = mat.getColumns();
		util::Matrix<T> result(matRows,matCols);
		result = -1; // give all pixel value which is -1
		// go throw the whole matrix and try to find edges
		for(int i=1;i<matRows-1;i++)
			for(int j=1;j<matCols-1;j++)
				{ // check that there are values
					if(mat(i,j)!=0 && result(i,j)==-1)
						{
							//inform the pixel that it has been checked
							result(i,j)=mat(i,j);
							util::SortedList<util::Pair<double,util::Pair<int,int> > > pairList(9);
							//add the pair to the list if there are values in pixels
							// first add the center pixel therefor that it could be just
							// the weakest pixel
							util::Pair<int,int> coordinates(i,j);// make coordinates.
							// add the value of gradient and coordinates to new util::Pair.
							util::Pair<double,util::Pair<int,int> > center(mat(i,j),coordinates);
							pairList.addElement(center);// and add them to list.

							// then add all others if there are something
							if(mat(i-1,j-1)!=0)
								{
									util::Pair<int,int> coordinates(i-1,j-1);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > northwest(mat(i-1,j-1),coordinates);
									pairList.addElement(northwest);// and add them to list.
								}
							if(mat(i-1,j)!=0)
								{
									util::Pair<int,int> coordinates(i-1,j);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > north(mat(i-1,j),coordinates);
									pairList.addElement(north);// and add them to list.
								}
							if(mat(i-1,j+1)!=0)
								{
									util::Pair<int,int> coordinates(i-1,j+1);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > northeast(mat(i-1,j+1),coordinates);
									pairList.addElement(northeast);// and add them to list.
								}
							if(mat(i,j-1)!=0)
								{
									util::Pair<int,int> coordinates(i,j-1);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > west(mat(i,j-1),coordinates);
									pairList.addElement(west);// and add them to list.
								}
							if(mat(i,j+1)!=0)
								{
									util::Pair<int,int> coordinates(i,j+1);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > east(mat(i,j+1),coordinates);
									pairList.addElement(east);// and add them to list.
								}
							if(mat(i+1,j-1)!=0)
								{
									util::Pair<int,int> coordinates(i+1,j-1);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > southwest(mat(i+1,j-1),coordinates);
									pairList.addElement(southwest);// and add them to list.
								}
							if(mat(i+1,j)!=0)
								{
									util::Pair<int,int> coordinates(i+1,j);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > south(mat(i+1,j),coordinates);
									pairList.addElement(south);// and add them to list.
								}
							if(mat(i+1,j+1)!=0)
								{
									util::Pair<int,int> coordinates(i+1,j+1);// make coordinates.
									// add the value of gradient and coordinates to new util::Pair.
									util::Pair<double,util::Pair<int,int> > southeast(mat(i+1,j+1),coordinates);
									pairList.addElement(southeast);// and add them to list.
								}

							// now removing pixels which differ the most from median. 
							// pixel will until there are only two pixels left.
							// (note using 8 conectivity)
							int len = pairList.getLength();
							while(pairList.getLength() > 2)
								{
									int last = len-1;
									int median = len >> 1;
									double medianValue=pairList[median].getFirst();
									// check if the first or last diifer the most from median.
									if((medianValue-pairList[0].getFirst()) > (pairList[last].getFirst()-medianValue))
										{ // set the value of pixel to zero to indicate that we have visited
											// here and removed pixel
											result(pairList[0].getSecond().getFirst(),pairList[0].getSecond().getSecond())=0;
											// then remove element from the list.
											pairList.removeElementAt(0);
										}
									// have to remove the last element.
									else
										{
											// set the value of pixel to zero to indicate that we have visited
											// here and removed pixel
											result(pairList[last].getSecond().getFirst(),pairList[last].getSecond().getSecond())=0;
											// then remove element from the list.
											pairList.removeElementAt(last);	
										}
									// and then new length again
									len = pairList.getLength();
								}// while loop							
						}// if mat(i,j)!=0...
					// we have to change the -1 value to zero if the value of gradient was zero.
					else result(i,j)=0;
				}//for loop
		return result;
	}
}}
#endif