TriangleMesh.cpp

		//constants are parameters of Andor iXon camera provided from vendor:
		//photon shot noise: dark current
		*p+=(float)GetRandomPoisson(0.06f);

		//read-out noise:
		//  variance up to 25.f (old camera on ILBIT)
		//  variance about 1.f (for the new camera on ILBIT)
		*p+=GetRandomGauss(480.f,20.f);
		//*p+=530.f;
	}
#ifdef SAVE_INTERMEDIATE_IMAGES
	dt.SaveImage("6_texture_filtered_resampled_finalized.ics");
#endif

	//obtain final GRAY16 image
	i3d::FloatToGrayNoWeight(dt,texture);
}


void ActiveMesh::CenterMesh(const Vector3F& newCentre)
{
	//calc geom. centre
	double x=0.,y=0.,z=0.;
	for (unsigned int i=0; i < Pos.size(); ++i)
	{
		x+=Pos[i].x;
		y+=Pos[i].y;
		z+=Pos[i].z;
	}
	x/=double(Pos.size());
	y/=double(Pos.size());
	z/=double(Pos.size());

	//std::cout << "mesh centre is: " << x << "," << y << "," << z << "\n";

	x-=newCentre.x;
	y-=newCentre.y;
	z-=newCentre.z;

	//shift the centre to point (0,0,0)
	for (unsigned int i=0; i < Pos.size(); ++i)
	{
		Pos[i].x-=float(x);
		Pos[i].y-=float(y);
		Pos[i].z-=float(z);
	}
	for (unsigned int i=0; i < fPoints.size(); ++i)
	{
		fPoints[i].x-=float(x);
		fPoints[i].y-=float(y);
		fPoints[i].z-=float(z);
	}
}


void ActiveMesh::ScaleMesh(const Vector3F& scale)
{
	for (unsigned int i=0; i < Pos.size(); ++i)
	{
		Pos[i].x*=scale.x;
		Pos[i].y*=scale.y;
		Pos[i].z*=scale.z;
	}
}


void ActiveMesh::TranslateMesh(const Vector3F& shift)
{
	for (unsigned int i=0; i < Pos.size(); ++i)
	{
		Pos[i].x+=shift.x;
		Pos[i].y+=shift.y;
		Pos[i].z+=shift.z;
	}
}


