symb.cpp
#include "symb.h"
#include <iomanip>
Interval EuclDist(const IVector & x, const IVector & y){
Interval res(0.0);
if (x.dimension()!=y.dimension()) return 1e100;
for (int i=0;i<x.dimension();++i) res+=sqr(x[i]-y[i]);
return sqrt(res);
}
string ColorString(int i){
if (i>=5) i=i%5;
string s;
if (i==0) s="1 0 0";
if (i==1) s="0 0 1";
if (i==2) s="1 0 1";
if (i==3) s="1 1 0";
if (i==4) s="0 1 1";
return s;
}
IVector GetRange(const Polygon & P) {
IVector Range;
for (int k=0;k<P.size();k++){
if (k==0) Range=P[k];
else Range=Hull(Range,P[k]);
}
return Range;
}
void PolygonSetSave(const PolygonSet & Polygons, string FileName, int prec){
ofstream OutFile(FileName.c_str()) ;
cout << "Saving " << FileName << endl;
OutFile << setprecision(prec);
OutFile << "Num " << Polygons.size() << endl;
if (Polygons.size()>0 && Polygons[0].size()>0) OutFile << "Dim " << Polygons[0][0].dimension() << endl;
for (int i=0;i<Polygons.size();i++){
OutFile << "P " << i << " " << Polygons[i].size() << endl;
for (int j=0;j<Polygons[i].size();++j){
for (int d=0;d<Polygons[i][j].dimension();d++) OutFile << (d>0?" ":"") << setw(prec+2) << setiosflags(ios::fixed) << Polygons[i][j][d].leftBound();
OutFile << endl;
}
}
OutFile.close();
}
void PolygonSetLoad(string FileName, PolygonSet & Polygons){
int dim=0;
Polygons.clear();
ifstream InFile(FileName.c_str());
if (!InFile){ cerr << "File " << FileName << " not found." << endl; return; }
int PolyNum=0;
char ch;
while (InFile.get(ch)){
string WS;
if (ch=='#'){
getline(InFile, WS);
continue;
}
InFile.putback(ch);
InFile >> WS;
if (WS=="Num") InFile >> PolyNum;
if (WS=="Dim") InFile >> dim;
if (WS=="P"){
Polygon Poly;
int pos;
InFile >> pos;
if (pos<0 || pos>=PolyNum) {cerr << "Read Error\n"; return;}
int PNum;
InFile >> PNum;
for (int k=0;k<PNum;k++) {
IVector P(dim);
for (int d=0;d<dim;d++){
double tmp;
InFile >> tmp;
P[d]=Interval(tmp);
}
Poly.push_back(P);
}
Polygons.push_back(Poly);
}
getline(InFile, WS);
}
InFile.close();
}
PolygonSet PolygonSetSymmetric(const PolygonSet & Polygons){
PolygonSet PolygonsSymm=Polygons;
for (int i=0;i<Polygons.size();i++){
Polygon PolySymm;
for (int j=0;j<Polygons[i].size();++j){
PolySymm.push_back(-Polygons[i][j]);
}
PolygonsSymm.push_back(PolySymm);
}
return PolygonsSymm;
}
IVector Segment::GetRange() const {
IVector Range;
for (int k=0;k<size();k++){
if (k==0) Range=P[k];
else Range=Hull(Range,P[k]);
}
return Range;
}
double SegmentDist(const Segment & S1, const Segment & S2){
double dist=1e10;
for (int i=0;i<S1.size();++i){
for (int j=0;j<S2.size();++j){
Interval d=EuclDist(S1.P[i],S2.P[i]);
// if (d<=dist) dist=d;
if (d.leftBound()<=dist) dist=d.leftBound();
}
}
return dist;
}
IVector Quadr::GetRange(){
IVector Range;
int PNum=0;
for (int j=0;j<SegmentNum();j++){
for (int k=0;k<S[j].PointNum();k++){
if (PNum==0) Range=S[j].P[k];
else Range=Hull(Range,S[j].P[k]);
PNum++;
}
}
return Range;
}
void QuadrSet::Read(string FileName){
int dim=0;
Q.clear();
ifstream InFile(FileName.c_str());
if (!InFile){ cerr << "File " << FileName << " not found." << endl; return; }
cout << "Reading " << FileName << "...\n";
int QNum=0;
char ch;
while (InFile.get(ch)){
string WS;
if (ch=='#'){
getline(InFile, WS);
continue;
}
InFile.putback(ch);
InFile >> WS;
if (WS=="QuadrNum") InFile >> QNum;
if (WS=="Dim") InFile >> dim;
if (WS=="Q"){
Quadr QTmp;
int pos;
InFile >> pos;
if (pos<0 || pos>=QNum) {cerr << "Read Error\n"; return;}
int SNum;
InFile >> SNum;
for (int j=0;j<SNum;j++){
Segment STmp;
int PNum;
InFile >> PNum;
for (int k=0;k<PNum;k++) {
IVector P(dim);
for (int d=0;d<dim;d++){
double tmp;
InFile >> tmp;
P[d]=Interval(tmp);
}
STmp.AddPoint(P);
}
QTmp.AddSegment(STmp);
}
AddQuadr(QTmp);
}
getline(InFile, WS);
}
InFile.close();
Write();
}
void QuadrSet::Write(int prec){
cout << setprecision(prec);
for (int i=0;i<QuadrNum();i++){
cout << "Q " << i << " " << Q[i].SegmentNum() << endl;
for (int j=0;j<Q[i].SegmentNum();j++){
cout << Q[i].S[j].PointNum() << " ";
for (int k=0;k<Q[i].S[j].PointNum();k++){
for (int d=0;d<Q[i].S[j].P[k].dimension();d++) cout << Q[i].S[j].P[k][d] << " ";
}
cout << endl;
}
}
}
void QuadrSet::Save(string FileName, int prec){
ofstream OutFile(FileName.c_str()) ;
cout << "Saving " << FileName << endl;
OutFile << setprecision(prec);
OutFile << "QuadrNum " << QuadrNum() << endl;
if (QuadrNum()>0) OutFile << "Dim " << Q[0].S[0].P[0].dimension() << endl;
for (int i=0;i<QuadrNum();i++){
OutFile << "Q " << i << " " << Q[i].SegmentNum() << endl;
for (int j=0;j<Q[i].SegmentNum();j++){
OutFile << Q[i].S[j].PointNum();
for (int k=0;k<Q[i].S[j].PointNum();k++){
for (int d=0;d<Q[i].S[j].P[k].dimension();d++) OutFile << " " << setw(prec+2) << setiosflags(ios::fixed) << Q[i].S[j].P[k][d].leftBound();
}
OutFile << endl;
}
}
OutFile.close();
}
IVector QuadrSet::GetRange(){
IVector Range;
int QNum=0;
for (int i=0;i<QuadrNum();i++){
if (i==0) Range=Q[i].GetRange();
else Range=Hull(Range,Q[i].GetRange());
}
return Range;
}
PolygonSet QuadrSet::GetPolygonSet(bool MakeSymmetric){
PolygonSet Polygons;
for (int i=0;i<Q.size();++i){
vector <IVector> Poly;
vector <IVector> PolySymm;
for (int j=0;j<Q[i].SegmentNum();j++){
for (int k=0;k<Q[i].S[j].PointNum();k++){
if (j>0 && k==0) continue;
Poly.push_back(Q[i].S[j].P[k]);
PolySymm.push_back(-Q[i].S[j].P[k]);
}
}
Polygons.push_back(Poly);
if (MakeSymmetric) Polygons.push_back(PolySymm);
}
return Polygons;
}
void QuadrSet::PrintEPS(string FileName, EPSOptions & EPSOpt, int XDisp, int YDisp, bool DoSymmetric, bool PlotLegend, bool UseSingleColor, bool fill){
cout << "QuadrSet::PrintEPS" << endl;
EPSOpt.EPSFileCreator="cc (Z. Galias)";
ostringstream oss;
oss << "QuadrSet";
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
bool PlotQuadr=true;
if (PlotQuadr){
for (int i=0;i<QuadrNum();i++){
IVector Range=Q[i].GetRange();
if (Inf(Range(XDisp))>EPSOpt.XMax || Sup(Range(XDisp))<EPSOpt.XMin || Inf(Range(YDisp))>EPSOpt.YMax || Sup(Range(YDisp))<EPSOpt.YMin) continue;
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
for (int j=0;j<Q[i].SegmentNum();j++){
if (j==0) aEPS.OutFile << "n " << aEPS.FunX(Mid(Q[i].S[j].P[0](XDisp))) << " " << aEPS.FunY(Mid(Q[i].S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q[i].S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(Mid(Q[i].S[j].P[k](XDisp))) << " " << aEPS.FunY(Mid(Q[i].S[j].P[k](YDisp))) << " l ";
}
}
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
if (DoSymmetric){
for (int j=0;j<Q[i].SegmentNum();j++){
if (j==0) aEPS.OutFile << "n " << aEPS.FunX(-Mid(Q[i].S[j].P[0](XDisp))) << " " << aEPS.FunY(-Mid(Q[i].S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q[i].S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(-Mid(Q[i].S[j].P[k](XDisp))) << " " << aEPS.FunY(-Mid(Q[i].S[j].P[k](YDisp))) << " l ";
}
}
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
}
for (int j=0;j<Q[i].SegmentNum();j++){
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
aEPS.OutFile << "n " << aEPS.FunX(Mid(Q[i].S[j].P[0](XDisp))) << " " << aEPS.FunY(Mid(Q[i].S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q[i].S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(Mid(Q[i].S[j].P[k](XDisp))) << " " << aEPS.FunY(Mid(Q[i].S[j].P[k](YDisp))) << " l ";
}
aEPS.OutFile << "s" << endl;
if (DoSymmetric){
aEPS.OutFile << "n " << aEPS.FunX(-Mid(Q[i].S[j].P[0](XDisp))) << " " << aEPS.FunY(-Mid(Q[i].S[j].P[0](YDisp))) << " m ";
for (int k=1;k<Q[i].S[j].PointNum();k++){
aEPS.OutFile << aEPS.FunX(-Mid(Q[i].S[j].P[k](XDisp))) << " " << aEPS.FunY(-Mid(Q[i].S[j].P[k](YDisp))) << " l ";
}
aEPS.OutFile << "s" << endl;
}
}
}
}
aEPS.OutFile << "grestore\n";
if (PlotLegend){
aEPS.OutFile << "gsave\n";
aEPS.OutFile.precision(5);
aEPS.ChangeCoordinates(0,1.0*(EPSOpt.bbXMax-EPSOpt.bbXMin)/(EPSOpt.bbYMax-EPSOpt.bbYMin),0,1.0);
aEPS.SetFont(0.027,"Times-Roman");
for (int i=0;i<QuadrNum();i++){
bool top=true;
double xofset=0.02;
double xsize=0.02;
double xskip=0.01;
double yofset=xofset;
if (top) yofset=1-xofset-xsize;
double ysize=xsize;
double yskip=xskip;
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
double ypos=yofset+(ysize+yskip)*i;
if (top) ypos=yofset-(ysize+yskip)*i;
aEPS.FillRect(xofset,ypos,xofset+xsize,ypos+ysize);
aEPS.OutFile << "0 0 0 rgb\n";
// char WS[3];
// sprintf(WS,"%2d",i+1);
string s=to_string(i+1);
aEPS.DrawStringRight(xofset+2.5*xsize,ypos,s.c_str());
}
aEPS.OutFile << "grestore\n";
}
aEPS.Finish();
cout << "gv " << FileName << endl;
}
void PolygonSetSaveEPS(const PolygonSet & Polygons, string FileName, EPSOptions & EPSOpt, int XDisp, int YDisp, bool DoSymmetric, bool PlotLegend, bool UseSingleColor, bool fill){
cout << "PolygonSetSaveEPS" << endl;
EPSOpt.EPSFileCreator="cc (Z. Galias)";
ostringstream oss;
oss << "Polygons";
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
bool PlotPolygons=true;
if (PlotPolygons){
for (int i=0;i<Polygons.size();i++){
IVector Range=GetRange(Polygons[i]);
if (Inf(Range(XDisp))>EPSOpt.XMax || Sup(Range(XDisp))<EPSOpt.XMin || Inf(Range(YDisp))>EPSOpt.YMax || Sup(Range(YDisp))<EPSOpt.YMin) continue;
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
aEPS.OutFile << "n " << aEPS.FunX(Mid(Polygons[i][0](XDisp))) << " " << aEPS.FunY(Mid(Polygons[i][0](YDisp))) << " m ";
for (int k=1;k<Polygons[i].size();k++){
aEPS.OutFile << aEPS.FunX(Mid(Polygons[i][k](XDisp))) << " " << aEPS.FunY(Mid(Polygons[i][k](YDisp))) << " l ";
}
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
if (DoSymmetric){
aEPS.OutFile << "n " << aEPS.FunX(-Mid(Polygons[i][0](XDisp))) << " " << aEPS.FunY(-Mid(Polygons[i][0](YDisp))) << " m ";
for (int k=1;k<Polygons[i].size();k++){
aEPS.OutFile << aEPS.FunX(-Mid(Polygons[i][k](XDisp))) << " " << aEPS.FunY(-Mid(Polygons[i][k](YDisp))) << " l ";
}
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
}
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
aEPS.OutFile << "n " << aEPS.FunX(Mid(Polygons[i][0](XDisp))) << " " << aEPS.FunY(Mid(Polygons[i][0](YDisp))) << " m ";
for (int k=1;k<Polygons[i].size();k++){
aEPS.OutFile << aEPS.FunX(Mid(Polygons[i][k](XDisp))) << " " << aEPS.FunY(Mid(Polygons[i][k](YDisp))) << " l ";
}
aEPS.OutFile << "s" << endl;
if (DoSymmetric){
aEPS.OutFile << "n " << aEPS.FunX(-Mid(Polygons[i][0](XDisp))) << " " << aEPS.FunY(-Mid(Polygons[i][0](YDisp))) << " m ";
for (int k=1;k<Polygons[i].size();k++){
aEPS.OutFile << aEPS.FunX(-Mid(Polygons[i][k](XDisp))) << " " << aEPS.FunY(-Mid(Polygons[i][k](YDisp))) << " l ";
}
aEPS.OutFile << "s" << endl;
}
}
}
aEPS.OutFile << "grestore\n";
if (PlotLegend){
aEPS.OutFile << "gsave\n";
aEPS.OutFile.precision(5);
aEPS.ChangeCoordinates(0,1.0*(EPSOpt.bbXMax-EPSOpt.bbXMin)/(EPSOpt.bbYMax-EPSOpt.bbYMin),0,1.0);
aEPS.SetFont(0.027,"Times-Roman");
for (int i=0;i<Polygons.size();i++){
bool top=true;
double xofset=0.02;
double xsize=0.02;
double xskip=0.01;
double yofset=xofset;
if (top) yofset=1-xofset-xsize;
double ysize=xsize;
double yskip=xskip;
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
double ypos=yofset+(ysize+yskip)*i;
if (top) ypos=yofset-(ysize+yskip)*i;
aEPS.FillRect(xofset,ypos,xofset+xsize,ypos+ysize);
aEPS.OutFile << "0 0 0 rgb\n";
// char WS[3];
// sprintf(WS,"%2d",i+1);
string s=to_string(i+1);
aEPS.DrawStringRight(xofset+2.5*xsize,ypos,s.c_str());
}
aEPS.OutFile << "grestore\n";
}
aEPS.Finish();
cout << "gv " << FileName << endl;
}
void VectorSetSave(const DVectorSet & VSet, string FileName, int prec){
if (VSet.size()==0){ cerr << "VectorSetSave: emptyset" << endl; return;}
ofstream OutFile(FileName.c_str()) ;
cout << "Saving " << FileName << endl;
OutFile << setprecision(prec);
OutFile << "Num " << VSet.size() << endl;
OutFile << "Dim " << VSet[0].dimension() << endl;
for (int i=0;i<VSet.size();i++){
OutFile << "V";
for (int d=0;d<VSet[i].dimension();++d) OutFile << ' ' << VSet[i][d];
OutFile << endl;
}
OutFile.close();
}
void VectorSetLoad(string FileName, DVectorSet & VSet){
VSet.clear();
ifstream InFile(FileName.c_str());
if (!InFile){ cerr << "File " << FileName << " not found." << endl; return; }
cout << "Reading " << FileName << "...\n";
char ch;
int Num=0;
int dim=0;
while (InFile.get(ch)){
string WS;
if (ch=='#'){
getline(InFile, WS);
continue;
}
InFile.putback(ch);
InFile >> WS;
if (WS=="Num") InFile >> Num;
if (WS=="Dim") InFile >> dim;
if (WS=="V"){
DVector v(dim);
for (int d=0;d<dim;++d){
double x;
InFile >> x;
v[d]=x;
}
VSet.push_back(v);
}
}
if (VSet.size()!=Num) cerr << "VectorSetLoad read error.";
InFile.close();
}
void IVectorSetSave(const IVectorSet & IVSet, string FileName, int prec){
if (IVSet.size()==0){ cerr << "IVectorSetSave: emptyset" << endl; return;}
ofstream OutFile(FileName.c_str()) ;
cout << "Saving " << FileName << endl;
OutFile << setprecision(prec);
OutFile << "Num " << IVSet.size() << endl;
OutFile << "Dim " << IVSet[0].dimension() << endl;
for (int i=0;i<IVSet.size();i++){
OutFile << "V";
for (int d=0;d<IVSet[i].dimension();++d) OutFile << ' ' << Inf(IVSet[i][d]) << ' ' << Sup(IVSet[i][d]);
OutFile << endl;
}
OutFile.close();
}
void IVectorSetLoad(string FileName, IVectorSet & IVSet){
IVSet.clear();
ifstream InFile(FileName.c_str());
if (!InFile){ cerr << "File " << FileName << " not found." << endl; return; }
cout << "Reading " << FileName << "...\n";
char ch;
int Num=0;
int dim=0;
while (InFile.get(ch)){
string WS;
if (ch=='#'){
getline(InFile, WS);
continue;
}
InFile.putback(ch);
InFile >> WS;
if (WS=="Num") InFile >> Num;
if (WS=="Dim") InFile >> dim;
if (WS=="V"){
IVector v(dim);
for (int d=0;d<dim;++d){
double x1,x2;
InFile >> x1 >> x2;
v[d]=Interval(x1,x2);
}
IVSet.push_back(v);
}
}
if (IVSet.size()!=Num) cerr << "IVectorSetLoad read error.";
InFile.close();
}
void IVectorSetsSaveEPS(const IVectorSets & IVSets, string FileName, EPSOptions & EPSOpt, int XDisp, int YDisp, bool PlotLegend, bool UseSingleColor, bool fill){
cout << "IVectorSetsSaveEPS" << endl;
EPSOpt.EPSFileCreator="cc (Z. Galias)";
ostringstream oss;
oss << "IVectorSets";
EPSOpt.EPSFileTitle=oss.str();
cout << "Saving " << FileName << endl;
SaveEPS aEPS(FileName.c_str());
aEPS.Init(EPSOpt);
aEPS.OutFile << "grestore\n";
aEPS.OutFile.precision(10);
aEPS.OutFile << "gsave\n";
aEPS.RescaleToBoundingBox(EPSOpt.XMin,EPSOpt.XMax,EPSOpt.YMin,EPSOpt.YMax);
aEPS.SetClip(EPSOpt.bbXMin,EPSOpt.bbXMax,EPSOpt.bbYMin,EPSOpt.bbYMax);
aEPS.DefSizesAutomaticOut(EPSOpt.Diameter);
aEPS.DefLineWidthAutomaticOut(EPSOpt.LineWidth);
for (int i=0;i<IVSets.size();i++){
aEPS.OutFile << ColorString(UseSingleColor?0:i) << " rgb\n";
for (int j=0;j<IVSets[i].size();++j){
aEPS.OutFile << "n ";
aEPS.OutFile << aEPS.FunX(Inf(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Inf(IVSets[i][j](YDisp))) << " m ";
aEPS.OutFile << aEPS.FunX(Inf(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Sup(IVSets[i][j](YDisp))) << " l ";
aEPS.OutFile << aEPS.FunX(Sup(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Sup(IVSets[i][j](YDisp))) << " l ";
aEPS.OutFile << aEPS.FunX(Sup(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Inf(IVSets[i][j](YDisp))) << " l ";
if (fill) aEPS.OutFile << "c f\n";
else aEPS.OutFile << "c s\n";
}
}
if (fill){
for (int i=0;i<IVSets.size();i++){
aEPS.OutFile << "0 0 0" << " rgb\n";
for (int j=0;j<IVSets[i].size();++j){
aEPS.OutFile << "n ";
aEPS.OutFile << aEPS.FunX(Inf(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Inf(IVSets[i][j](YDisp))) << " m ";
aEPS.OutFile << aEPS.FunX(Inf(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Sup(IVSets[i][j](YDisp))) << " l ";
aEPS.OutFile << aEPS.FunX(Sup(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Sup(IVSets[i][j](YDisp))) << " l ";
aEPS.OutFile << aEPS.FunX(Sup(IVSets[i][j](XDisp))) << " " << aEPS.FunY(Inf(IVSets[i][j](YDisp))) << " l ";
aEPS.OutFile << "c s\n";
}
}
}
aEPS.OutFile << "grestore\n";
aEPS.Finish();
cout << "NumberOfBoxes: ";
int BoxNum=0;
for (int i=0;i<IVSets.size();i++){
cout << IVSets[i].size() << ' ';
BoxNum+=IVSets[i].size();
}
cout << "Total: " << BoxNum << endl;
cout << "gv " << FileName << endl;
}
// IntPointInPoly(): winding number test for a point in a polygon
// Input: P = a point,
// V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
// Return: wn = the winding number (=0 only when P is outside)
int IntPointInPoly(IntPoint P, const vector <IntPoint> & V){
int n=V.size()-1;
int wn = 0; // the winding number counter
// loop through all edges of the polygon
for (int i=0;i<n;i++){ // edge from V[i] to V[i+1]
if (V[i].y<=P.y){ // start y <= P.y
if (V[i+1].y>P.y) // an upward crossing
if (isLeft( V[i], V[i+1], P) > 0) ++wn; // P left of edge, have a valid up intersect
}
else { // start y > P.y (no test needed)
if (V[i+1].y<=P.y) // a downward crossing
if (isLeft(V[i],V[i+1],P)<0) --wn; // P right of edge, have a valid down intersect
}
}
return wn;
}
bool PointNotInSegment(const IVector & P0, const IVector & P1, const IVector & P, int x, int y){
if (P[x]<P0[x] && P[x]<P1[x]) return true;
if (P[x]>P0[x] && P[x]>P1[x]) return true;
if (P[y]<P0[y] && P[y]<P1[y]) return true;
if (P[y]>P0[y] && P[y]>P1[y]) return true;
Interval tmp=(P1[x]-P0[x])*(P[y]-P0[y])-(P[x]-P0[x])*(P1[y]-P0[y]);
if (!tmp.contains(0.0)) return true; // P does not belong to the line containing (P0,P1)
return false;
}
bool PointNotInBrokenLine(const IVector & P, const BrokenLine & BL, int x, int y){
for (int s=0;s<BL.size()-1;s++){ // process edge from BL[s] to BL[s+1]
if (!PointNotInSegment(BL[s],BL[s+1],P,x,y)) return false; // P has non-empty intersection with the edge from SS[s] to SS[s+1]
}
return true;
}
// PointInPoly(): winding number test for a point in a polygon
// Input: P = a point,
// Poly[] = vertex points of a polygon Poly[n+1] with Poly[n]=Poly[0]
// Return: (0 when P is outside, 1 when P is inside, -1, when do not know)
int PointInPolyInterior(const IVector & P, const Polygon & Poly, int x, int y){
int n=Poly.size()-1;
for (int i=0;i<n;i++){ // process edge from Poly[i] to Poly[i+1]
if (!PointNotInSegment(Poly[i],Poly[i+1],P,x,y)) return -1; // P has non-empty intersection with the edge from Poly[i] to Poly[i+1]
}
// now we know that P has non-empty intersection with the border of polygon Poly
// it is sufficient to verify conditions for a single point in P
IVector PP(P.dimension());
for (int d=0;d<P.dimension();++d) PP[d]=P[d].leftBound();
int wn = 0; // the winding number counter
// loop through all edges of the polygon
for (int i=0;i<n;i++){ // edge from Poly[i] to Poly[i+1]
if (Poly[i][y]<=PP[y]){ // start y <= P[y]
if (Poly[i+1][y]>PP[y]){ // an upward crossing
int res=isLeftOrRight(Poly[i],Poly[i+1],PP,x,y);
if (res==0) return -1; // do not know
if (res==+1) ++wn; // P left of edge, have a valid up intersect
}
}
else { // start y > P[y] (no test needed)
if (Poly[i+1][y]<=PP[y]){ // a downward crossing
int res=isLeftOrRight(Poly[i],Poly[i+1],PP,x,y);
if (res==0) return -1; // do not know
if (res==-1) --wn; // P right of edge, have a valid down intersect
}
}
}
if (wn==0) return 0; // point P is outside polygon Poly
return 1;
}
void TestIntPointInPoly(){
vector <IntPoint> Poly;
Poly.push_back(IntPoint(0,0));
Poly.push_back(IntPoint(0,2));
Poly.push_back(IntPoint(2,2));
Poly.push_back(IntPoint(2,0));
Poly.push_back(IntPoint(0,0));
cout << "IntPointInPoly: " << IntPointInPoly(IntPoint(1,1),Poly) << endl;
cout << "IntPointInPoly: " << IntPointInPoly(IntPoint(3,1),Poly) << endl;
}
void TestPointInPoly(){
vector <IVector> Poly;
IVector P(2);
P[0]=0; P[1]=0;
Poly.push_back(P);
P[0]=0; P[1]=2;
Poly.push_back(P);
P[0]=2; P[1]=2;
Poly.push_back(P);
P[0]=2; P[1]=0;
Poly.push_back(P);
P[0]=0; P[1]=0;
Poly.push_back(P);
P[0]=1;
P[1]=1;
cout << "Result should be 1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=3;
P[1]=1;
cout << "Result should be 0: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=2;
P[1]=3;
cout << "Result should be 0: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=2;
P[1]=1;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=2;
P[1]=0;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=1;
P[1]=0;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=3;
P[1]=0;
cout << "Result should be 0: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=2;
P[1]=2;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=0;
P[1]=0;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=2;
P[1]=0;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=0;
P[1]=2;
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=Interval(0.9,1.1);
P[1]=Interval(0.9,1.1);
cout << "Result should be 1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=Interval(2.9,3.1);
P[1]=Interval(0.9,1.1);
cout << "Result should be 0: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=Interval(1.9,2.1);
P[1]=Interval(2.9,3.1);
cout << "Result should be 0: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
P[0]=Interval(1.9,2.1);
P[1]=Interval(-0.1,0.1);
cout << "Result should be -1: PointInPoly: " << PointInPolyInterior(P,Poly,0,1) << endl;
IVector P0(2);
IVector P1(2);
P0[0]=0;
P0[1]=0;
P1[0]=2;
P1[1]=2;
P[0]=Interval(1.1,2.0);
P[1]=Interval(0.0,0.9);
cout << "Result should be true: PointNotInSegment: " << (PointNotInSegment(P0,P1,P,0,1)==true?"true":"false") << endl;
P[0]=Interval(2.1,2.3);
P[1]=Interval(1.7,1.9);
cout << "Result should be true: PointNotInSegment: " << (PointNotInSegment(P0,P1,P,0,1)==true?"true":"false") << endl;
P[0]=Interval(1.0);
P[1]=Interval(1.0);
cout << "Result should be false: PointNotInSegment: " << (PointNotInSegment(P0,P1,P,0,1)==true?"true":"false") << endl;
}