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

#ifndef _SparseImage3D_hpp
#define _SparseImage3D_hpp


#include<ImLib3D/Image3D.hpp>
#include<ImLib3D/BitImage3D.hpp>
#include<boost/type_traits/arithmetic_traits.hpp>

#define UintToVect3Di(v) Vect3Di(v&((1<<10)-1),(v>>10)&((1<<10)-1),(v>>20)&((1<<10)-1))
#define Vect3DiToUint(p) (p.x+(p.y<<10)+(p.z<<20))

class SparseImage3DBase;
//! defines the constant sparse structure (which voxels are zero or non zero) of a sparse image. 
/*! Non zero values of a sparse image are stored in an array indexed by valNb. This class
 provides fast valNb to 3D Pos and 3D pos to valNb conversion. These conversions require
special datastructures that are automatically computed only if needed. */
class SparseStructure3D : public Streamable
{
	friend void SparseImage3D_Test();
	//! number of non zero values in mask
	int size;
	//! one bit per voxel image describing zero and non zero values
	BitMask3D mask;

    //@{ \name Several SparseImage3D's may share the same SparseStructure3D, refence counting is used for auto deallocation
private:
	//! if this is true, then object automatically manges it's own destruction
	bool autoDestroy;
//  	//! number of times this object is used. when it goes down to zero, it should delete itself
//  	mutable int refCount;
	//! images referencing this  object. when it goes down to zero, it should delete itself
	mutable vector<SparseImage3DBase *> referencingImages;
	//! indicates whether that call to the destructor is not a user error
	mutable bool suicideInProgress;
public:
	void Reference(SparseImage3DBase *image)   const {referencingImages.push_back(image);}
	bool Dereference(SparseImage3DBase *image) const ;
    //@}
    //@{ \name Several SparseImage3D files may share the same SparseStructure3D file
private:
	//! keep track of existing structs, so that anybody that needs one can hook on it by name
	static vector<const SparseStructure3D *> allreadyExisting;
	mutable string name;
public:
	void SetName(string _name){name=_name;}
	string  Name() const {return name;}
	static const SparseStructure3D *Find(string sname);
    //@}
private:
    //@{ \name Data structures for rapid valnb (value index) to 3D position conversion
	mutable bool valNbToPosOk;
	mutable vector<uint> valNbToPosIndex;
	void MakeValNbToPosIndex() const
	{
		valNbToPosOk=true;
		valNbToPosIndex.resize(size);
		int valNb=0;
		for(BitMask3D::const_iteratorXYZ p=mask.begin();p!=mask.end();++p)
		{
			if(*p)
			{
				valNbToPosIndex[valNb]=Vect3DiToUint(p);
				valNb++;
			}
		}
	}
    //@}

    //@{ \name Data structures for rapid 3D position to valnb (value index) conversion
	mutable bool posToValNbOk;
	mutable Container3D<vector<uint> > xProjValNb;
	void MakePosToValNb() const
	{
		posToValNbOk=true;
		xProjValNb.Resize(Size3D(mask.Height(),mask.Depth(),1));
		if(!valNbToPosOk){MakeValNbToPosIndex();}
		// we want each vector in xProjValNb to be of exact good size
		// so first find exact sizes
		LabelImage3D sizes(mask.Height(),mask.Depth(),1);
		sizes.Fill(0);

		for(uint valNb=0;valNb<valNbToPosIndex.size();valNb++)
		{
			Vect3Di pos=UintToVect3Di(valNbToPosIndex[valNb]);
			sizes(pos.y,pos.z,0)++;
		}
		// now correctly resize each one
		for(LabelImage3D::iteratorXYZ p=sizes.begin();p!=sizes.end();++p)
		{
			xProjValNb(p.Pos()).resize(*p);
		}
		// finally fill in xProjValNb, using sizes as a counter
		sizes.Fill(0);
		for(uint valNb=0;valNb<valNbToPosIndex.size();valNb++)
		{
			Vect3Di pos=UintToVect3Di(valNbToPosIndex[valNb]);
			xProjValNb(pos.y,pos.z,0)[sizes(pos.y,pos.z,0)++]=valNb;
		}
	}
    //@}
public:
	void Read( ImLib3DFile *file,XMLNode *parentNode=NULL,XMLNode *node=NULL)       ;
	void Write(ImLib3DFile *file,XMLNode *parentNode=NULL,XMLNode *node=NULL) const ;
	void WriteToFile (const string &fname) const {name=fname;Streamable::WriteToFile (fname);}
	void ReadFromFile(const string &fname)       {name=fname;Streamable::ReadFromFile(fname);}

