rasterizeTriangles.cpp

///////////////////////////////////////////////////////////////////////////////
//
//  Name:        rasterizeTriangles.mex.cc
//  Purpose:     Draw triangles on a binary image
//  Author:      Daeyun Shin <dshin11@illinois.edu>
//  Created:     04.05.2015
//  Modified:    04.05.2015
//  Version:     0.1
//
//  This Source Code Form is subject to the terms of the Mozilla Public
//  License, v. 2.0. If a copy of the MPL was not distributed with this
//  file, You can obtain one at http://mozilla.org/MPL/2.0/.
//
///////////////////////////////////////////////////////////////////////////////

/**
 * Example usage - drawing a silhouette from isosurface
 *
 * fv = isosurface(X, Y, Z, V, 0.5);
 * xy = worldToImage(camera, fv.vertices);
 * triangles = round([xy(fv.faces(:,1),:) xy(fv.faces(:,2),:)
 *xy(fv.faces(:,3),:)]);
 *
 * sz = size(image);
 * silhouette = rasterizeTriangles(triangles, sz(1:2));
 */

#include "mex.h"
#include "mexutil.h"

struct Vertex {
  int x, y;
};

void DrawVertLine(const int x, const int y1, const int y2, int8_t* im,
                  const int* dims) {
  if (x < 0 || y2 < 0 || x >= dims[1] || y1 >= dims[0]) return;

  for (int y = std::max(y1, 0); y <= std::min(y2, dims[0] - 1); ++y) {
    im[x * dims[0] + y] = 1;
  }
}

void FillLeftFlatTriangle(const Vertex& v1, const Vertex& v2, const Vertex& v3,
                          int8_t* im, const int* dims) {
  float slope1 = (float)(v1.y - v2.y) / (float)(v1.x - v2.x);
  float slope2 = (float)(v1.y - v3.y) / (float)(v1.x - v3.x);

  float cury1 = v1.y;
  float cury2 = v1.y;

  for (int scanlineX = v1.x; scanlineX >= v2.x; scanlineX--) {
    DrawVertLine(scanlineX, cury1, cury2, im, dims);
    cury1 -= slope1;
    cury2 -= slope2;
  }
}

void FillRightFlatTriangle(const Vertex& v1, const Vertex& v2, const Vertex& v3,
                           int8_t* im, const int* dims) {
  float slope1 = (float)(v3.y - v1.y) / (float)(v3.x - v1.x);
  float slope2 = (float)(v3.y - v2.y) / (float)(v3.x - v2.x);

  float cury1 = v3.y;
  float cury2 = v3.y;

  for (int scanlineX = v3.x; scanlineX < v1.x; scanlineX++) {
    DrawVertLine(scanlineX, cury1, cury2, im, dims);
    cury1 += slope1;
    cury2 += slope2;
  }
}

void DrawTriangle(Vertex& v1, Vertex& v2, Vertex& v3, int8_t* im,
                  const int* dims) {
  if (v1.x < v2.x) std::swap(v1, v2);
  if (v1.x < v3.x) std::swap(v1, v3);
  if (v2.x < v3.x) std::swap(v2, v3);

  if (v2.x == v3.x) {
    if (v2.y > v3.y) {
      FillLeftFlatTriangle(v1, v3, v2, im, dims);
    } else {
      FillLeftFlatTriangle(v1, v2, v3, im, dims);
    }
  } else if (v1.x == v2.x) {
    if (v1.y > v2.y) {
      FillRightFlatTriangle(v2, v1, v3, im, dims);
    } else {
      FillRightFlatTriangle(v1, v2, v3, im, dims);
    }
  } else {
    Vertex v4;
    v4.x = v2.x;
    float slope = (float)(v3.y - v1.y) / (float)(v1.x - v3.x);
    v4.y = v3.y - (v2.x - v3.x) * slope;

    if (v2.y > v4.y) {
      FillLeftFlatTriangle(v1, v4, v2, im, dims);
      FillRightFlatTriangle(v4, v2, v3, im, dims);
    } else {
      FillLeftFlatTriangle(v1, v2, v4, im, dims);
      FillRightFlatTriangle(v2, v4, v3, im, dims);
    }
  }
}

void DrawTriangles(const int* triangle, int num_triangles, int8_t* im,
                   const int* dims) {
  for (int i = 0; i < num_triangles; ++i) {
    Vertex v1, v2, v3;
    v1.x = triangle[i * 6];
    v1.y = triangle[i * 6 + 1];
    v2.x = triangle[i * 6 + 2];
    v2.y = triangle[i * 6 + 3];
    v3.x = triangle[i * 6 + 4];
    v3.y = triangle[i * 6 + 5];
    DrawTriangle(v1, v2, v3, im, dims);
  }
}

void mexFunction(int nargout, mxArray* out[], int nargin, const mxArray* in[]) {
  using namespace mexutil;

  N_IN(2);
  N_OUT(1);

  M_ASSERT(mxGetNumberOfDimensions(in[0]) == 2);
  M_ASSERT(mxGetNumberOfDimensions(in[1]) == 2);

  const mwSize* dims0 = mxGetDimensions(in[0]);
  M_ASSERT(dims0[0] >= 1);
  M_ASSERT(dims0[1] == 6);

  const mwSize* dims1 = mxGetDimensions(in[1]);
  M_ASSERT(dims1[0] == 1);
  M_ASSERT(dims1[1] == 2);

  int tri_mat_size[2] = {(int)dims0[0], (int)dims0[1]};
  double* triangles = (double*)mxGetPr(in[0]);
  double* size1 = (double*)mxGetPr(in[1]);
  int imsize[2] = {(int)size1[0], (int)size1[1]};

  int* triangles_t =
      (int*)calloc(tri_mat_size[0] * tri_mat_size[1], sizeof(int));
  int count = 0;
  for (int i = 0; i < tri_mat_size[0]; ++i) {
    for (int j = 0; j < tri_mat_size[1]; ++j) {
      triangles_t[count++] = (int)(triangles[i + j * tri_mat_size[0]] + 0.5);
    }
  }

  int8_t* im = (int8_t*)calloc(imsize[0] * imsize[1], sizeof(int8_t));

  DrawTriangles(triangles_t, tri_mat_size[0], im, imsize);

  mwSize out_dims[2] = {(mwSize)imsize[0], (mwSize)imsize[1]};
  out[0] = mxCreateNumericArray(2, out_dims, mxDOUBLE_CLASS, mxREAL);
  double* out_im = (double*)mxGetPr(out[0]);

  for (int i = 0; i < imsize[1]; ++i) {
    for (int j = 0; j < imsize[0]; ++j) {
      out_im[i * imsize[0] + j] = (double)im[i * imsize[0] + j];
    }
  }
}