/* ImLib3D
 * Copyright (c) 2001, ULP-IPB Strasbourg.
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
/*! \file MorphologicalOperators.hxx */
#ifndef _MorphologicalOperators_hxx
#define _MorphologicalOperators_hxx

#include<ImLib3D/Image3D.hpp>
#include<ImLib3D/MorphologicalOperators.hpp>
#include<algorithm>
#include<ImLib3D/ImageProgress.hpp>

template <class ImageType>
void
IP3D::MedianFilter(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
	ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
	args.SameSize();
	ImageType &res=args.GetDest();

	// create vector for containing neighborValues
	vector<typename ImageType::value_type> neighborValues;
	neighborValues.assign(mask.size(),src(0,0,0));
	int imedian=mask.size()/2;

	typename ImageType::iteratorXYZ rp;
	ImageProgress tracker("MedianFilter",rp);
	for(rp=res.begin();rp!=res.end();rp++)
	{
		// find neighborValues
		neighborValues[0]=src.SafeValue(rp.Pos()+mask.GetPoint(0));
		for(uint i=1;i<mask.size();i++)
		{
			neighborValues[i]=src.SafeValue(rp.Pos()+mask[i]);
		}
		// now find median of neighborValues
		sort(neighborValues.begin(),neighborValues.end());
		*rp=neighborValues[imedian];
	}	

	//args.Finished();
}


template <class ImageType>
void
IP3D::SharpeningFilter(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
	typedef typename ImageType::value_type Im3DValue;
	ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
	args.SameSize();
	ImageType &res=args.GetDest();

	// create vector for containing neighborValues
	vector<Im3DValue> neighborValues;
	neighborValues.assign(mask.size(),src(0,0,0));

	typename ImageType::iteratorXYZ rp;
	ImageProgress tracker("SharpeningFilter",rp);
	for(rp=res.begin();rp!=res.end();rp++)
	{
		*rp=src(rp.Pos());
		// find neighborValues
		neighborValues[0]=src.SafeValue(rp.Pos()+mask.GetPoint(0));
		for(uint i=1;i<mask.size();i++)
		{
			neighborValues[i]=src.SafeValue(rp.Pos()+mask[i]);
		}
		Im3DValue min = neighborValues[0];
		Im3DValue max = neighborValues[0];
		int minindex=0;
		int maxindex=0; 
		// find min and max neighbor-value
		for(uint i=1;i<mask.size();i++)
		{
			if (neighborValues[i] > max)
			{max = neighborValues[i];maxindex=i;}
			if (neighborValues[i] < min)
			{min = neighborValues[i];minindex=i;}
		}
		// now assign current to min or max, whichever value is closest:
		if ( ((*rp)-min) < (max-(*rp)))
			*rp = min;
		else if ( ((*rp)-min) > (max-(*rp)))
			*rp = max;
		// or leave unchanged when distance the same
	}	
	//args.Finished();
}


template <class ImageType>
void
IP3D::Erosion(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
	ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
	args.SameSize();
	ImageType &res=args.GetDest();

	ImageNeighbors<const ImageType> neighbors(src,mask);
	uint npts=mask.size();
	typename ImageType::iteratorXYZ rp;
	ImageProgress tracker("Erosion",rp);
	for(rp=res.begin();rp!=res.end();rp++)
	{
		const Vect3Di &pos=rp.Pos();
		if(neighbors.IsSafe(pos))
		{
			neighbors.MoveTo(pos);
			*rp=neighbors[0];
			for(uint i=1;i<npts;i++)
			{
				*rp=min(*rp,neighbors[i]);
			}
		}
		else
		{
			*rp=src.SafeValue(rp.Pos()+mask.GetPoint(0));
			for(uint i=1;i<npts;i++)
			{
				*rp=min(*rp,src.SafeValue(rp.Pos()+mask[i]));
			}
		}
	}	
}

template <class ImageType>
void
IP3D::Dilation(const ImageType& src, const StructureElement &mask, ImageType& res0)
{
	ImageProcessorArgs<ImageType> args(ImArgs(src),res0);
	args.SameSize();
	ImageType &res=args.GetDest(0);
	typename ImageType::iteratorXYZ rp;
	ImageProgress tracker("Dilation",rp);
	for(rp=res.begin();rp!=res.end();rp++)
	{
		*rp=src.SafeValue(rp.Pos()+mask.GetPoint(0));
		for(uint i=1;i<mask.size();i++)
		{
			*rp=max(*rp,src.SafeValue(rp.Pos()+mask[i]));
		}
	}	
}


template <class ImageType>
void
IP3D::MakeValueList(const ImageType& src, ImageValues<typename ImageType::value_type>& res)
{
	typename ImageType::const_iteratorXYZ sp;
	for(sp=src.begin();sp!=src.end();sp++)
	{
		res.values[*sp].point=sp.Pos();
		res.values[*sp].size++;
	}	
}