// ============== surface fitting ============== 
int ActiveMesh::CalcQuadricSurface_Taubin(const int vertexID,
                                          float (&coeffs)[10])
{
	//determine some reasonable number of nearest neighbors
	std::vector< std::vector<size_t> > neigsLeveled;
	ulm::getVertexNeighbours(*this,vertexID,2,neigsLeveled);

	//make it flat...
	std::vector<size_t> neigs;
	for (unsigned int l=0; l < neigsLeveled.size(); ++l)
		for (unsigned int i=0; i < neigsLeveled[l].size(); ++i)
			neigs.push_back(neigsLeveled[l][i]);
	neigsLeveled.clear();

	//V be the vector of [x,y,z] combinations, a counterpart to coeffs
	float V[10],V1[10],V2[10],V3[10];

	//M be the matrix holding initially sum of (square matrices) V*V'
	//N is similar, just a sum of three different Vs is used
	float M[100],N[100];
	for (int i=0; i < 100; ++i) M[i]=N[i]=0.f;

	//over all neighbors (including the centre vertex itself),
	//
	//this order garuantees that array V will be relevant for input
	//vertexID after the cycles are over -- will become handy later
	for (signed int n=(int)neigs.size()-1; n >= 0; --n)
	{
		//shortcut to the current point
		const Vector3FC& v=Pos[neigs[n]];

		//get Vs for the given point 'v'
		V[0]=1.f;       V1[0]=0.f;      V2[0]=0.f;      V3[0]=0.f;
		V[1]=v.x;       V1[1]=1.f;      V2[1]=0.f;      V3[1]=0.f;
		V[2]=v.y;       V1[2]=0.f;      V2[2]=1.f;      V3[2]=0.f;
		V[3]=v.z;       V1[3]=0.f;      V2[3]=0.f;      V3[3]=1.f;
		V[4]=v.x*v.y;   V1[4]=v.y;      V2[4]=v.x;      V3[4]=0.f;
		V[5]=v.x*v.z;   V1[5]=v.z;      V2[5]=0.f;      V3[5]=v.x;
		V[6]=v.y*v.z;   V1[6]=0.f;      V2[6]=v.z;      V3[6]=v.y;
		V[7]=v.x*v.x;   V1[7]=2.f*v.x;  V2[7]=0.f;      V3[7]=0.f;
		V[8]=v.y*v.y;   V1[8]=0.f;      V2[8]=2.f*v.y;  V3[8]=0.f;
		V[9]=v.z*v.z;   V1[9]=0.f;      V2[9]=0.f;      V3[9]=2.f*v.z;

		//construct V*V' and add it to M and N
		for (int j=0; j < 9; ++j)  //for column
		 for (int i=0; i < 9; ++i) //for row;  note the order optimal for Fortran
		 {
			//C order                 (row-major):  M_i,j -> M[i][j] -> &M +i*STRIDE +j
			//Lapack/Fortran order (column-major):  M_i,j -> M[i][j] -> &M +j*STRIDE +i
			//const int off=i*10 +j; //C
			const int off=j*9 +i; //Fortran

			M[off]+= V[j+1]* V[i+1];
			N[off]+=V1[j+1]*V1[i+1];
			N[off]+=V2[j+1]*V2[i+1];
			N[off]+=V3[j+1]*V3[i+1];
		 }
	}

	//now, solve the generalized eigenvector of the matrix pair (M,N):
	//MC = nNC
	//
	//M,N are (reduced) 9x9 matrices constructed above,
	//C is vector (infact, the coeff), n is scalar Lagrange multiplier
	//
	//if M,N were (full-size) 10x10 matrices, the N would be singular
	//as the 1st row would contain only zeros,
	//it is therefore reduced to 9x9 sacrifing the first row
	//
	//according to netlib (Lapack) docs, http://www.netlib.org/lapack/lug/node34.html
	//type 1, Az=lBz -- A=M, B=N, z=C
	//function: SSYGV

/*

	lapack_int itype=1;
	char jobz='V';
	char uplo='U';
	lapack_int n=9;
	float w[10];
	float work[512];
	lapack_int lwork=512;
	lapack_int info;
	LAPACK_ssygv(&itype,&jobz,&uplo,&n,M,&n,N,&n,w,work,&lwork,&info);

	std::cout << "vertices considered: " << neigs.size() << "\n";
	std::cout << "info=" << info << " (0 is OK)\n";
	std::cout << "work(1)=" << work[0] << " (should be below 512)\n";
	
	//if some error, report it to the caller
	if (info != 0) return info;

	//M is now matrix of eigenvectors
	//it should hold (according to Lapack docs):
	//Z^T N Z = I where Z is one eigenvector, I is identity matrix
	//
	//w holds eigenvalues in ascending order

	//our result c[1]...c[9] is the eigenvector
	//corresponding to the smallest non-negative eigenvalue, so the j-th eigenvector
	int j=0;
	while (j < 9 && w[j] < 0.f) ++j;
	//have we found some non-negative eigenvalue?
	if (j == 9) return(-9999);

	//also:
	//the last missing coefficient c[0] we will determine by submitting
	//the given input vertex to the algebraic expresion of the surface
	//(given with coeffs) and equating it to zero:
	coeffs[0]=0.f;
	for (int i=0; i < 9; ++i)
	{
		coeffs[i+1]=M[j*9 +i];          //copy eigenvector
		coeffs[0]-=coeffs[i+1]*V[i+1];  //determine c[0]
	}

	std::cout << "w(j)=" << w[j] << ", j=" << j << "\n";

*/

	return 0;
}