	int ValSize() const {return size;}
	Size3D Size() const {return mask.Size();}
	const BitMask3D &Mask() const {return mask;}

	Vect3Di valNbToPos(int valNb)  const
	{
		if(!valNbToPosOk){MakeValNbToPosIndex();}
		return UintToVect3Di(valNbToPosIndex[valNb]);
	}
	int posToValNb(const Vect3Di &pos)  const
	{
		if(!mask(pos)){return -1;}
		if(!posToValNbOk){MakePosToValNb();}
		const vector<uint> &proj=xProjValNb(pos.y,pos.z,0);
		vector<uint>::const_iterator p=proj.begin();
		for(p=proj.begin();p!=proj.end();++p){if((valNbToPosIndex[*p]&((1<<10)-1))==(uint)pos.x){return *p;}}
		ThrowError("SparseStructure3D::posToValNb, point in mask, but not xProjValNb");
		return -1;
	}
	
public:
	virtual ~SparseStructure3D();

public:
	static SparseStructure3D *Create(const string &name);
	static SparseStructure3D *Create(ImLib3DFile *file,XMLNode *parentNode=NULL,XMLNode *node=NULL);
	static SparseStructure3D *Create(const BitMask3D &_mask,bool _autoDestroy);
private:
	void ComputeValSize()
	{
		size=0;
		for(BitMask3D::iteratorXYZ p=mask.begin();p!=mask.end();++p)
		{
			if(*p){size++;}
		}
	}
	SparseStructure3D();
	SparseStructure3D(BitMask3D _mask,bool _autoDestroy=false,int _size=-1);
};

//! SparseImage3D iterator, this just encapulates a float pointer
template<class ImageType,class _Ptr,class _Ref>
class SIM3D_iteratorFast : public random_access_iterator<float,size_t>
{
  	typedef SIM3D_iteratorFast<ImageType,_Ptr,_Ref> _Self;
public:
	ImageType *ima;
	_Ptr pos;
public:
	void AssertBounds()
	{
		if(pos<ima->begin().pos || pos>=ima->end().pos)
		{
			ThrowError("iteratorFast: out of bounds:pos:%x begin:%x end:%x\n",pos,ima->begin().pos,ima->end().pos);
		}
	}
	inline _Ref   operator*(){return *pos;}
	double Progress() const {if(!this){return 0;}return (_Self(ima,pos)-ima->begin())/(float)(ima->end()-ima->begin());}
	inline _Self  operator++(int)     {_Self __tmp=*this;pos++;return __tmp;}
	inline _Self &operator++()        {pos++;return *this;}
	inline void   operator+=(size_t m){pos+=m;}
	inline _Self  operator+ (size_t m){_Self __tmp=*this;__tmp.pos+=m;return __tmp;}
	inline _Self  operator- (size_t m){_Self __tmp=*this;__tmp.pos-=m;return __tmp;}
	inline _Self  operator--(int)     {_Self __tmp=*this;pos--;return __tmp;}
	inline _Self &operator--()        {pos--;return *this;}
	inline void   operator-=(size_t m){pos-=m;}
	inline int    operator- (const _Self &other){return pos-other.pos;}
	inline bool   operator< (const _Self &other){return pos< other.pos;}
	inline bool   operator<=(const _Self &other){return pos<=other.pos;}
	inline bool   operator> (const _Self &other){return pos> other.pos;}
	inline bool   operator>=(const _Self &other){return pos>=other.pos;}

