diff --git a/src/Classic/Utilities/CMakeLists.txt b/src/Classic/Utilities/CMakeLists.txt index b88b5026cbd6578e6b3d6bc96c7837b7e43c0e2a..8cb66c78e82c9103093ef6a7f05e22c29fc6466a 100644 --- a/src/Classic/Utilities/CMakeLists.txt +++ b/src/Classic/Utilities/CMakeLists.txt @@ -25,6 +25,7 @@ set (_SRCS SingularMatrixError.cpp SizeError.cpp SwitcherError.cpp + Mesher.cpp Util.cpp ) @@ -61,6 +62,7 @@ set (HDRS SingularMatrixError.h SizeError.h SwitcherError.h + Mesher.h Util.h ) diff --git a/src/Classic/Utilities/MSLang.cpp b/src/Classic/Utilities/MSLang.cpp index 9ffe5bb384ca0d55eb5864e206e554805f5e039a..e60cfde1d9f939df387b6b15d4f040655f586530 100644 --- a/src/Classic/Utilities/MSLang.cpp +++ b/src/Classic/Utilities/MSLang.cpp @@ -1,5 +1,6 @@ #include "Utilities/MSLang.h" #include "Utilities/PortableBitmapReader.h" +#include "Utilities/Mesher.h" #include "Algorithms/Quaternion.h" #include "Physics/Physics.h" @@ -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 { Function* func_m; unsigned int N_m; @@ -648,6 +733,67 @@ namespace mslang { 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): level_m(right.level_m), objects_m(right.objects_m.begin(), @@ -801,58 +947,51 @@ namespace mslang { if (identifier == "rectangle") { fun = new Rectangle; - /*iterator it2 = */it += shift; + it += shift; if (!Rectangle::parse_detail(it, end, fun)) return false; - // it = it2; return true; } else if (identifier == "ellipse") { fun = new Ellipse; - /*iterator it2 = */it += shift; + it += shift; 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; } else if (identifier == "repeat") { fun = new Repeat; - /*iterator it2 = */it += shift; + it += shift; if (!Repeat::parse_detail(it, end, fun)) return false; - // it = it2; - return true; } else if (identifier == "rotate") { fun = new Rotate; - // iterator it2 = it += shift; if (!Rotate::parse_detail(it, end, fun)) return false; - // // it = it2; - return true; } else if (identifier == "translate") { fun = new Translate; - /*iterator it2 = */it += shift; + it += shift; if (!Translate::parse_detail(it, end, fun)) return false; - // it = it2; - return true; } else if (identifier == "shear") { fun = new Shear; - /*iterator it2 = */it += shift; + it += shift; if (!Shear::parse_detail(it, end, fun)) return false; - // it = it2; - return true; } else if (identifier == "union") { fun = new Union; - /*iterator it2 = */it += shift; + it += shift; if (!Union::parse_detail(it, end, fun)) return false; - // it = it2; - return true; } else if (identifier == "mask") { fun = new Mask; diff --git a/src/Classic/Utilities/MSLang.h b/src/Classic/Utilities/MSLang.h index bf005f461e658f51ddc14c09157af915f4257d35..445aa34c6c81879444c6fcefae0b8e05629a8146 100644 --- a/src/Classic/Utilities/MSLang.h +++ b/src/Classic/Utilities/MSLang.h @@ -1,3 +1,6 @@ +#ifndef MSLANG_H +#define MSLANG_H + #include "Algorithms/Vektor.h" #include "AppTypes/Tenzor.h" @@ -33,6 +36,13 @@ namespace mslang { 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 { if (2 * std::abs(X[0] - center_m[0]) <= width_m && 2 * std::abs(X[1] - center_m[1]) <= height_m) @@ -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 { int level_m; std::list<Base*> objects_m; @@ -210,4 +245,6 @@ namespace mslang { }; bool parse(std::string str, Function* &fun); -} \ No newline at end of file +} + +#endif \ No newline at end of file diff --git a/src/Classic/Utilities/Mesher.cpp b/src/Classic/Utilities/Mesher.cpp new file mode 100644 index 0000000000000000000000000000000000000000..f95caa4277bda9b0ecdaeca38a537e80cdc76ae0 --- /dev/null +++ b/src/Classic/Utilities/Mesher.cpp @@ -0,0 +1,218 @@ +#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 diff --git a/src/Classic/Utilities/Mesher.h b/src/Classic/Utilities/Mesher.h new file mode 100644 index 0000000000000000000000000000000000000000..696bf484c5d75dc03ddae4b96734b10fe9a7801d --- /dev/null +++ b/src/Classic/Utilities/Mesher.h @@ -0,0 +1,43 @@ +#ifndef MESHER_H +#define MESHER_H + +#include "Algorithms/Vektor.h" +#include "Utilities/MSLang.h" + +class Mesher { +public: + Mesher(std::vector<Vector_t> &vertices); + + std::vector<mslang::Triangle> getTriangles() const; + +private: + void orientVerticesCCW(); + bool isConvex(unsigned int i) const; + double crossProduct(const Vector_t &a, + const Vector_t &b) const; + bool isPointOnLine(unsigned int i, + unsigned int j, + const Vector_t &pt) const; + bool isPointRightOfLine(unsigned int i, + unsigned int j, + const Vector_t &pt) const; + bool lineSegmentTouchesOrCrossesLine(unsigned int i, + unsigned int j, + unsigned int k, + unsigned int l) const; + bool isPotentialEdgeIntersected(unsigned int i) const; + bool isPointInsideCone(unsigned int i, + unsigned int j, + unsigned int jPlusOne, + unsigned int jMinusOne) const; + bool isEar(unsigned int i) const; + std::vector<unsigned int> findAllEars() const; + double computeMinimumAngle(unsigned int i) const; + unsigned int selectBestEar(std::vector<unsigned int> &ears) const; + void apply(); + + std::vector<mslang::Triangle> triangles_m; + std::vector<Vector_t> vertices_m; +}; + +#endif \ No newline at end of file diff --git a/tools/mslang/mslang.cpp b/tools/mslang/mslang.cpp index fdd289a7c83dfa72c772bef1cb395bb7d741404c..16247f0886a5cb72c338dd574f141a01ff5a00f3 100644 --- a/tools/mslang/mslang.cpp +++ b/tools/mslang/mslang.cpp @@ -50,7 +50,6 @@ int main(int argc, char *argv[]) std::ofstream out("test.gpl"); for (mslang::Base* bfun: baseBlocks) { // bfun->print(0); - // std::cout << std::endl; bfun->computeBoundingBox(); bfun->writeGnuplot(out); }