/* 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   VoxelCoding.hpp*/
#ifndef _VoxelCoding_hpp
#define _VoxelCoding_hpp
#include<ImLib3D/MorphologicalOperators.hpp>
#include<ImLib3D/Threshold.hpp>


namespace IP3D
{
/*! #&VoxelCode&#     Sort of distance transform of a shape, starting from seed points, and using an 
    integer metric
        m0 is the "distance" associated with F-Neighbors (Cross7). 
        m1 is the "distance" associated with e-Neighbors (Planar19). 
        m2 is the "distance" associated with v-Neighbors (Cubic27). 
        Example values are m0=1,m1=2,m2=3, or m0=3,m1=4,m2=5. 
     <UL>Args:
      <LI>Shape            type:Mask3D
      <LI>Seed             type:Mask3D
      <LI>m0               type:int
      <LI>m1               type:int
      <LI>m2               type:int
      <LI>resulting codes  type:Image3Dreal<int>
     </UL>
~*/
		void VoxelCode(const Mask3D &shape,const Mask3D &seeds,int m0,int m1,int m2,LabelImage3D& res,
					  const Vect3Di *stopPoint=NULL);
};

namespace IP3D
{
/*! #&SSCode&#     VoxelCode starting from a single seed point
        For example, this is usefull for shortest path extraction
     <UL>Args:
      <LI>Shape            type:Mask3D
      <LI>p0               type:Vect3Di
      <LI>m0               type:int
      <LI>m1               type:int
      <LI>m2               type:int
      <LI>resulting codes  type:Image3Dreal<int>
      <LI>start            type:Vect3Di
     </UL>
~*/
	inline void SSCode(const Mask3D &shape,const Vect3Di &sourcePoint,int m0,int m1,int m2,LabelImage3D& res,
					  const Vect3Di *stopPoint=NULL)
	{
		Mask3D seeds(shape.Size());
		seeds.Fill(0);
		seeds(sourcePoint)=1;
		VoxelCode(shape,seeds,m0,m1,m2,res,stopPoint);
	}
};
namespace IP3D
{
/*! #&BSCode&#     VoxelCode starting from boundary points of shape
        This is pretty much like a distance transform
     <UL>Args:
      <LI>Shape            type:Mask3D
      <LI>m0               type:int
      <LI>m1               type:int
      <LI>m2               type:int
      <LI>resulting codes  type:Image3Dreal<int>
     </UL>
~*/
	inline void BSCode(const Mask3D &shape,int m0,int m1,int m2,LabelImage3D& res)
	{
		Mask3D seeds=shape;
		Mask3D t;
		Erosion(shape,MORPHO_Cross7,t);
		seeds-=t;
		VoxelCode(shape,seeds,m0,m1,m2,res);
	}
};
class CodeCluster;

namespace IP3D
{
namespace ShortestPathExtractionHelpers
{
	//! argh, call fct. FIXME, this would be prettier with libsigc++
	typedef bool (*StopFctType)(void *dat,const Vect3Di &p);
	//! argh, call fct + data. FIXME, this would be prettier with libsigc++
	typedef pair<StopFctType,void*> StopCBType;
}

/*! #&ShortestPathExtraction&#     Shortest path between two points using VoxelCoding
     <UL>Args:
      <LI>Shape      type:Mask3D
      <LI>start      type:Vect3Di
      <LI>end        type:Vect3Di
      <LI>path       type:vector<Vect3Di>
      <LI>res image  type:Mask3D
      <LI>sel        type:StructureElement
      <LI>ssCode     type:Image3Dreal<int>
     </UL>
~*/
	void ShortestPathExtraction(const Mask3D &shape,const Vect3Di &start,const Vect3Di &end,vector<Vect3Di> &path,
								Mask3D *pres=NULL,StructureElement *psel=NULL,
								const LabelImage3D *pssCode0=NULL,
								ShortestPathExtractionHelpers::StopCBType stopCB=
								ShortestPathExtractionHelpers::StopCBType(NULL,NULL));


};


