ラスタ→STLを目指せ #1 3点からSTLを書き出す

ここでいうSTLはStandard Triangulated Languageです。
ラスタの特定のバンドの値を高さ方向にとってSTLにしてしまおう、というもの。普通はラスタと言ってもDEMです。

STLの仕様は http://ja.wikipedia.org/wiki/Standard_Triangulated_Language 参照。

今回は、"facet"ブロックをひとつ出力してみましょう。

ふとどうすればいいのか分からないのが法線ベクトル。
三角形が仕様通り右ネジにそろっているとします。

右ネジにそろっているなら法線ベクトルは計算できます。たとえば http://zahyou.6.ql.bz/3db/gai.htm 参照。

/*
 * p = { x1, y1, z1, x2, y2, z2, x3, y3, z3};
 * で
 * (x1,y1,z1), (x2,y2,z2), (x3,y3,z3) が左回りになっているとします。
 */
void print_triangle(FILE *fp, double *p) {
  double v1x, v1y, v1z, v2x, v2y, v2z;
  double nx, ny, nz, an;
  int n;
  v1x = p[3] - p[0];
  v1y = p[4] - p[1];
  v1z = p[5] - p[2];
  v2x = p[6] - p[0];
  v2y = p[7] - p[1];
  v2z = p[8] - p[2];
  /* nx,ny,nzが法線ベクトル */
  nx = v1y*v2z - v1z*v2y;
  ny = v1z*v2x - v1x*v2z;
  nz = v1x*v2y - v1y*v2x;
  /* 正規化 (%eならいらないかも知れない) */
  an = sqrt(nx*nx+ny*ny+nz*nz);
  fprintf(fp,"facet normal %.14e %.14e %.14e\nouter loop\n", nx/an, ny/an, nz/an);
  for( n = 0; n < 3; n++ ) {
    fprintf(fp,"  vertex %.14e %.14e %.14e\n", p[0],p[1],p[2]);
    p += 3;
  }
  fprintf(fp,"endloop\nendfacet\n");
}