FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2012 Imetric 3D GmbH                                    *
3
 *                                                                         *
4
 *   This file is part of the FreeCAD CAx development system.              *
5
 *                                                                         *
6
 *   This library is free software; you can redistribute it and/or         *
7
 *   modify it under the terms of the GNU Library General Public           *
8
 *   License as published by the Free Software Foundation; either          *
9
 *   version 2 of the License, or (at your option) any later version.      *
10
 *                                                                         *
11
 *   This library  is distributed in the hope that it will be useful,      *
12
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
13
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
14
 *   GNU Library General Public License for more details.                  *
15
 *                                                                         *
16
 *   You should have received a copy of the GNU Library General Public     *
17
 *   License along with this library; see the file COPYING.LIB. If not,    *
18
 *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
19
 *   Suite 330, Boston, MA  02111-1307, USA                                *
20
 *                                                                         *
21
 ***************************************************************************/
22

23
#include "PreCompiled.h"
24
#ifndef _PreComp_
25
#include <algorithm>
26
#include <functional>
27
#endif
28

29
#include <QFuture>
30
#include <QFutureWatcher>
31
#include <QtConcurrentMap>
32

33
#include <Base/Sequencer.h>
34
#include <Base/Tools.h>
35

36
// #define OPTIMIZE_CURVATURE
37
#ifdef OPTIMIZE_CURVATURE
38
#include <Eigen/Eigenvalues>
39
#else
40
#include <Mod/Mesh/App/WildMagic4/Wm4MeshCurvature.h>
41
#endif
42

43
#include "Approximation.h"
44
#include "Curvature.h"
45
#include "Iterator.h"
46
#include "MeshKernel.h"
47
#include "Tools.h"
48

49

50
using namespace MeshCore;
51
namespace sp = std::placeholders;
52

53
MeshCurvature::MeshCurvature(const MeshKernel& kernel)
54
    : myKernel(kernel)
55
    , myMinPoints(20)
56
    , myRadius(0.5f)
57
{
58
    mySegment.resize(kernel.CountFacets());
59
    std::generate(mySegment.begin(), mySegment.end(), Base::iotaGen<FacetIndex>(0));
60
}
61

62
MeshCurvature::MeshCurvature(const MeshKernel& kernel, std::vector<FacetIndex> segm)
63
    : myKernel(kernel)
64
    , myMinPoints(20)
65
    , myRadius(0.5f)
66
    , mySegment(std::move(segm))
67
{}
68

69
void MeshCurvature::ComputePerFace(bool parallel)
70
{
71
    myCurvature.clear();
72
    MeshRefPointToFacets search(myKernel);
73
    FacetCurvature face(myKernel, search, myRadius, myMinPoints);
74

75
    if (!parallel) {
76
        Base::SequencerLauncher seq("Curvature estimation", mySegment.size());
77
        for (FacetIndex it : mySegment) {
78
            CurvatureInfo info = face.Compute(it);
79
            myCurvature.push_back(info);
80
            seq.next();
81
        }
82
    }
83
    else {
84
        // NOLINTBEGIN
85
        QFuture<CurvatureInfo> future =
86
            QtConcurrent::mapped(mySegment, std::bind(&FacetCurvature::Compute, &face, sp::_1));
87
        // NOLINTEND
88
        QFutureWatcher<CurvatureInfo> watcher;
89
        watcher.setFuture(future);
90
        watcher.waitForFinished();
91
        for (const auto& it : future) {
92
            myCurvature.push_back(it);
93
        }
94
    }
95
}
96

