Off3DFormatTriangularPicture.cpp

#include "Off3DFormatPicture.h"
#include <string>
#include <fstream>
// Draw triangle with vertices {a, a+u, a+v} on grid
// For every entry in grid (except the null vector -1) add a +1 to features
// vector.
void drawTriangleOFF(SparseGrid &grid, int inputFieldSize,
                     std::vector<float> &features, int &nSpatialSites, float a0,
                     float a1, float a2, float u0, float u1, float u2, float v0,
                     float v1, float v2) {
  float base = powf(u0 * u0 + u1 * u1 + u2 * u2, 0.5);
  u0 /= base;
  u1 /= base;
  u2 /= base;                                 // scale u to a unit vector
  float offset = u0 * v0 + u1 * v1 + u2 * v2; // u dot v
  v0 -= offset * u0;
  v1 -= offset * u1;
  v2 -= offset * u2; // make v orthogonal to u
  float height = powf(v0 * v0 + v1 * v1 + v2 * v2, 0.5);
  v0 /= height;
  v1 /= height;
  v2 /= height; // scale v to be a unit vector
  // u and v are now orthogonal
  // The triangle now has points {a, a+base*u, a+offset*u+height*v}

  for (float h = 0; h <= height; h = std::min(h + 1, height) + (h == height)) {
    float l = base * (1 - h / height);
    for (float b = 0; b <= l; b = std::min(b + 1, l) + (b == l)) {
      float p0 = a0 + (b + offset * h / height) * u0 + h * v0,
            p1 = a1 + (b + offset * h / height) * u1 + h * v1,
            p2 = a2 + (b + offset * h / height) * u2 + h * v2;
      int i0 = p0 + 0.25 * inputFieldSize, //(0.25,0.25,0.25) for the pyramid
          // center in cubic coord system,
          // corresponds to
          //(1/2,1/2/sqrt(3),1/sqrt(24))
          i1 = p1 + 0.25 * inputFieldSize, i2 = p2 + 0.25 * inputFieldSize;
      if (i0 >= 0 && i1 >= 0 && i2 >= 0 && i0 + i1 + i2 < inputFieldSize) {
        int n = i0 * inputFieldSize * inputFieldSize + i1 * inputFieldSize + i2;
        if (grid.mp.find(n) == grid.mp.end()) {
          grid.mp[n] = nSpatialSites++;
          features.push_back(1);
        }
      }
      // else
      //   std::cout << "!"<<std::flush;
    }
  }
}

OffSurfaceModelPicture::OffSurfaceModelPicture(std::string filename,
                                               int renderSize, int label)
    : Picture(label), renderSize(renderSize) {
  std::ifstream file(filename.c_str());
  std::string off;
  getline(file, off);
  int nPoints, nTriangles, nLines;
  file >> nPoints >> nTriangles >> nLines;
  points.set_size(nPoints, 3);
  surfaces.resize(nTriangles);
  for (int i = 0; i < nPoints; i++)
    file >> points(i, 0) >> points(i, 1) >> points(i, 2);
  for (int i = 0; i < nTriangles; i++) {
    surfaces[i].resize(3);
    int three;
    file >> three >> surfaces[i][0] >> surfaces[i][1] >> surfaces[i][2];
  }
}

OffSurfaceModelPicture::~OffSurfaceModelPicture() {}
void OffSurfaceModelPicture::normalize() { // Fit centrally in the cube
                                           // [-renderSize/2,renderSize/2]^3
  arma::mat pointsm = arma::min(points, 0);
  arma::mat pointsM = arma::max(points, 0);
  float scale = arma::mat(pointsM - pointsm).max();
  assert(scale > 0);
  points = points - arma::repmat(0.5 * (pointsm + pointsM), points.n_rows, 1);
  points *= renderSize / scale;
}

void OffSurfaceModelPicture::random_rotation(RNG &rng) {
  arma::mat L, Q, R;
  L.set_size(3, 3);
  for (int i = 0; i < 3; i++)
    for (int j = 0; j < 3; j++)
      L(i, j) = rng.uniform();
  arma::qr(Q, R, L);
  points = points * Q;
}

void OffSurfaceModelPicture::jiggle(RNG &rng, float alpha) {
  for (int i = 0; i < 3; i++)
    points.col(i) += renderSize * rng.uniform(-alpha, alpha);
}

void OffSurfaceModelPicture::affineTransform(RNG &rng, float alpha) {
  arma::mat L = arma::eye<arma::mat>(3, 3);
  for (int i = 0; i < 3; i++)
    for (int j = 0; j < 3; j++)
      L(i, j) += rng.uniform(-alpha, alpha);
  points = points * L;
}

void OffSurfaceModelPicture::codifyInputData(SparseGrid &grid,
                                             std::vector<float> &features,
                                             int &nSpatialSites,
                                             int spatialSize) {
  features.push_back(0); // Background feature
  grid.backgroundCol = nSpatialSites++;
  for (int i = 0; i < surfaces.size(); ++i) {
    // assume triangles
    drawTriangleOFF(grid, spatialSize, features, nSpatialSites,
                    points(surfaces[i][0], 0), points(surfaces[i][0], 1),
                    points(surfaces[i][0], 2),
                    points(surfaces[i][1], 0) - points(surfaces[i][0], 0),
                    points(surfaces[i][1], 1) - points(surfaces[i][0], 1),
                    points(surfaces[i][1], 2) - points(surfaces[i][0], 2),
                    points(surfaces[i][2], 0) - points(surfaces[i][0], 0),
                    points(surfaces[i][2], 1) - points(surfaces[i][0], 1),
                    points(surfaces[i][2], 2) - points(surfaces[i][0], 2));
  }
}

// Map tetrahedron to simplex
arma::mat convPyrCub =
    "1 0 0;  -0.57735   1.1547  0;  -0.40825  -0.40825 1.2247";
// (0,0,0) -> (0.0.0)
// (1,0,0) -> (1,0,0)
// (0.5,0.866,0) -> (0,1,0)
// (0.5,0.289,0.82) -> (0,0,10

Picture *OffSurfaceModelPicture::distort(RNG &rng, batchType type) {
  OffSurfaceModelPicture *pic = new OffSurfaceModelPicture(*this);
  pic->random_rotation(rng);
  pic->normalize();
  if (type == TRAINBATCH) {
    pic->affineTransform(rng, 0.2);
    pic->jiggle(rng, 0.2);
  }
  pic->points = pic->points * convPyrCub;
  return pic;
}