src/libengrid/reducedpolydatareader.cpp

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