Commit f0af5788 authored by kraus's avatar kraus
Browse files

adding polygons to multi-slit language; first step towards #238

parent fae0d63b
...@@ -25,6 +25,7 @@ set (_SRCS ...@@ -25,6 +25,7 @@ set (_SRCS
SingularMatrixError.cpp SingularMatrixError.cpp
SizeError.cpp SizeError.cpp
SwitcherError.cpp SwitcherError.cpp
Mesher.cpp
Util.cpp Util.cpp
) )
...@@ -61,6 +62,7 @@ set (HDRS ...@@ -61,6 +62,7 @@ set (HDRS
SingularMatrixError.h SingularMatrixError.h
SizeError.h SizeError.h
SwitcherError.h SwitcherError.h
Mesher.h
Util.h Util.h
) )
......
#include "Utilities/MSLang.h" #include "Utilities/MSLang.h"
#include "Utilities/PortableBitmapReader.h" #include "Utilities/PortableBitmapReader.h"
#include "Utilities/Mesher.h"
#include "Algorithms/Quaternion.h" #include "Algorithms/Quaternion.h"
#include "Physics/Physics.h" #include "Physics/Physics.h"
...@@ -273,6 +274,90 @@ namespace mslang { ...@@ -273,6 +274,90 @@ namespace mslang {
} }
}; };
void Triangle::print(int indentwidth) {
std::string indent(indentwidth, ' ');
std::string indentbase(4, ' ');
Vector_t origin = trafo_m.getOrigin();
double angle = trafo_m.getAngle() * Physics::rad2deg;
std::cout << indent << "triangle, \n";
for (unsigned int i = 0; i < 3u; ++ i) {
std::cout << indent << indentbase << "node " << i << ": " << nodes_m[i] << "\n";
}
std::cout << indent << indentbase << "origin: " << origin[0] << ", " << origin[1] << ",\n"
<< indent << indentbase << "angle: " << angle << "\n"
<< indent << indentbase << trafo_m(0, 0) << "\t" << trafo_m(0, 1) << "\t" << trafo_m(0, 2) << "\n"
<< indent << indentbase << trafo_m(1, 0) << "\t" << trafo_m(1, 1) << "\t" << trafo_m(1, 2) << "\n"
<< indent << indentbase << trafo_m(2, 0) << "\t" << trafo_m(2, 1) << "\t" << trafo_m(2, 2) << std::endl;
}
void Triangle::apply(std::vector<Base*> &bfuncs) {
bfuncs.push_back(this->clone());
}
Base* Triangle::clone() const {
Triangle *tri = new Triangle(*this);
return tri;
}
void Triangle::writeGnuplot(std::ofstream &out) const {
unsigned int width = out.precision() + 8;
for (unsigned int i = 0; i < 4u; ++ i) {
Vector_t corner = trafo_m.transformFrom(nodes_m[i % 3u]);
out << std::setw(width) << corner[0]
<< std::setw(width) << corner[1]
<< std::endl;
}
out << std::endl;
bb_m.writeGnuplot(out);
}
void Triangle::computeBoundingBox() {
std::vector<Vector_t> corners;
for (unsigned int i = 0; i < 3u; ++ i) {
corners.push_back(trafo_m.transformFrom(nodes_m[i]));
}
Vector_t llc = corners[0], urc = corners[0];
for (unsigned int i = 1u; i < 3u; ++ i) {
if (corners[i][0] < llc[0]) llc[0] = corners[i][0];
else if (corners[i][0] > urc[0]) urc[0] = corners[i][0];
if (corners[i][1] < llc[1]) llc[1] = corners[i][1];
else if (corners[i][1] > urc[1]) urc[1] = corners[i][1];
}
bb_m = BoundingBox(llc, urc);
}
double Triangle::crossProduct(const Vector_t &pt, unsigned int nodeNum) const {
nodeNum = nodeNum % 3u;
unsigned int nextNode = (nodeNum + 1) % 3u;
Vector_t nodeToPt = pt - nodes_m[nodeNum];
Vector_t nodeToNext = nodes_m[nextNode] - nodes_m[nodeNum];
return nodeToPt[0] * nodeToNext[1] - nodeToPt[1] * nodeToNext[0];
}
bool Triangle::isInside(const Vector_t &R) const {
bool test0 = (crossProduct(R, 0) <= 0.0);
bool test1 = (crossProduct(R, 1) <= 0.0);
bool test2 = (crossProduct(R, 2) <= 0.0);
return test0 && test1 && test2;
}
void Triangle::orientNodesCCW() {
if (crossProduct(nodes_m[0], 1) > 0.0) {
std::swap(nodes_m[1], nodes_m[2]);
}
}
struct Repeat: public Function { struct Repeat: public Function {
Function* func_m; Function* func_m;
unsigned int N_m; unsigned int N_m;
...@@ -648,6 +733,67 @@ namespace mslang { ...@@ -648,6 +733,67 @@ namespace mslang {
std::vector<Rectangle> pixels_m; std::vector<Rectangle> pixels_m;
}; };
struct Polygon: public Function {
std::vector<Triangle> triangles_m;
void triangulize(std::vector<Vector_t> &nodes) {
Mesher mesher(nodes);
triangles_m = mesher.getTriangles();
}
static
bool parse_detail(iterator &it, const iterator &end, Function* &fun) {
Polygon *poly = static_cast<Polygon*>(fun);
std::vector<Vector_t> nodes;
std::string str(it, end);
boost::regex argument(Double + "," + Double + ";(.*)");
boost::regex argumentEnd(Double + "," + Double + "(\\).*)");
boost::smatch what;
while (boost::regex_match(str, what, argument)) {
Vector_t p(atof(std::string(what[1]).c_str()),
atof(std::string(what[3]).c_str()),
1.0);
nodes.push_back(p);
std::string fullMatch = what[0];
std::string rest = what[5];
it += (fullMatch.size() - rest.size());
str = std::string(it, end);
}
if (!boost::regex_match(str, what, argumentEnd) ||
nodes.size() < 2) return false;
Vector_t p(atof(std::string(what[1]).c_str()),
atof(std::string(what[3]).c_str()),
1.0);
nodes.push_back(p);
std::string fullMatch = what[0];
std::string rest = what[5];
it += (fullMatch.size() - rest.size() + 1);
str = std::string(it, end);
poly->triangulize(nodes);
return true;
}
virtual void print(int ident) {
// for (auto pix: pixels_m) pix.print(ident);
}
virtual void apply(std::vector<Base*> &bfuncs) {
for (Triangle &tri: triangles_m)
bfuncs.push_back(tri.clone());
}
};
QuadTree::QuadTree(const QuadTree &right): QuadTree::QuadTree(const QuadTree &right):
level_m(right.level_m), level_m(right.level_m),
objects_m(right.objects_m.begin(), objects_m(right.objects_m.begin(),
...@@ -801,58 +947,51 @@ namespace mslang { ...@@ -801,58 +947,51 @@ namespace mslang {
if (identifier == "rectangle") { if (identifier == "rectangle") {
fun = new Rectangle; fun = new Rectangle;
/*iterator it2 = */it += shift; it += shift;
if (!Rectangle::parse_detail(it, end, fun)) return false; if (!Rectangle::parse_detail(it, end, fun)) return false;
// it = it2;
return true; return true;
} else if (identifier == "ellipse") { } else if (identifier == "ellipse") {
fun = new Ellipse; fun = new Ellipse;
/*iterator it2 = */it += shift; it += shift;
if (!Ellipse::parse_detail(it, end, fun)) return false; if (!Ellipse::parse_detail(it, end, fun)) return false;
// it = it2; return true;
} else if (identifier == "polygon") {
fun = new Polygon;
it += shift;
if (!Polygon::parse_detail(it, end, fun)) return false;
return true; return true;
} else if (identifier == "repeat") { } else if (identifier == "repeat") {
fun = new Repeat; fun = new Repeat;
/*iterator it2 = */it += shift; it += shift;
if (!Repeat::parse_detail(it, end, fun)) return false; if (!Repeat::parse_detail(it, end, fun)) return false;
// it = it2;
return true; return true;
} else if (identifier == "rotate") { } else if (identifier == "rotate") {
fun = new Rotate; fun = new Rotate;
// iterator it2 =
it += shift; it += shift;
if (!Rotate::parse_detail(it, end, fun)) return false; if (!Rotate::parse_detail(it, end, fun)) return false;
// // it = it2;
return true; return true;
} else if (identifier == "translate") { } else if (identifier == "translate") {
fun = new Translate; fun = new Translate;
/*iterator it2 = */it += shift; it += shift;
if (!Translate::parse_detail(it, end, fun)) return false; if (!Translate::parse_detail(it, end, fun)) return false;
// it = it2;
return true; return true;
} else if (identifier == "shear") { } else if (identifier == "shear") {
fun = new Shear; fun = new Shear;
/*iterator it2 = */it += shift; it += shift;
if (!Shear::parse_detail(it, end, fun)) return false; if (!Shear::parse_detail(it, end, fun)) return false;
// it = it2;
return true; return true;
} else if (identifier == "union") { } else if (identifier == "union") {
fun = new Union; fun = new Union;
/*iterator it2 = */it += shift; it += shift;
if (!Union::parse_detail(it, end, fun)) return false; if (!Union::parse_detail(it, end, fun)) return false;
// it = it2;
return true; return true;
} else if (identifier == "mask") { } else if (identifier == "mask") {
fun = new Mask; fun = new Mask;
......
#ifndef MSLANG_H
#define MSLANG_H
#include "Algorithms/Vektor.h" #include "Algorithms/Vektor.h"
#include "AppTypes/Tenzor.h" #include "AppTypes/Tenzor.h"
...@@ -33,6 +36,13 @@ namespace mslang { ...@@ -33,6 +36,13 @@ namespace mslang {
height_m(urc[1] - llc[1]) height_m(urc[1] - llc[1])
{ } { }
bool doesIntersect(const BoundingBox &bb) {
return ((center_m[0] - 0.5 * width_m <= bb.center_m[0] + 0.5 * bb.width_m) &&
(center_m[0] + 0.5 * width_m >= bb.center_m[0] - 0.5 * bb.width_m) &&
(center_m[1] - 0.5 * height_m <= bb.center_m[1] + 0.5 * bb.height_m) &&
(center_m[1] + 0.5 * height_m >= bb.center_m[1] - 0.5 * bb.height_m));
}
bool isInside(const Vector_t &X) const { bool isInside(const Vector_t &X) const {
if (2 * std::abs(X[0] - center_m[0]) <= width_m && if (2 * std::abs(X[0] - center_m[0]) <= width_m &&
2 * std::abs(X[1] - center_m[1]) <= height_m) 2 * std::abs(X[1] - center_m[1]) <= height_m)
...@@ -162,6 +172,31 @@ namespace mslang { ...@@ -162,6 +172,31 @@ namespace mslang {
}; };
struct Triangle: public Base {
std::vector<Vector_t> nodes_m;
Triangle():
Base(),
nodes_m(std::vector<Vector_t>(3, Vector_t(0, 0, 1)))
{ }
Triangle(const Triangle &right):
Base(right),
nodes_m(right.nodes_m)
{ }
virtual ~Triangle()
{ }
virtual void print(int indentwidth);
virtual void apply(std::vector<Base*> &bfuncs);
virtual Base* clone() const;
virtual void writeGnuplot(std::ofstream &out) const;
virtual void computeBoundingBox();
double crossProduct(const Vector_t &pt, unsigned int nodeNum) const;
virtual bool isInside(const Vector_t &R) const;
void orientNodesCCW();
};
struct QuadTree { struct QuadTree {
int level_m; int level_m;
std::list<Base*> objects_m; std::list<Base*> objects_m;
...@@ -210,4 +245,6 @@ namespace mslang { ...@@ -210,4 +245,6 @@ namespace mslang {
}; };
bool parse(std::string str, Function* &fun); bool parse(std::string str, Function* &fun);
} }
\ No newline at end of file
#endif
\ No newline at end of file
#include "Utilities/Mesher.h"
Mesher::Mesher(std::vector<Vector_t> &vertices):
vertices_m(vertices)
{
apply();
}
std::vector<mslang::Triangle> Mesher::getTriangles() const {
return triangles_m;
}
void Mesher::orientVerticesCCW() {
const unsigned int size = vertices_m.size();
double sum = 0.0;
for (unsigned int i = 0; i < size; ++ i) {
unsigned int iPlusOne = (i + 1) % size;
Vector_t edge(vertices_m[iPlusOne][0] - vertices_m[i][0],
vertices_m[iPlusOne][1] + vertices_m[i][1],
0.0);
sum += edge[0] * edge[1];
}
if (sum <= 0.0) return;
std::vector<Vector_t> reoriented(vertices_m.rbegin(), vertices_m.rend());
vertices_m.swap(reoriented);
}
bool Mesher::isConvex(unsigned int i) const {
unsigned int numVertices = vertices_m.size();
i = i % numVertices;
unsigned int iPlusOne = (i + 1) % numVertices;
unsigned int iMinusOne = (i + numVertices - 1) % numVertices;
Vector_t edge0 = vertices_m[iMinusOne] - vertices_m[i];
Vector_t edge1 = vertices_m[iPlusOne] - vertices_m[i];
double vectorProduct = edge0[0] * edge1[1] - edge0[1] * edge1[0];
return (vectorProduct < 0.0);
}
double Mesher::crossProduct(const Vector_t &a,
const Vector_t &b) const {
return a[0] * b[1] - a[1] * b[0];
}
bool Mesher::isPointOnLine(unsigned int i,
unsigned int j,
const Vector_t &pt) const {
Vector_t aTmp = vertices_m[j] - vertices_m[i];
Vector_t bTmp = pt - vertices_m[i];
double r = crossProduct(aTmp, bTmp);
return std::abs(r) < 1e-10;
}
bool Mesher::isPointRightOfLine(unsigned int i,
unsigned int j,
const Vector_t &pt) const {
Vector_t aTmp = vertices_m[j] - vertices_m[i];
Vector_t bTmp = pt - vertices_m[i];
return crossProduct(aTmp, bTmp) < 0.0;
}
bool Mesher::lineSegmentTouchesOrCrossesLine(unsigned int i,
unsigned int j,
unsigned int k,
unsigned int l) const {
return (isPointOnLine(i, j, vertices_m[k]) ||
isPointOnLine(i, j, vertices_m[l]) ||
(isPointRightOfLine(i, j, vertices_m[k]) ^
isPointRightOfLine(i, j, vertices_m[l])));
}
bool Mesher::isPotentialEdgeIntersected(unsigned int i) const {
unsigned int numVertices = vertices_m.size();
if (numVertices < 5) return false;
i = i % numVertices;
unsigned int iPlusOne = (i + 1) % numVertices;
unsigned int iMinusOne = (i + numVertices - 1) % numVertices;
mslang::BoundingBox bbPotentialNewEdge;
bbPotentialNewEdge.center_m = 0.5 * (vertices_m[iMinusOne] + vertices_m[iPlusOne]);
bbPotentialNewEdge.width_m = std::abs(vertices_m[iMinusOne][0] - vertices_m[iPlusOne][0]);
bbPotentialNewEdge.height_m = std::abs(vertices_m[iMinusOne][1] - vertices_m[iPlusOne][1]);
for (unsigned int j = iPlusOne + 1; j < iPlusOne + numVertices - 3; ++ j) {
unsigned int k = (j % numVertices);
unsigned int kPlusOne = ((k + 1) % numVertices);
mslang::BoundingBox bbThisEdge;
bbThisEdge.center_m = 0.5 * (vertices_m[k] + vertices_m[kPlusOne]);
bbThisEdge.width_m = std::abs(vertices_m[k][0] - vertices_m[kPlusOne][0]);
bbThisEdge.height_m = std::abs(vertices_m[k][1] - vertices_m[kPlusOne][1]);
if (bbPotentialNewEdge.doesIntersect(bbThisEdge) &&
lineSegmentTouchesOrCrossesLine(iPlusOne, iMinusOne,
k, kPlusOne) &&
lineSegmentTouchesOrCrossesLine(k, kPlusOne,
iPlusOne, iMinusOne))
return true;
}
return false;
}
bool Mesher::isPointInsideCone(unsigned int i,
unsigned int j,
unsigned int jPlusOne,
unsigned int jMinusOne) const {
const Vector_t &pt = vertices_m[i];
return !(isPointOnLine(jMinusOne, j, pt) ||
isPointRightOfLine(jMinusOne, j, pt) ||
isPointOnLine(j, jPlusOne, pt) ||
isPointRightOfLine(j, jPlusOne, pt));
}
bool Mesher::isEar(unsigned int i) const {
unsigned int size = vertices_m.size();
unsigned int iMinusTwo = (i + size - 2) % size;
unsigned int iMinusOne = (i + size - 1) % size;
unsigned int iPlusOne = (i + 1) % size;
unsigned int iPlusTwo = (i + 2) % size;
bool convex = isConvex(i);
bool isInsideCone1 = isPointInsideCone(iMinusOne, iPlusOne, iPlusTwo, i);
bool isInsideCone2 = isPointInsideCone(iPlusOne, iMinusOne, i, iMinusTwo);
bool notCrossed = !isPotentialEdgeIntersected(i);
return (convex &&
isInsideCone1 &&
isInsideCone2 &&
notCrossed);
}
std::vector<unsigned int> Mesher::findAllEars() const {
unsigned int size = vertices_m.size();
std::vector<unsigned int> ears;
for (unsigned int i = 0; i < size; ++ i) {
if (isEar(i)) {
ears.push_back(i);
}
}
return ears;
}
double Mesher::computeMinimumAngle(unsigned int i) const {
unsigned int numVertices = vertices_m.size();
unsigned int previous = (i + numVertices - 1) % numVertices;
unsigned int next = (i + 1) % numVertices;
Vector_t edge0 = vertices_m[i] - vertices_m[previous];
Vector_t edge1 = vertices_m[next] - vertices_m[i];
Vector_t edge2 = vertices_m[previous] - vertices_m[next];
double length0 = euclidean_norm(edge0);
double length1 = euclidean_norm(edge1);
double length2 = euclidean_norm(edge2);
double angle01 = std::asin((edge0[0] * edge1[1] - edge0[1] * edge1[0]) / length0 / length1);
double angle12 = std::asin((edge1[0] * edge2[1] - edge1[1] * edge2[0]) / length1 / length2);
double angle20 = M_PI - angle01 - angle12;
return std::min(std::min(angle01, angle12), angle20);
}
unsigned int Mesher::selectBestEar(std::vector<unsigned int> &ears) const {
unsigned int numEars = ears.size();
double maxMinAngle = computeMinimumAngle(ears[0]);
unsigned int earWithMaxMinAngle = 0;
for (unsigned int i = 1; i < numEars; ++ i) {
double angle = computeMinimumAngle(ears[i]);
if (angle > maxMinAngle) {
maxMinAngle = angle;
earWithMaxMinAngle = i;
}
}
return ears[earWithMaxMinAngle];
}
void Mesher::apply() {
orientVerticesCCW();
unsigned int numVertices = vertices_m.size();
while (numVertices > 3) {
std::vector<unsigned int> allEars = findAllEars();
unsigned int bestEar = selectBestEar(allEars);
unsigned int next = (bestEar + 1) % numVertices;
unsigned int previous = (bestEar + numVertices - 1) % numVertices;
mslang::Triangle tri;
tri.nodes_m[0] = vertices_m[previous];
tri.nodes_m[1] = vertices_m[bestEar];
tri.nodes_m[2] = vertices_m[next];
triangles_m.push_back(tri);
vertices_m.erase(vertices_m.begin() + bestEar);
-- numVertices;
}
mslang::Triangle tri;
tri.nodes_m = vertices_m;
triangles_m.push_back(tri);
}
\ No newline at end of file