bool ActiveMesh::GetPointOnQuadricSurface(const float x,const float y,
                                          float &z1, float &z2,
                                          const float (&coeffs)[10])
{
	const float a=coeffs[9];
	const float b=coeffs[3] +coeffs[5]*x +coeffs[6]*y;
	const float c=coeffs[0] +coeffs[1]*x +coeffs[2]*y
	             +coeffs[4]*x*y +coeffs[7]*x*x +coeffs[8]*y*y;

	const float sqArg=b*b - 4*a*c;
	if (sqArg < 0.f) return false;
	if (a == 0.f) return false;

	z1=(-b + sqrtf(sqArg)) / (2.f*a);
	z2=(-b - sqrtf(sqArg)) / (2.f*a);

	return true;
}


float ActiveMesh::GetClosestPointOnQuadricSurface(Vector3F& point,
                                                  const float (&coeffs)[10])
{
	//backup original input coordinate
	const float x=point.x;
	const float y=point.y;
	const float z=point.z;
	float tmp1,tmp2;

	//list of possible coordinates
	std::vector<Vector3F> pointAdepts;

	//took a pair of coordinates, calculate the third one
	//and make it an adept...
	if (GetPointOnQuadricSurface(x,y,tmp1,tmp2,coeffs))
	{
		pointAdepts.push_back(Vector3F(x,y,tmp1));
		pointAdepts.push_back(Vector3F(x,y,tmp2));
	}
	if (GetPointOnQuadricSurface(x,z,tmp1,tmp2,coeffs))
	{
		pointAdepts.push_back(Vector3F(x,tmp1,z));
		pointAdepts.push_back(Vector3F(x,tmp2,z));
	}
	if (GetPointOnQuadricSurface(y,z,tmp1,tmp2,coeffs))
	{
		pointAdepts.push_back(Vector3F(tmp1,y,z));
		pointAdepts.push_back(Vector3F(tmp2,y,z));
	}

	//are we doomed?
	if (pointAdepts.size() == 0)
		return (-999999.f);

	//find the closest
	int closestIndex=ChooseClosestPoint(pointAdepts,point);
	
	//calc distance to it
	point-=pointAdepts[closestIndex];
	tmp1=point.Len();

	//adjust the input/output point
	point=pointAdepts[closestIndex];

	return (tmp1);
}


int ActiveMesh::ChooseClosestPoint(const std::vector<Vector3F>& points,
                                   const Vector3F& point)
{
	int minIndex=-1;
	float minSqDist=9999999999999.f;

	Vector3F p;
	for (unsigned int i=0; i < points.size(); ++i)
	{
		p=point;
		p-=points[i];
		if (p.LenQ() < minSqDist)
		{
			minIndex=i;
			minSqDist=p.LenQ();
		}
	}

	return minIndex;
}