97
#ifdef OPTIMIZE_CURVATURE
98
namespace MeshCore
99
{
100
void GenerateComplementBasis(Eigen::Vector3f& rkU, Eigen::Vector3f& rkV, const Eigen::Vector3f& rkW)
101
{
102
    float fInvLength;
103

104
    if (fabs(rkW[0]) >= fabs(rkW[1])) {
105
        // W.x or W.z is the largest magnitude component, swap them
106
        fInvLength = 1.0 / sqrt(rkW[0] * rkW[0] + rkW[2] * rkW[2]);
107
        rkU[0] = -rkW[2] * fInvLength;
108
        rkU[1] = 0.0;
109
        rkU[2] = +rkW[0] * fInvLength;
110
        rkV[0] = rkW[1] * rkU[2];
111
        rkV[1] = rkW[2] * rkU[0] - rkW[0] * rkU[2];
112
        rkV[2] = -rkW[1] * rkU[0];
113
    }
114
    else {
115
        // W.y or W.z is the largest magnitude component, swap them
116
        fInvLength = 1.0 / sqrt(rkW[1] * rkW[1] + rkW[2] * rkW[2]);
117
        rkU[0] = 0.0;
118
        rkU[1] = +rkW[2] * fInvLength;
119
        rkU[2] = -rkW[1] * fInvLength;
120
        rkV[0] = rkW[1] * rkU[2] - rkW[2] * rkU[1];
121
        rkV[1] = -rkW[0] * rkU[2];
122
        rkV[2] = rkW[0] * rkU[1];
123
    }
124
}
125
}  // namespace MeshCore
126

