object.C

/*
 * This file is part of the "Archon" framework.
 * (http://files3d.sourceforge.net)
 *
 * Copyright © 2002 by Kristian Spangsege and Brian Kristiansen.
 *
 * Permission to use, copy, modify, and distribute this software and
 * its documentation under the terms of the GNU General Public License is
 * hereby granted. No representations are made about the suitability of
 * this software for any purpose. It is provided "as is" without express
 * or implied warranty. See the GNU General Public License
 * (http://www.gnu.org/copyleft/gpl.html) for more details.
 *
 * The characters in this file are ISO8859-1 encoded.
 *
 * The documentation in this file is in "Doxygen" style
 * (http://www.doxygen.org).
 */

#include <float.h>
#include <math.h>
#include <vector>
#include <float.h>

#include <archon/math/geometry.H>
#include <archon/math/quartic_solve.H>
#include <archon/util/miniball.H>

#include "object.H"

#if ! defined M_PI
#define M_E           2.7182818284590452353602874713526625  /* e */
#define M_LOG2E       1.4426950408889634073599246810018922  /* log_2 e */
#define M_LOG10E      0.4342944819032518276511289189166051  /* log_10 e */
#define M_LN2         0.6931471805599453094172321214581766  /* log_e 2 */
#define M_LN10        2.3025850929940456840179914546843642  /* log_e 10 */
#define M_PI          3.1415926535897932384626433832795029  /* pi */
#define M_PI_2        1.5707963267948966192313216916397514  /* pi/2 */
#define M_PI_4        0.7853981633974483096156608458198757  /* pi/4 */
#define M_1_PI        0.3183098861837906715377675267450287  /* 1/pi */
#define M_2_PI        0.6366197723675813430755350534900574  /* 2/pi */
#define M_2_SQRTPI    1.1283791670955125738961589031215452  /* 2/sqrt(pi) */
#define M_SQRT2       1.4142135623730950488016887242096981  /* sqrt(2) */
#define M_SQRT1_2     0.7071067811865475244008443621048490  /* 1/sqrt(2) */
#endif

namespace Archon
{
  namespace Raytracer
  {
    using namespace Utilities;
    using namespace Math;

    inline bool leastPositiveDist(double t1, double t2, double &dist)
    {
      if(t1 > t2)
      {
        if(t1 >= 0)
        {
          dist = t2 >= 0 ? t2 : t1;
          return true;
        }
      }
      else
      {
        if(t2 >= 0)
        {
          dist = t1 >= 0 ? t1 : t2;
          return true;
        }
      }
      return false;
    }

    bool Plane::intersect(const Line3 &ray, const Object *originObject,
                          const Object::Part *, double &dist,
                          const Object::Part **p) const
    {
      /*
       * A ray that originates from a plane and leaves
       * it can never intersect it again.
       */
      if(originObject == this) return false;

      const double cosine = dot(normal, ray.direction);
      if(cosine > 0) return false; // backface culling!!
      dist = dot(normal, origin - ray.point) / cosine;
      if(p) *p = this;
      return true;
    }

    void Plane::map(const Vector3 &point, Vector3 &normal) const
    {
      normal = this->normal;
    }

    void Plane::map(const Vector3 &point, Vector3 &normal,
                    Vector2 &texturePoint) const
    {
      Vector3 p = invMap(point);
      texturePoint[0] = 0.5 + p[0]/2;
      texturePoint[1] = 0.5 - p[2]/2;
      normal = this->normal;
    }

    Cylinder::Cylinder(const Surface *surface, const CoordSystem3x3 &s):
      Object::Part(surface)
    {
      invMap.setInverseOf(s);
      pointToNormalMap = invMap;
      pointToNormalMap.basis.transpose();
      pointToNormalMap.basis.map(pointToNormalMap.origin);
      pointToNormalMap.basis *= Matrix3x3(Vector3(1, 0, 0), Vector3(0, 0, 0), Vector3(0, 0, 1)) * invMap.basis;
      {
        CoordSystem3x3 s2 = s;
        s2.translate(Vector3(0, 1, 0));
        topCap = new Plane(surface, s2);
      }
      {
        CoordSystem3x3 s2 = s;
        s2.basis.scale(Vector3(1, -1, -1));
        s2.translate(Vector3(0, 1, 0));
        bottomCap = new Plane(surface, s2);
      }
    }

    Cylinder::~Cylinder()
    {
      delete topCap;
      delete bottomCap;
    }

