plane_uv_calculation_using_intersection_point.c

plane_uv_calculation_using_intersection_point.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "../raytracing/helpers.h"
#include "../raytracing/vector.h"
#include "../raytracing/matrix.h"

int main(void) {

  #define QUAD 0
  #define RECT 0
  #define TRIA 1
  #if QUAD
  Vector p00 = Vector_New(1,1,0);
  Vector p10 = Vector_New(5,2,0);
  Vector p11 = Vector_New(5,4,0);
  Vector p01 = Vector_New(2,4,0);
  #elif RECT
  Vector p00 = Vector_New(1,1,0);
  Vector p10 = Vector_New(5,1,0);
  Vector p11 = Vector_New(5,4,0);
  Vector p01 = Vector_New(1,4,0);
  #elif TRIA
  Vector p00 = Vector_New(0,0,0);
  Vector p10 = Vector_New(10,0,0);
  Vector p11 = Vector_New(5,5,0);
  Vector p01 = p11;
  #endif

  // F(u,v) = p00 + u(p10-p00) + v(p01-p00) + uv(p11-p01 + p00-p10);

  Vector tmp1 = Vector_SubVector(p11,p01);
  Vector tmp2 = Vector_SubVector(p00,p10);
  Vector pa = Vector_AddVector(tmp1,tmp2);
  Vector pb = Vector_SubVector(p10,p00);
  Vector pc = Vector_SubVector(p01,p00);
  Vector pd = p00;

  // F(u,v) = pd + u(pb) + v(pc) + uv(pa);

  Vector pn = Vector_Cross(pb,pc);

  // calc a normal in the F(u,0) to F(u,1) line
  // using this line as vector that lies in perpendicular
  // plane with respect the plane polygon (pn)
  // N_u = pn X (F(u,1) - F(u,0))
  // N_u = pn X ((pd + u(pb) + pc + u(pa)) - (pd + u(pb)))
  // N_u = pn X (pd + u(pb) + pc + u(pa) - pd - u(pb))
  // N_u = pn X (pc + u(pa))
  // N_u = pn X pc + u(pn X pa)
  Vector nc = Vector_Cross(pn,pc);
  Vector na = Vector_Cross(pn,pa);
  // N_u = nc + u(na)

  Vector ri = Vector_DivScalar(
    Vector_AddVector(p00,Vector_AddVector(p10,Vector_AddVector(p11,p01))),
    4
  );
  // equation of the plane with u is (ri is the intersection point)
  // N_u(F(u,0) - ri) = 0
  // (nc + u(na))(p00 + u(p10-p00) - ri) = 0
  // (nc + u(na))(pd + u(pb) - ri) = 0
  // (nc + u(na))pd + (nc + u(na))u(pb) - (nc + u(na))ri = 0
  // * is the dot product
  // nc*pd + u(na*pd) + u(nc*pb) + u^2(na*pb) - nc*ri - u(na*ri) = 0
  // u^2(na*pb) + u(na*pd + nc*pb - na*ri) + (nc*pd - nc*ri) = 0
  // this is quadratic in u
  // A = na*pb
  // B = na*pd + nc*pb - na*ri
  // C = nc*pd - nc*ri

  double A = Vector_Dot(na,pb);
  double B = Vector_Dot(na,pd) + Vector_Dot(nc,pb) - Vector_Dot(na,ri);
  double C = Vector_Dot(nc,pd) - Vector_Dot(nc,ri);
  double u;
  if(fabs(A) < EPSILON) {
    u = -C/B;
    printf("A is zero so u is %f\n",u);

  } else {
    double disc = B*B - 4*A*C;
    printf("disc: %f\n",disc);
    double radicand = sqrt(disc);
    double deno = 2*A;
    double u1 = (-B - radicand) / deno;
    double u2 = (-B + radicand) / deno;
    printf("u1: %f, u2: %f\n",u1,u2);
    u = MIN(u1,u2);
  }

  if(u < 0 || u > 1) {
    printf("not intersection possible\n");
    return 0;
  }

  // calc normal for N_v
  // N_v = pn X (F(1,v) - F(0,v))
  // N_v = pn X ((pd + pb + v(pc) + v(pa)) - (pd + v(pc)))
  // N_v = pn X (pd + pb + v(pc) + v(pa) - pd - v(pc))
  // N_v = pn X (pb + v(pa)
  // N_v = pn X pb + v(pn X pa)
  Vector nb = Vector_Cross(pn,pb);
  //Vector na = Vector_Cross(pn,pa);
  // N_v = nb + v(na)

  // equation of the plane with respect to v is
  // N_v(F(0,v) - ri) = 0
  // (nb + v(na))(pd + v(pc) - ri) = 0
  // (nb + v(na))pd + (nb + v(nc))v(pc) - (nb + v(nc))ri = 0
  // nb*pd + v(na*pd) + v(nb*pc) + v^2(na*pc) - nb*ri - v(na*ri) = 0
  // v^2(na*pc) + v(na*pd + nb*pc - na*ri) + (nb*pd - nb*ri) = 0

  A = Vector_Dot(na,pc);
  B = Vector_Dot(na,pd) + Vector_Dot(nb,pc) - Vector_Dot(na,ri);
  C = Vector_Dot(nb,pd) - Vector_Dot(nb,ri);
  double v;
  if(fabs(A) < EPSILON) {
    v = -C/B;
    printf("A is zero so v is %f\n",v);

  } else {
    double disc = B*B - 4*A*C;
    printf("disc: %f\n",disc);
    double radicand = sqrt(disc);
    double deno = 2*A;
    double v1 = (-B - radicand) / deno;
    double v2 = (-B + radicand) / deno;
    printf("v1: %f, v2: %f\n",v1,v2);
    v = MIN(v1,v2);
  }
  if(v < 0 || v > 1) {
    printf("not intersection possible\n");
    return 0;
  }

  printf("u,v -> %f,%f\n",u,v);

  return 0;
}