Merge branch 'ISBI2016_TextureRenderer' of...

#undef max

//some debug-code enabling triggers

//#define SAVE_INTERMEDIATE_IMAGES

#define SAVE_INTERMEDIATE_IMAGES

#define TEXTURE_QUANTIZATION	64

#define TEXTURE_QUANTIZATION 8	

#define SQR(x) ((x)*(x))

#define SKIP GetRandomUniform(0.0, 1.0)

//'multiple' should be ideally 10^desired_decimal_accuracy

int inline RoundTo(const float val, const float multiple=1000.f)

{

	i3d::Image3d<float> perlinInner,perlinOutside;

	perlinInner.CopyMetaData(mask);

//	DoPerlin3D(perlinInner,5.0,0.8*1.5,0.7*1.5,6);

	DoPerlin3D(perlinInner,1.0,0.8*1.3,0.7*1.5,6);

	DoPerlin3D(perlinInner,1.0,0.3,1.0,6);

	Stretch(perlinInner, 0.0f, 1.0f, 4.0f);

#ifdef SAVE_INTERMEDIATE_IMAGES

	//perlinInner.SaveImage("2_PerlinAlone_Inner.ics");

	perlinInner.SaveImage("2_PerlinAlone_Inner.ics");

#endif

	perlinOutside.CopyMetaData(mask);

	DoPerlin3D(perlinOutside,1.4,0.8*1.1,0.7*1.5,6);

	DoPerlin3D(perlinOutside,2.5,0.2,1.0,6);

	Stretch(perlinOutside, 0.0f, 1.0f, 2.0f);

#ifdef SAVE_INTERMEDIATE_IMAGES

	//perlinOutside.SaveImage("2_PerlinAlone_Outside.ics");

	perlinOutside.SaveImage("2_PerlinAlone_Outside.ics");

#endif

	i3d::Image3d<i3d::GRAY16> erroded;

	const float* pI=perlinInner.GetFirstVoxelAddr();

	const float* pO=perlinOutside.GetFirstVoxelAddr();

	const i3d::GRAY16* er=erroded.GetFirstVoxelAddr();

	// intensity values

	const float innerBg = 0.0;

	const float innerPerlinStrength = 20.0f;

	const float coronaBg = 1.0f;

	const float coronaPerlinStrength = 330.0f;

	// distances (in microns)

	const float coronaWidth = 0.1f;

	const float interfaceCoronaInterior = 0.05f;

	while (p != pL)

	{

		//are we within the mask?

		if (*p > 0.f)

		{

			//close to the surface?

			if (*p < 0.2f || *er == 0) 

				 *p=1000.f + 5000.f*(*pO); //corona

			else if (*p < 0.4f)

			if (*p < coronaWidth || *er == 0) 

				 *p=coronaBg + coronaPerlinStrength*(*pO); // corona

			else if (*p < (coronaWidth + interfaceCoronaInterior))

			{

				 float factor = (*p - 0.2f)/(0.2f);

				 *p=(300.f  + 600.f*(*pI))*factor +  (1000.f + 5000.f*(*pO))*(1.0f-factor);

				 // linear interpolation factor (between corona and interior)

				 float factor = (*p - coronaWidth)/(interfaceCoronaInterior); 

				 *p=(innerBg + innerPerlinStrength*(*pI))*factor +  

				 	 (coronaBg + coronaPerlinStrength*(*pO))*(1.0f-factor);

			}

			else

				 *p=300.f  + 600.f*(*pI);  //inside

				 *p=innerBg  + innerPerlinStrength*(*pI);  //inside

			if (*er == 0) *p=1500.f; //corona

			if (*er == 0) *p=coronaBg; //corona envelope (was 1500)

			if (*p < 0.f) *p=0.f;

		}

	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+=TEXTURE_QUANTIZATION)

   for (int v=0; v < floorf(*p); v+=TEXTURE_QUANTIZATION)

	{

		//convert px coords into um, use centres of voxels

		const float X=(float(x)+0.5f)/xRes + xOff;

void ActiveMesh::InitDots_filo(const i3d::Image3d<i3d::GRAY16>& mask,const unsigned int idx)

{

   // get the resolution of image

	i3d::Vector3d<float> res = mask.GetResolution().GetRes();

	float max_res = std::max(res.x, std::max(res.y, res.z)); // ... pixels/micron

	const int density = 10; // ... how may dots should be pressed between two neighbouring pixel positions

	//syntactic short-cut to the proper filopodium geometry data

	mesh_t& filoMesh=Ftree[idx];

		const Vector3F& fP=filoMesh.fPoints[filoMesh.segFromPoint[i]];

		const Vector3F& tP=filoMesh.fPoints[filoMesh.segToPoint[i]];

		const float length=(tP-fP).Len(); //full segment has length = 0.005

		const float length=(tP-fP).Len(); // segment length in microns

		// discrete stepping

		const int S = (int) ceil(length * max_res * density);

		/// OBSOLETE - with magic factor

		//stepping along the axis 

		const int S=(int)ceil(5.f*length/0.005); //"fraction factor"

		// const int S=(int)ceil(5.f*length/0.005); //"fraction factor"

		for (int s=0; s < S; ++s)

		{

			Vector3F p =fS*fP;

			         p+=tS*tP;

			// OBSOLETE - magic factor

			//insert fl. dots

			float dotsCount=400000.f*r;

//			float dotsCount=400000.f*r;

			// the number of dots should correspond the width (radius) of filopodium

			// multiplied by given density. The density should be a large number as the

			// radius is typically very small

			const int dotDensity = 5;

		   float dotsCount = dotDensity*r;

			// prevent data overflow (5 - multiplicity of addition)

			int upper_limit = std::numeric_limits<i3d::GRAY16>::max() / (5 * TEXTURE_QUANTIZATION);

			Mul(mainDirection, vector2, vector4);

			// make them as small as the size of the radius in the given position

			vector1 *= r / vector1.Len();

			vector2 *= r / vector2.Len();

			vector3 *= r / vector3.Len();

			vector4 *= r / vector4.Len();

			vector1 /= vector1.Len();

			vector2 /= vector2.Len();

			vector3 /= vector3.Len();

			vector4 /= vector4.Len();

		   int radiusInPixels = (int) ceil(r*max_res);

			for (int q=0; q < qM; ++q)

			{

				dots_filo.push_back(p);

				//dots_filo.push_back(p+Vector3F(0.0f,0.0f,+100*r));

				//dots_filo.push_back(p+Vector3F(0.0f,0.0f,-100*r));

				dots_filo.push_back(p+vector1);

				dots_filo.push_back(p+vector2);

				dots_filo.push_back(p+vector3);

				dots_filo.push_back(p+vector4);

				// The following cycle guarantees, that the dots will fill the filopodium

				// smoothly - from the center to the outer envelope. 

				for (int step=0; step<=radiusInPixels; step++)

				{

				   // how much the given vector should enlarge

				   float scale = (r*(step+1.0f))/(radiusInPixels+1.0f);

					dots_filo.push_back(p+scale*vector1);

					dots_filo.push_back(p+scale*vector2);

					dots_filo.push_back(p+scale*vector3);

					dots_filo.push_back(p+scale*vector4);

				}

			}

		}

	}

			{

				const int X=int(x-0.5f+float(ix)/float(xn));

				++T[X];

				//DAVID: tady nepricitat +1, ale +GetRandomGauss(40.f,40.f/2.5f)

				//DAVID: ( za predpokladu, ze GetRandomGauss(mean,variance) )

				//DAVID: variance nastavena zamerne s malym prekryvem

				//DAVID: pozor na pretekani intervalu GRAY16

				//std::cout << "populating at [" << X << "," << Y << "," << Z << "]\n";

			}

		 for (int y=0; y < (signed)blurred_texture.GetSizeY(); ++y)

			  for (int x=0; x < (signed)blurred_texture.GetSizeX(); ++x, ++p)

			  {

				// background signal:

			   // max distance

				float maxDist = SQR((float)xC) + SQR((float)yC);

				// background signal (inverted parabola)

				float distSq = 

						  -(SQR((float)(x-xC)) + SQR((float)(y-yC)))/24200.f + 1.f;

						  -(SQR((float)(x-xC)) + SQR((float)(y-yC)))/maxDist + 1.f;

				*p *= distSq;

			  }

#ifdef SAVE_INTERMEDIATE_IMAGES

void ActiveMesh::PhaseIII(i3d::Image3d<float>& blurred,

	                       i3d::Image3d<i3d::GRAY16>& texture)

{

	// NONSPECIFIC BACKGROUND

   i3d::Image3d<float> bgImg;

	bgImg.CopyMetaData(blurred);

	DoPerlin3D(bgImg,10.0,7.0,1.0,10); // very smooth a wide coherent noise

	// reset the maximum and minimin intesity levels of the background

	// to the expected values

	const float oldMaxI = bgImg.GetVoxelData().max();

	const float oldMinI = bgImg.GetVoxelData().min();

	const float newMaxI = 2.2f;

	const float newMinI = 2.0f;

	const float scale = (newMaxI-newMinI)/(oldMaxI-oldMinI);

	for (size_t i=0; i<bgImg.GetImageSize(); i++)

	{

		float value = bgImg.GetVoxel(i);

		value = scale*(value - oldMinI) + newMinI;

	   bgImg.SetVoxel(i, value);

	}

   // given gain of EMCCD camera

   const int EMCCDgain = 1000;

	float* p=blurred.GetFirstVoxelAddr();

	for (size_t i=0; i<blurred.GetImageSize(); i++, p++)

	{

		// shift the signal (simulates non-ideal black background)

		// ?reflection of medium?

		*p+=10.0f;

		*p+=bgImg.GetVoxel(i);

		// ENF ... Excess Noise Factor (stochasticity of EMCCD gain)

		// source of information:

		// https://www.qimaging.com/resources/pdfs/emccd_technote.pdf

		const float ENF = GetRandomUniform(1.0f, 1.4f);

		// PHOTON NOISE 

		// uncertainty in the number of incoming photons

		// from statistics: shot noise = sqrt(signal) 

		// from statistics: shot noise mean = sqrt(signal) 

		const float noiseMean = sqrtf(*p);

		// for Poisson distribution E(X) = D(X)

		const float noiseVar = noiseMean; 

		// SKIP avoid strong discretization of intenzity histogram due

		// to the Poisson probabilty function

 		*p += ENF * ((float)GetRandomPoisson(noiseMean) - noiseMean + SKIP);

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

		//DAVID: tady odstranit 40.0f*, cili uz nenasobit

		//DAVID: bude asi nutne preladit vsechny Bulharske konstanty

		// EMCCD GAIN

      // amplification of signal (and inevitably also the noise)

	   *p *= EMCCDgain; 

		// DARK CURRENT

		// constants are parameters of Andor iXon camera provided from vendor:

		*p+=(float)GetRandomPoisson(0.06f);

		*p += ENF*EMCCDgain*(float)GetRandomPoisson(0.06f);

		// BASELINE (camera electron level)

		*p += 400.0f;

		// READOUT NOISE

		// variance up to 25.f (old camera on ILBIT)

		// variance about 1.f (for the new camera on ILBIT)

		*p+=GetRandomGauss(450.f,10.f);

		*p+=GetRandomGauss(0.0f,1.0f);

		// ADC (analogue-digital converter)

		// ADCgain ... how many electrons correcpond one intensity value

		// ADCoffset ... how much intensity levels are globally added to each pixel

		const float ADCgain = 28.0f; 

		*p/=ADCgain;

		const float ADCoffset = 400.0f;

		*p+=ADCoffset;

	}

#ifdef SAVE_INTERMEDIATE_IMAGES

	blurred.SaveImage("6_texture_finalized.ics");

	void PopulateSurfTriangles_Ftree(const unsigned int idx,const bool enforce=false);

	///renders filopodium first aside, then overlays the given mask

	void RenderMask_filo(i3d::Image3d<i3d::GRAY16>& mask,const unsigned int idx,const bool showTriangles=false);

	void RenderMask_filo(i3d::Image3d<i3d::GRAY16>& mask,const unsigned int idx,

		const unsigned short color=100,const bool showTriangles=false);

	/// render directly into mask using volumetric mesh

	void RenderMask_bodyVol(i3d::Image3d<i3d::GRAY16>& mask,const bool showTriangles=false);

	//translating all mesh points, to assure the meshes appear within

	//the produced images

	mesh.meshShift=minBB -Vector3F(1.0f);

	mesh.meshPlaygroundSize=(maxBB-minBB) +Vector3F(2.0f);

	mesh.meshPlaygroundSize=(maxBB-minBB) +Vector3F(5.0f);

	//NB: We also introduce 1um thick border around the mesh

	std::cout << "corners of a detected bounding box (with 1um thick border):\n";

	std::cout << "[" << minBB.x << "," << minBB.y << "," << minBB.z << "] <-> ["