    bool Cylinder::intersect(const Line3 &ray, const Object *originObject,
                             const Object::Part *, double &dist,
                             const Object::Part **p) const
    {
      /*
       * A cylinder is convex so:
       * A ray that originates from the surface of a cylinder
       * and leaves the cylinder can never intersect its surface again.
       */
      if(originObject == this) return false;

      const Vector3 rp = invMap(ray.point);
      const Vector3 rd = invMap.basis(ray.direction);

      const double a = 2 * (rd[0] * rd[0] + rd[2] * rd[2]);
      const double b = -2 * (rd[0] * rp[0] + rd[2] * rp[2]);
      const double c = (rp[0] * rp[0] + rp[2] * rp[2]) - 1;
      const double d = b * b - 2 * a * c;

      if(d < 0) return false; // No solutions

      const double sqrt_d = sqrt(d);
      const double t = (b - sqrt_d) / a;
      // Note that division by zero is fine
      // Note that if b < 0 then t < 0

      const double y = rp[1] + rd[1] * t;

      if(y > 1)
      {
        // Intersection with top cap
        const double t2 = (b + sqrt_d) / a;
        if(t2 < 0) return false;
        const double y2 = rp[1] + t2 * rd[1];
        if(y2 > 1) return false;
        dist = (1 - rp[1]) / rd[1];
        if(p) *p = topCap;
        return true;
      }

      if(y < -1)
      {
        // Intersection with bottom cap
        const double t2 = (b + sqrt_d) / a;
        if(t2 < 0) return false;
        const double y2 = rp[1] + t2 * rd[1];
        if(y2 < -1) return false;
        dist = (-1 - rp[1]) / rd[1];
        if(p) *p = bottomCap;
        return true;
      }

      if(b <= 0) return false;

      dist = t;
      if(p) *p = this;
      return true;
    }

    void Cylinder::interiorIntersect(const Line3 &ray,
                                     const Object::Part *,
                                     double &dist,
                                     const Object::Part *&p) const
    {
      const Vector3 rp = invMap(ray.point);
      const Vector3 rd = invMap.basis(ray.direction);

      const double a = 2 * (rd[0] * rd[0] + rd[2] * rd[2]);
      const double b = -2 * (rd[0] * rp[0] + rd[2] * rp[2]);
      const double c = (rp[0] * rp[0] + rp[2] * rp[2]) - 1;
      const double d = b * b - 2 * a * c;

      if(d < 0)
      {
        dist = 0;
        p = this;
        return;
      }

      const double sqrt_d = sqrt(d);

      // Because we are inside the cylinder we may assume that there is
      // one positive and one negative distance, we only want the
      // positive one which of cause is the highest one.
      // Note that division by zero is fine
      // Note that if b < 0 then t < 0
      const double t = (b + sqrt_d) / a;

      const double y = rp[1] + rd[1] * t;

      if(y > 1)
      {
        // Intersection with top cap
        dist = (1 - rp[1]) / rd[1];
        p = topCap;
        return;
      }

      if(y < -1)
      {
        // Intersection with bottom cap
        dist = (-1 - rp[1]) / rd[1];
        p = bottomCap;
        return;
      }

      dist = t;
      p = this;
    }

    void Cylinder::map(const Vector3 &point, Vector3 &normal,
                       Vector2 &texturePoint) const
    {
      map(point, normal);

      Vector3 p = invMap(point);
    
      texturePoint[0] = p[0] == 0 ? p[2] < 0 ? 0 : 0.5 :
        (p[0] < 0 ? 0.25 : 0.75) - atan(p[2]/p[0]) / (2 * M_PI);

      texturePoint[1] = 0.5 + p[1]/2;
    }

    bool Sphere::intersect(const Line3 &ray, const Object *originObject,
                           const Object::Part *, double &dist,
                           const Object::Part **p) const
    {
      /*
       * A ray that originates from the surface of a sphere
       * and leaves the sphere can never intersect its surface again.
       */
      if(originObject == this) return false;

      const Vector3 rp = invMap(ray.point);
      const Vector3 rd = invMap.basis(ray.direction);
      const double b = -2 * dot(rd, rp);
      if(b <= 0) return false;
      const double a = 2 * rd.squareSum();
      const double c = rp.squareSum() - 1;
      const double d = b * b - 2 * a * c;
      if(d < 0) return false; // No solutions
      dist = (b - sqrt(d)) / a;
      if(p) *p = this;
      return true;
    }

