main.cpp

#include <i3d/image3d.h>
#include <i3d/transform.h>

#include <stdlib.h> //for rnd generator seeding
#include <time.h>
#include <unistd.h>
#include <string.h>

#include "params.h"
#include "graphics.h"
#include "../cmath3d/TriangleMesh.h"

//create some shared objects...
ParamsClass params;
void ParamsSetup(void);
//
ActiveMesh mesh;

int LoadNewMesh(int fileNo);
int RenderMesh(int fileNo);

int main(int argc,char **argv)
{
	if (strstr(argv[0],"Gen") != NULL)
	{
		//I'm called as generator/simulator
		ParamsSetup();

		//ID53: mesh #14 looks nice

		//ISBI:
		//for (int i=0; i < 50; ++i)
		for (int i=0; i < 25; ++i)
		{
			LoadNewMesh(i);

			std::cout << "mesh: #" << i << "\n";
			std::cout << "vertices  #: " << mesh.Pos.size() << "\n";
			std::cout << "triangles #: " << mesh.ID.size()/3 << "\n";
			std::cout << "normals   #: " << mesh.norm.size() << "\n";

			std::cout << "filopodia #: " << mesh.Ftree.size() << "\n";
			for (size_t f=0; f < mesh.Ftree.size(); ++f)
			{
				std::cout << f << ": vertices  #: " << mesh.Ftree[f].Pos.size() << "\n";
				std::cout << f << ": triangles #: " << mesh.Ftree[f].ID.size()/3 << "\n";
				std::cout << f << ": normals   #: " << mesh.Ftree[f].norm.size() << "\n";
			}

			RenderMesh(i);
		}

		/*
		initializeGL();
		loopGL();
		closeGL();
		*/

		return (0);
	}

	//I'm called as Compactor or Visualizator
	if (argc != 6 && argc != 7)
	{
		std::cout << "wrong number of parameters, should be:\n";
		std::cout << "in_nuclVol.vtk in_nuclSurf.vtk in_filoA.vtk in_filoB.vtk out_mesh.stl [filoStretcher]\n";
		std::cout << "the last filoStretcher is optional parameter\n";
		return(99);
	}

	float stretcher=1.f;
	if (argc == 7) stretcher=(float)atof(argv[6]);
	std::cout << "using filopodia radius stretcher = " << stretcher << "\n";

	//int retval=mesh.ImportVTK_Volumetric(argv[1]);
	int retval=mesh.ImportVTK(argv[2]);
	if (retval) 
	{
		std::cout << "some error reading nucleus: " << retval << "\n";
		return(1);
	}

	retval=mesh.ImportVTK_Ftree(argv[3],0,stretcher);
	if (retval) 
	{
		std::cout << "some error reading f-tree: " << retval << "\n";
		return(2);
	}

	retval=mesh.ImportVTK_Ftree(argv[4],1,stretcher);
	if (retval) 
	{
		std::cout << "some error reading f-tree: " << retval << "\n";
		return(2);
	}

	std::cout << argv[0] << std::endl;
	std::cout << argv[0][9] <<  std::endl;
	if (strstr(argv[0],"Comp") != NULL)
		//I am called as Compactor
		mesh.ExportSTL(argv[5]);
	else
	{
		//I am visualizer
		ParamsSetup();

		std::cout << "vertices  #: " << mesh.Pos.size() << "\n";
		std::cout << "triangles #: " << mesh.ID.size()/3 << "\n";
		std::cout << "normals   #: " << mesh.norm.size() << "\n";
		mesh.CenterMesh(params.sceneCentre);

		initializeGL();
		loopGL();
		closeGL();
	}
	return(0);
}


