From 53b73dec4d6e71899d183335509591b7c89a419b Mon Sep 17 00:00:00 2001 From: mutantturkey Date: Thu, 14 Jun 2012 18:34:23 -0400 Subject: pushing the fly-tools into this repository --- fly-tools/FlyObject.cpp | 118 ++ fly-tools/FlyObject.h | 49 + fly-tools/FlyTrackingFilter.cpp | 666 ++++++ fly-tools/FlyTrackingIterativeFilter.cpp | 707 +++++++ fly-tools/FlyTrackingMain.cpp | 3181 +++++++++++++++++++++++++++++ fly-tools/FlyTrackingNoImages | Bin 0 -> 97256 bytes fly-tools/FrameInfo.cpp | 72 + fly-tools/FrameInfo.h | 35 + fly-tools/FurthestPointAlongMajorAxis.cpp | 3084 ++++++++++++++++++++++++++++ fly-tools/MyPair.h | 30 + fly-tools/StandardDeviation.cpp | 90 + fly-tools/misc/linesweep.cpp | 693 +++++++ fly-tools/misc/makeimage.cpp | 51 + fly-tools/misc/test.cpp | 92 + 14 files changed, 8868 insertions(+) create mode 100644 fly-tools/FlyObject.cpp create mode 100644 fly-tools/FlyObject.h create mode 100644 fly-tools/FlyTrackingFilter.cpp create mode 100644 fly-tools/FlyTrackingIterativeFilter.cpp create mode 100644 fly-tools/FlyTrackingMain.cpp create mode 100755 fly-tools/FlyTrackingNoImages create mode 100644 fly-tools/FrameInfo.cpp create mode 100644 fly-tools/FrameInfo.h create mode 100644 fly-tools/FurthestPointAlongMajorAxis.cpp create mode 100644 fly-tools/MyPair.h create mode 100644 fly-tools/StandardDeviation.cpp create mode 100644 fly-tools/misc/linesweep.cpp create mode 100644 fly-tools/misc/makeimage.cpp create mode 100644 fly-tools/misc/test.cpp (limited to 'fly-tools') diff --git a/fly-tools/FlyObject.cpp b/fly-tools/FlyObject.cpp new file mode 100644 index 0000000..d6cbee3 --- /dev/null +++ b/fly-tools/FlyObject.cpp @@ -0,0 +1,118 @@ +/* + * FlyObject.cpp + * + * + * Created by Md. Alimoor Reza on 6/26/10. + * Copyright 2010 Drexel University. All rights reserved. + * + */ + +#include "FlyObject.h" + +FlyObject::FlyObject(int area, pair centroid, pair majorAxisEV, pair velocityV,bool headIsInDirectionMAEV, pair head, double speed) { +//FlyObject::FlyObject(int area, pair centroid, pair majorAxisEV, pair velocityV, vector > areaCoord) { + this->area = area; + this->centroid = centroid; + this->majorAxisEV = majorAxisEV; + this->velocityV = velocityV; + this->headIsInDirectionMAEV = headIsInDirectionMAEV; + this->head = head; + this->speed = speed; +// this->areaCoord = areaCoord; +} + +FlyObject::FlyObject(const FlyObject &f){ + this->area = f.getArea(); + this->centroid=f.getCentroid(); + this->majorAxisEV =f.getMajorAxisEV(); + this->velocityV = f.getVelocityV(); +// this->areaCoord = f.getAreaCoord(); + this->headIsInDirectionMAEV = f.getHeadIsInDirectionMAEV(); + this->head = f.getHead(); + this->speed = f.getSpeed(); +} + +int FlyObject::getArea() const { + return area; +} + +pair FlyObject::getCentroid() const { + return this->centroid; +} + +pair FlyObject::getMajorAxisEV() const { + return this->majorAxisEV; +} + +pair FlyObject::getVelocityV() const { + return this->velocityV; +} +bool FlyObject::getHeadIsInDirectionMAEV() const { + return this->headIsInDirectionMAEV; +} +pair FlyObject::getHead() const { + return this->head; +} +//vector > FlyObject::getAreaCoord() const { +// return this->areaCoord; +//} +double FlyObject::getSpeed() const { + return this->speed; +} +void FlyObject::setArea(int area) { + this->area = area; +} + +void FlyObject::setCentroid(pair centroid) { + this->centroid = centroid; +} + +void FlyObject::setMajorAxisEV(pair majorAxisEV) { + this->majorAxisEV = majorAxisEV; +} + +void FlyObject::setVelocityV(pair velocityV) { + this->velocityV = velocityV; +// cout << "velocityV is set to "<velocityV.first<<","<velocityV.second< head) { + this->head = head; + cout << "new head is set"< > areaCoord){ +// this->areaCoord = areaCoord; +//} +void FlyObject::setSpeed(double speed) { + this->speed = speed; +} + +void FlyObject::normalizeVelocity() { + double temp = velocityV.first*velocityV.first + velocityV.second*velocityV.second; + temp = sqrt(temp); + cout << "sum = "<headIsInDirectionMAEV = headIsInDirectionMAEV; + cout << "head flag set to "<headIsInDirectionMAEV<headIsInDirectionMAEV<head.first<<","< +#include +#include +#include +using namespace std; + +class FlyObject { +public: + //FlyObject(int area, pair centroid, pair majorAxisEV, pair velocityV, vector > areaCoord); + FlyObject(int area, pair centroid, pair majorAxisEV, pair velocityV,bool headIsInDirectionMAEV, pair head, double speed); + FlyObject(const FlyObject &f); + int getArea() const; + pair getCentroid() const; + pair getMajorAxisEV() const; + pair getVelocityV() const; + bool getHeadIsInDirectionMAEV() const; + pair getHead() const; + double getSpeed() const; + //vector > getAreaCoord() const; + void setArea(int area); + void setCentroid(pair); + void setMajorAxisEV(pair); + void setVelocityV(pair); + void normalizeVelocity(); + void setHeadIsInDirectionMAEV(bool); + void setHead(pair head); + void setSpeed(double speed); + //void setAreaCoord(vector >); + void output(ostream &out); + + +private: + int area; + //vector > areaCoord; + pair centroid; + pair majorAxisEV; + pair velocityV; + bool headIsInDirectionMAEV; + pair head; + double speed; +}; diff --git a/fly-tools/FlyTrackingFilter.cpp b/fly-tools/FlyTrackingFilter.cpp new file mode 100644 index 0000000..9020522 --- /dev/null +++ b/fly-tools/FlyTrackingFilter.cpp @@ -0,0 +1,666 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "FrameInfo.h" + +using namespace Magick; +using namespace std; +void findObj(Image* img, int x, int y, vector > & shape ,bool eightCon=true, bool colorLookingFor=true); +void eightConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +void fourConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +vector covariantDecomposition(vector > & points); +pair getCentroid(vector > & 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 > > shapeVectors; +vector > shape; +vector > sizeNIndexVector; + +void bubbleSort() { + + for(int i=1; i a = sizeNIndexVector[j]; + pair b = sizeNIndexVector[j+1]; + + if (a.first < b.first) { + pair c = sizeNIndexVector[j]; + sizeNIndexVector[j] = sizeNIndexVector[j+1]; + sizeNIndexVector[j+1] = c; + } + } + } + +} + +void fillResidualWithObj(vector > & obj, ColorRGB c) +{ + for (unsigned int i = 0; ipixelColor(obj[i].first, obj[i].second, c); +} + +void writeHist(const char* filename, map & len) +{ + map::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 hist(last-first, 0); + vector hist(last+1, 0); + + cout << "hist size: " << hist.size() << endl; + try{ + for(unsigned int j = 0; j= 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 " << 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/"<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 "< 0) + cout << "black object size is "<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<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; x0) { + shapeVectors.push_back(shape); + cout << "new object pushed back at position "< si(s, objectCounter); + sizeNIndexVector.push_back(si); + objectCounter++; + } + + } + } + + // sort the sizes and take the largest to find average largest + cout<<"shapeVectors size = "<(sizeNIndexVector[0].first); + + cout << "Current total size is "< (totalSize/frameCounter); + cout << "Average largest size is "<>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 "< 0) { + + shapeVectors.push_back(shape); + cout << "New object found at position ("< si(s, objectCounter); + sizeNIndexVector.push_back(si); + objectCounter++; + + } + + } + } + + cout << "Total object found "<(sizeNIndexVector[0].first); + + double secondLargest = 0; + double ratio = 0; + if (sizeNIndexVector.size() > 1) { + + secondLargest = static_cast(sizeNIndexVector[1].first); + ratio = secondLargest/currentLargestSize; + cout << "Ratio is "<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 = "<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 > & 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 > & 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 p; + p.first = x; + p.second = y; + obj.push_back(p); + +// if (obj.size() > barrier) { +// cout<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 > & 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 p; + p.first = x; + p.second = y; + obj.push_back(p); + +// if (obj.size() > barrier) { +// //cout<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 getCentroid(vector > & points) +{ + pair centroid; + centroid.first = 0; + centroid.second = 0; + + for (unsigned int i = 0; i covariantDecomposition(vector > & points) +{ + unsigned int i,j,k; + pair centroid = getCentroid(points); + vector retval; + + gsl_matrix* matrice = gsl_matrix_alloc(2, 2); + + double sumX2 = 0, sumXY = 0, sumY2 = 0; + for (k = 0; ksize; i++) + retval.push_back(gsl_vector_get(eigenVal, i)); + + for (j = 0; jsize2; j++) + for (i = 0; isize1; i++) + retval.push_back(gsl_matrix_get(eigenVec, i, j)); + + retval.push_back(static_cast(centroid.first)); + retval.push_back(static_cast (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 diff --git a/fly-tools/FlyTrackingIterativeFilter.cpp b/fly-tools/FlyTrackingIterativeFilter.cpp new file mode 100644 index 0000000..8ecb55b --- /dev/null +++ b/fly-tools/FlyTrackingIterativeFilter.cpp @@ -0,0 +1,707 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "FrameInfo.h" +#include "MyPair.h" + +using namespace Magick; +using namespace std; +void findObj(Image* img, int x, int y, vector > & shape ,bool eightCon=true, bool colorLookingFor=true); +void eightConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +void fourConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +void findObjIterative(Image* img, int x, int y, vector > & shape, bool eightCon, double colorLookingFor); +void fourConnObjIterative(Image* img, int x, int y, vector > & obj, double colorLookingFor); +vector covariantDecomposition(vector > & points); +pair getCentroid(vector > & 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; +Image* imgForFilter; + +ostream &operator<<(ostream &out, FlyObject & fO) { + fO.output(out); + return out; +} + +ostream &operator<<(ostream &out, FrameInfo & fI) { + fI.output(out); + return out; +} + +vector > > shapeVectors; +vector > shape; +vector > sizeNIndexVector; + +void bubbleSort() { + + for(int i=1; i a = sizeNIndexVector[j]; + pair b = sizeNIndexVector[j+1]; + + if (a.first < b.first) { + pair c = sizeNIndexVector[j]; + sizeNIndexVector[j] = sizeNIndexVector[j+1]; + sizeNIndexVector[j+1] = c; + } + } + } + +} + +void fillResidualWithObj(vector > & obj, ColorRGB c) +{ + for (unsigned int i = 0; ipixelColor(obj[i].first, obj[i].second, c); +} + +void writeHist(const char* filename, map & len) +{ + map::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 hist(last-first, 0); + vector hist(last+1, 0); + + cout << "hist size: " << hist.size() << endl; + try{ + for(unsigned int j = 0; j= 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 " << endl; // input file contains name of the + // input image files + return -1; + } + + //MagickCore::SetMagickResourceLimit(MagickCore::MemoryResource, 1536); + //MagickCore::SetMagickResourceLimit(MagickCore::MapResource, 2048); + + + // Arg 1 is a file name, this is the actual name (First10MinSet48_0000001.png) + string fileName = argv[1]; + // Arg 2 is the ratio of the second largest to largest. Use 15. + double ratioSecondLargestToLargest = atof(argv[2]); + // Arg 3 is the location of the Masks + string inputMaskFileLocation(argv[3]); + // Arg 4 is the output folder for the Filtered images. + string outputFileLocation(argv[4]); + + ratioSecondLargestToLargest = 1/ratioSecondLargestToLargest; + + char buffer[100]; + + string savedFileName = fileName; + string finalImageName = outputFileLocation + "final/"+ savedFileName; + + // get the input mask file + fileName = inputMaskFileLocation + fileName; + + Image* original = new Image(fileName.c_str()); + int width = original->columns(); + int height = original->rows(); + sprintf(buffer,"%ix%i",width,height); + + imgForFilter = new Image(buffer, "white"); + shape.clear(); + + // find the black background from location (0,0) + ColorMono topLeftColor = ColorMono(original->pixelColor(0,0)); + ColorMono topRightColor = ColorMono(original->pixelColor(width-1,0)); + ColorMono bottomRightColor = ColorMono(original->pixelColor(width-1,height-1)); + ColorMono bottomLeftColor = ColorMono(original->pixelColor(0,height-1)); + + if (topLeftColor.mono() == false) { + + findObjIterative(original, 0, 0, shape, false, 0.0); + + } else if (topRightColor.mono() == false) { + + cout << "Top left is not black pixel for FLOOD FILLING so flood filling from top right\n"; + findObjIterative(original, width-1, 0, shape, false, 0.0); + + } else if (bottomRightColor.mono() == false) { + + cout << "Top left/Top right are not black pixel for FLOOD FILLING so flood filling from bottom right\n"; + findObjIterative(original, width-1, height-1, shape, false, 0.0); + + } else { + + cout << "Top left/top right/bottom right are not black pixel for FLOOD FILLING so flood filling from the bottomleft\n"; + findObjIterative(original, 0, height-1, shape, false, 0.0); + + } + + string inputFileName = outputFileLocation + "temp/" + savedFileName; + + Image* final_image = imgForFilter; + + + sprintf(buffer,"%ix%i",width,height); + + // residual image is initialized with black representing not visited. + residual = new Image(buffer, "black"); + + + shapeVectors.clear(); + sizeNIndexVector.clear(); + + // find the objects and sort according to size + int objectCounter = 0; + for (int x=0; x 0) { + shapeVectors.push_back(shape); + pair si(s, objectCounter); + sizeNIndexVector.push_back(si); + objectCounter++; + } + } + } + + + bubbleSort(); + + // take the largest object + double currentLargestSize = static_cast(sizeNIndexVector[0].first); + double secondLargest = 0; + double ratio = 0; + + if (sizeNIndexVector.size() > 1) { + + secondLargest = static_cast(sizeNIndexVector[1].first); + ratio = secondLargest/currentLargestSize; + + } + + // 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) { + numberOfObjects = 1; + } + else if (ratio <= ratioSecondLargestToLargest ) { + numberOfObjects = 1; + } else { + numberOfObjects = 2; + } + + Image* imgFinal = new Image(buffer, "black"); + + for (int n=0; npixelColor(shapeVectors[ sizeNIndexVector[n].second ][i].first, shapeVectors[ sizeNIndexVector[n].second ][i].second, "white"); + } + + } + + // write final image + cout << finalImageName << " \r"; + 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 = "<pixelColor(shapeVectors[ sizeNIndexVector[0].second ][i].first, shapeVectors[ sizeNIndexVector[0].second ][i].second, "red"); +// } +// +// //string singleImageName = "output/filtered/single/"+fileName; +// string singleImageName = outputFileLocation + "single/"+ savedFileName; +// +// singleObjectFinal->write(singleImageName.c_str()); +// +// } +// +} + + return 0; +} + +void findObjIterative(Image* img, int x, int y, vector > & shape, bool eightCon, double colorLookingFor) { + + assert(imgForFilter != NULL); + + //if (eightCon == true) + // eightConnObjIterative(img, x, y, shape, colorLookingFor); + //else { + fourConnObjIterative(img, x, y, shape, colorLookingFor); + + //} + + +} + +void fourConnObjIterative(Image* img, int x, int y, vector > & obj, double colorLookingFor) { + + /* + Flood-fill (node, target-color, replacement-color): + 1. Set Q to the empty queue. + 2. If the color of node is not equal to target-color, return. + 3. Add node to Q. + 4. For each element n of Q: + 5. If the color of n is equal to target-color: + 6. Set w and e equal to n. + 7. Move w to the west until the color of the node to the west of w no longer matches target-color. + 8. Move e to the east until the color of the node to the east of e no longer matches target-color. + 9. Set the color of nodes between w and e to replacement-color. + 10. For each node n between w and e: + 11. If the color of the node to the north of n is target-color, add that node to Q. + If the color of the node to the south of n is target-color, add that node to Q. + 12. Continue looping until Q is exhausted. + 13. Return. + + */ + + + queue< MyPair > Q; + + ColorRGB imgpixel = ColorRGB(img->pixelColor(x,y)); + + if ( (imgpixel.red() != colorLookingFor) and (imgpixel.green() != colorLookingFor) and (imgpixel.blue() != colorLookingFor)) { + + cout << "Returning without floodfilling because the first pixel is not the colorLookingFor"<columns(),height = img->rows(); + + while (Q.empty() != true) { + + MyPair n = Q.front(); + Q.pop(); + + ColorRGB westColor; + ColorRGB eastColor; + ColorRGB nColor; + MyPair i,j; + + nColor = ColorRGB(img->pixelColor(n.first, n.second)); + + if ( (nColor.red() == colorLookingFor) and (nColor.green() == colorLookingFor) and (nColor.blue() == colorLookingFor)) { + + //cout << "Current pixel is of the black color ("<pixelColor(e.first, e.second)); + else { + //cout << "outside of the image boundary in x direction so break the while loop for east"<pixelColor(i.first, i.second, ColorRGB(colorLookingFor, colorLookingFor, colorLookingFor) ); + + // change the color to green to indicate that it is visited + img->pixelColor(i.first, i.second, ColorRGB(0.0, 1.0, 0.0)); + + //cout << "Current pixel visited "<=0 ) { + + MyPair n(i.first, i.second-1); + ColorRGB northColor = ColorRGB(img->pixelColor(n.first, n.second)); + if ((northColor.red() == colorLookingFor) and (northColor.green() == colorLookingFor) and (northColor.blue() == colorLookingFor)) { + Q.push(n); + //cout << "North pixel not visited so pushed "<pixelColor(s.first, s.second)); + //cout << "South color "<size; i++) + retval.push_back(gsl_vector_get(eigenVal, i)); + + for (j = 0; jsize2; j++) + for (i = 0; isize1; i++) + retval.push_back(gsl_matrix_get(eigenVec, i, j)); + + retval.push_back(static_cast(centroid.first)); + retval.push_back(static_cast (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) +{ + 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; + +} diff --git a/fly-tools/FlyTrackingMain.cpp b/fly-tools/FlyTrackingMain.cpp new file mode 100644 index 0000000..4114a08 --- /dev/null +++ b/fly-tools/FlyTrackingMain.cpp @@ -0,0 +1,3181 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "FrameInfo.h" + +using namespace Magick; +using namespace std; + +const double PI = atan(1.0)*4.0; +const double FACTOR_EIGEN = 100; +const int STUCKING_TO_A_SINGLE_BLOB = 1; +const int SEPARATING_FROM_SINGLE_BLOB = 2; + + +Image* residual; + +vector fOVector; +vector fIVector; +// temporary store of the frames for finding initial head direction +vector fIForHeadVector; +// since for the first time the head is automatically need to be set +int startIndexToFindNewHead=0; +int endIndexToFindNewHead; +double initLargeLocationX=-1; +double initLargeLocationY=-1; +double initSmallLocationX=-1; +double initSmallLocationY=-1; +pair largeCollisionBeforeDirection; +pair smallCollisionBeforeDirection; + +// start location in the vector where two object get together +int startIndexSingleBlob=-1; +int endIndexSingleBlob =-1; + +// indices of the new set of sequences +int startOfATrackSequence = 0; // fix it later on inside main +int endOfATrackSequence=-1; +int sequenceSize=1; + +// indices for printing the one object frames +int startOfAOneObject = -1; +int endOfAOneObject = -1; + +vector fnVector; +string inputFileName; + +// GLOBAL PATHS +string maskImagePath; +string origImagePath; +string finalOutputPath; +string outputFilePrefix; + +vector > velocityDirectionsF; +vector > velocityDirectionsS; + +pair avgVelocityF; +pair avgVelocityS; + +pair overAllVelocityF; +pair overAllVelocityS; + +vector > evDirectionF; +vector > evDirectionS; + +int totalMaleLookingAtFemale = 0; +int totalFemaleLookingAtMale = 0; +int totalSingleBlob = 0; +int totalUnprocessedFrame = 0; +int totalSeparated = 0; +map centroidDistanceMap; +map headDirAngleMap; +map speedMap; + +void initSequence(){ + + startOfATrackSequence = -1; + endOfATrackSequence = -1; + sequenceSize = 1; + + +} + +ostream &operator<<(ostream &out, FlyObject & fO) { + fO.output(out); + return out; +} + +ostream &operator<<(ostream &out, FrameInfo & fI) { + fI.output(out); + return out; +} + +void bubbleSort(vector & fov) { + + //FlyObject a,b,c; + for(int i=1; i dataMap) +{ + cout << "In the beginning of the write hist\n"; + cout << "dataMap size "<::iterator front = dataMap.begin(), + back = dataMap.end(); + back--; + + + unsigned int first = front->first, last = back->first; + cout << "Min: " << first << endl<< "Max: " << last << endl<< "Count: " << last-first << endl; + //vector hist(last-first, 0); + vector hist(last+1, 0); + + cout << "hist size: " << hist.size() << endl; + try{ + for(unsigned int j = 0; j= 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] = dataMap[j]; + } + } + catch (...) + { cerr << "Bad histogram bucketing" << endl; } + + try + { + ofstream fout(filename); + for (unsigned int i = 0; i > & shape ,bool eightCon=true, bool colorLookingFor=true); +void eightConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +void fourConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +vector covariantDecomposition(vector > & points); +pair getCentroid(vector > & points); +bool isInterface(Image* orig, unsigned int x, unsigned int y); +void writeFrameImage(int fn, string imS); +void drawTheFlyObject(FrameInfo currentFI, string fileName, int isFirst, bool singleBlob=false,bool unprocessed = false); +void drawTheSequence(int startIndex, int endIndex, int isFirst, bool singleBlob = false, bool unprocessed = false); +double euclideanDist(FlyObject a, FlyObject b); +bool identifyFlyObjectFromFrameToFrame(FrameInfo prevFI, FrameInfo& currentFI, bool gotRidOfSingleBlob=false) ; +int roundT(double v) {return int(v+0.5);} +void determineHeadDirection(int fileCounter); + + +void normalizeVector(pair &a); +double calculateDotProduct(pair v, pair eV); +void calculateHeadVector(FlyObject fO, pair &headDirection); + +void normalizeVector(pair &a) { + double temp = a.first*a.first + a.second*a.second; + temp = sqrt(temp); + if (temp != 0) { + a.first = a.first/temp; + a.second = a.second/temp; + } +} +double calculateDotProduct(pair v, pair eV) { + return (v.first*eV.first + v.second*eV.second); +} +void calculateHeadVector(FlyObject fO, pair &headDirection) { + bool headDirectionBool = fO.getHeadIsInDirectionMAEV(); + pair fOMajorAxis = fO.getMajorAxisEV(); + if (headDirectionBool == true) { + headDirection.first = fOMajorAxis.first; + headDirection.second = fOMajorAxis.second; + } else { + headDirection.first = -fOMajorAxis.first; + headDirection.second = -fOMajorAxis.second; + } +} + +void determineHeadDirection(int fileCounter) { + + FrameInfo prevFI = fIVector[fileCounter-1]; + FrameInfo currentFI = fIVector[fileCounter]; + + vector pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + FlyObject pFLargeFO = pFOVector[0]; + FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + // calculate velocity + pair pFLCentroid = pFLargeFO.getCentroid(); + pair pFSCentroid = pFSmallFO.getCentroid(); + + pair cFLCentroid = cFLargeFO.getCentroid(); + pair cFSCentroid = cFSmallFO.getCentroid(); + + double velocityXLarge = static_cast (cFLCentroid.first) - static_cast (pFLCentroid.first); + double velocityYLarge = static_cast (cFLCentroid.second) - static_cast (pFLCentroid.second); + + double velocityXSmall = static_cast (cFSCentroid.first) - static_cast (pFSCentroid.first); + double velocityYSmall = static_cast (cFSCentroid.second) - static_cast (pFSCentroid.second); + + pair cLVV(velocityXLarge,velocityYLarge); + pair cSVV(velocityXSmall,velocityYSmall); + + cFLargeFO.setVelocityV(cLVV); + cFSmallFO.setVelocityV(cSVV); + + // normalize the velocity + cFLargeFO.normalizeVelocity(); + cFSmallFO.normalizeVelocity(); + + // determine the head direction for larger object + pair cLEV = cFLargeFO.getMajorAxisEV(); + // normalize the eigenvector + normalizeVector(cLEV); + pair identifiedCLEV; + + // the head from the previous frame + pair pLH = pFLargeFO.getHead(); + + // find the dot product with the previous frame head with current major axis both direction + double largerDotProdPrevHeadAndCLEV = calculateDotProduct(pLH, cLEV); + pair cLREV(-cLEV.first, -cLEV.second); + //double largerDotProdPrevHeadAndCLREV = calculateDotProduct(pLH, clREV); + if (largerDotProdPrevHeadAndCLEV >=0) { + cout << "Current eigen is in direction with the historical head"< newHead(newHeadFirst,newHeadSecond); + cFLargeFO.setHead(newHead); + cFLargeFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + identifiedCLEV = cLEV; + + } else { + cout << "Current eigen is in opposite direction of the historical head\n"; + // record this vector + double newHeadFirst = 0.1*cLREV.first + 0.9*pLH.first; + double newHeadSecond =0.1*cLREV.second+0.9*pLH.second; + pair newHead(newHeadFirst,newHeadSecond); + cFLargeFO.setHead(newHead); + cFLargeFO.setHeadIsInDirectionMAEV(false); + // set the identified EV to reverse direction of the current EV + identifiedCLEV = cLREV; + + } + + + // tracking the forward and backward movement of the fly object + double largeDotProd = calculateDotProduct(cLVV, identifiedCLEV); + cout << "largerDotProd "<= 0) { + cout<<"Larger Dot prod is positive"< cSEV = cFSmallFO.getMajorAxisEV(); + normalizeVector(cSEV); + pair identifiedCSEV; + // head from the previous frame + pair pSH = pFSmallFO.getHead(); + + // find the dot product with the previous frame head with current major axis both direction + double smallerDotProdPrevHeadAndCSEV = calculateDotProduct(pSH, cSEV); + pair cSREV(-cSEV.first, -cSEV.second); + if (smallerDotProdPrevHeadAndCSEV >=0) { + cout << "Current eigen is in direction with the historical head for the smaller fly object"< newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + cFSmallFO.setHeadIsInDirectionMAEV(true); + // set the identified EV to current EV + identifiedCSEV = cSEV; + } else { + cout << "Current eigen is in direction with the historical head for the smaller fly object"< newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + cFSmallFO.setHeadIsInDirectionMAEV(false); + // set the identified EV to reverse direction of current EV + identifiedCSEV = cSREV; + + } + + // to detect the swift change of head direction with 2 frames + // find previous frame's HD + pair previousHeadDirection; + calculateHeadVector(pFSmallFO, previousHeadDirection); + // find the current frame's HD + pair currentHeadDirection; + calculateHeadVector(cFSmallFO, currentHeadDirection); + + // calculate the dot product of previous head direction and current head direction. + double previousHDDotCurrentHD = calculateDotProduct(previousHeadDirection, currentHeadDirection); + if (previousHDDotCurrentHD < 0 and fileCounter > 1) { + + // now check if velocity direction( for the future historical velocity might be considered too) and + // current head direction dot product is also less than zero + // asssumption that current velocity and new head direction are towards same direction + // pair currentVelocity = cFSmallFO.getVelocityV(); + // double currentVVDotCurrentHD = calculateDotProduct(currentVelocity, currentHeadDirection); + // if (currentVVDotCurrentHD < 0) + { + // toggle current head direction + bool currentHeadDirectionBool = cFSmallFO.getHeadIsInDirectionMAEV(); + if (currentHeadDirectionBool == true) + cFSmallFO.setHeadIsInDirectionMAEV(false); + else { + cFSmallFO.setHeadIsInDirectionMAEV(true); + + } + + // update the historical head + // reset historical head to conform to the swift change in head direction + bool cFSmallFOFinalHeadDirectionBool = cFSmallFO.getHeadIsInDirectionMAEV(); + if (cFSmallFOFinalHeadDirectionBool == true) { + // record this vector for history + double newHeadFirst = cSEV.first; + double newHeadSecond = cSEV.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } else { + // record this vector for history + double newHeadFirst = cSREV.first; + double newHeadSecond = cSREV.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } + + + } + } + + /* retain historical head + bool cFSmallFOFinalHeadDirectionBool = cFSmallFO.getHeadIsInDirectionMAEV(); + if (cFSmallFOFinalHeadDirectionBool == true) { + // record this vector for history + double newHeadFirst = 0.1*cSEV.first + 0.9*pSH.first; + double newHeadSecond = 0.1*cSEV.second + 0.9*pSH.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } else { + // record this vector for history + double newHeadFirst = 0.1*cSREV.first + 0.9*pSH.first; + double newHeadSecond = 0.1*cSREV.second + 0.9*pSH.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } + */ + + double smallDotProd = calculateDotProduct(cSVV, cSEV); + + if (smallDotProd >= 0) { + cout << "Smaller dot product is positive"< updatedFOVector; + updatedFOVector.push_back(cFLargeFO); + updatedFOVector.push_back(cFSmallFO); + + currentFI.setFOVector(updatedFOVector); + + fIVector[fileCounter] = currentFI; + + // cout << "checking the update"< tempFIFOVector = tempFI.getFOVector(); + // FlyObject tempFILFO = tempFIFOVector[0]; + // FlyObject tempFISFO = tempFIFOVector[1]; + // pair tempFILFOVV = tempFILFO.getVelocityV(); + // cout << "Large object velocity vector "< tempFISFOVV = tempFISFO.getVelocityV(); + // cout << "Small object velocity vector "< x1 + if (x0 > x1) + { + int temp = x0; + x0 = x1; + x1 = temp; + temp = y0; + y0 = y1; + y1 = temp; + } + + int dx, dy, d, x, y, incrE, incrNE; + dx = x1 - x0; + dy = y1 - y0; + y = y0; + x = x0; + + // if they are on a vertical line + if (dx == 0) + { + if (y0 < y1) + for (int i = y0; i<=y1; i++) + lookAt(x,i,img); + else + for (int i = y1; i<=y0; i++) + lookAt(x,i,img); + return; + } + + // if they are on a horizontal line + if (dy == 0) + { + for (int i = x0; i<=x1; i++) + lookAt(i,y,img); + return; + } + + int dir = 0; + double m = double(dy)/double(dx); + if (m >= 1.0) + dir = 1; + else if ( (m < 0.0) && (m > -1.0) ) + dir = 2; + else if (m <= -1.0) + dir = 3; + + switch(dir) + { + // when slope m: 0< m <1 + case 0: + d = dy*2 - dx; + incrE = dy*2; + incrNE = (dy - dx)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (x= 1 + case 1: + d = dx*2 - dy; + incrE = dx*2; + incrNE = (dx - dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y= 0) + { + d+= incrE; + x++; + } + else + { + d+= incrNE; + x++; + y--; + } + + lookAt(x,y,img); + } + break; + // since initially we swapped the P0 with P1 so that P0 always holds values with smaller x cooridinate + // so we are sure that m<0 is the result of y0 being greater that + case 3: + d = dx*2 + dy; + incrE = dx*2; + incrNE = (dx + dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y>y1) + { + if (d <= 0) + { + d+= incrE; + y--; + } + else + { + d+= incrNE; + x++; + y--; + } + lookAt(x,y,img); + } + break; + } + + + if (inWhite) + { + // it is a fraction of the blob so dont count it + //len[dist(x_init,y_init,x,y)+1]++; + inWhite = false; + startedWhite = false; + } +}*/ + +bool isInFemaleBlob; +int maskImageHeight; +int maskImageWidth; +vector > bresenhamLine; + + +void putPixel(Image* maskImage, int x, int y, int color) { + + + // ignore the pixels outside of the image boundary + // if outside maximum y + if (y >= maskImageHeight) return; + if (y < 0) return; + if (x >= maskImageWidth) return; + if (x < 0) return; + + ColorMono currPixelColor = ColorMono(maskImage->pixelColor(x,y)); + + if (currPixelColor.mono() == true) { + isInFemaleBlob = true; + //isInBlackZone = false; + cout << "Hit the male object at "< temp(x,y); + bresenhamLine.push_back(temp); + } + + +} + + +int draw_line_bm(Image* maskImage, int x0, int y0, int x1, int y1) { + + isInFemaleBlob = false; + bresenhamLine.clear(); + + maskImageHeight = maskImage->rows(); + maskImageWidth = maskImage->columns(); + + int x, y; + int dx, dy; + dx = x1 - x0; + dy = y1 - y0; + double m = static_cast (dy)/static_cast (dx); + + int octant = -1; + if (( m >= 0 and m <= 1) and x0 < x1 ) { + cout << "Octant 1"< 1) and (y0 < y1)) { + cout << "Octant 2"<= -1) and (x0 > x1)) { + cout << "Octant 4"< 0 and m <=1) and (x0 > x1) ) { + cout << "Octant 5"< 1) and (y0 > y1) ) { + cout << "Octant 6"< y1) ) { + cout << "Octant 7"<= -1) and (x0 < x1) ) { + cout << "Octant 8"< x1 ) { + + if (d <= 0) { + + d = d + delW; + + x = x - 1; + + + } else { + d = d + delNW; + + + x = x - 1; + y = y + 1; + + } + + putPixel(maskImage,x, y, 4); + if (isInFemaleBlob == true) + return 1; + + } + + break; + + /*case 4: + + d = dx - 2*dy; + cout << "dx - 2*dy = "<<(dx-2*dy)< 0 "< x1) { + + if (d<=0) { + + d = d + delW; + + + x = x - 1; + + } else { + + d = d + delSW; + + + x = x - 1; + + y = y - 1; + + + } + + putPixel(maskImage,x, y, 5); + if (isInFemaleBlob == true) + return 1; + + + } + + break; + + + case 6: + + d = 2*dx - dy; + delS = 2*dx; + delSW = 2*(dx-dy); + + x = x0; + y = y0; + + + putPixel(maskImage,x,y,6); + if (isInFemaleBlob == true) + return 1; + + while (y > y1) { + + if (d<=0) { + + d = d + delS; + + y = y -1; + } + else { + + d = d + delSW; + + y = y -1; + x = x -1; + } + + putPixel(maskImage,x, y, 6); + if (isInFemaleBlob == true) + return 1; + + } + + break; + + case 7: + + d = 2*dx - dy; + delS = 2*dx; + delSE = 2*(dx-dy); + + x = x0; + y = y0; + + putPixel(maskImage,x,y,7); + if (isInFemaleBlob == true) + return 1; + + while (y > y1) { + + if (d<=0) { + d = d + delS; + y = y -1; + } + else { + d = d + delSE; + + + y = y - 1; + x = x + 1; + } + + putPixel(maskImage,x, y, 7); + if (isInFemaleBlob == true) + return 1; + + } + + break; + + case 8: + + d = 2*dy - dx; + delE = 2*dy; + delSE = 2*(dy - dx); + + x = x0; + y = y0; + + putPixel(maskImage,x,y,8); + if (isInFemaleBlob == true) + return 1; + + + while (x < x1) { + + if (d<=0) { + d = d + delE; + + x = x + 1; + } + else { + d = d + delSE; + + + y = y - 1; + x = x + 1; + } + cout << "putpixel "< newLocation, pair initLocation) { + + double temp = pow((newLocation.first - initLocation.first), 2.0) + pow((newLocation.second - initLocation.second), 2.0); + temp = sqrt(temp); + return temp; + +} + +int sequenceCondition(FrameInfo prevFI,FrameInfo currentFI) { + bool prevFIsSingleBlob = prevFI.getIsSingleBlob(); + bool currentFIsSingleBlob = currentFI.getIsSingleBlob(); + + if (prevFIsSingleBlob == false and currentFIsSingleBlob == true) + return STUCKING_TO_A_SINGLE_BLOB; + else if (prevFIsSingleBlob == true and currentFIsSingleBlob == false) + return SEPARATING_FROM_SINGLE_BLOB; + else + return -1; + +} + + + + +void drawTheSequence(int startIndex, int endIndex, int isFirst, bool singleBlob, bool unprocessed) { + + cout << "Should draw "<>fileName; + //inputFileName = "output/identified/"+fileName; + FrameInfo prevFI = fIVector[startIndex]; + cout << "Extracting information for image "<< fnVector[startIndex] << endl; + cout << "----------------------------------------\n"; + cout<>fileName; + fileName = fnVector[i]; + // inputFileName = "output/identified/"+fileName; + drawTheFlyObject(nextFI, fileName, isFirst, singleBlob, unprocessed); + } + + //inputFile.close(); +} + + + +void objectHeadDirection(FlyObject prevFO, FlyObject ¤tFO) { + + // take the head direction from the previous frame + pair pFHH = prevFO.getHead(); + + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdEVAndHH = calculateDotProduct(cFOEV, pFHH); + if (dotProdEVAndHH > 0 ) { + cout << "Current head is in direction with the Previous frame Head(which is the historical head from the maxDistIndex)"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(true); + + } else { + cout << "Current head is in direction with the previous frame head"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(false); + + } + + // copy the velocity vector from the previous frame + pair pFVV = prevFO.getVelocityV(); + currentFO.setVelocityV(pFVV); + +} + +void objectHeadDirection(FlyObject prevFO, FlyObject ¤tFO, bool prevFFOHD) { + + // take the head direction from the previous frame + //pair pFHH = prevFO.getHead(); + + pair pFHeadDir; + calculateHeadVector(prevFO, pFHeadDir); + + + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdEVAndHD = calculateDotProduct(cFOEV, pFHeadDir); + if (dotProdEVAndHD > 0 ) { + cout << "Current head is in direction with the Previous frame Head direction"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(true); + */ + currentFO.setHead(cFOEV); + currentFO.setHeadIsInDirectionMAEV(true); + + } else { + cout << "Current head is in reverse direction with the previous frame head direction"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(false); + */ + currentFO.setHead(cFOREV); + currentFO.setHeadIsInDirectionMAEV(false); + + } + + // copy the velocity vector from the previous frame + if (prevFFOHD == true) { + pair pFVV = prevFO.getVelocityV(); + currentFO.setVelocityV(pFVV); + } + +} + +void objectHeadDirection(FlyObject ¤tFO, int saveEV=0) { + + // get the velocity vector + pair cFOVV = currentFO.getVelocityV(); + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdVVAndEV = calculateDotProduct(cFOEV, cFOVV); + if (dotProdVVAndEV > 0 ) { + cout << "Current head is in direction with the velocity vector\n"; + // set the head to the eigen vector + currentFO.setHead(cFOEV); + currentFO.setHeadIsInDirectionMAEV(true); + + if (saveEV == 1) { + cout << "saving the eigen vector for the first object"< zero(0.0,0.0); + + if (saveEV == 1) { + cout << "saving the zero eigen vector for the first object"< cFV) { + +// // get the velocity vector +// pair cFOVV = currentFO.getVelocityV(); + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdVVAndEV = calculateDotProduct(cFOEV, cFV); + if (dotProdVVAndEV > 0 ) { + cout << "Current head is in direction with the vector used\n"; + // set the head to the eigen vector + currentFO.setHead(cFOEV); + currentFO.setHeadIsInDirectionMAEV(true); + } else { + cout << "Current head is in reverse direction with the vector used"< &velDirectionF, pair&velDirectionS) { + + // find the average velocity vector + cout << "Finding average velocity vector from "< pFOVector = prevFI.getFOVector(); + FlyObject pFFirstFO = pFOVector[0]; + FlyObject pFSecondFO = pFOVector[1]; + + + FrameInfo currentFI = fIVector[end]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + // get the summation of (velocity direction from start frame to half of the frames in this interval) + pair cFFCentroid = cFFirstFO.getCentroid(); + pair cFSCentroid = cFSecondFO.getCentroid(); + + + pair pFFCentroid = pFFirstFO.getCentroid(); + pair pFSCentroid = pFSecondFO.getCentroid(); + + + int velXFirst = cFFCentroid.first - pFFCentroid.first; + int velYFirst = cFFCentroid.second - pFFCentroid.second; + + int velXSecond = cFSCentroid.first - pFSCentroid.first; + int velYSecond = cFSCentroid.second - pFSCentroid.second; + + cout << "Velocity vector"< cFVV(static_cast (velXFirst), static_cast (velYFirst)); + pair cSVV(static_cast (velXSecond), static_cast (velYSecond)); + + velDirectionF = cFVV; + velDirectionS = cSVV; + + +} + +double getSpeed(pair vector) { + double value = vector.first*vector.first + vector.second*vector.second; + value = sqrt(value); + + return value; + + +} + +void velocityDirections(int stIndex, int endIndex) { + + + velocityDirectionsF.clear(); + velocityDirectionsS.clear(); + + //int intervalLength = 5; + cout << "Initial velocity direction calculation"< velDirectionF; + pair velDirectionS; + + // get the velocity + velocityDirection(i,i+1, velDirectionF, velDirectionS); + + // extract the fly object for update + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + // get the speed of the velocity + double speedF = getSpeed(velDirectionF); + cFFirstFO.setSpeed(speedF); + double speedS = getSpeed(velDirectionS); + cFSecondFO.setSpeed(speedS); + + // save the velocity direction + normalizeVector(velDirectionF); + velocityDirectionsF.push_back(velDirectionF); + normalizeVector(velDirectionS); + velocityDirectionsS.push_back(velDirectionS); + + cout << "Normalized velocity"< updatedFOVector; + updatedFOVector.push_back(cFFirstFO); + updatedFOVector.push_back(cFSecondFO); + currentFI.setFOVector(updatedFOVector); + cout << "Updating the frame "< > velocityDirs, int &startIndex, int &endIndex) { + + int positiveVelSeqSize = 0; + int flag = false; + int maxSeqSize = 0; + int st = 0; + for (int j=0; j prevVel = velocityDirs[j]; + pair currVel = velocityDirs[j+1]; + + double dotProd = calculateDotProduct(prevVel, currVel); + + if( dotProd > 0 && flag == false) { + st = j; + positiveVelSeqSize++; + flag = true; + //cout << "In first if positiveSize "< 0 && flag == true) { + positiveVelSeqSize++; + //cout << "In second if positive "< maxSeqSize) { + maxSeqSize = positiveVelSeqSize; + startIndex = st; + endIndex = st+positiveVelSeqSize; + //cout << "maxseq updated \npositiveSize "< prevVel = velocityDirs[j]; + if (prevVel.first != 0 || prevVel.second != 0) { + startIndex = j; + endIndex = j; + zero = true; + break; + } + } + + if (zero != true) { + cout << "All directions zero"< cFOVector = currentFI.getFOVector(); + FlyObject cFFO = cFOVector[object-1]; + pair cFVV = cFFO.getVelocityV(); + //cout << "Velocity before normalization "< cFEV = cFFO.getMajorAxisEV(); + cout << "Eigen vector before normalization "< maxVelValue) { + maxVelValue = velValue; + maxVelValIndex = i; + } + + + } + + cout << "Maximum speed is chosen for for frame "< cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + cout << "Setting representative head direction of frame "< tempVelocity = cFSecondFO.getVelocityV(); + //cout << "Velocity was "< intervalLength) { + + cout << "Propagating in the longest positive frame interval starting the direction from the frame "< pFOVector = prevFI.getFOVector(); + FlyObject pFFirstFO = pFOVector[0]; + FlyObject pFSecondFO = pFOVector[1]; + + for (int i=t+1; i<=e; i++) { + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< tempVelocity = cFSecondFO.getVelocityV(); + //cout << "Velocity was "<=s; i--) { + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< pFOVector = prevFI.getFOVector(); + FlyObject pFFirstFO = pFOVector[0]; + FlyObject pFSecondFO = pFOVector[1]; + + if (e < origEnd) { + cout << "Propagating front from "< cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< origStart ) { + cout << "Propagating down wards from "<=origStart; i--) { + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< z = velocityDirectionsF[a]; + pair zEV = evDirectionF[a]; + cout <<" First object velocity = "< w = velocityDirectionsS[a]; + pair wEV = evDirectionS[a]; + + cout << "Second object velocity = "< pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + // bool currentlySingleBlob = currentFI.getIsSingleBlob(); + + // if (currentlySingleBlob == false) { + + FlyObject pFLargeFO = pFOVector[0]; + FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + double distLTL = euclideanDist(pFLargeFO, cFLargeFO); + cout << "previousFirst to currentFirst "< aCentroid = a.getCentroid(); + pair bCentroid = b.getCentroid(); + + //cout << "Distance from "< (aCentroid.first)-static_cast (bCentroid.first)), 2.0); + double y2 = pow((static_cast (aCentroid.second)-static_cast (bCentroid.second)), 2.0); + double dist = sqrt((x2+y2)); + cout << dist< cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + double currDist = euclideanDist(cFFirstFO, cFSecondFO); + + if (currDist > maxDist) { + maxDist = currDist; + maxDistIndex = i; + cout << "New max distance" << maxDist << " at the frame number "<=startOfATrackSequence; i--) { + FrameInfo currentFI = fIVector[i]; + objectCorrespondence(prevFI, currentFI); + + // update the fIVector + fIVector[i] = currentFI; + prevFI = currentFI; + + } + cout << "Assigning object correspondance between maxDistIndex+1 "<<(maxDistIndex+1)<<" to endIndex "< fOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = fOVector[0]; + FlyObject cFSecondFO = fOVector[1]; + + sequenceFirstAverage += cFFirstFO.getArea(); + sequenceSecondAverage += cFSecondFO.getArea(); + + cout << "For frame number "< sequenceSecondAverage) { + foutLPS << "First object is the female and avg size is "< sequenceSecondAverage) { + cout << "First sequence is the object A"< " << endl; // input file contains name of the + // input image files + return -1; + } + + MagickCore::SetMagickResourceLimit(MagickCore::MemoryResource, 1536); + MagickCore::SetMagickResourceLimit(MagickCore::MapResource, 4092); + + string ifns(argv[1]); + inputFileName = ifns; +// cout << "input file name is "< > shape; + vector tempFOV; + + // to find the new head direction + bool currentlyCalculatingHead = true; + + + + while (inputFile>>fileName) { + +// Image* img = new Image(argv[1]);// = new Image(argv[1]); + + int fi = fileName.find("_"); + // current sequence numbers spans from 0 - 18019, so 5 digits are needed + int span = 5; + string tempString = fileName.substr(fi+1,span); + int frameCounter = atoi(tempString.c_str()); + //cout << frameCounter<columns(),height = img->rows(); + diagLength= static_cast ( sqrt( (height*height) + (width*width) ) ); + //cout << "Diagonal length is "< 0 ) + { + + // cout << "size of the object is: " << s < eigenVal = covariantDecomposition(shape); + + { +// objCount++; + + double velocity_x=0.0, velocity_y=0.0; + // save the object information + FlyObject tempFO(s, + pair (eigenVal[6], eigenVal[7]), + pair (eigenVal[4], eigenVal[5]), + pair (velocity_x, velocity_y), + false, + pair (eigenVal[4], eigenVal[5]), + 0.0); + tempFOV.push_back(tempFO); + + } + } + + } + } + + + delete img; + + delete residual; + +// cout<<"Sorting the objects according to size"< 1) { + + FrameInfo prevFI = fIVector[fileCounter-1]; + + int seqCond = sequenceCondition(prevFI, currentFI); + + if (seqCond == STUCKING_TO_A_SINGLE_BLOB) { + + endOfATrackSequence = fileCounter-1; + + // save the index for printing the one object later on + startOfAOneObject = fileCounter; + + } else if (seqCond == SEPARATING_FROM_SINGLE_BLOB) { + + startOfATrackSequence = fileCounter; + + // draw the single blob sequence + endOfAOneObject = fileCounter - 1; + cout << "Only one object StartIndex "<= 15 to be a single blob state, so pass false parameter. + drawTheSequence(startOfAOneObject, endOfAOneObject,0, true, false); + startOfAOneObject = -1; + endOfAOneObject = -1; + //startOfATrackSequence = 1; + sequenceSize = 1; + } + + + if(seqCond == STUCKING_TO_A_SINGLE_BLOB) { + cout << "StartIndex "< 15 when the size of the sequence is > 16. + if ((endOfATrackSequence - startOfATrackSequence) >=15 ) { + processASequence(startOfATrackSequence, endOfATrackSequence); + cout << "Done processing"< 15 ) { + + cout << "Last sequence that does not stick to a single blob status startIndex "<(totalMaleLookingAtFemale)/static_cast (fileCounter-totalUnprocessedFrame); + //percentageLookingAt *= 100.0; + + //double percentageSingleBlob = static_cast(totalSingleBlob)/static_cast(fileCounter); + //percentageSingleBlob *= 100.0; + + + // new calculation of percentage look at should consider only those frames where the flies are separated + double percentageLookingAt = static_cast(totalMaleLookingAtFemale+totalFemaleLookingAtMale)/static_cast(totalSeparated); + percentageLookingAt *= 100.0; + + double percentageSingleBlob = static_cast(totalSingleBlob)/static_cast(fileCounter); + percentageSingleBlob *= 100.0; + + foutSt<<"Total number of single blob "<=0; i--) { + + pair tempP = bresenhamLine[i]; + + int x = tempP.first; + + int y = tempP.second; + + ColorMono currPixelColor = ColorMono(maskImage->pixelColor(x,y)); + + if (currPixelColor.mono() == true) { + cout << "Hit the target fly"<<" "<columns(); + int height = image->rows(); + + char buffer[100]; + sprintf(buffer,"%ix%i",width,height); + + // the residual image should be newed + residual = new Image(buffer, "black"); + + bool found = false; + vector > foundShape; + vector > shape; + cout<<"Detecting the male object for finding the start point"< 0 ) + { + + cout << "size of the object is: " << s < tmpCentroid = getCentroid(shape); + foutDebugCen<<"tmpCentroid.first, tmpCentroid.second = ("< point = shape[l]; + cout << "point.first, point.second "< eigenVal = covariantDecomposition(foundShape); + + int x0, y0, x1, y1; + + x0 = cen_x; + y0 = cen_y; + + double ev_x; + double ev_y; + + if (eVDirection == true) { + ev_x = static_cast (x0) + static_cast (diagLength)*eigenVal[4]; + ev_y = static_cast (y0) + static_cast (diagLength)*eigenVal[5]; + } + else { + ev_x = static_cast (x0) - static_cast (diagLength)*eigenVal[4]; + ev_y = static_cast (y0) - static_cast (diagLength)*eigenVal[5]; + + } + + x1 = static_cast (ev_x); + y1 = static_cast (ev_y); + + cout<<"Endpoint: centroid (x0,y0)==("< point = foundShape[i]; + maskImage->pixelColor(point.first, point.second,"white"); + } + + + int hits = draw_line_bm(maskImage, x1, y1, x0, y0); + + + //maskImage->strokeColor("red"); + //maskImage->draw(DrawableLine(x1, y1, x0, y0)); + //maskImage->write("test.png"); + + cout<<"BresenhamLine size is "< 0) { + pair temp = bresenhamLine[bresenhamLine.size()-1]; + cout << "Finding the starting point: Hits source at "<pixelColor(temp.first, temp.second); + //maleSP_x = prev_x; + //maleSP_y = prev_y; + + if (c.mono() == true) { + cout << "start point from the source object should be black"< 0) { + pair temp = bresenhamLine[bresenhamLine.size()-1]; + cout << "Finding the starting point: Hits source at "<pixelColor(temp.first, temp.second); + //maleSP_x = prev_x; + //maleSP_y = prev_y; + + if (c.mono() == true) { + cout << "start point from the source object should be black"< fOVector = currentFI.getFOVector(); + + if (singleBlob == false and unprocessed == false) { + + FlyObject maleFO = fOVector[isFirst]; + FlyObject femaleFO = fOVector[1-isFirst]; + + pair maleCentroid = maleFO.getCentroid(); + + pair maleMajorAxisEV = maleFO.getMajorAxisEV(); + + bool maleEVDir = maleFO.getHeadIsInDirectionMAEV(); + + pair femaleCentroid = femaleFO.getCentroid(); + + pair femaleMajorAxisEV = femaleFO.getMajorAxisEV(); + + bool femaleEVDir = femaleFO.getHeadIsInDirectionMAEV(); + + + // 1. finding the distance between the centroids + double tempDist = pow(static_cast(maleCentroid.first - femaleCentroid.first),2) + pow(static_cast(maleCentroid.second-femaleCentroid.second),2); + + tempDist = sqrt(tempDist); + + // round the function + unsigned int dist = roundT(tempDist); + + centroidDistanceMap[dist] = centroidDistanceMap[dist] + 1; + + foutSt << "Centroid distance "< maleHeadDir; + pair femaleHeadDir; + + if (maleEVDir == true) { + maleHeadDir = maleMajorAxisEV; + foutSt<<"Male Head In direction of ev"<(dp)); + cout<<"Angle in radian "<(PI); + cout<<"Angle in deg "<(roundT(static_cast(deg))); + cout<<"Angle after rounding "< fOVector = currentFI.getFOVector(); + + cout << "While drawing it found objects = "< femaleCentroid = fO.getCentroid(); + bool eVDirectionFemale = fO.getHeadIsInDirectionMAEV(); + /////////// + + FlyObject currentFO = fOVector[isFirst]; + + pair centroid = currentFO.getCentroid(); + + bool eVDirection = currentFO.getHeadIsInDirectionMAEV(); + + // debug: + cout<<"Calling the findTheStartPoint() function"< centroid = currentFO.getCentroid(); + pair majorAxisEV = currentFO.getMajorAxisEV(); + + bool eVDirection = currentFO.getHeadIsInDirectionMAEV(); + + double ev_x, ev_y; + + // draw the female when tracked by the male fly + if (isHitting == 1) { + if (n != isFirst and n == 0) { + ///img->strokeColor("red"); + + //img->draw(DrawableRectangle(centroid.first - 6, centroid.second - 6, centroid.first + 6, centroid.second + 6)); + } else if (n != isFirst and n==1) { + + // img->strokeColor("red"); + + //img->draw(DrawableRectangle(centroid.first - 6, centroid.second - 6, centroid.first + 6, centroid.second + 6)); + } + } + + + // draw the male when tracked by the female fly + // this situation will draw the YELLOW circle as well as the ORANGE rectangle around the + // male object. This double color ensure that this situation is incorrectly detects the male + // as female. Because our assumption is that female never tracks the male. So the actual female + // tracking male will occur very insignificant times. + if (isHittingFemaleToMale == 1) { + if (n == isFirst and n == 0) { + // img->strokeColor("Red"); + + // img->draw(DrawableRectangle(centroid.first - 6, centroid.second - 6, centroid.first + 6, centroid.second + 6)); + } else if (n == isFirst and n==1) { + + //img->strokeColor("Red"); + + //img->draw(DrawableRectangle(centroid.first - 6, centroid.second - 6, centroid.first + 6, centroid.second + 6)); + } + } + + + // draw the Axis direction + if (eVDirection == true) { + //ev_x = static_cast(centroid.first) + 50.0 * majorAxisEV.first; + //ev_y = static_cast(centroid.second) + 50.0 * majorAxisEV.second; + //img->strokeColor("green"); + //img->draw( DrawableLine( centroid.first, centroid.second, static_cast(ev_x), static_cast(ev_y) )); + } else { + //ev_x = static_cast(centroid.first) - 50.0 * majorAxisEV.first; + //ev_y = static_cast(centroid.second) - 50.0 * majorAxisEV.second; + //img->strokeColor("green"); + //img->draw( DrawableLine( centroid.first, centroid.second, static_cast(ev_x), static_cast(ev_y) )); + } + + // draw the velocity vector + + /*img->strokeColor("blue"); + pair velocityV = currentFO.getVelocityV(); + ev_x = static_cast(centroid.first) + 30.0 * velocityV.first; + ev_y = static_cast(centroid.second) + 30.0 * velocityV.second; + img->draw(DrawableLine( centroid.first, centroid.second, static_cast(ev_x), static_cast(ev_y) )); + */ + + /* + // draw the historical head vector + img->strokeColor("white"); + pair headV = currentFO.getHead(); + ev_x = static_cast (centroid.first) + 25.0*headV.first; + ev_y = static_cast (centroid.second) + 25.0*headV.second; + img->draw( DrawableLine(centroid.first, centroid.second, static_cast (ev_x), static_cast (ev_y)) ); + */ + + //draw the object tracking circle + if (n == isFirst and n==0) { + cout << "Tracking the n = "<strokeColor("yellow"); + //img->draw(DrawableCircle(centroid.first, centroid.second, centroid.first+5, centroid.second)); + //img->pixelColor(prev_x, prev_y, "red"); + + } else if ( n == isFirst and n==1) { + + cout << "Tracking the "<strokeColor("yellow"); + //img->fillColor("none"); + //img->pixelColor(prev_x, prev_y, "red"); + //img->draw(DrawableCircle(centroid.first, centroid.second, centroid.first+5, centroid.second)); + + + } + + + } + } + + // overwrite the file now with axis + // img->write(inputFileName.c_str()); + + // when do not want to identify on the original comment below line and uncomment the above one + //img->write(outputFileName.c_str()); + delete maskImage; + if (isHitting == 1 || isHittingFemaleToMale == 0) + calculateStatistics(currentFI, fileName, isFirst, singleBlob, true, false, unprocessed); + else if (isHitting == 0 || isHittingFemaleToMale == 1) + calculateStatistics(currentFI, fileName, isFirst, singleBlob, false, true, unprocessed); + else if (isHitting == 1 || isHittingFemaleToMale == 1) + calculateStatistics(currentFI, fileName, isFirst, singleBlob, true, true, unprocessed); + else + calculateStatistics(currentFI, fileName, isFirst, singleBlob, false, false, unprocessed); + +} + +void findObj(Image* img, int x, int y, vector > & 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 > & 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 p; + p.first = x; + p.second = y; + obj.push_back(p); + +// if (obj.size() > barrier) { +// //cout<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 > & 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 p; + p.first = x; + p.second = y; + obj.push_back(p); + +// if (obj.size() > barrier) { +// //cout<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 getCentroid(vector > & points) +{ + pair centroid; + centroid.first = 0; + centroid.second = 0; + + for (unsigned int i = 0; i covariantDecomposition(vector > & points) +{ + unsigned int i,j,k; + pair centroid = getCentroid(points); + vector retval; + + gsl_matrix* matrice = gsl_matrix_alloc(2, 2); + + double sumX2 = 0, sumXY = 0, sumY2 = 0; + for (k = 0; ksize; i++) + retval.push_back(gsl_vector_get(eigenVal, i)); + + for (j = 0; jsize2; j++) + for (i = 0; isize1; i++) + retval.push_back(gsl_matrix_get(eigenVec, i, j)); + + retval.push_back(static_cast(centroid.first)); + retval.push_back(static_cast (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; + +} diff --git a/fly-tools/FlyTrackingNoImages b/fly-tools/FlyTrackingNoImages new file mode 100755 index 0000000..b61fcd9 Binary files /dev/null and b/fly-tools/FlyTrackingNoImages differ diff --git a/fly-tools/FrameInfo.cpp b/fly-tools/FrameInfo.cpp new file mode 100644 index 0000000..d00e4cd --- /dev/null +++ b/fly-tools/FrameInfo.cpp @@ -0,0 +1,72 @@ +/* + * FrameInfo.cpp + * + * + * Created by Md. Alimoor Reza on 6/26/10. + * Copyright 2010 Drexel University. All rights reserved. + * + */ + +#include "FrameInfo.h" + +FrameInfo::FrameInfo(int frameNo, vector fOVector, bool isSingleBlob) { + this->frameNo = frameNo; + this->fOVector = fOVector; + this->isSingleBlob = isSingleBlob; + +} +FrameInfo::FrameInfo(const FrameInfo &f) { + this->frameNo = f.getFrameNo(); + this->fOVector = f.getFOVector(); + this->isSingleBlob = f.getIsSingleBlob(); + +} +int FrameInfo::getFrameNo() const { + return frameNo; +} +bool FrameInfo::getIsSingleBlob() const { + return isSingleBlob; +} +vector FrameInfo::getFOVector() const{ + return fOVector; +} + +void FrameInfo::setFrameNo(int fn) { + this->frameNo = fn; +} +void FrameInfo::setIsSingleBlob(bool isSingleBlob) { + this->isSingleBlob = isSingleBlob; +} +void FrameInfo::setFOVector(vector fov) { + this->fOVector = fov; + +// cout << "setting fov \n"; +// FlyObject a = fOVector[0]; +// pair av = a.getVelocityV(); +// cout << "velocity "< bv = b.getVelocityV(); +// cout << "velocity "< 1) { + cout << "swapping\n"; + FlyObject a = fOVector[0]; + fOVector[0] = fOVector[1]; + fOVector[1] = a; + } +} +void FrameInfo::output(ostream &out) { + out<<"FrameNo : "< +#include +#include +#include "FlyObject.h" +using namespace std; + +class FrameInfo { +public: + FrameInfo(int frameNo, vector fOVector, bool isSingleBlob); + FrameInfo(const FrameInfo &f); + int getFrameNo() const; + bool getIsSingleBlob() const; + vector getFOVector() const; + + void setFrameNo(int fn); + void setIsSingleBlob(bool isSingleBlob); + void setFOVector(vector fov); + void swapTheFlyObject(); + void output(ostream &out); + + +private: + bool isSingleBlob; + int frameNo; + vector fOVector; +}; diff --git a/fly-tools/FurthestPointAlongMajorAxis.cpp b/fly-tools/FurthestPointAlongMajorAxis.cpp new file mode 100644 index 0000000..8c8a086 --- /dev/null +++ b/fly-tools/FurthestPointAlongMajorAxis.cpp @@ -0,0 +1,3084 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "FrameInfo.h" + +using namespace Magick; +using namespace std; + +const double PI = atan(1.0)*4.0; +const double FACTOR_EIGEN = 100; +const int STUCKING_TO_A_SINGLE_BLOB = 1; +const int SEPARATING_FROM_SINGLE_BLOB = 2; + + +Image* residual; + +vector fOVector; +vector fIVector; +// temporary store of the frames for finding initial head direction +vector fIForHeadVector; +// since for the first time the head is automatically need to be set +int startIndexToFindNewHead=0; +int endIndexToFindNewHead; +double initLargeLocationX=-1; +double initLargeLocationY=-1; +double initSmallLocationX=-1; +double initSmallLocationY=-1; +pair largeCollisionBeforeDirection; +pair smallCollisionBeforeDirection; + +// start location in the vector where two object get together +int startIndexSingleBlob=-1; +int endIndexSingleBlob =-1; + +// indices of the new set of sequences +int startOfATrackSequence = 0; // fix it later on inside main +int endOfATrackSequence=-1; +int sequenceSize=1; + +// indices for printing the one object frames +int startOfAOneObject = -1; +int endOfAOneObject = -1; + +vector fnVector; +string inputFileName; + +// GLOBAL PATHS +string maskImagePath; +string origImagePath; +string finalOutputPath; + +vector > velocityDirectionsF; +vector > velocityDirectionsS; + +pair avgVelocityF; +pair avgVelocityS; + +pair overAllVelocityF; +pair overAllVelocityS; + +vector > evDirectionF; +vector > evDirectionS; + +void initSequence(){ + + startOfATrackSequence = -1; + endOfATrackSequence = -1; + sequenceSize = 1; + + +} + +ostream &operator<<(ostream &out, FlyObject & fO) { + fO.output(out); + return out; +} + +ostream &operator<<(ostream &out, FrameInfo & fI) { + fI.output(out); + return out; +} + +void bubbleSort(vector & fov) { + + //FlyObject a,b,c; + for(int i=1; i > & shape ,bool eightCon=true, bool colorLookingFor=true); +void eightConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +void fourConnObj(Image* img, int x, int y, vector > & obj, bool color=true); +vector covariantDecomposition(vector > & points); +pair getCentroid(vector > & points); +bool isInterface(Image* orig, unsigned int x, unsigned int y); +void writeFrameImage(int fn, string imS); +void drawTheFlyObject(FrameInfo currentFI, string fileName, int isFirst, bool singleBlob=false); +void drawTheSequence(int startIndex, int endIndex, int isFirst, bool singleBlob = false); +double euclideanDist(FlyObject a, FlyObject b); +bool identifyFlyObjectFromFrameToFrame(FrameInfo prevFI, FrameInfo& currentFI, bool gotRidOfSingleBlob=false) ; +int roundT(double v) {return int(v+0.5);} +void determineHeadDirection(int fileCounter); + + +void normalizeVector(pair &a); +double calculateDotProduct(pair v, pair eV); +void calculateHeadVector(FlyObject fO, pair &headDirection); + +/* +void lookAt(int x, int y, jzImage& img) +{ + int imageWidth =img.width(); + int imageHeight = img.height(); + // if current pixel is white + if (img.pixelColor(x,y) == 1) { + // if it was the first white pixel then + if (!inWhite) + { + // check whether it started as a white pixel; if it started as white pixel then this line segment should + // not be counted. because it is part of the larger white blob + if ( (x == 0 || x == (imageWidth -1) ) || ( y == 0 || y == (imageHeight - 1))) + startedWhite = true; + inWhite = true; + x_init = x; + y_init = y; + } + // if we are on a white region + //else { + + //} + + } + + if (img.pixelColor(x,y) == 0) { + // if we are going through a white line and reached the black pixel + if (inWhite) + { + // if the line started as a white pixel then discard the line fragment + if (startedWhite == false) { + int val = roundT(dist(x_init, y_init, x, y)); + //cout<<" val is = " << val << " at (x0,y0) = " << x_init << "," << y_init << " to (x1,y1) = " << x << "," << y < x1 + if (x0 > x1) + { + int temp = x0; + x0 = x1; + x1 = temp; + temp = y0; + y0 = y1; + y1 = temp; + } + + int dx, dy, d, x, y, incrE, incrNE; + dx = x1 - x0; + dy = y1 - y0; + y = y0; + x = x0; + + // if they are on a vertical line + if (dx == 0) + { + if (y0 < y1) + for (int i = y0; i<=y1; i++) + lookAt(x,i,img); + else + for (int i = y1; i<=y0; i++) + lookAt(x,i,img); + return; + } + + // if they are on a horizontal line + if (dy == 0) + { + for (int i = x0; i<=x1; i++) + lookAt(i,y,img); + return; + } + + int dir = 0; + double m = double(dy)/double(dx); + if (m >= 1.0) + dir = 1; + else if ( (m < 0.0) && (m > -1.0) ) + dir = 2; + else if (m <= -1.0) + dir = 3; + + switch(dir) + { + // when slope m: 0< m <1 + case 0: + d = dy*2 - dx; + incrE = dy*2; + incrNE = (dy - dx)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (x= 1 + case 1: + d = dx*2 - dy; + incrE = dx*2; + incrNE = (dx - dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y= 0) + { + d+= incrE; + x++; + } + else + { + d+= incrNE; + x++; + y--; + } + + lookAt(x,y,img); + } + break; + // since initially we swapped the P0 with P1 so that P0 always holds values with smaller x cooridinate + // so we are sure that m<0 is the result of y0 being greater that + case 3: + d = dx*2 + dy; + incrE = dx*2; + incrNE = (dx + dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y>y1) + { + if (d <= 0) + { + d+= incrE; + y--; + } + else + { + d+= incrNE; + x++; + y--; + } + lookAt(x,y,img); + } + break; + } + + + if (inWhite) + { + // it is a fraction of the blob so dont count it + //len[dist(x_init,y_init,x,y)+1]++; + inWhite = false; + startedWhite = false; + } +}*/ + +bool isInMaleBlob; +bool isInBlackZone; +bool isInFemaleBlob; +bool isAtTheBoundary; +int maskImageHeight; +int maskImageWidth; +void putPixel(Image* maskImage, int x, int y, int color) { + + ColorMono currPixelColor = ColorMono(maskImage->pixelColor(x,y)); + + // if current pixel is still inside the male fly do nothing + // if current pixel is black and prior to that it was inside the male fly then enter the black zone + if (currPixelColor.mono() == false and isInMaleBlob == true) { + isInBlackZone = true; + isInMaleBlob = false; + cout << "Enters black zone "<rows(); + maskImageWidth = maskImage->columns(); + + int x, y; + int dx, dy; + dx = x1 - x0; + dy = y1 - y0; + double m = static_cast (dy)/static_cast (dx); + + int octant = -1; + if (( m >= 0 and m <= 1) and x0 < x1 ) { + cout << "Octant 1"< 1) and (y0 < y1)) { + cout << "Octant 2"<= -1) and (x0 > x1)) { + cout << "Octant 4"< 0 and m <=1) and (x0 > x1) ) { + cout << "Octant 5"< 1) and (y0 > y1) ) { + cout << "Octant 6"< y1) ) { + cout << "Octant 7"<= -1) and (x0 < x1) ) { + cout << "Octant 8"< 0\n"; + y = y + 1; + x = x - 1; + } + putPixel(maskImage,x, y, 3); + if (isInFemaleBlob == true) + return 1; + else if (isAtTheBoundary == true) + return 2; + + + } + + break; + + case 4: + + d = -dx + 2*dy; + //cout << "-dx + 2*dy = "<<(-dx+2*dy)< 0 "< 0 "<=x1) { + + if (d<=0) { + + //cout << "Going W since d<0"< 0"<=y1) { + + if (d<=0) { + d = d + delS; + y = y -1; + } + else { + d = d + delSW; + y = y -1; + x = x -1; + } + + putPixel(maskImage,x, y, 6); + if (isInFemaleBlob == true) + return 1; + else if (isAtTheBoundary == true) + return 2; + + + } + + break; + + case 7: + + d = 2*dx - dy; + delS = 2*dx; + delSE = 2*(dx-dy); + + x = x0; + y = y0; + + putPixel(maskImage,x,y,7); + if (isInFemaleBlob == true) + return 1; + else if (isAtTheBoundary == true) + return 2; + + + while (y>=y1) { + + if (d<=0) { + d = d + delS; + y = y -1; + } + else { + d = d + delSE; + y = y - 1; + x = x + 1; + } + + putPixel(maskImage,x, y, 7); + if (isInFemaleBlob == true) + return 1; + else if (isAtTheBoundary == true) + return 2; + + } + + break; + + case 8: + + d = 2*dy - dx; + delE = 2*dy; + delSE = 2*(dy - dx); + + x = x0; + y = y0; + + putPixel(maskImage,x,y,8); + if (isInFemaleBlob == true) + return 1; + else if (isAtTheBoundary == true) + return 2; + + + while (x<=x1) { + + if (d<=0) { + d = d + delE; + x = x + 1; + } + else { + d = d + delSE; + y = y - 1; + x = x + 1; + } + + putPixel(maskImage,x, y, 8); + if (isInFemaleBlob == true) + return 1; + else if (isAtTheBoundary == true) + return 2; + + + } + + break; + + default: + cout << "No octant which should be a bug\n"; + exit(0); + break; + + + } + + return 0; + +} + + +void fillResidualWithObj(vector > & obj, ColorRGB c) +{ + for (unsigned int i = 0; ipixelColor(obj[i].first, obj[i].second, c); +} +double euclideanDist(pair newLocation, pair initLocation) { + + double temp = pow((newLocation.first - initLocation.first), 2.0) + pow((newLocation.second - initLocation.second), 2.0); + temp = sqrt(temp); + return temp; + +} +bool calculateDisplacement(FrameInfo prevFI, FrameInfo currentFI) { + + vector pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + FlyObject pFLargeFO = pFOVector[0]; + FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + pair cFLargeCentroid = cFLargeFO.getCentroid(); + pair initLarge(initLargeLocationX, initLargeLocationY); + + double largeDist = euclideanDist(cFLargeCentroid, initLarge); + cout << "Large object current position "< cFSmallCentroid = cFSmallFO.getCentroid(); + pair initSmall(initSmallLocationX, initSmallLocationY); + + double smallDist = euclideanDist(cFSmallCentroid, initSmall); + cout << "Small object current position "<= 20 || smallDist >=20) +// return true; +// else +// return false; + if (largeDist >= 5 || smallDist >=5) + return true; + else + return false; + + +} + +int sequenceCondition(FrameInfo prevFI,FrameInfo currentFI) { + bool prevFIsSingleBlob = prevFI.getIsSingleBlob(); + bool currentFIsSingleBlob = currentFI.getIsSingleBlob(); + + if (prevFIsSingleBlob == false and currentFIsSingleBlob == true) + return STUCKING_TO_A_SINGLE_BLOB; + else if (prevFIsSingleBlob == true and currentFIsSingleBlob == false) + return SEPARATING_FROM_SINGLE_BLOB; + else + return -1; + +} + + +bool isStuckToSingleBlob(FrameInfo prevFI,FrameInfo currentFI) { + bool prevFIsSingleBlob = prevFI.getIsSingleBlob(); + bool currentFIsSingleBlob = currentFI.getIsSingleBlob(); + + if (prevFIsSingleBlob == false and currentFIsSingleBlob == true) + return true; + else + return false; + +} + + +bool isStartOfNewHeadDirectionSearch(FrameInfo prevFI, FrameInfo currentFI) { + + bool prevFIsSingleBlob = prevFI.getIsSingleBlob(); + bool currentFIsSingleBlob = currentFI.getIsSingleBlob(); + + if (prevFIsSingleBlob == true and currentFIsSingleBlob == false) + return true; + else + return false; + +} + +void separateObjectAndInitForNewHeadCalculation(FrameInfo prevFI, FrameInfo currentFI) { + + // ASSUMPTION : the bigger object retains the previous frame head direction. And the smaller frame head calculation starts + vector pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + // previous frame should have only one frame when this function is called + if (pFOVector.size() > 1) { + cerr<<"Previous frame does contain more than one fly object in the funciton separateObjectAndInitForNewHeadCalculation()"< currentFLargeCentroid = cFLargeFO.getCentroid(); + initLargeLocationX = currentFLargeCentroid.first; + initLargeLocationY = currentFLargeCentroid.second; + + cout << "Init location of the large object when get departed first ("< currentFSmallCentroid = cFSmallFO.getCentroid(); + initSmallLocationX = currentFSmallCentroid.first; + initSmallLocationY = currentFSmallCentroid.second; + + + cout << "Init location of the small object when get departed first ("< windowLargeHeadDirection; + + +void computeHeadBackward(FrameInfo nextFI, FrameInfo ¤tFI) { + vector nextFOVector = nextFI.getFOVector(); + vector currentFOVector = currentFI.getFOVector(); + + // next frame objects + FlyObject nextFLargeFO = nextFOVector[0]; + FlyObject nextFSmallFO = nextFOVector[1]; + + // current frame objects + FlyObject currentFLargeFO = currentFOVector[0]; + FlyObject currentFSmallFO = currentFOVector[1]; + + // larger + pair cLEV = currentFLargeFO.getMajorAxisEV(); + normalizeVector(cLEV); + // historical head + pair nLHH = nextFLargeFO.getHead(); + // next frame head direction +// pair nLHD; +// calculateHeadVector(nextFLargeFO, nLHD); + // take the minimum angle with the historical head + pair nCentroid = nextFLargeFO.getCentroid(); + cout << "Next Large was at "< cLREV(-cLEV.first, -cLEV.second); + if (largeDotProdNHWCEV >=0) { + cout << "Current eigen is in direction with the historical head"< newHead(newHeadFirst,newHeadSecond); + currentFLargeFO.setHead(newHead); + currentFLargeFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + //identifiedELEV = eLEV; + + } else { + cout << "Current eigen is in opposite direction of the historical head\n"; + // record this vector + double newHeadFirst = 0.1*cLREV.first+0.9*nLHH.first; + double newHeadSecond = 0.1*cLREV.second+0.9*nLHH.second; + pair newHead(newHeadFirst,newHeadSecond); + currentFLargeFO.setHead(newHead); + currentFLargeFO.setHeadIsInDirectionMAEV(false); + // set the identified EV + //identifiedELEV = eLEV; + + } + + + // small + + pair cSEV = currentFSmallFO.getMajorAxisEV(); + normalizeVector(cSEV); + // small historical head + pair nSHH = nextFSmallFO.getHead(); + // next frame head direction +// pair nSHD; +// calculateHeadVector(nextFSmallFO, nSHD); + // take the minimum angle with the historical head + pair nSCentroid = nextFSmallFO.getCentroid(); + cout << "Next small centroid at "< cSREV(-cSEV.first, -cSEV.second); + if (smallDotProdNHWCEV >=0) { + cout << "Current eigen is in direction with the historical head"< newHead(newHeadFirst,newHeadSecond); + currentFSmallFO.setHead(newHead); + currentFSmallFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + //identifiedELEV = eLEV; + + } else { + cout << "Current eigen is in opposite direction of the historical head\n"; + // record this vector + double newHeadFirst = 0.1*cSREV.first+0.9*nSHH.first; + double newHeadSecond = 0.1*cSREV.second+0.9*nSHH.second; + pair newHead(newHeadFirst,newHeadSecond); + currentFSmallFO.setHead(newHead); + currentFSmallFO.setHeadIsInDirectionMAEV(false); + // set the identified EV + //identifiedELEV = eLEV; + + } + + + // update the currentFrame + // update the fIForHeadVector update + vector updatedFOVector; + updatedFOVector.push_back(currentFLargeFO); + updatedFOVector.push_back(currentFSmallFO); + + currentFI.setFOVector(updatedFOVector); + + +} + +void calculateNewHeadDirection() { + + cout << "Start calculating new head direction for large,small fly objects"< firstFOVector = firstFI.getFOVector(); + vector endFOVector = endFI.getFOVector(); + + // init frame objects + FlyObject firstFLargeFO = firstFOVector[0]; + FlyObject firstFSmallFO = firstFOVector[1]; + + // end frame objects + FlyObject endFLargeFO = endFOVector[0]; + FlyObject endFSmallFO = endFOVector[1]; + + // large + pair firstFLargeFOCentroid = firstFLargeFO.getCentroid(); + pair endFLargeFOCentroid = endFLargeFO.getCentroid(); + double largeDirectionX = endFLargeFOCentroid.first - firstFLargeFOCentroid.first; + double largeDirectionY = endFLargeFOCentroid.second - firstFLargeFOCentroid.second; + pair largeDirection(largeDirectionX, largeDirectionY); + normalizeVector(largeDirection); + + + // find for the large object + cout << "For the larger end frame head finding \n"; + pair eLEV = endFLargeFO.getMajorAxisEV(); + // normalize the eigenvector + normalizeVector(eLEV); + pair identifiedELEV; + +// // the head from the previous frame +// pair pLH = pFLargeFO.getHead(); + + // find the whether the object is moving forward or backward + double largerDotProdMovDirAndLCBD = calculateDotProduct(largeDirection, largeCollisionBeforeDirection); + + // find the dot product with the large object movement direction with end frames major axis both direction + double largerDotProdMovDirAndELEV = calculateDotProduct(largeDirection, eLEV); + pair eLREV(-eLEV.first, -eLEV.second); + //double largerDotProdMovDirAndELEV = calculateDotProduct(pLH, clREV); + if (largerDotProdMovDirAndELEV >=0 and largerDotProdMovDirAndLCBD >= 0) { + cout << "Current eigen is in direction with the movement direction and moving forward \n"; + // record this vector for history + double newHeadFirst = 0.1*eLEV.first+0.9*largeDirection.first; + double newHeadSecond = 0.1*eLEV.second+0.9*largeDirection.second; + pair newHead(newHeadFirst,newHeadSecond); + endFLargeFO.setHead(newHead); + endFLargeFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + identifiedELEV = eLEV; + + } else if (largerDotProdMovDirAndELEV <0 and largerDotProdMovDirAndLCBD >=0){ + cout << "Current eigen is in opposite direction of the movement direction and moving forward\n"; + // record this vector + double newHeadFirst = 0.1*eLREV.first+0.9*largeDirection.first; + double newHeadSecond = 0.1*eLREV.second+0.9*largeDirection.second; + pair newHead(newHeadFirst,newHeadSecond); + endFLargeFO.setHead(newHead); + endFLargeFO.setHeadIsInDirectionMAEV(false); + // set the identified EV to reverse direction of the current EV + identifiedELEV = eLREV; + + } + + // moving backward + else if (largerDotProdMovDirAndELEV <0 and largerDotProdMovDirAndLCBD <0) { + cout << "Current eigen is in direction with the movement direction and moving backward"< newHead(newHeadFirst,newHeadSecond); + endFLargeFO.setHead(newHead); + endFLargeFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + identifiedELEV = eLEV; + + } else if (largerDotProdMovDirAndELEV >=0 and largerDotProdMovDirAndLCBD <0 ) { + cout << "Current eigen is in opposite direction of movement direction and moving backward\n"; + // record this vector + double newHeadFirst = 0.1*eLREV.first+0.9*largeDirection.first; + double newHeadSecond = 0.1*eLREV.second+0.9*largeDirection.second; + pair newHead(newHeadFirst,newHeadSecond); + endFLargeFO.setHead(newHead); + endFLargeFO.setHeadIsInDirectionMAEV(false); + // set the identified EV to reverse direction of the current EV + identifiedELEV = eLREV; + + } + + + + + // small + pair firstFSmallFOCentroid = firstFSmallFO.getCentroid(); + pair endFSmallFOCentroid = endFSmallFO.getCentroid(); + double smallDirectionX = endFSmallFOCentroid.first - firstFSmallFOCentroid.first; + double smallDirectionY = endFSmallFOCentroid.second - firstFSmallFOCentroid.second; + pair smallDirection(smallDirectionX, smallDirectionY); + cout << "Small Direction "< smallRevDirection(-smallDirection.first, -smallDirection.second); + cout << "Normalized small Reverse Direction "< eSEV = endFSmallFO.getMajorAxisEV(); + // normalize the eigenvector + normalizeVector(eSEV); + pair identifiedESEV; + + // find the whether the object is moving forward or backward + double smallerDotProdMovDirAndSCBD = calculateDotProduct(smallDirection, smallCollisionBeforeDirection); + + // find the dot product with the small object movement direction with end frames major axis both direction + double smallerDotProdMovDirAndESEV = calculateDotProduct(smallDirection, eSEV); + pair eSREV(-eSEV.first, -eSEV.second); + + if (smallerDotProdMovDirAndESEV >=0 && smallerDotProdMovDirAndSCBD >=0) { + cout << "Current eigen is in direction with the movement direction and moving forward "< newHead(newHeadFirst,newHeadSecond); + endFSmallFO.setHead(newHead); + endFSmallFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + identifiedESEV = eSEV; + + } else if (smallerDotProdMovDirAndESEV <0 and smallerDotProdMovDirAndSCBD >=0) { + cout << "Current eigen is in opposite direction of the movement direction and moving forward \n"; + // record this vector + double newHeadFirst = 0.1*eSREV.first+0.9*smallDirection.first; + double newHeadSecond = 0.1*eSREV.second+0.9*smallDirection.second; + pair newHead(newHeadFirst,newHeadSecond); + endFSmallFO.setHead(newHead); + endFSmallFO.setHeadIsInDirectionMAEV(false); + // set the identified EV + identifiedESEV = eSREV; + + } else if (smallerDotProdMovDirAndESEV <0 and smallerDotProdMovDirAndSCBD <0 ) { + cout << "Current eigen is in direction with the movement direction and moving backward"< newHead(newHeadFirst,newHeadSecond); + endFSmallFO.setHead(newHead); + endFSmallFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + identifiedESEV = eSEV; + + } else if (smallerDotProdMovDirAndESEV >=0 and smallerDotProdMovDirAndSCBD <0 ) { + cout << "Current eigen is in opposite direction with the movement direction and moving backward\n"; + // record this vector + double newHeadFirst = 0.1*eSREV.first+0.9*smallRevDirection.first; + double newHeadSecond = 0.1*eSREV.second+0.9*smallRevDirection.second; + pair newHead(newHeadFirst,newHeadSecond); + endFSmallFO.setHead(newHead); + endFSmallFO.setHeadIsInDirectionMAEV(false); + // set the identified EV + identifiedESEV = eSREV; + + } + + + + // update the fIForHeadVector update + vector updatedFOVector; + updatedFOVector.push_back(endFLargeFO); + updatedFOVector.push_back(endFSmallFO); + + endFI.setFOVector(updatedFOVector); + + // update the vector with the correct info + fIForHeadVector[endIndexToFindNewHead-startIndexToFindNewHead] = endFI; + + // calculate the head directions from end to first backwards + FrameInfo nextFI = fIForHeadVector[endIndexToFindNewHead-startIndexToFindNewHead]; + + cout << "fIForHeadVector size "<=0; i--) { + + cout << "Calculating backwards head direction i "<>fileName; + //inputFileName = "output/identified/"+fileName; + FrameInfo prevFI = fIVector[startIndex]; + cout << "Extracting information for image "<< prevFI.getFrameNo() << endl; + cout << "----------------------------------------\n"; + cout<>fileName; + fileName = fnVector[i]; + // inputFileName = "output/identified/"+fileName; + drawTheFlyObject(nextFI, fileName, isFirst, singleBlob); + } + + //inputFile.close(); + +} + +void normalizeVector(pair &a) { + double temp = a.first*a.first + a.second*a.second; + temp = sqrt(temp); + if (temp != 0) { + a.first = a.first/temp; + a.second = a.second/temp; + } +} +double calculateDotProduct(pair v, pair eV) { + return (v.first*eV.first + v.second*eV.second); +} +void calculateHeadVector(FlyObject fO, pair &headDirection) { + bool headDirectionBool = fO.getHeadIsInDirectionMAEV(); + pair fOMajorAxis = fO.getMajorAxisEV(); + if (headDirectionBool == true) { + headDirection.first = fOMajorAxis.first; + headDirection.second = fOMajorAxis.second; + } else { + headDirection.first = -fOMajorAxis.first; + headDirection.second = -fOMajorAxis.second; + } +} +void determineHeadDirection(int fileCounter) { + + FrameInfo prevFI = fIVector[fileCounter-1]; + FrameInfo currentFI = fIVector[fileCounter]; + + vector pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + FlyObject pFLargeFO = pFOVector[0]; + FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + // calculate velocity + pair pFLCentroid = pFLargeFO.getCentroid(); + pair pFSCentroid = pFSmallFO.getCentroid(); + + pair cFLCentroid = cFLargeFO.getCentroid(); + pair cFSCentroid = cFSmallFO.getCentroid(); + + double velocityXLarge = static_cast (cFLCentroid.first) - static_cast (pFLCentroid.first); + double velocityYLarge = static_cast (cFLCentroid.second) - static_cast (pFLCentroid.second); + + double velocityXSmall = static_cast (cFSCentroid.first) - static_cast (pFSCentroid.first); + double velocityYSmall = static_cast (cFSCentroid.second) - static_cast (pFSCentroid.second); + + pair cLVV(velocityXLarge,velocityYLarge); + pair cSVV(velocityXSmall,velocityYSmall); + + cFLargeFO.setVelocityV(cLVV); + cFSmallFO.setVelocityV(cSVV); + + // normalize the velocity + cFLargeFO.normalizeVelocity(); + cFSmallFO.normalizeVelocity(); + + // determine the head direction for larger object + pair cLEV = cFLargeFO.getMajorAxisEV(); + // normalize the eigenvector + normalizeVector(cLEV); + pair identifiedCLEV; + + // the head from the previous frame + pair pLH = pFLargeFO.getHead(); + + // find the dot product with the previous frame head with current major axis both direction + double largerDotProdPrevHeadAndCLEV = calculateDotProduct(pLH, cLEV); + pair cLREV(-cLEV.first, -cLEV.second); + //double largerDotProdPrevHeadAndCLREV = calculateDotProduct(pLH, clREV); + if (largerDotProdPrevHeadAndCLEV >=0) { + cout << "Current eigen is in direction with the historical head"< newHead(newHeadFirst,newHeadSecond); + cFLargeFO.setHead(newHead); + cFLargeFO.setHeadIsInDirectionMAEV(true); + // set the identified EV + identifiedCLEV = cLEV; + + } else { + cout << "Current eigen is in opposite direction of the historical head\n"; + // record this vector + double newHeadFirst = 0.1*cLREV.first + 0.9*pLH.first; + double newHeadSecond =0.1*cLREV.second+0.9*pLH.second; + pair newHead(newHeadFirst,newHeadSecond); + cFLargeFO.setHead(newHead); + cFLargeFO.setHeadIsInDirectionMAEV(false); + // set the identified EV to reverse direction of the current EV + identifiedCLEV = cLREV; + + } + + + // tracking the forward and backward movement of the fly object + double largeDotProd = calculateDotProduct(cLVV, identifiedCLEV); + cout << "largerDotProd "<= 0) { + cout<<"Larger Dot prod is positive"< cSEV = cFSmallFO.getMajorAxisEV(); + normalizeVector(cSEV); + pair identifiedCSEV; + // head from the previous frame + pair pSH = pFSmallFO.getHead(); + + // find the dot product with the previous frame head with current major axis both direction + double smallerDotProdPrevHeadAndCSEV = calculateDotProduct(pSH, cSEV); + pair cSREV(-cSEV.first, -cSEV.second); + if (smallerDotProdPrevHeadAndCSEV >=0) { + cout << "Current eigen is in direction with the historical head for the smaller fly object"< newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + cFSmallFO.setHeadIsInDirectionMAEV(true); + // set the identified EV to current EV + identifiedCSEV = cSEV; + } else { + cout << "Current eigen is in direction with the historical head for the smaller fly object"< newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + cFSmallFO.setHeadIsInDirectionMAEV(false); + // set the identified EV to reverse direction of current EV + identifiedCSEV = cSREV; + + } + + // to detect the swift change of head direction with 2 frames + // find previous frame's HD + pair previousHeadDirection; + calculateHeadVector(pFSmallFO, previousHeadDirection); + // find the current frame's HD + pair currentHeadDirection; + calculateHeadVector(cFSmallFO, currentHeadDirection); + + // calculate the dot product of previous head direction and current head direction. + double previousHDDotCurrentHD = calculateDotProduct(previousHeadDirection, currentHeadDirection); + if (previousHDDotCurrentHD < 0 and fileCounter > 1) { + + // now check if velocity direction( for the future historical velocity might be considered too) and + // current head direction dot product is also less than zero + // asssumption that current velocity and new head direction are towards same direction + // pair currentVelocity = cFSmallFO.getVelocityV(); + // double currentVVDotCurrentHD = calculateDotProduct(currentVelocity, currentHeadDirection); + // if (currentVVDotCurrentHD < 0) + { + // toggle current head direction + bool currentHeadDirectionBool = cFSmallFO.getHeadIsInDirectionMAEV(); + if (currentHeadDirectionBool == true) + cFSmallFO.setHeadIsInDirectionMAEV(false); + else { + cFSmallFO.setHeadIsInDirectionMAEV(true); + + } + + // update the historical head + // reset historical head to conform to the swift change in head direction + bool cFSmallFOFinalHeadDirectionBool = cFSmallFO.getHeadIsInDirectionMAEV(); + if (cFSmallFOFinalHeadDirectionBool == true) { + // record this vector for history + double newHeadFirst = cSEV.first; + double newHeadSecond = cSEV.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } else { + // record this vector for history + double newHeadFirst = cSREV.first; + double newHeadSecond = cSREV.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } + + + } + } + + /* retain historical head + bool cFSmallFOFinalHeadDirectionBool = cFSmallFO.getHeadIsInDirectionMAEV(); + if (cFSmallFOFinalHeadDirectionBool == true) { + // record this vector for history + double newHeadFirst = 0.1*cSEV.first + 0.9*pSH.first; + double newHeadSecond = 0.1*cSEV.second + 0.9*pSH.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } else { + // record this vector for history + double newHeadFirst = 0.1*cSREV.first + 0.9*pSH.first; + double newHeadSecond = 0.1*cSREV.second + 0.9*pSH.second; + pair newHead(newHeadFirst, newHeadSecond); + cFSmallFO.setHead(newHead); + + } + */ + + double smallDotProd = calculateDotProduct(cSVV, cSEV); + + if (smallDotProd >= 0) { + cout << "Smaller dot product is positive"< updatedFOVector; + updatedFOVector.push_back(cFLargeFO); + updatedFOVector.push_back(cFSmallFO); + + currentFI.setFOVector(updatedFOVector); + + fIVector[fileCounter] = currentFI; + + // cout << "checking the update"< tempFIFOVector = tempFI.getFOVector(); + // FlyObject tempFILFO = tempFIFOVector[0]; + // FlyObject tempFISFO = tempFIFOVector[1]; + // pair tempFILFOVV = tempFILFO.getVelocityV(); + // cout << "Large object velocity vector "< tempFISFOVV = tempFISFO.getVelocityV(); + // cout << "Small object velocity vector "< pFHH = prevFO.getHead(); + + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdEVAndHH = calculateDotProduct(cFOEV, pFHH); + if (dotProdEVAndHH > 0 ) { + cout << "Current head is in direction with the Previous frame Head(which is the historical head from the maxDistIndex)"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(true); + + } else { + cout << "Current head is in direction with the previous frame head"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(false); + + } + + // copy the velocity vector from the previous frame + pair pFVV = prevFO.getVelocityV(); + currentFO.setVelocityV(pFVV); + +} + +void objectHeadDirection(FlyObject prevFO, FlyObject ¤tFO, bool prevFFOHD) { + + // take the head direction from the previous frame + //pair pFHH = prevFO.getHead(); + + pair pFHeadDir; + calculateHeadVector(prevFO, pFHeadDir); + + + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdEVAndHD = calculateDotProduct(cFOEV, pFHeadDir); + if (dotProdEVAndHD > 0 ) { + cout << "Current head is in direction with the Previous frame Head direction"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(true); + */ + currentFO.setHead(cFOEV); + currentFO.setHeadIsInDirectionMAEV(true); + + } else { + cout << "Current head is in reverse direction with the previous frame head direction"< newHead(newHeadX, newHeadY); + currentFO.setHead(newHead); + currentFO.setHeadIsInDirectionMAEV(false); + */ + currentFO.setHead(cFOREV); + currentFO.setHeadIsInDirectionMAEV(false); + + } + + // copy the velocity vector from the previous frame + if (prevFFOHD == true) { + pair pFVV = prevFO.getVelocityV(); + currentFO.setVelocityV(pFVV); + } + +} + +void objectHeadDirection(FlyObject ¤tFO, int saveEV=0) { + + // get the velocity vector + pair cFOVV = currentFO.getVelocityV(); + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdVVAndEV = calculateDotProduct(cFOEV, cFOVV); + if (dotProdVVAndEV > 0 ) { + cout << "Current head is in direction with the velocity vector\n"; + // set the head to the eigen vector + currentFO.setHead(cFOEV); + currentFO.setHeadIsInDirectionMAEV(true); + + if (saveEV == 1) { + cout << "saving the eigen vector for the first object"< zero(0.0,0.0); + + if (saveEV == 1) { + cout << "saving the zero eigen vector for the first object"< cFV) { + +// // get the velocity vector +// pair cFOVV = currentFO.getVelocityV(); + pair cFOEV = currentFO.getMajorAxisEV(); + normalizeVector(cFOEV); + pair cFOREV(-cFOEV.first, -cFOEV.second); + + double dotProdVVAndEV = calculateDotProduct(cFOEV, cFV); + if (dotProdVVAndEV > 0 ) { + cout << "Current head is in direction with the vector used\n"; + // set the head to the eigen vector + currentFO.setHead(cFOEV); + currentFO.setHeadIsInDirectionMAEV(true); + } else { + cout << "Current head is in reverse direction with the vector used"< mFOVector = midFI.getFOVector(); + FlyObject mFFirstFO = mFOVector[0]; + FlyObject mFSecondFO = mFOVector[1]; + + pair midFCentroid = mFFirstFO.getCentroid(); + pair midSCentroid = mFSecondFO.getCentroid(); + + cout << "mid index is "< &velDirectionF, pair&velDirectionS) { + + // find the average velocity vector + cout << "Finding average velocity vector from "< pFOVector = prevFI.getFOVector(); + FlyObject pFFirstFO = pFOVector[0]; + FlyObject pFSecondFO = pFOVector[1]; + + + FrameInfo currentFI = fIVector[end]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + // get the summation of (velocity direction from start frame to half of the frames in this interval) + pair cFFCentroid = cFFirstFO.getCentroid(); + pair cFSCentroid = cFSecondFO.getCentroid(); + + + pair pFFCentroid = pFFirstFO.getCentroid(); + pair pFSCentroid = pFSecondFO.getCentroid(); + + + int velXFirst = cFFCentroid.first - pFFCentroid.first; + int velYFirst = cFFCentroid.second - pFFCentroid.second; + + int velXSecond = cFSCentroid.first - pFSCentroid.first; + int velYSecond = cFSCentroid.second - pFSCentroid.second; + + cout << "Velocity vector"< cFVV(static_cast (velXFirst), static_cast (velYFirst)); + pair cSVV(static_cast (velXSecond), static_cast (velYSecond)); + + velDirectionF = cFVV; + velDirectionS = cSVV; + + +} +double getSpeed(pair vector) { + double value = vector.first*vector.first + vector.second*vector.second; + value = sqrt(value); + + return value; + + +} + +void velocityDirections(int stIndex, int endIndex) { + + + velocityDirectionsF.clear(); + velocityDirectionsS.clear(); + + //int intervalLength = 5; + cout << "Initial velocity direction calculation"< velDirectionF; + pair velDirectionS; + + // get the velocity + velocityDirection(i,i+1, velDirectionF, velDirectionS); + + // extract the fly object for update + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + // get the speed of the velocity + double speedF = getSpeed(velDirectionF); + cFFirstFO.setSpeed(speedF); + double speedS = getSpeed(velDirectionS); + cFSecondFO.setSpeed(speedS); + + // save the velocity direction + normalizeVector(velDirectionF); + velocityDirectionsF.push_back(velDirectionF); + normalizeVector(velDirectionS); + velocityDirectionsS.push_back(velDirectionS); + + cout << "Normalized velocity"< updatedFOVector; + updatedFOVector.push_back(cFFirstFO); + updatedFOVector.push_back(cFSecondFO); + currentFI.setFOVector(updatedFOVector); + cout << "Updating the frame "< > velocityDirs, int &startIndex, int &endIndex) { + + int positiveVelSeqSize = 0; + int flag = false; + int maxSeqSize = 0; + int st = 0; + for (int j=0; j prevVel = velocityDirs[j]; + pair currVel = velocityDirs[j+1]; + + double dotProd = calculateDotProduct(prevVel, currVel); + + if( dotProd > 0 && flag == false) { + st = j; + positiveVelSeqSize++; + flag = true; + //cout << "In first if positiveSize "< 0 && flag == true) { + positiveVelSeqSize++; + //cout << "In second if positive "< maxSeqSize) { + maxSeqSize = positiveVelSeqSize; + startIndex = st; + endIndex = st+positiveVelSeqSize; + //cout << "maxseq updated \npositiveSize "< prevVel = velocityDirs[j]; + if (prevVel.first != 0 || prevVel.second != 0) { + startIndex = j; + endIndex = j; + zero = true; + break; + } + } + + if (zero != true) { + cout << "All directions zero"< cFOVector = currentFI.getFOVector(); + FlyObject cFFO = cFOVector[object-1]; + pair cFVV = cFFO.getVelocityV(); + //cout << "Velocity before normalization "< cFEV = cFFO.getMajorAxisEV(); + cout << "Eigen vector before normalization "< maxVelValue) { + maxVelValue = velValue; + maxVelValIndex = i; + } + + + } + + cout << "Maximum speed is chosen for for frame "< cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + cout << "Setting representative head direction of frame "< tempVelocity = cFSecondFO.getVelocityV(); + //cout << "Velocity was "< intervalLength) { + + cout << "Propagating in the longest positive frame interval starting the direction from the frame "< pFOVector = prevFI.getFOVector(); + FlyObject pFFirstFO = pFOVector[0]; + FlyObject pFSecondFO = pFOVector[1]; + + for (int i=t+1; i<=e; i++) { + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< tempVelocity = cFSecondFO.getVelocityV(); + //cout << "Velocity was "<=s; i--) { + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< pFOVector = prevFI.getFOVector(); + FlyObject pFFirstFO = pFOVector[0]; + FlyObject pFSecondFO = pFOVector[1]; + + if (e < origEnd) { + cout << "Propagating front from "< cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< origStart ) { + cout << "Propagating down wards from "<=origStart; i--) { + FrameInfo currentFI = fIVector[i]; + vector cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + if (object == 1) { + //cout << "First object extract"< z = velocityDirectionsF[a]; + pair zEV = evDirectionF[a]; + cout <<" First object velocity = "< w = velocityDirectionsS[a]; + pair wEV = evDirectionS[a]; + + cout << "Second object velocity = "< pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + // bool currentlySingleBlob = currentFI.getIsSingleBlob(); + + // if (currentlySingleBlob == false) { + + FlyObject pFLargeFO = pFOVector[0]; + FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + double distLTL = euclideanDist(pFLargeFO, cFLargeFO); + cout << "previousFirst to currentFirst "< aCentroid = a.getCentroid(); + pair bCentroid = b.getCentroid(); + + //cout << "Distance from "< (aCentroid.first)-static_cast (bCentroid.first)), 2.0); + double y2 = pow((static_cast (aCentroid.second)-static_cast (bCentroid.second)), 2.0); + double dist = sqrt((x2+y2)); + cout << dist< cFOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = cFOVector[0]; + FlyObject cFSecondFO = cFOVector[1]; + + double currDist = euclideanDist(cFFirstFO, cFSecondFO); + + if (currDist > maxDist) { + maxDist = currDist; + maxDistIndex = i; + cout << "New max distance" << maxDist << " at the frame number "<=startOfATrackSequence; i--) { + FrameInfo currentFI = fIVector[i]; + objectCorrespondence(prevFI, currentFI); + + // update the fIVector + fIVector[i] = currentFI; + prevFI = currentFI; + + } + cout << "Assigning object correspondance between maxDistIndex+1 "<<(maxDistIndex+1)<<" to endIndex "< fOVector = currentFI.getFOVector(); + FlyObject cFFirstFO = fOVector[0]; + FlyObject cFSecondFO = fOVector[1]; + + sequenceFirstAverage += cFFirstFO.getArea(); + sequenceSecondAverage += cFSecondFO.getArea(); + + cout << "For frame number "< sequenceSecondAverage) { + fout << "First object is the female"< sequenceSecondAverage) { + cout << "First sequence is the object A"< " << endl; // input file contains name of the + // input image files + return -1; + } + + MagickCore::SetMagickResourceLimit(MagickCore::MemoryResource, 1536); + MagickCore::SetMagickResourceLimit(MagickCore::MapResource, 4092); + + string fileName; + ifstream inputFile(argv[1]); + + // save the input file name + string ifns(argv[1]); + inputFileName = ifns; + + if (inputFile.fail() ) { + cout << "cannot open the input file that contains name of the input images\n"; + exit(1); + } + + // set the global paths + string tempOIP(argv[2]); + origImagePath = tempOIP; + + string tempFOP(argv[3]); + finalOutputPath = tempFOP; + + string tempMIP(argv[4]); + maskImagePath = tempMIP; + + unsigned int objCount = 0; + int i; + + int frameCounter = 0; + int fileCounter=0; + + char buffer[100]; + string imgSize; +// FlyObject a,b; + vector > shape; + vector tempFOV; + + // to find the new head direction + bool currentlyCalculatingHead = true; + + + + while (inputFile>>fileName) { + +// Image* img = new Image(argv[1]);// = new Image(argv[1]); + + int fi = fileName.find("Test"); + // current sequence numbers spans from 0 - 18019, so 5 digits are needed + int span = 5; + string tempString = fileName.substr(fi+4,span); + int frameCounter = atoi(tempString.c_str()); + cout << frameCounter<columns(),height = img->rows(); + diagLength= static_cast ( sqrt( (height*height) + (width*width) ) ); + cout << "Diagonal length is "< 0 ) + { + + cout << "size of the object is: " << s < eigenVal = covariantDecomposition(shape); + + { +// objCount++; + + double velocity_x=0.0, velocity_y=0.0; + // save the object information + FlyObject tempFO(s, + pair (eigenVal[6], eigenVal[7]), + pair (eigenVal[4], eigenVal[5]), + pair (velocity_x, velocity_y), + false, + pair (eigenVal[4], eigenVal[5]), + 0.0); + tempFOV.push_back(tempFO); + + + + + } + } + + } + } + +// cout<<"Sorting the objects according to size"< pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + bool currentlySingleBlob = currentFI.getIsSingleBlob(); + + if (currentlySingleBlob == false) { + + FlyObject pFLargeFO = pFOVector[0]; + //FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + // if just got rid off the single blob then the distance from the centroid of the large object before collision + // which is saved to the new two centroids needs to be calculated. This will give two direction + // dot product with the saved large object direction and directions are calculated the positive direction + // gives the new large object + /// ASSUMPTION is that the large object tends to go toward the direction before collision + // the smaller object will go either opposite or will be behind the object + // find the minimum distance among the centroids and chose the close one + // find minimum distance from the larger object = min1 + double distLTL = euclideanDist(pFLargeFO, cFLargeFO); + cout << "previousL to currentL "< pFOVector = prevFI.getFOVector(); + vector cFOVector = currentFI.getFOVector(); + + FlyObject pFLargeFO = pFOVector[0]; + FlyObject pFSmallFO = pFOVector[1]; + + FlyObject cFLargeFO = cFOVector[0]; + FlyObject cFSmallFO = cFOVector[1]; + + // find the area ratio of the current two objects. + // if they are within 5 percent, calculate the minimum distance from previous larger to current larger and smaller + double curr_larger = cFLargeFO.getArea(); + double curr_smaller = cFSmallFO.getArea(); + double ratio = (curr_smaller/curr_larger) * 100.0; + if (ratio >= 95.0) { + cout <<"Current ratio greater than equal 95.0 percent. Areas are : "<columns()<<" and height "<rows()< fOVector = currentFI.getFOVector(); + cout << "While drawing it found objects = "< centroid = currentFO.getCentroid(); + pair majorAxisEV = currentFO.getMajorAxisEV(); + + bool eVDirection = currentFO.getHeadIsInDirectionMAEV(); + + double ev_x, ev_y; + + if (eVDirection == true) { + + x0 = centroid.first; + y0 = centroid.second; + ev_x = static_cast(centroid.first) + static_cast(diagLength)* majorAxisEV.first; + ev_y = static_cast(centroid.second) + static_cast(diagLength)* majorAxisEV.second; + + x1 = static_cast (ev_x); + y1 = static_cast (ev_y); + + } else { + + x0 = centroid.first; + y0 = centroid.second; + ev_x = static_cast(centroid.first) - static_cast(diagLength)* majorAxisEV.first; + ev_y = static_cast(centroid.second) - static_cast(diagLength)* majorAxisEV.second; + + x1 = static_cast (ev_x); + y1 = static_cast (ev_y); + + + } + + int isHitting = draw_line_bm(maskImage, x0, y0, x1, y1 ); + + + + + + for (int n=0; n centroid = currentFO.getCentroid(); + pair majorAxisEV = currentFO.getMajorAxisEV(); + + bool eVDirection = currentFO.getHeadIsInDirectionMAEV(); + + double ev_x, ev_y; + //draw the object tracking circle + if (n == isFirst and n==0) { + cout << "Tracking the n = "<strokeColor("yellow"); + img->draw(DrawableCircle(centroid.first, centroid.second, centroid.first+5, centroid.second)); + + } else if ( n == isFirst and n==1) { + + cout << "Tracking the "<strokeColor("yellow"); + //img->fillColor("none"); + img->draw(DrawableCircle(centroid.first, centroid.second, centroid.first+5, centroid.second)); + + + } + + + // draw the female when tracked by the male fly + if (isHitting == 1) { + if (n != isFirst and n == 0) { + img->strokeColor("red"); + + /*img->draw(DrawableLine(centroid.first - 5, centroid.second - 5, centroid.first + 5, centroid.second - 5)); + img->draw(DrawableLine(centroid.first + 5, centroid.second - 5, centroid.first + 5, centroid.second + 5)); + img->draw(DrawableLine(centroid.first + 5, centroid.second + 5, centroid.first - 5, centroid.second + 5)); + img->draw(DrawableLine(centroid.first - 5, centroid.second + 5, centroid.first - 5, centroid.second - 5)); + */ + img->draw(DrawableRectangle(centroid.first - 5, centroid.second - 5, centroid.first + 5, centroid.second + 5)); + } else if (n != isFirst and n==1) { + + img->strokeColor("red"); + + /*img->draw(DrawableLine(centroid.first - 5, centroid.second - 5, centroid.first + 5, centroid.second - 5)); + img->draw(DrawableLine(centroid.first + 5, centroid.second - 5, centroid.first + 5, centroid.second + 5)); + img->draw(DrawableLine(centroid.first + 5, centroid.second + 5, centroid.first - 5, centroid.second + 5)); + img->draw(DrawableLine(centroid.first - 5, centroid.second + 5, centroid.first - 5, centroid.second - 5));*/ + + img->draw(DrawableRectangle(centroid.first - 5, centroid.second - 5, centroid.first + 5, centroid.second + 5)); + } + } + + // draw the Axis direction + if (eVDirection == true) { + ev_x = static_cast(centroid.first) + 50.0 * majorAxisEV.first; + ev_y = static_cast(centroid.second) + 50.0 * majorAxisEV.second; + img->strokeColor("green"); + img->draw( DrawableLine( centroid.first, centroid.second, static_cast(ev_x), static_cast(ev_y) )); + } else { + ev_x = static_cast(centroid.first) - 50.0 * majorAxisEV.first; + ev_y = static_cast(centroid.second) - 50.0 * majorAxisEV.second; + img->strokeColor("green"); + img->draw( DrawableLine( centroid.first, centroid.second, static_cast(ev_x), static_cast(ev_y) )); + } + + // draw the velocity vector + + /*img->strokeColor("blue"); + pair velocityV = currentFO.getVelocityV(); + ev_x = static_cast(centroid.first) + 30.0 * velocityV.first; + ev_y = static_cast(centroid.second) + 30.0 * velocityV.second; + img->draw(DrawableLine( centroid.first, centroid.second, static_cast(ev_x), static_cast(ev_y) )); + */ + + /* + // draw the historical head vector + img->strokeColor("white"); + pair headV = currentFO.getHead(); + ev_x = static_cast (centroid.first) + 25.0*headV.first; + ev_y = static_cast (centroid.second) + 25.0*headV.second; + img->draw( DrawableLine(centroid.first, centroid.second, static_cast (ev_x), static_cast (ev_y)) ); + */ + + /* + // overall velocity direction + if (n == 0) { + ev_x = static_cast (centroid.first) + 10*overAllVelocityF.first; + ev_y = static_cast (centroid.second)+ 10*overAllVelocityF.second; + img->strokeColor("cyan"); + img->draw(DrawableLine(centroid.first, centroid.second, static_cast (ev_x), static_cast (ev_y) ) ); + } else { + ev_x = static_cast (centroid.first) + 10*overAllVelocityS.first; + ev_y = static_cast (centroid.second)+ 10*overAllVelocityS.second; + img->strokeColor("cyan"); + img->draw(DrawableLine(centroid.first, centroid.second, static_cast (ev_x), static_cast (ev_y) ) ); + + } + + */ + // average velocity direction + /*if (n == 0) { + ev_x = static_cast (centroid.first) + 10*avgVelocityF.first; + ev_y = static_cast (centroid.second)+ 10*avgVelocityF.second; + img->strokeColor("cyan"); + img->draw(DrawableLine(centroid.first, centroid.second, static_cast (ev_x), static_cast (ev_y) ) ); + } else { + ev_x = static_cast (centroid.first) + 10*avgVelocityS.first; + ev_y = static_cast (centroid.second)+ 10*avgVelocityS.second; + img->strokeColor("cyan"); + img->draw(DrawableLine(centroid.first, centroid.second, static_cast (ev_x), static_cast (ev_y) ) ); + + }*/ + + + } + } + // overwrite the file now with axis + // img->write(inputFileName.c_str()); + + // when do not want to identify on the original comment below line and uncomment the above one + img->write(outputFileName.c_str()); + delete img; + delete maskImage; + +} + +void findObj(Image* img, int x, int y, vector > & 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 > & 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 p; + p.first = x; + p.second = y; + obj.push_back(p); + +// if (obj.size() > barrier) { +// //cout<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 > & 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 p; + p.first = x; + p.second = y; + obj.push_back(p); + +// if (obj.size() > barrier) { +// //cout<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 getCentroid(vector > & points) +{ + pair centroid; + centroid.first = 0; + centroid.second = 0; + + for (unsigned int i = 0; i covariantDecomposition(vector > & points) +{ + unsigned int i,j,k; + pair centroid = getCentroid(points); + vector retval; + + gsl_matrix* matrice = gsl_matrix_alloc(2, 2); + + double sumX2 = 0, sumXY = 0, sumY2 = 0; + for (k = 0; ksize; i++) + retval.push_back(gsl_vector_get(eigenVal, i)); + + for (j = 0; jsize2; j++) + for (i = 0; isize1; i++) + retval.push_back(gsl_matrix_get(eigenVec, i, j)); + + retval.push_back(static_cast(centroid.first)); + retval.push_back(static_cast (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; + +} diff --git a/fly-tools/MyPair.h b/fly-tools/MyPair.h new file mode 100644 index 0000000..6a2a4fa --- /dev/null +++ b/fly-tools/MyPair.h @@ -0,0 +1,30 @@ +class MyPair { + +public: + int first, second; + MyPair(int x, int y) { + + this->first = x; + this->second = y; + + } + MyPair() { + + this->first = -1; + this->second = -1; + + } + MyPair(const MyPair &f) { + this->first = f.getX(); + this->second = f.getY(); + + } + int getX() const { + return this->first; + } + + int getY() const { + return this->second; + } + +}; diff --git a/fly-tools/StandardDeviation.cpp b/fly-tools/StandardDeviation.cpp new file mode 100644 index 0000000..a157d11 --- /dev/null +++ b/fly-tools/StandardDeviation.cpp @@ -0,0 +1,90 @@ +#include +#include +#include +#include +#include +using namespace std; + +vector currentHistogramValues; + +int main(int argc, char* argv[]) { + + // argv[0] name of the executable + // argv[1] input file containing the list of files + // argv[2] output file containing the standard deviation and file name pair + // argv[3] data file path + // argv[4] data file postfix of metricname + + if (argc < 3) { + cout<<"Please provide the parameters ./executable iputfilename outputfilename"<>currentFileName) { + + string currentFileWithExtension = prefixPath + currentFileName + "/" + currentFileName + "_"+ metricName +".txt"; + ifstream currentFile(currentFileWithExtension.c_str()); + if (currentFile.fail() == true) { + cout << currentFileName + " cannot be opened"<> currentValueOfHistogram) { + sumOfValues = sumOfValues + i*currentValueOfHistogram; + currentHistogramValues.push_back(currentValueOfHistogram); + + N = N + currentValueOfHistogram; + i=i+1.0; + } + + M = currentHistogramValues.size(); + // mean + double mean = sumOfValues/N; + + + // sigma^2 = (sum( (i-mean)^2*H_i ) )/(N-1) for i = 0 to M-1 + double standardDev = 0.0; + double sumSquaredResults = 0.0; + int j = 0; + for (i=0.0; i +#include +#include +#include +#include +#include +#include +#include "jzImage.h" + +using namespace Magick; +using namespace std; + +unsigned int jzImage::_h; +unsigned int jzImage::_w; +char* jzImage::data = NULL; + +// colors in debugging +const int C = 10; + +static map len; +bool inWhite; +bool startedWhite; +int x_init, y_init; + +inline int roundT(double v) { return int(v+0.5); } +inline int dist(int x0, int y0, int x1, int y1) {return roundT(sqrt(pow(double(x1)-x0, 2.0)+pow(double(y1)-y0, 2.0)));} + +inline void drawLine(int x0, int y0, int x1, int y1, jzImage& img); +inline void lookAt(int x, int y, jzImage& img); + + +void writeHist(const char* filename) +{ + map::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 hist(last-first, 0); + vector hist(last+1, 0); + + cout << "hist size: " << hist.size() << endl; + try{ + for(unsigned int j = 0; j= 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 len) +{ + map::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 hist(last-first, 0); + vector hist(last+1, 0); + + cout << "hist size: " << hist.size() << endl; + try{ + for(unsigned int j = 0; j= 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 " << endl; + return -1; + } + + int i, j; + + inWhite = false; + len.clear(); + + MagickCore::SetMagickResourceLimit(MagickCore::MemoryResource, 1536); + MagickCore::SetMagickResourceLimit(MagickCore::MapResource, 4096); + Image* img = new Image(argv[1]); + + int width = img->columns(), + height = img->rows(); + + jzImage &myImg = jzImage::Instance(img); + + delete img; + + /* + // doesn't take pixel(0,0) and also takes pixel(width-1, height-1) twice + //map, pair >, bool> hitMap; + const int boundSize = 2*(width+height-2); + static pair* boundary = new pair[boundSize]; + + for (i = 0; i, pair >, bool> hitMap; + const int boundSize = 2*(width+height-2); + static pair* boundary = new pair[boundSize]; + + for (i = 0; i<=width-1; i++) + { + boundary[i].first = i; + boundary[i].second = 0; + boundary[width+(height-2)+i].first = i; + boundary[width+(height-2)+i].second = height-1; + } + + for (i = 0; i= 1 + //drawLine(3,2,5,5,myImg); + // 0<=m<1 + //drawLine(5,5,1,4,myImg); + // -1< m < 0 + //drawLine(7,0,2,3,myImg); + // m<-1 + //drawLine(1,7,3,2,myImg); + // dx = 0 + //drawLine(5,0,5,5,myImg); + // dy = 0 + //drawLine(0,3,7,3,myImg); + // multiple length + //drawLine(1,7,7,1,myImg); + + + + for (i = 0; i lineSweepValLog; + map::iterator front = len.begin(), back = len.end(); + + for( front; front != back; front++ ) { + double logvalue = log(front->second); + lineSweepValLog[front->first] = logvalue; + } + + string fileNameLog(argv[2]); + fileNameLog = fileNameLog + "_log_hist.txt"; + writeHistLog(fileNameLog.c_str(), lineSweepValLog); + + + delete[] boundary; + + return 0; +} + +/* +void lookAt(int x, int y, jzImage& img) +{ + if (img.pixelColor(x,y) == 1) + // if it was the first white pixel then + if (!inWhite) + { + inWhite = true; + x_init = x; + y_init = y; + } + + if (img.pixelColor(x,y) == 0) + if (inWhite) + { + int val = roundT(dist(x_init, y_init, x, y)); + len[val]++; + inWhite = false; + } +} +*/ + +void lookAt(int x, int y, jzImage& img) +{ + int imageWidth =img.width(); + int imageHeight = img.height(); + // if current pixel is white + if (img.pixelColor(x,y) == 1) { + // if it was the first white pixel then + if (!inWhite) + { + // check whether it started as a white pixel; if it started as white pixel then this line segment should + // not be counted. because it is part of the larger white blob + if ( (x == 0 || x == (imageWidth -1) ) || ( y == 0 || y == (imageHeight - 1))) + startedWhite = true; + inWhite = true; + x_init = x; + y_init = y; + } + // if we are on a white region + /*else { + + } + */ + } + + if (img.pixelColor(x,y) == 0) { + // if we are going through a white line and reached the black pixel + if (inWhite) + { + // if the line started as a white pixel then discard the line fragment + if (startedWhite == false) { + int val = roundT(dist(x_init, y_init, x, y)); + //cout<<" val is = " << val << " at (x0,y0) = " << x_init << "," << y_init << " to (x1,y1) = " << x << "," << y < x1 + if (x0 > x1) + { + int temp = x0; + x0 = x1; + x1 = temp; + temp = y0; + y0 = y1; + y1 = temp; + } + + int dx, dy, d, x, y, incrE, incrNE; + dx = x1 - x0; + dy = y1 - y0; + y = y0; + x = x0; + + // if they are on a vertical line + if (dx == 0) + { + if (y0 < y1) + for (int i = y0; i<=y1; i++) + lookAt(x,i,img); + else + for (int i = y1; i<=y0; i++) + lookAt(x,i,img); + return; + } + + // if they are on a horizontal line + if (dy == 0) + { + for (int i = x0; i<=x1; i++) + lookAt(i,y,img); + return; + } + + int dir = 0; + double m = double(dy)/double(dx); + if (m >= 1.0) + dir = 1; + else if ( (m < 0.0) && (m > -1.0) ) + dir = 2; + else if (m <= -1.0) + dir = 3; + + switch(dir) + { + // when slope m: 0< m <1 + case 0: + d = dy*2 - dx; + incrE = dy*2; + incrNE = (dy - dx)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (x= 1 + case 1: + d = dx*2 - dy; + incrE = dx*2; + incrNE = (dx - dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y= 0) + { + d+= incrE; + x++; + } + else + { + d+= incrNE; + x++; + y--; + } + + lookAt(x,y,img); + } + break; + // since initially we swapped the P0 with P1 so that P0 always holds values with smaller x cooridinate + // so we are sure that m<0 is the result of y0 being greater that + case 3: + d = dx*2 + dy; + incrE = dx*2; + incrNE = (dx + dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y>y1) + { + if (d <= 0) + { + d+= incrE; + y--; + } + else + { + d+= incrNE; + x++; + y--; + } + lookAt(x,y,img); + } + break; + } + + + if (inWhite) + { + // it is a fraction of the blob so dont count it + //len[dist(x_init,y_init,x,y)+1]++; + inWhite = false; + startedWhite = false; + } +} +*/ + +void drawLine(int x0, int y0, int x1, int y1, jzImage& img) +{ + inWhite = false; + startedWhite = false; + + // always go from x0 -> x1 + if (x0 > x1) + { + int temp = x0; + x0 = x1; + x1 = temp; + temp = y0; + y0 = y1; + y1 = temp; + } + + int dx, dy, d, x, y, incrE, incrNE; + dx = x1 - x0; + dy = y1 - y0; + y = y0; + x = x0; + + // if they are on a vertical line + if (dx == 0) + { + if (y0 < y1) + for (int i = y0; i<=y1; i++) + lookAt(x,i,img); + else + for (int i = y1; i<=y0; i++) + lookAt(x,i,img); + return; + } + + // if they are on a horizontal line + if (dy == 0) + { + for (int i = x0; i<=x1; i++) + lookAt(i,y,img); + return; + } + + int dir = 0; + double m = double(dy)/double(dx); + if (m >= 1.0) + dir = 1; + else if ( (m < 0.0) && (m > -1.0) ) + dir = 2; + else if (m <= -1.0) + dir = 3; + + switch(dir) + { + // when slope m: 0< m <1 + case 0: + d = dy*2 - dx; + incrE = dy*2; + incrNE = (dy - dx)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (x= 1 + case 1: + d = dx*2 - dy; + incrE = dx*2; + incrNE = (dx - dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y= 0) + { + d+= incrE; + x++; + } + else + { + d+= incrNE; + x++; + y--; + } + + lookAt(x,y,img); + } + break; + // since initially we swapped the P0 with P1 so that P0 always holds values with smaller x cooridinate + // so we are sure that m<0 is the result of y0 being greater that + case 3: + d = dx*2 + dy; + incrE = dx*2; + incrNE = (dx + dy)*2; + x = x0; + y = y0; + + lookAt(x,y,img); + + while (y>y1) + { + if (d <= 0) + { + d+= incrE; + y--; + } + else + { + d+= incrNE; + x++; + y--; + } + lookAt(x,y,img); + } + break; + } + + + if (inWhite) + { + // it is a fraction of the blob so dont count it + //len[dist(x_init,y_init,x,y)+1]++; + inWhite = false; + startedWhite = false; + } +} diff --git a/fly-tools/misc/makeimage.cpp b/fly-tools/misc/makeimage.cpp new file mode 100644 index 0000000..06e2897 --- /dev/null +++ b/fly-tools/misc/makeimage.cpp @@ -0,0 +1,51 @@ +#include +#include +#include +#include +using namespace std; +using namespace Magick; + +int main(int argc,char **argv) + +{ + // Construct the image object. Seperating image construction from the + // the read operation ensures that a failure to read the image file + // doesn't render the image object useless. + + Image* img; + char buffer[100]; + sprintf(buffer,"%ix%i",7,7); + + // residual image is initialized with black representing not visited. + //residual = new Image(buffer, "black"); + + img = new Image(buffer, "white"); + + for (int j=0; j<=3;j++) { + for (int i=0; i<=3-j; i++) { + img->pixelColor(i,j, "black"); + img->pixelColor(6-i,j, "black"); + } + } + + int k; + for (int j=4; j<=6;j++) { + for (int i=0; i<=j-3; i++) { + img->pixelColor(i,j, "black"); + img->pixelColor(6-i,j, "black"); + } + } + + + //img->pixelColor(0,3, "red"); + + + string namei = "7x7.png"; + img->write(namei.c_str()); + + delete img; + + return 0; + +} + diff --git a/fly-tools/misc/test.cpp b/fly-tools/misc/test.cpp new file mode 100644 index 0000000..a709482 --- /dev/null +++ b/fly-tools/misc/test.cpp @@ -0,0 +1,92 @@ +#include +#include + +using namespace std; + +void largestIncreasingPositiveDotProductSeq(vector velocityDirs, int &startIndex, int &endIndex) { + + int positiveVelSeqSize = 0; + int flag = false; + int maxSeqSize = 0; + int st = 0; + for (int j=0; j 0 && flag == false) { + st = j; + positiveVelSeqSize++; + flag = true; + cout << "In first if positiveSize "< 0 && flag == true) { + positiveVelSeqSize++; + cout << "In second if positive "< maxSeqSize) { + maxSeqSize = positiveVelSeqSize; + startIndex = st; + endIndex = st+positiveVelSeqSize; + cout << "maxseq updated \npositiveSize "< test(4); + test[0]=0; + test[1]=7; + test[2]=0; + test[3]=7; + /*test[4]=5; + test[5]=2; + test[6]=-3; + test[7]=-4; + test[8]=2; + test[9]=3; + */ + int st; + int endIndex; + + largestIncreasingPositiveDotProductSeq(test,st, endIndex ); + for (int i=st; i<=endIndex; i++) { + cout << test[i]<=st; i=i-5) { + c = i; + cout << c << endl; + size++; + } + if ((c-st) != 0) { + cout <<"additional " <