---

gunu_revolve.cpp

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */


/************************************************************************
  CIRCLE, SURFACES OF REVOLUTION, SPHERE
*************************************************************************/
 
#include "stdafx.h"

#include <iostream>

#include "gul_types.h"
#include "gul_vector.h"
#include "guma_lineq.h"
#include "gunu_revolve.h"

namespace gunu {

using gul::Ptr;
using gul::rtr;
using gul::point;
using gul::hpoint;
using gul::point2;
using gul::hpoint2;
using gul::ProjectToLine;
using guma::SVDIntersectLines;
using gul::weight;
using gul::euclid;
using gul::normalize;
using gul::cross_product;

/*---------------------------------------------------------------------*//**
  Create a circle as a NURBS curve (see "The NURBS book")                 */
/*------------------------------------------------------------------------*/  
template< class T >
void Circle(
            point<T> O,  /* center (origin of local coordinate system) */
            point<T> X,  /* X-axis of loc.-coord.-syst., |X| = |Y| = 1 */
            point<T> Y,  /* Y-axis of loc.-coord.-syst., perpendicular to X */
            T r,           /* radius */
            T theta_start, /* start angle (in radians)(in regard to X) */
            T theta_end,   /* end angle (in radians)(in regard to X) */
            int    *ret_n, /* => number of controlpoints - 1 */
            int *ret_p,    /* => degree */
            Ptr<T> *retU,             /* => knot vector */
            Ptr< hpoint<T> > *retPw   /* => control points */  )  
{
  Ptr< hpoint<T> > Pw;
  point<T> P0, T0, P1, P2, T2;
  T ths,the,theta,dtheta, w1, angle, dummy;
  Ptr<T> U;
  int narcs, n, index, i, j;

  ths = theta_start;
  the = theta_end;
  
  if( the < ths ) the += rtr<T>::pi() * (T)2;
  
  theta = the - ths;
  
  if( theta <= rtr<T>::pi() / (T)2 ) narcs = 1;
  else if( theta <= rtr<T>::pi() ) narcs = 2;
  else if( theta <= (T)1.5 * rtr<T>::pi()  ) narcs = 3;
  else narcs = 4;
  
  dtheta = theta / ((T)narcs);
  
  n = 2 * narcs;      /* n = number of control points - 1 */

  Pw.reserve_pool(n+1); 
  U.reserve_pool(n+4);

  w1 = rtr<T>::cos( dtheta/(T)2 );

  P0.x = O.x + r*rtr<T>::cos(ths)*X.x + r*rtr<T>::sin(ths)*Y.x;
  P0.y = O.y + r*rtr<T>::cos(ths)*X.y + r*rtr<T>::sin(ths)*Y.y;
  P0.z = O.z + r*rtr<T>::cos(ths)*X.z + r*rtr<T>::sin(ths)*Y.z;

  T0.x = -rtr<T>::sin(ths)*X.x + rtr<T>::cos(ths)*Y.x;
  T0.y = -rtr<T>::sin(ths)*X.y + rtr<T>::cos(ths)*Y.y;
  T0.z = -rtr<T>::sin(ths)*X.z + rtr<T>::cos(ths)*Y.z;

  Pw[0].x = P0.x;
  Pw[0].y = P0.y;
  Pw[0].z = P0.z;
  Pw[0].w = (T)1;
  
  index = 0;
  angle = ths;
  
  for( i = 1; i <= narcs; i++ )
  {
    angle += dtheta;

    P2.x = O.x + r*rtr<T>::cos(angle)*X.x + r*rtr<T>::sin(angle)*Y.x;
    P2.y = O.y + r*rtr<T>::cos(angle)*X.y + r*rtr<T>::sin(angle)*Y.y;
    P2.z = O.z + r*rtr<T>::cos(angle)*X.z + r*rtr<T>::sin(angle)*Y.z;    
        
    Pw[index+2].x = P2.x;
    Pw[index+2].y = P2.y;
    Pw[index+2].z = P2.z;
    Pw[index+2].w = (T)1;

    T2.x = -rtr<T>::sin(angle)*X.x + rtr<T>::cos(angle)*Y.x;
    T2.y = -rtr<T>::sin(angle)*X.y + rtr<T>::cos(angle)*Y.y;
    T2.z = -rtr<T>::sin(angle)*X.z + rtr<T>::cos(angle)*Y.z;

    SVDIntersectLines( P0, T0, P2, T2, &dummy, &dummy, &P1 );
  
    Pw[index+1].x = w1 * P1.x;
    Pw[index+1].y = w1 * P1.y;
    Pw[index+1].z = w1 * P1.z;
    Pw[index+1].w = w1;            

    index += 2;
    
    if( i < narcs )
    {
      P0 = P2;
      T0 = T2;
    }  
  }  
  j = 2 * narcs + 1;
  
  for( i = 0; i < 3; i++ )
  {
    U[i] = (T)0;
    U[i+j] = (T)1;
  }  
  switch( narcs )
  {
    case 2:
      U[3] = U[4] = (T)0.5;
      break;
    case 3:
      U[3] = U[4] = (T)1 / (T)3;
      U[5] = U[6] = (T)2 / (T)3;
      break;
    case 4:
      U[3] = U[4] = (T)0.25;  
      U[5] = U[6] = (T)0.5;  
      U[7] = U[8] = (T)0.75;  
      break;
  }    
  *ret_n = n;
  *ret_p = 2;
  *retPw = Pw;
  *retU = U;
}
template void Circle(
            point<float> O,  /* center (origin of local coordinate system) */
            point<float> X,  /* X-axis of loc.-coord.-syst., |X| = |Y| = 1 */
            point<float> Y,  /* Y-axis of loc.-coord.-syst., perpendicular to X */
            float r,           /* radius */
            float theta_start, /* start angle (in regard to X) */
            float theta_end,   /* end angle (in regard to X) */
            int    *ret_n, /* => number of controlpoints - 1 */
            int *ret_p,    /* => degree */
            Ptr<float> *retU,             /* => knot vector */
            Ptr< hpoint<float> > *retPw   /* => control points */  );
template  void Circle(
            point<double> O,  /* center (origin of local coordinate system) */
            point<double> X,  /* X-axis of loc.-coord.-syst., |X| = |Y| = 1 */
            point<double> Y,/* Y-axis of loc.-coord.-syst., perpendicular to X*/
            double r,           /* radius */
            double theta_start, /* start angle (in regard to X) */
            double theta_end,   /* end angle (in regard to X) */
            int    *ret_n, /* => number of controlpoints - 1 */
            int *ret_p,    /* => degree */
            Ptr<double> *retU,             /* => knot vector */
            Ptr< hpoint<double> > *retPw   /* => control points */  );



/*----------------------------------------------------------------------*//**
  Creates a surface of revolution                                          */
/*-------------------------------------------------------------------------*/  
template< class T, class CHP >
void GULAPI Revolution(
           point<T> S,  /* point on rotation axis */
           point<T> A,  /* direction of axis, |A| = 1 */
           T theta,     /* rotation angle (in radians) */
           int m,       /* number of control points -1 of rotated curve */
           Ptr< CHP >& Pwj,  /* control points of rotated curve */
           int *ret_nu, /* number of control points -1  per row of the surface*/
           int *ret_pu,  /* degree of surface of revolution (for rows) */
           Ptr<T> *retU, /* knotvector of surface of revolution (for rows)*/
           Ptr< Ptr< hpoint<T> > >*retPw) /*control points of surface of rev. */
{
  Ptr<T> U;
  T dtheta, wm, wj, angle, sines[5], cosines[5], r, dummy;
  int narcs, i, j, n, index;
  point<T> P0, T0, P1, P2, T2, X, Y, O;
  Ptr< Ptr< hpoint<T> > > Pwij;

  if( theta <= rtr<T>::pi() / 2.0 ) narcs = 1;
  else if( theta <= rtr<T>::pi() ) narcs = 2;
  else if( theta <= (T)1.5 * rtr<T>::pi() ) narcs = 3;
  else narcs = 4;
  
  dtheta = theta / ((T)narcs);

  n = 2 * narcs;              /* number of control points per row */
                              /* (rows = circles) */
  U.reserve_pool(n+4);

  Pwij.reserve_pool(m+1);
  for( i = 0; i <= m; i++ )
    Pwij[i].reserve_pool(n+1);

  wm = rtr<T>::cos( dtheta / (T)2 ); /* medium weight for circle segment */
  
  angle = (T)0;
  
  for( i = 1; i <= narcs; i++ )
  {
    angle += dtheta;
    cosines[i] = rtr<T>::cos(angle);
    sines[i] = rtr<T>::sin(angle);
  }  
  for( j = 0; j <= m; j++ )       /* calc rows */  
  {
    wj = weight(Pwj[j]);  /* j. control point of curve to P0 */
    P0 = euclid(Pwj[j]);

    O = ProjectToLine( S, A, P0 );

    X.x = P0.x - O.x;
    X.y = P0.y - O.y;
    X.z = P0.z - O.z;  

    r = normalize( &X, X );
    Y = cross_product( A, X );
    
    Pwij[j][0] = Pwj[j];        /* first control point of row */

    T0 = Y;
    index = 0;
    angle = (T)0;
    
    for( i = 1; i <= narcs; i++ )  /* remaining  control points of row */     
    {
      P2.x = O.x + r*cosines[i]*X.x + r*sines[i]*Y.x;
      P2.y = O.y + r*cosines[i]*X.y + r*sines[i]*Y.y;
      P2.z = O.z + r*cosines[i]*X.z + r*sines[i]*Y.z;

      Pwij[j][index+2].x = wj * P2.x;
      Pwij[j][index+2].y = wj * P2.y;
      Pwij[j][index+2].z = wj * P2.z;      
      Pwij[j][index+2].w = wj;

      T2.x = -sines[i]*X.x + cosines[i]*Y.x;
      T2.y = -sines[i]*X.y + cosines[i]*Y.y;
      T2.z = -sines[i]*X.z + cosines[i]*Y.z;
      
      SVDIntersectLines( P0, T0, P2, T2, &dummy, &dummy, &P1 );

      Pwij[j][index+1].x = wm*wj*P1.x;
      Pwij[j][index+1].y = wm*wj*P1.y;
      Pwij[j][index+1].z = wm*wj*P1.z;
      Pwij[j][index+1].w = wm*wj;

      index += 2;
      
      if( i < narcs )
      {
        P0 = P2;
        T0 = T2;
      }         
    }
  }
  /* --- construct knot vector for rows ---------------------------- */

  j = 2 * narcs + 1;
  
  for( i = 0; i < 3; i++ )
  {
    U[i] = (T)0;
    U[i+j] = (T)1;
  }  
  switch( narcs )
  {
    case 2:
      U[3] = U[4] = (T)0.5;
      break;
      
    case 3:
      U[3] = U[4] = (T)1 / (T)3;
      U[5] = U[6] = (T)2 / (T)3;
      break;
      
    case 4:
      U[3] = U[4] = (T)0.25;  
      U[5] = U[6] = (T)0.5;  
      U[7] = U[8] = (T)0.75;  
      break;
  }
  *ret_nu = n;
  *ret_pu = 2;
  *retU = U;
  *retPw = Pwij;
}
template GULAPI void Revolution(
        point<float> S,  /* point on rotation axis */
        point<float> A,  /* direction of axis, |A| = 1 */
        float theta,     /* rotation angle (in radians) */
        int m,       /* number of control points -1 of rotated curve */
        Ptr< hpoint<float> >& Pwj,   /* control points of rotated curve */
        int *ret_nu, /* number of control points -1  per row of the surface*/
        int *ret_pu,  /* degree of surface of revolution (for rows) */
        Ptr<float> *retU, /* knotvector of surface of revolution (for rows)*/
        Ptr< Ptr< hpoint<float> > >*retPw);/*control points of surface of rev.*/
template GULAPI void Revolution(
       point<double> S,  /* point on rotation axis */
       point<double> A,  /* direction of axis, |A| = 1 */
       double theta,     /* rotation angle (in radians) */
       int m,       /* number of control points -1 of rotated curve */
       Ptr< hpoint<double> >& Pwj,   /* control points of rotated curve */
       int *ret_nu, /* number of control points -1  per row of the surface*/
       int *ret_pu,  /* degree of surface of revolution (for rows) */
       Ptr<double> *retU, /* knotvector of surface of revolution (for rows)*/
       Ptr< Ptr< hpoint<double> > >*retPw);/*control points of surface of rev.*/
template GULAPI void Revolution(
        point<float> S,  /* point on rotation axis */
        point<float> A,  /* direction of axis, |A| = 1 */
        float theta,     /* rotation angle (in radians) */
        int m,       /* number of control points -1 of rotated curve */
        Ptr< point<float> >& Pwj,   /* control points of rotated curve */
        int *ret_nu, /* number of control points -1  per row of the surface*/
        int *ret_pu,  /* degree of surface of revolution (for rows) */
        Ptr<float> *retU, /* knotvector of surface of revolution (for rows)*/
        Ptr< Ptr< hpoint<float> > >*retPw);/*control points of surface of rev.*/
template GULAPI void Revolution(
       point<double> S,  /* point on rotation axis */
       point<double> A,  /* direction of axis, |A| = 1 */
       double theta,     /* rotation angle (in radians) */
       int m,       /* number of control points -1 of rotated curve */
       Ptr< point<double> >& Pwj,   /* control points of rotated curve */
       int *ret_nu, /* number of control points -1  per row of the surface*/
       int *ret_pu,  /* degree of surface of revolution (for rows) */
       Ptr<double> *retU, /* knotvector of surface of revolution (for rows)*/
       Ptr< Ptr< hpoint<double> > >*retPw);/*control points of surface of rev.*/

template GULAPI void Revolution(
        point<float> S,  /* point on rotation axis */
        point<float> A,  /* direction of axis, |A| = 1 */
        float theta,     /* rotation angle (in radians) */
        int m,       /* number of control points -1 of rotated curve */
        Ptr< hpoint2<float> >& Pwj,   /* control points of rotated curve */
        int *ret_nu, /* number of control points -1  per row of the surface*/
        int *ret_pu,  /* degree of surface of revolution (for rows) */
        Ptr<float> *retU, /* knotvector of surface of revolution (for rows)*/
        Ptr< Ptr< hpoint<float> > >*retPw);/*control points of surface of rev.*/
template GULAPI void Revolution(
       point<double> S,  /* point on rotation axis */
       point<double> A,  /* direction of axis, |A| = 1 */
       double theta,     /* rotation angle (in radians) */
       int m,       /* number of control points -1 of rotated curve */
       Ptr< hpoint2<double> >& Pwj,   /* control points of rotated curve */
       int *ret_nu, /* number of control points -1  per row of the surface*/
       int *ret_pu,  /* degree of surface of revolution (for rows) */
       Ptr<double> *retU, /* knotvector of surface of revolution (for rows)*/
       Ptr< Ptr< hpoint<double> > >*retPw);/*control points of surface of rev.*/
template GULAPI void Revolution(
        point<float> S,  /* point on rotation axis */
        point<float> A,  /* direction of axis, |A| = 1 */
        float theta,     /* rotation angle (in radians) */
        int m,       /* number of control points -1 of rotated curve */
        Ptr< point2<float> >& Pwj,   /* control points of rotated curve */
        int *ret_nu, /* number of control points -1  per row of the surface*/
        int *ret_pu,  /* degree of surface of revolution (for rows) */
        Ptr<float> *retU, /* knotvector of surface of revolution (for rows)*/
        Ptr< Ptr< hpoint<float> > >*retPw);/*control points of surface of rev.*/
template GULAPI void Revolution(
       point<double> S,  /* point on rotation axis */
       point<double> A,  /* direction of axis, |A| = 1 */
       double theta,     /* rotation angle (in radians) */
       int m,       /* number of control points -1 of rotated curve */
       Ptr< point2<double> >& Pwj,   /* control points of rotated curve */
       int *ret_nu, /* number of control points -1  per row of the surface*/
       int *ret_pu,  /* degree of surface of revolution (for rows) */
       Ptr<double> *retU, /* knotvector of surface of revolution (for rows)*/
       Ptr< Ptr< hpoint<double> > >*retPw);/*control points of surface of rev.*/



/*-----------------------------------------------------------------------*//**
  Creates a unit sphere (centered at origin)                                */
/*--------------------------------------------------------------------------*/
template< class T >
GULAPI void UnitSphere( 
            int *ret_nu, int *ret_pu, Ptr<T> *retU,
            int *ret_nv, int *ret_pv, Ptr<T> *retV,
            Ptr< Ptr< hpoint<T> > > *retPw ) 
{
  point<T> O,X,Y,S,A;
  T r;
  Ptr<T> U,V;
  Ptr< hpoint<T> > Cw;
  Ptr< Ptr< hpoint<T> > > Pw;
  int nu,nv,pu,pv;

  O.x = (T)0.0;  O.y = (T)0.0;  O.z = (T)0.0;
  Y.x = (T)0.0;  Y.y = (T)1.0;  Y.z = (T)0.0;
  X.x = (T)1.0;  X.y = (T)0.0;  X.z = (T)0.0;
  r = (T)1.0;

  Circle( O, X, Y, r, (T)0, rtr<T>::pi(), &nv, &pv, &V, &Cw );

  /* create rotation surface: */
  S.x = (T)0.0;  S.y = (T)0.0;  S.z = (T)0.0;
  A.x = (T)1.0;  A.y = (T)0.0;  A.z = (T)0.0;

  Revolution( S, A, (T)2*rtr<T>::pi(), nv, Cw, &nu, &pu, &U, &Pw );

  *ret_nu = nu;
  *ret_pu = pu;
  *retU = U;
  *ret_nv = nv;
  *ret_pv = pv;
  *retV = V;
  *retPw = Pw;
}
template GULAPI void UnitSphere( 
            int *ret_nu, int *ret_pu, Ptr<float> *retU,
            int *ret_nv, int *ret_pv, Ptr<float> *retV,
            Ptr< Ptr< hpoint<float> > > *retPw );
template GULAPI void UnitSphere( 
            int *ret_nu, int *ret_pu, Ptr<double> *retU,
            int *ret_nv, int *ret_pv, Ptr<double> *retV,
            Ptr< Ptr< hpoint<double> > > *retPw );

}
 

---