1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Img3D/Mesh/Sphere.cpp
//! @brief Implements utility functions in ba3d namespace.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Base/Util/Assert.h"
#include "Img3D/Model/Geometry.h"
#include <numbers>
using std::numbers::pi;
namespace Img3D {
// cut: 0..1 - how much is cut off from the bottom
// cutFromTop - how much fraction of the radius is removed from the top
Geometry::Mesh Geometry::meshSphere(float cut, float baseShift, float cutFromTop)
{
if (1 <= cut)
return {};
cut = qMax(0.f, cut);
ASSERT(0 <= cut && cut < 1);
// 'rings' are the # of horizontal cross-sections ranging from bottom to top of the sphere
// 'slices' are the # of vertices in a given ring
int rings;
int slices = SLICES;
float minPh, maxPh, phRge;
if (cut > 0) // South pole absent
minPh = asinf(2 * cut - 1);
else // South pole present
minPh = -float((pi / 2)) + float(pi) / RINGS;
if (cutFromTop > 0) // North pole absent
maxPh = asinf(1 - 2 * cutFromTop);
else // North pole present
maxPh = float((pi / 2)) - float(pi) / RINGS;
phRge = maxPh - minPh;
rings = qMax(2, qCeil(qreal(RINGS * phRge) / pi)); // At least 2 rings (incl. lowest ring)
ASSERT(qAbs(minPh) < float((pi / 2)));
ASSERT(2 <= rings && 2 <= slices);
// meshes of vertices and normals, without poles, _[ring][slice]
QVector<Vertices> vs_(rings);
QVector<Vertices> ns_(rings);
for (auto& ring : vs_)
ring.resize(slices);
for (auto& ring : ns_)
ring.resize(slices);
float const R = .5f;
for (int r = 0; r < rings; ++r) {
float ph = minPh + phRge * r / (rings - 1);
float cp = cosf(ph), sp = sinf(ph);
for (int s = 0; s < slices; ++s) {
auto th = float((2 * pi) * s / slices);
F3 v(R * cp * cosf(th), R * cp * sinf(th), R * sp + baseShift);
// baseShift is used for shifting the bottom of the spherical shape to z=0 plane
vs_[r][s] = v;
ns_[r][s] = v.normalized();
}
}
// make into triangles
int const nv = 6 * (rings)*slices;
Vertices vs;
vs.reserve(nv);
Vertices ns;
ns.reserve(nv);
for (int r = 0; r < rings; ++r) {
auto &vr = vs_.at(r), &nr = ns_.at(r);
for (int s = 0; s < slices; ++s) {
int s0 = s, s1 = (s + 1) % slices;
auto &v0 = vr.at(s0), &v1 = vr.at(s1);
auto &n0 = nr.at(s0), &n1 = nr.at(s1);
if (r == 0) { // bottom most ring
F3 vp, n0, n1, np(F3(0, 0, -1));
if (cut > 0) { // South pole absent
vp = F3(0, 0, v0.z());
n0 = n1 = np;
} else { // South pole present
vp = F3(0, 0, -R + baseShift);
n0 = nr.at(s0);
n1 = nr.at(s1);
}
vs.addTriangle(v0, vp, v1);
ns.addTriangle(n0, np, n1);
}
if (r + 1 == rings) { // top most ring
F3 vp, n0, n1;
F3 np(0, 0, 1);
if (cutFromTop > 0) { // North pole absent
vp = F3(0, 0, v0.z());
n0 = n1 = np;
} else { // North pole present
vp = F3(0, 0, +R + baseShift);
n0 = nr.at(s0);
n1 = nr.at(s1);
}
vs.addTriangle(v0, v1, vp);
ns.addTriangle(n0, n1, np);
} else { // in between poles
auto &vr1 = vs_.at(r + 1), &nr1 = ns_.at(r + 1);
auto &n2 = nr1.at(s1), &n3 = nr1.at(s0);
vs.addQuad(v0, v1, vr1.at(s1), vr1.at(s0));
ns.addQuad(n0, n1, n2, n3);
}
}
}
return makeMesh(vs, ns);
}
} // namespace Img3D
|