void ActiveMesh::InitDots(const i3d::Image3d<i3d::GRAY16>& mask)
{
	i3d::Image3d<float> dt;
	i3d::GrayToFloat(mask,dt);

	i3d::EDM(dt,0,100.f,false);
#ifdef SAVE_INTERMEDIATE_IMAGES
	dt.SaveImage("1_DTAlone.ics");
#endif

	i3d::Image3d<float> perlinInner,perlinOutside;
	perlinInner.CopyMetaData(mask);
	DoPerlin3D(perlinInner,5.0,0.8*1.5,0.7*1.5,6);
#ifdef SAVE_INTERMEDIATE_IMAGES
	//perlinInner.SaveImage("2_PerlinAlone_Inner.ics");
#endif

	perlinOutside.CopyMetaData(mask);
	DoPerlin3D(perlinOutside,1.8,0.8*1.5,0.7*1.5,6);
#ifdef SAVE_INTERMEDIATE_IMAGES
	//perlinOutside.SaveImage("2_PerlinAlone_Outside.ics");
#endif

	i3d::Image3d<i3d::GRAY16> erroded;
	i3d::ErosionO(mask,erroded,1);

	//initial object intensity levels
	float* p=dt.GetFirstVoxelAddr();
	float* const pL=p+dt.GetImageSize();
	const float* pI=perlinInner.GetFirstVoxelAddr();
	const float* pO=perlinOutside.GetFirstVoxelAddr();
	const i3d::GRAY16* er=erroded.GetFirstVoxelAddr();
	while (p != pL)
	{
		//are we within the mask?
		if (*p > 0.f)
		{
			//close to the surface?
			if (*p < 0.3f || *er == 0) *p=2000.f + 5000.f*(*pO); //corona
			else                       *p=600.f  + 600.f*(*pI);  //inside
			if (*er == 0) *p=2000.f; //std::max(*p,2000.f);                //corona

			if (*p < 0.f) *p=0.f;
		}
		++p; ++pI; ++pO; ++er;
	}
#ifdef SAVE_INTERMEDIATE_IMAGES
	dt.SaveImage("3_texture.ics");
#endif

	perlinInner.DisposeData();
	perlinOutside.DisposeData();
	erroded.DisposeData();

	//now, read the "molecules"
	dots.clear();
	dots.reserve(1<<23);

	//time savers...
	const float xRes=mask.GetResolution().GetX();
	const float yRes=mask.GetResolution().GetY();
	const float zRes=mask.GetResolution().GetZ();
	const float xOff=mask.GetOffset().x;
	const float yOff=mask.GetOffset().y;
	const float zOff=mask.GetOffset().z;

	p=dt.GetFirstVoxelAddr();
	for (int z=0; z < (signed)dt.GetSizeZ(); ++z)
	for (int y=0; y < (signed)dt.GetSizeY(); ++y)
	for (int x=0; x < (signed)dt.GetSizeX(); ++x, ++p)
	for (int v=0; v < *p; v+=50)
	{
		//convert px coords into um
		const float X=(float)x/xRes + xOff;
		const float Y=(float)y/yRes + yOff;
		const float Z=(float)z/zRes + zOff;

		dots.push_back(Vector3F(X,Y,Z));
		/*
		if (dots.size() == dots.capacity())
		{
			std::cout << "reserving more at position: " << x << "," << y << "," << z << "\n";
			dots.reserve(dots.size()+(1<<20));
		}
		*/
	}
//#ifdef SAVE_INTERMEDIATE_IMAGES
	std::cout << "intiated " << dots.size() << " fl. molecules (capacity is for "
	          << dots.capacity() << ")\n";
//#endif
}


void ActiveMesh::BrownDots(const i3d::Image3d<i3d::GRAY16>& mask)
{
	//TODO REMOVE ME
	for (size_t i=0; i < dots.size(); ++i)
		dots[i].x+=1.0f;
}


