moved old filter cpp and updated readme

1 files changed, 666 insertions, 0 deletions

diff --git a/fly-tools/misc/FlyTrackingFilter.cpp b/fly-tools/misc/FlyTrackingFilter.cpp
new file mode 100644
index 0000000..9020522
--- /dev/null
+++ b/fly-tools/misc/FlyTrackingFilter.cpp
@@ -0,0 +1,666 @@
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <cmath>
+#include <vector>
+#include <list>
+#include <fstream>
+#include <cassert>
+#include <cstdlib>
+
+#include <ImageMagick/Magick++.h>
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_eigen.h>
+
+#include "FrameInfo.h"
+
+using namespace Magick;
+using namespace std;
+void findObj(Image* img, int x, int y, vector<pair<int,int> > & shape ,bool eightCon=true, bool colorLookingFor=true);
+void eightConnObj(Image* img, int x, int y, vector<pair<int, int> > & obj, bool color=true);
+void fourConnObj(Image* img, int x, int y, vector<pair<int, int> > & obj, bool color=true);
+vector<double> covariantDecomposition(vector<pair<int,int> > & points);
+pair<int,int> getCentroid(vector<pair<int,int> > & points);
+bool isInterface(Image* orig, unsigned int x, unsigned int y);
+void writeFrameImage(int fn, string imS);
+int roundT(double v) {return int(v+0.5);}
+
+
+const double PI = atan(1.0)*4.0;
+const double FACTOR_EIGEN = 100;
+
+Image* residual;
+
+ostream &operator<<(ostream &out, FlyObject & fO) {
+	fO.output(out);
+	return out;
+}
+
+ostream &operator<<(ostream &out, FrameInfo & fI) {
+	fI.output(out);
+	return out;
+}
+
+vector<vector<pair<int, int> > > shapeVectors;
+vector<pair<int,int> > shape;
+vector<pair<int, int> > sizeNIndexVector;
+
+void bubbleSort() {
+	
+	for(int i=1; i<sizeNIndexVector.size(); i++) {
+		for(int j=0; j<sizeNIndexVector.size()-i; j++) {
+			pair<int, int> a = sizeNIndexVector[j];
+			pair<int, int> b = sizeNIndexVector[j+1];
+			
+			if (a.first < b.first) {
+				pair<int, int> c = sizeNIndexVector[j];
+				sizeNIndexVector[j] = sizeNIndexVector[j+1];
+				sizeNIndexVector[j+1] = c;
+			}
+		}		
+	}
+
+}
+
+void fillResidualWithObj(vector<pair<int, int> > & obj, ColorRGB c)
+{
+	for (unsigned int i = 0; i<obj.size(); i++)
+		residual->pixelColor(obj[i].first, obj[i].second, c);
+}
+
+void writeHist(const char* filename, map<unsigned int, unsigned int> & len)
+{
+	map<unsigned int,unsigned int>::iterator front = len.begin(),
+	back = len.end();
+	back--;
+	
+	
+	unsigned int first = front->first, last = back->first;
+	/*if (cutoff != -1 && cutoff < int(last))
+	 last = cutoff;
+	 */
+	cout << "Min: " << first << endl
+	<< "Max: " << last << endl
+	<< "Count: " << last-first << endl;
+	//vector<unsigned int> hist(last-first, 0);
+	vector<unsigned int> hist(last+1, 0);
+	
+	cout << "hist size: " << hist.size() << endl;
+	try{
+		for(unsigned int j = 0; j<first; j++) {
+			hist[j] = 0;
+		}
+		for (unsigned int j = first; j<=last; j++)
+		{
+			/*if ( roundT(j-first) >= int(hist.size()) )
+			 hist.resize(j-first,0);
+			 hist[roundT(j-first)] = len[j];
+			 */
+			
+			/*if ( roundT(j) >= int(hist.size()) )
+			 hist.resize(j,0);
+			 hist[roundT(j)] = len[j];
+			 */
+			hist[j] = len[j];
+		}
+	}
+	catch (...)
+	{ cerr << "Bad histogram bucketing" << endl; }
+	
+	/*if ( (cutoff >= 0) && (cutoff<int(hist.size())) )
+	 hist.resize(cutoff);
+	 */
+	len.clear();
+	try
+	{
+		ofstream fout(filename);
+		for (unsigned int i = 0; i<hist.size(); i++) {
+			fout << hist[i] << endl;
+		}
+		fout << first << " " << last << " " << hist.size() << endl;
+		fout.close();
+	}
+	catch (...)
+	{ cerr << "Bad memory loc for opening file" << endl; }
+}
+
+int main(int argc, char* argv[])
+{
+	if (argc < 5)
+	{
+		cerr << "Usage: executablename <inputFile.txt> <ratio_largest_to_second_largest> <InputLocationOfMaskImage> <outputFolderName>" << endl; // input file contains name of the
+		// input image files
+		return -1;
+	}
+	
+	MagickCore::SetMagickResourceLimit(MagickCore::MemoryResource, 1536);
+	MagickCore::SetMagickResourceLimit(MagickCore::MapResource, 2048);
+	
+	string fileName;
+	ifstream inputFile(argv[1]);
+	if (inputFile.fail() ) {
+		cout << "cannot open the input file that contains name of the input images\n";
+		exit(1);
+	}
+
+	// get the output file name along with the location from argv[4]
+	string outputFileLocation(argv[4]);
+	string outputFileName = outputFileLocation + "final/outputFile.txt";
+	ofstream outputFile(outputFileName.c_str());
+	
+	// get the location of the input mask file
+	string inputMaskFileLocation(argv[3]);
+	
+	int frameCounter = 0;
+	// use 15, 20 or 10
+	// 15 found to be working correct
+	double ratioSecondLargestToLargest = atof(argv[2]);
+	
+	cout << "Ratio is 1/"<<ratioSecondLargestToLargest<< " = "<<(1/ratioSecondLargestToLargest)<<endl;
+	
+	outputFile<<"Ratio given 1/"<<ratioSecondLargestToLargest<<" = "<<(1/ratioSecondLargestToLargest)<<endl;
+	
+	ratioSecondLargestToLargest = 1/ratioSecondLargestToLargest;
+	// to find the largest object avg size
+	
+	long totalSize = 0;
+	
+	char buffer[100];
+	
+	while (inputFile>>fileName) {
+		
+		string savedFileName = fileName;
+	
+		// get the input mask file
+//		fileName = "input/"+fileName;
+		fileName =  inputMaskFileLocation + fileName;
+
+		Image* img = new Image(fileName.c_str());
+		int width = img->columns(),height = img->rows();
+		Image* imgForFilter;
+		sprintf(buffer,"%ix%i",width,height);
+		
+		// residual image is initialized with black representing not visited.
+		residual = new Image(buffer, "black");
+		
+		imgForFilter = new Image(buffer, "white");
+		cout << "Reading "<<savedFileName<<endl;
+		cout << "Filter wxh "<<width<<","<<height<<endl;
+		shape.clear();
+		// find the black background from location (0,0)
+		findObj(img, 0, 0, shape, false, false);
+		int s = shape.size();
+		if (s > 0)
+			cout << "black object size is "<<s;
+		for (int i=0; i<s; i++) {
+			imgForFilter->pixelColor(shape[i].first, shape[i].second, "black");
+		}
+		
+		// store the intermediate file in temp folder under the output location
+//		string oFilteredFileName = "output/filtered/temp/"+outputFile;
+		string oFilteredFileName =  outputFileLocation +"temp/"+savedFileName;
+		
+		cout << "Saving the filtered image "<< oFilteredFileName<<endl;
+		imgForFilter->write(oFilteredFileName.c_str());
+
+		delete residual;
+		/*
+		residual = new Image(buffer, "black");
+		
+		shapeVectors.clear();
+		sizeNIndexVector.clear();
+		
+		// find the two largest object
+		int objectCounter = 0;
+		for (int x=0; x<width; x++) {
+			for (int y=0; y<height; y++) {
+				//find the white object now using eight connected neighbours
+				shape.clear();
+				findObj(imgForFilter, x, y, shape, true, true);
+				int s = shape.size();
+				if (s>0) {
+					shapeVectors.push_back(shape);
+					cout << "new object pushed back at position "<<x<<","<<y<<" of size "<<s<<endl;
+					pair<int, int> si(s, objectCounter);
+					sizeNIndexVector.push_back(si);
+					objectCounter++;
+				}
+
+			}
+		}
+		
+		// sort the sizes and take the largest to find average largest
+		cout<<"shapeVectors size = "<<shapeVectors.size()<<endl;
+		
+		bubbleSort();
+		
+		cout <<"Largest object size is "<<sizeNIndexVector[0].first<<endl;
+		
+		// add the largest size to sum
+		totalSize += static_cast<long>(sizeNIndexVector[0].first);
+		
+		cout << "Current total size is "<<totalSize<<" for object "<<(frameCounter+1)<<endl;
+		
+		cout << "-----------------------------------------------------------"<<endl;
+		*/
+		
+		frameCounter++;
+		
+		delete imgForFilter;
+
+//		delete residual;
+		
+		delete img;
+	
+
+	}
+	
+	inputFile.close();
+	
+	/*
+	avgLargestSize = static_cast<long> (totalSize/frameCounter);
+	cout << "Average largest size is "<<avgLargestSize<<endl;
+	*/
+	
+	// previous loop calculates the average largest size
+	 
+	inputFile.open(argv[1]);
+	
+	if (inputFile.fail() == true) {
+		cout << "Cannot open the input file again"<<endl;
+		exit(1);
+	}
+	
+	
+	while (inputFile>>fileName) {
+		
+		//string inputFileName = "output/filtered/temp/"+fileName;
+		string inputFileName =  outputFileLocation + "temp/"+fileName;
+		
+		Image* img = new Image(inputFileName.c_str());
+		int width = img->columns();
+		int height = img->rows();
+		
+		outputFile<<"File name is "<<inputFileName<<endl;
+		outputFile<<"----------------------------------------------\n";
+		
+		sprintf(buffer,"%ix%i",width,height);
+		
+		// residual image is initialized with black representing not visited.
+		residual = new Image(buffer, "black");
+		
+		cout << "Reading "<<inputFileName<<endl;
+		cout << "Filter wxh "<<width<<","<<height<<endl;
+		
+		shapeVectors.clear();
+		sizeNIndexVector.clear();
+
+		// find the objects and sort according to size
+		int objectCounter = 0;
+		for (int x=0; x<width; x++) {
+			for (int y=0; y<height; y++) {
+				// find the white object using eight connected
+				shape.clear();
+				findObj(img, x, y, shape, true, true);
+				int s = shape.size();
+				
+				if (s > 0) {
+					
+					shapeVectors.push_back(shape);
+					cout << "New object found at position ("<<x<<","<<y<<") of size "<<s<<endl;
+					pair<int, int> si(s, objectCounter);
+					sizeNIndexVector.push_back(si);
+					objectCounter++;
+				
+				}
+					
+			}
+		}
+		
+		cout << "Total object found "<<sizeNIndexVector.size()<<endl;
+		
+		outputFile<<"Total objects found "<<sizeNIndexVector.size()<<endl;
+		
+		bubbleSort();
+		
+		// take the largest object
+		double currentLargestSize = static_cast<double>(sizeNIndexVector[0].first);
+		
+		double secondLargest = 0;
+		double ratio = 0;
+		if (sizeNIndexVector.size() > 1) {
+			
+			secondLargest = static_cast<double>(sizeNIndexVector[1].first);
+			ratio = secondLargest/currentLargestSize;
+			cout << "Ratio is "<<secondLargest<<"/"<<currentLargestSize<<" = "<<ratio<<endl;
+			outputFile<<"secondLargest = "<<secondLargest<<"\ncurrentLargest = "<<currentLargestSize<<"\nRatio = "<<ratio<<endl;
+			
+		}
+		
+		// find the largest to second largest ratio if it is less than the defined ratio then
+		// the objects are single object
+				
+		int numberOfObjects = 0;
+		
+		if (sizeNIndexVector.size() == 1) {
+			cout << "Frame contains one object\n";
+			outputFile << "Frame contains one object\n";
+			numberOfObjects = 1;
+		}
+		else if (ratio <= ratioSecondLargestToLargest ) {
+			cout << "Single object calculated in the frame because current ratio = "<<ratio<<" is less than defined ratio = "<<ratioSecondLargestToLargest<<endl;
+			outputFile<< "Single object calculated in the frame because current ratio = "<<ratio<<" is less than defined ratio = "<<ratioSecondLargestToLargest<<endl;
+			numberOfObjects = 1;
+		} else { 
+			cout << "Two objects in the frame\n";
+			outputFile<<"Two objects in the frame\n";
+			numberOfObjects = 2;
+		}
+		
+		cout << "Total object found "<<numberOfObjects<<endl;
+		outputFile << "Total object found "<<numberOfObjects<<endl;
+		
+		Image* imgFinal = new Image(buffer, "black");
+				
+		for (int n=0; n<numberOfObjects; n++) {
+			
+			int totalPoints = sizeNIndexVector[n].first;
+			
+			for (int i=0; i<totalPoints; i++) {
+				
+				imgFinal->pixelColor(shapeVectors[ sizeNIndexVector[n].second ][i].first, shapeVectors[ sizeNIndexVector[n].second ][i].second, "white");
+			}
+			
+
+		}
+		
+		//string finalImageName = "output/filtered/final/"+fileName;
+		
+		string finalImageName = outputFileLocation + "final/"+fileName;
+		
+		imgFinal->write( finalImageName.c_str() );
+		
+		
+		
+		
+		// writing the single in red
+		if (numberOfObjects == 1) {
+			
+			Image* singleObjectFinal = new Image(buffer, "black");
+		
+			int totalPoints = sizeNIndexVector[0].first;
+			
+			cout << "Output the single object of size = "<<totalPoints<<endl;
+			
+			for (int i=0; i<totalPoints; i++) {
+				
+				singleObjectFinal->pixelColor(shapeVectors[ sizeNIndexVector[0].second ][i].first, shapeVectors[ sizeNIndexVector[0].second ][i].second, "red");
+			}
+			
+			//string singleImageName = "output/filtered/single/"+fileName;
+			string singleImageName = outputFileLocation + "single/"+fileName;
+			
+			singleObjectFinal->write(singleImageName.c_str());
+			
+			delete singleObjectFinal;
+			
+		}
+		
+		outputFile<<"----------------------------------------------------\n";
+		
+			
+		delete img;
+		
+		delete residual;
+		
+		delete imgFinal;
+		
+		
+	}
+	
+	inputFile.close();
+	outputFile.close();
+	
+	return 0;
+}
+
+
+void findObj(Image* img, int x, int y, vector<pair<int,int> > & shape ,bool eightCon, bool colorLookingFor)
+{
+	assert(residual != NULL);
+
+	if (eightCon == true)
+		eightConnObj(img, x, y, shape, colorLookingFor);
+	else {
+		fourConnObj(img, x, y, shape, colorLookingFor);
+	}
+}
+
+int barrier = 1000;
+void fourConnObj(Image* img, int x, int y, vector<pair<int, int> > & obj, bool color)
+{
+	int width = img->columns(),height = img->rows();
+	
+	// boundary violation check
+	if ( (x >= (width)) || (x < 0) || (y >= (height) ) || (y < 0) )
+		return;
+	
+	// residualpixel.mono() == true implies it is visited. Otherwise not visited.
+	ColorMono residualpixel = ColorMono(residual->pixelColor(x,y));
+	// originalpixel.mono() == true implies it is an object pixel. Otherwise it is blank region pixel.
+	ColorMono originalpixel = ColorMono(img->pixelColor(x,y));
+	
+	// If the current pixel is already visited then return
+	if (residualpixel.mono() == true)
+		return;
+	
+	// Else if current pixel is not visited and it is black, which means it is not an object pixel; so return
+	else if (residualpixel.mono() == false && originalpixel.mono() != color)
+		return;
+	// If current pixel is not visited and its value is white, which means a new object is found.
+	else if (residualpixel.mono() == false && originalpixel.mono() == color) {
+		// Save the coordinates of the current pixel into the vector and make the pixel visited in the residual image
+		pair<int,int> p;
+		p.first = x;
+		p.second = y;
+		obj.push_back(p);
+		
+//		if (obj.size() > barrier) {
+//			cout<<obj.size()<<endl;
+//			barrier = barrier + 1000;
+//		}
+		// setting the residual image at pixel(x,y) to white.
+		residual->pixelColor(x,y, ColorMono(true));
+		
+		// Recursively call all of it's eight neighbours.
+		fourConnObj(img, x+1, y, obj, color);
+		fourConnObj(img, x, y-1, obj, color);
+		
+		fourConnObj(img, x-1, y, obj, color);
+		fourConnObj(img, x, y+1, obj, color);
+	}
+	
+}
+
+void eightConnObj(Image* img, int x, int y, vector<pair<int, int> > & obj, bool color)
+{
+	int width = img->columns(),height = img->rows();
+	
+	// boundary violation check
+	if ( (x >= (width)) || (x < 0) || (y >= (height) ) || (y < 0) )
+		return;
+	
+	// residualpixel.mono() == true implies it is visited. Otherwise not visited.
+	ColorMono residualpixel = ColorMono(residual->pixelColor(x,y));
+	// originalpixel.mono() == true implies it is an object pixel. Otherwise it is blank region pixel.
+	ColorMono originalpixel = ColorMono(img->pixelColor(x,y));
+	
+	// If the current pixel is already visited then return
+	if (residualpixel.mono() == true)
+		return;
+	
+	// Else if current pixel is not visited and it is black, which means it is not an object pixel; so return
+	else if (residualpixel.mono() == false && originalpixel.mono() != color)
+		return;
+	// If current pixel is not visited and its value is white, which means a new object is found.
+	else if (residualpixel.mono() == false && originalpixel.mono() == color) {
+		// Save the coordinates of the current pixel into the vector and make the pixel visited in the residual image
+		pair<int,int> p;
+		p.first = x;
+		p.second = y;
+		obj.push_back(p);
+
+//		if (obj.size() > barrier) {
+//			//cout<<obj.size()<<endl;
+//			barrier = barrier + 1000;
+//		}
+		// setting the residual image at pixel(x,y) to white.
+		residual->pixelColor(x,y, ColorMono(true));
+		
+		// Recursively call all of it's eight neighbours.
+		eightConnObj(img, x+1, y, obj, color);
+		eightConnObj(img, x+1, y-1, obj, color);
+		eightConnObj(img, x, y-1, obj, color);
+		eightConnObj(img, x-1, y-1, obj, color);
+		
+		eightConnObj(img, x-1, y, obj, color);
+		eightConnObj(img, x-1, y+1, obj, color);
+		eightConnObj(img, x, y+1, obj, color);
+		eightConnObj(img, x+1, y+1, obj, color);		
+
+	}
+	
+}
+
+// Aspect Ratio
+pair<int,int> getCentroid(vector<pair<int,int> > & points)
+{
+	pair<int,int> centroid;
+	centroid.first = 0;
+	centroid.second = 0;
+	
+	for (unsigned int i = 0; i<points.size(); i++)
+	{
+		centroid.first += points[i].first;
+		centroid.second += points[i].second;
+	}
+	
+	centroid.first = roundT(double(centroid.first)/points.size());
+	centroid.second = roundT(double(centroid.second)/points.size());
+	
+	return centroid;
+}
+
+
+vector<double> covariantDecomposition(vector<pair<int,int> > & points)
+{
+	unsigned int i,j,k;
+	pair<int,int> centroid = getCentroid(points);
+	vector<double> retval;
+	
+	gsl_matrix* matrice = gsl_matrix_alloc(2, 2);
+	
+	double sumX2 = 0, sumXY = 0, sumY2 = 0;
+	for (k = 0; k<points.size(); k++)
+	{
+		sumX2 += pow(double(points[k].first - centroid.first),2.0);
+		sumY2 += pow(double(points[k].second - centroid.second),2.0);
+		// should we take the absolute value of X*Y
+		sumXY += (points[k].first - centroid.first) * (points[k].second - centroid.second);
+	}
+	gsl_matrix_set(matrice, 0, 0, roundT(sumX2/points.size()));
+	gsl_matrix_set(matrice, 0, 1, roundT(sumXY/points.size()));
+	gsl_matrix_set(matrice, 1, 0, roundT(sumXY/points.size()));
+	gsl_matrix_set(matrice, 1, 1, roundT(sumY2/points.size()));
+	
+	//	outputMatrix("Covariant", matrice);
+	
+	// This function allocates a workspace for computing eigenvalues of n-by-n
+	// real symmetric matrices. The size of the workspace is O(2n).
+	gsl_eigen_symmv_workspace* eigenSpace = gsl_eigen_symmv_alloc(2);
+	gsl_vector* eigenVal = gsl_vector_alloc(2);
+	gsl_matrix* eigenVec = gsl_matrix_alloc(2, 2);
+	// This function computes the eigenvalues and eigenvectors of the real
+	// symmetric matrix A. Additional workspace of the appropriate size must be
+	// provided in w. The diagonal and lower triangular part of A are destroyed
+	// during the computation, but the strict upper triangular part is not
+	// referenced. The eigenvalues are stored in the vector eval and are unordered.
+	// The corresponding eigenvectors are stored in the columns of the matrix evec.
+	// For example, the eigenvector in the first column corresponds to the first
+	// eigenvalue. The eigenvectors are guaranteed to be mutually orthogonal and
+	// normalised to unit magnitude.
+	gsl_eigen_symmv (matrice, eigenVal, eigenVec, eigenSpace);
+	gsl_eigen_symmv_free (eigenSpace);
+	
+	gsl_eigen_symmv_sort(eigenVal, eigenVec, GSL_EIGEN_SORT_VAL_ASC);
+	
+	for (i = 0; i<eigenVal->size; i++)
+		retval.push_back(gsl_vector_get(eigenVal, i));
+	
+	for (j = 0; j<eigenVec->size2; j++)
+		for (i = 0; i<eigenVec->size1; i++)
+			retval.push_back(gsl_matrix_get(eigenVec, i, j));
+	
+	retval.push_back(static_cast<double>(centroid.first));
+	retval.push_back(static_cast<double> (centroid.second));
+	
+//	for (i=0; i<2; i++) {
+//		gsl_vector_view evec_i = gsl_matrix_column (eigenVec, i);
+//		//printf ("eigenvalue = %g\n", eval_i);
+//		cout<<"eigenvector = \n";
+//		gsl_vector_fprintf (stdout, &evec_i.vector, "%g");
+//	}
+	
+	gsl_vector_free(eigenVal);
+	gsl_matrix_free(matrice);
+	gsl_matrix_free(eigenVec);
+	
+	return retval;
+}
+
+// isInterface for binary image
+bool isInterface(Image* orig, unsigned int x, unsigned int y)
+{
+	// Get the current pixel's color
+	ColorMono currentpixel = (ColorMono)orig->pixelColor(x,y);
+	// If the current pixel is black pixel then it is not boundary pixel
+	// error check
+	if (currentpixel.mono() == false)
+		return false;
+	
+	// If the current pixel is not black then it is white. So, now we need
+	// to check whether any four of its neighbor pixels (left, top, right,
+	// bottom ) is black. If any of this neighbor is black then current
+	// pixel is a neighbor pixel. Otherwise current pixel is not neighbor
+	// pixel.
+	
+	ColorMono leftneighborpixel = (ColorMono)orig->pixelColor(x-1,y);
+	ColorMono topneighborpixel = (ColorMono)orig->pixelColor(x,y-1);
+	ColorMono rightneighborpixel = (ColorMono)orig->pixelColor(x+1,y);
+	ColorMono bottomneighborpixel = (ColorMono)orig->pixelColor(x,y+1);
+	
+	// If leftneighborpixel is black and currentpixel is white then it is
+	// boundary pixel
+	if ( leftneighborpixel.mono() != currentpixel.mono())
+		return true;
+	// If topneighborpixel is black and currentpixel is white then it is
+	// boundary pixel
+	else if (topneighborpixel.mono() != currentpixel.mono())
+		return true;
+	// If rightneighborpixel is black and currentpixel is white then it
+	// is boundary pixel
+	else if (rightneighborpixel.mono() != currentpixel.mono())
+		return true;
+	// If bottomneighborpixel is black and currentpixel is white then it
+	// is boundary pixel
+	else if (bottomneighborpixel.mono() != currentpixel.mono())
+		return true;
+	// Else all of its neighbor pixels are white so it can not be a
+	// boundary pixel
+	else
+		return false;
+	
+}
\ No newline at end of file