FreeCAD

1
/***************************************************************************
2
 *   Copyright (c) 2009 Werner Mayer <wmayer[at]users.sourceforge.net>     *
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

25
#include <Base/Tools.h>
26

27
#include "Algorithm.h"
28
#include "Approximation.h"
29
#include "Iterator.h"
30
#include "MeshKernel.h"
31
#include "Smoothing.h"
32

33

34
using namespace MeshCore;
35

36

37
AbstractSmoothing::AbstractSmoothing(MeshKernel& m)
38
    : kernel(m)
39
{}
40

41
AbstractSmoothing::~AbstractSmoothing() = default;
42

43
void AbstractSmoothing::initialize(Component comp, Continuity cont)
44
{
45
    this->component = comp;
46
    this->continuity = cont;
47
}
48

49
PlaneFitSmoothing::PlaneFitSmoothing(MeshKernel& m)
50
    : AbstractSmoothing(m)
51
{}
52

53
void PlaneFitSmoothing::Smooth(unsigned int iterations)
54
{
55
    MeshCore::MeshPoint center;
56
    MeshCore::MeshPointArray PointArray = kernel.GetPoints();
57

58
    MeshCore::MeshPointIterator v_it(kernel);
59
    MeshCore::MeshRefPointToPoints vv_it(kernel);
60
    MeshCore::MeshPointArray::_TConstIterator v_beg = kernel.GetPoints().begin();
61

62
    for (unsigned int i = 0; i < iterations; i++) {
63
        Base::Vector3f N, L;
64
        for (v_it.Begin(); v_it.More(); v_it.Next()) {
65
            MeshCore::PlaneFit pf;
66
            pf.AddPoint(*v_it);
67
            center = *v_it;
68
            const std::set<PointIndex>& cv = vv_it[v_it.Position()];
69
            if (cv.size() < 3) {
70
                continue;
71
            }
72

73
            std::set<PointIndex>::const_iterator cv_it;
74
            for (cv_it = cv.begin(); cv_it != cv.end(); ++cv_it) {
75
                pf.AddPoint(v_beg[*cv_it]);
76
                center += v_beg[*cv_it];
77
            }
78

79
            float scale = 1.0f / (static_cast<float>(cv.size()) + 1.0f);
80
            center.Scale(scale, scale, scale);
81

82
            // get the mean plane of the current vertex with the surrounding vertices
83
            pf.Fit();
84
            N = pf.GetNormal();
85
            N.Normalize();
86

87
            // look in which direction we should move the vertex
88
            L.Set(v_it->x - center.x, v_it->y - center.y, v_it->z - center.z);
89
            if (N * L < 0.0f) {
90
                N.Scale(-1.0, -1.0, -1.0);
91
            }
92

93
            // maximum value to move is distance to mean plane
94
            float d = std::min<float>(fabs(this->maximum), fabs(N * L));
95
            N.Scale(d, d, d);
96

97
            PointArray[v_it.Position()].Set(v_it->x - N.x, v_it->y - N.y, v_it->z - N.z);
98
        }
99

100
        // assign values without affecting iterators
101
        PointIndex count = kernel.CountPoints();
102
        for (PointIndex idx = 0; idx < count; idx++) {
103
            kernel.SetPoint(idx, PointArray[idx]);
104
        }
105
    }
106
}
107

108
void PlaneFitSmoothing::SmoothPoints(unsigned int iterations,
109
                                     const std::vector<PointIndex>& point_indices)
110
{
111
    MeshCore::MeshPoint center;
112
    MeshCore::MeshPointArray PointArray = kernel.GetPoints();
113

114
    MeshCore::MeshPointIterator v_it(kernel);
115
    MeshCore::MeshRefPointToPoints vv_it(kernel);
116
    MeshCore::MeshPointArray::_TConstIterator v_beg = kernel.GetPoints().begin();
117

118
    for (unsigned int i = 0; i < iterations; i++) {
119
        Base::Vector3f N, L;
120
        for (PointIndex it : point_indices) {
121
            v_it.Set(it);
122
            MeshCore::PlaneFit pf;
123
            pf.AddPoint(*v_it);
124
            center = *v_it;
125
            const std::set<PointIndex>& cv = vv_it[v_it.Position()];
126
            if (cv.size() < 3) {
127
                continue;
128
            }
129

130
            std::set<PointIndex>::const_iterator cv_it;
131
            for (cv_it = cv.begin(); cv_it != cv.end(); ++cv_it) {
132
                pf.AddPoint(v_beg[*cv_it]);
133
                center += v_beg[*cv_it];
134
            }
135

136
            float scale = 1.0f / (static_cast<float>(cv.size()) + 1.0f);
137
            center.Scale(scale, scale, scale);
138

139
            // get the mean plane of the current vertex with the surrounding vertices
140
            pf.Fit();
141
            N = pf.GetNormal();
142
            N.Normalize();
143

144
            // look in which direction we should move the vertex
145
            L.Set(v_it->x - center.x, v_it->y - center.y, v_it->z - center.z);
146
            if (N * L < 0.0f) {
147
                N.Scale(-1.0, -1.0, -1.0);
148
            }
149

150
            // maximum value to move is distance to mean plane
151
            float d = std::min<float>(fabs(this->maximum), fabs(N * L));
152
            N.Scale(d, d, d);
153

154
            PointArray[v_it.Position()].Set(v_it->x - N.x, v_it->y - N.y, v_it->z - N.z);
155
        }
156

157
        // assign values without affecting iterators
158
        PointIndex count = kernel.CountPoints();
159
        for (PointIndex idx = 0; idx < count; idx++) {
160
            kernel.SetPoint(idx, PointArray[idx]);
161
        }
162
    }
163
}
164

165
LaplaceSmoothing::LaplaceSmoothing(MeshKernel& m)
166
    : AbstractSmoothing(m)
167
{}
168

169
void LaplaceSmoothing::Umbrella(const MeshRefPointToPoints& vv_it,
170
                                const MeshRefPointToFacets& vf_it,
171
                                double stepsize)
172
{
173
    const MeshCore::MeshPointArray& points = kernel.GetPoints();
174
    MeshCore::MeshPointArray::_TConstIterator v_it, v_beg = points.begin(), v_end = points.end();
175

176
    PointIndex pos = 0;
177
    for (v_it = points.begin(); v_it != v_end; ++v_it, ++pos) {
178
        const std::set<PointIndex>& cv = vv_it[pos];
179
        if (cv.size() < 3) {
180
            continue;
181
        }
182
        if (cv.size() != vf_it[pos].size()) {
183
            // do nothing for border points
184
            continue;
185
        }
186

187
        size_t n_count = cv.size();
188
        double w {};
189
        w = 1.0 / double(n_count);
190

191
        double delx = 0.0, dely = 0.0, delz = 0.0;
192
        std::set<PointIndex>::const_iterator cv_it;
193
        for (cv_it = cv.begin(); cv_it != cv.end(); ++cv_it) {
194
            delx += w * static_cast<double>((v_beg[*cv_it]).x - v_it->x);
195
            dely += w * static_cast<double>((v_beg[*cv_it]).y - v_it->y);
196
            delz += w * static_cast<double>((v_beg[*cv_it]).z - v_it->z);
197
        }
198

199
        float x = static_cast<float>(static_cast<double>(v_it->x) + stepsize * delx);
200
        float y = static_cast<float>(static_cast<double>(v_it->y) + stepsize * dely);
201
        float z = static_cast<float>(static_cast<double>(v_it->z) + stepsize * delz);
202
        kernel.SetPoint(pos, x, y, z);
203
    }
204
}
205

206
void LaplaceSmoothing::Umbrella(const MeshRefPointToPoints& vv_it,
207
                                const MeshRefPointToFacets& vf_it,
208
                                double stepsize,
209
                                const std::vector<PointIndex>& point_indices)
210
{
211
    const MeshCore::MeshPointArray& points = kernel.GetPoints();
212
    MeshCore::MeshPointArray::_TConstIterator v_beg = points.begin();
213

214
    for (PointIndex it : point_indices) {
215
        const std::set<PointIndex>& cv = vv_it[it];
216
        if (cv.size() < 3) {
217
            continue;
218
        }
219
        if (cv.size() != vf_it[it].size()) {
220
            // do nothing for border points
221
            continue;
222
        }
223

224
        size_t n_count = cv.size();
225
        double w {};
226
        w = 1.0 / double(n_count);
227

228
        double delx = 0.0, dely = 0.0, delz = 0.0;
229
        std::set<PointIndex>::const_iterator cv_it;
230
        for (cv_it = cv.begin(); cv_it != cv.end(); ++cv_it) {
231
            delx += w * static_cast<double>((v_beg[*cv_it]).x - (v_beg[it]).x);
232
            dely += w * static_cast<double>((v_beg[*cv_it]).y - (v_beg[it]).y);
233
            delz += w * static_cast<double>((v_beg[*cv_it]).z - (v_beg[it]).z);
234
        }
235

236
        float x = static_cast<float>(static_cast<double>((v_beg[it]).x) + stepsize * delx);
237
        float y = static_cast<float>(static_cast<double>((v_beg[it]).y) + stepsize * dely);
238
        float z = static_cast<float>(static_cast<double>((v_beg[it]).z) + stepsize * delz);
239
        kernel.SetPoint(it, x, y, z);
240
    }
241
}
242

243
void LaplaceSmoothing::Smooth(unsigned int iterations)
244
{
245
    MeshCore::MeshRefPointToPoints vv_it(kernel);
246
    MeshCore::MeshRefPointToFacets vf_it(kernel);
247

248
    for (unsigned int i = 0; i < iterations; i++) {
249
        Umbrella(vv_it, vf_it, lambda);
250
    }
251
}
252

253
void LaplaceSmoothing::SmoothPoints(unsigned int iterations,
254
                                    const std::vector<PointIndex>& point_indices)
255
{
256
    MeshCore::MeshRefPointToPoints vv_it(kernel);
257
    MeshCore::MeshRefPointToFacets vf_it(kernel);
258

259
    for (unsigned int i = 0; i < iterations; i++) {
260
        Umbrella(vv_it, vf_it, lambda, point_indices);
261
    }
262
}
263

264
TaubinSmoothing::TaubinSmoothing(MeshKernel& m)
265
    : LaplaceSmoothing(m)
266
{}
267

268
void TaubinSmoothing::Smooth(unsigned int iterations)
269
{
270
    MeshCore::MeshRefPointToPoints vv_it(kernel);
271
    MeshCore::MeshRefPointToFacets vf_it(kernel);
272

273
    // Theoretically Taubin does not shrink the surface
274
    iterations = (iterations + 1) / 2;  // two steps per iteration
275
    for (unsigned int i = 0; i < iterations; i++) {
276
        Umbrella(vv_it, vf_it, GetLambda());
277
        Umbrella(vv_it, vf_it, -(GetLambda() + micro));
278
    }
279
}
280

281
void TaubinSmoothing::SmoothPoints(unsigned int iterations,
282
                                   const std::vector<PointIndex>& point_indices)
283
{
284
    MeshCore::MeshRefPointToPoints vv_it(kernel);
285
    MeshCore::MeshRefPointToFacets vf_it(kernel);
286

287
    // Theoretically Taubin does not shrink the surface
288
    iterations = (iterations + 1) / 2;  // two steps per iteration
289
    for (unsigned int i = 0; i < iterations; i++) {
290
        Umbrella(vv_it, vf_it, GetLambda(), point_indices);
291
        Umbrella(vv_it, vf_it, -(GetLambda() + micro), point_indices);
292
    }
293
}
294

295
namespace
296
{
297
using AngleNormal = std::pair<double, Base::Vector3d>;
298
inline Base::Vector3d find_median(std::vector<AngleNormal>& container)
299
{
300
    auto compare_angle_normal = [](const AngleNormal& an1, const AngleNormal& an2) {
301
        return an1.first < an2.first;
302
    };
303
    size_t n = container.size() / 2;
304
    std::nth_element(container.begin(),
305
                     container.begin() + n,
306
                     container.end(),
307
                     compare_angle_normal);
308

309
    if ((container.size() % 2) == 1) {
310
        return container[n].second;
311
    }
312
    else {
313
        // even sized vector -> average the two middle values
314
        auto max_it =
315
            std::max_element(container.begin(), container.begin() + n, compare_angle_normal);
316
        Base::Vector3d vec = (max_it->second + container[n].second) / 2.0;
317
        vec.Normalize();
318
        return vec;
319
    }
320
}
321
}  // namespace
322

323
MedianFilterSmoothing::MedianFilterSmoothing(MeshKernel& m)
324
    : AbstractSmoothing(m)
325
{}
326

327
void MedianFilterSmoothing::Smooth(unsigned int iterations)
328
{
329
    std::vector<unsigned long> point_indices(kernel.CountPoints());
330
    std::generate(point_indices.begin(), point_indices.end(), Base::iotaGen<unsigned long>(0));
331
    MeshCore::MeshRefFacetToFacets ff_it(kernel);
332
    MeshCore::MeshRefPointToFacets vf_it(kernel);
333

334
    for (unsigned int i = 0; i < iterations; i++) {
335
        UpdatePoints(ff_it, vf_it, point_indices);
336
    }
337
}
338

339
void MedianFilterSmoothing::SmoothPoints(unsigned int iterations,
340
                                         const std::vector<PointIndex>& point_indices)
341
{
342
    MeshCore::MeshRefFacetToFacets ff_it(kernel);
343
    MeshCore::MeshRefPointToFacets vf_it(kernel);
344

345
    for (unsigned int i = 0; i < iterations; i++) {
346
        UpdatePoints(ff_it, vf_it, point_indices);
347
    }
348
}
349

350
void MedianFilterSmoothing::UpdatePoints(const MeshRefFacetToFacets& ff_it,
351
                                         const MeshRefPointToFacets& vf_it,
352
                                         const std::vector<PointIndex>& point_indices)
353
{
354
    const MeshCore::MeshPointArray& points = kernel.GetPoints();
355
    const MeshCore::MeshFacetArray& facets = kernel.GetFacets();
356

357
    // Initialize the array with the real normals
358
    std::vector<Base::Vector3d> faceNormals;
359
    faceNormals.reserve(facets.size());
360
    MeshCore::MeshFacetIterator iter(kernel);
361
    for (iter.Init(); iter.More(); iter.Next()) {
362
        faceNormals.emplace_back(Base::toVector<double>(iter->GetNormal()));
363
    }
364

365
    // Step 1: determine face normals
366
    for (FacetIndex pos = 0; pos < facets.size(); pos++) {
367
        iter.Set(pos);
368
        Base::Vector3d refNormal = Base::toVector<double>(iter->GetNormal());
369
        const std::set<FacetIndex>& cv = ff_it[pos];
370
        const MeshCore::MeshFacet& facet = facets[pos];
371

372
        std::vector<AngleNormal> anglesWithFaces;
373
        for (auto fi : cv) {
374
            iter.Set(fi);
375
            Base::Vector3d faceNormal = Base::toVector<double>(iter->GetNormal());
376
            double angle = refNormal.GetAngle(faceNormal);
377

378
            int absWeight = std::abs(weights);
379
            if (absWeight > 1 && facet.IsNeighbour(fi)) {
380
                if (weights < 0) {
381
                    angle = -angle;
382
                }
383
                for (int i = 0; i < absWeight; i++) {
384
                    anglesWithFaces.emplace_back(angle, faceNormal);
385
                }
386
            }
387
            else {
388
                anglesWithFaces.emplace_back(angle, faceNormal);
389
            }
390
        }
391

392
        faceNormals[pos] = find_median(anglesWithFaces);
393
    }
394

395
    // Step 2: move vertices
396
    for (auto pos : point_indices) {
397
        Base::Vector3d P = Base::toVector<double>(points[pos]);
398
        const std::set<FacetIndex>& cv = vf_it[pos];
399

400
        double totalArea = 0.0;
401
        Base::Vector3d totalvT;
402
        for (auto it : cv) {
403
            iter.Set(it);
404

405
            double faceArea = iter->Area();
406
            totalArea += faceArea;
407

408
            Base::Vector3d C = Base::toVector<double>(iter->GetGravityPoint());
409

410
            Base::Vector3d PC = C - P;
411
            Base::Vector3d mT = faceNormals[it];
412
            Base::Vector3d vT = (PC * mT) * mT;
413
            totalvT += vT * faceArea;
414
        }
415

416
        P = P + totalvT / totalArea;
417
        kernel.SetPoint(pos, Base::toVector<float>(P));
418
    }
419
}
420