lib/creaVtk/Utilities.cxx

   1 /*
   2 Copyright 2012-2022 Ronald Römer
   3
   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
   7
   8     http://www.apache.org/licenses/LICENSE-2.0
   9
  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.
  15 */
  16
  17 #include "Utilities.h"
  18
  19 #include <cmath>
  20
  21 #include <vtkPoints.h>
  22 #include <vtkIdList.h>
  23 #include <vtkMath.h>
  24 #include <vtkPolyData.h>
  25 #include <vtkDataWriter.h>
  26 #include <vtkSmartPointer.h>
  27
  28 #include <vtkPolyDataWriter.h>
  29
  30 void ComputeNormal (vtkPoints *pts, double *n, vtkIdType num, const vtkIdType *poly) {
  31     n[0] = 0; n[1] = 0; n[2] = 0;
  32
  33     if (num == 3) {
  34         double pA[3], pB[3], pC[3], a[3], b[3];
  35
  36         pts->GetPoint(poly[0], pA);
  37         pts->GetPoint(poly[1], pB);
  38         pts->GetPoint(poly[2], pC);
  39
  40         vtkMath::Subtract(pB, pA, a);
  41         vtkMath::Subtract(pC, pA, b);
  42
  43         vtkMath::Cross(a, b, n);
  44     } else {
  45         double pA[3], pB[3];
  46
  47         vtkIdType i, a, b;
  48
  49         for (i = 0; i < num; i++) {
  50             a = poly[i];
  51             b = poly[i+1 == num ? 0 : i+1];
  52
  53             pts->GetPoint(a, pA);
  54             pts->GetPoint(b, pB);
  55
  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]);
  59         }
  60     }
  61
  62     vtkMath::Normalize(n);
  63 }
  64
  65 void FindPoints (vtkKdTreePointLocator *pl, const double *pt, vtkIdList *pts, double tol) {
  66     pts->Reset();
  67
  68     vtkPolyData *pd = vtkPolyData::SafeDownCast(pl->GetDataSet());
  69
  70     vtkIdList *closest = vtkIdList::New();
  71
  72     // vtkKdTree.cxx#L2505
  73     // arbeitet mit single-precision
  74     pl->FindPointsWithinRadius(std::max(1e-3, tol), pt, closest);
  75
  76     vtkIdType i, numPts = closest->GetNumberOfIds();
  77
  78     double c[3], v[3];
  79
  80     for (i = 0; i < numPts; i++) {
  81         pd->GetPoint(closest->GetId(i), c);
  82         vtkMath::Subtract(pt, c, v);
  83
  84         if (vtkMath::Norm(v) < tol) {
  85             pts->InsertNextId(closest->GetId(i));
  86         }
  87     }
  88
  89     closest->Delete();
  90 }
  91
  92 void WriteVTK (const char *name, vtkPolyData *pd) {
  93     vtkPolyDataWriter *w = vtkPolyDataWriter::New();
  94     w->SetInputData(pd);
  95     w->SetFileName(name);
  96     w->Update();
  97     w->Delete();
  98 }
  99
 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
 102
 103     double _vA[3];
 104
 105     vtkMath::Cross(n, vA, _vA);
 106     double ang = std::atan2(vtkMath::Dot(_vA, vB), vtkMath::Dot(vA, vB));
 107
 108     if (ang < 0) {
 109         ang += 2*M_PI;
 110     }
 111
 112     return ang;
 113 }
 114
 115 Base::Base (vtkPoints *pts, vtkIdType num, const vtkIdType *poly) {
 116     ComputeNormal(pts, n, num, poly);
 117
 118     double ptA[3],
 119         ptB[3];
 120
 121     pts->GetPoint(poly[0], ptA);
 122     pts->GetPoint(poly[1], ptB);
 123
 124     ei[0] = ptB[0]-ptA[0];
 125     ei[1] = ptB[1]-ptA[1];
 126     ei[2] = ptB[2]-ptA[2];
 127
 128     vtkMath::Normalize(ei);
 129
 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];
 133
 134     vtkMath::Normalize(ej);
 135
 136     d = n[0]*ptA[0]+n[1]*ptA[1]+n[2]*ptA[2];
 137 }
 138
 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];
 142
 143     out[0] = x;
 144     out[1] = y;
 145 }
 146
 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];
 151
 152     out[0] = x;
 153     out[1] = y;
 154     out[2] = z;
 155 }
 156
 157 void ComputeNormal (const Poly &poly, double *n) {
 158     n[0] = 0; n[1] = 0; n[2] = 0;
 159
 160     Poly::const_iterator itrA, itrB;
 161
 162     for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
 163         itrB = itrA+1;
 164
 165         if (itrB == poly.end()) {
 166             itrB = poly.begin();
 167         }
 168
 169         const Point3d &ptA = *itrA,
 170             &ptB = *itrB;
 171
 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);
 175     }
 176
 177     vtkMath::Normalize(n);
 178 }
 179
 180 bool PointInPoly (const Poly &poly, const Point3d &p) {
 181     bool in = false;
 182
 183     Poly::const_iterator itrA, itrB;
 184
 185     for (itrA = poly.begin(); itrA != poly.end(); ++itrA) {
 186         itrB = itrA+1;
 187
 188         if (itrB == poly.end()) {
 189             itrB = poly.begin();
 190         }
 191
 192         const Point3d &ptA = *itrA,
 193             &ptB = *itrB;
 194
 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))) {
 198
 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) {
 201                 in = !in;
 202             }
 203         }
 204     }
 205
 206     return in;
 207 }
 208
 209 void WritePolys (const char *name, const PolysType &polys) {
 210     vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
 211
 212     vtkSmartPointer<vtkPolyData> pd = vtkSmartPointer<vtkPolyData>::New();
 213     pd->SetPoints(pts);
 214     pd->Allocate(1);
 215
 216     for (auto &poly : polys) {
 217         vtkSmartPointer<vtkIdList> cell = vtkSmartPointer<vtkIdList>::New();
 218
 219         for (auto &p : poly) {
 220             cell->InsertNextId(pts->InsertNextPoint(p.x, p.y, p.z));
 221         }
 222
 223         pd->InsertNextCell(VTK_POLYGON, cell);
 224     }
 225
 226     WriteVTK(name, pd);
 227 }
 228
 229 void GetPolys (const ReferencedPointsType &pts, const IndexedPolysType &indexedPolys, PolysType &polys) {
 230     for (const auto &poly : indexedPolys) {
 231         Poly newPoly;
 232
 233         for (auto &id : poly) {
 234             newPoly.push_back(pts.at(id));
 235         }
 236
 237         polys.push_back(std::move(newPoly));
 238     }
 239 }