/*
 * Cylinder.C - methods for cylinder manipulations.
 *
 * Copyright (C) 1992, Christoph Streit (streit@iam.unibe.ch)
 *                     University of Berne, Switzerland
 * Copyright (C) 1989, 1991, Craig E. Kolb
 * 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 "Cylinder.h"
#include "transform.h"

//___________________________________________________________ Cylinder

Cylinder::Cylinder(real radius, const Vector& bottom, const Vector& top)
{
  Vector axis = top - bottom;
  real len = axis.normalize();

  /*
   * matrix trans transforms unit cylinder (1, [0,0,0], [0,0,1]) to
   * the given cylinder (radius, bottom, top).
   */
  trans  = new TransMatrix(coordSystemTransform(bottom, axis, 
						radius, len));
  itrans = new TransMatrix(*trans);
  itrans->invert();
}

GeoObject* Cylinder::create(real radius, const Vector& bottom, const Vector& top)
{
  if (radius < EPSILON)
    return NULL;

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

  return new Cylinder(radius, bottom, top);
}

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

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

  intersectionTests++;

  a = ray.getDir()[0]*ray.getDir()[0] + ray.getDir()[1]*ray.getDir()[1];
  if (a < EPSILON*EPSILON)
    /* |ray.getDir()[2]| == 1. */
    return FALSE;
  
  b = ray.getDir()[0]*ray.getOrig()[0]
    + ray.getDir()[1]*ray.getOrig()[1];
  c = ray.getOrig()[0]*ray.getOrig()[0] 
    + ray.getOrig()[1]*ray.getOrig()[1] - 1;
  disc = b*b - a*c;

  if(disc < 0.)
    return FALSE;

  disc = sqrt(disc);
  t1 = (-b + disc) / a;
  t2 = (-b - disc) / a;
  if (t1 < minDist && t2 < minDist)
    return FALSE;

  zpos1 = ray.getOrig()[2] + t1 * ray.getDir()[2];
  zpos2 = ray.getOrig()[2] + t2 * ray.getDir()[2];
  
  if (t1 < minDist || zpos1 < 0. || zpos1 > 1.) {
    if (t2 < minDist || zpos2 < 0. || zpos2 > 1.)
      return FALSE;
    else
      t1 = t2;
  } 
  else if (t2 >= minDist && zpos2 >= 0. && zpos2 <= 1.)
      t1 = (t1 < t2) ? t1 : t2;
  
  if (t1 < maxDist) {
    maxDist = t1;
    intersections++;
    return TRUE;
  }
  return FALSE;
}

/*
 * Compute cylinder normal at the given point.
 */

Vector Cylinder::normal(const Vector& p) const
{
  /*
   * The normal is equal to the point of intersection in cylinder
   * space, but with Z = 0.;
   */
  return Vector(p[0], p[1], 0);
}

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

PolygonList* Cylinder::tesselate(const BoundingBox&)
{
  const int resolution = 10;
  real delta = 2*M_PI/(real)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, lastY, 0), Vector(x, y, 0),
		    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 cylinder, postmultiply
 * it to the already computed transformations.
 */

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

  return GeoObject::setTransform(matMul);
}