Photo CD Demo 1/Demo.bin/graphtal/cone.c

C/C++ Source or Header | 1992-10-22 | 5KB | 228 lines

/*
 * Cone.C - methods for cone manipulations.
 *
 * Copyright (C) 1992, Christoph Streit ([e-mail])
 * All rights reserved.
 *
 * This software may be freely copied, modified, and redistributed
 * provided that this copyright notice is preserved on all copies.
 *
 * You may not distribute this software, in whole or in part, as part of
 * any commercial product without the express consent of the authors.
 *
 * There is no warranty or other guarantee of fitness of this software
 * for any purpose.  It is provided solely "as is".
 *
 */

#include "Cone.h"
#include "Cylinder.h"
#include "transform.h"

//___________________________________________________________ Cone

Cone::Cone(real r1, const Vector& bottom,
       real r2, const Vector& top)
{
  /*
   * Find the axis and axis length.
   */
  Vector axis = top - bottom;
  real len = axis.normalize();

  /*
   * To render a cone, we transform the desired cone into
   * a canonical, Z-axis aligned, unit length, unit radius
   * at the apex cone.
   */

  /*
   * "tantheta" is the change in radius per unit length along
   * the cone axis.
   */
  real tantheta = (r2-r1)/len;

  /*
   * lprime defines the distance along the axis where the cone starts
   */
  real lprime = r1/tantheta;

  /*
   * Find the true base (origin) of the cone.
   */
  Vector base = (-lprime)*axis + bottom;

  /*
   * tlen is the total length of the cone.
   */
  real totalLengthOfCone = lprime+len;

  startDistance = lprime/totalLengthOfCone;

  trans = new TransMatrix(coordSystemTransform(base, axis,
                     r2, totalLengthOfCone));
  itrans = new TransMatrix(*trans);
  itrans->invert();
}

GeoObject* Cone::create(real r1, const Vector& bottom,
            real r2, const Vector& top)
{
  /*
   * Check the arguments.
   */
  if (equal(r1, r2))
    return Cylinder::create(r1, bottom, top);

  if ((top-bottom).length() < EPSILON)
    return NULL;

  if ((r1 < 0) || r2 < 0)
    return NULL;

  /*
   * The passed basepoint must be closer to the origin of the
   * cone than the apex point, implying that the base radius
   * must be smaller than the apex radius.  If the values passed
   * reflect the opposite, we switch everything.
   */
  if (r2 < r1)
    return new Cone(r2, top, r1, bottom);
  else
    return new Cone(r1, bottom, r2, top);
}

/*
 * Ray cone intersection.
 *
 * The code is adapted from Craig Kolbs rayshade.
 */

int Cone::intersect(const Ray& ray, real minDist, real& maxDist)
{
  real t1, t2, a, b, c, disc, zpos;

  intersectionTests++;

  a = ray.getDir()[0]*ray.getDir()[0] + ray.getDir()[1]*ray.getDir()[1]
    - ray.getDir()[2]*ray.getDir()[2];
  b = ray.getDir()[0]*ray.getOrig()[0] + ray.getDir()[1]*ray.getOrig()[1]
    - ray.getDir()[2]*ray.getOrig()[2];
  c = ray.getOrig()[0]*ray.getOrig()[0] + ray.getOrig()[1]*ray.getOrig()[1]
    - ray.getOrig()[2]*ray.getOrig()[2];

  if (fabs(a) < EPSILON) {
    /*
     * Only one intersection point...
     */
    t1 = -0.5*c / b;
    zpos = ray.getOrig()[2] + t1*ray.getDir()[2];
    if (t1 < minDist || zpos < startDistance || zpos > 1.)
      return FALSE;

    if (t1 < maxDist) {
      maxDist = t1;
      intersections++;
      return TRUE;
    }
    return FALSE;
  }
  else {
    disc = b*b - a*c;
    if (disc < 0.)
      return FALSE;        /* No possible intersection */

    disc = sqrt(disc);
    t1 = (-b + disc) / a;
    t2 = (-b - disc) / a;

    /*
     * Clip intersection points.
     */
    zpos = ray.getOrig()[2] + t1*ray.getDir()[2];
    if (t1 < minDist || zpos < startDistance || zpos > 1.) {
      zpos = ray.getOrig()[2] + t2*ray.getDir()[2];
      if (t2 < minDist || zpos < startDistance || zpos > 1.)
    return FALSE;
      else
    t1 = t2;
    }
    else {
      zpos = ray.getOrig()[2] + t2*ray.getDir()[2];
      if (t2 >= minDist && zpos >= startDistance && zpos <= 1. && t2 < t1)
    t1 = t2;
    }
    if (t1 < maxDist) {
      maxDist = t1;
      intersections++;
      return TRUE;
    }
    return FALSE;
  }
}

/*
 * Compute cone normal at point p.
 */

Vector Cone::normal(const Vector& p) const
{
  /*
   * If p is on the z-axis, return the cone axis (0,0,1)
   */
  if (equal(p[0], 0) && equal(p[1], 0))
    return Vector(0,0,1);

  /*
   * The following is equal to (p*(0, 0, 1))*p
   */
  return Vector(p[0]*p[2], p[1]*p[2], -p[0]*p[0] - p[1]*p[1]);
}


/*
 * Split the cone into polygon (rectangles)
 */

PolygonList* Cone::tesselate(const BoundingBox&)
{
  const real resolution = 10;
  real delta = 2*M_PI/resolution;
  real lastX = 1;
  real lastY = 0;
  real x, y;

  Polygon* p;
  PolygonList* polys = new PolygonList(resolution);

  /*
   * Generate points on the edge of the disc.
   */
  for (real alpha=delta; alpha <= 2*M_PI; alpha += delta) {
    x = cos(alpha); y = sin(alpha);
    p = new Polygon(Vector(lastX*startDistance,
               lastY*startDistance, startDistance),
            Vector(x*startDistance, y*startDistance, startDistance),
            Vector(x, y, 1), Vector(lastX, lastY, 1));
    p->transform(*trans);
    polys->append(p);
    lastX = x; lastY = y;
  }

  return polys;
}

/*
 * If a transformation is given for the defined cone, postmultiply
 * it to the already computed transformations.
 */

int Cone::setTransform(TransMatrix* tmat)
{
  TransMatrix* matMul = new TransMatrix((*trans) * (*tmat));
  delete itrans;
  delete trans;

  return GeoObject::setTransform(matMul);
}