/* 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.hpp */

#ifndef _MorphologicalOperators_hpp
#define _MorphologicalOperators_hpp

#include<ImLib3D/ImageProcessor.hpp>
#include<map>

//! Neighborhoods for morphological operations
/*! These are several values used for creating neighborhoods for morphological operations 
 */
enum StructureElementType
{
	MORPHO_CUSTOM,  //!< User defined neighborhood
	MORPHO_FORWARD,  //!< For internal dist-transform use (all points in mask before reaching center, in scan order)
	MORPHO_BACKWARD, //!< For internal dist-transform use (all points in mask after reaching center, in scan order)
	MORPHO_ALLDIRECTIONS, //!< For internal dist-transform use 
	MORPHO_Cubic27, //!< Cas du voisinage cube 
	MORPHO_Cross7,  //!< Cas du voisinage croix 
	MORPHO_Planar19,//!< Cas du voisinage plans 
	MORPHO_Corner9 ,//!< Cas du voisinage coin 
	MORPHO_Pts81,		
	MORPHO_Pts14,
	MORPHO_Pts4,
	MORPHO_Pts10,
	MORPHO_Pts5,
	MORPHO_Pts41,
	MORPHO_XYCross5, //! Planar xy
	MORPHO_XZCross5, //! Planar xz
	MORPHO_YZCross5, //! Planar yz
	MORPHO_SphereR3, //! Sphere with radius 3
	MORPHO_SphereR4  //! Sphere with radius 4
};

//! Structure Element for Morphological Operators (also called neighborhood or mask)
/*! Coefficients may be associated to neighborhood points, this is used
    for computing distance transforms .  
 */
class StructureElement
{
protected:
	//! Original mask type
	StructureElementType maskType;
	//! Create mask from maskPoints
	void CreateMask();
public:
	//! Computes the smallest and largest positions in maskPoints 
	void ComputeMaskPointsSize(Vect3Di &minv,Vect3Di &maxv) const;
protected:
	//! Creates maskPoints from maskType
	void CreateMaskPoints(StructureElementType maskType);
	//! Creates maskPoints from mask
	void CreateMaskPoints();

	//! Set of points that define this mask
	vector<Vect3Di> maskPoints;
	//! Coefficients associated to these points (used by distance transform)
	vector<float> coeffs;
	//! Mask image, same points as in maskPoints
	/*! The mask's iCenter gives the center of the mask
	 */
	Mask3D *mask;
public:
	static void CreateSpherePointList(vector<Vect3Di>& pl, int radius);
	void ComputeCoeffs();
	//! Display this mask (for debugining)
	void Show(const string &name) const;
	//@{ \name Accesors to specific data positions that will return Zero for position outside of scope
	inline const Vect3Di &GetPoint(int i) const {return maskPoints[i];}
	inline const Vect3Di &operator[](uint i) const {return maskPoints[i];}
	inline float GetCoeff(int i) const {return coeffs[i];}
    //! zone in an image where it is safe to consider all neighbors
	RectZone3Di SafeZone(const Size3D &imageSize) const
	{
		Vect3Di vmin,vmax;
		ComputeMaskPointsSize(vmin,vmax);
		return RectZone3Di(Vect3Di(0,0,0),imageSize).Intersection(
			RectZone3Di(-vmin,imageSize.GetV()-vmax-Vect3Di(1,1,1)));
	}
	//! Find point number given it's position
	int   FindPoint(const Vect3Di &P) const 
	{
		for(uint i=0;i<maskPoints.size();i++)
		{if(P==maskPoints[i]){return i;}}
		return -1;
	}
	//! Find point's coefficient given it's position
	float GetCoeff(const Vect3Di &P) const 
	{
		int i=FindPoint(P);
		if(i<0)ThrowError("GetCoeff(P=%d %d %d) point not found",P.x,P.y,P.z);
		return coeffs[i];
	}
	const Mask3D  &Mask() const {return *mask;}
	inline int  GetNbPoints() const {return maskPoints.size();}
	inline uint size() const {return maskPoints.size();}
	//@}
	
