?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
feature magic is an amplification factor for edge sharpening (not corners)
[engrid.git] / src / libengrid / createcadtesselation.cpp
blob59f2a87b551f7854845914e6f7b7daf8747b2a26
1 //
2 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 // + +
4 // + This file is part of enGrid. +
5 // + +
6 // + Copyright 2008-2013 enGits GmbH +
7 // + +
8 // + enGrid is free software: you can redistribute it and/or modify +
9 // + it under the terms of the GNU General Public License as published by +
10 // + the Free Software Foundation, either version 3 of the License, or +
11 // + (at your option) any later version. +
12 // + +
13 // + enGrid is distributed in the hope that it will be useful, +
14 // + but WITHOUT ANY WARRANTY; without even the implied warranty of +
15 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +
16 // + GNU General Public License for more details. +
17 // + +
18 // + You should have received a copy of the GNU General Public License +
19 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>. +
20 // + +
21 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22 //
23 #include "createcadtesselation.h"
24 #include "guimainwindow.h"
25
26 #include <vtkImageData.h>
27 #include <vtkXMLImageDataWriter.h>
28 #include <vtkFloatArray.h>
29 #include <vtkTriangleFilter.h>
30 #include <vtkContourFilter.h>
31 #include <vtkDecimatePro.h>
32
33 CreateCadTesselation::CreateCadTesselation(CadInterface *cad_interface)
34 {
35 m_CadInterface = cad_interface;
36 m_ScanMemory = 0.5;
37 m_NumIterations = 0;
38 m_PreservationType = 0;
39 m_SmallestFeatureSize = 1e99;
40 m_SmallestResolution = 0;
41 m_TargetReduction = 0;
42 }
43
44 bool CreateCadTesselation::shootRay(vec3_t x, vec3_t v, vec3_t &x_in, vec3_t &x_out, vec3_t &n_in, vec3_t &n_out)
45 {
46 v.normalise();
47 double r_hit;
48 vec3_t x_hit, n_hit;
49 CadInterface::HitType hit_type = m_CadInterface->shootRay(x, v, x_hit, n_hit, r_hit);
50 if (hit_type == CadInterface::Miss || hit_type == CadInterface::HitOut) {
51 return false;
52 }
53 x_in = x_hit;
54 n_in = n_hit;
55 x = x_in;
56 do {
57 x += 1e-10*v;
58 hit_type = m_CadInterface->shootRay(x, v, x_hit, n_hit, r_hit);
59 if (hit_type == CadInterface::HitOut) {
60 x_out = x_hit;
61 n_out = n_hit;
62 }
63 x = x_hit;
64 } while (hit_type == CadInterface::HitIn);
65 return true;
66 }
67
68 void CreateCadTesselation::scan(bool create_grid, int interlaces)
69 {
70 m_Dx = (m_X2[0] - m_X1[0])/(m_Ni-1);
71 m_Dy = (m_X2[1] - m_X1[1])/(m_Nj-1);
72 m_Dz = (m_X2[2] - m_X1[2])/(m_Nk-1);
73
74 EG_VTKSP(vtkImageData, gdata);
75 gdata->SetDimensions(m_Ni, m_Nj, m_Nk);
76 EG_VTKSP(vtkFloatArray, g);
77 g->SetName("g");
78 g->SetNumberOfValues(m_Ni*m_Nj*m_Nk);
79 gdata->GetPointData()->AddArray(g);
80 for (int i = 0; i < m_Ni; ++i) {
81 for (int j = 0; j < m_Nj; ++j) {
82 for (int k = 0; k < m_Nk; ++k) {
83 if (preserveFluid()) {
84 g->SetValue(getIdx(i,j,k), 1);
85 } else {
86 g->SetValue(getIdx(i,j,k), 0);
87 }
88 }
89 }
90 }
91 m_XScan1 = m_X2;
92 m_XScan2 = m_X1;
93
94 vec3_t x_in, x_out, n_in, n_out;
95
96 int progress = 0;
97 m_GeometryFound = false;
98
99 double dxi = m_Dx/(interlaces + 1);
100 double dyi = m_Dy/(interlaces + 1);
101 double dzi = m_Dz/(interlaces + 1);
102
103 double shift = 1e-6;
104
105 // scan ij plane -> k direction
106 // xy plane -> z direction
107 for (int i = 0; i < m_Ni; ++i) {
108 for (int j = 0; j < m_Nj; ++j) {
109 vec3_t x0 = getX(i,j,0);
110 for (int i_il = 0; i_il <= interlaces; ++i_il) {
111 for (int j_il = 0; j_il <= interlaces; ++j_il) {
112 vec3_t x = x0;
113 x[0] += i_il*dxi;
114 x[1] += j_il*dyi;
115 int k_last = 0;
116 while (shootRay(x, vec3_t(0,0,1), x_in, x_out, n_in, n_out)) {
117 m_GeometryFound = true;
118 m_XScan1[2] = min(m_XScan1[2], x_in[2]);
119 m_XScan2[2] = max(m_XScan2[2], x_out[2]);
120 int k1 = int((x_in[2] - m_X1[2])/m_Dz) + 1;
121 int k2 = int((x_out[2] - m_X1[2])/m_Dz);
122 if (preserveFluid()) {
123 for (int k = k_last; k < min(m_Nk-1,k1); ++k) {
124 g->SetValue(getIdx(i,j,k), 0);
125 }
126 } else {
127 for (int k = max(0,k1); k <= min(m_Nk-1,k2); ++k) {
128 g->SetValue(getIdx(i,j,k), 1);
129 }
130 }
131 x = x_out;
132 x[2] += shift*m_Dz;
133 k_last = min(m_Nk-1,k2+1);
134 }
135 if (preserveFluid()) {
136 for (int k = k_last; k <= m_Nk-1; ++k) {
137 g->SetValue(getIdx(i,j,k), 0);
138 }
139 }
140 }
141 }
142 }
143 progress += m_Nj;
144 GuiMainWindow::pointer()->setProgress(progress);
145 }
146
147 // scan ik plane -> j direction
148 // xz plane -> y direction
149 for (int i = 0; i < m_Ni; ++i) {
150 for (int k = 0; k < m_Nk; ++k) {
151 vec3_t x0 = getX(i,0,k);
152 for (int i_il = 0; i_il <= interlaces; ++i_il) {
153 for (int k_il = 0; k_il <= interlaces; ++k_il) {
154 vec3_t x = x0;
155 x[0] += i_il*dxi;
156 x[2] += k_il*dzi;
157 int j_last = 0;
158 /*
159 if (fabs(x0[0]) < 1.5 && fabs(x0[2]) < 1.5 && create_grid) {
160 cout << x0 << endl;
161 }
162 */
163 while (shootRay(x, vec3_t(0,1,0), x_in, x_out, n_in, n_out)) {
164 m_GeometryFound = true;
165 m_XScan1[1] = min(m_XScan1[1], x_in[1]);
166 m_XScan2[1] = max(m_XScan2[1], x_out[1]);
167 int j1 = int((x_in[1]-m_X1[1])/m_Dy) + 1;
168 int j2 = int((x_out[1]-m_X1[1])/m_Dy);
169 if (preserveFluid()) {
170 for (int j = j_last; j < min(m_Nj-1,j1); ++j) {
171 g->SetValue(getIdx(i,j,k), 0);
172 }
173 } else {
174 for (int j = max(0,j1); j <= min(m_Nj-1,j2); ++j) {
175 g->SetValue(getIdx(i,j,k), 1);
176 }
177 }
178 x = x_out;
179 x[1] += shift*m_Dy;
180 j_last = min(m_Nj-1,j2+1);
181 }
182 if (preserveFluid()) {
183 for (int j = j_last; j <= m_Nj-1; ++j) {
184 g->SetValue(getIdx(i,j,k), 0);
185 }
186 }
187 }
188 }
189 }
190 progress += m_Nk;
191 GuiMainWindow::pointer()->setProgress(progress);
192 }
193
194 // scan jk plane -> i direction
195 // yz plane -> x direction
196 for (int j = 0; j < m_Nj; ++j) {
197 for (int k = 0; k < m_Nk; ++k) {
198 vec3_t x0 = getX(0,j,k);
199 for (int j_il = 0; j_il <= interlaces; ++j_il) {
200 for (int k_il = 0; k_il <= interlaces; ++k_il) {
201 vec3_t x = x0;
202 x[1] += j_il*dyi;
203 x[2] += k_il*dzi;
204 int i_last = 0;
205 while (shootRay(x, vec3_t(1,0,0), x_in, x_out, n_in, n_out)) {
206 m_GeometryFound = true;
207 m_XScan1[0] = min(m_XScan1[0], x_in[0]);
208 m_XScan2[0] = max(m_XScan2[0], x_out[0]);
209 int i1 = int((x_in[0]-m_X1[0])/m_Dx) + 1;
210 int i2 = int((x_out[0]-m_X1[0])/m_Dx);
211 if (preserveFluid()) {
212 for (int i = i_last; i < min(m_Ni-1,i1); ++i) {
213 g->SetValue(getIdx(i,j,k), 0);
214 }
215 } else {
216 for (int i = max(0,i1); i <= min(m_Ni-1,i2); ++i) {
217 g->SetValue(getIdx(i,j,k), 1);
218 }
219 }
220 x = x_out;
221 x[0] += shift*m_Dx;
222 i_last = min(m_Ni-1,i2+1);
223 }
224 if (preserveFluid()) {
225 for (int i = i_last; i <= m_Ni-1; ++i) {
226 g->SetValue(getIdx(i,j,k), 0);
227 }
228 }
229 }
230 }
231 }
232 progress += m_Nk;
233 GuiMainWindow::pointer()->setProgress(progress);
234 }
235
236 if (create_grid) {
237
238 GuiMainWindow::pointer()->resetProgress("smoothing", m_Ni*m_Nj*m_Nk*m_NumIterations);
239
240 // smooth image data
241 int progress = 0;
242 for (int iter = 1; iter <= m_NumIterations; ++iter) {
243 for (int i = 1; i < m_Ni-1; ++i) {
244 for (int j = 1; j < m_Nj-1; ++j) {
245 for (int k = 1; k < m_Nk-1; ++k) {
246 int idx = getIdx(i,j,k);
247 bool preserve = false;
248 if (m_PreservationType == 1 && g->GetValue(idx) > 0.999) {
249 preserve = true;
250 }
251 if (m_PreservationType == 2 && g->GetValue(idx) < 0.001) {
252 preserve = true;
253 }
254 if (!preserve) {
255 double g_new = 0;
256 g_new += g->GetValue(getIdx(i+1,j,k));
257 g_new += g->GetValue(getIdx(i-1,j,k));
258 g_new += g->GetValue(getIdx(i,j+1,k));
259 g_new += g->GetValue(getIdx(i,j-1,k));
260 g_new += g->GetValue(getIdx(i,j,k+1));
261 g_new += g->GetValue(getIdx(i,j,k-1));
262 g_new *= 1.0/6.0;
263 g->SetValue(idx, g_new);
264 }
265 }
266 }
267 progress += m_Nj*m_Nk;
268 GuiMainWindow::pointer()->setProgress(progress);
269 }
270 }
271
272 // write image data for testing and reporting purposes
273 EG_VTKSP(vtkXMLImageDataWriter, vti);
274 QString file_name = GuiMainWindow::pointer()->getCwd() + "/g.vti";
275 vti->SetFileName(qPrintable(file_name));
276 vti->SetDataModeToBinary();
277 vti->SetInput(gdata);
278 vti->Write();
279
280 gdata->GetPointData()->SetActiveScalars("g");
281 EG_VTKSP(vtkContourFilter, contour);
282 contour->SetInput(gdata);
283 contour->SetNumberOfContours(1);
284 double g_level = 0.5;
285 if (m_PreservationType == 1) {
286 g_level = 0.5;
287 }
288 if (m_PreservationType == 2) {
289 g_level = 0.5;
290 }
291 contour->SetValue(0, g_level);
292 EG_VTKSP(vtkTriangleFilter, tri);
293 tri->SetInput(contour->GetOutput());
294 EG_VTKSP(vtkDecimatePro, decimate);
295 decimate->PreserveTopologyOn();
296 decimate->SetTargetReduction(m_TargetReduction);
297 decimate->SetFeatureAngle(GeometryTools::deg2rad(45));
298 decimate->SetInput(tri->GetOutput());
299 decimate->Update();
300
301 allocateGrid(m_Grid, decimate->GetOutput()->GetNumberOfPolys(), decimate->GetOutput()->GetNumberOfPoints());
302 for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
303 vec3_t x;
304 decimate->GetOutput()->GetPoint(id_node, x.data());
305 x[0] = m_X1[0] + x[0]*m_Dx;
306 x[1] = m_X1[1] + x[1]*m_Dx;
307 x[2] = m_X1[2] + x[2]*m_Dx;
308 m_Grid->GetPoints()->SetPoint(id_node, x.data());
309 }
310 EG_VTKDCC(vtkIntArray, cell_code, m_Grid, "cell_code");
311 for (vtkIdType id_cell = 0; id_cell < decimate->GetOutput()->GetNumberOfPolys(); ++id_cell) {
312 EG_GET_CELL(id_cell, decimate->GetOutput());
313 vtkIdType id_new_cell = m_Grid->InsertNextCell(type_cell, num_pts, pts);
314 cell_code->SetValue(id_new_cell, 1);
315 }
316 }
317 }
318
319 void CreateCadTesselation::operate()
320 {
321 double preview_memory = min(1024*1024*100.0, m_ScanMemory);
322 int N = int(pow(preview_memory/sizeof(float), 1.0/3.0));
323 m_Ni = N;
324 m_Nj = N;
325 m_Nk = N;
326 m_X1 = vec3_t(-1e6,-1e6,-1e6);
327 m_X2 = vec3_t( 1e6, 1e6, 1e6);
328 m_Dx = (m_X2[0] - m_X1[0])/(m_Ni-1);
329 m_Dy = (m_X2[1] - m_X1[1])/(m_Nj-1);
330 m_Dz = (m_X2[2] - m_X1[2])/(m_Nk-1);
331 double new_vol, old_vol;
332 int count = 0;
333 do {
334 old_vol = (m_X2[0]-m_X1[0])*(m_X2[1]-m_X1[1])*(m_X2[2]-m_X1[2]);
335 QString num;
336 QString text = "scanning (V=";
337 num.setNum(old_vol);
338 text += num + ")";
339
340 GuiMainWindow::pointer()->resetProgress(text, 3*N*N);
341 scan(false);
342 if (m_GeometryFound) {
343 m_X1 = m_XScan1 - 2*vec3_t(m_Dx, m_Dy, m_Dz);
344 m_X2 = m_XScan2 + 2*vec3_t(m_Dx, m_Dy, m_Dz);
345 } else {
346 m_X1 *= 0.1;
347 m_X2 *= 0.1;
348 }
349 new_vol = (m_X2[0]-m_X1[0])*(m_X2[1]-m_X1[1])*(m_X2[2]-m_X1[2]);
350 } while (count < 20 && (old_vol - new_vol)/old_vol > 0.05);
351
352 // bounding box should now be established
353 // last scan run with the full resoluion (if required)
354 double Lx = m_X2[0] - m_X1[0];
355 double Ly = m_X2[1] - m_X1[1];
356 double Lz = m_X2[2] - m_X1[2];
357 double max_size = m_ScanMemory/sizeof(float);
358 double delta = max(m_SmallestResolution, pow(Lx*Ly*Lz/max_size, 1.0/3.0));
359 int interlaces = 0;
360 if (preserveFluid() || preserveSolid()) {
361 interlaces = int(2*delta/m_SmallestFeatureSize);
362 }
363 m_Ni = int(Lx/delta) + 1;
364 m_Nj = int(Ly/delta) + 1;
365 m_Nk = int(Lz/delta) + 1;
366 QString num;
367 QString text = "scanning (h=";
368 num.setNum(delta/(interlaces+1));
369 text += num + ")";
370 GuiMainWindow::pointer()->resetProgress(text, m_Ni*m_Nj + m_Ni*m_Nk + m_Nj*m_Nk);
371 scan(true, interlaces);
372
373 UpdateNodeIndex(m_Grid);
374 UpdateCellIndex(m_Grid);
375 GuiMainWindow::pointer()->resetXmlDoc();
376 GuiMainWindow::pointer()->clearBCs();
377 GuiMainWindow::pointer()->addBC(1, BoundaryCondition("imported", "patch"));
378
379 GuiMainWindow::pointer()->resetProgress(" ", 100);
380 }
Clone of the official engrid repository