    void Sphere::interiorIntersect(const Line3 &ray,
                                   const Object::Part *,
                                   double &dist,
                                   const Object::Part *&p) const
    {
      const Vector3 rp = invMap(ray.point);
      const Vector3 rd = invMap.basis(ray.direction);
      const double b = -2 * dot(rd, rp);
      const double a = 2 * rd.squareSum();
      const double c = rp.squareSum() - 1;
      const double d = b * b - 2 * a * c;
      dist = d < 0 ? 0 : (b + sqrt(d)) / a;
      p = this;
    }

    void Sphere::map(const Vector3 &point, Vector3 &normal,
                     Vector2 &texturePoint) const
    {
      map(point, normal);

      Vector3 p = invMap(point);
    
      /*
       * Determine the longitude angle of the point as a value between 0
       * and 1 (starting and ending at the prime meridian and rising to
       * the east.)
       * Note that the longitude is actually the shortest angle between the
       * prime meridian and the meridian of the point.
       */
      texturePoint[0] = p[0] == 0 ? p[2] < 0 ? 0 : 0.5 :
        (p[0] < 0 ? 0.25 : 0.75) - atan(p[2]/p[0]) / (2 * M_PI);

      /*
       * Determine the latitude angle of the point as a value between 0
       * and 1 (0 beeing at the south pole and 1 beeing at the north
       * pole.)
       * Note that the latitude is actually the shortest angle between the
       * intersection direction and the equator plane.
       */
      texturePoint[1] = acos(-p[1]/p.length()) / M_PI;
    }


    bool Torus::intersect(const Line3 &ray, const Object *originObject,
                          const Object::Part *, double &dist,
                          const Object::Part **part) const
    {
      /*
       * Map the ray into the coordinate system of the torus.
       */
      Vector3 p = ray.point;
      invMap.map(p);
      Vector3 d = ray.direction;
      invMap.basis.map(d);
      const double l = d.length();
      d /= l;

      const double spx = p[0]*p[0];
      const double spy = p[1]*p[1];
      const double spz = p[2]*p[2];

      const double dpx = d[0]*p[0];
      const double dpy = d[1]*p[1];
      const double dpz = d[2]*p[2];

      const double dp  = dpx + dpy + dpz;

      const double sma = majorRadius*majorRadius;

      const double a = spx + spy + spz - minorRadius*minorRadius;
      const double b = a + sma;

      const double k3 = 4*dp;
      const double k2 = 2*(2*dp*dp + a - sma*(1 - 2*d[1]*d[1]));
      const double k1 = 4*(dp*a - sma*(dp - 2*dpy));
      const double k0 = b*b - 4*sma*(spx + spz);

      // Solve: u^4 + k3 u^3 + k2 u^2 + k1 u + k0 = 0

      double roots[4];
      const int numberOfRealRoots =
        Math::quarticSolve(k3, k2, k1, k0, roots);

      if(numberOfRealRoots == 0) return false;

      if(originObject == this)
      {
        int minAbs = -1;
        for(int i=0; i<numberOfRealRoots; ++i)
          if(minAbs == -1 || fabs(roots[i]) < fabs(roots[minAbs])) minAbs = i;
        roots[minAbs] = -1;
      }

      double minDist = 0;
      bool any = false;
      for(int i=0; i<numberOfRealRoots; ++i)
      {
        if(roots[i] > 0 && (!any || roots[i] < minDist))
        {
          minDist = roots[i];
          any = true;
        }
      }
      if(!any) return false;
      dist = minDist / l;
      if(part) *part = this;
      return true;
    }

    void Torus::interiorIntersect(const Line3 &ray,
                                  const Object::Part *,
                                  double &dist,
                                  const Object::Part *&p) const
    {
      if(!intersect(ray, this, this, dist, &p)) p = this, dist = 0;
    }

    void Torus::map(const Vector3 &point, Vector3 &normal) const
    {
      normal = point;
      invMap.map(normal);
      const double a = 1-majorRadius/sqrt(normal[0]*normal[0] + normal[2]*normal[2]);
      normal[0] *= a;
      normal[2] *= a;
      normalMap.map(normal);
      normal.normalize();
    }

