SHARE
TWEET

3D envelope

a guest Mar 21st, 2014 110 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <unordered_set>
  2. #include <unordered_map>
  3. #include <fstream>
  4. #include <iostream>
  5. #include <chrono>
  6. #include <memory>
  7.  
  8. using namespace std;
  9.  
  10. static string IN_FILE = "points_uniq.txt";
  11. static double NORM_ZERO = 1e-9;
  12. static double ZERO = 1e-15;
  13.  
  14. class Point3D {
  15.  
  16. public:
  17.         int id;
  18.  
  19. private:
  20.         double x;
  21.         double y;
  22.         double z;
  23.  
  24. public:
  25.         Point3D() : Point3D(0, 0, 0) {}
  26.         Point3D(int iid, double ix, double iy, double iz) : id(iid), x(ix), y(iy), z(iz) {}
  27.  
  28.         Point3D(double x, double y, double z) :Point3D(-1, x, y, z) {}
  29.  
  30.         bool operator== (const Point3D & x) const {
  31.                 return id == x.id;
  32.         }
  33.  
  34.         Point3D neg() const {
  35.                 return Point3D(-x, -y, -z);
  36.         }
  37.  
  38.         static Point3D add(initializer_list<Point3D> pts) {
  39.                 double x = 0;
  40.                 double y = 0;
  41.                 double z = 0;
  42.                 for (Point3D pt : pts) {
  43.                         x += pt.x;
  44.                         y += pt.y;
  45.                         z += pt.z;
  46.                 }
  47.                 return Point3D(x, y, z);
  48.         }
  49.  
  50.         static Point3D sub(const Point3D a, const Point3D b) {
  51.                 return Point3D(a.x - b.x, a.y - b.y, a.z - b.z);
  52.         }
  53.  
  54.         Point3D mul(double d) const {
  55.                 return Point3D(x * d, y * d, z * d);
  56.         }
  57.  
  58.         static double scal(const Point3D a, const Point3D b) {
  59.                 return a.x * b.x + a.y * b.y + a.z * b.z;
  60.         }
  61.  
  62.         static Point3D prod(const Point3D a, const Point3D b) {
  63.                 double x = a.y * b.z - a.z * b.y;
  64.                 double y = a.z * b.x - a.x * b.z;
  65.                 double z = a.x * b.y - a.y * b.x;
  66.  
  67.                 return Point3D(x, y, z);
  68.         }
  69.  
  70.         double size() const {
  71.                 double d = scal(*this, *this);
  72.                 return sqrt(d);
  73.         }
  74.  
  75.         Point3D norm() const {
  76.                 double d = size();
  77.                 if (d < ZERO)
  78.                         return *this;
  79.                 return mul(1 / d);
  80.         }
  81. };
  82.  
  83. class Line {
  84.  
  85.         Point3D a;
  86.         Point3D b;
  87.         bool mark;
  88.  
  89. public:
  90.         Line(Point3D ia, Point3D ib): a(ia), b(ib), mark(false) {}
  91.  
  92.         void doMark() {
  93.                 mark = !mark;
  94.         }
  95.  
  96.         bool isMark() const {
  97.                 return mark;
  98.         }
  99.  
  100.         Point3D beg() const {
  101.                 return a;
  102.         }
  103.  
  104.         Point3D end() const {
  105.                 return b;
  106.         }
  107.  
  108.         Point3D beg(bool rev) const {
  109.                 return rev ? b : a;
  110.         }
  111.  
  112.         Point3D end(bool rev) const {
  113.                 return rev ? a : b;
  114.         }
  115.  
  116.         Point3D dir(bool rev) const {
  117.                 Point3D d = Point3D::sub(b, a);
  118.                 return rev ? d.neg() : d;
  119.         }
  120.  
  121.         void addIDs(unordered_set<int> & pts) const {
  122.                 pts.insert(a.id);
  123.                 pts.insert(b.id);
  124.         }
  125. };
  126.  
  127. class Triangle {
  128.  
  129.         Point3D norm;
  130.         shared_ptr<Line> a;
  131.         shared_ptr<Line> b;
  132.         shared_ptr<Line> c;
  133.         bool arev;
  134.         bool brev;
  135.         bool crev;
  136.  
  137. public:
  138.         Triangle(Point3D a, Point3D b, Point3D c, Point3D * p) : Triangle(make_shared<Line>(a, b), make_shared<Line>(b, c), make_shared<Line>(c, a), p) {}
  139.  
  140.         Triangle(shared_ptr<Line> a, shared_ptr<Line> b, shared_ptr<Line> c, Point3D * p) : Triangle(a, b, c, p, false, false) {}
  141.  
  142.         Triangle(shared_ptr<Line> ia, shared_ptr<Line> ib, shared_ptr<Line> ic, Point3D * p, bool iarev, bool icrev) : a(ia), b(ib), c(ic), arev(iarev), crev(icrev), norm(0, 0, 0) {
  143.                 brev = a->end(arev) == b->end();
  144.  
  145.                 bool ba = a->beg(arev) == c->end(crev);
  146.                 bool bb = b->beg(brev) == a->end(arev);
  147.                 bool bc = c->beg(crev) == b->end(brev);
  148.                 if (!(ba && bb && bc))
  149.                         cout << "failed ori\n";
  150.  
  151.                 Point3D da = a->dir(arev);
  152.                 Point3D dc = c->dir(!crev);
  153.  
  154.                 Point3D n = Point3D::prod(da, dc);
  155.                 if (p != nullptr) {
  156.                         double d = Point3D::scal(n, Point3D::sub(*p, a->beg()));
  157.                         if (d < 0)
  158.                                 n = n.mul(-1);
  159.                         if (fabs(d) < NORM_ZERO)
  160.                                 cout << "inv norm " << d << "\n";
  161.                 }
  162.                 norm = n.norm();
  163.         }
  164.  
  165.         bool isValid() const {
  166.                 return Point3D::scal(norm, norm) > 1e-5;
  167.         }
  168.  
  169.         double getOri(const Point3D p) const {
  170.                 auto diff = Point3D::sub(p, a->beg());
  171.                 auto sc = Point3D::scal(norm, diff);
  172.                 return sc;
  173.         }
  174.  
  175.         void marklines() {
  176.                 a->doMark();
  177.                 b->doMark();
  178.                 c->doMark();
  179.         }
  180.  
  181.         void addMarked(list<shared_ptr<Line>> & bord) {
  182.                 if (a->isMark())
  183.                         bord.push_back(a);
  184.                 if (b->isMark())
  185.                         bord.push_back(b);
  186.                 if (c->isMark())
  187.                         bord.push_back(c);
  188.         }
  189.  
  190.         void addIDs(unordered_set<int> & pts) {
  191.                 a->addIDs(pts);
  192.                 b->addIDs(pts);
  193.                 c->addIDs(pts);
  194.         }
  195. };
  196.  
  197. class Envelope {
  198.  
  199.         Point3D in;
  200.         list<shared_ptr<Triangle>> tri;
  201.  
  202. public:
  203.         Envelope(Point3D iin, shared_ptr<Triangle> tr) : in(iin) {
  204.                 tri.push_back(tr);
  205.         }
  206.  
  207.         void addPoint(Point3D pt) {
  208.                 for (auto tr = tri.begin(); tr != tri.end();)
  209.                 {
  210.                         if ((*tr)->getOri(pt) < NORM_ZERO) {
  211.                                 tr = tri.erase(tr);
  212.                                 continue;
  213.                         }
  214.                         (*tr)->marklines();
  215.                         tr++;
  216.                 }
  217.  
  218.                 list<shared_ptr<Line>> bord;
  219.                 for (auto tr : tri)
  220.                         tr->addMarked(bord);
  221.  
  222.                 unordered_map<int, shared_ptr<Line>> map;
  223.                 for (auto lin : bord) {
  224.                         Point3D b = lin->beg(false);
  225.                         Point3D e = lin->end(false);
  226.  
  227.                         shared_ptr<Line> x = map[b.id];
  228.                         if (x == nullptr) {
  229.                                 x = make_shared<Line>(pt, b);
  230.                                 map[b.id] = x;
  231.                         }
  232.                         bool xrev = x->end() == pt;
  233.  
  234.                         shared_ptr<Line> z = map[e.id];
  235.                         if (z == nullptr) {
  236.                                 z = make_shared<Line>(e, pt);
  237.                                 map[e.id] = z;
  238.                         }
  239.                         bool zrev = z->beg() == pt;
  240.  
  241.                         shared_ptr<Triangle> tr = make_shared<Triangle>(x, lin, z, &in, xrev, zrev);
  242.                         if (tr->isValid())
  243.                                 tri.push_back(tr);
  244.                         else
  245.                                 cout << "inv triangle" << "\n";
  246.                         lin->doMark();
  247.                 }
  248.         }
  249.  
  250.         size_t getPtCnt() {
  251.                 unordered_set<int> pts;
  252.                 for (auto tr : tri)
  253.                         tr->addIDs(pts);
  254.  
  255.                 return pts.size();
  256.         }
  257.  
  258.         size_t getTriCnt() {
  259.                 return tri.size();
  260.         }
  261.  
  262. };
  263.  
  264. static int ptCnt;
  265.  
  266. Point3D getPoint(istream & sc) {
  267.         double x, y, z;
  268.         sc >> x;
  269.         sc >> y;
  270.         sc >> z;
  271.         if (sc.eof())
  272.                 return Point3D();
  273.         return Point3D(++ptCnt, x, y, z);
  274. }
  275.  
  276. Envelope * initEnvelope(istream & sc) {
  277.         Point3D a = getPoint(sc);
  278.         Point3D b = getPoint(sc);
  279.  
  280.         vector<Point3D> pts;
  281.  
  282.         Point3D c(0,0,0);
  283.  
  284.         unique_ptr<Triangle> tr = nullptr;
  285.         while (!sc.eof()) {
  286.                 c = getPoint(sc);
  287.                 tr = make_unique<Triangle>(a, b, c, nullptr);
  288.                 if (tr->isValid())
  289.                         break;
  290.                 pts.push_back(c);
  291.         }
  292.  
  293.         // TODO: if !tr valid
  294.         Point3D d(0,0,0);
  295.         while (!sc.eof()) {
  296.                 d = getPoint(sc);
  297.                 double n = tr->getOri(d);
  298.                 if (fabs(n) > NORM_ZERO)
  299.                         break;
  300.                 pts.push_back(d);
  301.         }
  302.  
  303.         Point3D in = Point3D::add({ a, b, c, d }).mul(1 / 4.0);
  304.         Envelope * env = new Envelope(in, make_shared<Triangle>(a, b, c, &in));
  305.         env->addPoint(d);
  306.  
  307.         for each (Point3D pt in pts)
  308.                 env->addPoint(pt);
  309.         return env;
  310. }
  311.  
  312. using namespace std::chrono;
  313.  
  314. int main(char ** arg) {
  315.         ptCnt = 0;
  316.  
  317.         for (int i = 0; i < 20; i++) {
  318.                 ifstream sc(IN_FILE);
  319.  
  320.                 auto start = high_resolution_clock::now();
  321.                 Envelope * env = initEnvelope(sc);
  322.  
  323.                 while (!sc.eof()) {
  324.                         Point3D pt = getPoint(sc);
  325.                         if (pt.id >= 0)
  326.                                 env->addPoint(pt);
  327.                 }
  328.                 auto end = high_resolution_clock::now();
  329.                 sc.close();
  330.  
  331.                 auto diff = duration_cast<duration<double>>(end - start);
  332.  
  333.                 cout << "time: " << diff.count() << "\n";
  334.                 cout << "poly cnt: " << env->getTriCnt() << " pt cnt: " << ptCnt << " pts: " << env->getPtCnt() << "\n";
  335.  
  336.                 delete env;
  337.         }
  338.     return 0;
  339. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top