	//! remove the center point of the mask 
	/*! Both in mask and in maskPoints
	 */
	void RemoveCenter();
	
	//! Build a  mask given a description string
	StructureElement(const string &s);

	//! Build a weighted mask (for distance transform) from an original std mask
	StructureElement(const StructureElement &other,StructureElementType direction);

	//! Build a  mask given it's description
	StructureElement(StructureElementType _maskType):
		maskType(_maskType)
	{
		CreateMaskPoints(maskType);
		CreateMask();
	}
	//! Build a  mask given it's points
	StructureElement(const vector<Vect3Di> &_maskPoints):
		maskType(MORPHO_CUSTOM)
	{
		maskPoints=_maskPoints;
		CreateMask();
	}
};

//! This class helps fast iterating through all neighbors of each point of an image
/*! Two situations are distinguished, if all neighbors are inside image (safe)
and if some neighbors are outside image (unsafe) */
template<class ImageType,class Inst=typename HasPointerIterators<ImageType>::value >
class ImageNeighbors
{
	typedef typename ImageType::value_type Im3DValue;
	typedef typename ImageType::iteratorFast Iterator;
	ImageType &image;
	const StructureElement &sel;
	/*! using this instead of iterators is an important
	  speedup when compiling in debug, (no difference in optimized..)*/
	Im3DValue ** iterators;//vector<Iterator> iterators;
	typename ImageType::iteratorXYZ p;
	uint npts;
	RectZone3Di safezone;
	int lastY;
	int lastX;
public:
	inline bool IsSafe(const Vect3Di &pos){return safezone.IsInside(pos);}
	inline void MoveTo(const Vect3Di &pos)
	{
		if(lastY!=pos.y || lastX!=pos.x-1)
		{// if we've changed line we have to recompute iterators (slow)
			for(uint i=0;i<npts;i++)
			{
				iterators[i]=(Im3DValue *)&image(pos+sel[i]);
//				iterators[i].MoveTo(pos+sel[i]);
			}
			lastY=pos.y;
			lastX=pos.x;
		}
		else
		{// fast!, just increment iterators
			for(uint i=0;i<npts;i++){iterators[i]++;}
			lastX++;
		}
	}
	inline void MoveBackwardsTo(const Vect3Di &pos)
	{
		if(lastY!=pos.y || lastX!=pos.x+1)
		{// if we've changed line we have to recompute iterators (slow)
			for(uint i=0;i<npts;i++)
			{
				iterators[i]=(Im3DValue *)&image(pos+sel[i]);
//				iterators[i].MoveTo(pos+sel[i]);
			}
			lastY=pos.y;
			lastX=pos.x;
		}
		else
		{// fast!, just decrement iterators
			for(uint i=0;i<npts;i++){iterators[i]--;}
			lastX--;
		}
	}
	inline Im3DValue &operator[](int i){return *(iterators[i]);}
	~ImageNeighbors(){delete [] iterators;}
	ImageNeighbors(ImageType &_image,const StructureElement &_sel):
		image(_image),sel(_sel),
		iterators(new Im3DValue*[sel.size()]),//iterators(sel.size()),
		npts(sel.size()),
		safezone(sel.SafeZone(image.Size())),lastY(-1),lastX(-2)
	{
        //for(uint i=0;i<npts;i++){iterators[i]=image.begin();}
	}
};

template<class ImageType>
struct IteratorType
{
	typedef typename ImageType::iterator type;
};
template<class ImageType>
struct IteratorType<const ImageType>
{
	typedef typename ImageType::const_iterator type;
};

