1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include <stdio.h>
7 #include "SmoothData.h"
8 #include "Numeric.h"
9 #include "OS.h"
10
11 // Basic coordinate-based floting point value averager
12
13 double xyzv::eps = 1.e-12;
14
xyzv(const xyzv & other)15 xyzv::xyzv(const xyzv &other)
16 {
17 x = other.x;
18 y = other.y;
19 z = other.z;
20 scaleValue = other.scaleValue; // Added by Trevor Strickler 07/10/2013
21 scale_numvals = other.scale_numvals; // Added by Trevor Strickler 07/10/2013
22 nbvals = other.nbvals;
23 nboccurrences = other.nboccurrences;
24 if(other.vals && other.nbvals) {
25 vals = new double[other.nbvals];
26 for(int i = 0; i < nbvals; i++) vals[i] = other.vals[i];
27 }
28 else
29 vals = nullptr;
30 }
31
operator =(const xyzv & other)32 xyzv &xyzv::operator=(const xyzv &other)
33 {
34 if(this != &other) {
35 x = other.x;
36 y = other.y;
37 z = other.z;
38 scaleValue = other.scaleValue; // Added by Trevor Strickler 07/10/2013
39 scale_numvals = other.scale_numvals; // Added by Trevor Strickler 07/10/2013
40 nbvals = other.nbvals;
41 nboccurrences = other.nboccurrences;
42 if(other.vals && other.nbvals) {
43 vals = new double[other.nbvals];
44 for(int i = 0; i < nbvals; i++) vals[i] = other.vals[i];
45 }
46 }
47 return *this;
48 }
49
update(int n,double * v)50 void xyzv::update(int n, double *v)
51 {
52 if(!vals) {
53 vals = new double[n];
54 for(int i = 0; i < n; i++) vals[i] = 0.0;
55 nbvals = n;
56 nboccurrences = 0;
57 }
58 else if(nbvals != n)
59 return; // error
60 double x1 = (double)(nboccurrences) / (double)(nboccurrences + 1);
61 double x2 = 1. / (double)(nboccurrences + 1);
62 for(int i = 0; i < nbvals; i++) vals[i] = (x1 * vals[i] + x2 * v[i]);
63 nboccurrences++;
64 }
65
66 // Added by Trevor Strickler
scale_update(double scale_inp)67 void xyzv::scale_update(double scale_inp)
68 {
69 if(std::abs(1.0 - scale_inp) <= eps) scale_inp = 1.0;
70 if(scale_inp != 1.0 || scaleValue != 1.0) {
71 double x1 = (double)(scale_numvals) / (double)(scale_numvals + 1);
72 double x2 = 1.0 / (double)(scale_numvals + 1);
73 scaleValue = (x1 * scaleValue + x2 * scale_inp);
74 }
75 if(std::abs(1.0 - scaleValue) <= eps) scaleValue = 1.0;
76 scale_numvals++;
77 }
78
add(double x,double y,double z,int n,double * vals)79 void smooth_data::add(double x, double y, double z, int n, double *vals)
80 {
81 xyzv xyz(x, y, z);
82 auto it = c.find(xyz);
83 if(it == c.end()) {
84 xyz.update(n, vals);
85 c.insert(xyz);
86 }
87 else {
88 // we can do this because we know that it will not destroy the set
89 // ordering
90 xyzv *p = (xyzv *)&(*it);
91 p->update(n, vals);
92 }
93 }
94
95 // added by Trevor Strickler
add_scale(double x,double y,double z,double scale_val)96 void smooth_data::add_scale(double x, double y, double z, double scale_val)
97 {
98 xyzv xyz(x, y, z);
99 auto it = c.find(xyz);
100 if(it == c.end()) {
101 xyz.scale_update(scale_val);
102 c.insert(xyz);
103 }
104 else {
105 // we can do this because we know that it will not destroy the set
106 // ordering
107 xyzv *p = (xyzv *)&(*it);
108 p->scale_update(scale_val);
109 }
110 }
111
get(double x,double y,double z,int n,double * vals) const112 bool smooth_data::get(double x, double y, double z, int n, double *vals) const
113 {
114 auto it = c.find(xyzv(x, y, z));
115 if(it == c.end()) return false;
116 for(int k = 0; k < n; k++) vals[k] = it->vals[k];
117 return true;
118 }
119
120 // added by Trevor Strickler
get_scale(double x,double y,double z,double * scale_val) const121 bool smooth_data::get_scale(double x, double y, double z, double *scale_val) const
122 {
123 auto it = c.find(xyzv(x, y, z));
124 if(it == c.end()) return false;
125 (*scale_val) = it->scaleValue;
126 return true;
127 }
128
normalize()129 void smooth_data::normalize()
130 {
131 auto it = c.begin();
132 while(it != c.end()) {
133 if(it->nbvals == 3) norme(it->vals);
134 it++;
135 }
136 }
137
exportview(const std::string & filename) const138 bool smooth_data::exportview(const std::string &filename) const
139 {
140 FILE *fp = Fopen(filename.c_str(), "w");
141 if(!fp) return false;
142 fprintf(fp, "View \"data\" {\n");
143 auto it = c.begin();
144 while(it != c.end()) {
145 switch(it->nbvals) {
146 case 1:
147 fprintf(fp, "SP(%.16g,%.16g,%.16g){%.16g};\n", it->x, it->y, it->z,
148 it->vals[0]);
149 break;
150 case 3:
151 fprintf(fp, "VP(%.16g,%.16g,%.16g){%.16g,%.16g,%.16g};\n", it->x, it->y,
152 it->z, it->vals[0], it->vals[1], it->vals[2]);
153 break;
154 }
155 it++;
156 }
157 fprintf(fp, "};\n");
158 fclose(fp);
159 return true;
160 }
161
162 // Normal smoother
163
164 float xyzn::eps = 1.e-6F;
165
angle(int i,char nx,char ny,char nz)166 float xyzn::angle(int i, char nx, char ny, char nz)
167 {
168 // returns the angle (in [-180,180]) between the ith normal stored
169 // at point xyz and the new normal nx,ny,nz
170 double a[3] = {char2float(n[i].nx), char2float(n[i].ny), char2float(n[i].nz)};
171 double b[3] = {char2float(nx), char2float(ny), char2float(nz)};
172 norme(a);
173 norme(b);
174
175 double c[3];
176 prodve(a, b, c);
177
178 double const cosc = prosca(a, b);
179 double const sinc = std::sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
180 double const angplan = myatan2(sinc, cosc);
181 return (float)(angplan * 180. / M_PI);
182 }
183
update(char nx,char ny,char nz,float tol)184 void xyzn::update(char nx, char ny, char nz, float tol)
185 {
186 // just ignore it if we have more than 100 clusters
187 if(n.size() > 100) return;
188
189 // we average by clusters of normals separated by tol; the result of
190 // the averaging depends on the order in which we average (since we
191 // store the average value as the cluster center as we go), but it
192 // seems to work very nicely in practice (and it's faster than
193 // storing everyting and averaging at the end)
194 for(std::size_t i = 0; i < n.size(); i++) {
195 if(tol >= 180. || std::abs(angle(i, nx, ny, nz)) < tol) {
196 // just ignore it if we have more than 100 contributions to a
197 // single point...
198 if(n[i].nb < 100) {
199 float c1 = (float)(n[i].nb) / (float)(n[i].nb + 1);
200 float c2 = 1.0F / (float)(n[i].nb + 1);
201 n[i].nx = (char)(c1 * n[i].nx + c2 * nx);
202 n[i].ny = (char)(c1 * n[i].ny + c2 * ny);
203 n[i].nz = (char)(c1 * n[i].nz + c2 * nz);
204 n[i].nb++;
205 }
206 return;
207 }
208 }
209
210 // create a new cluster
211 nnb nn = {nx, ny, nz, 0};
212 n.push_back(nn);
213 }
214
add(double x,double y,double z,double nx,double ny,double nz)215 void smooth_normals::add(double x, double y, double z, double nx, double ny,
216 double nz)
217 {
218 xyzn xyz((float)x, (float)y, (float)z);
219 auto it = c.find(xyz);
220 if(it == c.end()) {
221 xyz.update(float2char((float)nx), float2char((float)ny),
222 float2char((float)nz), tol);
223 c.insert(xyz);
224 }
225 else {
226 // we can do this because we know that it will not destroy the set
227 // ordering
228 xyzn *p = (xyzn *)&(*it);
229 p->update(float2char((float)nx), float2char((float)ny),
230 float2char((float)nz), tol);
231 }
232 }
233
get(double x,double y,double z,double & nx,double & ny,double & nz) const234 bool smooth_normals::get(double x, double y, double z, double &nx, double &ny,
235 double &nz) const
236 {
237 auto it =
238 c.find(xyzn((float)x, (float)y, (float)z));
239
240 if(it == c.end()) return false;
241
242 xyzn *p = (xyzn *)&(*it);
243 for(std::size_t i = 0; i < p->n.size(); i++) {
244 if(std::abs(p->angle(i, float2char((float)nx), float2char((float)ny),
245 float2char((float)nz))) < tol) {
246 nx = char2float(p->n[i].nx);
247 ny = char2float(p->n[i].ny);
248 nz = char2float(p->n[i].nz);
249 break;
250 }
251 }
252 return true;
253 }
254