template <class VT>
VT GetPixel(i3d::Image3d<VT> const &img,const float x,const float y,const float z)
{
	/*
	//nearest neighbor:
	int X=static_cast<int>(roundf(x));
	int Y=static_cast<int>(roundf(y));
	int Z=static_cast<int>(roundf(z));

	if (img.Include(X,Y,Z)) return(img.GetVoxel(X,Y,Z));
	else return(0);
	*/

	//nearest not-greater integer coordinate, "o" in the picture in docs
	//X,Y,Z will be coordinate of the voxel no. 2
	const int X=static_cast<int>(floorf(x));
	const int Y=static_cast<int>(floorf(y));
	const int Z=static_cast<int>(floorf(z));
	//now we can write only to pixels at [X or X+1,Y or Y+1,Z or Z+1]

	//quit if too far from the "left" borders of the image
	//as we wouldn't be able to draw into the image anyway
	if ((X < -1) || (Y < -1) || (Z < -1)) return (0);

	//residual fraction of the input coordinate
	const float Xfrac=x - static_cast<float>(X);
	const float Yfrac=y - static_cast<float>(Y);
	const float Zfrac=z - static_cast<float>(Z);

	//the weights
	float A=0.0f,B=0.0f,C=0.0f,D=0.0f; //for 2D

	//x axis:
	A=D=Xfrac;
	B=C=1.0f - Xfrac;

	//y axis:
	A*=1.0f - Yfrac;
	B*=1.0f - Yfrac;
	C*=Yfrac;
	D*=Yfrac;

	//z axis:
	float A_=A,B_=B,C_=C,D_=D;
	A*=1.0f - Zfrac;
	B*=1.0f - Zfrac;
	C*=1.0f - Zfrac;
	D*=1.0f - Zfrac;
	A_*=Zfrac;
	B_*=Zfrac;
	C_*=Zfrac;
	D_*=Zfrac;

	//portions of the value in a bit more organized form, w[z][y][x]
	const float w[2][2][2]={{{B ,A },{C ,D }},
				{{B_,A_},{C_,D_}}};

	//the return value
	float v=0;

	//reading from the input image,
	//for (int zi=0; zi < 2; ++zi) if (Z+zi < (signed)img.GetSizeZ()) { //shortcut for 2D cases to avoid some computations...
	for (int zi=0; zi < 2; ++zi)
	  for (int yi=0; yi < 2; ++yi)
	    for (int xi=0; xi < 2; ++xi)
		if (img.Include(X+xi,Y+yi,Z+zi)) {
			//if we got here then we can safely change coordinate types
			v+=(float)img.GetVoxel((size_t)X+xi,(size_t)Y+yi,(size_t)Z+zi) * w[zi][yi][xi];
		}
	//}

	return ( static_cast<VT>(v) );
}

void ActiveMesh::FFDots(const i3d::Image3d<i3d::GRAY16>& mask,
                        const FlowField<float> &FF)
{
	//TODO: tests: FF consistency, same size as mask?

	//time savers...
	const float xRes=mask.GetResolution().GetX();
	const float yRes=mask.GetResolution().GetY();
	const float zRes=mask.GetResolution().GetZ();
	const float xOff=mask.GetOffset().x;
	const float yOff=mask.GetOffset().y;
	const float zOff=mask.GetOffset().z;

	//apply FF on the this->dots (no boundary checking)
	for (size_t i=0; i < dots.size(); ++i)
	{
		//turn micron position into pixel one
		const float X=(dots[i].x -xOff) *xRes;
		const float Y=(dots[i].y -yOff) *yRes;
		const float Z=(dots[i].z -zOff) *zRes;

		//note: GetPixel() returns 0 in case we ask for value outside the image
		//TODO: check against mask
		dots[i].x += GetPixel(*FF.x, X,Y,Z);
		dots[i].y += GetPixel(*FF.y, X,Y,Z);
		dots[i].z += GetPixel(*FF.z, X,Y,Z);
	}
}


void ActiveMesh::RenderDots(const i3d::Image3d<i3d::GRAY16>& mask,
                            i3d::Image3d<i3d::GRAY16>& texture)
{
	texture.CopyMetaData(mask);
	texture.GetVoxelData()=0;

	//time savers...
	const float xRes=texture.GetResolution().GetX();
	const float yRes=texture.GetResolution().GetY();
	const float zRes=texture.GetResolution().GetZ();
	const float xOff=texture.GetOffset().x;
	const float yOff=texture.GetOffset().y;
	const float zOff=texture.GetOffset().z;

	//time-savers for boundary checking...
	const int maxX=(int)texture.GetSizeX()-1;
	const int maxY=(int)texture.GetSizeY()-1;
	const int maxZ=(int)texture.GetSizeZ()-1;

	//time-savers for accessing neigbors...
	const size_t xLine=texture.GetSizeX();
	const size_t Slice=texture.GetSizeY() *xLine;
	i3d::GRAY16* const T=texture.GetFirstVoxelAddr();

	//render the points into the texture image
	for (size_t i=0; i < dots.size(); ++i)
	{
		const int x=(int)roundf( (dots[i].x-xOff) *xRes);
		const int y=(int)roundf( (dots[i].y-yOff) *yRes);
		const int z=(int)roundf( (dots[i].z-zOff) *zRes);

		if ((x > 0) && (y > 0) && (z > 0)
		   && (x < maxX) && (y < maxY) && (z < maxZ)) T[z*Slice +y*xLine +x]+=i3d::GRAY16(50);
	}
}

