Off3DFormatPicture.cpp

#include "Off3DFormatPicture.h"
#include <string>
#include <fstream>

int mapToGridOFF(float coord, int inputFieldSize) {
  return std::max(
      0, std::min(inputFieldSize - 1, (int)(coord + 0.5 * inputFieldSize)));
}
// 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 n =
          mapToGridOFF(p0, inputFieldSize) * inputFieldSize * inputFieldSize +
          mapToGridOFF(p1, inputFieldSize) * inputFieldSize +
          mapToGridOFF(p2, inputFieldSize);
      if (grid.mp.find(n) == grid.mp.end()) {
        grid.mp[n] = nSpatialSites++;
        features.push_back(1);
      }
    }
  }
}

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));
  }
}

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);
  }
  return pic;
}