src/libengrid/reducedpolydatareader.cpp

   1 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   2 // +                                                                      +
   3 // + This file is part of enGrid.                                         +
   4 // +                                                                      +
   5 // + Copyright 2008-2014 enGits GmbH                                      +
   6 // +                                                                      +
   7 // + enGrid is free software: you can redistribute it and/or modify       +
   8 // + it under the terms of the GNU General Public License as published by +
   9 // + the Free Software Foundation, either version 3 of the License, or    +
  10 // + (at your option) any later version.                                  +
  11 // +                                                                      +
  12 // + enGrid is distributed in the hope that it will be useful,            +
  13 // + but WITHOUT ANY WARRANTY; without even the implied warranty of       +
  14 // + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        +
  15 // + GNU General Public License for more details.                         +
  16 // +                                                                      +
  17 // + You should have received a copy of the GNU General Public License    +
  18 // + along with enGrid. If not, see <http://www.gnu.org/licenses/>.       +
  19 // +                                                                      +
  20 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  21
  22 #include "reducedpolydatareader.h"
  23 #include "engrid.h"
  24 #include "guimainwindow.h"
  25 #include "octree.h"
  26
  27 #include <QFileInfo>
  28
  29 #include <vtkPolyDataReader.h>
  30 #include <vtkEgPolyDataToUnstructuredGridFilter.h>
  31 #include <vtkDelaunay3D.h>
  32 #include <vtkClipVolume.h>
  33 #include <vtkContourFilter.h>
  34 #include <vtkGeometryFilter.h>
  35 #include <vtkDecimatePro.h>
  36 #include <vtkSmoothPolyDataFilter.h>
  37
  38 ReducedPolyDataReader::ReducedPolyDataReader()
  39 {
  40   setFormat("legacy VTK files(*.vtk)");
  41   setExtension(".vtk");
  42   m_MaxEdgeLength = 10;
  43 }
  44
  45 void ReducedPolyDataReader::computeLevelSet(vtkUnstructuredGrid* m_Grid, vtkPolyData* poly)
  46 {
  47   // create triangles
  48   poly->BuildCells();
  49   QVector<Triangle> triangles(poly->GetNumberOfPolys());
  50   for (vtkIdType id_poly = 0; id_poly < poly->GetNumberOfPolys(); ++id_poly) {
  51     EG_GET_CELL(id_poly, poly);
  52     if (num_pts == 3) {
  53       poly->GetPoint(pts[0], triangles[id_poly].a.data());
  54       poly->GetPoint(pts[1], triangles[id_poly].b.data());
  55       poly->GetPoint(pts[2], triangles[id_poly].c.data());
  56       triangles[id_poly].id_a = pts[0];
  57       triangles[id_poly].id_b = pts[1];
  58       triangles[id_poly].id_c = pts[2];
  59       triangles[id_poly].g1 = triangles[id_poly].b - triangles[id_poly].a;
  60       triangles[id_poly].g2 = triangles[id_poly].c - triangles[id_poly].a;
  61       triangles[id_poly].g3 = triangles[id_poly].g1.cross(triangles[id_poly].g2);
  62       triangles[id_poly].A  = 0.5*triangles[id_poly].g3.abs();
  63       triangles[id_poly].g3.normalise();
  64       triangles[id_poly].G.column(0, triangles[id_poly].g1);
  65       triangles[id_poly].G.column(1, triangles[id_poly].g2);
  66       triangles[id_poly].G.column(2, triangles[id_poly].g3);
  67       triangles[id_poly].GI = triangles[id_poly].G.inverse();
  68       triangles[id_poly].smallest_length = (triangles[id_poly].b - triangles[id_poly].a).abs();
  69       triangles[id_poly].smallest_length = min(triangles[id_poly].smallest_length, (triangles[id_poly].c - triangles[id_poly].b).abs());
  70       triangles[id_poly].smallest_length = min(triangles[id_poly].smallest_length, (triangles[id_poly].a - triangles[id_poly].c).abs());
  71     } else {
  72       EG_BUG;
  73     }
  74   }
  75
  76   //vtkDoubleArray *scalars = vtkDoubleArray::New();
  77   EG_VTKSP(vtkDoubleArray, scalars);
  78   scalars->SetNumberOfComponents(1);
  79   scalars->SetName("g");
  80   scalars->Allocate(m_Grid->GetNumberOfPoints());
  81
  82   QProgressDialog progress("Reducing triangulation", "Abort", 0, m_Grid->GetNumberOfPoints());
  83   progress.setWindowModality(Qt::ApplicationModal);
  84   for (vtkIdType id_node = 0; id_node < m_Grid->GetNumberOfPoints(); ++id_node) {
  85     progress.setValue(id_node);
  86     QApplication::processEvents();
  87     if (progress.wasCanceled()) {
  88       EG_ERR_RETURN("interrupted by user");
  89     }
  90     double g_levelset = 1e99;
  91     vec3_t xp;
  92     m_Grid->GetPoint(id_node, xp.data());
  93     foreach (Triangle T, triangles) {
  94       vec3_t xi(1e99,1e99,1e99);
  95       vec3_t ri;
  96       double scal = (xp - T.a)*T.g3;
  97       vec3_t x1, x2;
  98       if (scal > 0) {
  99         x1 = xp + T.g3;
 100         x2 = xp - scal*T.g3 - T.g3;
 101       } else {
 102         x1 = xp - T.g3;
 103         x2 = xp - scal*T.g3 + T.g3;
 104       }
 105       double d = 1e99;
 106
 107       bool intersects_face = GeometryTools::intersectEdgeAndTriangle(T.a, T.b, T.c, x1, x2, xi, ri);
 108       if (intersects_face) {
 109         vec3_t dx = xp - xi;
 110         d = fabs(dx*T.g3);
 111       } else {
 112         double kab = GeometryTools::intersection(T.a, T.b - T.a, xp, T.b - T.a);
 113         double kac = GeometryTools::intersection(T.a, T.c - T.a, xp, T.c - T.a);
 114         double kbc = GeometryTools::intersection(T.b, T.c - T.b, xp, T.c - T.b);
 115         double dab = (T.a + kab*(T.b-T.a) - xp).abs();
 116         double dac = (T.a + kac*(T.c-T.a) - xp).abs();
 117         double dbc = (T.b + kbc*(T.c-T.b) - xp).abs();
 118         bool set = false;
 119         if ((kab >= 0) && (kab <= 1)) {
 120           if (dab < d) {
 121             xi = T.a + kab*(T.b-T.a);
 122             d = dab;
 123             set = true;
 124           }
 125         }
 126         if ((kac >= 0) && (kac <= 1)) {
 127           if (dac < d) {
 128             xi = T.a + kac*(T.c-T.a);
 129             d = dac;
 130             set = true;
 131           }
 132         }
 133         if ((kbc >= 0) && (kbc <= 1)) {
 134           if (dbc < d) {
 135             xi = T.b + kbc*(T.c-T.b);
 136             d = dbc;
 137             set = true;
 138           }
 139         }
 140         double da = (T.a - xp).abs();
 141         double db = (T.b - xp).abs();
 142         double dc = (T.c - xp).abs();
 143         if (da < d) {
 144           xi = T.a;
 145           d = da;
 146           set = true;
 147         }
 148         if (db < d) {
 149           xi = T.b;
 150           d = db;
 151         }
 152         if (dc < d) {
 153           xi = T.c;
 154           d = dc;
 155           set = true;
 156         }
 157         if (!set) {
 158           EG_BUG;
 159         }
 160       }
 161       if (xi[0] > 1e98) {
 162         EG_BUG;
 163       }
 164       vec3_t dx = xp - xi;
 165       if (d < fabs(g_levelset)) {
 166         if (dx*T.g3 > 0) {
 167           g_levelset = d;
 168         } else {
 169           g_levelset = -d;
 170         }
 171       }
 172     }
 173     scalars->InsertTuple1(id_node, g_levelset);
 174   }
 175   progress.setValue(m_Grid->GetNumberOfPoints());
 176   m_Grid->GetPointData()->SetScalars(scalars);
 177 }
 178
 179
 180 void ReducedPolyDataReader::operate()
 181 {
 182   try {
 183     QFileInfo file_info(GuiMainWindow::pointer()->getFilename());
 184     readInputFileName(file_info.completeBaseName() + ".vtk");
 185     if (isValid()) {
 186       EG_VTKSP(vtkPolyDataReader, vtk);
 187       double bounds[6];
 188       vtk->SetFileName(qPrintable(getFileName()));
 189       vtk->Update();
 190       vtk->GetOutput()->GetBounds(bounds);
 191       vec3_t x1(bounds[0], bounds[2], bounds[4]);
 192       vec3_t x2(bounds[1], bounds[3], bounds[5]);
 193       double m_MaxEdgeLength = (x1-x2).abs()/20;
 194 //      double m_MinEdgeLength = (x1-x2).abs()/100;
 195       vec3_t x12 = 0.5*(x1 + x2);
 196       vec3_t dx = x2 - x1;
 197       double d = max(dx[0], max(dx[1], dx[2]));
 198       dx = vec3_t(d,d,d);
 199       x1 = x12 - 0.75*dx;
 200       x2 = x12 + 0.75*dx;
 201       Octree octree;
 202       octree.setBounds(x1, x2);
 203       octree.setSmoothTransitionOn();
 204
 205       EG_VTKSP(vtkSmoothPolyDataFilter, smooth1);
 206       smooth1->SetInputConnection(vtk->GetOutputPort());
 207       smooth1->FeatureEdgeSmoothingOn();
 208       smooth1->BoundarySmoothingOn();
 209       smooth1->SetNumberOfIterations(200);
 210       smooth1->Update();
 211       smooth1->GetOutput()->BuildCells();
 212
 213       /*
 214       QVector<double> max_length(smooth->GetOutput()->GetNumberOfPoints(), m_MaxEdgeLength);
 215       for (vtkIdType id_node = 0; id_node < smooth->GetOutput()->GetNumberOfPolys(); ++id_node) {
 216         vtkIdType N_pts, *pts;
 217         smooth1->GetOutput()->GetCellPoints(id_node, N_pts, pts);
 218         vec3_t x[N_pts];
 219         for (int i = 0; i < N_pts; ++i) {
 220           smooth1->GetOutput()->GetPoint(pts[i], x[i].data());
 221         }
 222         for (int i = 0; i < N_pts-1; ++i) {
 223           int j = i+1;
 224           if (j >= N_pts) {
 225             j = 0;
 226           }
 227           double L = (x[i]-x[j]).abs();
 228           max_length[pts[i]] = max(max_length[pts[i]], L);
 229           max_length[pts[j]] = max(max_length[pts[j]], L);
 230         }
 231       }
 232       */
 233
 234       EG_VTKSP(vtkDecimatePro, decimate);
 235       decimate->SetFeatureAngle(20.0);
 236       decimate->PreserveTopologyOn();
 237       decimate->SetInputConnection(smooth1->GetOutputPort());
 238       decimate->SetTargetReduction(0.9);
 239       EG_VTKSP(vtkSmoothPolyDataFilter, smooth2);
 240       smooth2->SetInputConnection(decimate->GetOutputPort());
 241       smooth2->FeatureEdgeSmoothingOn();
 242       smooth2->BoundarySmoothingOn();
 243       smooth2->SetNumberOfIterations(200);
 244       smooth2->Update();
 245
 246       int N;
 247       do {
 248         for (vtkIdType id_node = 0; id_node < smooth2->GetOutput()->GetNumberOfPoints(); ++id_node) {
 249           vec3_t x;
 250           smooth2->GetOutput()->GetPoint(id_node, x.data());
 251           int cell = octree.findCell(x);
 252           if (octree.getDx(cell) >m_MaxEdgeLength) {
 253             octree.markToRefine(cell);
 254           }
 255         }
 256         N = octree.refineAll();
 257       } while (N > 0);
 258
 259       EG_VTKSP(vtkUnstructuredGrid, oct_grid);
 260       octree.toVtkGridHangingNodes(oct_grid);
 261
 262       computeLevelSet(oct_grid, smooth2->GetOutput());
 263
 264       EG_VTKSP(vtkDelaunay3D, delaunay);
 265       delaunay->SetInputData(oct_grid);
 266
 267       EG_VTKSP(vtkContourFilter, contour);
 268       contour->SetInputConnection(delaunay->GetOutputPort());
 269       contour->GenerateValues(1, 0, 0);
 270       contour->SetInputData(oct_grid);
 271       contour->Update();
 272
 273       {
 274         EG_VTKSP(vtkXMLUnstructuredGridWriter,vtu);
 275         QString file_name = GuiMainWindow::pointer()->getCwd() + "/oct_grid.vtu";
 276         vtu->SetFileName(qPrintable(file_name));
 277         vtu->SetDataModeToBinary();
 278         vtu->SetInputData(oct_grid);
 279         vtu->Write();
 280       }
 281       if (contour->GetOutput()->GetNumberOfPolys() == 0) {
 282         EG_ERR_RETURN("unable to extract reduced surface");
 283       }
 284
 285       EG_VTKSP(vtkEgPolyDataToUnstructuredGridFilter, poly2ug);
 286       poly2ug->SetInputConnection(contour->GetOutputPort());
 287       poly2ug->Update();
 288
 289       makeCopy(poly2ug->GetOutput(), m_Grid);
 290       createBasicFields(m_Grid, m_Grid->GetNumberOfCells(), m_Grid->GetNumberOfPoints());
 291       updateNodeIndex(m_Grid);
 292       updateCellIndex(m_Grid);
 293       EG_VTKDCC(vtkIntArray, bc, m_Grid, "cell_code");
 294       for (vtkIdType id_cell = 0; id_cell < m_Grid->GetNumberOfCells(); ++id_cell) {
 295         bc->SetValue(id_cell, 0);
 296       }
 297     }
 298   } catch (Error err) {
 299     err.display();
 300   }
 301 }
 302