void ParamsSetup(void)
{
	//set up the environment
   params.sceneOffset=Vector3d<float>(0.f);
   params.sceneSize=Vector3d<float>(27.5f,27.5f,22.0f); //for sample cell

	params.sceneCentre=params.sceneSize;
	params.sceneCentre/=2.0f;
	params.sceneCentre+=params.sceneOffset;

	params.sceneOuterBorder=Vector3d<float>(2.f);
	params.sceneBorderColour.r=0.7f;
	params.sceneBorderColour.g=0.7f;
	params.sceneBorderColour.b=0.0f;

	params.inputCellsFilename="cells/cell%d.txt";
	params.numberOfAgents=2;

	params.friendshipDuration=10.f;
	params.maxCellSpeed=0.20f;

	//adapted for cell cycle length of 24 hours
	params.cellCycleLength=14.5*60; //[min]

	params.currTime=0.0f; //all three are in units of minutes
	params.incrTime=1.0f;
	params.stopTime=50.0f;

	params.imgSizeX=500; //pixels
	params.imgSizeY=500;
	params.imgResX=1.0f; //pixels per micrometer
	params.imgResY=1.0f;

	params.imgOutlineFilename="Outline%05d.tif";
	params.imgPhantomFilename="Phantom%05d.tif";
	params.imgMaskFilename="Mask%05d.tif";
	params.imgFluoFilename="FluoImg%05d.tif";
	params.imgPhCFilename="PhCImg%05d.tif";
	params.tracksFilename="tracks.txt";


	//re-seed the standard rnd generator
	//(which is used for Perlin texture)
	unsigned int seed=-1 * (int)time(NULL) * (int)getpid();
	srand(seed);
}


int LoadNewMesh(int fileNo)
{
	char VTK_nucleusSurf[1024];
	char VTK_nucleus[1024];
	char VTK_fTreeA[1024];
	char VTK_fTreeB[1024];

	sprintf(VTK_nucleusSurf,"../sample_sequence/ID64/ID64_t%03d_nucleusSurf.vtk",fileNo*2);
	sprintf(VTK_nucleus,    "../sample_sequence/ID64/ID64_t%03d_nucleus.vtk",fileNo*2);
	sprintf(VTK_fTreeA,     "../sample_sequence/ID64/ID64_t%03d_fTree01.vtk",fileNo*2);
	sprintf(VTK_fTreeB,     "../sample_sequence/ID64/ID64_t%03d_fTree02.vtk",fileNo*2);

	//first, obtain a fresh volumetric mesh of the cell body
	//(before it gets overwritten in the following code)
	//
	//keep the volumetric configuration of the first load
	bool saveAlsoTetrahedra=(mesh.GetVolIDsize() == 0);
	int retval=mesh.ImportVTK_Volumetric(VTK_nucleus,saveAlsoTetrahedra);
	if (retval) 
	{
		std::cout << "some error reading nucleus: " << retval << "\n";
		return(1);
	}
	//position the mesh somewhat inside the scene
	mesh.ScaleMesh(Vector3F(180.f,180.f,250.f));
	mesh.TranslateMesh(params.sceneCentre);
	//and rotate both arrays
	mesh.RotateVolPos();

	//continue reading surface meshes and filopodia
	retval=mesh.ImportVTK(VTK_nucleusSurf);
	if (retval) 
	{
		std::cout << "some error reading nucleus: " << retval << "\n";
		return(1);
	}

	mesh.dots_filo.clear();
	mesh.dots_filo.reserve(300000);
	retval=mesh.ImportVTK_Ftree(VTK_fTreeA,0);
	if (retval) 
	{
		std::cout << "some error reading f-tree: " << retval << "\n";
		return(2);
	}

	retval=mesh.ImportVTK_Ftree(VTK_fTreeB,1);
	if (retval) 
	{
		std::cout << "some error reading f-tree: " << retval << "\n";
		return(2);
	}

	std::cout << "loaded mesh #" << fileNo*2 << "\n";

	//position the mesh somewhat inside the scene
	mesh.ScaleMesh(Vector3F(180.f,180.f,250.f));
	mesh.TranslateMesh(params.sceneCentre);

	return(0);
}


