2 Copyright 2012-2022 Ronald Römer
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
8 http://www.apache.org/licenses/LICENSE-2.0
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
17 #include "Utilities.h"
21 #include <vtkPoints.h>
22 #include <vtkIdList.h>
24 #include <vtkPolyData.h>
25 #include <vtkDataWriter.h>
26 #include <vtkSmartPointer.h>
28 #include <vtkPolyDataWriter.h>
30 void ComputeNormal (vtkPoints *pts, double *n, vtkIdType num, const vtkIdType *poly) {
31 n[0] = 0; n[1] = 0; n[2] = 0;
34 double pA[3], pB[3], pC[3], a[3], b[3];
36 pts->GetPoint(poly[0], pA);
37 pts->GetPoint(poly[1], pB);
38 pts->GetPoint(poly[2], pC);
40 vtkMath::Subtract(pB, pA, a);
41 vtkMath::Subtract(pC, pA, b);
43 vtkMath::Cross(a, b, n);
49 for (i = 0; i < num; i++) {
51 b = poly[i+1 == num ? 0 : i+1];
56 n[0] += (pA[1]-pB[1])*(pA[2]+pB[2]);
57 n[1] += (pA[2]-pB[2])*(pA[0]+pB[0]);
58 n[2] += (pA[0]-pB[0])*(pA[1]+pB[1]);
62 vtkMath::Normalize(n);
65 void FindPoints (vtkKdTreePointLocator *pl, const double *pt, vtkIdList *pts, double tol) {
68 vtkPolyData *pd = vtkPolyData::SafeDownCast(pl->GetDataSet());
70 vtkIdList *closest = vtkIdList::New();
72 // vtkKdTree.cxx#L2505
73 // arbeitet mit single-precision
74 pl->FindPointsWithinRadius(std::max(1e-3, tol), pt, closest);
76 vtkIdType i, numPts = closest->GetNumberOfIds();
80 for (i = 0; i < numPts; i++) {
81 pd->GetPoint(closest->GetId(i), c);
82 vtkMath::Subtract(pt, c, v);
84 if (vtkMath::Norm(v) < tol) {
85 pts->InsertNextId(closest->GetId(i));
92 void WriteVTK (const char *name, vtkPolyData *pd) {
93 vtkPolyDataWriter *w = vtkPolyDataWriter::New();
100 double GetAngle (double *vA, double *vB, double *n) {
101 // http://math.stackexchange.com/questions/878785/how-to-find-an-angle-in-range0-360-between-2-vectors
105 vtkMath::Cross(n, vA, _vA);
106 double ang = std::atan2(vtkMath::Dot(_vA, vB), vtkMath::Dot(vA, vB));
115 Base::Base (vtkPoints *pts, vtkIdType num, const vtkIdType *poly) {
116 ComputeNormal(pts, n, num, poly);
121 pts->GetPoint(poly[0], ptA);
122 pts->GetPoint(poly[1], ptB);
124 ei[0] = ptB[0]-ptA[0];
125 ei[1] = ptB[1]-ptA[1];
126 ei[2] = ptB[2]-ptA[2];
128 vtkMath::Normalize(ei);
130 ej[0] = n[1]*ei[2]-n[2]*ei[1];
131 ej[1] = -n[0]*ei[2]+n[2]*ei[0];
132 ej[2] = n[0]*ei[1]-n[1]*ei[0];
134 vtkMath::Normalize(ej);
136 d = n[0]*ptA[0]+n[1]*ptA[1]+n[2]*ptA[2];
139 void Transform (const double *in, double *out, const Base &base) {
140 double x = base.ei[0]*in[0]+base.ei[1]*in[1]+base.ei[2]*in[2],
141 y = base.ej[0]*in[0]+base.ej[1]*in[1]+base.ej[2]*in[2];
147 void BackTransform (const double *in, double *out, const Base &base) {
148 double x = in[0]*base.ei[0]+in[1]*base.ej[0]+base.d*base.n[0],
149 y = in[0]*base.ei[1]+in[1]*base.ej[1]+base.d*base.n[1],
150 z = in[0]*base.ei[2]+in[1]*base.ej[2]+base.d*base.n[2];
157 void ComputeNormal (const Poly &poly, double *n) {
158 n[0] = 0; n[1] = 0; n[2] = 0;
160 Poly::const_iterator itrA, itrB;
162 for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
165 if (itrB == poly.end()) {
169 const Point3d &ptA = *itrA,
172 n[0] += (ptA.y-ptB.y)*(ptA.z+ptB.z);
173 n[1] += (ptA.z-ptB.z)*(ptA.x+ptB.x);
174 n[2] += (ptA.x-ptB.x)*(ptA.y+ptB.y);
177 vtkMath::Normalize(n);
180 bool PointInPoly (const Poly &poly, const Point3d &p) {
183 Poly::const_iterator itrA, itrB;
185 for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
188 if (itrB == poly.end()) {
192 const Point3d &ptA = *itrA,
195 if ((ptA.x <= p.x || ptB.x <= p.x)
196 && ((ptA.y < p.y && ptB.y >= p.y)
197 || (ptB.y < p.y && ptA.y >= p.y))) {
199 // schnittpunkt mit bounding box und strahlensatz
200 if (ptA.x+(p.y-ptA.y)*(ptB.x-ptA.x)/(ptB.y-ptA.y) < p.x) {
209 void WritePolys (const char *name, const PolysType &polys) {
210 vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
212 vtkSmartPointer<vtkPolyData> pd = vtkSmartPointer<vtkPolyData>::New();
216 for (auto &poly : polys) {
217 vtkSmartPointer<vtkIdList> cell = vtkSmartPointer<vtkIdList>::New();
219 for (auto &p : poly) {
220 cell->InsertNextId(pts->InsertNextPoint(p.x, p.y, p.z));
223 pd->InsertNextCell(VTK_POLYGON, cell);
229 void GetPolys (const ReferencedPointsType &pts, const IndexedPolysType &indexedPolys, PolysType &polys) {
230 for (const auto &poly : indexedPolys) {
233 for (auto &id : poly) {
234 newPoly.push_back(pts.at(id));
237 polys.push_back(std::move(newPoly));