void ActiveMesh::PhaseII(const i3d::Image3d<i3d::GRAY16>& texture,
	                      i3d::Image3d<float>& intermediate)
{
	i3d::GrayToFloat(texture,intermediate);

	i3d::GaussIIR(intermediate,3.f,3.f,2.5f);
#ifdef SAVE_INTERMEDIATE_IMAGES
	intermediate.SaveImage("4_texture_filtered.ics");
#endif
}

void ActiveMesh::PhaseIII(i3d::Image3d<float>& intermediate,
	                       i3d::Image3d<i3d::GRAY16>& texture)
{
#ifdef SAVE_INTERMEDIATE_IMAGES
	intermediate.SaveImage("5_texture_filtered_resampled.ics");
#endif
	float* p=intermediate.GetFirstVoxelAddr();
	for (int z=0; z < (signed)intermediate.GetSizeZ(); ++z)
	for (int y=0; y < (signed)intermediate.GetSizeY(); ++y)
	for (int x=0; x < (signed)intermediate.GetSizeX(); ++x, ++p)
	{
		//background signal:
		float distSq=((float)x-110.f)*((float)x-110.f) + ((float)y-110.f)*((float)y-110.f);
		*p+=150.f*expf(-0.5f * distSq / 2500.f);

		//uncertainty in the number of incoming photons
		const float noiseMean = sqrtf(*p), // from statistics: shot noise = sqrt(signal) 
				      noiseVar = noiseMean;  // for Poisson distribution E(X) = D(X)

		*p+=40.f*((float)GetRandomPoisson(noiseMean) - noiseVar);

		//constants are parameters of Andor iXon camera provided from vendor:
		//photon shot noise: dark current
		*p+=(float)GetRandomPoisson(0.06f);

		//read-out noise:
		//  variance up to 25.f (old camera on ILBIT)
		//  variance about 1.f (for the new camera on ILBIT)
		*p+=GetRandomGauss(480.f,20.f);
		//*p+=530.f;
	}
#ifdef SAVE_INTERMEDIATE_IMAGES
	intermediate.SaveImage("6_texture_filtered_resampled_finalized.ics");
#endif

	//obtain final GRAY16 image
	i3d::FloatToGrayNoWeight(intermediate,texture);
}


