From 68eed608949b3b2ecb3d87a72886ae917ff7bc0f Mon Sep 17 00:00:00 2001 From: mutantturkey Date: Thu, 14 Jun 2012 23:36:25 -0400 Subject: initial makefile for fly-tools, still broken --- fly-tools/misc/FurthestPointAlongMajorAxis.cpp | 3084 ++++++++++++++++++++++++ 1 file changed, 3084 insertions(+) create mode 100644 fly-tools/misc/FurthestPointAlongMajorAxis.cpp (limited to 'fly-tools/misc') diff --git a/fly-tools/misc/FurthestPointAlongMajorAxis.cpp b/fly-tools/misc/FurthestPointAlongMajorAxis.cpp new file mode 100644 index 0000000..8c8a086 --- /dev/null +++ b/fly-tools/misc/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; + +} -- cgit v1.2.3