    void Torus::map(const Vector3 &point, Vector3 &normal,
                    Vector2 &texturePoint) const
    {
      normal = point;
      invMap.map(normal);
      const double a = sqrt(normal[0]*normal[0] + normal[2]*normal[2]);
      const double b = a - majorRadius;

      texturePoint[0] = normal[0] == 0 ? normal[2] < 0 ? 0 : 0.5 :
        (normal[0] < 0 ? 0.25 : 0.75) - atan(normal[2]/normal[0]) / (2 * M_PI);

      texturePoint[1] = normal[1] == 0 ? b < 0 ? 0 : 0.5 :
        (normal[1] < 0 ? 0.25 : 0.75) - atan(b/normal[1]) / (2 * M_PI);

      const double c = b/a;
      normal[0] *= c;
      normal[2] *= c;
      normalMap.map(normal);
      normal.normalize();
    }


    Polygon::Polygon(const Surface *surface,
                     const std::vector<Vector3 *> &vertexList,
                     const Vector2 &textureVertex1,
                     const Vector2 &textureVertex2,
                     const Vector2 &textureVertex3):
      Object::Part(surface), vertexList(vertexList)
    {
      if(vertexList.size() < 3) ARCHON_THROW1(ArgumentException, "Too few vertices");
      Vector3 x = *vertexList[2];
      x -= *vertexList[1];
      Vector3 y = *vertexList[0];
      y -= *vertexList[1];
      normal = x;
      normal *= y;
      normal.normalize();

      if(!dynamic_cast<const Texture *>(surface)) return;

      // Polygon to texture map
      Matrix3x3 m;
      m.x = Vector3(textureVertex3);
      m.x -= Vector3(textureVertex2);
      m.y = Vector3(textureVertex1);
      m.y -= Vector3(textureVertex2);
      m.z = m.x;
      m.z *= m.y;

      // World to texture map
      m /= Matrix3x3(x, y, normal);

      texMapX = m.getRow1();
      texMapY = m.getRow2();
      texMapOrigin = textureVertex2;
      texMapOrigin -= Vector2(m(*vertexList[1]));
    }

    Polygon::Polygon(const Surface *surface,
                     const std::vector<Vector3 *> &vertexList,
                     const Vector2 &textureVertex1,
                     const Vector2 &textureVertex2):
      Object::Part(surface), vertexList(vertexList)
    {
      if(vertexList.size() < 3) ARCHON_THROW1(ArgumentException, "Too few vertices");
      normal = *vertexList[1];
      normal -= *vertexList[0];
      Vector3 v = *vertexList[2];
      v -= *vertexList[1];
      normal *= v;
      normal.normalize();
      /*
      // Texture
      Vector3 wX = *vertexList[1];
      wX -= *vertexList[0];

      Vector2 tX = textureVertex2;
      tX -= textureVertex1;

      scaling = tX.length() / wX.length();

      wX.normalize();

      Vector3 wY = normal;
      wY *= wX;

      tX.normalize();

      Vector2 tY = tX;
      tY.cross();

      Vector2 w;

      w = -textureVertex1;
      textureOrigin = (wX * dot(w, tX) + wY * dot(w, tY))/scaling +
      *vertexList[0];

      w = Vector2::unitX - textureVertex1;
      textureAxisX = (wX * dot(w, tX) + wY * dot(w, tY))/scaling +
      *vertexList[0] - textureOrigin;

      w = Vector2::unitY - textureVertex1;
      textureAxisY = (wX * dot(w, tX) + wY * dot(w, tY))/scaling +
      *vertexList[0] - textureOrigin;

      textureAxisX.normalize();
      textureAxisY.normalize();
      */
    }

    bool Polygon::intersect(const Line3 &ray, const Object *originObject,
                            const Object::Part *, double &dist,
                            const Object::Part **part) const
    {
      /*
       * A ray that originates from the surface of a polygon and leaves
       * the polygon can never intersect its polygon again.
       */
      if(originObject == this) return false;

      const double cosine = dot(normal, ray.direction);
      if(cosine > 0) return false; // backface culling!!
      Vector3 v, w = *vertexList.back();
      for(std::vector<Vector3 *>::const_iterator i=vertexList.begin();
          i != vertexList.end(); ++i)
      {
        v = **i;
        w -= v;
        v -= ray.point;
        w *= v;
        /*
         * 'w' is now the normal to the plain containing both the eye
         * point and the polygon edge between the current and the
         * previous vertex.  The direction of the normal from the plane
         * is away from the polygon interior.
         */
        if(dot(w, ray.direction) > 0) return false;
        w = **i;
      }

      if(part) *part = this;
      dist = dot(normal, v) / cosine;
      return true;
    }