127
void MeshCurvature::ComputePerVertex()
128
{
129
    // get all points
130
    const MeshPointArray& pts = myKernel.GetPoints();
131

132
    MeshCore::MeshRefPointToFacets pt2f(myKernel);
133
    MeshCore::MeshRefPointToPoints pt2p(myKernel);
134
    unsigned long numPoints = myKernel.CountPoints();
135

136
    myCurvature.clear();
137
    myCurvature.reserve(numPoints);
138

139
    std::vector<Eigen::Vector3f> akNormal(numPoints);
140
    std::vector<Eigen::Vector3f> akVertex(numPoints);
141
    for (unsigned long i = 0; i < numPoints; i++) {
142
        Base::Vector3f n = pt2f.GetNormal(i);
143
        akNormal[i][0] = n.x;
144
        akNormal[i][1] = n.y;
145
        akNormal[i][2] = n.z;
146
        const Base::Vector3f& p = pts[i];
147
        akVertex[i][0] = p.x;
148
        akVertex[i][1] = p.y;
149
        akVertex[i][2] = p.z;
150
    }
151

152
    // One could iterate over the triangles and then for each vertex of a triangle compute the
153
    // derivates. One could also iterate over the points and then for each adjacent point calculate
154
    // the derivates. Both methods must lead to the same values in the above matrices.
155
    //
156
    // Iterate over the vertexes
157
    for (unsigned long i = 0; i < numPoints; i++) {
158
        Eigen::Matrix3f akDNormal;
159
        akDNormal.setZero();
160
        Eigen::Matrix3f akWWTrn;
161
        akWWTrn.setZero();
162
        Eigen::Matrix3f akDWTrn;
163
        akDWTrn.setZero();
164

165
        int iV0 = i;
166
        int iV1;
167
        const std::set<unsigned long>& nb = pt2p[i];
168
        for (std::set<unsigned long>::const_iterator it = nb.begin(); it != nb.end(); ++it) {
169
            iV1 = *it;
170

171
            // Compute edge from V0 to V1, project to tangent plane of vertex,
172
            // and compute difference of adjacent normals.
173
            Eigen::Vector3f kE = akVertex[iV1] - akVertex[iV0];
174
            Eigen::Vector3f kW = kE - (kE.dot(akNormal[iV0])) * akNormal[iV0];
175
            Eigen::Vector3f kD = akNormal[iV1] - akNormal[iV0];
176
            for (int iRow = 0; iRow < 3; iRow++) {
177
                for (int iCol = 0; iCol < 3; iCol++) {
178
                    akWWTrn(iRow, iCol) += 2 * kW[iRow] * kW[iCol];
179
                    akDWTrn(iRow, iCol) += 2 * kD[iRow] * kW[iCol];
180
                }
181
            }
182
        }
183

184
        // Add in N*N^T to W*W^T for numerical stability.  In theory 0*0^T gets
185
        // added to D*W^T, but of course no update needed in the implementation.
186
        // Compute the matrix of normal derivatives.
187
        for (int iRow = 0; iRow < 3; iRow++) {
188
            for (int iCol = 0; iCol < 3; iCol++) {
189
                akWWTrn(iRow, iCol) =
190
                    0.5 * akWWTrn(iRow, iCol) + akNormal[i][iRow] * akNormal[i][iCol];
191
                akDWTrn(iRow, iCol) *= 0.5;
192
            }
193
        }
194

195
        akDNormal = akDWTrn * akWWTrn.inverse();
196

197
        // If N is a unit-length normal at a vertex, let U and V be unit-length
198
        // tangents so that {U, V, N} is an orthonormal set.  Define the matrix
199
        // J = [U | V], a 3-by-2 matrix whose columns are U and V.  Define J^T
200
        // to be the transpose of J, a 2-by-3 matrix.  Let dN/dX denote the
201
        // matrix of first-order derivatives of the normal vector field.  The
202
        // shape matrix is
203
        //   S = (J^T * J)^{-1} * J^T * dN/dX * J = J^T * dN/dX * J
204
        // where the superscript of -1 denotes the inverse.  (The formula allows
205
        // for J built from non-perpendicular vectors.) The matrix S is 2-by-2.
206
        // The principal curvatures are the eigenvalues of S.  If k is a principal
207
        // curvature and W is the 2-by-1 eigenvector corresponding to it, then
208
        // S*W = k*W (by definition).  The corresponding 3-by-1 tangent vector at
209
        // the vertex is called the principal direction for k, and is J*W.
210
        // compute U and V given N
211
        float minCurvature;
212
        float maxCurvature;
213
        Base::Vector3f minDirection;
214
        Base::Vector3f maxDirection;
215

216
        Eigen::Vector3f kU, kV;
217
        Eigen::Vector3f kN = akNormal[i];
218
        float len = kN.squaredNorm();
219
        if (len == 0) {
220
            continue;  // skip
221
        }
222
        MeshCore::GenerateComplementBasis(kU, kV, kN);
223

224
        // Compute S = J^T * dN/dX * J.  In theory S is symmetric, but
225
        // because we have estimated dN/dX, we must slightly adjust our
226
        // calculations to make sure S is symmetric.
227
        float fS01 = kU.dot(akDNormal * kV);
228
        float fS10 = kV.dot(akDNormal * kU);
229
        float fSAvr = 0.5 * (fS01 + fS10);
230
        Eigen::Matrix2f kS;
231
        kS(0, 0) = kU.dot(akDNormal * kU);
232
        kS(0, 1) = fSAvr;
233
        kS(1, 0) = fSAvr;
234
        kS(1, 1) = kV.dot(akDNormal * kV);
235

236
        // compute the eigenvalues of S (min and max curvatures)
237
        float fTrace = kS(0, 0) + kS(1, 1);
238
        float fDet = kS(0, 0) * kS(1, 1) - kS(0, 1) * kS(1, 0);
239
        float fDiscr = fTrace * fTrace - (4.0) * fDet;
240
        float fRootDiscr = sqrt(fabs(fDiscr));
241
        minCurvature = (0.5) * (fTrace - fRootDiscr);
242
        maxCurvature = (0.5) * (fTrace + fRootDiscr);
243

244
        // compute the eigenvectors of S
245
        Eigen::Vector2f kW0(kS(0, 1), minCurvature - kS(0, 0));
246
        Eigen::Vector2f kW1(minCurvature - kS(1, 1), kS(1, 0));
247
        if (kW0.squaredNorm() >= kW1.squaredNorm()) {
248
            float len = kW0.squaredNorm();
249
            if (len > 0 && len != 1) {
250
                kW0.normalize();
251
            }
252
            Eigen::Vector3f v = kU * kW0[0] + kV * kW0[1];
253
            minDirection.Set(v[0], v[1], v[2]);
254
        }
255
        else {
256
            float len = kW1.squaredNorm();
257
            if (len > 0 && len != 1) {
258
                kW1.normalize();
259
            }
260
            Eigen::Vector3f v = kU * kW1[0] + kV * kW1[1];
261
            minDirection.Set(v[0], v[1], v[2]);
262
        }
263

264
        kW0 = Eigen::Vector2f(kS(0, 1), maxCurvature - kS(0, 0));
265
        kW1 = Eigen::Vector2f(maxCurvature - kS(1, 1), kS(1, 0));
266
        if (kW0.squaredNorm() >= kW1.squaredNorm()) {
267
            float len = kW0.squaredNorm();
268
            if (len > 0 && len != 1) {
269
                kW0.normalize();
270
            }
271
            Eigen::Vector3f v = kU * kW0[0] + kV * kW0[1];
272
            maxDirection.Set(v[0], v[1], v[2]);
273
        }
274
        else {
275
            float len = kW1.squaredNorm();
276
            if (len > 0 && len != 1) {
277
                kW1.normalize();
278
            }
279
            Eigen::Vector3f v = kU * kW1[0] + kV * kW1[1];
280
            maxDirection.Set(v[0], v[1], v[2]);
281
        }
282

283
        CurvatureInfo ci;
284
        ci.fMaxCurvature = maxCurvature;
285
        ci.cMaxCurvDir = maxDirection;
286
        ci.fMinCurvature = minCurvature;
287
        ci.cMinCurvDir = minDirection;
288
        myCurvature.push_back(ci);
289
    }
290
}
291
#else
292
void MeshCurvature::ComputePerVertex()
293
{
294
    myCurvature.clear();
295

296
    // get all points
297
    std::vector<Wm4::Vector3<double>> aPnts;
298
    aPnts.reserve(myKernel.CountPoints());
299
    MeshPointIterator cPIt(myKernel);
300
    for (cPIt.Init(); cPIt.More(); cPIt.Next()) {
301
        Wm4::Vector3<double> cP(cPIt->x, cPIt->y, cPIt->z);
302
        aPnts.push_back(cP);
303
    }
304

305
    // get all point connections
306
    std::vector<int> aIdx;
307
    aIdx.reserve(3 * myKernel.CountFacets());
308
    const MeshFacetArray& raFts = myKernel.GetFacets();
309
    for (const auto& it : raFts) {
310
        for (PointIndex point : it._aulPoints) {
311
            aIdx.push_back(int(point));
312
        }
313
    }
314

315
    // in case of an empty mesh no curvature can be calculated
316
    if (myKernel.CountPoints() == 0 || myKernel.CountFacets() == 0) {
317
        return;
318
    }
319

320
    // compute vertex based curvatures
321
    Wm4::MeshCurvature<double> meshCurv(myKernel.CountPoints(),
322
                                        &(aPnts[0]),
323
                                        myKernel.CountFacets(),
324
                                        &(aIdx[0]));
325

326
    // get curvature information now
327
    const Wm4::Vector3<double>* aMaxCurvDir = meshCurv.GetMaxDirections();
328
    const Wm4::Vector3<double>* aMinCurvDir = meshCurv.GetMinDirections();
329
    const double* aMaxCurv = meshCurv.GetMaxCurvatures();
330
    const double* aMinCurv = meshCurv.GetMinCurvatures();
331

332
    myCurvature.reserve(myKernel.CountPoints());
333
    for (unsigned long i = 0; i < myKernel.CountPoints(); i++) {
334
        CurvatureInfo ci;
335
        ci.cMaxCurvDir = Base::Vector3f((float)aMaxCurvDir[i].X(),
336
                                        (float)aMaxCurvDir[i].Y(),
337
                                        (float)aMaxCurvDir[i].Z());
338
        ci.cMinCurvDir = Base::Vector3f((float)aMinCurvDir[i].X(),
339
                                        (float)aMinCurvDir[i].Y(),
340
                                        (float)aMinCurvDir[i].Z());
341
        ci.fMaxCurvature = (float)aMaxCurv[i];
342
        ci.fMinCurvature = (float)aMinCurv[i];
343
        myCurvature.push_back(ci);
344
    }
345
}
346
#endif  // OPTIMIZE_CURVATURE
347

