IMOD - code - icontextra
Jump to navigation
Jump to search
About
NOTE: This page is a daughter page of: IMOD
I wrote a bunch of plugins for IMOD (a tomography and image segmentation program) back in my PhD. Some of the code files I wrote actually solved a lot of polygon problems, so I thought it would be useful to paste my old code here.
Disclaimer: This is GNU code, but, as with most academic code, it's a mess and many of the functions may not work 100%. Also, the polygon (contour/icont) classes and calls I make are specific to IMOD, but may still be helpful if you can translate to your own structures.
icontextra.h
#ifndef _ICONTEXTRA_H
#define _ICONTEXTRA_H
//############################################################
#include "_common_functions.h"
#include "imodplugin.h"
#include <vector>
using namespace std;
//############################################################
//-------------------------------
//## CONSTANTS
enum textalign { TA_LEFT, TA_CENTER, TA_RIGHT };
enum textalignv { TV_TOP, TV_CENTER, TV_BOTTOM, TV_LEADING };
const float DEFAULT_FONT_HEIGHT = 12.0f;
const float DEFAULT_FONT_SPACING = 10.0f;
//############################################################
//-------------------------------
//## SPECIAL DATA STRUCTURES USED IN FUCTIONS:
struct IcontPtr // used as wrapper to store a pointer to a contour
{ // so I can have vectors of "IcontPtr"
Icont *cont; // pointer to a contour
float area; // used in interoplator methods to stre and sort vector by area
IcontPtr( Icont* cont_ );
IcontPtr();
//~IcontPtr();
void deleteContour();
};
struct PtConnection // used especially for getIntersectingPolygon function
{
Ipoint intercept; // stores intercept point
int cont1Idx; // stores index for cont1 where intercept point exists
int cont2Idx; // stores index for cont1 where intercept point exists
bool included; // wether pts included yet in generating intersecting polygon(s)
PtConnection( Ipoint _intercept );
PtConnection( Ipoint _intercept, int _cont1Idx );
};
bool operator< (const PtConnection &a, const PtConnection &b);
bool operator== (const PtConnection lhs, const PtConnection rhs);
struct IdxAndFloat // used especially for getIntersectingPolygon function
{
int idx; // stores an idx reference to the intercept point (in conn)
float dist; // stores distance from a pt on a contour to its interecpt pt
IdxAndFloat(int _idx, float _dist);
};
bool operator<(const IdxAndFloat &a, const IdxAndFloat &b);
struct IntAndInt // used to connect point indexes between two contours
{
int idx1; // idx reference to a point in contour 1
int idx2; // idx reference to a point in contour 2
IntAndInt(int _idx1, int _idx2);
};
bool operator< (const IntAndInt &a, const IntAndInt &b);
bool operator== (const IntAndInt lhs, const IntAndInt rhs);
struct IdxToSort // used especially for getIntersectingPolygon function
{
int idx; // stores and idx reference to the intercept point (in conn)
float float1; // stores distance from a pt on a contour to its interecpt pt
float float2; // used as a tie breaker
IdxToSort();
IdxToSort( int _idx, float _float1 );
IdxToSort( int _idx, float _float1, float _float2 );
};
bool operator<(const IdxToSort &a, const IdxToSort &b);
//############################################################
//-------------------------------
//## INLINE FUNCTION DECLARATIONS: (DEFINED AT THE END OF THIS FILE)
inline bool isContValid(Icont *cont);
inline bool isEmpty(Icont *cont);
inline bool isInterpolated(Icont *cont);
inline void setInterpolated(Icont *cont, int on);
inline bool isOpenFlag(Icont *cont);
inline void setOpenFlag(Icont *cont, int on);
inline bool isDeleteFlag(Icont *cont);
inline void setDeleteFlag(Icont *cont, int on);
inline int psize(Icont *cont);
inline Ipoint *getPt(Icont *cont, int idx);
inline int getPtZInt(Icont *cont, int idx);
inline Ipoint *getPtNoWrap(Icont *cont, int idx);
inline Ipoint *getLastPt(Icont *cont );
inline Ipoint *getFirstPt(Icont *cont );
inline void setPt(Ipoint *pt, float x, float y, float z);
inline void copyPt(Ipoint *to, Ipoint *from);
inline Ipoint newPt( float x, float y, float z );
inline bool ptsEqual( Ipoint *pt1, Ipoint *pt2 );
inline bool ptsApproxEqual( Ipoint *pt1, Ipoint *pt2, float prec );
inline void removePtsSize( Icont *cont );
inline void removePtSize( Icont *cont, int idx );
inline bool isDefaultSize( Iobj *obj, Icont *cont, int idx );
inline void printPt( Ipoint *pt );
inline void printCont( Icont *cont );
inline void printContPts( Icont *cont );
inline void deleteAllPts( Icont *cont );
inline float getZ( Icont *cont );
inline int getZInt( Icont *cont );
inline float getZRange( Icont *cont );
inline void setZValue( Icont *cont, int newZValue );
inline void cont_copyPts( Icont *from, Icont *to, bool clearToCont );
inline void deleteContours( vector<IcontPtr> &conts );
inline void eraseContour( vector<IcontPtr> &conts, int idx );
inline Icont* getCont( Iobj *obj, int idx );
inline Icont* getLastCont( Iobj *obj );
inline int csize( Iobj *obj );
inline int osize( Imod *imod );
inline Iobj* getObj( Imod *imod, int idx );
inline bool isObjClosed(Iobj *obj);
inline bool isContClosed(Iobj *obj, Icont *cont);
inline bool isObjectValidAndShown(Iobj *obj);
inline void setObjColor( Iobj *obj, int red, int green, int blue );
//-------------------------------
//## POINT RELATED FUNCTIONS:
void point_rotatePointAroundPoint2D( Ipoint *pt, Ipoint *center, float theta );
void point_scalePtAboutPt( Ipoint *pt, Ipoint *center, float scaleX, float scaleY, float scaleZ );
void point_scalePtAboutPt2D( Ipoint *pt, Ipoint *center, float scaleX, float scaleY );
float point_findAvgDistToPt( Icont *pts, Ipoint *startPt, Ipoint *scalePt );
float point_maxDiffInDistToPoint( Icont *pts, Ipoint *startPt, Ipoint *scalePt );
float point_stdDevDistToPoint( Icont *pts, Ipoint *startPt, Ipoint *scalePt );
Ipoint point_findCenterOfPts( Iobj *obj, Icont *pts, Ipoint *startPt, float &avgRadius,
float zScale, float &avgResidual, int maxIts );
Ipoint point_findCenterOfPtsX( Iobj *obj, Icont *pts, Ipoint *startPt, float avgRadius, float zScale, int maxIts );
Ipoint point_findOptimalZScaleAndCenterOfPts( Iobj *obj, Icont *pts, Ipoint *startPt, float &avgRadius, float zScale,
float &bestZScale, float &bestAvgResidual, float startChangeZ,
float accuracyZScale, int maxIts );
float getValCardinalSpline( float fract, float p0, float p1, float p2, float p3, float tensileFract );
Ipoint getPtCardinalSpline( float fract, Ipoint p0, Ipoint p1, Ipoint p2, Ipoint p3, float tensileFract );
//-------------------------------
//## MINIMUM BOUNDING RECTANGLE (MBR) RELATED FUNCTIONS:
void mbr_reset( Ipoint *ll, Ipoint *ur );
void mbr_addPt( Ipoint *pt, Ipoint *ll, Ipoint *ur );
float mbr_distToNearestEdge(float val, float min, float max);
float mbr_distToNearestEdge(float min1, float max1, float min2, float max2);
float mbr_distBetweenBBoxes2D(Ipoint *ll1, Ipoint *ur1, Ipoint *ll2, Ipoint *ur2);
float mbr_distBetweenBBoxes3D(Ipoint *ll1, Ipoint *ur1, Ipoint *ll2, Ipoint *ur2, float zScale);
bool mbr_distToBBox2D(Ipoint *pt, Ipoint *ll, Ipoint *ur);
bool mbr_isPtInsideBBox(Ipoint *pt, Ipoint *ll, Ipoint *ur);
bool mbr_isPtInsideBBox2D(Ipoint *pt, Ipoint *ll, Ipoint *ur);
bool mbr_doEdgesOverlap(float min1, float max1, float min2, float max2);
bool mbr_doBBoxesOverlap2D(Ipoint *p1ll, Ipoint *p1ur, Ipoint *p2ll, Ipoint *p2ur);
bool mbr_isBBoxInsideBBox(Ipoint *ll1, Ipoint *ur1, Ipoint *ll2, Ipoint *ur2);
//-------------------------------
//## LINE RELATED FUNCTIONS:
float line_getAngle2D ( Ipoint *linept1, Ipoint *linept2 );
float line_getAngle2DPos ( Ipoint *pt1, Ipoint *pt2 );
float line_getRadians2D ( Ipoint *linept1, Ipoint *linept2 );
float line_radiansFormed3Pts( Ipoint *pt1, Ipoint *pt2, Ipoint *pt3 );
Ipoint line_getPtHalfwayBetween(Ipoint *pt1, Ipoint *pt2);
Ipoint line_findPtFractBetweenPts2D( const Ipoint *pt1, const Ipoint *pt2, float fractBetweenPts );
Ipoint line_findPtFractBetweenPts( const Ipoint *pt1, const Ipoint *pt2, const float fractBetweenPts );
float line_sqDistBetweenPts2D( Ipoint *pt1, Ipoint *pt2 );
float line_distBetweenPts2D( Ipoint *pt1, Ipoint *pt2 );
bool line_isPointOnLine( Ipoint *pt, Ipoint *lineStart, Ipoint *lineEnd );
bool line_getLineEquation( Ipoint *pt1, Ipoint *pt2, float *gradient, float *offset );
float line_crossProduct3Points( Ipoint *pt1, Ipoint *pt2, Ipoint *pt3);
float line_angleFormed3Pts( Ipoint *pt1, Ipoint *pt2, Ipoint *pt3 );
Ipoint line_getPtRelativeToEnd( Ipoint *start, Ipoint *end, float distFromEnd, float angleFromStraight );
bool line_doLinesCrossAndWhere( Ipoint *line1pt1, Ipoint *line1pt2, Ipoint *line2pt1, Ipoint *line2pt2, Ipoint *intercept );
bool line_getInterceptWhereRayCross( Ipoint *line1pt1, Ipoint *line1pt2, Ipoint *line2pt1, Ipoint *line2pt2, Ipoint *intercept );
double line_getFractPtBetweenTwoRays( Ipoint *ray1pt1, Ipoint *ray1pt2, Ipoint *ray2pt1, Ipoint *ray2pt2, Ipoint *pt );
double line_getFractPtBetweenTwoLines( Ipoint *l1p1, Ipoint *l1p2, Ipoint *l2p1, Ipoint *l2p2, Ipoint *pt, bool checkOutsideNegative );
Ipoint point_findPtInQuad1InMatchingQuad2( Ipoint *pt, Ipoint *q1BL, Ipoint *q1TL, Ipoint *q1TR, Ipoint *q1BR, Ipoint *q2BL, Ipoint *q2TL, Ipoint *q2TR, Ipoint *q2BR );
bool line_isKiss( Ipoint *pMid, Ipoint *p1,Ipoint *p2, Ipoint *p3, Ipoint *p4 );
bool line_twoKiss( Ipoint *a1, Ipoint *a2, Ipoint *a3, Ipoint *b1, Ipoint *b2, Ipoint *b3 );
//-------------------------------
//## EXTRA CONTOUR FUNCTIONS:
int cont_generateDigitUsing16SegDisplay( Iobj *obj, char ch, int z );
int cont_generateDigitUsing7SegDisplay( Iobj *obj, int number, int z );
int cont_generateTextAsConts( Iobj *obj, string text, Ipoint pos, float fontSize, int textAlign=TA_LEFT, bool smallCaps=false );
int cont_generateTextAreaAsConts( Iobj *obj, string text, Ipoint pos, float fontSize, int alignHoriz=TA_LEFT, int alignVert=TV_LEADING, bool smallCaps=false, float lineSpacing=4 );
int cont_addTwoPointContourToObj( Iobj *obj, float p1x, float p1y, float p2x, float p2y, float z, int open=1 );
int cont_addTwoPointContourToObj( Iobj *obj, Ipoint p1, Ipoint p2, int open=1 );
int cont_isEqual( Icont *cont1, Icont *cont2 );
int cont_doesPtExistInCont( Icont *cont, Ipoint *pt );
bool cont_getMBR( Icont *cont, Ipoint *ll, Ipoint *ur );
bool cont_getCenterOfMBR( Icont *cont, Ipoint *rpt );
bool cont_getCentroid( Icont *cont, Ipoint *rpt );
bool cont_insideCrude( Icont *cont1, Icont *cont2 );
float cont_getRadius( Icont *c );
void cont_findClosestPtInContToGivenPt( Ipoint *pt, Icont *cont, float *closestDist, Ipoint *closestPt, int *closestPtIdx );
void cont_findMinMaxAndAvgDistFromPt( Ipoint *pt, Icont *cont, float zScale, float &minDist, float &maxDist, float &avgDist );
bool cont_doContsTouch( Icont *cont1, Icont *cont2 );
float cont_minDistPtAndContourPts2D( Ipoint *pt, Icont *cont, bool returnZeroIfPtInside );
float cont_minDistBetweenContPts2D( Icont *cont1, Icont *cont2, bool returnZeroIfTouch );
float cont_minDistBetweenContPts3D( Icont *cont1, Icont *cont2, float zScale, Ipoint *pt1, Ipoint *pt2 );
void cont_reorderPtsToStartAtIdx( Icont *c, int idxNewFirstPt );
int cont_removePointsInCircle( Icont *cont, Ipoint *center, float radius, bool checkZEachPoint );
void cont_translate( Icont *cont, Ipoint *translate );
void cont_translate( Icont *cont, float x, float y );
void cont_rotateAroundPoint2D( Icont *cont, Ipoint *center, float angle );
void cont_scaleAboutPtXY( Icont *cont, Ipoint *center, float scaleX, float scaleY );
void cont_scaleAboutPt3D( Icont *cont, Ipoint *center, float scaleX, float scaleY, float scaleFactorZ );
void cont_scaleAboutPt( Icont *cont, Ipoint *pt, const float scaleFactor, bool ignoreZ );
void cont_stretchAlongAngle( Icont *cont, Ipoint *center, float angle, float stretchFactor );
void cont_generateCircle( Icont *cont, float radius, int numPoints, Ipoint center, bool addEndPt );
void cont_generateBox( Icont *cont, float llX, float llY, float width, float height, float z, bool repeatFirstPt=true );
void cont_generateBox( Icont *cont, Ipoint ll, Ipoint ur, bool repeatFirstPt=true );
int cont_numTimesCountsCross( Icont *cont1, Icont *cont2 );
bool cont_doCountoursCross( Icont *cont1, Icont *cont2, bool cont1Closed, bool cont2Closed );
bool cont_doCountoursCrossAndWhere( Icont *cont1, Icont *cont2, bool cont1Closed, bool cont2Closed, int *pt1BeforeCross, int *pt2BeforeCross );
bool cont_doesPtTouchContLine( Ipoint *pt, Icont *cont );
float cont_findClosestPts2D( Icont *cont1, Icont *cont2, int *closestPtIdxInCont1, int *closestPtIdxInCont2 );
void cont_reversePts( Icont *c );
void cont_addChamferPts( Icont *cont, Ipoint *ptPrev, Ipoint *ptCurr, Ipoint *ptNext, float distOffset, float minAngle );
void cont_expandOpenCont( Icont *contOrig, Icont *contReturn, float thickness, float minAngleForChamfers, bool roundEnds );
void cont_expandClosedCont( Icont *contOrig, Icont *innerCont, Icont *outerCont, float thickness, float minAngleForChamfers );
void cont_concat( Icont *contNew, Icont *cont1, Icont *cont2, bool matchClosestEnds );
int cont_addPtsCrude( Icont *cont, float maxDist, bool closed ); // MODIFIED
int cont_addPtsSmoothIteration( Icont *cont, float maxDist, float tensileFract, bool closed );
int cont_addPtsSmooth( Icont *cont, float maxDist, float tensileFract, bool closed,
bool roundZOpenPts=true, bool addPtEveryZ=true );
int cont_avgPtsPos( Icont *cont, float moveFract, float minDistToMove, bool closed, bool rescale );
int cont_reducePtsCrude( Icont *cont, float minDist, bool closed );
int cont_reducePtsTol(Icont *cont, float tol);
int cont_reducePtsMinArea( Icont *cont, float minArea, bool closed ); // MODIFIED
int cont_removeRedundantPts( Icont *cont, bool removeStraightLinePts, bool closed, bool removePts );
bool cont_isSimple( Icont *cont, bool closed=true );
bool cont_isSimpleSeg( Icont *cont, bool closed, int *ptCross );
void cont_makeSimple( Icont *cont );
int cont_breakIntoSimple( vector<IcontPtr> &conts, Icont *cont );
int cont_killVertAndHorzSegments( Icont *cont );
void cont_nudgeAllPointsRandomly( Icont *cont );
int cont_nudgeAnyPointsLyingOnOtherContour( Icont *cont1, Icont *cont2 );
bool cont_isConvex( Icont *cont );
int cont_makeConvexCrude( Icont *cont );
int cont_makeConvex( Icont *cont );
int cont_markConvexPtsNegOne( Icont *cont );
void cont_calcConvexProperties( Icont *cont, bool closed, int *numConvexPts, float *convexLen, float *hullLen, float *hullArea );
bool cont_breakContourEitherSide( Icont *cont, Icont *contBreak1, Icont *contBreak2, int idxPt1, int idxPt2, bool shareEdge );
bool cont_breakContourByLine( Icont *cont, Icont *contBreak1, Icont *contBreak2, Ipoint *linePt1, Ipoint *linePt2, Ipoint expectedPt, bool useExpectedPtInsteadOfMaxAreaSmallerSide );
int cont_breakContourByContour( vector<IcontPtr> &contSegs, Icont *contO, Icont *contBreakO, float minDist );
void cont_joinContsAtClosestApproach( Icont *newCont, Icont *cont1, Icont *cont2, bool addPtInMiddle );
void cont_joinContsAtClosestApproachVector( Icont *newCont, vector<IcontPtr> conts, bool addPtInMiddle );
void cont_getIntersectingConvexPolygon( Icont *newCont, Icont *cont1, Icont *cont2 );
int cont_breakContByZValue( Icont *contOrig, vector<IcontPtr> &contSegs, int zValue, bool removeOffSegments );
int cont_breakOpenContAtZValue( Icont *contOrig, vector<IcontPtr> &contSegs, int zValueToBreak );
int cont_breakContByCircle( Icont *contOrig, vector<IcontPtr> &contSegs, Ipoint *center, float radius );
int cont_addPtsAtIntersection( Icont *cont1, Icont *cont2 );
int cont_getIntersectingSegments( Icont *cont1, Icont *cont2, vector<IcontPtr> &cont1Segs, vector<IcontPtr> &cont2Segs );
int cont_getIntersectingSegmentsSafe( Icont *cont1, Icont *cont2, vector<IcontPtr> &cont1Segs, vector<IcontPtr> &cont2Segs );
int cont_getIntersectingPolygons( vector<IcontPtr> &finalConts, Icont *cont1, Icont *cont2 );
int cont_getUnionPolygons( vector<IcontPtr> &finalConts, Icont *cont1, Icont *cont2 );
bool cont_getOuterUnionPolygon( Icont *newCont, Icont *cont1O, Icont *cont2O );
Ipoint cont_getPtDistAlongLength( Icont *cont, float dist, bool closed, int startPt );
Ipoint cont_getPtFractAlongLength( Icont *cont, float fract, bool closed, int startPt );
vector<float> cont_getFractPtsAlongLength( Icont *cont, bool closed, int startPt );
int cont_addPtsFractsAlongLength( Icont *cont, Icont *contNew, vector<float> fractsAlongLen, bool closed, bool keepExistingPts, int startPt );
//void cont_calcPtsizeInfo( Iobj *obj, Icont *cont, const float zScale, float &openLength, float &fullLength, float &openVol, float &fullVol, float &avgRadius, float &minRadius, float &maxRadius, float &minMidRadius, float &surfaceArea );
void cont_calcPtsizeInfo( Iobj *obj, Icont *cont, const float zScale, float pixelSize,
float &openLength, float &fullLength,
float &avgR, float &firstR, float &lastR,
float &minR, float &maxR, float &minMidR,
float &openVol, float &fullVol, float &surfaceArea );
//############################################################
//-------------------------------
//## INLINE FUNCTION DEFINITIONS:
//------------------------
//-- Returns true if the contour exists and has points.
inline bool isContValid(Icont *cont)
{
return ( cont != NULL && imodContourGetMaxPoint(cont) > 0 );
}
//------------------------
//-- Returns true if the contour's interpolated flag is set.
inline bool isInterpolated(Icont *cont)
{
return ( imodContourGetFlag( cont, ICONT_STIPPLED ) != 0 );
}
//------------------------
//-- Sets the contour's interpolated flag to on or off.
inline void setInterpolated(Icont *cont, int on)
{
imodContourSetFlag( cont, ICONT_STIPPLED, on );
}
//------------------------
//-- Returns true if the contour's delete flag is set.
inline bool isOpenFlag(Icont *cont)
{
return ( imodContourGetFlag( cont, ICONT_OPEN ) != 0 );
}
//------------------------
//-- Sets the contour's delete flag ("ICONT_TEMPUSE") to on or off.
inline void setOpenFlag(Icont *cont, int on )
{
imodContourSetFlag( cont, ICONT_OPEN, on );
}
//------------------------
//-- Returns true if the contour's delete flag is set.
inline bool isDeleteFlag(Icont *cont)
{
return ( imodContourGetFlag( cont, (unsigned int)ICONT_TEMPUSE ) != 0 );
}
//------------------------
//-- Sets the contour's delete flag ("ICONT_TEMPUSE") to on or off.
inline void setDeleteFlag(Icont *cont, int on )
{
imodContourSetFlag( cont, ICONT_TEMPUSE, on );
}
//------------------------
//-- Returns true is the contour has no points or is marked for deletion.
inline bool isEmpty(Icont *cont)
{
return ( imodContourGetMaxPoint( cont ) == 0 || isDeleteFlag( cont ) );
}
//------------------------
//-- Shorter function name for "imodContourGetMaxPoint()"
inline int psize(Icont *cont)
{
return imodContourGetMaxPoint( cont );
}
//------------------------
//-- Returns a pointer to the point at the given index in the contour...
//-- wrapping around to the first point if idx is >= imodContourGetMaxPoint().
inline Ipoint *getPt(Icont *cont, int idx)
{
int contPts = imodContourGetMaxPoint( cont );
idx = (idx+contPts) % contPts;
return imodContourGetPoint(cont, idx);
//Ipoint *pts = imodContourGetPoints( cont );
//return &pts[ (idx+contPts) % contPts ];
}
//------------------------
//-- Returns the Z value of the point at the given index in the contour
//-- rounded to the nearest integer
inline int getPtZInt(Icont *cont, int idx)
{
return int( getPt(cont, idx)->z + 0.5 );
}
//------------------------
//-- Returns a pointer to the point at the given index in the contour...
//-- but does not wrap around.
inline Ipoint *getPtNoWrap(Icont *cont, int idx)
{
int contPts = imodContourGetMaxPoint( cont );
keepWithinRange( idx, 0, contPts-1 );
return imodContourGetPoint(cont, idx);
}
//------------------------
//-- Returns a pointer to the last point in the contour.
inline Ipoint *getLastPt(Icont *cont )
{
if(!cont)
return (NULL);
int contPts = imodContourGetMaxPoint( cont );
Ipoint *pts = imodContourGetPoints( cont );
return &pts[ contPts-1 ];
}
//------------------------
//-- Returns a pointer to the first point in the contour.
inline Ipoint *getFirstPt(Icont *cont )
{
if(!cont)
return (NULL);
Ipoint *pts = imodContourGetPoints( cont );
return &pts[ 0 ];
}
//------------------------
//-- Used to set the three coordinates of a point in one call.
inline void setPt(Ipoint *pt, float x, float y, float z)
{
pt->x = x;
pt->y = y;
pt->z = z;
}
//------------------------
//-- Used to set the three coordinates of a point in one call.
inline void copyPt(Ipoint *to, Ipoint *from)
{
to->x = from->x;
to->y = from->y;
to->z = from->z;
}
//------------------------
//-- Used to create a new point in one call.
inline Ipoint newPt( float x, float y, float z )
{
Ipoint pt;
pt.x = x;
pt.y = y;
pt.z = z;
return pt;
}
//------------------------
//-- Returns true if two points are equal... but unlike imodPointIsEqual
//-- only checks x and y values
inline bool ptsEqual( Ipoint *pt1, Ipoint *pt2 )
{
return (pt1->x == pt2->x) && (pt1->y == pt2->y);
}
//------------------------
//-- Returns true if two points are with plus or minus "prec"
//-- of each other in X and Y
inline bool ptsApproxEqual( Ipoint *pt1, Ipoint *pt2, float prec )
{
return (pt1->x >= pt2->x-prec && pt1->x <= pt2->x+prec )
&& (pt1->y >= pt2->y-prec && pt1->y <= pt2->y+prec );
}
//------------------------
//-- Resets the point size and fine grain of all points in the current contour
//-- to the default.
inline void removeExtraInfo( Icont *cont )
{
int openFlag = imodContourGetFlag( cont, ICONT_OPEN );
int stippledFlag = imodContourGetFlag( cont, ICONT_STIPPLED );
Icont *tempCont = imodContourDup( cont );
imodContourDefault( cont );
imodContourCopy( tempCont, cont );
cont_copyPts( tempCont, cont, true );
imodContourSetFlag( cont, ICONT_OPEN, openFlag );
imodContourSetFlag( cont, ICONT_STIPPLED, stippledFlag );
imodContourDelete( tempCont );
}
//------------------------
//-- Resets the point size and fine grain of all points in the current contour
//-- to the default.
inline void removePtsSize( Icont *cont )
{
for (int p=0; p<psize(cont); p++)
removePtSize( cont, p );
}
//------------------------
//-- Resets the point size and fine grain info of the given point in the current
//-- contour to the default.
inline void removePtSize( Icont *cont, int idx )
{
Ipoint *pt = getPt(cont, idx);
Ipoint newPt;
setPt( &newPt, pt->x, pt->y, pt->z );
imodPointAdd( cont, &newPt, idx );
imodPointDelete( cont, idx+1 );
}
//------------------------
//-- Returns true if the size of the point equals the default sphere size for the object
//-- or returns false if different.
inline bool isDefaultSize( Iobj *obj, Icont *cont, int idx )
{
float ptSize = imodPointGetSize(obj, cont, idx);
float objSphereSize = imodObjectGetValue(obj,IobjPointSize);
return (ptSize == objSphereSize);
}
//------------------------
//-- Prints some simple information about point - used in debugging.
inline void printPt( Ipoint *pt )
{
wprint("pt = %f,%f,%f\n", pt->x, pt->y, pt->z );
}
//------------------------
//-- Prints some simple information about the contour - used in debugging.
inline void printCont( Icont *cont )
{
wprint("cont pts=%d z=%d \n", psize(cont), imodContourZValue(cont) );
}
//------------------------
//-- Prints all points in the contour.
inline void printContPts( Icont *cont )
{
wprint("cont: " );
for(int p=0; p<psize(cont); p++)
wprint(" { %d,%d,%d }",
(int)getPt(cont,p)->x, (int)getPt(cont,p)->y, (int)getPt(cont,p)->z );
wprint("\n");
}
//------------------------
//-- Useful for deleting all points without calling "imodContourDefault",
//-- which will reset other contour properties.
inline void deleteAllPts( Icont *cont )
{
int nPoints = psize( cont );
for (int i=nPoints-1; i>=0; i--)
imodPointDelete( cont, i );
}
//------------------------
//-- Shorter function name for "imodContourZValue()"
inline float getZ( Icont *cont )
{
return imodContourZValue( cont );
}
//------------------------
//-- Returns Z value of contour rounded to nearest integer
inline int getZInt( Icont *cont )
{
return floor(imodContourZValue( cont ) + 0.5);
}
//------------------------
//-- Returns the span of the contour in Z (will be 0 if contour on one slice)
inline float getZRange( Icont *cont )
{
if( isEmpty(cont) )
return 0;
float maxZ = INT_MIN;
float minZ = INT_MAX;
for(int p=0; p<psize(cont); p++)
{
float z = getPt(cont,p)->z;
updateMin( minZ, z );
updateMax( maxZ, z );
}
return (maxZ - minZ);
}
//------------------------
//-- Changes the z value of all the points to the given value.
inline void setZValue( Icont *cont, int newZValue )
{
for (int i=0; i<psize( cont ); i++)
getPt( cont, i )->z = (float)newZValue;
imodContourSetFlag( cont, ICONT_WILD, 0 );
//imodel_contour_check_wild(cont)
}
//------------------------
//-- Copies all the points from the first contour to the second.
inline void cont_copyPts( Icont *from, Icont *to, bool clearToCont )
{
if(clearToCont)
imodContourDefault(to);
for (int i=0; i<psize(from); i++)
imodPointAppend( to, getPt(from,i) );
}
//------------------------
//-- Deletes all contours in this array by calling
//-- imodContourDelete for any non-null contour pointers
inline void deleteContours( vector<IcontPtr> &conts )
{
for (int iCont=0; iCont<(int)conts.size(); iCont++)
conts[iCont].deleteContour();
}
//------------------------
//-- Erases the specified contours in this array by calling
//-- imodContourDelete and then removing it from the array
inline void eraseContour( vector<IcontPtr> &conts, int idx )
{
conts[idx].deleteContour();
conts.erase( conts.begin() + idx );
}
//------------------------
//-- Shorter function name for "imodObjectGetContour()"
inline Icont* getCont( Iobj *obj, int idx )
{
return imodObjectGetContour(obj, idx);
}
//------------------------
//-- Returns the last contour in the object
inline Icont* getLastCont( Iobj *obj )
{
if(!obj)
return (NULL);
return imodObjectGetContour(obj, csize(obj)-1 );
}
//------------------------
//-- Shorter function name for "imodObjectGetMaxContour()"
inline int csize( Iobj *obj )
{
return imodObjectGetMaxContour(obj);
}
//------------------------
//-- Shorter function name for "imodGetMaxObject()"
inline int osize( Imod *imod )
{
return imodGetMaxObject(imod);
}
//---------------------------------
//-- Returns the object at the specified index (or NULL if there is none).
inline Iobj* getObj( Imod *imod, int idx )
{
int objIdx, contIdx, ptIdx;
imodGetIndex(imod, &objIdx, &contIdx, &ptIdx);
imodSetIndex( imod, idx, 0, 0);
Iobj *obj = imodObjectGet(imod);
imodSetIndex( imod, objIdx, contIdx, ptIdx );
return obj;
}
//---------------------------------
//-- Returns true if the object is closed.
inline bool isObjClosed(Iobj *obj)
{
return (imodObjectGetValue(obj, IobjFlagClosed) == 1);
}
//---------------------------------
//-- Returns true if the object is valid and has it's draw flag on.
inline bool isContClosed(Iobj *obj, Icont *cont)
{
bool objClosed = isObjClosed( obj );
bool contClosed = !isOpenFlag( cont );
return ( objClosed && contClosed );
}
//---------------------------------
//-- Returns true if the object is valid and has it's draw flag on.
inline bool isObjectValidAndShown(Iobj *obj)
{
// imodObjectGetValue returns 1 when object is hidden
int objHidden = (imodObjectGetValue(obj, IobjFlagDraw) );
return ( obj != NULL && !objHidden );
}
//---------------------------------
//-- Wrapper function for "imodObjectSetColor" but using integers between 0 & 255
//-- intead of a float beteen 0 & 1.
//-- Sets the objects color using red, green, blue integers with values between
//-- 0 and 255. I created this function because I had problems with
//-- "imodObjectSetColor( obj, red/255, green/255, blue/255 );" sometimes
//-- typecasting some colors incorrectly.
inline void setObjColor( Iobj *obj, int red, int green, int blue )
{
float redF = (float)red / 255.0f;
float greenF = (float)green / 255.0f;
float blueF = (float)blue / 255.0f;
imodObjectSetColor( obj, redF,greenF,blueF );
}
//############################################################
#endif
icontextra.cpp
/*
* icontextra.c -- Special file for contour functions I need to use in my
* special plugins and do not exist in icont.h/icont.cpp
*/
/*****************************************************************************
* Copyright (C) 2007 by Andrew Noske from the Institute for Molecular *
* Bioscience at the University of Queensland (Australia) *
*****************************************************************************/
/* include needed Qt headers and imod headers */
#include "icontextra.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "_common_functions.h"
#include "imodplugin.h"
#include <vector>
using namespace std;
//############################################################
//----------------------------------------------------------------------------
//## SPECIAL DATA STRUCTURES USED IN FUNCTIONS:
//----------------------------------------------------------------------------
//## IcontPtr
IcontPtr::IcontPtr( Icont* cont_ )
{
cont = imodContourDup( cont_ );
}
IcontPtr::IcontPtr()
{
cont = imodContourNew();
}
void IcontPtr::deleteContour()
{
if( cont != NULL )
{
imodContourDelete( cont );
cont = NULL;
}
}
//## PtConnection // used especially for getIntersectingPolygon function
PtConnection::PtConnection( Ipoint _intercept ) {
intercept = _intercept;
cont1Idx = -1;
cont2Idx = -1;
included = false;
}
PtConnection::PtConnection( Ipoint _intercept, int _cont1Idx ) {
intercept = _intercept;
cont1Idx = _cont1Idx;
cont2Idx = -1;
included = false;
}
//-- Used to sort points by y value - using x as a tie-breaker.
bool operator< (const PtConnection &a, const PtConnection &b) {
if ( a.intercept.y == b.intercept.y )
return ( a.intercept.x <= b.intercept.x );
return (a.intercept.y < b.intercept.y);
}
//-- Used to eliminate duplicates
bool operator== (const PtConnection lhs, const PtConnection rhs) {
return ( lhs.intercept.x == rhs.intercept.x );
}
//## IdxAndFloat // used especially for getIntersectingPolygon function
IdxAndFloat::IdxAndFloat(int _idx, float _dist) {
idx = _idx;
dist = _dist;
}
//-- Test if distance is "less than"
bool operator<(const IdxAndFloat &a, const IdxAndFloat &b) {
return (a.dist < b.dist);
}
//## IntAndInt // used to connect point indexes between two contours
IntAndInt::IntAndInt(int _idx1, int _idx2) {
idx1 = _idx1;
idx2 = _idx2;
}
//-- Test if indexes are "less than"
bool operator<(const IntAndInt &a, const IntAndInt &b) {
return (a.idx1 < b.idx1) || (a.idx1 == b.idx1 && a.idx2 < b.idx2);
}
//-- Test if both indexes are "equal"
bool operator==(const IntAndInt lhs, const IntAndInt rhs) {
return (lhs.idx1 == rhs.idx1 && lhs.idx2 == rhs.idx2);
}
//## IdxToSort // used especially for getIntersectingPolygon function
IdxToSort::IdxToSort() {
idx = 0;
float1 = 0;
float2 = 0;
}
IdxToSort::IdxToSort( int _idx, float _float1 ) {
idx = _idx;
float1 = _float1;
float2 = 0;
}
IdxToSort::IdxToSort( int _idx, float _float1, float _float2 ) {
idx = _idx;
float1 = _float1;
float2 = _float2;
}
//-- Test if distance is "less than"
bool operator<(const IdxToSort &a, const IdxToSort &b) {
if ( a.float1 == b.float1 )
return ( a.float2 <= b.float2 );
return (a.float1 < b.float1);
}
//############################################################
//----------------------------------------------------------------------------
//## POINT RELATED FUNCTIONS:
//----------------------------------------------------------------------------
//------------------------
//-- Rotates the point by "theta" around the "center" in X and Y
void point_rotatePointAroundPoint2D( Ipoint *pt, Ipoint *center, float theta )
{
pt->x -= center->x; //|-- translate point (so center is at origin)
pt->y -= center->y; //|
float newX = pt->x*cos(theta) - pt->y*sin(theta); //|-- determine new position
float newY = pt->x*sin(theta) + pt->y*cos(theta); //| after rotation
pt->x = newX;
pt->y = newY;
pt->x += center->x; //|-- translate back point
pt->y += center->y; //| (so center returns to its original position)
}
//------------------------
//-- "Scales" the point by shifting it away from the center along each axis
void point_scalePtAboutPt( Ipoint *pt, Ipoint *center,
float scaleX, float scaleY, float scaleZ )
{
pt->x -= center->x; //|-- translate point (so center is at origin)
pt->y -= center->y; //|
pt->z -= center->z; //|
pt->x *= scaleX; // "scale"/shift point in x
pt->y *= scaleY; // "scale"/shift point in y
pt->z *= scaleZ; // "scale"/shift point in y
pt->x += center->x; //|-- translate back point
pt->y += center->y; //| (so center returns to its original position)
pt->z += center->z; //|
}
//------------------------
//-- "Scales" the point by shifting it away from the center by
//-- scaleX in X and scaleY in Y
void point_scalePtAboutPt2D( Ipoint *pt, Ipoint *center, float scaleX, float scaleY )
{
pt->x -= center->x; //|-- translate point (so center is at origin)
pt->y -= center->y; //|
pt->x *= scaleX; // "scale"/shift point in x
pt->y *= scaleY; // "scale"/shift point in y
pt->x += center->x; //|-- translate back point
pt->y += center->y; //| (so center returns to its original position)
}
//------------------------
//-- Takes a series of points, and return their average distance from the startPt,
//-- and also records the distance of each in the radius value of the points.
float point_findAvgDistToPt( Icont *pts, Ipoint *startPt, Ipoint *scalePt )
{
float totalDist = 0.0; // tallys the total distance from each point to the start point
//## FIND AVERAGE DISTANCE FROM THE POINTS TO THE START POINT:
for(int p=0; p<psize(pts); p++ )
{
double dist = imodPoint3DScaleDistance( getPt(pts,p), startPt, scalePt );
totalDist += dist;
imodPointSetSize( pts, p, dist ); // store distance to point in point radius
//pts[p].radius = dist; // NOTE: stores distance as points radius so we don't have to calculate a second time.
}
return ( fDiv( totalDist, psize(pts) ) );
}
//------------------------
//-- Takes a series of points, and return the "maxDeviation" as the furthest distance
//-- minus the closest distance, divided by the furthest distance.
float point_maxDiffInDistToPoint( Icont *pts, Ipoint *startPt, Ipoint *scalePt )
{
float minDist = FLOAT_MAX;
float maxDist = 0.0;
//## FIND CLOSEST AND FURTHEST DISTANCE FROM A POINT:
for(int p=0; p<psize(pts); p++ )
{
float dist = imodPoint3DScaleDistance( getPt(pts,p), startPt, scalePt );
minDist = MIN( minDist, dist );
maxDist = MAX( maxDist, dist );
}
return ( fDiv( (maxDist-minDist), maxDist ) );
}
//------------------------
//-- Takes a series of points, and return the standard devaition these are
//-- from the "startPt".
float point_stdDevDistToPoint( Icont *pts, Ipoint *startPt, Ipoint *scalePt )
{
float avgDist = point_findAvgDistToPt( pts, startPt, scalePt );
float sumDifferenceSquared = 0.0;
//## FIND CLOSEST AND FURTHEST DISTANCE FROM A POINT:
for(int p=0; p<psize(pts); p++ )
{
float dist = imodPoint3DScaleDistance( getPt(pts,p), startPt, scalePt );
sumDifferenceSquared += SQ(dist - avgDist);
}
return sqrt( fDiv( sumDifferenceSquared, (psize(pts) - 1) ) );
}
//------------------------
//-- Takes a series of points, and a starting points, and tries to find the
//-- "equidistant center" such that all points are an equal distance from it.
//-- In other words, it tries to find the best "center" of a sphere/circle,
//-- where the set of points are points around the edge of the sphere but they
//-- DON'T have to be evenly distributed. It also sets the radius to the average
//-- distance to the points. The best starting point is one that is already
//-- inside the bounding polygon and/or close to the center.
//--
//-- NOTE: The illustration below demonstrates the difference between a center of points
//-- and the centroid using the same set of 4 points.
//--
//-- CENTER OF POINTS: | CENTER OF AREA: |
//-- | |
//-- 2 | 2 |
//-- 1 o 3 | 1 __-o-__ 3 |
//-- o | o | o- -o |
//-- \ | / | \__ \ center of area |
//-- \ | / | \ X \ <-- centroid |
//-- startPt -> X \ | / | \__ \ (center of area) |
//-- \|/ | \___ \ |
//-- center -> X---------o 4 | \__o 4 |
//-- of points | |
//-- | |
Ipoint point_findCenterOfPts( Iobj *obj, Icont *pts, Ipoint *startPt, float &avgRadius,
float zScale, float &avgResidual, int maxIts )
{
int numPoints = psize(pts);
if( numPoints <= 0 )
return (*startPt);
Ipoint scalePt;
setPt( &scalePt, 1, 1, zScale );
float prevAvgDist = FLOAT_MAX;
Ipoint prevStartPt;
setPt( &prevStartPt, startPt->x, startPt->y, startPt->z);
for(int i=0; i<maxIts; i++)
{
//## FIND AVERAGE DISTANCE FROM THE POINTS TO THE START POINT:
float totalDist = 0.0; // tallys the total distance from each point to the start point
for(int p=0; p<psize(pts); p++ )
{
float dist = imodPoint3DScaleDistance( getPt(pts,p), startPt, &scalePt );
totalDist += dist;
imodPointSetSize( pts, p, dist );
//pts[p].radius = dist;
// NOTE: stores the distance as the points radius so we don't have to calculate a second time.
}
float avgDist = fDiv( totalDist, numPoints );
//## CHECK IF IMPROVEMENT HAS BEEN MADE:
float improvementDist = prevAvgDist - avgDist;
if( improvementDist == 0 ) // if no improvement has been made: early exit
{
break;
}
else if( improvementDist < 0 ) // if we are going the WRONG way from the center (which can happen if the start point is too far outside the sphereical profile):
{
//startPt = line_findPtFractBetweenPts( startPt, pts[0], 0.5 ); // put the start point on the other side
//continue;
}
prevAvgDist = avgDist;
prevStartPt = *startPt;
//## PROJECT LINES FROM EACH POINT THROUGH THE START POINT, AND AVERAGE THE PROJECTED POINTS TO FIND A NEW START POINT:
Ipoint totalProjectedVals; // tally's the total distance in X, Y and Z of the projected points, so that a averaged "middle point" can be found.
setPt( &totalProjectedVals, 0,0,0 );
float totalResidual = 0;
for(int p=0; p<numPoints; p++ )
{
float radius = imodPointGetSize(obj, pts, p);
float fractBetweenPts = fDivCustom( avgDist, radius, 1 );
Ipoint projectedPt = line_findPtFractBetweenPts( getPt(pts,p), startPt, fractBetweenPts );
totalProjectedVals.x += projectedPt.x;
totalProjectedVals.y += projectedPt.y;
totalProjectedVals.z += projectedPt.z;
float residual = ABS( avgDist - radius );
totalResidual += residual;
}
Ipoint newStartPt;
newStartPt.x = totalProjectedVals.x / numPoints;
newStartPt.y = totalProjectedVals.y / numPoints;
newStartPt.z = totalProjectedVals.z / numPoints;
avgRadius = avgDist;
avgResidual = totalResidual / numPoints;
*startPt = newStartPt;
}
return *startPt;
}
Ipoint point_findCenterOfPtsX( Iobj *obj, Icont *pts, Ipoint *startPt, float avgRadius, float zScale, int maxIts )
{
float avgResidual;
return point_findCenterOfPts( obj, pts, startPt, avgRadius, zScale, avgResidual, maxIts );
}
//------------------------
//-- Makes multiple calls to "point_findCenterOfPts" in an attempt to find the zScale
//-- and center point (for this zScale) which give the lowest average residual.
//--
//-- NOTE: This uses a system where five evenly spaced zScales (ranging from "low" to "high")
//-- are tested, and the value with the lowest residual becomes the new middle,
//-- and the increment reduces until the "best Z scale" is found.
//-- The diagram below shows how the guess gets closer and closer to the local minimum.
//-- Each round two more values have to be calculate.
//--
//--
//-- Low MidLow Mid MidHigh High
//-- round 1: | | | | |
//-- *
//-- round 2: | | | | |
//-- *
//-- round 3: | | | | |
//-- *
//-- \___ ____/
//-- \______ ____/ ^ avgResidual (Y AXIS)
//-- \________. ______/
//-- \_ _./ > zScale (X AXIS)
//-- \._/
//-- ^
int MAX_GUESSES = 50;
Ipoint point_findOptimalZScaleAndCenterOfPts( Iobj *obj, Icont *pts, Ipoint *startPt, float &avgR, float zScale,
float &bestZScale, float &bestAvgResidual, float startChangeZ,
float accuracyZScale, int maxIts )
{
float zScaleLow = zScale - startChangeZ;
float zScaleHigh = zScale + startChangeZ;
float zScaleMid = getFractBetween( zScaleLow, zScaleHigh, 0.5);
float zScaleMidLow = getFractBetween( zScaleLow, zScaleHigh, 0.25);
float zScaleMidHigh = getFractBetween( zScaleLow, zScaleHigh, 0.75);
float zScaleLow_resid = 0;
float zScaleHigh_resid = 0;
float zScaleMid_resid = 0;
float zScaleMidLow_resid = 0;
float zScaleMidHigh_resid = 0;
Ipoint centerMid = point_findCenterOfPts( obj, pts, startPt, avgR, zScaleMid, zScaleMid_resid, maxIts );
Ipoint centerLow = point_findCenterOfPts( obj, pts, ¢erMid, avgR, zScaleLow, zScaleLow_resid, maxIts );
Ipoint centerHigh = point_findCenterOfPts( obj, pts, ¢erMid, avgR, zScaleHigh, zScaleHigh_resid, maxIts );
Ipoint centerMidHigh = centerMid;
Ipoint centerMidLow = centerMid;
for (int i=0; i<MAX_GUESSES; i++)
{
//## CALCULATE THE NEW MID POINTS:
zScaleMidLow = getFractBetween( zScaleLow, zScaleHigh, 0.25);
zScaleMidHigh = getFractBetween( zScaleLow, zScaleHigh, 0.75);
centerMidHigh = point_findCenterOfPts( obj, pts, ¢erMid, avgR, zScaleMidLow, zScaleMidLow_resid, maxIts );
centerMidLow = point_findCenterOfPts( obj, pts, ¢erMid, avgR, zScaleMidHigh, zScaleMidHigh_resid, maxIts );
//## DETERMINE MINIMUM AND MAXIMUM AVERAGE RESIDUAL OVER FIVE Z-SCALES:
float min_resid = MIN(MIN(MIN(MIN( zScaleLow_resid, zScaleHigh_resid ), zScaleMid_resid), zScaleMidLow_resid), zScaleMidHigh_resid);
float max_resid = MAX(MAX(MAX(MAX( zScaleLow_resid, zScaleHigh_resid ), zScaleMid_resid), zScaleMidLow_resid), zScaleMidHigh_resid);
//## ADJUST FIVE POINTS SO THAT ZSCALE WITH THE LOWEST VALUE BECOMES THE NEW MID POINT:
if( min_resid == zScaleLow_resid )
{
zScaleHigh = zScaleMid; zScaleHigh_resid = zScaleMid_resid; centerHigh = centerMid;
zScaleMid = zScaleLow; zScaleMid_resid = zScaleLow_resid; centerMid = centerLow;
zScaleLow = getFractBetween( zScaleHigh, zScaleMid, 2.0 ); centerLow = point_findCenterOfPts( obj, pts, ¢erMid, avgR, zScaleLow, zScaleLow_resid, maxIts );
}
else if( min_resid == zScaleHigh_resid )
{
zScaleLow = zScaleMid; zScaleLow_resid = zScaleMid_resid; centerLow = centerMid;
zScaleMid = zScaleHigh; zScaleMid_resid = zScaleHigh_resid; centerMid = centerHigh;
zScaleHigh = getFractBetween( zScaleLow, zScaleMid, 2.0 ); centerHigh = point_findCenterOfPts( obj, pts, ¢erMid, avgR, zScaleHigh, zScaleHigh_resid, maxIts );
}
else if( min_resid == zScaleMid_resid )
{
zScaleLow = zScaleMidLow; zScaleLow_resid = zScaleMidLow_resid; centerLow = centerMidLow;
zScaleHigh = zScaleMidHigh; zScaleHigh_resid = zScaleMidHigh_resid; centerHigh = centerMidHigh;
}
else if( min_resid == zScaleMidLow_resid )
{
zScaleHigh = zScaleMid; zScaleHigh_resid = zScaleMid_resid; centerHigh = centerMid;
zScaleMid = zScaleMidLow; zScaleMid_resid = zScaleMidLow_resid; centerMid = centerMidLow;
}
else if( min_resid == zScaleMidHigh_resid )
{
zScaleLow = zScaleMid; zScaleLow_resid = zScaleMid_resid; centerLow = centerMid;
zScaleMid = zScaleMidHigh; zScaleMid_resid = zScaleMidHigh_resid; centerMid = centerMidHigh;
}
//cout << " i=" << i << " zScaleMid=" << zScaleMid << "low=" << zScaleLow << "high=" << zScaleHigh << " min_resid=" << min_resid << " diff=" << (max_resid-min_resid) << endl; //%%%%%%
//## CHECK IF SUFFICIENT ACCURACY TO EXIT EARLY:
if( (zScaleHigh - zScaleLow) < accuracyZScale )
break;
}
bestZScale = zScaleMid;
bestAvgResidual = zScaleMid_resid;
return centerMid;
}
//------------------------
//-- Calculates a value "fract" of the distance between "p1" and "p2" key values
//-- according to the cardinal spline algorithm.
//-- It does this by first calculating the gradients at these two key values.
//-- See: http://www.mvps.org/directx/articles/catmull/
//--
//-- NOTE: This function is used in smooth interpolation.
//--
//-- p0 p1 p2 p3
//-- o------------o-----#------o---------------------o
//-- |
//-- desired point
//--
//-- ACTUAL EQUATION:
//--
//-- q(t) = 0.5 * [ ( -p0 + 3*p1 - 3*p2 + p3 ) * t^3
//-- + ( 2*p0 - 5*p1 + 4*p2 - p3 ) * t^2
//-- + ( -p0 + p2 ) * t
//-- + ( 2*p1 ) ]
float getValCardinalSpline( float fract, float p0, float p1, float p2, float p3,
float tensileFract )
{
//## CALCULATE "FRACTION INTO KEYFRAME" SQUARED AND CUBED:
float fract_2 = fract * fract;
float fract_3 = fract * fract * fract;
//## CALCULATE VALUES:
//calculate gradient at previous keyframe:
float gPrev = tensileFract*(p1 - p0)/1.0
+ tensileFract*(p2 - p1)/1.0;
//calculate gradient at current keyframe:
float gCurr = tensileFract*(p2 - p1)/1.0
+ tensileFract*(p3 - p2)/1.0;
//calculate value for current frame:
float curr = p1*(2*fract_3 - 3*fract_2 + 1)
+ p2*(3*fract_2 - 2*fract_3)
+ gPrev*(fract_3 - 2*fract_2 + fract)
+ gCurr*(fract_3 - fract_2);
return curr;
// return value on spline:
return 0.5 * ( ( 2 * p1 )
+ (-p0 + p2) * fract
+ (2*p0 - 5*p1 + 4*p2 - p3) * fract_2
+ (-p0 + 3*p1 - 3*p2 + p3) * fract_3 );
}
//------------------------
//-- Calculates point "fract" of the distance between points "p1" and "p2"
//-- according to the catmull-rom spline algorithm.
//-- See: getValCardinalSpline()
Ipoint getPtCardinalSpline( float fract, Ipoint p0, Ipoint p1, Ipoint p2,
Ipoint p3, float tensileFract )
{
Ipoint returnPt;
returnPt.x = getValCardinalSpline( fract, p0.x,p1.x,p2.x,p3.x, tensileFract );
returnPt.y = getValCardinalSpline( fract, p0.y,p1.y,p2.y,p3.y, tensileFract );
returnPt.z = getValCardinalSpline( fract, p0.z,p1.z,p2.z,p3.z, tensileFract );
return (returnPt);
}
//############################################################
//----------------------------------------------------------------------------
//## MINIMUM BOUNDING RECTANGLE (MBR) RELATED FUNCTIONS:
//----------------------------------------------------------------------------
//------------------------
//-- Sets MBR to default - with LL point all max value and UR point min value
//-- ready to add points with mbr_addPt.
void mbr_reset( Ipoint *ll, Ipoint *ur )
{
setPt( ll, FLOAT_MAX, FLOAT_MAX, FLOAT_MAX );
setPt( ur, FLOAT_MIN, FLOAT_MIN, FLOAT_MIN );
}
//------------------------
//-- Adjusts the MBR points so they include the point "pt" provided.
void mbr_addPt( Ipoint *pt, Ipoint *ll, Ipoint *ur )
{
if( pt->x < ll->x ) ll->x = pt->x;
if( pt->y < ll->y ) ll->y = pt->y;
if( pt->z < ll->z ) ll->z = pt->z;
if( pt->x > ur->x ) ur->x = pt->x;
if( pt->y > ur->y ) ur->y = pt->y;
if( pt->z > ur->z ) ur->z = pt->z;
}
//------------------------
//-- Used to calculate distance between a point and 2 edges.
//-- (but does NOT check wrap around)
float mbr_distToNearestEdge(float val, float min, float max)
{
if ( val < min ) return (min - val);
else if ( val > max ) return (max - val);
else return (0.0);
}
//------------------------
//-- Used to calculate minimum distance between two edges (which may overlap).
float mbr_distToNearestEdge(float min1, float max1, float min2, float max2)
{
if ( max1 < min2 ) return (min2 - max1);
else if ( min1 > max2 ) return (max2 - min1);
else return (0.0);
}
//------------------------
//-- Used to calculate minimum distance between two bounding boxes (which may overlap)
//-- in X and Y
float mbr_distBetweenBBoxes2D(Ipoint *ll1, Ipoint *ur1, Ipoint *ll2, Ipoint *ur2)
{
float distX = mbr_distToNearestEdge(ll1->x,ur1->x,ll2->x,ur2->x);
float distY = mbr_distToNearestEdge(ll1->y,ur1->y,ll2->y,ur2->y);
return sqrt( distX*distX + distY*distY );
}
//------------------------
//-- Used to calculate minimum distance between two bounding boxes (which may overlap)
//-- in 3D
float mbr_distBetweenBBoxes3D(Ipoint *ll1, Ipoint *ur1, Ipoint *ll2, Ipoint *ur2, float zScale)
{
float distX = mbr_distToNearestEdge(ll1->x,ur1->x,ll2->x,ur2->x);
float distY = mbr_distToNearestEdge(ll1->y,ur1->y,ll2->y,ur2->y);
float distZ = mbr_distToNearestEdge(ll1->z,ur1->z,ll2->z,ur2->z) * zScale;
float distXY = sqrt( distX*distX + distY*distY );
return sqrt( distXY*distXY + distZ*distZ );
}
//------------------------
//-- Used to calculate minimum distance between a point and a bounding boxing X and Y
bool mbr_distToBBox2D(Ipoint *pt, Ipoint *ll, Ipoint *ur)
{
float distX = mbr_distToNearestEdge( pt->x, ll->x, ur->x );
float distY = mbr_distToNearestEdge( pt->y, ll->y, ur->y );
return sqrt( distX*distX + distY*distY );
}
//------------------------
//-- Checks if point is inside the given bounding box
bool mbr_isPtInsideBBox(Ipoint *pt, Ipoint *ll, Ipoint *ur)
{
return ( isBetween(ll->z,pt->z,ur->z)
&& isBetween(ll->x,pt->x,ur->x)
&& isBetween(ll->y,pt->y,ur->y) );
}
//------------------------
//-- Checks if point is inside the given bounding box
//-- without considering Z
bool mbr_isPtInsideBBox2D(Ipoint *pt, Ipoint *ll, Ipoint *ur)
{
return ( isBetween(ll->x,pt->x,ur->x) && isBetween(ll->y,pt->y,ur->y));
}
//------------------------
//-- Used to calculate MINDIST between two edges.
bool mbr_doEdgesOverlap(float min1, float max1, float min2, float max2)
{
return !( (max1 < min2) || (min1 > max2) ); // says: if edge1 to the LEFT or RIGHT of edge2: it does not overlap (otherwise it does)
}
//------------------------
//-- Returns true if two bounding boxes overlap in 2D
bool mbr_doBBoxesOverlap2D(Ipoint *p1ll, Ipoint *p1ur, Ipoint *p2ll, Ipoint *p2ur)
{
return ( mbr_doEdgesOverlap( p1ll->x, p1ur->x, p2ll->x, p2ur->x )
&& mbr_doEdgesOverlap( p1ll->y, p1ur->y, p2ll->y, p2ur->y ) );
}
//------------------------
//-- Returns true if the first box is completely inside the second box
//-- (along all dimensions)
bool mbr_isBBoxInsideBBox(Ipoint *ll1, Ipoint *ur1, Ipoint *ll2, Ipoint *ur2)
{
return ( ( ll1->x >= ll2->x && ur1->x <= ur2->x )
&& ( ll1->y >= ll2->y && ur1->y <= ur2->y )
&& ( ll1->z >= ll2->z && ur1->z <= ur2->z ) );
}
//############################################################
//----------------------------------------------------------------------------
//## LINE RELATED FUNCTIONS:
//----------------------------------------------------------------------------
//------------------------
//-- Calculates the angle (in degrees) of a line from the horizontal
//-- (the angle formed when a ray projected left of the first point)
//-- NOTE: Returned angle will be between -180 and 180 degrees.
float line_getAngle2D ( Ipoint *linept1, Ipoint *linept2 )
{
float oppY = linept2->y - linept1->y;
float adjX = linept2->x - linept1->x;
return (float)(RADS_TO_DEGS * atan2( oppY , adjX ) );
}
//------------------------
//-- Returns the angle, between 0 and 360 degrees, of a line from the horizontal.
float line_getAngle2DPos ( Ipoint *linept1, Ipoint *linept2 )
{
return (float)fMod( (line_getAngle2D(linept1,linept2)+360.0f), 360.0f );
}
//------------------------
//-- Calculates the angle (in RADIANS) of a line from the horizontal
//-- (the angle formed when a ray projected left of the first point)
//-- NOTE: Returned angle will be between -PI and PI.
float line_getRadians2D ( Ipoint *linept1, Ipoint *linept2 )
{
float oppY = linept2->y - linept1->y;
float adjX = linept2->x - linept1->x;
return atan2( oppY , adjX );
}
//------------------------
//-- Calculates the angle theta (in RADIANS) between two connected lines defined by the points: {(x1,y1), (x2,y2)} and {(x2,y2), (x3,y3)}.
//-- (i.e. two lines connected with pt2 in the middle).
float line_radiansFormed3Pts( Ipoint *pt1, Ipoint *pt2, Ipoint *pt3 )
{
float line1Rads = line_getRadians2D ( pt1, pt2 );
float line2Rads = line_getRadians2D ( pt3, pt2 );
return fModWithinRange(line2Rads-line1Rads , -PI, PI );
}
//------------------------
//-- Returns the point halfway between the two given points.
Ipoint line_getPtHalfwayBetween(Ipoint *pt1, Ipoint *pt2)
{
Ipoint pt;
pt.x = avg( pt1->x, pt2->x );
pt.y = avg( pt1->y, pt2->y );
pt.z = avg( pt1->z, pt2->z );
return ( pt );
}
//------------------------
//-- Finds a point some fraction of the distance between pt1 and pt2 (ignoring Z)
//-- NOTE: Can also be used to return a point BEYOND the line
//-- (if the fractBetweenPts value is > 1 or < 0)
Ipoint line_findPtFractBetweenPts2D( const Ipoint *pt1, const Ipoint *pt2,
float fractBetweenPts )
{
float diffX = pt2->x - pt1->x;
float diffY = pt2->y - pt1->y;
Ipoint pt;
pt.x = (fractBetweenPts*diffX) + pt1->x;
pt.y = (fractBetweenPts*diffY) + pt1->y;
pt.z = pt1->z;
return pt;
}
//------------------------
//-- Finds a point some fraction of the distance between pt1 and pt2
//-- NOTE: Can also be used to return a point BEYOND the line
//-- (if the fractBetweenPts value is > 1 or < 0)
Ipoint line_findPtFractBetweenPts( const Ipoint *pt1, const Ipoint *pt2,
const float fractBetweenPts )
{
float diffX = pt2->x - pt1->x;
float diffY = pt2->y - pt1->y;
float diffZ = pt2->z - pt1->z;
Ipoint pt;
pt.x = (fractBetweenPts*diffX) + pt1->x;
pt.y = (fractBetweenPts*diffY) + pt1->y;
pt.z = (fractBetweenPts*diffZ) + pt1->z;
return pt;
}
//------------------------
//-- Calculates the distance SQUARED between two points in x and y
float line_sqDistBetweenPts2D( Ipoint *pt1, Ipoint *pt2 )
{
float diffX = pt2->x - pt1->x; // |-- calculate distance along x & y axis
float diffY = pt2->y - pt1->y; // |
return ( diffX*diffX + diffY*diffY ); // calculate distance squared
}
//------------------------
//-- Calculates the distance between two points in 2D space (z values are ignored)
float line_distBetweenPts2D( Ipoint *pt1, Ipoint *pt2 )
{
float diffX = pt2->x - pt1->x; // |-- calculate distance along x & y axis
float diffY = pt2->y - pt1->y; // |
return sqrt( diffX*diffX + diffY*diffY ); // calculate distance
}
//------------------------
//-- Returns true if the point lines on the line
//-- NOTE: This is done by calculating angles, which is NOT the most efficient way.
bool line_isPointOnLine( Ipoint *pt, Ipoint *lineStart, Ipoint *lineEnd )
{
return ( isBetween( lineStart->x, pt->x, lineEnd->x )
&& isBetween( lineStart->y, pt->y, lineEnd->y )
&& ( line_getAngle2D(lineStart,lineEnd) == line_getAngle2D(lineStart,pt) ) );
}
//------------------------
//-- Calculates the equation of a line through the two points provided and returns
//-- true if successful or false if points have same x value (or line is vertical)
//-- y = gradient*x + offset
bool line_getLineEquation( Ipoint *pt1, Ipoint *pt2, float *gradient, float *offset )
{
float diffY = pt2->y - pt1->y;
float diffX = pt2->x - pt1->x;
if(diffX == 0)
{
*gradient = FLOAT_MAX;
*offset = pt1->y;
return false;
}
*gradient = diffY/diffX;
*offset = pt1->x*(*gradient) + pt1->y;
return true;
}
//------------------------
//-- Computes the cross product of the two vectors defined by points
//-- {(x1,y1), (x2,y2)} and {(x1,y1), (x3,y3)}.
//-- (i.e. two lines connected with pt1 in the middle).
//--
//-- NOTE: If the cross product is < 0 then
//-- pt3 o_ the line formed by {pt1-pt2},{pt2,pt3}
//-- -_ (with pt2 in the middle) forms a "left turn"
//-- -_
//-- pt2 o---------o pt1 < "middle"
float line_crossProduct3Points( Ipoint *pt1, Ipoint *pt2, Ipoint *pt3)
{
return ((pt2->x - pt1->x)*(pt3->y - pt1->y)) - ((pt2->y - pt1->y)*(pt3->x - pt1->x));
}
//------------------------
//-- Calculates the angle theta (in DEGREES) between two connected lines
//-- defined by the points: {(x1,y1), (x2,y2)} and {(x2,y2), (x3,y3)}.
//-- (i.e. two lines connected with pt2 in the middle).
//-- The result will be between 0 and 360.
//--
//-- pt3 o
//-- /
//-- / theta
//-- pt2 o---------o pt1
float line_angleFormed3Pts( Ipoint *pt1, Ipoint *pt2, Ipoint *pt3 )
{
float line1Ang = line_getAngle2D ( pt1, pt2 );
float line2Ang = line_getAngle2D ( pt3, pt2 );
return fMod( line2Ang-line1Ang, 360.0 );
}
//------------------------
//-- Returns the position of a point placed "distFromEnd" away from the "end" point
//-- of a line in the direction "angleFromStraight" degrees away from the direction
//-- the line is pointing.
//--
//-- eg: if distFromEnd=2 and angleFromStraight=90:
//--
//-- o <-- the point returned would be here
//--
//-- start --> o-----------o <-- end
Ipoint line_getPtRelativeToEnd( Ipoint *start, Ipoint *end,
float distFromEnd, float angleFromStraight )
{
float angleFromHorz = line_getAngle2D( start, end );
float theta = (angleFromHorz + angleFromStraight) * DEGS_TO_RADS;
float offsetX = distFromEnd*cos(theta);
float offsetY = distFromEnd*sin(theta);
Ipoint returnPt;
returnPt.x = offsetX + end->x;
returnPt.y = offsetY + end->y;
returnPt.z = end->z;
return returnPt;
}
//------------------------
//-- Takes two lines (four coordinates) and returns true, plus point of
//-- intersection if the lines cross each other - otherwise returns false.
//-- NOTE: doubles are used here to increase precision of final value.
bool line_doLinesCrossAndWhere( Ipoint *line1pt1, Ipoint *line1pt2,
Ipoint *line2pt1, Ipoint *line2pt2, Ipoint *intercept )
{
if( !imodPointIntersect( line1pt1, line1pt2, line2pt1, line2pt2 ) )
return false;
double a1 = line1pt2->y-line1pt1->y; // rise (difference Y)
double b1 = line1pt1->x-line1pt2->x; // run (difference X)
double c1 = line1pt2->x*line1pt1->y - line1pt1->x*line1pt2->y;
//{ a1*x + b1*y + c1 = 0 is line 1 }
double a2 = line2pt2->y-line2pt1->y;
double b2 = line2pt1->x-line2pt2->x;
double c2 = line2pt2->x*line2pt1->y - line2pt1->x*line2pt2->y;
//{ a2*x + b2*y + c2 = 0 is line 2 }
double denom = (a1*b2) - (a2*b1);
if (denom == 0) // lines are parallel
return false;
intercept->x = (b1*c2 - b2*c1) / denom;
intercept->y = (a2*c1 - a1*c2) / denom;
intercept->z = line1pt1->z;
return true;
}
//------------------------
//-- Takes two ray (defined by two points) and determines the point they intersect.
//-- NOTE: If the rays are parallen (thus don't intersect) returns false.
bool line_getInterceptWhereRayCross( Ipoint *line1pt1, Ipoint *line1pt2,
Ipoint *line2pt1, Ipoint *line2pt2,
Ipoint *intercept )
{
double a1 = line1pt2->y-line1pt1->y; // rise (difference Y)
double b1 = line1pt1->x-line1pt2->x; // run (difference X)
double c1 = line1pt2->x*line1pt1->y - line1pt1->x*line1pt2->y;
//{ a1*x + b1*y + c1 = 0 is line 1 }
double a2 = line2pt2->y-line2pt1->y;
double b2 = line2pt1->x-line2pt2->x;
double c2 = line2pt2->x*line2pt1->y - line2pt1->x*line2pt2->y;
//{ a2*x + b2*y + c2 = 0 is line 2 }
double denom = (a1*b2) - (a2*b1);
if (denom == 0) // lines are parallel
return false;
intercept->x = (b1*c2 - b2*c1) / denom;
intercept->y = (a2*c1 - a1*c2) / denom;
intercept->z = line1pt1->z;
return true;
}
//------------------------
//-- Takes two ray (each defined by two points) and a given point between them and calculates how close
//-- the point is to the first ray, compared to the second ray.
//--
//-- ===================================================
//-- X intercept
//--
//-- o
//-- | o
//-- | \
//-- ray1 | \ ray2
//-- o X \
//-- pt \
//-- o
//--
//-- In this example ray1-intercept-pt forms 15 degrees
//-- and ray1-intercept-ray2 forms 30 degrees, so value
//-- returned would be 0.5
//-- ===================================================
//--
//-- o o
//-- \ \
//-- \ pt \
//-- ray1X-- X X X --- ray2X
//-- \ \
//-- o o
//--
//-- In this example the rays are parallel so it's calculated differently.
//-- The pt is 1/4 the distance between the points ray1X and ray2X
//-- so 0.25 is returned.
//--
double line_getFractPtBetweenTwoRays( Ipoint *ray1pt1, Ipoint *ray1pt2,
Ipoint *ray2pt1, Ipoint *ray2pt2, Ipoint *pt )
{
double a1 = ray1pt2->y-ray1pt1->y;
double b1 = ray1pt1->x-ray1pt2->x;
double c1 = ray1pt2->x*ray1pt1->y - ray1pt1->x*ray1pt2->y; //{ a1*x + b1*y + c1 = 0 is line 1 }
double a2 = ray2pt2->y-ray2pt1->y;
double b2 = ray2pt1->x-ray2pt2->x;
double c2 = ray2pt2->x*ray2pt1->y - ray2pt1->x*ray2pt2->y; //{ a2*x + b2*y + c2 = 0 is line 2 }
double denom = a1*b2 - a2*b1;
if (denom == 0) // if rays are parallel.
{
if( a1 == 0 || a2 == 0 ) // if both lines are horizontal (no x-intercept) calculate y intercept.
{
double ray1Y = -(a1*pt->x + c1)/b1;
double ray2Y = -(a2*pt->x + c2)/b2;
//cout << " HORIZONTAL! ray1Y=" << ray1Y << " ray2Y=" << ray2Y << endl;
return fDiv( pt->y - ray1Y, ray2Y - ray1Y );
}
else
{
double ray1X = -(b1*pt->y + c1)/a1;
double ray2X = -(b2*pt->y + c2)/a2;
//cout << " PARALLEL! ray1X=" << ray1X << " ray2X=" << ray2X << endl;
return fDiv( pt->x - ray1X , ray2X - ray1X );
}
}
else
{
//## CALCULATE INTERCEPT POINT:
Ipoint intercept;
intercept.x = (b1*c2 - b2*c1) / denom;
intercept.y = (a2*c1 - a1*c2) / denom;
intercept.z = ray1pt1->z;
double angleBetweenRays = line_radiansFormed3Pts( ray1pt1, &intercept, ray2pt1 );
double angleRayOneAndInterceptToPoint = line_radiansFormed3Pts( ray1pt1, &intercept, pt );
//cout << " angleBetweenRays=" << angleBetweenRays << " angleRayOneAndInterceptToPoint=" << angleRayOneAndInterceptToPoint << endl; //%%%%%%
return fDiv( angleRayOneAndInterceptToPoint, angleBetweenRays );
}
}
//------------------------
//-- Takes two lines (each defined by two points) and a given point between them and calculates how close
//-- far along the lines the point lies, imagining the lines are joined to form a quadrilatral.
//-- NOTE: If you are sure the point is inside the quadrilateral you can set checkOutsideNegative to false.
//--
//-- ===================================================
//-- l1pt1 l1pt2
//-- o--------o
//-- ` `
//-- X pt `
//-- ` X `
//-- o__ X
//-- l2pt1 --__ `
//-- --__ `
//-- --__ `
//-- --o l2pt2
//--
//-- In this example the fractBetweenLines = 0.5
//-- ===================================================
//--
//-- o--------------o
//-- ` `
//-- ` pt `
//-- X X X
//-- ` `
//-- o--------------o
//--
//-- In this example the fractBetweenLines = 0.9
//--
//-- ===================================================
//--
//-- o--------------o
//-- ` `
//-- ` `
//-- X X X
//-- ` `
//-- o--------------o
//--
//-- In this example the fractBetweenLines = -0.3
//-- (but only if checkOutsideNegative==true)
//--
double line_getFractPtBetweenTwoLines( Ipoint *l1p1, Ipoint *l1p2,
Ipoint *l2p1, Ipoint *l2p2,
Ipoint *pt, bool checkOutsideNegative )
{
//## CALCULATE FRACTRIGHT:
double a1 = l1p2->y-l1p1->y;
double b1 = l1p1->x-l1p2->x;
double c1 = l1p2->x*l1p1->y - l1p1->x*l1p2->y; //{ a1*x + b1*y + c1 = 0 is line 1 }
double a2 = l2p2->y-l2p1->y;
double b2 = l2p1->x-l2p2->x;
double c2 = l2p2->x*l2p1->y - l2p1->x*l2p2->y; //{ a2*x + b2*y + c2 = 0 is line 2 }
double denom = a1*b2 - a2*b1;
//## CALCULATE INTERCEPT POINT OF TOP AND BOTTOM LINES:
Ipoint lineIntercept; // where the two lines intercept each other
if (denom == 0) // if lines are parallel: make "lineIntercept" same offset/angle as first line
{
lineIntercept.x = pt->x + (l1p2->x - l1p1->x);
lineIntercept.y = pt->y + (l1p2->y - l1p1->y);
lineIntercept.z = l1p1->z;
}
else // if lines not parallel: calculate "lineIntercept" where they cross
{
lineIntercept.x = (b1*c2 - b2*c1) / denom;
lineIntercept.y = (a2*c1 - a1*c2) / denom;
lineIntercept.z = l1p1->z;
}
Ipoint incptBetweenPt1s;
line_getInterceptWhereRayCross( &lineIntercept, pt, l1p1, l2p1, &incptBetweenPt1s );
Ipoint incptBetweenPt2s;
line_getInterceptWhereRayCross( &lineIntercept, pt, l1p2, l2p2, &incptBetweenPt2s );
double fractBetweenLines = fDiv( line_distBetweenPts2D( &incptBetweenPt1s, pt ), line_distBetweenPts2D( &incptBetweenPt1s, &incptBetweenPt2s ) );
if (checkOutsideNegative)
if( ABS( line_radiansFormed3Pts(pt,&incptBetweenPt1s,&incptBetweenPt2s) ) > 3 ) // if point occurs behind the line connecting l1p1 and l2p1: make fractBetweenLines negative // if (incptBetweenPt1s < incptBetweenPt2s && pt < incptBetweenPt1s || incptBetweenPt1s > incptBetweenPt2s && pt > incptBetweenPt1s)
fractBetweenLines = -fractBetweenLines;
return (fractBetweenLines);
}
//------------------------
//-- Takes a given point inside quadrilateral 1 (defined by four points) and returns
//-- the matching point inside quadrilateral 2.
//--
//-- ===================================================
//--
//-- quad1 quad2
//-- o---X---o o----X----o
//-- | \ \ \
//-- X pt \ \ ptR \
//-- | X \ X X X
//-- o__ X \ \
//-- --__ \ \ \
//-- -X__ \ o----X----o
//-- --__\
//-- -o
//--
//-- In this example the pt is fractionUp = 0.5 & fractionRight = 0.5
//-- inside quad1 and cooresponding point in quad2 is shown as (ptR).
//--
Ipoint point_findPtInQuad1InMatchingQuad2( Ipoint *pt,
Ipoint *q1BL, Ipoint *q1TL,
Ipoint *q1TR, Ipoint *q1BR,
Ipoint *q2BL, Ipoint *q2TL,
Ipoint *q2TR, Ipoint *q2BR )
{
double q1FractUp = line_getFractPtBetweenTwoLines( q1BL, q1TL, q1BR, q1TR, pt, true );
double q1FractRight = line_getFractPtBetweenTwoLines( q1BL, q1BR, q1TL, q1TR, pt, true );
Ipoint q2PartwayAlongBottom = line_findPtFractBetweenPts2D( q2BL, q2BR, q1FractRight );
Ipoint q2PartwayAlongTop = line_findPtFractBetweenPts2D( q2TL, q2TR, q1FractRight );
Ipoint ptR = line_findPtFractBetweenPts2D( &q2PartwayAlongBottom, &q2PartwayAlongTop, q1FractUp );
return (ptR);
}
//------------------------
//-- Takes five points which form two 3-point lines:
//-- {p1-pMid-p2} AND {p3-pMid-p4}
//-- Returns true if the first line appears to "kiss" the second -
//-- if it touches the other line but does NOT pass through it (example below).
//-- Returns false if the first line goes through the other.
//--
//--
//-- p3 p4
//-- o o
//-- | /
//-- |/
//-- o======o======o <-- this is a "kiss" because the line from p1-pMid-p2
//-- p2 pMid p1 does not pass through the other line.
bool line_isKiss( Ipoint *pMid, Ipoint *p1,Ipoint *p2, Ipoint *p3, Ipoint *p4 )
{
float ang1 = line_getAngle2D( pMid, p1 ); // angle each line relative to EAST
float ang2 = line_getAngle2D( pMid, p2 );
float ang3 = line_getAngle2D( pMid, p3 );
float ang4 = line_getAngle2D( pMid, p4 );
float relAng1 = 0;
float relAng2 = fMod( ang2 - ang1, 360.0 ); //|-- calculates angle of each
float relAng3 = fMod( ang3 - ang1, 360.0 ); //| line relative to the
float relAng4 = fMod( ang4 - ang1, 360.0 ); //| first line
return ( (relAng3 < relAng2 && relAng4 < relAng2 )
|| (relAng3 > relAng2 && relAng4 > relAng2 ) ); // determines if lines kiss.
}
//------------------------
//-- Takes six points which form two 3-point lines: {a1-a2-a3} AND {b1-b2-b3}
//-- Returns true if the first line appears to "kiss" the second
bool line_twoKiss( Ipoint *a1, Ipoint *a2, Ipoint *a3,
Ipoint *b1, Ipoint *b2, Ipoint *b3 )
{
if (a2 != b2)
return false;
return line_isKiss( a2, a1,a3, b1,b3 );
}
//############################################################
//----------------------------------------------------------------------------
//## EXTRA CONTOUR FUNCTIONS:
//----------------------------------------------------------------------------
//------------------------
//--
//-- __0___ ___1__ a _______b_______c
//-- |\_ | _/| |\_ | _/|
//-- | \12 10 13 | | \_ | _/ |
//-- 7| \_ | _/ |2 | \_ | _/ |
//-- | __8__\|/__9__ | h| _____\i/_____ |d
//-- | _/|\_ | | _/|\_ |
//-- 6| _/ | \_ |3 | _/ | \_ |
//-- | _/15 11 14_ | | _/ | \_ |
//-- |/______|______\| |/______|______\|
//-- 5 4 g f e
void cont_gen16SegDisplay( Iobj *obj,int s0,int s1,int s2,int s3,int s4,int s5,int s6,
int s7,int s8,int s9,int s10,int s11,int s12,int s13,int s14,int s15, int z )
{
Ipoint a; setPt(&a, 0, 12, z);
Ipoint b; setPt(&b, 3, 12, z);
Ipoint c; setPt(&c, 6, 12, z);
Ipoint d; setPt(&d, 6, 6, z);
Ipoint e; setPt(&e, 6, 0, z);
Ipoint f; setPt(&f, 3, 0, z);
Ipoint g; setPt(&g, 0, 0, z);
Ipoint h; setPt(&h, 0, 6, z);
Ipoint i; setPt(&i, 3, 6, z);
if( s0 ) cont_addTwoPointContourToObj(obj, a, b); // seg 0
if( s1 ) cont_addTwoPointContourToObj(obj, b, c); // seg 1
if( s2 ) cont_addTwoPointContourToObj(obj, c, d); // seg 2
if( s3 ) cont_addTwoPointContourToObj(obj, d, e); // seg 3
if( s4 ) cont_addTwoPointContourToObj(obj, e, f); // seg 4
if( s5 ) cont_addTwoPointContourToObj(obj, f, g); // seg 5
if( s6 ) cont_addTwoPointContourToObj(obj, g, h); // seg 6
if( s7 ) cont_addTwoPointContourToObj(obj, h, a); // seg 7
if( s8 ) cont_addTwoPointContourToObj(obj, h, i); // seg 8
if( s9 ) cont_addTwoPointContourToObj(obj, i, d); // seg 9
if( s10) cont_addTwoPointContourToObj(obj, b, i); // seg 12
if( s11) cont_addTwoPointContourToObj(obj, f, i); // seg 13
if( s12) cont_addTwoPointContourToObj(obj, a, i); // seg 10
if( s13) cont_addTwoPointContourToObj(obj, c, i); // seg 11
if( s14) cont_addTwoPointContourToObj(obj, e, i); // seg 12
if( s15) cont_addTwoPointContourToObj(obj, g, i); // seg 13
}
int cont_generateDigitUsing16SegDisplay( Iobj *obj, char ch, int z )
{
//int segDisp[14];
switch (ch)
{ // seg: 0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5
case('0'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,1,1,0,0,0,0,0,1,0,1,z); break;
case('1'): cont_gen16SegDisplay(obj,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,z); break;
case('2'): cont_gen16SegDisplay(obj,1,1,1,0,1,1,1,0,1,1,0,0,0,0,0,0,z); break;
case('3'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,0,0,1,1,0,0,0,0,0,0,z); break;
case('4'): cont_gen16SegDisplay(obj,0,0,1,1,0,0,0,1,1,1,0,0,0,0,0,0,z); break;
case('5'): cont_gen16SegDisplay(obj,1,1,0,1,1,1,0,1,1,1,0,0,0,0,0,0,z); break;
case('6'): cont_gen16SegDisplay(obj,1,1,0,1,1,1,1,1,1,1,0,0,0,0,0,0,z); break;
case('7'): cont_gen16SegDisplay(obj,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,z); break;
case('8'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,z); break;
case('9'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,0,1,1,1,0,0,0,0,0,0,z); break;
// seg: 0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5
case('a'): case('A'): cont_gen16SegDisplay(obj,1,1,1,1,0,0,1,1,1,1,0,0,0,0,0,0,z); break;
case('b'): case('B'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,0,0,0,1,1,1,0,0,0,0,z); break;
case('c'): case('C'): cont_gen16SegDisplay(obj,1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,z); break;
case('d'): case('D'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,0,0,0,0,1,1,0,0,0,0,z); break;
case('e'): case('E'): cont_gen16SegDisplay(obj,1,1,0,0,1,1,1,1,1,1,0,0,0,0,0,0,z); break;
case('f'): case('F'): cont_gen16SegDisplay(obj,1,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,z); break;
case('g'): case('G'): cont_gen16SegDisplay(obj,1,1,0,1,1,1,1,1,0,1,0,0,0,0,0,0,z); break;
case('h'): case('H'): cont_gen16SegDisplay(obj,0,0,1,1,0,0,1,1,1,1,0,0,0,0,0,0,z); break;
case('i'): case('I'): cont_gen16SegDisplay(obj,1,1,0,0,1,1,0,0,0,0,1,1,0,0,0,0,z); break;
case('j'): case('J'): cont_gen16SegDisplay(obj,0,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,z); break;
case('k'): case('K'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,1,1,1,0,0,0,0,1,1,0,z); break;
case('l'): case('L'): cont_gen16SegDisplay(obj,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,z); break;
case('m'): case('M'): cont_gen16SegDisplay(obj,0,0,1,1,0,0,1,1,0,0,0,0,1,1,0,0,z); break;
case('n'): case('N'): cont_gen16SegDisplay(obj,0,0,1,1,0,0,1,1,0,0,0,0,1,0,1,0,z); break;
case('o'): case('O'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,z); break;
case('p'): case('P'): cont_gen16SegDisplay(obj,1,1,1,0,0,0,1,1,1,1,0,0,0,0,0,0,z); break;
case('q'): case('Q'): cont_gen16SegDisplay(obj,1,1,1,1,1,1,1,1,0,0,0,0,0,0,1,0,z); break;
case('r'): case('R'): cont_gen16SegDisplay(obj,1,1,1,0,0,0,1,1,1,1,0,0,0,0,1,0,z); break;
case('s'): case('S'): cont_gen16SegDisplay(obj,1,1,0,1,1,1,0,1,1,1,0,0,0,0,0,0,z); break;
case('t'): case('T'): cont_gen16SegDisplay(obj,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,0,z); break;
case('u'): case('U'): cont_gen16SegDisplay(obj,0,0,1,1,1,1,1,1,0,0,0,0,0,0,0,0,z); break;
case('v'): case('V'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,1,z); break;
case('w'): case('W'): cont_gen16SegDisplay(obj,0,0,1,1,0,0,1,1,0,0,0,0,0,0,1,1,z); break;
case('x'): case('X'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,z); break;
case('y'): case('Y'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,z); break;
case('z'): case('Z'): cont_gen16SegDisplay(obj,1,1,0,0,1,1,0,0,0,0,0,0,0,1,0,1,z); break;
case('.'): cont_addTwoPointContourToObj(obj, 2,0, 4,0, z);
cont_addTwoPointContourToObj(obj, 4,0, 4,2, z);
cont_addTwoPointContourToObj(obj, 4,2, 2,2, z);
cont_addTwoPointContourToObj(obj, 2,2, 2,0, z);
break;
case(':'): cont_addTwoPointContourToObj(obj, 2,2, 4,2, z);
cont_addTwoPointContourToObj(obj, 4,2, 4,4, z);
cont_addTwoPointContourToObj(obj, 4,4, 2,4, z);
cont_addTwoPointContourToObj(obj, 2,4, 2,2, z);
cont_addTwoPointContourToObj(obj, 2,8, 4,8, z);
cont_addTwoPointContourToObj(obj, 4,8, 4,10, z);
cont_addTwoPointContourToObj(obj, 4,10, 2,10, z);
cont_addTwoPointContourToObj(obj, 2,10, 2,8, z);
break;
case(','): cont_addTwoPointContourToObj(obj, 0,0, -0.5,-0.5, z);
break;
case('='): cont_addTwoPointContourToObj(obj, 0,7.5, 6,7.5, z);
cont_addTwoPointContourToObj(obj, 0,4.5, 6,4.5, z);
break;
// seg: 0,1,2,3,4,5,6,7,8,9,0,1,2,3,4,5
case('-'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,z); break;
case('_'): cont_gen16SegDisplay(obj,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,z); break;
case('['): cont_gen16SegDisplay(obj,0,1,0,0,1,0,0,0,0,0,1,1,0,0,0,0,z); break;
case(']'): cont_gen16SegDisplay(obj,1,0,0,0,0,1,0,0,0,0,1,1,0,0,0,0,z); break;
case('('): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,z); break;
case(')'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,z); break;
case('*'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,z); break;
case('/'): cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,z); break;
case('\\'):cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,z); break;
case('\''):cont_gen16SegDisplay(obj,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,z); break;
case('"'): cont_gen16SegDisplay(obj,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,z); break;
//case(''): cont_gen16SegDisplay(obj,,,,,,,,,,,,,,,,,z); break;
default: return 0;
}
return 1;
}
//------------------------
//--
int cont_generateDigitUsing7SegDisplay( Iobj *obj, int number, int z )
{
if( number < 0 || number > 9 )
return 0;
static int segDisp[10][7] = {
// seg: A,B,C,D,E,F,G
{1,1,1,1,1,1,0}, // "0"
{0,1,1,0,0,0,0}, // "1"
{1,1,0,1,1,0,1}, // "2"
{1,1,1,1,0,0,1}, // "3"
{0,1,1,0,0,1,1}, // "4"
{1,0,1,1,0,1,1}, // "5"
{1,0,1,1,1,1,1}, // "6"
{1,1,1,0,0,0,0}, // "7"
{1,1,1,1,1,1,1}, // "8"
{1,1,1,1,0,1,1} // "9"
};
if( segDisp[number][0] ) cont_addTwoPointContourToObj(obj,0,2, 1,2, z); // seg A
if( segDisp[number][1] ) cont_addTwoPointContourToObj(obj,1,2, 1,1, z); // seg B
if( segDisp[number][2] ) cont_addTwoPointContourToObj(obj,1,1, 1,0, z); // seg C
if( segDisp[number][3] ) cont_addTwoPointContourToObj(obj,1,0, 0,0, z); // seg D
if( segDisp[number][4] ) cont_addTwoPointContourToObj(obj,0,0, 0,1, z); // seg E
if( segDisp[number][5] ) cont_addTwoPointContourToObj(obj,0,1, 0,2, z); // seg F
if( segDisp[number][6] ) cont_addTwoPointContourToObj(obj,0,1, 1,1, z); // seg G
return 1;
}
//------------------------
//--
//--
int cont_generateTextAsConts( Iobj *obj, string text, Ipoint pos,
float fontSize, int textAlign, bool smallCaps )
{
if( text.length() <= 0 )
return 0;
Ipoint originPt; setPt( &originPt, 0, 0, 0 );
float scaleFactor = fontSize / DEFAULT_FONT_HEIGHT;
float finalWidth = text.length() * DEFAULT_FONT_SPACING * scaleFactor;
float translateX = pos.x;
if(textAlign == TA_CENTER)
translateX -= finalWidth/2.0;
if(textAlign == TA_RIGHT )
translateX -= finalWidth;
//## GENERATE EACH CHARACTER:
int nContsBeforeText = csize(obj);
for(int i=0; i<text.length(); i++ )
{
int nContsBeforeChar = csize(obj);
cont_generateDigitUsing16SegDisplay( obj, text[i], pos.z );
bool isLowerCase = ( int(text[i]) >= int('a') ) && ( int(text[i]) <= int('z') );
for (int c=nContsBeforeChar; c<csize(obj); c++)
{
if( smallCaps && isLowerCase )
cont_scaleAboutPtXY( getCont(obj,c), &originPt, 0.9, 0.7 );
cont_translate( getCont(obj,c), DEFAULT_FONT_SPACING*i, 0 );
}
}
//## FINAL SCALE AND TRANSLATE OF TEXT
for (int c=nContsBeforeText; c<csize(obj); c++)
{
cont_scaleAboutPtXY( getCont(obj,c), &originPt, scaleFactor, scaleFactor );
cont_translate( getCont(obj,c), translateX, pos.y );
}
return 1;
}
//------------------------
//--
//--
int cont_generateTextAreaAsConts( Iobj *obj, string text, Ipoint pos, float fontSize,
int alignHoriz, int alignVert, bool smallCaps, float lineSpacing )
{
if (text.length() == 0 )
return 0;
vector<string> lines = string_explode( text, "\n" );
int numLines = (int)lines.size();
float scaleFactor = fontSize / DEFAULT_FONT_HEIGHT;
int maxCharsWide = 0;
for (int i=0; i<(int)lines.size(); i++)
maxCharsWide = MAX( maxCharsWide, (int)lines[i].length() );
float textWidth = maxCharsWide * DEFAULT_FONT_SPACING * scaleFactor;
float textHeight = (numLines * fontSize) + ((numLines-1) * lineSpacing);
float yPosFirstLine = pos.y;
if( alignVert == TV_CENTER )
yPosFirstLine -= fontSize*0.5 - ( (numLines-1) * 0.5*(fontSize+lineSpacing) );
if( alignVert == TV_TOP )
yPosFirstLine += fontSize/2.0;
if( alignVert == TV_BOTTOM )
yPosFirstLine -= (textHeight - fontSize);
for (int i=0; i<(int)lines.size(); i++)
{
Ipoint linePos; setPt(&linePos, pos.x, yPosFirstLine, pos.z );
float lineWidth = lines[i].length() * DEFAULT_FONT_SPACING * scaleFactor;
if( alignHoriz == TA_CENTER )
linePos.x -= lineWidth/2.0;
if( alignHoriz == TA_RIGHT )
linePos.x -= lineWidth;
linePos.y = yPosFirstLine - i*(fontSize+lineSpacing);
cont_generateTextAsConts( obj, lines[i], linePos, fontSize, TA_LEFT, smallCaps );
}
return numLines;
}
//------------------------
//-- Adds a new contour to the specified object
int cont_addTwoPointContourToObj( Iobj *obj, float p1x,float p1y,
float p2x,float p2y, float z, int open )
{
Icont *newCont = imodContourNew(); // malloc new contour and don't delele it
imodPointAppendXYZ( newCont, p1x, p1y, z );
imodPointAppendXYZ( newCont, p2x, p2y, z );
setOpenFlag( newCont, open );
int numConts = csize(obj);
int newContPos = imodObjectAddContour( obj, newCont );
free(newCont);
return newContPos;
}
//------------------------
//-- Adds a new contour to the specified object
int cont_addTwoPointContourToObj( Iobj *obj, Ipoint p1, Ipoint p2, int open )
{
Icont *newCont = imodContourNew(); // malloc new contour and don't delele it
imodPointAppend( newCont, &p1 );
imodPointAppend( newCont, &p2 );
int numConts = csize(obj);
setOpenFlag( newCont, open );
int newContPos = imodObjectAddContour( obj, newCont );
free(newCont);
return newContPos;
}
//------------------------
//-- Returns 1 if the two contours have an identical set of points
//-- (i.e. same number of points in same position)
int cont_isEqual( Icont *cont1, Icont *cont2 )
{
if( psize(cont1) != psize(cont2) )
return false;
for (int i=0; i<psize(cont1) && i<<psize(cont2); i++)
if( !imodPointIsEqual( getPt(cont1,i), getPt(cont2,i) ) )
return false;
return true;
}
//------------------------
//-- Returns 1 if the contour has a point with the same x and y
//-- value as the given point
int cont_doesPtExistInCont( Icont *cont, Ipoint *pt )
{
for (int i=0; i<psize(cont); i++)
if( ptsEqual(pt, getPt(cont,i)) )
return 1;
return 0;
}
//------------------------
//-- Returns the lower left and upper right corners of the minimum bounding rectangle
//-- around a contour.
bool cont_getMBR( Icont *cont, Ipoint *ll, Ipoint *ur )
{
if( isEmpty(cont) || psize(cont) == 0 )
{
setPt( ll, 0, 0, 0 );
setPt( ur, 0, 0, 0 );
return false;
}
else if( psize(cont) == 1 )
{
setPt( ll, getFirstPt(cont)->x, getFirstPt(cont)->y, getFirstPt(cont)->z );
setPt( ur, getFirstPt(cont)->x, getFirstPt(cont)->y, getFirstPt(cont)->z );
}
setPt( ll, FLOAT_MAX, FLOAT_MAX, FLOAT_MAX );
setPt( ur, FLOAT_MIN, FLOAT_MIN, FLOAT_MIN );
for( int p=0; p<psize(cont); p++ )
{
Ipoint *pt = getPt(cont, p);
if( ll->x > pt->x ) ll->x = pt->x;
if( ll->y > pt->y ) ll->y = pt->y;
if( ll->z > pt->z ) ll->z = pt->z;
if( ur->x < pt->x ) ur->x = pt->x;
if( ur->y < pt->y ) ur->y = pt->y;
if( ur->z < pt->z ) ur->z = pt->z;
}
return true;
}
//------------------------
//-- Returns the center of the minimum bounding rectangle around a contour.
//-- NOTE: Is a bit less expensive than "cont_getCentroid"
bool cont_getCenterOfMBR( Icont *cont, Ipoint *rpt )
{
if( isEmpty(cont) || psize(cont) == 0 )
{
setPt( rpt, 0, 0, 0 );
return false;
}
else if( psize(cont) == 1 )
{
setPt( rpt, getPt(cont,0)->x, getPt(cont,0)->y, getPt(cont,0)->z );
return true;
}
Ipoint minPt;
Ipoint maxPt;
setPt( &minPt, FLOAT_MAX, FLOAT_MAX, FLOAT_MAX );
setPt( &maxPt, FLOAT_MIN, FLOAT_MIN, FLOAT_MIN );
for( int p=0; p<psize(cont); p++ )
{
Ipoint *pt = getPt(cont, p);
if( minPt.x > pt->x ) minPt.x = pt->x;
if( minPt.y > pt->y ) minPt.y = pt->y;
if( minPt.z > pt->z ) minPt.z = pt->z;
if( maxPt.x < pt->x ) maxPt.x = pt->x;
if( maxPt.y < pt->y ) maxPt.y = pt->y;
if( maxPt.z < pt->z ) maxPt.z = pt->z;
}
setPt( rpt, avg(minPt.x, maxPt.x), avg(minPt.y, maxPt.y), avg(minPt.z, maxPt.z) );
return true;
}
//------------------------
//-- Returns the centroid/center of mass of a contour AND the area.
//-- WARNING: Contour must be enclosed (must not cross it's own path).
//-- WARNING: If contour is drawn clockwise, area will be a negative value.
bool cont_getCentroid( Icont *cont, Ipoint *rpt )
{
if( isEmpty(cont) || psize(cont) == 0 )
{
setPt( rpt, 0, 0, 0 );
return false;
}
else if( psize(cont) == 1 )
{
setPt( rpt, getPt(cont,0)->x, getPt(cont,0)->y, getPt(cont,0)->z );
return true;
}
float totArea = 0;
float xCentroid = 0;
float yCentroid = 0;
for(int i=0; i<psize(cont); i++) {
float boundingRecOverLine = (getPt(cont,i)->x * getPt(cont,i+1)->y)
- (getPt(cont,i+1)->x * getPt(cont,i)->y);
totArea += boundingRecOverLine;
xCentroid += (( getPt(cont,i)->x + getPt(cont,i+1)->x ) * boundingRecOverLine);
yCentroid += (( getPt(cont,i)->y + getPt(cont,i+1)->y ) * boundingRecOverLine);
}
totArea = 0.5*totArea;
xCentroid = 1.0/(6.0*totArea) * xCentroid;
yCentroid = 1.0/(6.0*totArea) * yCentroid;
setPt( rpt, xCentroid, yCentroid, getPt(cont,0)->z );
return true;
}
//------------------------
//-- Returns true if the "middle-most" point (in sequence) in "cont1" is inside "cont2".
//-- WARNING: If the contour only has only one or two points, the first is used.
bool cont_insideCrude( Icont *cont1, Icont *cont2 )
{
if( psize(cont1)==0 )
return false;
int middlePtIdx = int( psize(cont1) / 2);
return imodPointInsideCont( cont2, getPt( cont1, middlePtIdx ) );
}
//------------------------
//-- Returns the radius which would give a circle of equal area to
//-- the contour provided.
//-- WARNING: Contour must be enclosed (must not cross it's own path).
float cont_getRadius( Icont *c )
{
float area = imodContourArea( c );
return sqrt(area/PI);
}
//------------------------
//-- Finds the closest point in the given contour to the given point
//-- and returns the distance, plus the index value and the point itself.
void cont_findClosestPtInContToGivenPt( Ipoint *pt, Icont *cont, float *closestDist,
Ipoint *closestPt, int *closestPtIdx )
{
float closestDist2 = FLOAT_MAX;
for(int i=0; i<psize( cont ); i++ )
{
float sqDistToPt = line_sqDistBetweenPts2D( getPt(cont,i) , pt );
if( closestDist2 > sqDistToPt ) {
closestDist2 = sqDistToPt;
*closestPtIdx = i;
}
}
*closestPt = *getPt(cont,*closestPtIdx);
*closestDist = sqrt(closestDist2);
}
//------------------------
//-- Find the minimum, average and maximum distance in 3D from the specified
//-- point to the points in the contour.
void cont_findMinMaxAndAvgDistFromPt( Ipoint *pt, Icont *cont, float zScale,
float &minDist, float &maxDist, float &avgDist )
{
minDist = FLOAT_MAX;
maxDist = 0;
avgDist = 0;
Ipoint scalePt;
setPt( &scalePt, 1,1,zScale );
for(int p=0; p<psize(cont); p++ )
{
float distToPt = imodPoint3DScaleDistance( getPt(cont,p), pt, &scalePt );
minDist = MIN( minDist, distToPt );
maxDist = MAX( maxDist, distToPt );
avgDist += distToPt;
}
if( psize(cont) > 0 )
avgDist = avgDist / psize(cont);
}
//------------------------
//-- Checks if two contours touch each other (or if one is contained in the other)
//-- First checks if minimum bounding rectangles around each overlap,
//-- then if contours cross, and finally if one is inside the other
bool cont_doContsTouch( Icont *cont1, Icont *cont2 )
{
if ( isEmpty(cont1) || isEmpty(cont2) )
return false;
Ipoint cont1ll, cont1ur, cont2ll, cont2ur;
Icont *cs1 = imodel_contour_scan( cont1 );
Icont *cs2 = imodel_contour_scan( cont2 );
imodContourGetBBox( cont1, &cont1ll, &cont1ur );
imodContourGetBBox( cont2, &cont2ll, &cont2ur );
bool result = imodel_scans_overlap( cs1,cont1ll,cont1ur, cs2,cont2ll,cont2ur);
imodContourDelete( cs1 );
imodContourDelete( cs2 );
return result;
}
//------------------------
//-- Finds the closest distance between the given point and any point
//-- in the contour.
//-- Note that if the point is INSIDE the contour and returnZeroIfPointIsInside
//-- is true the function will return 0.
float cont_minDistPtAndContourPts2D(Ipoint *pt, Icont *cont, bool returnZeroIfPtInside)
{
if( returnZeroIfPtInside && imodPointInsideCont( cont, pt ) )
return (0);
float closestDist2 = FLOAT_MAX;
for(int i=0; i<psize( cont ); i++ )
updateMin ( closestDist2, line_sqDistBetweenPts2D( getPt(cont,i), pt ) );
return sqrt(closestDist2);
}
//------------------------
//-- Finds the closest distance between any two points (ignoring Z axis)...
//-- but note that this is NOT necessarily the closest distance between
//-- the lines formed by both contours.
float cont_minDistBetweenContPts2D( Icont *cont1, Icont *cont2, bool returnZeroIfTouch )
{
if( returnZeroIfTouch && cont_doContsTouch( cont1, cont2) )
return 0;
float closestDist2 = FLOAT_MAX;
for(int i=0; i<psize(cont1); i++ )
for(int j=0; j<psize(cont2); j++ )
updateMin(closestDist2,line_sqDistBetweenPts2D(getPt(cont1,i),getPt(cont2,j)));
return sqrt(closestDist2);
}
//------------------------
//-- Returns the closest distance between any two points in 3D and also
//-- returns the corrdinates of those points in pt1 and pt2.
//-- but note that this is NOT necessarily the closest distance between
//-- the lines formed by both contours.
float cont_minDistBetweenContPts3D( Icont *cont1, Icont *cont2, float zScale,
Ipoint *pt1, Ipoint *pt2 )
{
Ipoint scalePt;
setPt( &scalePt, 1,1,zScale );
float closestDist2 = FLOAT_MAX;
for(int i=0; i<psize(cont1); i++ )
for(int j=0; j<psize(cont2); j++ )
{
float distX = getPt(cont1,i)->x - getPt(cont2,j)->x;
float distY = getPt(cont1,i)->y - getPt(cont2,j)->y;
float distZ = (getPt(cont1,i)->z - getPt(cont2,j)->z) * zScale;
float distXZ = sqrt( distX*distX + distY*distY );
float dist2 = (distXZ*distXZ + distZ*distZ);
if( dist2 < closestDist2 )
{
closestDist2 = dist2;
copyPt( pt1, getPt(cont1,i) );
copyPt( pt2, getPt(cont2,j) );
}
}
return sqrt(closestDist2);
}
//------------------------
//-- Reorder points so that the point at index "idxNewFirstPt" become the
//-- first point in the contour
void cont_reorderPtsToStartAtIdx( Icont *c, int idxNewFirstPt )
{
if( idxNewFirstPt==0 )
return;
int nPts = psize( c );
if( idxNewFirstPt < 0 || idxNewFirstPt > nPts ) { //%%%%%% SHOULD NEVER HAPPEN
wprint( "ERROR: cont_reorderPtsToStartAtIdx() - bad index" );
}
Icont *copy = imodContourDup( c );
deleteAllPts( c );
for (int i=idxNewFirstPt; i<(nPts + idxNewFirstPt); i++)
imodPointAppend( c, getPt( copy, i % nPts ) );
imodContourDelete( copy );
}
//------------------------
//-- Removes any points in cont that fall within the circle defined.
int cont_removePointsInCircle( Icont *cont, Ipoint *center,
float radius, bool checkZEachPoint )
{
int radiusSq = (radius*radius);
int numRemovedPoints = 0;
for(int p=0; p<psize(cont); p++ )
{
if( !checkZEachPoint || getPt(cont,p)->z == center->z )
{
float distSq = line_sqDistBetweenPts2D( center, getPt(cont,p) );
if( distSq < radiusSq )
{
imodPointDelete( cont, p );
p--;
numRemovedPoints++;
}
}
}
return (numRemovedPoints);
}
//------------------------
//-- Translates a contour by a given amout.
void cont_translate( Icont *cont, Ipoint *translate )
{
for(int i=0; i<psize( cont ); i++)
{
getPt(cont,i)->x += translate->x;
getPt(cont,i)->y += translate->y;
getPt(cont,i)->z += translate->z;
}
}
//------------------------
//-- Translates (moves) a contour by a given amount in x and y.
void cont_translate( Icont *cont, float x, float y )
{
Ipoint translate;
translate.x = x;
translate.y = y;
translate.z = 0;
cont_translate( cont, &translate );
}
//------------------------
//-- Rotates a contour by a given angle (in DEGREES) around specified
//-- rotate point (ignores Z axis).
void cont_rotateAroundPoint2D( Icont *cont, Ipoint *center, float angle )
{
float theta = angle * DEGS_TO_RADS;
for(int i=0; i<psize( cont ); i++)
point_rotatePointAroundPoint2D( getPt(cont,i), center, theta );
}
//------------------------
//-- Scales a contour about a given point by scaleX in X and scaleY in Y.
void cont_scaleAboutPtXY( Icont *cont, Ipoint *center, float scaleX, float scaleY )
{
for(int i=0; i<psize( cont ); i++)
point_scalePtAboutPt2D( getPt(cont,i), center, scaleX, scaleY );
}
//------------------------
//-- Scales a contour about a given point by the given factors in each dimension
void cont_scaleAboutPt3D( Icont *cont, Ipoint *center, float scaleX,
float scaleY, float scaleFactorZ )
{
for(int i=0; i<psize(cont); i++)
point_scalePtAboutPt( getPt(cont, i), center, scaleX, scaleY, scaleFactorZ );
}
//------------------------
//-- Returns the contour after it has been scaled outwards/inwards
void cont_scaleAboutPt( Icont *cont, Ipoint *pt, float scaleFactor, bool ignoreZ )
{
float zOrig = getZ( cont );
for(int i=0; i<psize(cont); i++)
*getPt(cont,i) = line_findPtFractBetweenPts( pt, getPt(cont,i), scaleFactor );
if(ignoreZ)
setZValue( cont, (int)zOrig );
}
//------------------------
//-- Stretches the contour by "stretchFactor" about "center" along
//-- the given "angle".
void cont_stretchAlongAngle( Icont *cont, Ipoint *center,
float angle, float stretchFactor )
{
cont_rotateAroundPoint2D( cont, center, -angle );
cont_scaleAboutPtXY( cont, center, stretchFactor, 1.0 );
cont_rotateAroundPoint2D( cont, center, angle );
}
//---------------------------------
//-- Generates a clockwise contour representing a circle with given center,
//-- radius and number of points.
void cont_generateCircle( Icont *cont, float radius, int numPoints,
Ipoint center, bool addEndPt )
{
Ipoint pt;
pt.z = center.z;
for (int i = 0; i < numPoints; i++) {
float angle = -2.0 * PI * ( (float)i / (float)numPoints );
pt.x = center.x + radius * cos(angle);
pt.y = center.y + radius * sin(angle);
imodPointAppend(cont, &pt);
}
if(addEndPt)
{
pt.x = center.x + radius;
pt.y = center.y;
imodPointAppend(cont, &pt);
}
}
//---------------------------------
//-- Generates a 2D rectangle of the given width and height with the bottom left corner
//-- at coordinates (llX, llY, z). The box is draw counterclockwise from the bottom
//-- left, and if "repeatFirstPt" is true it will "close" the box with a 5th point
void cont_generateBox( Icont *cont, float llX, float llY, float width, float height,
float z, bool repeatFirstPt )
{
imodPointAppendXYZ( cont, llX , llY , z );
imodPointAppendXYZ( cont, llX+width, llY , z );
imodPointAppendXYZ( cont, llX+width, llY+height, z );
imodPointAppendXYZ( cont, llX , llY+height, z );
if( repeatFirstPt )
imodPointAppendXYZ( cont, llX , llY , z );
}
//---------------------------------
//-- Generates a 2D rectangle between the two points given (lower left and upper right).
//-- The box is draw counterclockwise from the bottom left, at the z of the first point,
//-- and if "repeatFirstPt" is true it will "close" the box with a 5th point.
void cont_generateBox( Icont *cont, Ipoint ll, Ipoint ur, bool repeatFirstPt )
{
imodPointAppendXYZ( cont, ll.x, ll.y, ll.z );
imodPointAppendXYZ( cont, ur.x, ll.y, ll.z );
imodPointAppendXYZ( cont, ur.x, ur.y, ll.z );
imodPointAppendXYZ( cont, ll.x, ur.y, ll.z );
if( repeatFirstPt )
imodPointAppendXYZ( cont, ll.x, ll.y, ll.z );
}
//------------------------
//-- Counts the number of times two (closed) contours cross each other.
//-- NOTE: This number should ALWAYS be even.
int cont_numTimesCountsCross( Icont *cont1, Icont *cont2 )
{
int numTimesLinesCross = 0;
for (int i=0; i<=psize(cont1);i++)
for (int j=0; j<=psize(cont2);j++)
if( imodPointIntersect( getPt(cont1,i), getPt(cont1,i+1),
getPt(cont2,j), getPt(cont2,j+1) ) )
numTimesLinesCross++;
return numTimesLinesCross;
}
//------------------------
//-- Returns true if the two (closed) contours cross paths.
bool cont_doCountoursCross( Icont *cont1, Icont *cont2,
bool cont1Closed=true, bool cont2Closed=true )
{
Ipoint interceptPt;
int cont1PtsToCheck = (cont1Closed) ? psize(cont1): psize(cont1)-1;
int cont2PtsToCheck = (cont2Closed) ? psize(cont2): psize(cont2)-1;
for (int i=0; i<cont1PtsToCheck;i++)
for (int j=0; j<cont2PtsToCheck;j++)
if( imodPointIntersect( getPt(cont1,i), getPt(cont1,i+1),
getPt(cont2,j), getPt(cont2,j+1) ) )
return true;
return false;
}
//------------------------
//-- Returns true if the two (closed) contours cross paths.
bool cont_doCountoursCrossAndWhere( Icont *cont1, Icont *cont2,
bool cont1Closed, bool cont2Closed,
int *pt1BeforeCross, int *pt2BeforeCross )
{
Ipoint interceptPt;
int cont1PtsToCheck = (cont1Closed) ? psize(cont1): psize(cont1)-1;
int cont2PtsToCheck = (cont2Closed) ? psize(cont2): psize(cont2)-1;
for (int i=0; i<cont1PtsToCheck;i++)
for (int j=0; j<cont2PtsToCheck;j++)
if(imodPointIntersect(getPt(cont1,i),getPt(cont1,i+1),
getPt(cont2,j),getPt(cont2,j+1)))
{
*pt1BeforeCross = i;
*pt2BeforeCross = j;
return true;
}
return false;
}
//------------------------
//-- Determines if point touches the given contour LINE
//-- (i.e: if it is on an existing point or between two consecutive contour points)
bool cont_doesPtTouchContLine( Ipoint *pt, Icont *cont )
{
for (int i=0; i<psize(cont); i++)
if( pt == getPt(cont,i) )
return true;
for (int i=0; i<psize(cont); i++)
if( line_isPointOnLine( pt, getPt(cont,i), getPt(cont,i+1) ) )
return true;
return false;
}
//------------------------
//-- Finds and returns the closest two points (ignoring Z axis) and returns
//-- their distance apart.
float cont_findClosestPts2D( Icont *cont1, Icont *cont2,
int *closestPtIdxInCont1, int *closestPtIdxInCont2 )
{
float closestDist2 = FLOAT_MAX;
for(int i=0; i<psize(cont1); i++ )
for(int j=0; j<psize(cont2); j++ )
{
float dist2 = line_sqDistBetweenPts2D( getPt(cont1,i), getPt(cont2,j) );
if ( dist2 < closestDist2 ) {
closestDist2 = dist2;
*closestPtIdxInCont1 = i;
*closestPtIdxInCont2 = j;
}
}
return sqrt(closestDist2);
}
//------------------------
//-- Reverses the order of the points
void cont_reversePts( Icont *c )
{
imodel_contour_invert( c );
}
//------------------------
//-- Takes three points and and a contour adds these points to the end of the contour,
//-- but will add extra points between them to create a nice curve if the angle formed
//-- by the three points is large.
//-- This function is used by the "cont_expandOpenCont" function to
//-- avoid square ends and sharp turns.
void cont_addChamferPts( Icont *cont, Ipoint *ptPrev, Ipoint *ptCurr, Ipoint *ptNext,
float distOffset, float minAngle )
{
if(minAngle <= 0)
return;
float angleMadeByPoint = 360 - line_angleFormed3Pts( ptPrev, ptCurr, ptNext );
// the angle made by the corner on the "outside"
if( angleMadeByPoint > 180 ) // if angle is reflex: add multiple points
{
float angleCurveSpan = angleMadeByPoint-180;
int numMidPtsToAdd = (int)floor( angleCurveSpan / minAngle );
if( numMidPtsToAdd == 0 ) // if no extra points are needed to smooth the corner:
{ // add two points
Ipoint pt1 = line_getPtRelativeToEnd(ptPrev,ptCurr,distOffset, 90);
Ipoint pt2 = line_getPtRelativeToEnd(ptPrev,ptCurr,distOffset, 90);
Ipoint pt3 = line_getPtRelativeToEnd(ptNext,ptCurr,distOffset,-90);
imodPointAppend( cont, &pt2 );
imodPointAppend( cont, &pt3 );
}
else // if extra points are needed to smooth the corner:
{ // add multiple points in an arc around the current (central) pt
float angleBetweenPts = angleCurveSpan / (numMidPtsToAdd+1);
for(int i=0; i<numMidPtsToAdd+1; i++)
imodPointAppend( cont, &line_getPtRelativeToEnd( ptPrev, ptCurr, distOffset,
90-(i*angleBetweenPts) ) );
}
}
else // else if angle on outside is <= 180 degrees (obtuse or acute) then:
{ // add a SINGLE point in the middle
float angleInset = angleMadeByPoint/2.0f;
float distFromCorner = fDiv( distOffset, sin(angleInset*DEGS_TO_RADS) );
Ipoint midSegmentPt = line_getPtRelativeToEnd( ptPrev, ptCurr, distFromCorner,
180.0-angleInset );
imodPointAppend( cont, &midSegmentPt );
}
}
//------------------------
//-- Takes an open contour and expands it by a certain thickness and returns a
//-- closed contour around the perimeter.
//-- For all reflex angles, extra points will be added to make a smoother curve
//-- according to the value of minAngleForChamfers.
//-- NOTE: There is NO checking for overlap in this function - that may be called
//-- seperately by calling: cont_makeSimple() or cont_breakIntoSimple()
void cont_expandOpenCont( Icont *contOrig, Icont *contR,
float thickness, float minAngleForChamfers, bool roundEnds )
{
deleteAllPts(contR);
Icont *cont = imodContourDup(contOrig);
//## DEAL WITH POSSIBLILTY THAT CONTOUR HAS ONE OR NO POINTS:
if ( thickness == 0 )
return;
if( psize(cont) == 0 )
return;
if ( psize(cont) == 1 ) {
cont_generateCircle( contR, thickness, 12, *getFirstPt(cont), false );
return;
}
//## PREPARE CONTOUR:
imodContourUnique( cont );
imodContourMakeDirection( cont, IMOD_CONTOUR_CLOCKWISE );
int N = psize(cont);
//## ADD START:
if( roundEnds ) {
cont_addChamferPts( contR, getPt(cont,1), getPt(cont,0), getPt(cont,1),
thickness, minAngleForChamfers);
}
else {
imodPointAppend(contR,
&line_getPtRelativeToEnd( getPt(cont,0),getPt(cont,1),thickness,-90 ) );
imodPointAppend(contR,
&line_getPtRelativeToEnd( getPt(cont,0),getPt(cont,1),thickness, 90 ) );
}
//## ADD POINTS AROUND OUTSIDE:
for(int i=1; i<N-1; i++)
cont_addChamferPts( contR, getPt(cont,i-1), getPt(cont,i), getPt(cont,i+1),
thickness, minAngleForChamfers);
//## ADD END POINT:
if( roundEnds ) {
cont_addChamferPts( contR, getPt(cont,N-2),getPt(cont,N-1),getPt(cont,N-2),
thickness,minAngleForChamfers);
}
else {
imodPointAppend(contR,
&line_getPtRelativeToEnd( getPt(cont,N-1), getPt(cont,N-2), thickness, -90 ) );
imodPointAppend(contR,
&line_getPtRelativeToEnd( getPt(cont,N-1), getPt(cont,N-2), thickness, 90 ) );
}
//## ADD POINTS AROUND INSIDE:
for(int i=N-2; i>0; i--)
cont_addChamferPts( contR, getPt(cont,i+1), getPt(cont,i), getPt(cont,i-1),
thickness, minAngleForChamfers);
return;
}
//------------------------
//-- Takes an closed contour and expands it by a certain thickness and returns
//-- several closed contours around the perimeter. For all reflex angles, extra
//-- points will be added to make a smoother curve according to the value of
//-- minAngleForChamfers.
//-- NOTE: If the inner contour overlapps itself it is broken into several simple
//-- contours, but NONE of these are deleted - that must be done manually.
void cont_expandClosedCont( Icont *contOrig, Icont *innerCont, Icont *outerCont,
float thickness, float minAngleForChamfers )
{
deleteAllPts(innerCont);
deleteAllPts(outerCont);
Icont *cont = imodContourDup(contOrig);
//## DEAL WITH POSSIBLILTY THAT CONTOUR HAS ONE OR NO POINTS:
if ( thickness == 0 ) {
return;
}
if( psize(cont) == 0 ) {
return;
}
if ( psize(cont) == 1 ) {
cont_generateCircle( outerCont, thickness, 12, *getFirstPt(cont), false );
return;
}
//## PREPARE CONTOUR:
imodContourUnique( cont );
imodContourMakeDirection( cont, IMOD_CONTOUR_CLOCKWISE );
//## GET INNER CONTOURS:
for(int i=psize(cont); i>0; i--)
cont_addChamferPts( innerCont, getPt(cont,i+1), getPt(cont,i), getPt(cont,i-1),
thickness, minAngleForChamfers);
//## ADD OUTER CONTOUR:
for(int i=0; i<psize(cont); i++)
cont_addChamferPts( outerCont, getPt(cont,i-1), getPt(cont,i), getPt(cont,i+1),
thickness, minAngleForChamfers);
return;
}
//------------------------
//-- Joins two contours together to form a longer contour.
//-- If matchClosestEnds is true it will join the two ends of cont1 and cont2
//-- which are closest together... otherwise it will simply add the points
//-- from cont2 onto the end of cont1
void cont_concat( Icont *contNew, Icont *cont1, Icont *cont2Orig, bool matchClosestEnds )
{
imodContourDefault( contNew );
Icont *cont2 = imodContourDup( cont2Orig );
if(matchClosestEnds) {
float distCont1EndToCont2Beg = imodPointDistance(getLastPt(cont1),getFirstPt(cont2));
float distCont1EndToCont2End = imodPointDistance(getLastPt(cont1),getLastPt(cont2));
if( distCont1EndToCont2End < distCont1EndToCont2Beg )
cont_reversePts( cont2 );
}
for(int i=0; i<psize(cont1); i++ )
imodPointAppend( contNew, getPt( cont1, i) );
for(int i=0; i<psize(cont2); i++ )
imodPointAppend( contNew, getPt( cont2, i) );
imodContourDelete(cont2);
imodContourUnique(contNew);
}
//------------------------
//-- Adds points to the using a very simple technique - to ensure
//-- no two consecutive points are > maxDist apart - and return the number of
//-- points added (if any).
int cont_addPtsCrude( Icont *cont, float maxDist, bool closed )
{
int pointsBefore = psize(cont);
int extra = (closed) ? 0 : -1;
for(int i=0; i<psize(cont)+extra; i++ )
{
Ipoint *currPt = getPt( cont, i );
Ipoint *nextPt = getPt( cont, i+1 );
float distToNextPt = line_distBetweenPts2D( currPt, nextPt );
if( distToNextPt > maxDist )
{
Ipoint newPt = line_getPtHalfwayBetween( currPt, nextPt );
imodPointAdd( cont, &newPt, i+1 );
i--;
}
}
return ( psize(cont) - pointsBefore );
}
//------------------------
//-- Smooths the contour by finding any occurance where two consecutive
//-- points are > maxDist away from each other and adding a SINGLE extra
//-- point in the middle using catumull-rom spline and return the number of
//-- points added (if any).
//-- NOTE: At the end there may still be consecutive points > maxDist
//-- apart (but will be closer than originally)
int cont_addPtsSmoothIteration(Icont *cont, float maxDist,
float tensileFract, bool closed)
{
int pointsBefore = psize(cont);
float sqMaxDist = SQ( maxDist );
if( psize(cont) < 5 )
return 0;
if(closed)
{
for(int i=1; i<(psize(cont))+1; i++ ) {
float sqDistToNextPt = line_sqDistBetweenPts2D( getPt(cont,i), getPt(cont,i+1) );
if ( sqDistToNextPt > sqMaxDist )
{
Ipoint newPt = getPtCardinalSpline(0.5, *getPt(cont,i-1), *getPt(cont,i),
*getPt(cont,i+1), *getPt(cont,i+2), tensileFract);
int insertIdx = i % psize(cont);
imodPointAdd( cont, &newPt, (insertIdx)+1 );
i++; // causes point just added to be skipped
}
}
}
else
{
imodPointAdd( cont, getFirstPt(cont), 0 );
imodPointAppend( cont, getLastPt(cont) );
for(int i=1; i<(psize(cont))-1; i++ ) {
float sqDistToNextPt = line_sqDistBetweenPts2D( getPtNoWrap(cont,i),
getPtNoWrap(cont,i+1) );
if ( sqDistToNextPt > sqMaxDist )
{
Ipoint newPt = getPtCardinalSpline(0.5, *getPtNoWrap(cont,i-1), *getPtNoWrap(cont,i),
*getPtNoWrap(cont,i+1), *getPtNoWrap(cont,i+2),
tensileFract);
int insertIdx = i % psize(cont);
imodPointAdd( cont, &newPt, (insertIdx)+1 );
i++; // causes point just added to be skipped
}
}
imodContourUnique( cont );
}
return ( psize(cont) - pointsBefore );
}
//------------------------
//-- Smooths the contour by finding any occurance where two consecutive
//-- points are > maxDist away from each other and adding MULTIPLE
//-- extra points using catumull-rom spline and returns the
//-- number of points added.
//-- NOTE: Running this a second time on the same contour will usually
//-- add extra points due to the curvature of points added.
int cont_addPtsSmooth( Icont *cont, float maxDist, float tensileFract, bool closed,
bool roundZOpenPts, bool addPtEveryZ )
{
if(maxDist<=0.0f)
maxDist = 1.0f;
int pointsBefore = psize(cont);
int pointsAdded = 0;
if(closed)
{
Icont *contO = imodContourDup(cont);
for(int i=0; i<(psize(contO)); i++ )
{
float distToNextPt = line_distBetweenPts2D( getPt(contO,i), getPt(contO,i+1) );
int numPtsToAdd = ceil(distToNextPt / maxDist) - 1;
if( distToNextPt==0 || numPtsToAdd<=0 )
continue;
float fractBetweenPts = fDiv( 1.0f, numPtsToAdd+1.0f );
for(int j=1; j<=numPtsToAdd; j++)
{
float fracAlongSegment = j * fractBetweenPts;
Ipoint newPt = getPtCardinalSpline(fracAlongSegment,
*getPt(contO,i-1),
*getPt(contO,i),
*getPt(contO,i+1),
*getPt(contO,i+2),
tensileFract);
int insertIdx = (i+pointsAdded) % psize(cont);
imodPointAdd( cont, &newPt, (insertIdx)+1 );
pointsAdded++;
}
}
imodContourDelete(contO);
}
else
{
imodPointAdd( cont, getFirstPt(cont), 0 );
imodPointAppend( cont, getLastPt(cont) );
Icont *contO = imodContourDup(cont);
for(int i=1; i<(psize(contO))-1; i++ )
{
float distToNextPt = line_distBetweenPts2D( getPtNoWrap(contO,i),
getPtNoWrap(contO,i+1) );
int numPtsToAdd = ceil(distToNextPt / maxDist) - 1;
if( addPtEveryZ )
{
int zDiffPts = ABS( roundToInt( getPt(contO,i)->z ) -
roundToInt( getPt(contO,i+1)->z ) );
numPtsToAdd = MAX( numPtsToAdd, zDiffPts-1 );
}
if( numPtsToAdd<=0 )
continue;
float fractBetweenPts = fDiv( 1.0f, numPtsToAdd+1.0f );
for(int j=1; j<=numPtsToAdd; j++)
{
float fracAlongSegment = j * fractBetweenPts;
Ipoint newPt = getPtCardinalSpline(fracAlongSegment,
*getPtNoWrap(contO,i-1),
*getPtNoWrap(contO,i),
*getPtNoWrap(contO,i+1),
*getPtNoWrap(contO,i+2),
tensileFract);
int insertIdx = (i+pointsAdded) % psize(cont);
imodPointAdd( cont, &newPt, (insertIdx)+1 );
pointsAdded++;
}
}
if( roundZOpenPts )
{
for( int p=0; p<psize(cont); p++ )
getPt(cont,p)->z = roundToInt(getPt(cont,p)->z);
}
imodContourDelete(contO);
imodContourUnique( cont );
}
return ( psize(cont) - pointsBefore );
}
//------------------------
//-- Smooths the contour by moving each selected point "moveFract" towards
//-- the position half-way between the point before and after it.
//-- Only points > "minDistToMove" away from their "expected position" are moved.
//-- Returns the number of points moved.
//-- If "rescale" is true the contour is expanded so the area before and after
//-- points are moved stays constant (even though the shape changes).
//-- NOTE: To help preserve shape, this should be run on contours with dense points
//-- ... thus it's a good idea to run "cont_addPtsCrude()" first.
int cont_avgPtsPos( Icont *cont, float moveFract, float minDistToMove,
bool closed, bool rescale )
{
if(psize(cont) <= 5)
return 0;
int pointsMoved = 0;
int offset = (closed) ? 1 : -1;
float areaBefore = imodContourArea( cont );
for( int p=1; p<psize(cont)+offset; p++ )
{
Ipoint *currPt = getPt(cont,p);
Ipoint expectedPos = line_getPtHalfwayBetween( getPt(cont,p-1), getPt(cont,p+1) );
float distToExpected = line_distBetweenPts2D( currPt, &expectedPos );
if( distToExpected > minDistToMove )
{
*currPt = line_findPtFractBetweenPts2D( currPt, &expectedPos, moveFract );
pointsMoved++;
}
}
float areaAfter = imodContourArea( cont );
if(rescale && pointsMoved && areaBefore > areaAfter && areaAfter > 0 )
{
float scaleXY = sqrt(areaBefore) / sqrt(areaAfter);
Ipoint centroid;
cont_getCentroid( cont, ¢roid );
cont_scaleAboutPt( cont, ¢roid, scaleXY, true );
}
return ( pointsMoved );
}
//------------------------
//-- Reduces the number of points in a contour using a VERY simple technique -
//-- by removing any point which is < minDist from the previous point -
//-- and returns the number of points deleted.
int cont_reducePtsCrude( Icont *cont, float minDist, bool closed )
{
int pointsBefore = psize(cont);
float sqMinDist = SQ( minDist );
int extra = (closed) ? 1 : 0;
for(int i=1; i<psize(cont)+extra; i++ ) {
float sqDistFromPrevPoint = line_sqDistBetweenPts2D(getPt(cont,i-1),getPt(cont,i));
if ( sqDistFromPrevPoint < sqMinDist ) {
imodPointDelete( cont, i );
i--;
}
}
return ( pointsBefore - psize(cont) );
}
//------------------------
//-- Reduces the number of points in a contour using the imodContourReduction()
//-- function using the given tolerance value and then returns the number
//-- of points deleted.
int cont_reducePtsTol( Icont *cont, float tol )
{
int pointsBefore = psize(cont);
imodContourReduce(cont, tol);
return ( pointsBefore - psize(cont) );
}
//------------------------
//-- Reduces the number of points in a contour by removing points
//-- which form an area < minArea when a triangle if formed with that point,
//-- the point before it, and the point after.
//-- Returns the number of points deleted.
int cont_reducePtsMinArea( Icont *cont, float minArea, bool closed )
{
int pointsBefore = psize(cont);
int extra = (closed) ? 0 : -1;
for(int i=1; i<psize(cont)+extra && psize(cont)>3; i++ )
{
Ipoint *p1 = getPt(cont,i-1); //|-- construct a triangle with the current pt
Ipoint *p2 = getPt(cont,i); //| plus the pt before and after it.
Ipoint *p3 = getPt(cont,i+1); //|
if ( imodPointArea( p1, p2, p3 ) < minArea ) {
imodPointDelete( cont, i );
i--;
}
}
return ( pointsBefore - psize(cont) );
}
//------------------------
//-- Removes duplicate points and any point which form a straight line with the
//-- point before and after it. Returns the number of points removed
//-- If removePts is false, then the function only counts points without removing them.
int cont_removeRedundantPts( Icont *cont, bool removeStraightLinePts, bool closed,
bool removePts )
{
int redundantPts = 0;
for(int p=psize(cont)-2; p>=0; p-- )
{
if ( imodPointIsEqual( getPt(cont,p), getPt(cont,p+1) ) )
{
if( removePts )
imodPointDelete(cont,p);
redundantPts++;
}
}
if( !removeStraightLinePts )
return (redundantPts);
int startIdx = psize(cont) - (closed) ? 0 : 1;
for(int p=startIdx; p>=1 && psize(cont)>3; p-- )
{
//cout << imodPointCross( getPt(cont,p-1), getPt(cont,p), getPt(cont,p+1) ) << endl;
if ( line_crossProduct3Points( getPt(cont,p-1), getPt(cont,p), getPt(cont,p+1) ) == 0 )
{
if( removePts )
imodPointDelete(cont, p);
redundantPts++;
}
}
return ( redundantPts );
}
//------------------------
//-- Will return true if the contour is a "simple polygon"
//-- (none of it's lines overlap) or false if it is complex
//-- (i.e. two or more lines intersect)
bool cont_isSimple( Icont *cont, bool closed )
{
int ptsToCheck = (closed) ? psize(cont) : psize(cont)-1;
for(int i=0; i<ptsToCheck; i++ )
for(int j=i+2; j<ptsToCheck; j++ )
if(imodPointIntersect(getPt(cont,i),getPt(cont,i+1),getPt(cont,j),getPt(cont,j+1))
&& !( i == 0 && j == psize(cont)-1 ) )
return false;
return true;
}
//------------------------
//-- If the contour is not-simple, will return false
//-- and the points postion of the first intersecting segment.
bool cont_isSimpleSeg( Icont *cont, bool closed, int *ptCross )
{
int ptsToCheck = (closed) ? psize(cont) : psize(cont)-1;
for(int i=0; i<ptsToCheck; i++ )
for(int j=i+2; j<ptsToCheck; j++ )
if(imodPointIntersect(getPt(cont,i),getPt(cont,i+1),getPt(cont,j),getPt(cont,j+1))
&& !( i == 0 && j == psize(cont)-1 ) )
{
*ptCross = i;
return false;
}
return true;
}
//------------------------
//-- Takes a closed contour and (if not already) makes it simple.
//-- It does this by determine points where the contour intersects itself,
//-- and deleting the SMALLER enclosed area (as computed by the fewest
//-- number of points) until no more self-intersections are found.
//-- NOTE: This function is still not perfect, but does the right
//-- thing 90% of the time.
void cont_makeSimple( Icont *cont )
{
int numIntersects = 0;
imodContourUnique( cont );
findIntersect: // LABEL (see "goto findIntersect")
for(int i=0; i<psize(cont); i++ )
for(int j=i+2; j<psize(cont); j++ )
{
if( i == 0 && j == psize(cont)-1 ) continue;
Ipoint intersectPt;
bool intersectionFound =
line_doLinesCrossAndWhere( getPt(cont,i),getPt(cont,i+1),
getPt(cont,j),getPt(cont,j+1),&intersectPt );
if ( intersectionFound )
{
int numPtsToDelete = j-i;
// if the region between the intersecting lines has more points:
// delete points in middle and add intersection point
if ( numPtsToDelete < (psize(cont) / 2) )
{
for (int p=0; p<(numPtsToDelete); p++)
imodPointDelete( cont, i+1 );
imodPointAdd(cont, &intersectPt, i+1);
}
// else (if region on the ends has more points):
// delete points either end of intersecting pts and add intersection pt
else
{
int altNumPtsToDelete = psize(cont)-(j+1);
for (int p=0; p<(altNumPtsToDelete); p++)
imodPointDelete( cont, j+1 );
imodPointAdd(cont, &intersectPt, i+1);
for (int p=0; p<(i+1); p++)
imodPointDelete( cont, 0 );
}
numIntersects++; // |-- used to avoid rare cases where
if(numIntersects>50) { // | problem occurs (not yet sure why)
wprint( "ERROR: cont_makeSimple()\n" );
return;
}
goto findIntersect; // will go to beginning of double loop and
// see if there is another intersection.
}
}
return;
}
//------------------------
//-- Takes a contour and breaks it into a series of simple contours
//-- (i.e. contours which don't overlap themselves) and returns these.
//--
//-- NOTE: If the given contour is already simple it will return itself.
int cont_breakIntoSimple( vector<IcontPtr> &conts, Icont *cont )
{
conts.clear();
conts.push_back( IcontPtr(cont) );
if( cont_isSimple ( cont ) )
return (conts.size());
Ipoint intersectPt;
Icont *newCont1 = imodContourNew();
Icont *newCont2 = imodContourNew();
for(int c=0; c<(int)conts.size(); c++ )
{
findIntersect: // LABEL
imodContourUnique( conts[c].cont );
for(int i=0; i< psize(conts[c].cont); i++ )
{
for(int j=i+2; j< psize(conts[c].cont); j++ )
{
if( i == 0 && j == psize(conts[c].cont)-1 ) continue;
if ( line_doLinesCrossAndWhere( getPt(conts[c].cont,i), getPt(conts[c].cont,i+1),
getPt(conts[c].cont,j), getPt(conts[c].cont,j+1),
&intersectPt ) )
{
imodPointAdd( conts[c].cont, &intersectPt, j+1);
imodPointAdd( conts[c].cont, &intersectPt, i+1);
cont_breakContourEitherSide(conts[c].cont,newCont1,newCont2,i+1,j+2,true);
imodContourUnique( newCont1 );
imodContourUnique( newCont2 );
eraseContour( conts, c );
conts.push_back( IcontPtr(newCont1) );
conts.push_back( IcontPtr(newCont2) );
goto findIntersect; // will go to beginning of double loop
// and see if there is another intersection.
}
}
}
}
imodContourDelete(newCont1);
imodContourDelete(newCont2);
return (conts.size());
}
//------------------------
//-- Ensure no two consequtive points have the same x or y value
//-- NOTE: I use this function because I tend to have trouble with perfectly
//-- vertical and (to a lesser extent) horizontal contour segments
//-- where I need to find intersections and measure gradients.
int cont_killVertAndHorzSegments( Icont *cont )
{
int nudgeEvents = 0;
for( int i=1; i<psize(cont); i++)
{
if( getPt(cont,i)->y == getPt(cont,i-1)->y )
{
getPt(cont,i)->y += 0.0001f;
nudgeEvents++;
}
if( getPt(cont,i)->x == getPt(cont,i-1)->x )
{
getPt(cont,i)->x += 0.0001f;
nudgeEvents++;
}
}
return (nudgeEvents);
}
//------------------------
//-- Randomly shifts all points in the contour by a very small amount in X and Y
void cont_nudgeAllPointsRandomly( Icont *cont )
{
const double MAX_NUDGE = 0.01;
/*static bool seed = false;
if(!seed)
seedRandom();
seed = true;
//*/
for(int i=0; i<psize(cont); i++)
{
getPt(cont,i)->x += randDbl( -MAX_NUDGE, MAX_NUDGE );
getPt(cont,i)->y += randDbl( -MAX_NUDGE, MAX_NUDGE );
}
}
//------------------------
//-- Ensures no point on "cont1" lies on another point or line segment of "cont2"
//-- and visa versa, by shifting these points slightly in X and/or Y.
//-- Returns the number of points which have been shifted.
int cont_nudgeAnyPointsLyingOnOtherContour( Icont *cont1, Icont *cont2 )
{
int ptsNudged = 0;
const float NUDGE_AMOUNT = 0.001;
for(int i=0; i<psize(cont1); i++)
{
for(int j=0; j<psize(cont2); j++)
{
if( imodPointIsEqual(getPt(cont1,i),getPt(cont2,j)) )
{
getPt(cont1,i)->x += NUDGE_AMOUNT;
ptsNudged++;
}
if( line_isPointOnLine(getPt(cont1,i),getPt(cont2,j),getPt(cont2,j+1)) )
{
getPt(cont1,i)->x += NUDGE_AMOUNT;
ptsNudged++;
}
if( line_isPointOnLine(getPt(cont2,j),getPt(cont1,i),getPt(cont1,i+1)) )
{
getPt(cont2,j)->x += NUDGE_AMOUNT;
ptsNudged++;
}
}
}
for(int i=0; i<psize(cont1); i++)
{
for(int j=0; j<psize(cont2); j++)
{
if( imodPointIsEqual(getPt(cont1,i),getPt(cont2,j)) )
{
getPt(cont1,i)->y += NUDGE_AMOUNT;
ptsNudged++;
}
if( line_isPointOnLine(getPt(cont1,i),getPt(cont2,j),getPt(cont2,j+1)) )
{
getPt(cont1,i)->y += NUDGE_AMOUNT;
ptsNudged++;
}
if( line_isPointOnLine(getPt(cont2,j),getPt(cont1,i),getPt(cont1,i+1)) )
{
getPt(cont2,j)->y += NUDGE_AMOUNT;
ptsNudged++;
}
}
}
return( ptsNudged );
}
//------------------------
//-- Determines if a contour is convex by checking that all the "turns"
//-- (the angle formed by the lines either side of each point)
//-- are in the same direction.
//-- Note that this algorithm may fail in rare cases where the contour
//-- is non-simple and forms a "loop".
bool cont_isConvex( Icont *cont )
{
if( psize(cont) <= 3 )
return true;
int numRightTurns = 0;
int numLeftTurns = 0;
for (int i=1; i<(psize(cont)+1); i++ )
{
float crossProduct = line_crossProduct3Points(getPt(cont,i-1),
getPt(cont,i),getPt(cont,i+1));
// calculates cross product using the line {i-1,i}
// and line {i,i+1} to determine if a "left turn" is made
if (crossProduct > 0) // if turn is left: tally it
numLeftTurns++;
else if (crossProduct == 0) // else if turn is straight: skip to next
continue;
else // else (if turn is right): tally it
numRightTurns++;
if( numLeftTurns>0 && numRightTurns>0 ) // if NOT all turns are in same direction:
return false; // contour is not convex
}
return true; // all turns are in the same direction: contour must be convex
}
//------------------------
//-- Attempts to makes a contour convex by eliminating all points which form a
//-- left turn (with the lines either side of it) for an clockwise contour
//-- or a right turn for a anti-clockwise contour. Note that this crude method
//-- can often fail on more convoluted contours with a concavity tree greater
//-- than one deep (i.e. the concave regions have concave regions).
//-- NOTE: It's a good idea to run "imodContourUnique" before running this.
int cont_makeConvexCrude( Icont *cont )
{
if( psize(cont) < 3 )
return 0;
int pointsBefore = psize(cont);
bool isClockwise = (imodContZDirection(cont) == IMOD_CONTOUR_CLOCKWISE);
for (int p=0; p<(psize(cont)+3) && psize(cont)>3; p++ )
{
float crossProduct = line_crossProduct3Points(getPt(cont,p-1),
getPt(cont,p),getPt(cont,p+1));
// calculates cross product using the line {i-1,i}
// and line {i,i+1} to determine if a "left turn" is made
if ( (crossProduct > 0 && isClockwise) || // if turn is left and clockwise:
(crossProduct < 0 && !isClockwise) ) // or turn is right and anticlockwise:
{
int ptIdx = intMod( p,psize(cont) );
imodPointDelete(cont, ptIdx );
p-=2;
p = MAX(0,p);
}
}
return ( pointsBefore - psize(cont) );
}
//------------------------
//-- Computes and returns a clockwise convex hull around a given set of points
//-- using the "Graham Scan" (3-coins algorithm) algorithm.
//-- This involves finding the lowest point (in y), sorting all points radially
//-- relative this point, then removing points which form a "left turn".
//-- Returns the number of points removed.
//-- see: http://en.wikipedia.org/wiki/Graham_scan
//-- NOTE: JUST DISCOVERED SOME SERIOUS PROBLEMS - RECOMMEND :-(
int cont_makeConvex( Icont *cont )
{
if( cont_isConvex(cont) )
return 0;
//## MAKE CONTOUR CLOCKWISE:
imodContourStrip(cont);
imodContourMakeDirection( cont, IMOD_CONTOUR_CLOCKWISE );
//## FIND THE LOWEST POINT IN THE CONTOUR:
int pointsBefore = psize(cont);
int idxStartPt = 0;
float lowestYVal = getFirstPt(cont)->y;
for (int p=1; p<pointsBefore; p++)
{
Ipoint *currPt = getPt(cont,p);
if(currPt->y <= lowestYVal) // if this point is lowest so far:
{
if( currPt->y == lowestYVal // if point has same y value:
&& currPt->x > getPt(cont,idxStartPt)->x ) // use x as a tie-breaker.
continue;
lowestYVal = currPt->y; // update as new lowest point.
idxStartPt = p;
}
}
//## MAKE CONTOUR CLOCKWISE AND START AT LOWEST POINT:
cont_reorderPtsToStartAtIdx( cont, idxStartPt );
//## SORT POINTS RADIALLY FROM THE LOWEST POINT:
Ipoint *lowestPt = getFirstPt(cont);
vector<IdxToSort> idxAngles;
for (int p=1; p<pointsBefore; p++)
{
float angle = line_getAngle2DPos( lowestPt, getPt(cont,p) );
float sqDist = line_sqDistBetweenPts2D( lowestPt, getPt(cont,p) );
idxAngles.push_back( IdxToSort( p, angle, FLOAT_MAX-sqDist ) );
}
idxAngles = vector_sort( idxAngles );
Icont *contCopy = imodContourDup( cont );
deleteAllPts(cont);
imodPointAppend( cont, getFirstPt(contCopy) ); // add lowest point
for (int i=0; i<(int)idxAngles.size(); i++) // use idx values to populate cont
imodPointAppend( cont, getPt(contCopy,idxAngles[i].idx ) );
imodContourDelete(contCopy);
//## ITERATE THROUGH LIST AND REMOVE ANY POINTS WHICH MAKE A LEFT TURN:
for (int p=1; p<psize(cont); p++ )
{
float crossProduct = line_crossProduct3Points(getPt(cont,p-1),
getPt(cont,p),getPt(cont,p+1));
if( crossProduct < 0 ) // if this point makes a left turn:
{
imodPointDelete( cont, p ); // delete this point... and
p = MAX(1,p-1); // go back two points
}
}
return ( pointsBefore - psize(cont) );
}
//------------------------
//-- Sets the z value of any concave points to -1 and
//-- returns the number of concave points found.
int cont_markConvexPtsNegOne( Icont *cont )
{
if( cont_isConvex(cont) )
return 0;
//## CREATE A CONVEX VERSION OF THE CONTOUR:
Icont *convexCont = imodContourDup(cont);
cont_makeConvex( convexCont );
//## FOR EACH POINT: IF IT'S NOT IN THE CONVEX VERSION, MARK IT AS A CONCAVE POINT
int numPts = psize(cont);
int numConvexPts = psize(convexCont);
int numConcavePts = numPts - numConvexPts;
int concaveIdx = 0;
for( int p=0; p<numPts && concaveIdx<numConcavePts; p++ )
{
Ipoint *currPt = getPt( cont,p );
int ptIsInside = cont_doesPtExistInCont(convexCont, currPt );
if( ptIsInside )
currPt->z = -1;
}
return numConvexPts;
}
//------------------------
//-- Takes a contour and returns the number of convex points "numConvexPts",
//-- the total length of convex line segments "convexLen", and both the length "hullLen"
//-- and area "hullArea" of the convex hull around the contour.
void cont_calcConvexProperties( Icont *cont, bool closed, int *numConvexPts,
float *convexLen, float *hullLen, float *hullArea )
{
if( cont_isConvex(cont) )
{
float totalLen = imodContourLength( cont, closed );
*numConvexPts = psize( cont );
*convexLen = totalLen;
*hullLen = totalLen;
*hullArea = imodContourArea( cont );
return;
}
//## CREATE A CONVEX HULL OVER THE CONTOUR AND COUNT IT'S POINTS:
Icont *convexCont = imodContourDup(cont);
cont_makeConvex( convexCont );
*numConvexPts = psize( convexCont );
*hullLen = imodContourLength( convexCont, closed );
*hullArea = imodContourArea( convexCont );
//## FOR EACH LINE SEGMENT: IF IT'S IN THE CONVEX VERSION, ADD IT TO THE CONVEX LENGTH
float totConvexLength = 0;
bool currPtInside = cont_doesPtExistInCont(convexCont, getFirstPt(cont) );
int numPts = (closed) ? psize(cont) : psize(cont)-1;
for( int p=0; p<numPts; p++ )
{
Ipoint *currPt = getPt(cont,p);
Ipoint *nextPt = getPt(cont,p+1);
bool nextPtInside = cont_doesPtExistInCont(convexCont, nextPt );
if( currPtInside && nextPtInside )
totConvexLength += line_distBetweenPts2D( currPt, nextPt );
currPtInside = nextPtInside;
}
*convexLen = totConvexLength;
imodContourDelete( convexCont );
}
//------------------------
//-- Takes a single contours and breaks it into two contours at the
//-- specified index points (idxPt1 and idxPt2).
//--
//-- If shareEdge is false:
//-- the second contour will include the points between idxPt1 and idxPt2-1
//-- and the first contour will include the pts between idxPt2 and idxPt1-1
//-- If shareEdge is set to true:
//-- both contours will include the pt at idxPt1 and idxPt2 and hence
//-- share that as a common edge
bool cont_breakContourEitherSide( Icont *cont, Icont *contBreak1, Icont *contBreak2,
int idxPt1, int idxPt2, bool shareEdge=true )
{
imodContourDefault( contBreak1 );
imodContourDefault( contBreak2 );
if( ( idxPt1 == idxPt2 )
|| ( !isBetweenAsc(0, idxPt1, psize(cont)-1) )
|| ( !isBetweenAsc(0, idxPt2, psize(cont)-1) ) ) {
wprint( "ERROR: cont_breakContourEitherSide()" );
return false;
}
if( !(idxPt1 < idxPt2) ) //|-- ensures that idxPt1 < idxPt2
swapVals( idxPt1, idxPt2 ); //|
if( shareEdge )
{
//## CREATE CONTOUR 1:
for( int i=0; i<=idxPt1; i++ )
imodPointAppend( contBreak1, getPt(cont,i) );
for( int i=idxPt2; i<psize(cont); i++ )
imodPointAppend( contBreak1, getPt(cont,i) );
//## CREATE CONTOUR 2:
for( int i=idxPt1; i<=idxPt2; i++ )
imodPointAppend( contBreak2, getPt(cont,i) );
}
else
{
//## CREATE CONTOUR 1:
for( int i=0; i<idxPt1; i++ )
imodPointAppend( contBreak1, getPt(cont,i) );
for( int i=idxPt2; i<psize(cont); i++ )
imodPointAppend( contBreak1, getPt(cont,i) );
//## CREATE CONTOUR 2:
for( int i=idxPt1; i<idxPt2; i++ )
imodPointAppend( contBreak2, getPt(cont,i) );
}
return true;
}
//------------------------
//-- Breaks a contour either side of the specified line.
bool cont_breakContourByLine( Icont *cont, Icont *contBreak1, Icont *contBreak2,
Ipoint *linePt1, Ipoint *linePt2,
Ipoint expectedPt,
bool useExpectedPtInsteadOfMaxAreaSmallerSide=false )
{
imodContourDefault(contBreak1);
imodContourDefault(contBreak2);
Ipoint intercept; // temp value for "line_doLinesCrossAndWhere" function
vector<PtConnection> intercepts; // stores a list of pts where the two contours
// intersect, plus an index in cont where it crossed.
for (int i=0; i<psize(cont);i++) // for each point/line in cont:
if(line_doLinesCrossAndWhere(getPt(cont,i),getPt(cont,i+1),
linePt1,linePt2,&intercept ) )
{ // if lines cross: add intercept point
imodPointAdd( cont, &intercept, i+1 );
intercepts.push_back( PtConnection(intercept, i+1) );
i++; // we want to skip the intersection point now.
}
int numIntercepts = (int)intercepts.size();
if ( numIntercepts < 2 ) { // if there are > 2 intercept points: return false
return false;
}
else if( numIntercepts == 2 ) { // if there are 2 intercept points:
// break the contour between these points
cont_breakContourEitherSide( cont,contBreak1,contBreak2,
intercepts[0].cont1Idx,intercepts[1].cont1Idx,true );
return true;
}
else // else (if there are > 2 points):
{ // test sequential break pts to find which pair looks most likely
intercepts = vector_sort( intercepts );
Icont *contB1 = imodContourNew();
Icont *contB2 = imodContourNew();
float minDistToExpectedPt = FLOAT_MAX;
float maxAreaOnSmallerSide = 0;
for(int i=0; i<(int)intercepts.size()-1;i++)
{
Ipoint midWayPt = line_findPtFractBetweenPts2D( &intercepts[i].intercept,
&intercepts[i+1].intercept, 0.5 );
if ( !imodPointInsideCont( cont, &midWayPt ) ) {
continue;
}
if( useExpectedPtInsteadOfMaxAreaSmallerSide )
{
float distToExpectedPt = imodPointDistance( &expectedPt, &midWayPt );
if( distToExpectedPt < minDistToExpectedPt ) {
minDistToExpectedPt = distToExpectedPt;
cont_copyPts( contBreak1, contB1, true );
cont_copyPts( contBreak2, contB2, true );
}
}
else
{
cont_breakContourEitherSide(cont,contB1,contB2,
intercepts[i].cont1Idx,intercepts[i+1].cont1Idx,true);
float areaOnSmallerSide = MIN( ABS(imodContourArea(contB1)),
ABS(imodContourArea(contB2)) );
if( areaOnSmallerSide > maxAreaOnSmallerSide ) {
maxAreaOnSmallerSide = areaOnSmallerSide;
cont_copyPts( contBreak1, contB1, true );
cont_copyPts( contBreak2, contB2, true );
}
}
}
imodContourDelete(contB1);
imodContourDelete(contB2);
return (psize(contBreak1) > 0 && psize(contBreak2) > 0 );
}
return true;
}
//------------------------
//-- Takes a single contours and breaks it into two contours at the
//--
int cont_breakContourByContour( vector<IcontPtr> &contSegs,
Icont *contO, Icont *contBreak, float minDist )
{
if( isEmpty(contO) )
return 0;
Imod *cont = imodContourDup( contO );
imodContourUnique( cont );
imodContourMakeDirection( cont, IMOD_CONTOUR_CLOCKWISE );
imodContourUnique( contBreak );
cont_addPtsCrude( cont, minDist, true );
int z = getZInt( cont );
const int PT_INSIDE = -2;
for( int p=0; p<psize(cont); p++ )
{
Ipoint *pt = getPt(cont,p);
if( imodPointInsideCont( contBreak, pt ) )
pt->z = PT_INSIDE;
}
cont_breakContByZValue( cont, contSegs, z, true );
imodContourDelete(cont);
return contSegs.size();
}
//------------------------
//-- Takes two contours and joins them together at the two points (one in each contour)
//-- which are closest together. This is used in a branching strategy.
//-- | __ __ | __ __ |
//-- | / \ / \ | / \___/ \ |
//-- | \__/ \__/ | \__/ ^ \__/ | NOTE: will have point in
//-- | | | middle if
//-- | takes two contours | returns single contour | addPtInMiddle is true
//-- | | joined by a double line |
void cont_joinContsAtClosestApproach( Icont *newCont, Icont *cont1Orig,
Icont *cont2Orig, bool addPtInMiddle=true )
{
imodContourDefault( newCont );
if( isEmpty(cont1Orig) || isEmpty(cont2Orig) )
return;
//## FIND CLOSEST POINTS IN BOTH CONTOURS:
Icont *cont1 = imodContourDup(cont1Orig);
Icont *cont2 = imodContourDup(cont2Orig);
imodContourMakeDirection( cont1, IMOD_CONTOUR_CLOCKWISE );
imodContourMakeDirection( cont2, IMOD_CONTOUR_CLOCKWISE );
int closestPtInCont1, closestPtInCont2;
cont_findClosestPts2D( cont1, cont2, &closestPtInCont1, &closestPtInCont2);
Ipoint middlePt = line_findPtFractBetweenPts2D( getPt(cont1,closestPtInCont1),
getPt(cont2,closestPtInCont2), 0.5 );
//## CONSTRUCT JOINED CONTOUR:
cont_copyPts( cont1, newCont, true );
cont_reorderPtsToStartAtIdx( newCont, closestPtInCont1 );
imodPointAppend( newCont, getPt(newCont,0) );
if( addPtInMiddle )
imodPointAppend( newCont, &middlePt );
cont_reorderPtsToStartAtIdx( cont2, closestPtInCont2 );
for(int i=0; i<psize(cont2); i++)
imodPointAppend( newCont, getPt(cont2,i) );
imodPointAppend( newCont, getPt(cont2,0) );
if( addPtInMiddle )
imodPointAppend( newCont, &middlePt );
imodContourDelete(cont1);
imodContourDelete(cont2);
}
//------------------------
//-- Takes a vector of contours and joines them all together in the place
//-- of closest approach to form a single big contour as shown in diagram:
//-- | __ | __ |
//-- | __ __ / \ | __ __ / \ |
//-- | / \ / \ \__/ | / \__/ \__\__/ |
//-- | \__/ \__/ | \__/ \__/ |
//-- | | |
//-- | takes multiple contours | returns single contour |
//-- | | (joined with thin line) |
void cont_joinContsAtClosestApproachVector( Icont *newCont, vector<IcontPtr> conts,
bool addPtInMiddle )
{
imodContourDefault(newCont);
if( (int)conts.size()==0 ) {
wprint( "ERROR: cont_joinContsAtClosestApproach() - empty vector\n" );
return;
}
else if( (int)conts.size()==1 ) {
cont_copyPts( conts[0].cont, newCont, true );
return;
}
else if( (int)conts.size()==2 ) {
cont_joinContsAtClosestApproach(newCont,conts[0].cont,conts[1].cont,addPtInMiddle);
return;
}
else
{
while( (int)conts.size() > 1 )
{
float closestDist = FLOAT_MAX;
int closestContIdx = 1;
for( int i=1; i<(int)conts.size(); i++) {
float minDistThisCont = cont_minDistBetweenContPts2D(conts[0].cont,
conts[i].cont,false);
if( minDistThisCont < closestDist ) {
closestDist = minDistThisCont;
closestContIdx = i;
}
}
cont_joinContsAtClosestApproach(conts[0].cont,conts[0].cont,
conts[closestContIdx].cont,addPtInMiddle);
eraseContour( conts, closestContIdx );
}
cont_copyPts( conts[0].cont, newCont, true );
return;
}
}
//------------------------
//-- Takes two contours, makes them convex and then returns a polygons
//-- representing the intersection (i.e. overlapping area) the two (convex)
//-- contours.
//-- NOTE: The intersection of two convex polygons is ALWAYS a (single)
//-- convex polygon or none at all.
//-- This diagram shows how it works:
//--
//-- | INPUT | make both convex | cont_getIntersectingConvexPolygon
//-- | _________ | _________ |
//-- | ____/__ \ | ____/__ \ | __
//-- | / __|__| | | / | | | | | |
//-- | | / | | > | | | | | > | | |
//-- | | |__|__ | | | | | | | | |
//-- | \____|__| | | \____|__| | | |__|
//-- | cont1 \_________/ | cont1 \_________/ |
//-- | cont2 | cont2 | returns single "intersection"
//-- contour
//-- (cont 1 has changed)
void cont_getIntersectingConvexPolygon( Icont *newCont,
Icont *cont1Orig, Icont *cont2Orig )
{
imodContourDefault(newCont);
if( isEmpty(cont1Orig) || isEmpty(cont2Orig) ) // if either contours are empty:
return; // return empty set
Icont *cont1 = imodContourDup(cont1Orig);
Icont *cont2 = imodContourDup(cont2Orig);
Icont *contInters = imodContourNew(); // stores final overlapping region
// (the intersection of the two contours)
//## PREPARE CONTOURS AND DATA STRUCTURES:
imodContourUnique( cont1 );
imodContourUnique( cont2 );
imodContourMakeDirection( cont1, IMOD_CONTOUR_CLOCKWISE );
imodContourMakeDirection( cont2, IMOD_CONTOUR_CLOCKWISE );
cont_makeConvex( cont1 );
cont_makeConvex( cont2 );
vector<PtConnection> intercepts; // stores a list of pts where the two contours
// intersects, plus the index in both cont1P and
// cont2P where this intersect point was added.
Ipoint intercept; // used in "line_doLinesCrossAndWhere"
// setup "matrix of intersects" - is used to index a list of intercept pts
// (in intercepts) which occur after EACH point in cont1 plus
// the distance to that point (and same for cont2)
vector< vector<IdxAndFloat> > cont1Intercepts( psize(cont1) );
vector< vector<IdxAndFloat> > cont2Intercepts( psize(cont2) );
const int NOT_INTERSECT_PT = -1; // used to mark a pt which is not an intersection pt
// use the z value show that none of these points are intersection points
// (the are all original)
// for intersection points the z value is set to an idx in the intercepts vector.
setZValue( cont1, NOT_INTERSECT_PT );
setZValue( cont2, NOT_INTERSECT_PT );
//## FIND ALL INTERCEPTION POINTS WHERE CONTOURS CROSS AND
//## ADD THEM TO THE LIST OF INTERCEPTS:
for (int i=0; i<psize(cont1);i++) // for each point/line in cont1:
for (int j=0; j<psize(cont2);j++) // for each point/line in cont2:
if(line_doLinesCrossAndWhere( getPt(cont1,i), getPt(cont1,i+1),
getPt(cont2,j), getPt(cont2,j+1), &intercept ) )
{ // if lines cross: add intercept point & calculate distances
PtConnection newInt = PtConnection(intercept);
if( vector_doesElementExistInVector( intercepts, newInt ) )
continue;
intercepts.push_back( newInt );
cont1Intercepts[i].push_back(IdxAndFloat(int(intercepts.size())-1,
line_distBetweenPts2D( &intercept, getPt(cont1,i))));
cont2Intercepts[j].push_back(IdxAndFloat(int(intercepts.size())-1,
line_distBetweenPts2D(&intercept,getPt(cont2,j))));
}
//## IF CONTOURS NEVER CROSS: TEST IF ONE IS INSIDE THE OTHER, AND RETURN APPRIATE VALUE
if( intercepts.size() == 0 ) // if contours don't cross paths at all:
{
if (imodPointInsideCont(cont2, getPt(cont1,0))) //(CASE 1) cont1 inside cont2
{
cont_copyPts(cont1, newCont, true);
imodContourDelete(cont1);
imodContourDelete(cont2);
imodContourDelete(contInters);
return;
}
else if (imodPointInsideCont(cont1, getPt(cont2,0))) //(CASE 2) cont2 inside cont1
{
cont_copyPts(cont2, newCont, true);
imodContourDelete(cont1);
imodContourDelete(cont2);
imodContourDelete(contInters);
return;
}
else { // (CASE 3) don't touch/overlap at all: return empty set
cont_copyPts(contInters, newCont, true);
imodContourDelete(cont1);
imodContourDelete(cont2);
imodContourDelete(contInters);
return;
}
}
//## FOR BOTH CONTOURS: IF MORE THAN ONE INTERCEPT POINTS OCCURS AFTER THE SINGLE POINT:
//## SORT THESE IN ORDER OF THEIR DISTANCE FROM THE POINT
// (this ensure they are added to the contour in the correct order)
for (int i=0; i<int(cont1Intercepts.size()); i++) // for each point in cont1:
if ( int(cont1Intercepts[i].size()) > 1 )
cont1Intercepts[i] = vector_sort( cont1Intercepts[i] );
for (int i=0; i<int(cont2Intercepts.size()); i++) // for each point in cont2:
if ( int(cont2Intercepts[i].size()) > 1 )
cont2Intercepts[i] = vector_sort( cont2Intercepts[i] );
//## FOR BOTH CONTOURS: CREATE A NEW VERSION WHEREBY THE INTERCEPTS POINT ARE ADDED,
//## AND AND MAP THE CONNECTION OF THESE POINTS BETWEEN CONTOURS
Icont *cont1P = imodContourNew(); //|- a version of the contours where the
Icont *cont2P = imodContourNew(); //| intercept pts have been ADDED as extra pts
for (int i=0; i<int(cont1Intercepts.size()); i++) { // for each point in cont1:
imodPointAppend( cont1P, getPt(cont1,i) ); // add it to cont1P
for (int j=0; j<int(cont1Intercepts[i].size()); j++ )
{
int interceptsIdx = cont1Intercepts[i][j].idx;
Ipoint interceptPt = intercepts.at(interceptsIdx).intercept;
interceptPt.z = interceptsIdx;
imodPointAppend( cont1P, &interceptPt ); // add the intercept to cont1P
intercepts.at(interceptsIdx).cont1Idx = psize(cont1P)-1;
}
}
for (int i=0; i<int(cont2Intercepts.size()); i++) { // for each point in cont2:
imodPointAppend( cont2P, getPt(cont2,i) ); // add it to cont2P
for (int j=0; j<int(cont2Intercepts[i].size()); j++ )
{
int interceptsIdx = cont2Intercepts[i][j].idx;
Ipoint interceptPt = intercepts.at(interceptsIdx).intercept;
interceptPt.z = (int)interceptsIdx;
imodPointAppend( cont2P, &interceptPt ); // add the intercept to cont2P
intercepts.at(interceptsIdx).cont2Idx = psize(cont2P)-1;
}
}
//## TRAVERSE NEW CONTOURS AND CONNECTIONS BETWEEN INTERCEPT POINTS
//## TO GENERATE OVERLAPPING CONTOUR:
Ipoint startPoint = intercepts[0].intercept; // the starting point
int currIdx = intercepts[0].cont1Idx;
bool currentlyCont1 = true; // indicates which contour we are traversing
// (if true: then it's cont1, if false: cont2)
while (true)
{
if( currentlyCont1 ) // if we are currently traversing cont1:
{
// add this point (in cont1) to the overlapping polygon
imodPointAppend( contInters, getPt(cont1P,currIdx));
currIdx++; // go to next point in the cont1
if( getPt( cont1P, currIdx)->z != -1 ) { // if this is a intercept point:
int nextIterceptIdx = (int)getPt(cont1P,currIdx)->z;
currIdx = intercepts.at(nextIterceptIdx).cont2Idx;
// ensures a new contour is not started from
// this intercept point (we've alread included it)
currentlyCont1 = false; // flip to cont2 (next iteration it will traverse cont2)
// if we've come back around to our start point: overlapping polygon is complete
if( intercepts.at(nextIterceptIdx).intercept == startPoint )
break;
}
}
else // else (if we are currently traversing cont2):
{
// add this point (in cont2) to the overlapping polygon
imodPointAppend( contInters, getPt(cont2P,currIdx));
currIdx++; // go to next point in cont2
if( getPt( cont2P, currIdx)->z != -1 ) { // if this is a intercept point:
int nextIterceptIdx = (int)getPt(cont2P,currIdx)->z;
currIdx = intercepts.at(nextIterceptIdx).cont1Idx;
// ensures a new contour is not started from
// this intercept point (we've alread included it)
currentlyCont1 = true; // flip back to cont1
// if we've come back around to our start point: overlapping polygon is complete
if( intercepts.at(nextIterceptIdx).intercept == startPoint )
break;
}
}
}
cont_copyPts(contInters, newCont, true);
imodContourDelete(cont1);
imodContourDelete(cont2);
imodContourDelete(contInters);
imodContourDelete(cont1P);
imodContourDelete(cont2P);
return;
}
//------------------------
//-- Takes a single closed contour and breaks it into fragments either side
//-- of any point which has a z value not equal to "zValue" with the point
//-- inclusive. The broken contours are returned in "contSegs" after
//-- setting all their points to zValue and deleting any fragments
//-- with only one point.
int cont_breakContByZValue( Icont *cont, vector<IcontPtr> &contSegs, int zValue,
bool removeOffSegments )
{
contSegs.clear();
if( isEmpty(cont) )
return 0;
int numPts = psize(cont);
float z = (float)zValue;
//## DETERMINE THE FIRST POINT WHICH IS OFF:
int pOffset = 0; // the index of the first point NOT on z
for( ; pOffset<numPts; pOffset++)
if( getPt(cont,pOffset)->z != z )
break;
//## GO THROUGH CONT FROM FIRST OFF POINT AND CREATE SEGMENTS FROM
//## SEQENTIAL SERIES OF ON POINTS:
contSegs.push_back( IcontPtr() );
for (int p=0; p<numPts+1;p++)
{
int i = (p+pOffset) % numPts;
bool currPointOn = getPt(cont,i)->z == z;
if( currPointOn ) // if point is on: add it
{
imodPointAppend( contSegs.back().cont, getPt(cont,i) );
}
else
{
if( removeOffSegments )
{
bool prevPointOn = getPt(cont,i-1)->z == z;
bool nextPointOn = getPt(cont,i+1)->z == z;
if( prevPointOn || nextPointOn ) // if intersection point:
{
if( prevPointOn )
imodPointAppend( contSegs.back().cont, getPt(cont,i) );
contSegs.push_back( IcontPtr() );
if( nextPointOn )
imodPointAppend( contSegs.back().cont, getPt(cont,i) );
}
}
else
{
imodPointAppend( contSegs.back().cont, getPt(cont,i) );
contSegs.push_back( IcontPtr() );
imodPointAppend( contSegs.back().cont, getPt(cont,i) );
}
}
}
//cout << "T1" << endl; flush( cout );
//## CLEAN ALL SEGMENTS, REMOVE ANY WITH ONLY ONE POINT OR LESS,
//## AND ADD AN EXTRA POINT IN THE MIDDLE OF ANY TWO POINT SEGMENTS:
for(int i=(int)contSegs.size()-1; i>=0; i--)
{
Icont *seg = contSegs[i].cont;
imodContourUnique( seg );
setZValue( seg,zValue );
if( psize(seg) <= 1 )
eraseContour( contSegs, i );
else if( psize(seg) == 2 )
{
Ipoint midPt = line_getPtHalfwayBetween( getPt(seg,0), getPt(seg,1) );
imodPointAdd( seg, &midPt, 1 ); // add a point in the middle
}
}
return (contSegs.size());
}
//------------------------
//-- Takes a single open contour and breaks it into fragments either
//-- side of any point which has a z value equal to "zValueToBreak"
//-- with the point not inclusive. The broken contours are returned
//-- in "contSegs" after deleting any fragments with only one point.
int cont_breakOpenContAtZValue( Icont *contOrig, vector<IcontPtr> &contSegs,
int zValueToBreak )
{
contSegs.clear();
if( isEmpty(contOrig) )
return 0;
Icont *cont = imodContourDup(contOrig);
imodContourUnique( cont );
contSegs.push_back( IcontPtr());
for (int i=0; i<psize(cont);i++)
{
// if this is an intersection point: add it, then start a new contour
if( getPt(cont,i)->z == (float)zValueToBreak ) {
contSegs.push_back( IcontPtr());
}
else
imodPointAppend( contSegs.back().cont, getPt(cont,i));
}
for (int i=0; i<(int)contSegs.size(); i++)
{
imodContourUnique( contSegs[i].cont );
if(psize(contSegs[i].cont) <= 1 ) {
eraseContour( contSegs, i );
i--;
}
}
imodContourDelete(cont);
return (contSegs.size());
}
//------------------------
//-- Takes a single open contour and breaks it into fragments either side of
//-- a circle, by marking any points inside the circle and getting rid of them.
int cont_breakContByCircle( Icont *contOrig, vector<IcontPtr> &contSegs,
Ipoint *center, float radius )
{
int radiusSq = (radius*radius);
int numRemovedPoints = 0;
int REMOVE_POINT_Z = -2;
Icont *cont = imodContourDup(contOrig);
for(int p=0; p<psize(cont); p++ )
{
if( getPt(cont,p)->z == center->z )
{
float distSq = line_sqDistBetweenPts2D( center, getPt(cont,p));
if( distSq < radiusSq )
{
getPt(cont, p)->z = REMOVE_POINT_Z;
numRemovedPoints++;
}
}
}
cont_breakOpenContAtZValue( cont, contSegs, REMOVE_POINT_Z );
imodContourDelete(cont);
return (numRemovedPoints);
}
//------------------------
//-- Adds points to closed contours "cont1" and "cont2" anywhere they intersect.
//-- Returns the number of points added, which will usually be equal to 2*the
//-- number of intersections, unless there were already points lying on the exact
//-- location of intersection.
int cont_addPtsAtIntersection( Icont *cont1, Icont *cont2 )
{
int ptsAdded = 0;
Ipoint intercept; // fed into "line_doLinesCrossAndWhere" function
imodContourUnique( cont1 );
imodContourUnique( cont2 );
for (int i=0; i<psize(cont1);i++) // for each point/line in cont1:
{
for (int j=0; j<psize(cont2);j++) // for each point/line in cont2:
{
if( !ptsEqual( getPt(cont1,i), getPt(cont2,j) )
&& !ptsEqual( getPt(cont1,i), getPt(cont2,j+1) )
&& !ptsEqual( getPt(cont1,i+1), getPt(cont2,j) )
&& !ptsEqual( getPt(cont1,i+1), getPt(cont2,j+1) )
&&
line_doLinesCrossAndWhere( getPt(cont1,i), getPt(cont1,i+1),
getPt(cont2,j), getPt(cont2,j+1), &intercept ) )
{
if( ptsAdded > 10000 )
{
wprint("\aERROR: cont_addPtsAtIntersection() - precision error");
imodContourUnique( cont1 );
imodContourUnique( cont2 );
return (ptsAdded);
}
//intercept.x = roundPrec( intercept.x, 0.00001 );
//intercept.y = roundPrec( intercept.y, 0.00001 );
bool ptExistsInCont1 = ptsEqual(&intercept, getPt(cont1,i))
|| ptsEqual(&intercept, getPt(cont1,i+1));
bool ptExistsInCont2 = ptsEqual(&intercept, getPt(cont2,j))
|| ptsEqual(&intercept, getPt(cont2,j+1));
if( !ptExistsInCont1 )
{
imodPointAdd( cont1, &intercept, i+1 ); // add the intercept to cont1
ptsAdded++;
}
if( !ptExistsInCont2 )
{
imodPointAdd( cont2, &intercept, j+1 ); // add the intercept to cont2
ptsAdded++;
}
if( !ptExistsInCont1 || !ptExistsInCont1 ) // if first time intercept found:
j = -1; // go from start of cont2 again
}
}
}
imodContourUnique( cont1 );
imodContourUnique( cont2 );
return (ptsAdded);
}
//------------------------
//-- Takes two contours and breaks both lines apart into pieces/segments
//-- in every place they intersect and returns the number of intercepts found.
int cont_getIntersectingSegments( Icont *cont1Orig, Icont *cont2Orig,
vector<IcontPtr> &cont1Seg,
vector<IcontPtr> &cont2Seg )
{
cont1Seg.clear();
cont2Seg.clear();
if( psize(cont1Orig) <= 1 || psize(cont2Orig) <= 1 ) // if either contours are empty:
return 0; // return empty set
Icont *cont1 = imodContourDup(cont1Orig);
Icont *cont2 = imodContourDup(cont2Orig);
cont_addPtsAtIntersection(cont1, cont2);
int cont1ZVal = getZInt(cont1);
int cont2ZVal = getZInt(cont2);
setZValue( cont1, cont1ZVal );
setZValue( cont2, cont2ZVal );
const float INTERSECT_PT = -2.0f;
for (int i=0; i<psize(cont1);i++) // for each point/line in cont1:
for (int j=0; j<psize(cont2);j++) // for each point/line in cont2:
if( ptsEqual( getPt(cont1,i), getPt(cont2,j) ) )
{
getPt(cont1,i)->z = INTERSECT_PT;
getPt(cont2,j)->z = INTERSECT_PT;
}
cont_breakContByZValue( cont1, cont1Seg, cont1ZVal, false );
cont_breakContByZValue( cont2, cont2Seg, cont2ZVal, false );
imodContourDelete(cont1);
imodContourDelete(cont2);
return ((int)cont1Seg.size()); // returns the number of intercept points
}
//------------------------
//-- Takes two contours and breaks both lines apart into pieces/segments
//-- in every place they intersect and returns the number of intercepts found.
//-- NOTE: If cont1 and cont2 are different slices no conts will be returned.
int cont_getIntersectingSegmentsSafe( Icont *cont1Orig, Icont *cont2Orig,
vector<IcontPtr> &cont1Seg,
vector<IcontPtr> &cont2Seg )
{
cont1Seg.clear();
cont2Seg.clear();
if( psize(cont1Orig) <= 1 || psize(cont2Orig) <= 1 ) // if either contours are empty:
return 0; // return empty set
Icont *cont1 = imodContourDup(cont1Orig);
Icont *cont2 = imodContourDup(cont2Orig);
imodContourUnique( cont1 );
imodContourUnique( cont2 );
int cont1ZVal = getZ(cont1);
int cont2ZVal = getZ(cont2);
//if( cont1ZVal != cont2ZVal)
// wprint( "ERROR: cont_getIntersectingSegments() - different z vals\n" );
vector<PtConnection> intercepts; // stores a list of points where the two contours
// intersect, plus the index in both cont1P and
// cont2P where this intersect point was added.
Ipoint intercept; // fed into "line_doLinesCrossAndWhere" function.
// create "matrix of intersects" - is used to index a list of intercept
// points (in intercepts) which occur after EACH point in cont1 plus
// the distance to that point (and same for cont2)
vector< vector<IdxAndFloat> > cont1Intercepts( psize(cont1));
vector< vector<IdxAndFloat> > cont2Intercepts( psize(cont2));
const int INTERSECT_PT = -2;
//## FIND ALL INTERCEPTION POINTS WHERE CONTOURS CROSS AND
//## ADD THEM TO THE LIST OF INTERCEPTS:
for (int i=0; i<psize(cont1);i++) // for each point/line in cont1:
for (int j=0; j<psize(cont2);j++) // for each point/line in cont2:
if( line_doLinesCrossAndWhere( getPt(cont1,i), getPt(cont1,i+1),
getPt(cont2,j), getPt(cont2,j+1), &intercept ) )
// if lines cross: add intercept point & calculate distances
{
PtConnection newInt = PtConnection(intercept);
if( vector_doesElementExistInVector( intercepts, newInt )
|| ptsEqual( getPt(cont1,i+1), &intercept)
|| ptsEqual( getPt(cont2,j+1), &intercept ) )
continue;
intercepts.push_back( newInt );
cont1Intercepts[i].push_back(IdxAndFloat(int(intercepts.size())-1,
line_distBetweenPts2D( &intercept, getPt(cont1,i))));
cont2Intercepts[j].push_back(IdxAndFloat(int(intercepts.size())-1,
line_distBetweenPts2D( &intercept, getPt(cont2,j))));
}
//## IF CONTOURS NEVER CROSS: RETURN ZERO
if( intercepts.size() == 0 )
{
imodContourDelete(cont1);
imodContourDelete(cont2);
return 0;
}
//## FOR BOTH CONTOURS: IF MORE THAN ONE INTERCEPT POINT AFTER A SINGLE POINT:
//## SORT THESE IN ORDER OF THEIR DISTANCE FROM THE POINT
for (int i=0; i<int(cont1Intercepts.size()); i++) // for each point in cont1:
if ( (int)cont1Intercepts[i].size() > 1 )
cont1Intercepts[i] = vector_sort( cont1Intercepts[i] );
for (int i=0; i<int(cont2Intercepts.size()); i++) // for each point in cont2:
if ( (int)cont2Intercepts[i].size() > 1 )
cont2Intercepts[i] = vector_sort( cont2Intercepts[i] );
//## FOR BOTH CONTOURS: CREATE A NEW VERSION WHEREBY THE INTERCEPTS POINT ARE ADDED,
//## AND MAP THE CONNECTION OF THESE POINTS BETWEEN CONTOURS
Icont *cont1P = imodContourNew(); //|- stores a version of the contours where the
Icont *cont2P = imodContourNew(); //| intercept points have been ADDED as extra pts
for (int i=0; i<int(cont1Intercepts.size()); i++) { // for each point in cont1:
imodPointAppend( cont1P, getPt(cont1,i)); // add it to cont1P
for (int j=0; j<(int)cont1Intercepts[i].size(); j++ )
{
int interceptsIdx = cont1Intercepts[i][j].idx;
Ipoint interceptPt = intercepts.at(interceptsIdx).intercept;
interceptPt.z = (float)INTERSECT_PT;
imodPointAppend( cont1P, &interceptPt ); // add the intercept to cont1P
intercepts.at(interceptsIdx).cont1Idx = psize(cont1P)-1;
}
}
for (int i=0; i<int(cont2Intercepts.size()); i++) { // for each point in cont2:
imodPointAppend( cont2P, getPt(cont2,i)); // add it to cont2P
for (int j=0; j<(int)cont2Intercepts[i].size(); j++ )
{
int interceptsIdx = cont2Intercepts[i][j].idx;
Ipoint interceptPt = intercepts.at(interceptsIdx).intercept;
interceptPt.z = (float)INTERSECT_PT;
imodPointAppend( cont2P, &interceptPt ); // add the intercept to cont2P
intercepts.at(interceptsIdx).cont2Idx = psize(cont2P)-1;
}
}
//## CLEAN LINE SEGMENTS AND DELETE ANY EMPTY ONES:
cont_breakContByZValue( cont1P, cont1Seg, cont1ZVal, false );
cont_breakContByZValue( cont2P, cont2Seg, cont2ZVal, false );
imodContourDelete(cont1);
imodContourDelete(cont2);
imodContourDelete(cont1P);
imodContourDelete(cont2P);
return ((int)cont1Seg.size()); // returns the number of intercept points
}
//------------------------
//-- Returns a vector of polygons representing the intersection
//-- (i.e. overlapping area) of two contours.
//-- WARNING: This function give correct result most of the time,
//-- but sometimes gives sligtly wrong polygons.
//--
//-- | INPUT | cont_getIntersectingPolygons |
//-- | _________ | |
//-- | ____/__ \ | __ |
//-- | / __|__| | | |__| |
//-- | | / | | > | |
//-- | | |__|__ | | __ |
//-- | \____|__| | | |__| |
//-- | cont1 \_________/ | |
//-- | cont2 | returns 2 contour |
int cont_getIntersectingPolygons(vector<IcontPtr> &finConts, Icont *cont1, Icont *cont2)
{
finConts.clear(); // will store the final overlapping regions
// (the intersection of the two contours)
if( isEmpty(cont1) || isEmpty(cont2) ) // if either contours is empty:
return (finConts.size()); // return empty set
//## BREAK CONTOURS INTO INTERSECTING SEGMENTS:
vector<IcontPtr> cont1Seg;
vector<IcontPtr> cont2Seg;
int numIntersectPts =
cont_getIntersectingSegmentsSafe( cont1, cont2, cont1Seg, cont2Seg );
//## IF CONTOURS NEVER CROSS: TEST IF ONE IS INSIDE THE OTHER,
//## AND RETURN APPROPRIATE VALUE
if( numIntersectPts == 0 ) // if contours don't cross paths at all:
{
if (imodPointInsideCont(cont2, getPt(cont1,0))) //(CASE 1) cont1 inside cont2
finConts.push_back( IcontPtr(cont1));
else if (imodPointInsideCont(cont1, getPt(cont2,0))) //(CASE 2) cont2 inside cont1
finConts.push_back( IcontPtr(cont2));
else //(CASE 3) don't touch/overlap at all
finConts.empty();
//deleteContours( cont1Seg );
//deleteContours( cont2Seg );
return (finConts.size());
}
//## ERASE LINE SEGMENTS WHICH ARE NOT PART OF INTERSECTING AREA:
if( numIntersectPts%2 ==1 ) {
wprint( "WARNING: ODD NUMBER INTERSECTIONS !!!!\n" );
for (int i=(int)cont1Seg.size()-1; i>=0; i--) // for all intercept points
if(imodPointInsideCont( cont2, getPt( cont1Seg[i].cont,1) ) )
eraseContour(cont1Seg, i);
for (int i=(int)cont2Seg.size()-1; i>=0; i--) // for all intercept points
if(imodPointInsideCont( cont1, getPt( cont2Seg[0].cont,1) ) )
eraseContour(cont2Seg, i);
}
else // delete every second segment (falling outside)
{
bool firstSegmentC1Inside = imodPointInsideCont(cont2, getPt(cont1Seg[0].cont,1));
int offset = (firstSegmentC1Inside) ? 1 : 0 ;
for(int i=(int)cont1Seg.size()-1; i>=0; i--)
if( i>=offset && i%2 == offset )
eraseContour( cont1Seg, i );
bool firstSegmentC2Inside = imodPointInsideCont(cont1, getPt(cont2Seg[0].cont,1));
offset = (firstSegmentC2Inside) ? 1 : 0 ;
for(int i=(int)cont2Seg.size()-1; i>=0; i--)
if( i>=offset && i%2 == offset )
eraseContour( cont2Seg, i );
}
//## JOIN CONT1 INTERSECTING SEGMENTS WITH TOUCHING INTERSECTING SEGMENTS IN CONT2
for (int i=0; i<(int)cont1Seg.size(); i++) // for all intercept points
{
for (int j=0; j<(int)cont2Seg.size(); j++) // for all intercept points
{
if( ptsEqual(getLastPt(cont1Seg[i].cont),getLastPt(cont2Seg[j].cont))
|| ptsEqual(getLastPt(cont1Seg[i].cont),getFirstPt(cont2Seg[j].cont)))
{
Icont *newIntersectingCont = imodContourNew();
cont_concat( newIntersectingCont, cont1Seg[i].cont, cont2Seg[j].cont, true );
finConts.push_back( IcontPtr(newIntersectingCont));
imodContourDelete(newIntersectingCont);
}
}
}
//## JOIN ANY REMAINING/TOUCHING INTERSECTING SEGMENTS:
Icont *tempCont = imodContourNew();
for (unsigned i=0; i<finConts.size(); i++)
{
for (unsigned j=i+1; j<finConts.size(); j++)
{
if( ptsEqual( getLastPt( finConts[i].cont ), getLastPt( finConts[j].cont ) )
|| ptsEqual( getLastPt( finConts[i].cont ), getFirstPt( finConts[j].cont ) ) )
{
cont_copyPts( finConts[i].cont, tempCont, true );
cont_concat( finConts[i].cont, tempCont, finConts[j].cont, true );
eraseContour( finConts, j );
i--;
break;
}
}
}
imodContourDelete( tempCont );
deleteContours( cont1Seg );
deleteContours( cont2Seg );
return (finConts.size());
}
//------------------------
//-- Returns a vector of polygons representing the union (i.e. combined area)
//-- of two contours.
//-- NOTE: If the two contours intersect in > 2 places then multiple
//-- contours will be returned - one of them an "outer" contour
//-- and the other ones will actual be holes in the outer contour.
//-- (see "cont_getOuterUnionPolygon" for a diagram)
//--
int cont_getUnionPolygons( vector<IcontPtr> &finConts, Icont *cont1, Icont *cont2 )
{
finConts.clear(); // will store the final combined regions
// (the union of the two contours)
if( isEmpty(cont1) || isEmpty(cont2) ) // if either contours is empty:
return (finConts.size()); // return empty set
//## BREAK CONTOURS INTO INTERSECTING SEGMENTS:
vector<IcontPtr> cont1Seg; // list of segments in cont1
vector<IcontPtr> cont2Seg; // list of segments in cont2
int numIntersectPts = cont_getIntersectingSegmentsSafe(cont1, cont2, cont1Seg, cont2Seg);
//wprint( "cont1Seg=%d cont1Seg=%d\n", cont1Seg.size(), cont2Seg.size()); //%%%%%%
if( numIntersectPts > 1000 )
{
wprint("WARNING: Floating point problem - too many intersections\n");
return 0;
}
//## IF CONTOURS NEVER CROSS: TEST IF ONE IS INSIDE THE OTHER,
//## AND RETURN APPROPRIATE VALUE
if( numIntersectPts == 0 ) // if contours don't cross paths at all:
{
if (imodPointInsideCont(cont2, getPt(cont1,0))) //(CASE 1) cont1 inside cont2
finConts.push_back( IcontPtr(cont1));
else if (imodPointInsideCont(cont1, getPt(cont2,0))) //(CASE 2) cont2 inside cont1
finConts.push_back( IcontPtr(cont2));
else // (CASE 3) don't touch/overlap at all
finConts.empty();
deleteContours(cont1Seg);
deleteContours(cont2Seg);
return (finConts.size());
}
//## ERASE LINE SEGMENTS WHICH ARE NOT PART OF INTERSECTING AREA:
if( numIntersectPts%2 == 1 )
{
wprint( "WARNING: ODD NUMBER INTERSECTIONS !!!! \n");
for (int i=(int)cont1Seg.size()-1; i>=0; i--) // for all intercept points
if(imodPointInsideCont( cont2, getPt(cont1Seg[i].cont,1) ) )
eraseContour( cont1Seg, i );
for (int i=(int)cont2Seg.size()-1; i>=0; i--) // for all intercept points
if(imodPointInsideCont( cont1, getPt(cont2Seg[0].cont,1) ) )
eraseContour( cont2Seg, i );
}
else // delete every second segment (falling inside)
{
bool firstSegmentC1Inside = imodPointInsideCont(cont2, getPt(cont1Seg[0].cont,1));
int offset = (firstSegmentC1Inside) ? 0 : 1 ;
for(int i=(int)cont1Seg.size()-1; i>=0; i--)
if( i>=offset && i%2 == offset )
eraseContour( cont1Seg, i );
bool firstSegmentC2Inside = imodPointInsideCont(cont1, getPt(cont2Seg[0].cont,1));
offset = (firstSegmentC2Inside) ? 0 : 1 ;
for(int i=(int)cont2Seg.size()-1; i>=0; i--)
if( i>=offset && i%2 == offset )
eraseContour( cont2Seg, i );
}
//## JOIN CONT1 INTERSECTING SEGMENTS WITH TOUCHING INTERSECTING SEGMENTS IN CONT2
Icont *newIntersectingCont = imodContourNew();
for (int i=0; i<(int)cont1Seg.size(); i++) // for all intercept points
{
for (int j=0; j<(int)cont2Seg.size(); j++) // for all intercept points
{
if( ptsEqual(getLastPt(cont1Seg[i].cont), getLastPt(cont2Seg[j].cont))
|| ptsEqual(getLastPt(cont1Seg[i].cont), getFirstPt(cont2Seg[j].cont)))
{
cont_concat( newIntersectingCont, cont1Seg[i].cont, cont2Seg[j].cont, true );
finConts.push_back( IcontPtr(newIntersectingCont));
}
}
}
imodContourDelete(newIntersectingCont);
//## JOIN ANY REMAINING/TOUCHING INTERSECTING SEGMENTS:
Icont *tempCont = imodContourNew();
for (int i=0; i<(int)finConts.size(); i++)
{
for (int j=i+1; j<(int)finConts.size(); j++)
{
if( ptsEqual( getLastPt( finConts[i].cont ), getLastPt( finConts[j].cont ) )
|| ptsEqual( getLastPt( finConts[i].cont ), getFirstPt( finConts[j].cont ) ) )
{
cont_copyPts( finConts[i].cont, tempCont, true );
cont_concat( finConts[i].cont, tempCont, finConts[j].cont, true );
eraseContour( finConts, j );
i--;
break;
}
}
}
imodContourDelete( tempCont );
deleteContours( cont1Seg );
deleteContours( cont2Seg );
return (finConts.size());
}
//------------------------
//-- Returns true if c1 has a smaller area than c2
bool cont_smallerArea( IcontPtr c1, IcontPtr c2 ) {
return ( ABS(imodContourArea( c1.cont )) < ABS(imodContourArea( c2.cont )));
}
//------------------------
//-- Returns a single polygon "newCont" representing the outer union
//-- (i.e. combined outer area) of the two contours "cont1O" and "cont2O".
//-- Returns true if successful, or false if could not find outer union successfully.
//-- NOTE: If the polygons don't overlap an empty polygon will be returned.
//-- NOTE: If the polygons intersect at > 2 places then the union area
//-- will consist of one OUTER polygon, and several smaller polygons
//-- inside this representing HOLES and this function will only
//-- return the largest (outer) polygon.
//--
//-- This diagram shows the difference between cont_getUnionPolygons
//-- AND cont_getOuterUnionPolygon
//--
//-- | INPUT | cont_getUnionPolygons | cont_getOuterUnionPolygon |
//-- | _________ | _________ | _________ |
//-- | ____/__ \ | ____/__ \ | ____/ \ |
//-- | / __|__| | | / |__| | | / | |
//-- | | / | | > | | | > | | | |
//-- | | |__|__ | | | __ | | | | |
//-- | \____|__| | | \____|__| | | \____ | |
//-- | cont1 \_________/ | \_________/ | \_________/ |
//-- | cont2 | returns | returns only |
//-- THREE CONTOURS OUTER CONTOUR
//-- (including 2 holes)
bool cont_getOuterUnionPolygon( Icont *newCont, Icont *cont1O, Icont *cont2O )
{
vector<IcontPtr> joinedConts;
cont_getUnionPolygons( joinedConts, cont1O, cont2O );
if ( joinedConts.empty() ) {
cont_getUnionPolygons( joinedConts, cont2O, cont1O );
}
if ( joinedConts.empty() ) { // if contours don't touch: return empty contour
int numCrosses = cont_numTimesCountsCross( cont1O, cont2O );
//wprint("%d crosses detected", numCrosses); //%%%%%%%%%%%%%%%
wprint( "Could not join contours\n" );
return (false);
}
else if ( (int)joinedConts.size() == 1 ) { // if is only one union polygon: return it
//wprint("copying\n");
cont_copyPts( joinedConts[0].cont, newCont, true );
}
else // if there are > 1 union polygons: return the biggest one
{
float maxArea=0;
int maxIdx=0;
for(int i=0; i<(int)joinedConts.size(); i++) {
float areaCont = imodContourArea(joinedConts[i].cont);
if( areaCont > maxArea ) {
maxArea = areaCont;
maxIdx = i;
}
}
if( maxArea < imodContourArea(cont1O) || maxArea < imodContourArea(cont2O) )
{
wprint( "\aERROR: try smoothing contours more before joining\n" );
return (false);
}
IcontPtr contWithBiggestArea = joinedConts[maxIdx];
cont_copyPts( contWithBiggestArea.cont, newCont, true );
}
deleteContours( joinedConts );
return (true);
}
//------------------------
//-- Returns the point a given distance along the path of the contour
//-- from the specified start point
Ipoint cont_getPtDistAlongLength( Icont *cont, float dist, bool closed, int startPt )
{
Ipoint returnPt;
if( isEmpty(cont) )
return (returnPt);
float contLength = imodContourLength( cont, closed ); // length of the contour
if( contLength == 0 || dist == 0 )
return (*getFirstPt(cont));
if( dist == contLength )
return(*getLastPt(cont));
dist = fMod( dist, contLength );
float distCurrPt = 0; // distance along contour to the current point
float distNextPt = 0; // distance along contour to the next point
for(int p=0; p<psize(cont); p++)
{
Ipoint *ptCurr = getPt( cont, startPt+ p );
Ipoint *ptNext = getPt( cont, startPt+(p+1) );
distNextPt += imodPointDistance( ptCurr, ptNext );
if( dist <= distNextPt )
{
float fractAlongSeg = fDiv((dist-distCurrPt),(distNextPt-distCurrPt));
returnPt = line_findPtFractBetweenPts2D( ptCurr, ptNext, fractAlongSeg );
return (returnPt);
}
distCurrPt = distNextPt;
}
wprint("\aERROR: cont_getPtDistAlongLength()\n");
return (*getLastPt(cont));
}
//------------------------
//-- Returns the point a given percentage distance along the path of the contour
//-- from the specified start point
Ipoint cont_getPtFractAlongLength( Icont *cont, float fract, bool closed, int startPt )
{
float contLength = imodContourLength( cont, closed ); // length of the contour
float distAlong = fract * contLength; // corresponding distance along contour
return( cont_getPtDistAlongLength(cont,distAlong,closed,startPt) );
}
//------------------------
//-- Returns the length up until each point in the contour as a fraction of the
//-- contour's total length.
//-- Note that the first value in the vector will always be 0,
//-- and if the contour is OPEN, the last value will be 1.
vector<float> cont_getFractPtsAlongLength( Icont *cont, bool closed, int startPt )
{
float contLength = imodContourLength( cont, closed );
vector<float> fractAlongLengthV;
fractAlongLengthV.push_back(0);
float cumLengthToPt = 0;
for(int p=1; p<psize(cont); p++)
{
Ipoint *ptCurr = getPt( cont, startPt+ p );
Ipoint *ptPrev = getPt( cont, startPt+(p-1) );
cumLengthToPt += imodPointDistance( ptCurr, ptPrev );
fractAlongLengthV.push_back( fDiv(cumLengthToPt,contLength) );
}
return fractAlongLengthV;
}
//------------------------
//-- Creates a new contour by adding points along 'cont' where each added point is
//-- "fractAlongLength" percent along the total length of 'cont' from the given
//-- "startPt and returns the final number of points in the new contour.
//-- If "keepExistingPts" is true, all existing points in the contour are preserved.
//-- Note that "fractAlongLength" must be in ascending order and all values between 0 and 1
int cont_addPtsFractsAlongLength( Icont *cont, Icont *contNew,
vector<float> fractsAlongLen,
bool closed, bool keepExistingPts, int startPt )
{
imodContourDefault(contNew);
if( (int)fractsAlongLen.size() == 0 )
{
if( keepExistingPts )
contNew = imodContourDup(contNew);
return psize(contNew);
}
float contLength = imodContourLength( cont, closed ); // length of the contour
float cumLengthToNextPt = 0; // total length up to the next point
float fractCurrPt = 0; // fraction of contour length up to the current point
float fractNextPt = 0; // fraction of contour length up to the next point
bool allFractsFound = false; // becomes true when all "fractsAlongLen" are added
float i = 0; // current index within "fractsAlongLen"
for(int p=0; p<psize(cont); p++) // for each point:
{
Ipoint *ptCurr = getPt( cont, startPt+ p );
Ipoint *ptNext = getPt( cont, startPt+(p+1) );
if( keepExistingPts )
imodPointAppend( contNew, ptCurr );
cumLengthToNextPt += imodPointDistance( ptCurr, ptNext );
fractNextPt = fDiv( cumLengthToNextPt, contLength );
while( !allFractsFound )
{
if( fractsAlongLen[i] <= fractNextPt )
{
float fractAlongSeg = fDiv( (fractsAlongLen[i] - fractCurrPt),
(fractNextPt - fractCurrPt) );
Ipoint newPt = line_findPtFractBetweenPts2D( ptCurr, ptNext, fractAlongSeg );
imodPointAppend( contNew, &newPt );
i++;
if( i >= (int)fractsAlongLen.size() )
{
allFractsFound = true;
if( !keepExistingPts )
return psize(contNew);
}
}
else
{
break;
}
}
fractCurrPt = fractNextPt;
}
if( !allFractsFound ) // should not happen... but sometimes does
{
for( ; i<(int)fractsAlongLen.size() && fractsAlongLen[i] >= 1.0; i++ )
imodPointAppend( contNew, getLastPt(cont) );
//wprint("\aERROR: cont_addPtsFractsAlongLength()\n");
}
return psize(contNew);
}
//------------------------
//-- Calculate a mass of information about a contour regarding point size.
//-- Is mostly designed for "tubes of varying thickness" - contours with
//-- different point sizes. To get all results in pixels, set "pixelSize" to 1,
//-- otherwise all values will be scaled by "pixelSize".
//-- NOTE: the values "openVol", "fullVol" and "surfaceArea" are estimates only
//-- they don't take into account the fact tubes can fold
//-- over themselves.
void cont_calcPtsizeInfo( Iobj *obj, Icont *cont, const float zScale, float pixelSize,
float &openLength, float &fullLength,
float &avgR, float &firstR, float &lastR,
float &minR, float &maxR, float &minMidR,
float &openVol, float &fullVol, float &surfaceArea )
{
int numPts = psize(cont);
openLength = 0; // length of the contour
fullLength = 0; // contour length + the radius of the first and last pt
firstR = 0; // radius of the first point
lastR = 0; // radius of the last point
avgR = 0; // average radius along the open portion of the contour
minR = FLOAT_MAX; // min point radius
maxR = 0; // max point radius
minMidR = FLOAT_MAX; // min point raidus excluding the first and last points
openVol = 0; // a crude estimate of the tube's volume in pixels cubed
fullVol = 0; // tubes volume estimate + hemisphere at start and end
surfaceArea = 0; // crude surface area for hemisphere capped tube
if( numPts == 0 )
return;
firstR = imodPointGetSize(obj, cont, 0);
lastR = imodPointGetSize(obj, cont, numPts-1);
float weightedR = 0;
Ipoint scale;
scale.x = 1.0;
scale.y = 1.0;
scale.z = zScale;
//## FIND MAXIMUM AND MINIMUM POINT SIZE:
float totalPtSize = 0; // only used if contour has no length
for(int p=0; p<psize(cont); p++)
{
totalPtSize += imodPointGetSize(obj,cont,p);
updateMin( minR, imodPointGetSize(obj,cont,p) );
updateMax( maxR, imodPointGetSize(obj,cont,p) );
}
for(int p=1; p<(psize(cont)-1); p++) // for the second point to the second last:
updateMin( minMidR, imodPointGetSize(obj,cont,p) );
if (minMidR==FLOAT_MAX) // can occur if only two points
minMidR = MIN(firstR,lastR);
//## CALCULATE LENGTH OF EACH SEGMENT AND USE IT TO MAKE A VOLUME
//## AND SURFACE AREA ESTIMATE:
for(int p=0; p<psize(cont)-1; p++)
{
Ipoint *pt1 = getPt(cont,p);
Ipoint *pt2 = getPt(cont,p+1);
float pt1Size = imodPointGetSize(obj,cont,p);
float pt2Size = imodPointGetSize(obj,cont,p+1);
float segmentLen = imodPoint3DScaleDistance( pt1, pt2, &scale );
float avgRForTwoPts = ( pt1Size + pt2Size) / 2.0;
float volBetweenTwoPs = geom_volumeConicalFrustrum( pt1Size, pt2Size, segmentLen );
openLength += segmentLen;
openVol += volBetweenTwoPs;
avgR += pt1Size;
weightedR += segmentLen * avgRForTwoPts;
surfaceArea += geom_surfaceAreaConicalFrustrum( pt1Size, pt2Size, segmentLen);
}
if(openLength==0)
avgR = totalPtSize / numPts;
else
avgR = weightedR / openLength;
fullLength = openLength + (firstR + lastR);
fullVol = openVol + 0.5*geom_volumeSphere(firstR) + 0.5*geom_volumeSphere(lastR);
surfaceArea = surfaceArea + 0.5*geom_surfaceAreaSphere(firstR)
+ 0.5*geom_surfaceAreaSphere(lastR);
//## CONVERT VALUES TO APPROPRIATE UNITS
if(pixelSize != 1)
{
openLength *= pixelSize;
fullLength *= pixelSize;
avgR *= pixelSize;
firstR *= pixelSize;
lastR *= pixelSize;
minR *= pixelSize;
maxR *= pixelSize;
minMidR *= pixelSize;
openVol *= CUBE(pixelSize);
fullVol *= CUBE(pixelSize);
surfaceArea *= SQ(pixelSize);
}
}
Acknowledgements: Brad Marsh and David Mastronarde for introducing me to IMOD and encouraging me to add to it's development |