	inline Vect3Di ComputePosition() const
	{
		return ima->sparseStructure->valNbToPos(pos-ima->begin().pos);
	}
	inline Vect3Di Pos() const
	{
		return ima->sparseStructure->valNbToPos(pos-ima->begin().pos);
	}
	inline bool  operator!=(const _Self  &other){return pos!=other.pos; }
	inline bool  operator==(const _Self  &other){return pos==other.pos; }
	SIM3D_iteratorFast():ima(NULL),pos(0){}
	SIM3D_iteratorFast(ImageType *_ima,_Ptr _pos):ima(_ima),pos(_pos){}
};



class SparseImage3DBase : public Image3D
{
public:
	void Invalidate()
	{
		sparseStructure=NULL;
	}
protected:
	const SparseStructure3D *sparseStructure;
	bool shareSparseStructureFile;
	void SetSparseStructure(const SparseStructure3D *_sparseStructure)
	{
		sparseStructure=_sparseStructure;sparseStructure->Reference(this);
		Resize(sparseStructure->Size());
	}
  	SparseImage3DBase():Image3D(),sparseStructure(NULL),shareSparseStructureFile(false)
	{;}
  	SparseImage3DBase(const SparseImage3DBase &other):	
		Image3D(other),sparseStructure(NULL),shareSparseStructureFile(false)
	{
		SetSparseStructure(other.sparseStructure);
	}
  	SparseImage3DBase(const SparseStructure3D *_sparseStructure):
		Image3D(_sparseStructure->Size()),
		sparseStructure(NULL),
		shareSparseStructureFile(false)
	{
		SetSparseStructure(_sparseStructure);
	}
};

//! set to zero only if second template arg is true
template<class T,bool>
struct OptionallySetToZero
{
	static void Set(T &){}
};
template<class T>
struct OptionallySetToZero<T,true>
{
	static void Set(T &v){v=0;}
};


//! An image where most points are zero. 
/*! The points which are zero or non zero are determined once, and connot be changed.
Though the values of non zero points may be changed. The zero/non zero structure
is described by  SparseStructure3D. Several SparseImage3D's may share the same
SparseStructure3D to save memory. When saved into a file they may also share a common
SparseStructure3D file.*/
template<class Im3DValue>
class SparseImage3D : public SparseImage3DBase//, public ArithmeticOperators<SparseStructure3D,float>
{
	typedef SparseImage3D<Im3DValue> _Self;

	//! Initialize zero value only for types for which it has a meaning
	void InitZero()
	{
		OptionallySetToZero<Im3DValue,boost::is_arithmetic<Im3DValue>::value>::Set(zero);
	}

	Im3DValue *imageData;
	Im3DValue *imageDataEnd;
public:
	//! So that other templates know 
	typedef Im3DValue value_type;
public:
    //@{ \name Iterator Typedefs    

	typedef SIM3D_iteratorFast<      _Self,      value_type *,      value_type &>       iteratorFast;
	typedef SIM3D_iteratorFast<const _Self,const value_type *,const value_type &> const_iteratorFast;
	typedef iteratorFast iterator;
	typedef iterator iteratorMasked;
	typedef iterator iteratorXYZ;
	typedef const_iteratorFast const_iterator;
	typedef const_iterator const_iteratorMasked;
	typedef const_iterator const_iteratorXYZ;
	friend  class SIM3D_iteratorFast<      _Self,      value_type *,      value_type &>;
	friend  class SIM3D_iteratorFast<const _Self,const value_type *,const value_type &>;
public:
	//@{ \name Begin and end iterating through an image
	//! (stl-like) get an  iterator positioned at the begining of this container
	iteratorFast         begin()       {return       iteratorFast(this,imageData);}
  	const_iteratorFast   begin() const {return const_iteratorFast(this,imageData);}
	//! (stl-like) get an  iterator positioned at the end of this container
	inline iteratorFast         end()         {return       iteratorFast(this,imageDataEnd);}
	inline const_iteratorFast   end()   const {return const_iteratorFast(this,imageDataEnd);}

