validation / pointmets

#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <fstream>
#include <stim/math/vector.h>
#include <stim/parser/arguments.h>
#include <stim/parser/filename.h>
#include <stim/image/image.h>
#include <ANN/ANN.h>

//load a set of points from an ASCII file
std::vector< stim::vec<float> > points_from_list(std::string filename){

	std::vector< stim::vec<float> > result;		//create an array to store the result

	std::string str;

	//load the ground truth
	std::ifstream gfile(filename.c_str());
	while(std::getline(gfile, str)){
		stim::vec<float> v(str);
		result.push_back(v);
	}

	return result;			//return the list of points
}

//load a set of points from an image (nonzero values are point positions)
std::vector< stim::vec<float> > points_from_image(std::string filename){

	stim::image<unsigned char> I;		//create an image object

	I.load(filename);					//load the image file

	std::vector<unsigned long> idx = I.sparse_idx();		//get the 1D indices for each nonzero pixel

	std::vector< stim::vec<float> > result;					//generate an array to store the 2D pixel positions
	result.resize(idx.size());								//allocate the space to save time

	unsigned long w = I.width();							//get the image width (used to resolve the 2D pixel position from the 1D index)

	unsigned long x, y;
	for(unsigned long i = 0; i < idx.size(); i++){
		y = idx[i] / w;
		x = idx[i] - y * w;

		result[i] = stim::vec<float>(x, y);
	}


	return result;											//return the array of points

}

/// This function calculates the number of hits (
void calc_hitlist(std::vector< unsigned int >& hits,
					std::vector< float >& dists,
				    std::vector< stim::vec<float> > data,
				    std::vector< stim::vec<float> > query){

	int n_data = data.size();
	int n_query = query.size();

	//allocate space in the hit lists
	hits.resize(n_query);
	dists.resize(n_query);

	//calculate the spatial dimension
	unsigned int dim = 1;
	for(unsigned int i = 0; i < data.size(); i++)
		if(data[i].size() > dim) dim = data[i].size();

										//number of data points
	ANNpointArray dataPts = annAllocPts(n_data, dim);				//array of data points
	ANNpoint queryPt = annAllocPt(dim);							//query point
	ANNidxArray nnIdx = new ANNidx[1];							//indices of nearest neighbors
	ANNdistArray sq_dist = new ANNdist[1];						//array of squared distances to nearest neighbors

	//load the ground truth points into the ANN point array
	for(unsigned int i = 0; i < data.size(); i++){
		for(unsigned int d = 0; d < dim; d++){
			dataPts[i][d] = data[i][d];
		}
	}

	//create a KD-tree
	ANNkd_tree* kdTree;											//KD tree search structure
	kdTree = new ANNkd_tree(dataPts, n_data, dim);				//create the KD tree search structure


	//submit each query point to the KD tree to find the nearest neighbor
	for(unsigned int q = 0; q < n_query; q++){

		for(unsigned int d = 0; d < dim; d++){
			queryPt[d] = query[q][d];
		}

		//query the KD-tree
		kdTree->annkSearch(queryPt, 1, nnIdx, sq_dist);

		//put the result in the hit array
		hits[q] = nnIdx[0];
		dists[q] = sqrt(sq_dist[0]);
	}
}

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

	//output advertisement
	std::cout<<std::endl<<std::endl;
	std::cout<<"========================================================================="<<std::endl;
	std::cout<<"Thank you for using the PointMets validation tool!"<<std::endl;
	std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
	std::cout<<"Developers: David Mayerich"<<std::endl;
	std::cout<<"Source: https://github.com/stimlab/hsiproc.git"<<std::endl;
	std::cout<<"========================================================================="<<std::endl<<std::endl;


	stim::arglist args;

#ifdef _WIN32
	args.set_ansi(false);
#endif

	//add arguments
	args.add("help", "prints this help");
	args.add("r", "maximum possible distance to consider a point 'detected'", "10", "[positive value]");

	//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: pointmets ground_truth test_case --option [A B C ...]"<<std::endl;
		std::cout<<std::endl<<std::endl
				  << "examples: pointmets blue_gt.txt blue_t.txt"<<std::endl;
		std::cout << "          pointmets blue_gt.bmp blue_t.bmp" << 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 comparison, enter pointmets --help for options."<<std::endl<<std::endl;
		exit(1);
	}

	//get the ground truth and test file
	stim::filename gfilename(args.arg(0));
	stim::filename tfilename(args.arg(1));

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

	std::vector< stim::vec<float> > G;
	std::vector< stim::vec<float> > T;
	

	//if an ASCII list of points is specified
	if(gfilename.get_extension() == "bmp")
		G = points_from_image(gfilename.str());
	else
		G = points_from_list(gfilename.str());

	if(tfilename.get_extension() == "bmp")
		T = points_from_image(gfilename.str());
	else
		T = points_from_list(tfilename.str());

	//calculate the hit list and distances for the test array
	std::vector< unsigned int > T_hits;
	std::vector< float > T_dists;
	calc_hitlist(T_hits, T_dists, G, T);

	//calculate the hit list and distances for the ground truth array
	std::vector< unsigned int > G_hits;
	std::vector< float > G_dists;
	calc_hitlist(G_hits, G_dists, T, G);

	//count up metrics for each test point
	std::vector< unsigned int > TP;
	std::vector< unsigned int > FP;
	unsigned int g, t;
	for(t = 0; t < T.size(); t++){

		//a true positive has a matching hit in both hit lists, with a distance less than the radius
		g = T_hits[t];
		if(G_hits[g] == t && T_dists[t] <= r) TP.push_back(t);
		else FP.push_back(t);
	}

	std::vector< unsigned int > FN;
	for(g = 0; g < G.size(); g++){
		t = G_hits[g];
		if(T_hits[t] != g || G_dists[g] > r) FN.push_back(g);
	}

	std::cout<<"N: "<<G.size()<<std::endl;
	std::cout<<"TP: "<<TP.size()<<"   ( "<<(double)TP.size() / (double)G.size() * 100<<"% detected)"<<std::endl;
	std::cout<<"FN: "<<FN.size()<<"   ( "<<(double)FN.size() / (double)G.size() * 100<<"% undetected)"<<std::endl;
	std::cout<<"FP: "<<FP.size()<<"   ( "<<(double)FP.size() / (double)T.size() * 100<<"% incorrectly detected)"<<std::endl;



}