void ActiveMesh::ConstructFF(FlowField<float> &FF)
{
	//erase the flow field
	//iterate:
	//  put displacement vectors
	//  smooth

	//tests: TODO
	//FF must be consistent
	//oldVolPos and newVolPos must be of the same length

	//erase
	FF.x->GetVoxelData()=0;
	FF.y->GetVoxelData()=0;
	FF.z->GetVoxelData()=0;

	//time savers...
	const float xRes=FF.x->GetResolution().GetX();
	const float yRes=FF.x->GetResolution().GetY();
	const float zRes=FF.x->GetResolution().GetZ();
	const float xOff=FF.x->GetOffset().x;
	const float yOff=FF.x->GetOffset().y;
	const float zOff=FF.x->GetOffset().z;

	//time-savers for boundary checking...
	const int maxX=(int)FF.x->GetSizeX()-1;
	const int maxY=(int)FF.x->GetSizeY()-1;
	const int maxZ=(int)FF.x->GetSizeZ()-1;

	//time-savers for accessing neigbors...
	const size_t xLine=FF.x->GetSizeX();
	const size_t Slice=FF.x->GetSizeY() *xLine;
	float* const ffx=FF.x->GetFirstVoxelAddr();
	float* const ffy=FF.y->GetFirstVoxelAddr();
	float* const ffz=FF.z->GetFirstVoxelAddr();

	//inject displacements
	for (size_t i=0; i < oldVolPos.size(); ++i)
	{
		const int x=(int)roundf( (oldVolPos[i].x-xOff) *xRes);
		const int y=(int)roundf( (oldVolPos[i].y-yOff) *yRes);
		const int z=(int)roundf( (oldVolPos[i].z-zOff) *zRes);

		const float dx=newVolPos[i].x - oldVolPos[i].x;
		const float dy=newVolPos[i].y - oldVolPos[i].y;
		const float dz=newVolPos[i].z - oldVolPos[i].z;

		if ((x > 0) && (y > 0) && (z > 0)
			&& (x < maxX) && (y < maxY) && (z < maxZ))
		{
			ffx[z*Slice +y*xLine +x]=dx;
			ffy[z*Slice +y*xLine +x]=dy;
			ffz[z*Slice +y*xLine +x]=dz;
		}
	}

	//smooth
	i3d::GaussIIR(*FF.x,15.0f);
	i3d::GaussIIR(*FF.y,15.0f);
	i3d::GaussIIR(*FF.z,15.0f);

	//multiply (to "correct" after normalized smoothing)
	FF.x->GetVoxelData()*=1240.f;
	FF.y->GetVoxelData()*=1240.f;
	FF.z->GetVoxelData()*=1240.f;
}


