$search

mesh_region.cc

Go to the documentation of this file.

/*
  James R. Diebel
  Stanford University

  Started: 24 August 2004
  Last revised:

  mesh_vert.cc - class implementation of Mesh class (mesh.hh)

  Depends on:
  - Mesh class (mesh.hh)
*/

#include "bmtk/triangulate.h"
#include "bmtk/mesh.hh"
#include <iostream>
#include <list>
#include <vector>
#include <cstring>

namespace bmtk {

// Region Functions //
void Mesh::findRegionFaceLists() {
  for (int i=0;i<nf;i++) f[i].r->fi.push_back(f[i].i);
}

void Mesh::findRegionProps() {
  renumberRegions();
  clearRegionProps();
  findRegionNormals();
  findRegionVariances();
}

void Mesh::clearRegionProps() {
  for (int i=0;i<nr;i++) {
    r[i].p = 0;
    r[i].ssn2 = 0;
    r[i].vi.clear();
    r[i].tvi.clear();
    r[i].ni.clear();
    r[i].fi.clear();
    r[i].n = 0;
    if (nd) for (int k=0;k<nd;k++) r[i].ns[k] = 0;
    if (nd) r[i].n = r[i].ns[0];
  }
}

void Mesh::findRegionNormals() {
  for (int i=0;i<nf;i++) {
    f[i].r->p++;
    f[i].r->n += f[i].n;
    if (nd) for (int k=0;k<nd;k++) f[i].r->ns[k] += f[i].ns[k];
  }
  for (int i=0;i<nr;i++) {
    r[i].n.normalize();
    if (nd) for (int k=0;k<nd;k++) r[i].ns[k].normalize();
  }
}

void Mesh::findRegionVariances() {
  for (int i=0;i<nf;i++) {
    if (nd) for (int k=0;k<nd;k++)
      f[i].r->ssn2 += f[i].ns[k].dist2To(f[i].r->ns[k]);
    else f[i].r->ssn2 += f[i].n.dist2To(f[i].r->n);
  }
}

void Mesh::renumberRegions() {
  bool temp = po;
  po = false;
  resetFaceFlags();
  nr = 0;
  for (int i=0;i<nf;i++) if (f[i].flag == -1) {
    Region* ri = f[i].r;
    q[0] = i;
    for (int in=1,out=0;in>out;out++) {
      if (f[q[out]].flag == -1 && f[q[out]].r == ri) {
  f[q[out]].flag = nr;
  for (int j=0;j<f[q[out]].nf;j++,in++) {
    q[in] = f[q[out]].f[j]->i;
  }
      }
    }
    nr++;
  }
  for (int i=0;i<nf;i++) f[i].r = &r[f[i].flag];
  po = temp;
}

void Mesh::findRegionSmoothNormals() {
  if (po) cout << "- Performing linear regression on region normals..."
         << flush << endl;

  // create lists of all the boundary faces in a given region
  for (int i=0;i<nf;i++) if (f[i].flag != -1)
    f[i].r->fi.push_back(f[i].i);

  // allocate memory
  float** aa = new float *[9];
  for (int i=0;i<9;i++) aa[i] = new float[9];
  float* xx = new float[9];
  float* bb = new float[9];

  // fill optimization matrix
  for (int i=0;i<nr;i++) {
    // clear out working memory
    for (int j=0;j<9;j++) memset(aa[j],0,9*sizeof(float));
    memset(bb,0,9*sizeof(float));
    memset(xx,0,9*sizeof(float));

    // construct basis
    Vec3d e1 = v[r[i].vi[1]].x - v[r[i].vi[0]].x;
    e1 -= e1.dot(r[i].n)*r[i].n;
    e1.normalize();
    Vec3d e2 = e1.cross(r[i].n);
    e2.normalize();

    // fill problem matrix-vector system
    for (int k=0;k<3;k++) aa[k*3][k*3] = 3*r[i].fi.size();
    for (int j=0;j<r[i].fi.size();j++) for (int t=0;t<3;t++) {
      // project the current vertex onto basis
      Vec3d x2D = Vec3d(f[r[i].fi[j]].v[t]->x.dot(e1),
      f[r[i].fi[j]].v[t]->x.dot(e2),0.0f);
      for (int k=0;k<3;k++) {
  bb[3*k] += f[r[i].fi[j]].n[k];
  for (int l=0;l<2;l++) {
    aa[3*k+1+l][3*k] += x2D[l];
    aa[3*k][3*k+1+l] += x2D[l];
    for (int s=0;s<2;s++) aa[3*k+1+s][3*k+1+l] += x2D[l]*x2D[s];
    bb[3*k+1+l] += f[r[i].fi[j]].n[k]*x2D[l];
  }
      }
    }

    // compute linear regression
    gelimd2(aa,bb,xx,9);

    // grab parameters
    Vec3d n0(xx[0],xx[3],xx[6]);
    Mat3x3 A(xx[1],xx[2],0.0f,xx[4],xx[5],0.0f,xx[7],xx[8],0.0f);

    // define vertex normals according to fit
    for (int j=0;j<r[i].tvi.size();j++) {
      Vec3d x2D = Vec3d(v[r[i].tvi[j]].x.dot(e1),
      v[r[i].tvi[j]].x.dot(e2),0.0f);
      r[i].ni.push_back(n0 + A*x2D);
      r[i].ni.back().normalize();
    }

    // print some things
    /*cout << endl;
    for (int j=0;j<9;j++) {
      for (int k=0;k<9;k++) {
  printf("%8.1f",aa[j][k]);
      }
      printf(" | %8.1f | %8.1f\n",xx[j],bb[j]);
    }
    cout << endl << "region " << i << ": " << r[i].n << endl;
    cout << n0 << endl << A << endl;
    cout << (n0 + A*(f[r[i].fi[0]].x)) << endl;*/
  }

  // clean up
  for (int i=0;i<9;i++) delete [] aa[i];
  delete [] xx;
  delete [] bb;

  if (po) cout << "  Done." << flush << endl;
}

void Mesh::findRegionSmoothNormalsOld() {
  // create lists of all the boundary faces in a given region
  for (int i=0;i<nf;i++) if (f[i].flag != -1) f[i].r->fi.push_back(f[i].i);

  // allocate memory
  float** aa = new float *[12];
  for (int i=0;i<12;i++) aa[i] = new float[12];
  float* xx = new float[12];
  float* bb = new float[12];

  // fill optimization matrix
  for (int i=0;i<nr;i++) {
    // clear out working memory
    for (int j=0;j<12;j++) memset(aa[j],0,12*sizeof(float));
    memset(bb,0,12*sizeof(float));
    memset(xx,0,12*sizeof(float));

    // fill problem matrix-vector system
    for (int k=0;k<3;k++) aa[k*4][k*4] = r[i].p;
    for (int j=0;j<r[i].fi.size();j++) for (int t=0;t<3;t++) {
      for (int k=0;k<3;k++) {
  bb[4*k] += f[r[i].fi[j]].n[k];
  for (int l=0;l<3;l++) {
    aa[4*k+1+l][4*k] += f[r[i].fi[j]].v[t]->x[l];
    aa[4*k][4*k+1+l] += f[r[i].fi[j]].v[t]->x[l];
    for (int s=0;s<3;s++) {
      aa[4*k+1+s][4*k+1+l] += (f[r[i].fi[j]].v[t]->x[l]*
             f[r[i].fi[j]].v[t]->x[s]);
    }
    bb[4*k+1+l] += f[r[i].fi[j]].n[k]*f[r[i].fi[j]].v[t]->x[l];
  }
      }
    }

    // compute linear regression
    gelimd2(aa,bb,xx,12);

    // grab parameters
    Vec3d n0(xx[0],xx[4],xx[8]);
    Mat3x3 A(xx[1],xx[2],xx[3],xx[5],xx[6],xx[7],xx[9],xx[10],xx[11]);

    // define vertex normals according to fit
    for (int j=0;j<r[i].tvi.size();j++) {
      r[i].ni.push_back(n0 + A*(v[r[i].tvi[j]].x));
      r[i].ni.back().normalize();
    }

    // print some things
    /*cout << endl;
    for (int j=0;j<12;j++) {
      for (int k=0;k<12;k++) {
  printf("%8.1f",aa[j][k]);
      }
      printf(" | %8.1f | %8.1f\n",xx[j],bb[j]);
    }
    cout << endl << "region " << i << ": " << r[i].n << endl;
    cout << n0 << endl << A << endl;
    cout << (n0 + A*(f[r[i].fi[0]].x)) << endl;*/
  }

  // clean up
  for (int i=0;i<12;i++) delete [] aa[i];
  delete [] xx;
  delete [] bb;
}

void Mesh::doSegmentationStep() {
  for (int i=0;i<nf;i++) {
    for (int j=0;j<f[i].nf;j++) if (f[i].r != f[i].f[j]->r) {
      float cssn2 = 0;
      int pop = nd*(f[i].r->p + f[i].f[j]->r->p);
      Vec3d mn;
      if (nd) for (int k=0;k<nd;k++) {
  mn = (f[i].r->ns[k]*f[i].r->p
        + f[i].f[j]->r->ns[k]*f[i].f[j]->r->p);
  mn.normalize();
  cssn2 += f[i].r->ssn2 + f[i].r->p*mn.dist2To(f[i].r->ns[k]) +
    f[i].f[j]->r->ssn2 + f[i].f[j]->r->p*mn.dist2To(f[i].f[j]->r->ns[k]);
      } else {
  mn = (f[i].r->n*f[i].r->p + f[i].f[j]->r->n*f[i].f[j]->r->p);
  mn.normalize();
  cssn2 = f[i].r->ssn2 + f[i].r->p*mn.dist2To(f[i].r->n) +
    f[i].f[j]->r->ssn2 + f[i].f[j]->r->p*mn.dist2To(f[i].f[j]->r->n);
  pop = f[i].r->p + f[i].f[j]->r->p;
      }
      if (cssn2/float(pop) < segThresh*float(nd+1)) {
  if (f[i].r->p > f[i].f[j]->r->p) {
    mergeRegions(*(f[i].f[j]->r),*(f[i].r),*(f[i].f[j]));
    if (nd) for (int k=0;k<nd;k++) {
      f[i].r->ns[k] = (f[i].r->ns[k]*f[i].r->p
       + f[i].f[j]->r->ns[k]*f[i].f[j]->r->p);
      f[i].r->ns[k].normalize();;
    } else f[i].r->n = mn;
    f[i].r->ssn2 = cssn2;
  } else {
    mergeRegions(*(f[i].r),*(f[i].f[j]->r),f[i]);
    if (nd) for (int k=0;k<nd;k++) {
      f[i].f[j]->r->ns[k] = (f[i].r->ns[k]*f[i].r->p
           + f[i].f[j]->r->ns[k]*f[i].f[j]->r->p);
      f[i].f[j]->r->ns[k].normalize();
    } else f[i].f[j]->r->n = mn;
    f[i].f[j]->r->ssn2 = cssn2;
  }
  continue;
      }

      float diff = 0;
      if (nd) for (int k=0;k<nd;k++)
  diff += f[i].ns[k].dot(f[i].f[j]->r->ns[k]) -
    f[i].ns[k].dot(f[i].r->ns[k]);
      else diff = f[i].n.dot(f[i].f[j]->r->n) - f[i].n.dot(f[i].r->n);
      if (diff > 0) {
  swapFace(f[i],*(f[i].r),*(f[i].f[j]->r));
  break;
      }
    }
  }
  findRegionProps();
}

// performs segementation
void Mesh::runSegmentation() {
  if (po) cout << endl << "Performing mesh segmentation..." << flush << endl;
  for (int i=0;i<maxIter;i++) {
    cout << "- iteration " << i << " of " << maxIter << ": from "
   << nr << " to " << flush;
    doSegmentationStep();
    cout << nr << " regions." << flush << endl;
    update();
  }

  if (po) cout << "Done with Segmentation." << flush << endl;
  findVertNormals();
  findBoundaries();
  triangulateContours();
  findRegionSmoothNormalsOld();
}

// swaps Face fi from Region ri to Region rj
void Mesh::swapFace(Face& fi, Region& ri, Region& rj) {
  if (!nd) {
    ri.ssn2 -= fi.n.dist2To(ri.n);
    ri.n *= ri.p;
    ri.n -= fi.n;
    ri.n.normalize();

    rj.n *= rj.p;
    rj.n += fi.n;
    rj.n.normalize();
    rj.ssn2 += fi.n.dist2To(rj.n);
  }

  ri.p--;
  rj.p++;
  fi.r = &rj;
}

// merge Region ri into Region rj starting at Face fi
void Mesh::mergeRegions(Region& ri, Region& rj, Face& fi) {
  rj.p += ri.p;
  ri.p = 0;
  q[0] = fi.i;
  for (int in=1,out=0;in>out;out++) if (f[q[out]].r == &ri) {
    f[q[out]].r = &rj;
    for (int j=0;j<f[q[out]].nf;j++,in++) q[in] = f[q[out]].f[j]->i;
  }
  nr--;
}

// find boundaries between regions, keypoints
void Mesh::findBoundaries() {
  // merge single-triangle floaters with surrounding region
  for (int i=0;i<nf;i++) if (f[i].r->p == 1) {
    bool floater = false;
    if (f[i].nf > 1) {
      floater = true;
      for (int j=1;j<f[i].nf;j++) {
  floater = floater && (f[i].f[j]->r == f[i].f[j-1]->r);
      }
    }
    if (floater) mergeRegions(*(f[i].r),*(f[i].f[0]->r),f[i]);
  }

  // recalculate region properties
  findRegionProps();

  // Mark all keypoints with current regions
  markKeypoints();

  // scan through all faces, looking for edges to trace
  for (int i=0;i<nf;i++) if (f[i].flag != -1) {
    for (int j=0;j<3;j++) {
      // if the edge is a region boundary and has not been visited yet
      if ((f[i].e[j]->flag != -1) && (f[i].e[j]->flag != f[i].r->i) &&
    f[i].r->vi.size() == 0){
  traceBoundary(f[i].e[j]->v[0],f[i].e[j],&f[i],f[i].r);
      }
    }
  }
}

void Mesh::markKeypoints() {
  resetEdgeFlags();
  resetVertFlags();
  resetFaceFlags();

  // flag edges, faces as boundaries, vertices as keypoints
  for (int i=0;i<ne;i++) {
    switch (e[i].nf) {
    case 1:
      e[i].flag = -2;
      e[i].f[0]->flag = -2;
      e[i].v[1]->flag-=2;
      e[i].v[0]->flag-=2;
      break;
    case 2:
      if (e[i].f[0]->r != e[i].f[1]->r) {
  e[i].flag =  -3;
  e[i].f[0]->flag = e[i].f[1]->flag = -3;
  e[i].v[0]->flag-=1;
  e[i].v[1]->flag-=1;
      }
      break;
    }
  }
}

void Mesh::traceBoundary(Vert* vc, Edge* ec, Face* fc, Region* rc) {
  list<int> trace;
  list<int>::iterator it;
  bool done,found;

  do {
    // add the current contour
    trace.push_back(vc->i);

    // move current vertex along current edge to next vertex
    if (vc == ec->v[0]) vc = ec->v[1]; else vc = ec->v[0];
    ec->flag = rc->i;

    // if new vertex has been visited before, delete everything between first
    // and second visits to ensure simply connected region
    if (vc->marker == rc->i) while (trace.back() != vc->i) trace.pop_back();

    // mark new vertex as visited (note that we leave the first vertex
    // unmarked because we wish to return to it without causing a problem
    vc->marker = rc->i;

    // look for next edge by traversing neighboring faces
    do {
      // assume that we have nowhere else to go
      done = true;
      found = false;

      // check for valid edge in current face
      for (int i=0;i<3;i++) {
  //cout << "." << flush;
  if (fc->e[i]->flag != -1 &&
      fc->e[i]->flag != rc->i &&
      *vc < *(fc->e[i])) {
    ec = fc->e[i];
    found = true; done = false;
    break;
  }
      }
      if (found) break;

      // otherwise try a neighboring face that also touches current vertex
      // and is part of the current region
      fc->flag = vc->i;
      for (int i=0;i<fc->nf;i++) {
  if (fc->f[i]->flag != vc->i &&
      fc->f[i]->r == rc &&
      *vc < *(fc->f[i])) {
    fc = fc->f[i];
    done = false;
    break;
  }
      }
    } while (!done);
  } while (!done);

  // traverse trace, filling the region contour with keypoints only
  trace.unique();
  for (it = trace.begin();it!=trace.end();it++) if (v[*it].flag < -3)
    rc->vi.push_back(*it);

  // terminate the chain with a flag
  rc->vi.push_back(-1);

  // print the contour
  /*cout << "Region " << rc->i << ": " << flush;
    for (int i=0;i<rc->vi.size();i++)
    cout << rc->vi[i] << " ";
    cout << endl << flush;*/
}

// generates decimated mesh after call to findBoundaries()
void Mesh::triangulateContours() {
  if (po) cout << "Triangulating region contours..." << endl << flush;

  // some STL lists used by the triangulation routines
  Vector2dVector contour;
  VertexList mappings,triangles;

  // for each (polygon <-> region), perform tessellation
  for (int i=0;i<nr;i++) {
    if (r[i].vi.size() < 4) continue;

    // construct local 2D basis
    Vec3d e1 = v[r[i].vi[1]].x - v[r[i].vi[0]].x;
    e1 -= e1.dot(r[i].n)*r[i].n;
    e1.normalize();
    Vec3d e2 = e1.cross(r[i].n);

    // grab contours and tessellate
    contour.clear();
    mappings.clear();
    for (int j=0;j<r[i].vi.size();j++) {
      if (r[i].vi[j] != -1) {
  float x1 = v[r[i].vi[j]].x.dot(e1);
  float x2 = v[r[i].vi[j]].x.dot(e2);
  contour.push_back(Vector2d(x1,x2));
  mappings.push_back(r[i].vi[j]);
      }
      else {
  if (mappings.size() < 3) {
    contour.clear();
    mappings.clear();
    break;
  }
  // and triangulate that contour, storing results in triangles
  Triangulate::Process(contour,triangles);

  // finally, remap the recent additions to our triangles list to respect
  // the global definitions of the vertices
  r[i].tvi.clear();
  for (int k=0;k<triangles.size();k++) {
    r[i].tvi.push_back(mappings[triangles[k]]);
  }
  contour.clear();
  mappings.clear();
  triangles.clear();
      }
    }
  }

  if (po) cout << "Done triangulating contour." << endl << flush;
}

// generates decimated mesh after call to findBoundaries()
Mesh* Mesh::genDecimatedMesh() {
  if (po) cout << "Generating decimated mesh..." << endl << flush;
  // count vertices in decimated mesh and create links
  int numvertices = 0;
  for (int i=0;i<nv;i++) if (v[i].flag < -3) {
    v[i].flag = numvertices;
    numvertices++;
  } else v[i].flag = -1;

  // create new vertex list and fill it from source mesh
  float* vertices = new float[3*numvertices];
  for (int i=0;i<nv;i++) if (v[i].flag >= 0) {
    for (int j=0;j<3;j++) {
      vertices[3*v[i].flag+j] = v[i].x[j];
    }
  }

  // some STL lists used by the triangulation routines
  Vector2dVector contour;
  VertexList mappings,triangles,polygons;

  // for each (polygon <-> region), perform tessellation
  for (int i=0;i<nr;i++) {
    if (r[i].vi.size() < 4) continue;

    // construct local 2D basis
    Vec3d e1 = v[r[i].vi[1]].x - v[r[i].vi[0]].x;
    e1 -= e1.dot(r[i].n)*r[i].n;
    e1.normalize();
    Vec3d e2 = e1.cross(r[i].n);

    // grab contours and tessellate
    contour.clear();
    mappings.clear();
    for (int j=0;j<r[i].vi.size();j++) {
      if (r[i].vi[j] != -1) {
  float x1 = v[r[i].vi[j]].x.dot(e1);
  float x2 = v[r[i].vi[j]].x.dot(e2);
  contour.push_back(Vector2d(x1,x2));
  mappings.push_back(v[r[i].vi[j]].flag);
      }
      else {
  if (mappings.size() < 3) {
    contour.clear();
    mappings.clear();
    break;
  }
  // and triangulate that contour, storing results in triangles
  int n = triangles.size();
  Triangulate::Process(contour,triangles);

  // finally, remap the recent additions to our triangles list to respect
  // the global definitions of the vertices
  for (int k=n;k<triangles.size();k++) {
    triangles[k] = mappings[triangles[k]];
    polygons.push_back(i);
  }
  contour.clear();
  mappings.clear();
      }
    }
  }

  // create and fill up triangle list
  int numfaces = triangles.size()/3;
  int* faces = new int[3*numfaces];
  for (int i=0;i<numfaces;i++) {
    faces[3*i] = triangles[3*i+1];
    faces[3*i+1] = triangles[3*i];
    faces[3*i+2] = triangles[3*i+2];
  }

  // Build new mesh
  Mesh* dm = new Mesh();
  vector<int> badFaces = dm->buildFrom(vertices,numvertices,faces,numfaces);

  // Allocate space for other datasets, if needed
  if (nd) {
    dm->nd = nd;
    for (int i=0;i<dm->nv;i++) dm->v[i].xs = new Vec3d[nd];
    for (int i=0;i<dm->nf;i++) {
      dm->f[i].ns = new Vec3d[nd];
      dm->r[i].ns = new Vec3d[nd];
    }
  }

  // Go through the source vertices and fill new mesh
  if (nd) for (int i=0;i<nv;i++) if (v[i].flag >= 0) {
    for (int k=0;k<nd;k++) dm->v[v[i].flag].xs[k] = v[i].xs[k];
  }

  // Fix the normal vectors to be the same as that of the parent region
  // and add in data from other data sets if available
  if (badFaces.size() == 0) {
    for (int i=0;i<dm->nf;i++) {
      dm->f[i].n = r[polygons[i*3]].n;
      if (nd) for (int ii=0;ii<nd;ii++)
  dm->f[i].ns[ii] = r[polygons[i*3]].ns[ii];
    }
  } else {
    for (int i=0,j=0,k=0;i<numfaces;i++) {
      if (j<badFaces.size()) if (i==badFaces[j]) {
  j++; continue;
      }
      dm->f[k].n = r[polygons[i*3]].n;
      if (nd) for (int ii=0;ii<nd;ii++) {
  dm->f[k].ns[ii] = r[polygons[i*3]].ns[ii];
      }
      k++;
    }
  }

  // Clean up and get out of here!
  delete [] vertices;
  delete [] faces;

  if (po) cout << "Done building decimated mesh." << endl << flush;

  return dm;
}

} // namespace bmtk

---

bmtk
Author(s): Benjamin Pitzer
autogenerated on Mon Mar 4 11:05:51 2013