template <class ImageType>
void
IP3D::ConnectedComponentLabelling(const ImageType& src, const StructureElement &mask0, LabelImage3D& res,typename ImageType::value_type *background0)
{
	// ALGO: For each pixel 
	typedef typename ImageType::value_type Im3DValue;
	pair<bool,Im3DValue> backgroundInfo=DefaultBackground<ImageType>();
	bool hasBackground=backgroundInfo.first;
	Im3DValue background=backgroundInfo.second;
	vector<int> indx2region;
	indx2region.reserve(1000);
	// default indx=0 -> region=0 = background
	indx2region.push_back(0);
	int nbRegions=1;

	res.Resize(src.Size());
	res.Fill(0);
//  	for(LabelImage3D::iteratorFast fp=res.begin();fp!=res.end();fp++){if(*fp){*fp=1;}}

	StructureElement mask(mask0,MORPHO_FORWARD);
	LabelImage3D::iteratorXYZ rp;
	vector<int> neigborIndexes;
	neigborIndexes.reserve(mask.size());

	ImageProgress tracker("ConnectedComponentLabelling",rp);
	for(rp=res.begin();rp!=res.end();rp++)
	{
		int talk=0;
		if(rp.Pos()==Vect3Di(48,5,76)){talk=1;}
		const Im3DValue &srcval=src(rp.Pos());
		if(hasBackground && background==srcval){*rp=0;continue;}
		// let's see if this voxel has any neigbors
		neigborIndexes.clear();
		for(uint i=0;i<mask.size();i++)
		{
//			if((currIndex=res.SafeValue(rp.Pos()+mask[i])))
			Vect3Di npos=rp.Pos()+mask[i];
			if(src.IsInside(npos) && src(npos)==srcval)
			{// found a neigbor!
//  					mprintf("found neigbor index:%2d ",currIndex);
//  					(rp.Pos()+mask[i]).Show("pos:","\n");
				int neighborIndex=res(npos);
				neigborIndexes.push_back(neighborIndex);
			}
		}
		if(neigborIndexes.size())
		{// found neigbors!
			if(talk)
			{
				for(uint i=0;i<neigborIndexes.size();i++)
				{
					mprintf("found neigborIndexes:%d:%d\n",i,neigborIndexes[i]);
				}
			}
			// set this pixel to the same value as one of those neigbors
			int anIndexOfThisRegion=neigborIndexes[0];
			*rp=anIndexOfThisRegion;
			// make all neigbor indexes point to the same region
			int thisRegion=indx2region[anIndexOfThisRegion];
			if(talk)
			{
				mprintf("anIndexOfThisRegion:%d\n",anIndexOfThisRegion);
				mprintf("thisRegion:%d\n",thisRegion);
			}
			for(uint i=1;i<neigborIndexes.size();i++)
			{
				if(talk)
				{
					mprintf("setting neigborIndexes:%d:%d: previous region:%d\n",i,neigborIndexes[i],indx2region[neigborIndexes[i]]);
				}
				int oldNeighborRegion=indx2region[neigborIndexes[i]];
				if(oldNeighborRegion!=thisRegion)
				{// merge 2 regions
					for(uint j=0;j<indx2region.size();j++)
					{
						if(indx2region[j]==oldNeighborRegion)
						{
							indx2region[j]=thisRegion;
						}
					}
				}
			}
		}
		else
		{// sorry, no neighbors 
			// this lonesome pixel is a new region and index
			*rp=indx2region.size();
			int thisRegion=nbRegions++;
			indx2region.push_back(thisRegion);
		}
//  	}else
//  	{
//  		*rp=0;// background index
//  	}
	}
//  	mprintf("CConnectedComponentLabelling:: region image\n");
//  	LabelImage3D regions=res;
//    	for(LabelImage3D::iteratorFast fp=regions.begin();fp!=regions.end();fp++){*fp=indx2region[*fp];}	
//  	ImLib3DDisplay::Show(regions);

    // ***** now compute new unique regions
	uint i;
	map<int,int> old2newregion;
	int nbNewRegions=0;
    // add  background region
	old2newregion[0]=nbNewRegions++;
	for(i=0;i<indx2region.size();i++)
	{
		int old=indx2region[i];
		if(old2newregion.find(old)==old2newregion.end())
		{
			old2newregion[old]=nbNewRegions++;
		}
	}
    // switch to new regions
	for(i=0;i<indx2region.size();i++)
	{
		indx2region[i]=old2newregion[indx2region[i]];
	}	


    // ***** now let's switch from indexes to regions
	for(rp=res.begin();rp!=res.end();rp++)
	{
		*rp=indx2region[*rp];
	}

}
#endif // _MorphologicalOperators_hxx
