segmentation / ivote3

  #include <iostream>

  #include <string>

  #include <fstream>

  #include <cuda_runtime.h>
  #include <stim/math/vector.h>
  #include <stim/parser/arguments.h>
  #include <stim/parser/filename.h>

  #include <stim/visualization/colormap.h>

  #include <stim/image/image.h>

  #define pi	3.14159
  

  //#define M_PI	3.14159
  //#include <stim/math/circle.h>
  //#include <stim/math/vec3.h>
  //#include <stim/math/plane.h>
  //#include <stim/math/vector.h>
  //#include <stim/visualization/aabb3.h>
  

  
  
  /*void test_3(float* gpu_out, float* gpu_grad, float rmax, float phi, int n, int x, int y, int z);
  
  int main(){
  	
  	int  n=20;
  	float rmax;
  	float phi_deg;
  	float phi;
  	rmax=4;
  	phi_deg = 15;
  	phi = phi_deg * pi/180;
  	int x,y,z;
  	x=y=z=1;
  	unsigned int size = x*y*z*sizeof (float);
  	float* cpu_grad = (float*) malloc(3*size);
  	float* gpu_grad;
  	cudaMalloc(&gpu_grad, 3*size);
  	cpu_grad[0]=1;
  	cpu_grad[1]=0;
  	cpu_grad[2]=-0.5;
  	cudaMemcpy(gpu_grad, cpu_grad, 3*size, cudaMemcpyHostToDevice);
  	float* cpu_out = (float*) malloc(3*size*(n+1));
  	float* gpu_out;
  	cudaMalloc(&gpu_out, 3*size*(n+1));
  	test_3(gpu_out, gpu_grad, rmax, phi, n, x, y, z);
  	cudaMemcpy(cpu_out, gpu_out, 3*size*(n+1), cudaMemcpyDeviceToHost);
  	std::ofstream list("circle_check_cuda1.txt");
  	if (list.is_open()){
  		for (int j=0; j<=n; ++j)
  			list << cpu_out[3*j] << '\t' << cpu_out[3*j +1] << '\t' << cpu_out[3*j + 2] << '\n';
  	}
  	list.close();
  	*/

  /*
  	int main(){
  
  
  		stim::vec3<float> g(-44,-3.4,-0.005);   // form a vec3 variable for the gradient vector
  			stim::vec3<float> g_sph = g.cart2sph();			//convert cartesian coordinate to spherical for the gradient vector
  			int n =36;										//set the number of points to find the boundaries of the conical voting area
  			int xi = 105;
  			int yi = 17;
  			int zi = 23;
  			float xc = 12 * cos(g_sph[1]) * sin(g_sph[2]);			//calculate the center point of the surface of the voting area for the voter
  			float yc = 10 * sin(g_sph[1]) * sin(g_sph[2]) ;
  			float zc = 10 * cos(g_sph[2]) ;
  			float r = sqrt(xc*xc + yc*yc + zc*zc);
  			xc+=xi;
  			yc+=yi;
  			zc+=zi;
  			stim::vec3<float> center(xc,yc,zc);
  			
  			float d = 2 * r * tan(25*pi/180 );		//find the diameter of the conical voting area
  			stim::vec3<float> norm = g.norm();			//compute the normalize gradient vector
  			float step = 360.0/(float) n;
  			stim::circle<float>  cir(center, d, norm);
  			stim::aabb3<int> bb(xi,yi,zi);
  			bb.insert(xc,yc,zc);
  			for(float j = 0; j <360.0; j += step){
  				stim::vec3<float> out = cir.p(j);
  				bb.insert(out[0], out[1], out[2]);
  			}
  
  			bb.trim_low(0,0,0);
  			bb.trim_high(128-1, 128-1, 128-1);
  
  			std::cout<< bb.low[0] << '\t' << bb.low[1] << '\t' << bb.low[2] << '\n';
  			std::cout<< bb.high[0] << '\t' << bb.high[1] << '\t' << bb.high[2] << '\n';
  			std::cin >> n;
  */
  		/*int  n=10;
  		stim::circle<float>  cir;
  		float* c0= (float*) malloc(3*sizeof(float));
  		c0[0] =-4;
  		c0[1]=0;
  		c0[2] = 3;
  		stim::vec3<float> c(c0[0],c0[1],c0[2]);
  		float len = c.len();
  		stim::vec3<float> norm(c0[0]/len,c0[1]/len,c0[2]/len);
  		std::cout<< len << '\n';
  		std::cout<< norm << '\n';
  		cir.center(c);
  		cir.normal(norm);
  		cir.scale(2);
  		stim::vec3<float> out = cir.p(45);
  		std::vector<stim::vec3<float>> out2 = cir.getPoints(n);

  		std::cout<< out << '\n';
  		std::cout <<out[0] << '\t' << out[1] << '\t' << out[2] <<'\n';
  		std::cout<< c << '\n';

  		for (std::vector<stim::vec3<float>>::const_iterator i = out2.begin(); i != out2.end(); ++i)
  		std::cout << *i << '\n';
  		std::ofstream list("circle_check.txt");
  		if (list.is_open()){
  			for (std::vector<stim::vec3<float>>::const_iterator j = out2.begin(); j != out2.end(); ++j)
  				list << *j << '\n';
  		}
  		list.close();
  		std::cin >> n;

  
  }

  */

  void ivote3(float* img, float std[], float anisotropy, float phi, float d_phi, unsigned int r[], int iter, float t, unsigned int conn[], 

  			unsigned int x, unsigned int y, unsigned int z);

  void lmax(float* center, float* vote, float t1, unsigned int conn[], unsigned int x, unsigned int y, unsigned int z);

  
  void invert_data(float* cpuI, unsigned int x, unsigned int y, unsigned int z){
  		for(int ix = 0; ix < x; ix++){
  			for (int iy = 0; iy < y; iy++){
  				for (int iz = 0; iz < z; iz++){
  					int idx = iz * x * y + iy * x + ix;
  					cpuI[idx] = 255 - cpuI[idx];
  				}
  			}
  		}
  	}

  int main(int argc, char** argv){
  

  
  	cudaDeviceProp prop;
  	int count;
  	cudaGetDeviceCount(&count);
  	//printf("cudadevicecount: %i\n", count);
  	for (int i=0; i<count; i++){
  		cudaGetDeviceProperties(&prop, i);
  		printf("current device ID: %d\n", i);
  		printf("device name: %s\n", prop.name);
  		printf("total global mem: %lu\n", prop.totalGlobalMem);

  		printf("shared memory per block: %lu\n", prop.sharedMemPerBlock);

  	}

  	//output advertisement
  	std::cout<<std::endl<<std::endl;
  	std::cout<<"========================================================================="<<std::endl;
  	std::cout<<"Thank you for using the ivote3 segmentation tool!"<<std::endl;
  	std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
  	std::cout<<"Developers: Laila Saadatifard and David Mayerich"<<std::endl;
  	std::cout<<"Source: https://git.stim.ee.uh.edu/segmentation/ivote3"<<std::endl;
  	std::cout<<"========================================================================="<<std::endl<<std::endl;
  
  	stim::arglist args;
  
  #ifdef _WIN32
  	args.set_ansi(false);
  #endif
  
  	//add arduments
  	args.add("help", "prints this help");
  	args.add("x", "size of the dataset along X axis", "positive value");
  	args.add("y", "size of the dataset along Y axis", "positive value");
  	args.add("z", "size of the dataset along Z axis", "positive value");
  	args.add("t", "threshold value for the final result", "positive valu");

  	args.add("invert", "to invert the input data set", "string");

  	args.add("anisotropy", "anisotropy value of the imaging", "1");

  	//parse the command line arguments.
  	args.parse(argc, argv);
  
  	//display the help text if requested
  	if(args["help"].is_set()){				
  		std::cout<<std::endl<<"usage: ivote input_image output_list --option [A B C ...]"<<std::endl;
  		std::cout<<std::endl<<std::endl
  				  << "examples: ivote blue.bmp list.txt "<<std::endl;
  		
  		std::cout<<std::endl<<std::endl;
  		std::cout<<args.str()<<std::endl;
  		exit(1);
  	}
  
  	//if the input and output files aren't specified, throw an error and exit
  	if(args.nargs() < 2){
  		std::cout<<"ERROR: two files must be specified for segmentation, enter ivote --help for options."<<std::endl<<std::endl;
  		exit(1);
  	}
  
  	//get the input image file
  	stim::filename Ifilename(args.arg(0));
  
  	//get the output file name
  	stim::filename OutName(args.arg(1));
  
  	//set the x, y, z.
  	int x = args["x"].as_int();
  	int y = args["y"].as_int();
  	int z = args["z"].as_int();
  
  	//set the threshold.
  	float t = args["t"].as_float();

  	//set the anisotropy
  	float anisotropy =  args["anisotropy"].as_float();

  	unsigned int rmax = 10;

  	unsigned int r[3] = { 12, rmax, rmax};

  	float std = 5;
  	float sigma[3] = { std, std, std};

  	unsigned int nlmax = 5;

  	unsigned int conn[3] = { nlmax, nlmax, nlmax};
  	float phi_deg = 25.0;

  	float phi = phi_deg * pi /180;

  	int iter = 8;
  	float d_phi = phi/(iter+2);

  	
  	std::string filename = Ifilename.str();

  	unsigned int bytes = x*y*z*sizeof(float);
  

  	//allocate space on the cpu for the input data

  	float* cpuI = (float*) malloc(bytes);
  

  	//load the input file into the cpuI

  	std::ifstream nissl(filename, std::ios::in | std::ios::binary);
  	nissl.read((char*)cpuI, bytes);
  	nissl.close();

  	if(args["invert"].is_set())
  		invert_data(cpuI, x, y, z);

  	
  	//write a new file from the cpuI.

  	std::ofstream original("shared2D-v1/inv-128.vol", std::ofstream::out | std::ofstream::binary);

  	original.write((char*)cpuI, bytes);
  	original.close();

  	//allocate space on the cpu for the output result

  	float* cpu_out = (float*) malloc(bytes);

  	// call the ivote function

  	//ivote3(cpu_out, cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z);
  	ivote3(cpuI, sigma, anisotropy, phi, d_phi, r, iter, t, conn, x, y, z);

  	//write the blurred file from the cpuI.

  	std::ofstream fblur("shared2D-v8/vote8.vol", std::ofstream::out | std::ofstream::binary);

  	fblur.write((char*)cpuI, bytes);
  	fblur.close();

  	
  	//stim::image<float>imgrad3;
  	//imgrad3.set_interleaved3(cpu_out, 128,128,128,3);
  	//std::ofstream fgx("syn/gx-128.vol", std::ofstream::out | std::ofstream::binary);
  	//fgx.write((char*)imgrad3.channel(0).data(), bytes);
  	//fgx.close();
  	

  	//write the output file.

  	//for (int t0=0; t0<=5000; t0+=100){
  	//	float t1 = t0;

  		int t0 = t;
  			lmax(cpu_out, cpuI, t, conn, x, y, z);
  			//std::ofstream fo("shared2D-v8/" + OutName.str(), std::ofstream::out | std::ofstream::binary);

  			std::ofstream fo( OutName.str()+std::to_string(t0)+".vol", std::ofstream::out | std::ofstream::binary);

  			fo.write((char*)cpu_out, bytes);
  			fo.close();

  	// creat a file for saving the list centers

  		std::ofstream list(OutName.str()+std::to_string(t0)+".obj");
  		// set the number of detected cells to zero.
  		int nod = 0;
  		if (list.is_open()){
  
  				for (int iz=0; iz<z; iz++){
  					for (int iy=0; iy<y; iy++){
  						for (int ix=0; ix<x; ix++){
  
  							int idx = iz * x * y + iy * x + ix;
  							if (cpu_out[idx]>0){
  								nod++;
  								list << "v" << "\t" << ix << "\t" << iy << "\t"<< iz << "\t" << cpu_out[idx] << '\n' ;
  					
  							}
  						}
  					}
  				}
  				list << "p" << "\t";
  				for (unsigned int i_nod =1 ; i_nod <=nod; i_nod++){
  					list << i_nod << "\t";
  				}
  
  		list.close();
  		}
  
  

  	//}

  		cudaDeviceReset();       

  }