API tools faq
paste
anutka

Untitled

anutka
Nov 9th, 2020
86
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 17.83 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <cmath>
  3. #include <vector>
  4. #include <assert.h>
  5. #include <iomanip>
  6. #include <queue>
  7. #include<algorithm>
  8. const double pi = 3.141592653589793238462643;
  9. const double accuracy = 0.0;
  10. const double MAX = 10000000;
  11. class Vector;
  12. class Vector2d {
  13. public:
  14. long double x;
  15. long double y;
  16. Vector2d() = default;
  17. Vector2d(long double x,long double y) : x(x), y(y){}
  18. Vector2d(const Vector2d& vector) = default;
  19.  
  20. friend bool operator == (const Vector2d& left, const Vector2d& right);
  21. friend bool operator != (const Vector2d& left, const Vector2d& right);
  22. friend Vector2d operator + (const Vector2d& left, const Vector2d& right);
  23. friend bool operator < (const Vector2d& left, const Vector2d& right);
  24. friend bool operator > (const Vector2d& left, const Vector2d& right);
  25. friend bool operator >= (const Vector2d& left, const Vector2d& right);
  26. friend bool operator <= (const Vector2d& left, const Vector2d& right);
  27. friend Vector2d operator - (const Vector2d& left, const Vector2d& right);
  28. friend Vector2d operator / (const Vector2d& left, const long double & right);
  29.  
  30. void operator += (const Vector2d& other){
  31. x += other.x;
  32. y += other.y;
  33. }
  34. void operator -= (const Vector2d& other){
  35. x -= other.x;
  36. y -= other.y;
  37. }
  38. void operator /= (const long double & c){
  39. assert(c >0);
  40. x /= c;
  41. y /= c;
  42. }
  43. void operator *= (const long double & c) {
  44. x *= c;
  45. y *= c;
  46. }
  47. };
  48. long double distance (const Vector2d& first, const Vector2d& second) {
  49. return sqrtl(((first.x-second.x)*(first.x-second.x)
  50. + (first.y-second.y)*(first.y-second.y)));
  51. }
  52.  
  53. long double scalarMul (const Vector2d& first, const Vector2d& second) {
  54. return first.x * second.x + first.y * second.y;
  55. }
  56.  
  57. /*
  58. * Арифметические операции
  59. */
  60. bool operator < (const Vector2d& left, const Vector2d& right) {
  61. if (right == left) {
  62. return false;
  63. }
  64. if (right == Vector2d(0,0)) {
  65. return false;
  66. }
  67. if (left == Vector2d(0,0)) {
  68. return true;
  69. }
  70. Vector2d a(0, -1);
  71. long double x = scalarMul(a, left) * sqrtl(256*(right.x* right.x + right.y*right.y))/16 -
  72. scalarMul(a, right) * sqrtl(256*(left.x* left.x + left.y*left.y))/16;
  73. if (x == 0) {
  74. return left.x* left.x + left.y*left.y < right.x* right.x + right.y*right.y;
  75. } else {
  76. return x > -accuracy;
  77. }
  78. }
  79. bool operator <= (const Vector2d& left, const Vector2d& right) {
  80. return left < right || left == right;
  81. }
  82. bool operator > (const Vector2d& left, const Vector2d& right) {
  83. return !(left <= right);
  84. }
  85. bool operator >= (const Vector2d& left, const Vector2d& right) {
  86. return left > right;
  87. }
  88. Vector2d operator + (const Vector2d& left, const Vector2d& right){
  89. Vector2d data = left;
  90. data += right;
  91. return data;
  92. }
  93. Vector2d operator - (const Vector2d& left, const Vector2d& right){
  94. Vector2d data = left;
  95. data -= right;
  96. return data;
  97. }
  98. Vector2d operator / (const Vector2d& left, const long double & right){
  99. assert(right != 0);
  100. Vector2d data = left;
  101. data /= right;
  102. return data;
  103. }
  104. Vector2d operator * (const Vector2d& left, const long double right){
  105. Vector2d data = left;
  106. data *= right;
  107. return data;
  108. }
  109.  
  110. bool operator == (const Vector2d& left, const Vector2d& right) {
  111. return fabsl(left.x - right.x) <= accuracy && fabsl(left.y - right.y) <= accuracy;
  112. }
  113. bool operator != (const Vector2d& left, const Vector2d& right) {
  114. return !(left == right);
  115. }
  116. long double OrientalSquare (const Vector2d& first, const Vector2d& second) {
  117. return first.x*second.y - first.y*second.x;
  118. }
  119.  
  120. size_t FindMinVertice(std::vector<std::pair<Vector2d, size_t>> & vertices) {
  121. size_t minInd = 0;
  122. Vector2d & min = vertices[minInd].first;
  123. for (size_t i = 0; i < vertices.size(); ++i) {
  124. Vector2d & cur = vertices[i].first;
  125. if (cur.x < min.x || (cur.x == min.x && cur.y < min.y)) {
  126. min = cur;
  127. minInd = i;
  128. }
  129. }
  130. return minInd;
  131. }
  132. bool isOnTheSameSide(Vector2d & begin, Vector2d& end, Vector2d& p1, Vector2d& p2) {
  133. Vector2d side = begin - end;
  134. Vector2d v1 = p1 - end;
  135. Vector2d v2 = p2 - end;
  136. long double s1 = OrientalSquare(side, v1);
  137. long double s2 = OrientalSquare(side, v2);
  138. return ((s1 >= 0 && s2 >= 0) || (s1 <= 0 && s2 <= 0));
  139. }
  140.  
  141. bool isOnTheSameLine(Vector2d & x, Vector2d & y, Vector2d & z) {
  142. Vector2d v1 = y - x, v2 = z - x;
  143. return OrientalSquare(v1, v2) == 0;
  144. }
  145.  
  146. class Vector {
  147. size_t num;
  148. public:
  149. double x;
  150. double y;
  151. double z;
  152. Vector() = default;
  153. Vector(double x,double y, double z) : x(x), y(y), z(z), num(-1){}
  154. Vector(const Vector& vector) = default;
  155.  
  156. friend bool operator == (const Vector& left, const Vector& right);
  157. friend bool operator != (const Vector& left, const Vector& right);
  158. friend Vector operator + (const Vector& left, const Vector& right);
  159. friend bool operator < (const Vector& left, const Vector& right);
  160. friend bool operator > (const Vector& left, const Vector& right);
  161. friend bool operator >= (const Vector& left, const Vector& right);
  162. friend bool operator <= (const Vector& left, const Vector& right);
  163. friend Vector operator - (const Vector& left, const Vector& right);
  164. friend Vector operator / (const Vector& left, const double& right);
  165.  
  166. void operator += (const Vector& other){
  167. x += other.x;
  168. y += other.y;
  169. z += other.z;
  170. }
  171. void operator -= (const Vector& other){
  172. x -= other.x;
  173. y -= other.y;
  174. z -= other.z;
  175. }
  176. void operator /= (const double& c){
  177. assert(c >0);
  178. x /= c;
  179. y /= c;
  180. z /= c;
  181. }
  182. void operator *= (const double& c) {
  183. x *= c;
  184. y *= c;
  185. z *= c;
  186. }
  187. size_t getNum() {
  188. return num;
  189. }
  190. void setNum(size_t i) {
  191. num = i;
  192. }
  193. };
  194.  
  195. /*
  196. * Арифметические операции
  197. */
  198. Vector operator + (const Vector& left, const Vector& right){
  199. Vector data = left;
  200. data += right;
  201. return data;
  202. }
  203. Vector operator - (const Vector& left, const Vector& right){
  204. Vector data = left;
  205. data -= right;
  206. return data;
  207. }
  208. Vector operator / (const Vector& left, const double& right){
  209. assert(right != 0);
  210. Vector data = left;
  211. data /= right;
  212. return data;
  213. }
  214. Vector operator * (const Vector& left, const double right){
  215. Vector data = left;
  216. data *= right;
  217. return data;
  218. }
  219.  
  220. /*
  221. * Операции сравнения
  222. */
  223. bool operator > (const Vector& left, const Vector& right) {
  224. return (left.x>right.x) && (left.y > right.y) && (left.z > right.z);
  225. }
  226. bool operator < (const Vector& left, const Vector& right) {
  227. return (left.x<right.x) && (left.y < right.y) && (left.z < right.z);
  228. }
  229. bool operator <= (const Vector& left, const Vector& right) {
  230. return (left.x<=right.x) && (left.y <= right.y) && (left.z <= right.z);
  231. }
  232. bool operator >= (const Vector& left, const Vector& right) {
  233. return (left.x>=right.x) && (left.y >= right.y) && (left.z >= right.z);
  234. }
  235.  
  236. bool operator == (const Vector& left, const Vector& right) {
  237. return fabs(left.x - right.x) <= accuracy && fabs(left.y - right.y) <= accuracy;
  238. }
  239. bool operator != (const Vector& left, const Vector& right) {
  240. return !(left == right);
  241. }
  242.  
  243.  
  244. double scalarMul (const Vector& first, const Vector& second) {
  245. return first.x * second.x + first.y * second.y + first.z * second.z;
  246. }
  247. Vector vectorMul (const Vector& first, const Vector& second) {
  248. return Vector(first.y*second.z - first.z*second.y,
  249. first.z*second.x - first.x*second.z,
  250. first.x*second.y - first.y*second.x);
  251. }
  252. double distance (const Vector& first, const Vector& second) {
  253. return sqrt( 256 *((first.x-second.x)*(first.x-second.x)
  254. + (first.y-second.y)*(first.y-second.y) +
  255. (first.z-second.z)*(first.z-second.z)))/16;
  256. }
  257.  
  258. double Volume(const Vector& first, const Vector& second, const Vector& third) {
  259. return first.x*second.y*third.z + first.y*second.z*third.x + first.z*second.x*third.y -
  260. first.x*second.z*third.y + first.y*second.x*third.z + first.z*second.y*third.x;
  261. }
  262.  
  263. class Plane {
  264. public:
  265. Vector normal;
  266. double D;
  267. Plane (const Vector& firstP, const Vector& secondP, const Vector& thirdP):
  268. normal(vectorMul(firstP - thirdP, secondP - thirdP)),
  269. D(-scalarMul(normal, thirdP)) {}
  270.  
  271. };
  272.  
  273. class Face {
  274. private:
  275. std::vector<Vector> vertices;
  276. std::vector<Vector> ConvexHull;
  277.  
  278. Plane plane;
  279. bool isCounted = false;
  280. bool isXConst = false;
  281. bool isZConst = false;
  282. bool isYConst = false;
  283. void FindConvexHull() {
  284. std::vector<std::pair<Vector2d, size_t>> verticesProection;
  285. if (!isZConst) {
  286. for(size_t i =0; i < vertices.size(); i++) {
  287. verticesProection.emplace_back(Vector2d(vertices[i].x, vertices[i].y), i);
  288. }
  289. } else {
  290. if (!isXConst) {
  291. for(size_t i =0; i < vertices.size(); i++) {
  292. verticesProection.emplace_back(Vector2d(vertices[i].y, vertices[i].z), i);
  293. }
  294. } else {
  295. if(!isYConst) {
  296. for(size_t i =0; i < vertices.size(); i++) {
  297. verticesProection.emplace_back(Vector2d(vertices[i].x, vertices[i].z), i);
  298. }
  299. }
  300. }
  301. }
  302. std::vector<std::pair<Vector2d, size_t>> vectors(vertices.size());
  303. std::vector<Vector2d> result;
  304. size_t minInd = FindMinVertice(verticesProection);
  305. Vector2d min2D = verticesProection[minInd].first;
  306. Vector min3D = vertices[verticesProection[minInd].second];
  307. for (auto &i: verticesProection) {
  308. vectors.emplace_back(i.first - min2D, i.second);
  309. }
  310. std::sort(vectors.begin(), vectors.end());
  311. Vector2d prev = min2D;
  312. ConvexHull.push_back(min3D);
  313. for (auto& i:vectors) {
  314. Vector2d nextVertice2d = min2D + i.first;
  315. Vector nextVertice3d = min3D + vertices[i.second];
  316. if (nextVertice2d != prev) {
  317. if (result.size() < 2) {
  318. result.push_back(nextVertice2d);
  319. ConvexHull.push_back(nextVertice3d);
  320. } else {
  321. while(result.size() > 1 && isOnTheSameLine(min2D, result[result.size() - 1], nextVertice2d) ){
  322. result.pop_back();
  323. ConvexHull.pop_back();
  324. }
  325. size_t j = result.size() - 1;
  326. while(j > 0 && (isOnTheSameLine(result[j], nextVertice2d, result[(j-1)]) || !isOnTheSameSide(result[j], nextVertice2d,
  327. result[(j-1)], result[(j + 1)% result.size()]))) {
  328. --j;
  329. }
  330. for (size_t t = result.size() - 1; t > j; --t) {
  331. result.pop_back();
  332. ConvexHull.pop_back();
  333. }
  334. result.push_back(nextVertice2d);
  335. ConvexHull.push_back(nextVertice3d);
  336. }
  337. }
  338. prev = nextVertice2d;
  339. }
  340. return;
  341. }
  342. public:
  343. Face() : plane(Vector(0,0,0), Vector(0,0,0), Vector(0,0,0)), vertices(0), ConvexHull(0){};
  344. Face(const Vector& a, const Vector& b, const Vector& c) : plane(a, b, c) {
  345. vertices.push_back(a);
  346. vertices.push_back(b);
  347. vertices.push_back(c);
  348. if (a.x == b.x && b.x == c.x) {
  349. isXConst = true;
  350. }
  351. if (a.y == b.y && b.y == c.y) {
  352. isYConst = true;
  353. }
  354. if (a.z == b.z && b.z == c.z) {
  355. isZConst = true;
  356. }
  357. }
  358.  
  359. Face(std::vector<Vector> & v) : vertices(v) , plane(v[0], v[1%v.size()], v[2%v.size()]) {
  360. Vector a = v[0];
  361. Vector b = v[1];
  362. Vector c = v[2];
  363. vertices.push_back(a);
  364. vertices.push_back(b);
  365. vertices.push_back(c);
  366. if (a.x == b.x && b.x == c.x) {
  367. isXConst = true;
  368. }
  369. if (a.y == b.y && b.y == c.y) {
  370. isYConst = true;
  371. }
  372. if (a.z == b.z && b.z == c.z) {
  373. isZConst = true;
  374. }
  375. }
  376. Vector ThirdVertex(const Vector& a, const Vector& b) const {
  377. for (auto&i : vertices) {
  378. if (i != a && i != b) {
  379. return i;
  380. }
  381. }
  382. assert(vertices.size() >= 3);
  383. return a;
  384. }
  385. Vector getNormal() const {
  386. return plane.normal;
  387. }
  388. void push_back(const Vector& a) {
  389. vertices.push_back(a);
  390. }
  391. std::vector<Vector> getConvexHull() {
  392. if (isCounted) {
  393. return ConvexHull;
  394. }
  395. FindConvexHull();
  396. isCounted = true;
  397. return ConvexHull;
  398. }
  399.  
  400. };
  401.  
  402. class Edge {
  403. public:
  404. Vector begin;
  405. Vector end;
  406. const Face& face;
  407. Edge(Vector begin, Vector end, const Face& face) :begin(begin), end(end), face(face){
  408. }
  409. };
  410. bool isOnOnePlane(Vector& a, Vector& b, Vector& c, Vector &d) {
  411. return fabs(Volume(b - a, c - a, d - a)) < accuracy;
  412. }
  413.  
  414.  
  415. class PointsContainer {
  416. std::vector<Vector> points;
  417. std::vector<Face> convexHull;
  418. public:
  419. PointsContainer() : points(0), convexHull(0){};
  420.  
  421. void emplace_back (double x, double y, double z) {
  422. points.emplace_back(x, y, z);
  423. size_t last = points.size() - 1;
  424. points[last].setNum(last);
  425. }
  426. Vector FindMaxAngle(Edge& e) {
  427. const Face& face = e.face;
  428. Vector thirdVertex = face.ThirdVertex(e.begin, e.end);
  429. Vector curMin = points[0];
  430. Vector curNormal = Plane(curMin, e.begin, e.end).normal;
  431. for (auto& point: points) {
  432. if (point != e.begin && point != e.end && !isOnOnePlane(e.begin, e.end, point, thirdVertex)) {
  433. Plane cur = Plane(point, e.begin, e.end);
  434. if (scalarMul(cur.normal, face.getNormal()) * distance(curNormal, Vector(0,0,0)) -
  435. scalarMul(curNormal, face.getNormal()) * distance(cur.normal, Vector(0,0,0)) > -accuracy){
  436. curNormal = cur.normal;
  437. curMin = point;
  438. }
  439. }
  440. }
  441. return curMin;
  442. }
  443. Face constructFaceWithAllPoints(Vector& a, Vector& b, Vector& c) {
  444. Face res({a, b, c});
  445. for (auto& i:points) {
  446. if (i!=a && i!= b && i!=c) {
  447. if (isOnOnePlane(a, b, c, i)) {
  448. res.push_back(i);
  449. }
  450. }
  451. }
  452. return res;
  453. }
  454. void FindFirstFace() {
  455. Vector& first = points[0];
  456. for (size_t i = 0; i < points.size(); ++i) {
  457. Vector & cur = points[i];
  458. if (cur.x < first.x || (cur.x == first.x && cur.y < first.y) || (cur.x == first.x && cur.y == first.y && cur.z < first.z)) {
  459. first = cur;
  460. }
  461. }
  462. Vector& second = points[1];
  463. for (size_t i = 0; i < points.size(); ++i) {
  464. Vector & cur = points[i];
  465. if (cur != first) {
  466. if (scalarMul(cur - first, Vector(1,0,0)) * distance(second, Vector(0,0,0)) -
  467. scalarMul(second - first, Vector(1, 0, 0)) * distance(cur, Vector(0,0,0)) < accuracy) {
  468. second = cur;
  469. }
  470. }
  471. }
  472. Vector third = first;
  473. bool isInicilized = false;
  474. for (size_t i = 0; i < points.size(); ++i) {
  475. Vector & cur = points[i];
  476. if (cur != first && vectorMul(cur-first, cur-second) != Vector(0,0,0)) {
  477. if (!isInicilized ||(scalarMul(cur - first, Vector(1,0,0)) * distance(second, Vector(0,0,0)) -
  478. scalarMul(second - first, Vector(1, 0, 0)) * distance(cur, Vector(0,0,0)) < accuracy)) {
  479. third = cur;
  480. isInicilized = true;
  481. }
  482. }
  483. }
  484. convexHull.push_back(constructFaceWithAllPoints(first, second, third));
  485. }
  486. void FindConvexHull () {
  487. std::vector<std::vector<bool>> isUsed(points.size(), std::vector<bool>(points.size(), 0));
  488. std::queue<Edge> edgesQ;
  489. FindFirstFace();
  490. std::vector<Vector> FaceConvHull = convexHull[0].getConvexHull();
  491. for (size_t i = 0; i < FaceConvHull.size() + 1; i++) {
  492. edgesQ.push(Edge(FaceConvHull[i],FaceConvHull[(i+1)%FaceConvHull.size()], convexHull[0]));
  493. }
  494. while (!edgesQ.empty()) {
  495. Edge e = edgesQ.front();
  496. edgesQ.pop();
  497. if (!isUsed[e.begin.getNum()][e.end.getNum()]) {
  498. isUsed[e.begin.getNum()][e.end.getNum()] = true;
  499. isUsed[e.end.getNum()] [e.begin.getNum()]= true;
  500. Vector p3 = FindMaxAngle(e);
  501. convexHull.push_back(constructFaceWithAllPoints(e.begin, e.end, p3));
  502. size_t last = convexHull.size() - 1;
  503. std::vector<Vector> FaceConvHull = convexHull[last].getConvexHull();
  504. for (size_t i = 0; i < FaceConvHull.size() + 1; i++) {
  505. edgesQ.push(Edge(FaceConvHull[i],FaceConvHull[(i+1)%FaceConvHull.size()], convexHull[last]));
  506. }
  507. }
  508. }
  509. }
  510. };
  511.  
  512.  
  513.  
  514. int main() {
  515. PointsContainer v;
  516. double x, y, z;
  517. int N;
  518. std::cin >> N;
  519. for (int i = 0; i < N; i++) {
  520. std::cin >> x>> y>> z;
  521. v.emplace_back(x, y, z);
  522. }
  523. v.FindConvexHull();
  524. return 0;
  525. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
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!