src/libengrid/pointfinder.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 #include "pointfinder.h"
  24
  25 PointFinder::PointFinder()
  26 {
  27   m_MinSize   = 1.0;
  28   m_MaxPoints = 100;
  29 }
  30
  31 void PointFinder::setGrid(vtkUnstructuredGrid *grid)
  32 {
  33   QVector<vec3_t> points(grid->GetNumberOfPoints());
  34   for (vtkIdType id_node = 0; id_node < grid->GetNumberOfPoints(); ++id_node) {
  35     grid->GetPoint(id_node, points[id_node].data());
  36   }
  37   setPoints(points);
  38 }
  39
  40 void PointFinder::setPoints(const QVector<vec3_t> &points)
  41 {
  42   m_Points = points;
  43   {
  44     vec3_t x1(1e99, 1e99, 1e99);
  45     vec3_t x2(-1e99, -1e99, -1e99);
  46     foreach (vec3_t x, m_Points) {
  47       x1[0] = min(x[0], x1[0]);
  48       x1[1] = min(x[1], x1[1]);
  49       x1[2] = min(x[2], x1[2]);
  50       x2[0] = max(x[0], x2[0]);
  51       x2[1] = max(x[1], x2[1]);
  52       x2[2] = max(x[2], x2[2]);
  53     }
  54     vec3_t xc = 0.5*(x1 + x2);
  55     vec3_t Dx1 = xc - x1;
  56     vec3_t Dx2 = x2 - xc;
  57     if (fabs(Dx1[0]) >= fabs(Dx1[1]) && fabs(Dx1[0]) >= fabs(Dx1[2])) {
  58       Dx1 = vec3_t(Dx1[0], Dx1[0], Dx1[0]);
  59     } else if (fabs(Dx1[1]) >= fabs(Dx1[0]) && fabs(Dx1[1]) >= fabs(Dx1[2])) {
  60       Dx1 = vec3_t(Dx1[1], Dx1[1], Dx1[1]);
  61     } else {
  62       Dx1 = vec3_t(Dx1[2], Dx1[2], Dx1[2]);
  63     }
  64     if (fabs(Dx2[0]) >= fabs(Dx2[1]) && fabs(Dx2[0]) >= fabs(Dx2[2])) {
  65       Dx2 = vec3_t(Dx2[0], Dx2[0], Dx2[0]);
  66     } else if (fabs(Dx2[1]) >= fabs(Dx2[0]) && fabs(Dx2[1]) >= fabs(Dx2[2])) {
  67       Dx2 = vec3_t(Dx2[1], Dx2[1], Dx2[1]);
  68     } else {
  69       Dx2 = vec3_t(Dx2[2], Dx2[2], Dx2[2]);
  70     }
  71     vec3_t Dx = Dx1;
  72     if (Dx2.abs() > Dx1.abs()) {
  73       Dx = Dx2;
  74     }
  75     x1 = xc - 2*Dx;
  76     x2 = xc + 2*Dx;
  77     m_Octree.setBounds(x1, x2);
  78     m_MinSize = 0.0001*Dx[0];
  79   }
  80   m_Buckets.resize(1);
  81   for (int i_points = 0; i_points < points.size(); ++i_points) {
  82     m_Buckets[0].append(i_points);
  83   }
  84   int N;
  85   do {
  86     N = refine();
  87   } while (N > 0);
  88   m_MaxBucketSize = 0;
  89   m_MinBucketSize = points.size();
  90   for (int i = 0; i < m_Buckets.size(); ++i) {
  91     m_MinBucketSize = min(m_MinBucketSize, m_Buckets[i].size());
  92     m_MaxBucketSize = max(m_MaxBucketSize, m_Buckets[i].size());
  93   }
  94 }
  95
  96 int PointFinder::refine()
  97 {
  98   int N = m_Octree.getNumCells();
  99   m_Octree.setSmoothTransitionOff();
 100   int N_new = 0;
 101   m_Octree.resetRefineMarks();
 102   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 103     double dx = m_Octree.getDx(cell);
 104     if (!m_Octree.hasChildren(cell) && m_Buckets[cell].size() > m_MaxPoints  && dx > 2*m_MinSize) {
 105       m_Octree.markToRefine(cell);
 106       ++N_new;
 107     }
 108   }
 109   int N_mrk = 0;
 110   for (int cell = 0; cell < m_Octree.getNumCells(); ++cell) {
 111     if (m_Octree.markedForRefine(cell)) {
 112       ++N_mrk;
 113     }
 114   }
 115   N_new *= 8;
 116   m_Octree.refineAll();
 117   if (m_Octree.getNumCells() != N + N_new) {
 118     EG_BUG;
 119   }
 120   m_Buckets.insert(N, N_new, QList<int>());
 121   for (int cell = N; cell < m_Octree.getNumCells(); ++cell) {
 122     m_Timer << "  " << m_Octree.getNumCells() << " octree cells" << Timer::endl;
 123     if (m_Octree.getNumCells() != N + N_new) {
 124       EG_BUG;
 125     }
 126     int parent = m_Octree.getParent(cell);
 127     if (parent < 0) {
 128       EG_BUG;
 129     }
 130     foreach (int i_points, m_Buckets[parent]) {
 131       vec3_t xcell = m_Octree.getCellCentre(cell);
 132       vec3_t x = m_Points[i_points];
 133       bool append = true;
 134       {
 135         double Dx = m_Octree.getDx(cell);
 136         double Dy = m_Octree.getDy(cell);
 137         double Dz = m_Octree.getDz(cell);
 138         if      (x[0] < xcell[0] - 1.5*Dx || x[0] > xcell[0] + 1.5*Dx) append = false;
 139         else if (x[1] < xcell[1] - 1.5*Dy || x[1] > xcell[1] + 1.5*Dy) append = false;
 140         else if (x[2] < xcell[2] - 1.5*Dz || x[2] > xcell[2] + 1.5*Dz) append = false;
 141       }
 142       if (append) {
 143         m_Buckets[cell].append(i_points);
 144       }
 145     }
 146   }
 147   return N_new;
 148 }
 149
 150 void PointFinder::getClosePoints(vec3_t x, QVector<int> &points, double dist)
 151 {
 152   int cell = m_Octree.findCell(x);
 153   if (cell < 0) {
 154     EG_BUG;
 155   }
 156   if (m_Octree.hasChildren(cell)) {
 157     EG_BUG;
 158   }
 159   while ((m_Buckets[cell].size() == 0  || dist > m_Octree.getDx(cell)) && m_Octree.getParent(cell) >= 0) {
 160     cell = m_Octree.getParent(cell);
 161   }
 162   points.resize(m_Buckets[cell].size());
 163   qCopy(m_Buckets[cell].begin(), m_Buckets[cell].end(), points.begin());
 164
 165 }
 166
 167 void PointFinder::writeOctreeMesh(QString file_name)
 168 {
 169   EG_VTKSP(vtkUnstructuredGrid, otg);
 170   m_Octree.toVtkGrid(otg);
 171   writeGrid(otg, file_name);
 172 }