348
// --------------------------------------------------------
349

350
namespace MeshCore
351
{
352
class FitPointCollector: public MeshCollector
353
{
354
public:
355
    explicit FitPointCollector(std::set<PointIndex>& ind)
356
        : indices(ind)
357
    {}
358
    void Append(const MeshCore::MeshKernel& kernel, FacetIndex index) override
359
    {
360
        PointIndex ulP1 {}, ulP2 {}, ulP3 {};
361
        kernel.GetFacetPoints(index, ulP1, ulP2, ulP3);
362
        indices.insert(ulP1);
363
        indices.insert(ulP2);
364
        indices.insert(ulP3);
365
    }
366

367
private:
368
    std::set<PointIndex>& indices;
369
};
370
}  // namespace MeshCore
371

372
// --------------------------------------------------------
373

374
FacetCurvature::FacetCurvature(const MeshKernel& kernel,
375
                               const MeshRefPointToFacets& search,
376
                               float r,
377
                               unsigned long pt)
378
    : myKernel(kernel)
379
    , mySearch(search)
380
    , myMinPoints(pt)
381
    , myRadius(r)
382
{}
383

384
CurvatureInfo FacetCurvature::Compute(FacetIndex index) const
385
{
386
    Base::Vector3f rkDir0, rkDir1;
387
    Base::Vector3f rkNormal;
388

389
    MeshGeomFacet face = myKernel.GetFacet(index);
390
    Base::Vector3f face_gravity = face.GetGravityPoint();
391
    Base::Vector3f face_normal = face.GetNormal();
392
    std::set<PointIndex> point_indices;
393
    FitPointCollector collect(point_indices);
394

395
    float searchDist = myRadius;
396
    int attempts = 0;
397
    do {
398
        mySearch.Neighbours(index, searchDist, collect);
399
        if (point_indices.empty()) {
400
            break;
401
        }
402
        float min_points = myMinPoints;
403
        float use_points = point_indices.size();
404
        searchDist = searchDist * sqrt(min_points / use_points);
405
    } while ((point_indices.size() < myMinPoints) && (attempts++ < 3));
406

407
    std::vector<Base::Vector3f> fitPoints;
408
    const MeshPointArray& verts = myKernel.GetPoints();
409
    fitPoints.reserve(point_indices.size());
410
    for (PointIndex it : point_indices) {
411
        fitPoints.push_back(verts[it] - face_gravity);
412
    }
413

414
    float fMin {}, fMax {};
415
    if (fitPoints.size() >= myMinPoints) {
416
        SurfaceFit surf_fit;
417
        surf_fit.AddPoints(fitPoints);
418
        surf_fit.Fit();
419
        rkNormal = surf_fit.GetNormal();
420
        double dMin {}, dMax {}, dDistance {};
421
        if (surf_fit.GetCurvatureInfo(0.0, 0.0, 0.0, dMin, dMax, rkDir1, rkDir0, dDistance)) {
422
            fMin = (float)dMin;
423
            fMax = (float)dMax;
424
        }
425
        else {
426
            fMin = FLT_MAX;
427
            fMax = FLT_MAX;
428
        }
429
    }
430
    else {
431
        // too few points => cannot calc any properties
432
        fMin = FLT_MAX;
433
        fMax = FLT_MAX;
434
    }
435

436
    CurvatureInfo info;
437
    if (fMin < fMax) {
438
        info.fMaxCurvature = fMax;
439
        info.fMinCurvature = fMin;
440
        info.cMaxCurvDir = rkDir1;
441
        info.cMinCurvDir = rkDir0;
442
    }
443
    else {
444
        info.fMaxCurvature = fMin;
445
        info.fMinCurvature = fMax;
446
        info.cMaxCurvDir = rkDir0;
447
        info.cMinCurvDir = rkDir1;
448
    }
449

450
    // Reverse the direction of the normal vector if required
451
    // (Z component of "local" normal vectors should be opposite in sign to the "local" view vector)
452
    if (rkNormal * face_normal < 0.0) {
453
        // Note: Changing the normal directions is similar to flipping over the object.
454
        // In this case we must adjust the curvature information as well.
455
        std::swap(info.cMaxCurvDir, info.cMinCurvDir);
456
        std::swap(info.fMaxCurvature, info.fMinCurvature);
457
        info.fMaxCurvature *= (-1.0);
458
        info.fMinCurvature *= (-1.0);
459
    }
460

461
    return info;
462
}
463