	value_type zero;
	value_type &Get(const Vect3Di &pos)
	{
		int valNb=sparseStructure->posToValNb(pos);
		return valNb<0 ? zero : imageData[valNb];
	}
	value_type Get(const Vect3Di &pos) const
	{
		int valNb=sparseStructure->posToValNb(pos);
		return valNb<0 ? zero : imageData[valNb];
	}
public:
	inline value_type &operator()(int x,int y,int z)        {return Get(Vect3Di(x,y,z));}
	inline value_type  operator()(int x,int y,int z)  const {return Get(Vect3Di(x,y,z));}
	inline value_type &operator()(const Vect3Di &pos)       {return Get(pos);}
	inline value_type  operator()(const Vect3Di &pos) const {return Get(pos);}

	//! Sets all elements of the container to this value 
	void  Fill(const value_type &value){fill(begin(),end(),value);}
	template<class OtherImageType >
	void MaskedCopy(OtherImageType other){for(iterator p=begin();p!=end();++p){*p=other(p.Pos());}}

	//@{ \name File READ/WRITE methods
public:
	virtual string GetTypeName() const {return "SparseImage3D";}
	//! This is where all the information is written
	/*! This virtual function writes xml/bin of any unset/unwritten  information
        (type,image size,image data,,...)
	    If you override it in a child class, you should still call it (using
		_Base::Write...), you can even call it several times...
	 */
	virtual void Read( ImLib3DFile *file,XMLNode *parentNode=NULL,XMLNode *node=NULL)      ;
	virtual void Write(ImLib3DFile *file,XMLNode *parentNode=NULL,XMLNode *node=NULL) const;
public:
	//! Write the image to file fname: This is the fct the end user should call
	virtual void WriteToFile (const string &fname) const {Streamable::WriteToFile(fname);}
	virtual void WriteToFile (const string &fname,bool _shareSparseStructureFile) 
	{
		shareSparseStructureFile=_shareSparseStructureFile;
		Streamable::WriteToFile(fname);
	}
	//@}
public:
	void Resize(const Size3D &newSize);
	const value_type &GetZero() const {return zero;}
	void SetZero(const value_type& _zero) {zero=_zero;}

    //@{ \name Copy and compare operators
	_Self& operator= (const _Self & other);
  	bool   operator==(const _Self & other) const;
  	bool   operator!=(const _Self & other) const {return !(this->operator==(other));}
  	template<class OtherImage >
  	bool operator== (const OtherImage & other){return Size()==other.Size() && equal(other.begin(),other.end(),begin());}
  	template<class OtherImage >
  	bool operator!= (const OtherImage & other){return !((*this)==other);}
	//@}

//	value_type      SafeValue(int pos)     const {return IsInside(pos) ?           Get(pos) : 0;}
  	value_type      SafeValue(int x,int y,int z) const 	
	{return IsInside(x,y,z) ?           (*this)(x,y,z) : 0 ;}
  	value_type      SafeValue(const Vect3Di &v)  const 	
	{return IsInside(v    ) ?           (*this)(v    ) : 0 ;}

	// **************************************************
	// ****************** constructors/destructors*******
	// **************************************************
	//! \name Constructors and destructors
	//@{
public:
	//! Destructor
	virtual ~SparseImage3D();
private:
	void Allocate();
public:
  	SparseImage3D();
  	explicit SparseImage3D(const string& fname):SparseImage3DBase(),imageData(NULL),imageDataEnd(NULL)
	{InitZero();ReadFromFile(fname);}
  	explicit SparseImage3D(const char   *fname):	SparseImage3DBase(),imageData(NULL),imageDataEnd(NULL)
	{InitZero();ReadFromFile(fname);}
  	explicit SparseImage3D(const SparseStructure3D *_sparseStructure);
	/*! Copy constructor */
	SparseImage3D(const SparseImage3D & origCont);
	//@}
};


#include<ImLib3D/SparseImage3D.hxx>

typedef SparseImage3D<float> SparseImage3Df;
#endif // _SparseImage3D_hpp