#define DO_PHASE_II_AND_III
#define DO_ANISOTROPIC_OUTPUT
int RenderMesh(int fileNo)
{
	//setup isotropic(!) image for sample cell
	i3d::Image3d<i3d::GRAY16> mask;
	mask.SetOffset(i3d::Offset(0.0,0.0,0.0));
	mask.SetResolution(i3d::Resolution(8.0,8.0,8.0));
	mask.MakeRoom(220,220,176);
	mask.GetVoxelData()=0;

	//renders the isotropic mask
	std::cout << "rendering mask #" << fileNo << "\n";
	mesh.RenderMask_body(mask);
	mesh.RenderMask_filo(mask,0);
	mesh.RenderMask_filo(mask,1);

	//updates the texture fl. points
	if (mesh.DotsFirstRun())
	{
		std::cout << "rendering texture first run\n";
		mesh.InitDots_body(mask);
		//mesh.InitDots_filo(mask); -- already re-created with ImportVTK_Ftree()
		std::cout << "filo dots size=" << mesh.dots_filo.size() << "\n";
	}
	else
	{
		std::cout << "rendering texture iteration\n";
		//mesh.BrownDots(mask);

		//obtain fresh flow field, so that point positions can be updated
		std::cout << "oldVol.size=" << mesh.oldVolPos.size() << ", newVol.size=" << mesh.newVolPos.size() << "\n";

		FlowField<float> FF;
		FF.x=new i3d::Image3d<float>;
		FF.y=new i3d::Image3d<float>;
		FF.z=new i3d::Image3d<float>;
		FF.x->CopyMetaData(mask);
		FF.y->CopyMetaData(mask);
		FF.z->CopyMetaData(mask);

		//create a FF by "diffusing displacement seeds/vertices"
		mesh.ConstructFF(FF);
/*
		char n[1024];
		sprintf(n,"old_vx_%03d.ics",fileNo);
		FF.x->SaveImage(n);
		sprintf(n,"old_vy_%03d.ics",fileNo);
		FF.y->SaveImage(n);
		sprintf(n,"old_vz_%03d.ics",fileNo);
		FF.z->SaveImage(n);

		//create a FF by "rendering displacement triangles"
		mesh.ConstructFF_T(FF);
		sprintf(n,"new_vx_%03d.ics",fileNo);
		FF.x->SaveImage(n);
		sprintf(n,"new_vy_%03d.ics",fileNo);
		FF.y->SaveImage(n);
		sprintf(n,"new_vz_%03d.ics",fileNo);
		FF.z->SaveImage(n);
*/
		//apply the FF on the fl. dots
		mesh.FFDots(mask,FF);

		//FF destructor will delete dyn. allocated images...
	}

	//renders the isotropic texture
	std::cout << "rendering texture I   #" << fileNo << "\n";
	i3d::Image3d<i3d::GRAY16> texture;
	mesh.RenderDots(mask,texture);
	//
	//or always one shot:
	//i3d::Image3d<i3d::GRAY16> texture;
	//mesh.RenderOneTimeTexture(mask,texture);

	char fileName[1024];
	sprintf(fileName,"phantom_t%03d.tif",fileNo);
	texture.SaveImage(fileName);

#ifdef DO_PHASE_II_AND_III
	std::cout << "rendering texture II  #" << fileNo << "\n";
	i3d::Image3d<float> intermediate;
	mesh.PhaseII(texture,intermediate);
#ifdef DO_ANISOTROPIC_OUTPUT
	//downsample now:
	float factor[3]={1.f,1.f,0.125f};
	i3d::Image3d<float>* tmp=NULL;
	i3d::lanczos_resample(&intermediate,tmp,factor,2);
	intermediate=*tmp;
	delete tmp;
#endif
	std::cout << "rendering texture III #" << fileNo << "\n";
	mesh.PhaseIII(intermediate,texture);
#endif

	std::cout << "exporting texture #" << fileNo << "\n";
	sprintf(fileName,"texture_t%03d.tif",fileNo);
	texture.SaveImage(fileName);

	std::cout << "exporting mask #" << fileNo << "\n\n";
	sprintf(fileName,"mask_t%03d.tif",fileNo);
#ifdef DO_ANISOTROPIC_OUTPUT
	//resample also the mask before exporting
	i3d::Resample(mask,texture,
	              mask.GetSizeX(),mask.GetSizeY(),mask.GetSizeZ()/8);
	texture.SaveImage(fileName);
#else
	mask.SaveImage(fileName);
#endif

	return(0);
}