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/FurthestPointAlongMajorAxis.cpp | 3084 ------------------------ fly-tools/Makefile | 51 + fly-tools/README | 0 fly-tools/config.mk | 23 + fly-tools/misc/FurthestPointAlongMajorAxis.cpp | 3084 ++++++++++++++++++++++++ 5 files changed, 3158 insertions(+), 3084 deletions(-) delete mode 100644 fly-tools/FurthestPointAlongMajorAxis.cpp create mode 100644 fly-tools/Makefile create mode 100644 fly-tools/README create mode 100644 fly-tools/config.mk create mode 100644 fly-tools/misc/FurthestPointAlongMajorAxis.cpp diff --git a/fly-tools/FurthestPointAlongMajorAxis.cpp b/fly-tools/FurthestPointAlongMajorAxis.cpp deleted file mode 100644 index 8c8a086..0000000 --- a/fly-tools/FurthestPointAlongMajorAxis.cpp +++ /dev/null @@ -1,3084 +0,0 @@ -#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/Makefile b/fly-tools/Makefile new file mode 100644 index 0000000..94b0636 --- /dev/null +++ b/fly-tools/Makefile @@ -0,0 +1,51 @@ +include config.mk + +SRC = FlyObject.cpp FlyTrackingFilter.cpp FrameInfo.cpp FlyTrackingMain.cpp +OBJ = ${SRC:.c=.o} + +all: options FilterFlyMask FlyTracking + +options: + @echo fly-tools build options: + @echo "CFLAGS = ${CFLAGS}" + @echo "LDFLAGS = ${LDFLAGS}" + @echo "CC = ${CC}" + @echo "SRC = ${SRC}" +.c.o: + @echo CC $< + @${CC} -c ${CFLAGS} $< + +${OBJ}: config.mk + +svte: ${OBJ} + @echo CC -o $@ + @${CC} -o $@ svte.o ${LDFLAGS} + +FlyTracking: ${OBJ} + @echo CC -o $@ + @${CC} -o $@ FlyTracking.o ${LDFLAGS} + +clean: + @echo cleaning + @rm -f FlyTracking FilterFlyMask ${OBJ} fly-${VERSION}.tar.gz + +dist: clean + @echo creating dist tarball + @mkdir -p fly-${VERSION} + @cp -R Makefile config.mk ${SRC} svte-${VERSION} + @tar -cf svte-${VERSION}.tar mt-${VERSION} + @gzip svte-${VERSION}.tar + @rm -rf svte-${VERSION} + +install: all + @echo installing executable file to ${DESTDIR}${PREFIX}/bin + @mkdir -p ${DESTDIR}${PREFIX}/bin + @cp -f svte ${DESTDIR}${PREFIX}/bin + @chmod 755 ${DESTDIR}${PREFIX}/bin/svte + +uninstall: + @echo removing executable file from ${DESTDIR}${PREFIX}/bin + @rm -f ${DESTDIR}${PREFIX}/bin/svte + + +.PHONY: all options clean dist install uninstall diff --git a/fly-tools/README b/fly-tools/README new file mode 100644 index 0000000..e69de29 diff --git a/fly-tools/config.mk b/fly-tools/config.mk new file mode 100644 index 0000000..e41bf88 --- /dev/null +++ b/fly-tools/config.mk @@ -0,0 +1,23 @@ +VERSION = 0.0.1 + +# Customize below to fit your system + +# paths +PREFIX = /usr +MANPREFIX = ${PREFIX}/share/man + + +# includes and libs + +GTKINC=$(shell pkg-config --cflags gtk+-2.0 vte ) +GTKLIB=$(shell pkg-config --libs gtk+-2.0 vte ) + +INCS = -I. -I/usr/include $(pkg-config --cflags Magick++ gsl) +LIBS = -L/usr/lib -lc $(pkg-config --libs Magick++ gsl) +# flags +CPPFLAGS = -DVERSION=\"${VERSION}\" +CFLAGS = -mtune=native -std=gnu99 -O2 -s ${INCS} ${CPPFLAGS} +LDFLAGS = -s ${LIBS} + +# compiler and linker +CC = g++ 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