//! cluster of conected voxels having same "code" value
class CodeCluster
{
	friend class SkeletonExtraction;
public:
	//! code value of all voxels of this cluster
	int value;
	//! geometric center of points in cluster
	Vect3Df center;
	//! neighboring clusters
	vector<CodeCluster *> neighbors;
	//! all points of this cluster
	vector<Vect3Di> points;
    //@{ \name flags describing cluster type
	bool isLM;
	bool isMerging;
	bool isDividing;
	bool isBranching;
	//@}
	//! user flag
	int stop;
private:
	CodeCluster(int _value):
		value(_value),stop(0)
	{
	}
	void Init()
	{
		center=points[0];
		for(uint i=1;i<points.size();i++){center+=points[i];}
		center/=points.size();
	}
	void FindType()
	{
		isLM=true;
		for(uint i=0;i<neighbors.size();i++)
		{
			if(neighbors[i]->value>value){isLM=false;break;}
		}
		int ct=0;
		for(uint i=0;i<neighbors.size();i++)
		{
			if(neighbors[i]->value==value-1){ct++;}
		}
		if(ct>=2){isMerging=true;}else{isMerging=false;}
		ct=0;
		for(uint i=0;i<neighbors.size();i++)
		{
			if(neighbors[i]->value==value+1){ct++;}
		}
		if(ct>=2){isDividing=true;}else{isDividing=false;}
		if(isDividing || isMerging){isBranching=true;}else{isBranching=false;}
	}
public:
	//! get neighboring clusters that have a lower code value
	void GetSuccessors(vector<CodeCluster *> &successors)
	{
		for(uint i=0;i<neighbors.size();i++)
		{
			if(neighbors[i]->value<value){successors.push_back(neighbors[i]);}
		}
	} 
	//! points in cluster that have maximal bsCode
	vector<Vect3Di> medialPoints;
	//! compute points in cluster that have maximal bsCode
	void FindMedialPoints(const LabelImage3D &bsCode)
	{
		int mx=bsCode(points[0]);
		for(uint i=1;i<points.size();i++){mx=max(mx,bsCode(points[i]));}
		for(uint i=0;i<points.size();i++){if(bsCode(points[i])==mx){medialPoints.push_back(points[i]);}}
	}
	//!  Get point in cluster that has maximal bsCode and is closest to ref point
	Vect3Di MedialPoint(const Vect3Df &ref=Vect3Di(0,0,0))
	{
		int imin=0;
		float d2=ref.Dist2(medialPoints[0]);
		for(uint i=1;i<medialPoints.size();i++)
		{
			float d=ref.Dist2(medialPoints[i]);
			if(d<d2){d2=d;imin=i;}
		}
		return medialPoints[imin];
	}
	//! determines neighboring clusters, using a position to cluster map
	void FindNeighbors(Container3D<CodeCluster *> &clustersIma)
	{
		// for each point in this cluster, check all if each of its neighbors is in another cluster
		// and add it if it's not there
		StructureElement sel(MORPHO_Cubic27);
		sel.RemoveCenter();
		for(uint pt=0;pt<points.size();pt++)
		{
			for(uint i=0;i<sel.size();i++)
			{
				Vect3Di npos=points[pt]+sel[i];
				if(!clustersIma.IsInside(npos)){continue;}
				CodeCluster *neighbor=clustersIma(npos);
				if(neighbor && neighbor->value!=value && abs(neighbor->value-value)<2)
				{
					if(find(neighbors.begin(),neighbors.end(),neighbor)==neighbors.end())
					{
						neighbors.push_back(neighbor);
					}
				}
			}
				
		}
//    		printf("neighbors for cluster %x %03d\n",this,value);
//    		for(uint i=0;i<neighbors.size();i++)
//    		{
//    			printf("cluster %x %03d:neighbor %2d: %x %03d\n",this,value,i,neighbors[i],neighbors[i]->value);
//    		}
	}
	
};
//! series of nodes linked by paths of points
class PathGraph
{
public:
	class Node;
	class Path
	{
		friend class PathGraph;
	public:
		Node *first,*last;
		vector<Vect3Di> points;
	private:
		Path(Node *_first,Node *_last,const vector<Vect3Di> &_points);
	};
	class Node
	{
		friend class PathGraph;
	public:
		Vect3Di point;
		vector<Path *> paths;
	private:
		Node(){;}
	};
	vector<Node *> nodes;
	vector<Path *> paths;
	Node *AddNode(){Node *node=new Node();nodes.push_back(node);return node;}
	RectZone3Di BoundingBox() const;
	Path *AddPath(Node *_first,Node *_last,const vector<Vect3Di> &_points)
	{
		Path *path=new Path(_first,_last,_points);paths.push_back(path);return path;
	}
};