void ActiveMesh::ConstructFF_T(FlowField<float> &FF)
{
	//erase the flow field
	//add displacement vectors by "rendering" displacement tetrahedra
	//smooth

	//tests: TODO
	//FF must be consistent
	//oldVolPos and newVolPos must be of the same length
	//length of oldVolPos*4 and length of VolID must be the same

	//erase
	FF.x->GetVoxelData()=0;
	FF.y->GetVoxelData()=0;
	FF.z->GetVoxelData()=0;

	i3d::Image3d<i3d::GRAY16> imgCounts;
	imgCounts.CopyMetaData(*FF.x);
	imgCounts.GetVoxelData()=0;

	//time savers...
	const float xRes=FF.x->GetResolution().GetX();
	const float yRes=FF.x->GetResolution().GetY();
	const float zRes=FF.x->GetResolution().GetZ();
	const float xOff=FF.x->GetOffset().x;
	const float yOff=FF.x->GetOffset().y;
	const float zOff=FF.x->GetOffset().z;

	//time-savers for boundary checking...
	const int maxX=(int)FF.x->GetSizeX()-1;
	const int maxY=(int)FF.x->GetSizeY()-1;
	const int maxZ=(int)FF.x->GetSizeZ()-1;

	//time-savers for accessing neigbors...
	const size_t xLine=FF.x->GetSizeX();
	const size_t Slice=FF.x->GetSizeY() *xLine;
	float* const ffx=FF.x->GetFirstVoxelAddr();
	float* const ffy=FF.y->GetFirstVoxelAddr();
	float* const ffz=FF.z->GetFirstVoxelAddr();
	i3d::GRAY16* const ffC=imgCounts.GetFirstVoxelAddr();

	//over all tetrahedra
	for (size_t i=0; i < VolID.size(); i+=4)
	{
		//tetrahedron to drive positions in the FF
		//(positions tetrahedron)
		const Vector3F& v1=oldVolPos[VolID[i+0]];
		const Vector3F& v2=oldVolPos[VolID[i+1]];
		const Vector3F& v3=oldVolPos[VolID[i+2]];
		const Vector3F& v4=oldVolPos[VolID[i+3]];

		//now iterate over the tetrahedra:
		for (float d=0.f; d <= 1.0f; d += 0.04f)
		for (float c=0.f; c <= (1.0f-d); c += 0.04f)
		for (float b=0.f; b <= (1.0f-d-c); b += 0.04f)
		{
			float a=1.0f -b -c -d;

			//float-point coordinate:
			Vector3F tmp;
			tmp =a*v1;
			tmp+=b*v2;
			tmp+=c*v3;
			tmp+=d*v4;

			//pixel (integer) coordinate:
			const int x=(int)roundf( (tmp.x-xOff) *xRes);
			const int y=(int)roundf( (tmp.y-yOff) *yRes);
			const int z=(int)roundf( (tmp.z-zOff) *zRes);

			if ((x > 0) && (y > 0) && (z > 0)
				&& (x < maxX) && (y < maxY) && (z < maxZ))
			{
				//tetrahedron to drive values to place at these positions
				//(values tetrahedron)
				const Vector3F dv1=newVolPos[VolID[i+0]] - v1;
				const Vector3F dv2=newVolPos[VolID[i+1]] - v2;
				const Vector3F dv3=newVolPos[VolID[i+2]] - v3;
				const Vector3F dv4=newVolPos[VolID[i+3]] - v4;

				//value
				tmp =a*dv1;
				tmp+=b*dv2;
				tmp+=c*dv3;
				tmp+=d*dv4;

				ffx[z*Slice +y*xLine +x]+=tmp.x;
				ffy[z*Slice +y*xLine +x]+=tmp.y;
				ffz[z*Slice +y*xLine +x]+=tmp.z;
				ffC[z*Slice +y*xLine +x]+=1;
			}
		}
	}

	//finish the averaging of FF
	for (size_t i=0; i < imgCounts.GetImageSize(); ++i)
	if (*(ffC+i))
	{
		*(ffx+i)/=float(*(ffC+i));
		*(ffy+i)/=float(*(ffC+i));
		*(ffz+i)/=float(*(ffC+i));
	}
	//imgCounts.SaveImage("counts.ics"); //TODO REMOVE

	//a bit more of fine-tunning:
	// make 2 rounds of 1px dilations into a copy image
	// extract the added 2px wide shell
	// widen the shell into another image with 2 rounds of 1px dilations
	// smooth sigma=1px this wider shell
	// mask/narrow the smoothed shell according to the previous shell
	// and combine with the original FF

	i3d::Image3d<i3d::GRAY16> tmp,shell,shellSmooth;
	i3d::Image3d<float> tmpF;
	tmpF=*FF.x;
	tmpF.GetVoxelData()*=1000.f;
	//tmpF.GetVoxelData()+=32768.f;
	tmpF.GetVoxelData()+=30000.f;
	i3d::FloatToGrayNoWeight(tmpF,shell);
	tmpF.SaveImage("shell0_F.ics");
	shell.SaveImage("shell0.ics");

	// make 2 rounds of 1px dilations into a copy image
	i3d::Dilation(shell,tmp,i3d::nb3D_o18);
	i3d::Dilation(tmp,shell,i3d::nb3D_o18);
	shell.SaveImage("shell1.ics");

	// extract the added 2px wide shell
	i3d::Dilation(shell,tmp,i3d::nb3D_o18);
	i3d::Dilation(tmp,shell,i3d::nb3D_o18);
	shell.SaveImage("shell2.ics");

	GrayToFloat(shell,tmpF);
	tmpF.SaveImage("shell2_F.ics");
	tmpF.GetVoxelData()-=30000.f;
	tmpF.GetVoxelData()/=1000.f;
	tmpF.SaveImage("shell2_normalValues_F.ics");

	// widen the shell into another image with 2 rounds of 1px dilations
	// smooth sigma=1px this wider shell
	// mask/narrow the smoothed shell according to the previous shell
	// and combine with the original FF


}