    bool Polygon::backfaceIntersect(const Line3 &ray,
                                    const Object *originObject,
                                    const Object::Part *, double &dist,
                                    const Object::Part **part) const
    {
      /*
       * A ray that originates from the surface of a polygon and leaves
       * the polygon can never intersect its polygon again.
       */
      if(originObject == this) return false;

      const double cosine = dot(normal, ray.direction);
      if(cosine < 0) return false; // frontface culling!!

      Vector3 v, w = *vertexList.front();
      for(std::vector<Vector3 *>::const_iterator i=vertexList.end();
          i != vertexList.begin(); --i)
      {
        v = **(i-1);
        w -= v;
        v -= ray.point;
        w *= v;
        /*
         * 'w' is now the normal to the plain containing both the eye
         * point and the polygon edge between the current and the
         * previous vertex.  The direction of the normal from the plane
         * is away from the polygon interior.
         */
        if(dot(w, ray.direction) > 0) return false;
        w = **(i-1);
      }

      if(part) *part = this;
      dist = dot(normal, v) / cosine;
      return true;
    }

    void Polygon::map(const Vector3 &point, Vector3 &normal,
                      Vector2 &texturePoint) const
    {
      normal = this->normal;

      /*
        Vector3 p = point;
        p -= textureOrigin;
        texturePoint[0] = dot(p, textureAxisX) * scaling;
        texturePoint[1] = dot(p, textureAxisY) * scaling;
      */

      texturePoint[0] = dot(point, texMapX);
      texturePoint[1] = dot(point, texMapY);
      texturePoint += texMapOrigin;
    }

    ConvexPolyhedron::
    ConvexPolyhedron(const std::vector<Polygon *> &polygonList):
      polygonList(polygonList) 
    {
      // Determine the minimal bounding sphere
      MinimalBoundingBall<3> ball;
      for(std::vector<Polygon *>::const_iterator i=polygonList.begin();
          i != polygonList.end(); ++i)
      {
        for(std::vector<Vector3 *>::iterator j=(*i)->vertexList.begin();
            j != (*i)->vertexList.end(); ++j)
        {
          const Vector3 *v = *j;
          MinimalBoundingBall<3>::Point p;
          p[0]=(*v)[0];
          p[1]=(*v)[1];
          p[2]=(*v)[2];
          ball.check_in(p);
        }
      }
      ball.build();
      const MinimalBoundingBall<3>::Point c = ball.center();
      boundingVolume.center.set(c[0], c[1], c[2]);
      boundingVolume.radius = sqrt(ball.squared_radius());
    }

    ConvexPolyhedron::~ConvexPolyhedron()
    {
      for(std::vector<Polygon *>::iterator i=polygonList.begin();
          i != polygonList.end(); ++i) delete *i;
    }

    bool ConvexPolyhedron::intersect(const Line3 &ray,
                                     const Object *originObject,
                                     const Part *originPart, double &dist,
                                     const Object::Part **part) const
    {
      /*
       * A ray that originates from the surface of a convex polyhedron
       * and leaves the polyhedron can never intersect its surface
       * again. 
       */
      if(originObject == this) return false;

      double d1 = 0, d2 = 0;
      if(!Math::intersect(Ray3(ray), boundingVolume, d1))
        return false;

      bool hit = false;
      const Part *p1, *p2 = 0;

      for(std::vector<Polygon *>::const_iterator i=polygonList.begin();
          i != polygonList.end(); ++i)
        if((*i)->intersect(ray, originObject, originPart, d1, &p1) &&
           (!hit || d1 < d2))
        {
          d2 = d1;
          p2 = p1;
          hit = true;
        }

      if(hit)
      {
        dist = d2;
        if(part) *part = p2;
        return true;
      }

      return false;
    }

    void ConvexPolyhedron::interiorIntersect(const Line3 &ray,
                                             const Object::Part *originPart,
                                             double &dist,
                                             const Object::Part *&p) const
    {
      for(std::vector<Polygon *>::const_iterator i=polygonList.begin();
          i != polygonList.end(); ++i)
      {
        if(originPart == *i) continue;
        if((*i)->backfaceIntersect(ray, this, originPart, dist, &p)) return;
      }

      // In case of numeric inaccuracy.
      p = originPart;
      dist = 0;
    }
  }
}

---