template<class ImageType>
class ImageNeighbors<ImageType,__false_type>
{
	typedef typename ImageType::value_type Im3DValue;
	typedef typename IteratorType<ImageType>::type Iterator;
	ImageType &image;
	const StructureElement &sel;
	/*! using this instead of iterators is an important
	  speedup when compiling in debug, (no difference in optimized..)*/
	vector<Iterator> iterators;
	typename ImageType::iteratorXYZ p;
	uint npts;
	RectZone3Di safezone;
	int lastY;
	int lastX;
public:
	inline bool IsSafe(const Vect3Di &pos){return safezone.IsInside(pos);}
	inline void MoveTo(const Vect3Di &pos)
	{
		if(lastY!=pos.y || lastX!=pos.x-1)
		{// if we've changed line we have to recompute iterators (slow)
			for(uint i=0;i<npts;i++)
			{
				iterators[i]=pos+sel[i];
			}
			lastY=pos.y;
			lastX=pos.x;
		}
		else
		{// fast!, just increment iterators
			for(uint i=0;i<npts;i++){iterators[i]++;}
			lastX++;
		}
	}
	inline void MoveBackwardsTo(const Vect3Di &pos)
	{
		if(lastY!=pos.y || lastX!=pos.x+1)
		{// if we've changed line we have to recompute iterators (slow)
			for(uint i=0;i<npts;i++)
			{
//				iterators[i]=(Im3DValue *)&image(pos+sel[i]);
				iterators[i]=pos+sel[i];
			}
			lastY=pos.y;
			lastX=pos.x;
		}
		else
		{// fast!, just decrement iterators
			for(uint i=0;i<npts;i++){iterators[i]--;}
			lastX--;
		}
	}
	inline Im3DValue operator[](int i){return *(iterators[i]);}
	~ImageNeighbors(){}
	ImageNeighbors(ImageType &_image,const StructureElement &_sel):
		image(_image),sel(_sel),
//		iterators(new Im3DValue*[sel.size()]),
		iterators(sel.size()),
		npts(sel.size()),
		safezone(sel.SafeZone(image.Size())),lastY(-1),lastX(-2)
	{
        for(uint i=0;i<npts;i++){iterators[i]=image.begin();}
	}
};



namespace IP3D
{
/*! #&MedianFilter&#     Median filter.
        Replaces each voxel by the median value of all of it's neighbors. Neighbors 
        are defined by structure element argument. Ignores masks.
     <UL>Args:
      <LI>Source image      type:ImageType
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:ImageType
     </UL>
~*/
	template <class ImageType>
	void MedianFilter(const ImageType& src, const StructureElement &mask, ImageType& res);
};

namespace IP3D
{
/*! #&SharpeningFilter&#     Sharpening filter.
        Replaces each voxel by the min or max value of all of it's neighbors, 
        whichever value is closer. If both are equivalent, leaves voxel untouched. 
        Neighbors are defined by structure element argument. Ignores masks.
     <UL>Args:
      <LI>Source image      type:ImageType
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:ImageType
     </UL>
~*/
	template <class ImageType>
	void SharpeningFilter(const ImageType& src, const StructureElement &mask, ImageType& res);
};


namespace IP3D
{
/*! #&Erosion&#     Morphological erosion.
        Replaces each voxel by the minimum value of all of it's neighbors. Ignores 
        masks.
     <UL>Args:
      <LI>Source image      type:ImageType
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:ImageType
     </UL>
~*/
	template <class ImageType>
	void Erosion(const ImageType& src, const StructureElement &mask, ImageType& res);
};

namespace IP3D
{
/*! #&Dilation&#     Morphological dilation.
        Replaces each voxel by the maximum value of all of it's neighbors. Ignores 
        masks.
     <UL>Args:
      <LI>Source image      type:ImageType
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:ImageType
     </UL>
~*/
	template <class ImageType>
	void Dilation(const ImageType& src, const StructureElement &mask, ImageType& res);
};

namespace IP3D
{
/*! #&Opening&#     Morphological opening.
        Succesively applies an erosion and a dilation.
     <UL>Args:
      <LI>Source image      type:ImageType
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:ImageType
     </UL>
~*/
	template <class ImageType>
	void Opening(const ImageType& src, const StructureElement &mask, ImageType& res)
	{
		IP3D::Erosion(src,mask,res);
		IP3D::Dilation(res,mask,res);
	}
};