//! Skeleton extraction 
class SkeletonExtraction
{
	static bool SpeStopCB(void *dat,const Vect3Di &p)
	{
		return ((SkeletonExtraction *)dat)->clustersIma(p)->stop!=-1;
	}
	class MedPathFind
	{
	public:
		CodeCluster *parentCluster;
		int id;
		int id1,id2;
		vector<Vect3Di> path;
		vector<CodeCluster *> cpath;
		vector<vector<int> > branches;
		vector<PathGraph::Node *> branchNodes;
		void ExtractPath(SkeletonExtraction *se);
		int FindCluster(CodeCluster *cluster)
		{
			int res=-1;
			for(uint i=0;i<cpath.size();i++){if(cpath[i]==cluster){res=i;break;}}
			return res;
		}
		int operator<(const MedPathFind &other) const {return cpath[0]->value<other.cpath[0]->value;}
		MedPathFind(CodeCluster *startCluster):
			id1(-1),id2(-1)
		{
			cpath.push_back(startCluster);
		}
	};
	friend class SkeletonExtraction::MedPathFind;
	//! shape from which to extract skeleton
	const Mask3D &shape;
	//! ssCode, "distance" field from seed point
	LabelImage3D ssCode;
	//! bsCode, "distance" field shape boundary
	LabelImage3D bsCode;
	//! each distinct value of ssCode that forms a set of conected points is grouped into a cluster
	vector<CodeCluster *> clusters;
	//! the cluster associated to each point
	Container3D<CodeCluster *> clustersIma;
	//! creates all clusters from ssCode
	void CreateClusters();
	//! initialises clusters, after all of them are created.
	void InitClusters();

	void PathGraphFromMedPaths(vector<MedPathFind> &medPaths,PathGraph &pathGraph);
	void ExtractAllMedPaths(vector<MedPathFind> &medPaths);
	void ExtractPaths(PathGraph &pathGraph);

	SkeletonExtraction(const Mask3D &_shape,const Vect3Di *psourcePoint=NULL);
	void MedPathFromPath(const vector<Vect3Di> &path,vector<Vect3Di> &medPath)
	{
		medPath.push_back(path[0]);
		for(uint i=1;i<path.size();i++)
		{
			medPath.push_back(clustersIma(path[i])->MedialPoint(medPath[i-1]));
		}
	}
public:
	static void ExtractAll(const Mask3D &shape,PathGraph &pathGraph);
	static void ExtractLinkingMedPath(const Mask3D &shape,const Vect3Di &start,const Vect3Di &end,vector<Vect3Di> &res);
	static void DrawPathGraph(const PathGraph &pathGraph,Mask3D &res);
};

namespace IP3D
{
/*! #&Skeleton&#     Centerline (1-D Skeleton) of a 3D shape
        The algorithm is described in Zhou and Toga 99. 
        WARNING: this implementation is still a bit buggy, and the algo 
        doesn't seem to be very good anyways.
     <UL>Args:
      <LI>Shape     type:Mask3D
      <LI>Skeleton  type:Mask3D
     </UL>
~*/
	inline void Skeleton(const Mask3D &shape,Mask3D& res)
	{
		PathGraph pathGraph;
		SkeletonExtraction::ExtractAll(shape,pathGraph);
		res.Resize(shape.Size());
		res.Fill(0);
		SkeletonExtraction::DrawPathGraph(pathGraph,res);
	}
};

#endif// _VoxelCoding_hpp