namespace IP3D
{
/*! #&Closing&#     Morphological closing.
        Succesively applies an dilation and a erosion.
     <UL>Args:
      <LI>Source image 1    type:ImageType
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:ImageType
     </UL>
~*/
	template <class ImageType>
	void Closing(const ImageType& src, const StructureElement &mask, ImageType& res)
	{
		IP3D::Dilation(src,mask,res);
		IP3D::Erosion(res,mask,res);
	}
};

namespace IP3D
{
/*! #&DistanceTransform&#     Distance to closest non zero point.
        Approximate distance within set, from inside outside of set. This is a fast, 
        two pass algorithm. The StructureElement contains coefficients that define 
        the distance to each neighbor. FIXME: These coeficients are currently not 
        optimal, and give a slightly biased distance. By default outside of image 
        is considered outside of shape. If border extend is true, outside of image 
        ignored
     <UL>Args:
      <LI>Source image      type:Mask3D
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:Image3Df
      <LI>border extend     type:bool
     </UL>
~*/
		void DistanceTransform(const Mask3D& src, const StructureElement &mask, Image3Df& res,bool borderExtend=false);
};

namespace IP3D
{
/*! #&ConnectedComponentLabelling&#     Creates an integer valued image where each set of connected neighbors of source 
    image has a seperate value.
        Result values go from 0 to n. Where n is the total number of separate 
        connected components in the image. Definition of neighbors uses 
        StructureElement.
     <UL>Args:
      <LI>Source image      type:Mask3D
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:Image3Dreal<int>
     </UL>
~*/
	template <class ImageType>
	void ConnectedComponentLabelling(const ImageType& src, const StructureElement &mask0, LabelImage3D& res,typename ImageType::value_type *background0=NULL);
};
template <class ImageType> inline pair<bool,typename ImageType::value_type> DefaultBackground(){return pair<bool,typename ImageType::value_type>(false,typename ImageType::value_type());}
template <> inline pair<bool,byte> DefaultBackground<Mask3D>(){return pair<bool,byte>(true,0);}

namespace IP3D
{
/*! #&FillHoles&#     Sets to 1 all points that are not connected to the outside of the source image.
        Definition of neighbors uses StructureElement. Algorithm simply uses 
        ConnectedComponentLabelling of an enlarged image.
     <UL>Args:
      <LI>Source image      type:Mask3D
      <LI>StructureElement  type:StructureElement
      <LI>Result image      type:Mask3D
     </UL>
~*/
		void FillHoles(const Mask3D& src, const StructureElement &mask, Mask3D& res);
};

class CMakeValueList;


//! Container for information on distinct image values
template<class Im3DValue>
class ImageValues
{
	friend class CMakeValueList;
public:
	class ImageValue
	{
	public:
		int size;// number of voxels having this value
		Vect3Di point;// position of a voxel having this value
		ImageValue(){size=0;point.x=-1;point.y=-1;}
	};
	typedef map<Im3DValue,ImageValue> MapType;
	map<Im3DValue,ImageValue> values;
	ImageValue operator[](const Im3DValue &v)
	{
		return values[v];
	}

};
namespace IP3D
{
/*! #&MakeValueList&#     Maps value to size, used for getting largest connected component... 
        Needs name cleanup, will be renamed. Ignores masks.
     <UL>Args:
      <LI>Source image        type:ImageType
      <LI>Maps value to size  type:ImageValues
     </UL>
~*/
	template <class ImageType>
	void MakeValueList(const ImageType& src, ImageValues<typename ImageType::value_type>& res);
};

namespace IP3D
{
/*! #&LargestConnectedComponent&#     Selects the largest connected component. Usefull after thresholding.
        Uses ConnectedComponentLabelling
     <UL>Args:
      <LI>Source image  type:Mask3D
      <LI>Result image  type:Mask3D
     </UL>
~*/
		void LargestConnectedComponent(const Mask3D &src,Mask3D &res);
};


#include<ImLib3D/MorphologicalOperators.hxx>

#endif // _MorphologicalOperators_hpp
