Text archives Help

---

Author: aek
Date: Tue Jul 12 12:22:38 2005
New Revision: 429

Modified:
   trunk/Model/Primitives/SuperEllipsoid.cc
   trunk/Model/Primitives/SuperEllipsoid.h
   trunk/scenes/primtest.cc
Log:
M    scenes/primtest.cc
    - Changed to match changes to SuperEllipsoid.cc/h
    - Changed to allow a single prim (-num 1x1) 

M    Model/Primitives/SuperEllipsoid.cc
M    Model/Primitives/SuperEllipsoid.h
    - Fixed and sped up.  Still a bit of self-shadowing though



Modified: trunk/Model/Primitives/SuperEllipsoid.cc
==============================================================================
--- trunk/Model/Primitives/SuperEllipsoid.cc    (original)
+++ trunk/Model/Primitives/SuperEllipsoid.cc    Tue Jul 12 12:22:38 2005
@@ -1,31 +1,28 @@
+#include <iostream>
+
 #include <Model/Primitives/SuperEllipsoid.h>
 #include <Interface/RayPacket.h>
 #include <Core/Geometry/BBox.h>
 #include <Core/Math/MinMax.h>
 #include <Core/Math/MiscMath.h>
 
-using namespace Manta;
 using namespace std;
-using SCIRun::Clamp;
+using namespace SCIRun;
+using namespace Manta;
 
 #define GOLDENMEAN 0.61803398874989484820
-#define BRACKETWIDTH 1.e-6
-#define MAXNEWTONITER 25
+#define BRACKETWIDTH 1.e-3
+#define MAXNEWTONITER 10
 
 SuperEllipsoid::SuperEllipsoid(
   Material *material,
-  Point const &center,
-  double radius,
   double alpha,
   double beta )
   : PrimitiveCommon( material, this ),
-    center( center ),
-    radius( radius )
+    two_over_alpha( 2. / alpha ),
+    two_over_beta( 2. / beta ),
+    alpha_over_beta( alpha / beta )
 {
-  inv_radius = 1. / radius;
-  two_over_alpha = 2. / alpha;
-  two_over_beta = 2. / beta;
-  alpha_over_beta = alpha / beta;
 }
 
 SuperEllipsoid::~SuperEllipsoid()
@@ -33,7 +30,7 @@
 }
 
 inline double SuperEllipsoid::functionValue(
-  Vector const &location ) const
+  Point const &location ) const
 {
   double x_power = pow( fabs( location.x() ), two_over_alpha );
   double y_power = pow( fabs( location.y() ), two_over_alpha );
@@ -42,7 +39,7 @@
 }
 
 inline Vector const SuperEllipsoid::functionGradient(
-  Vector const &location,
+  Point const &location,
   double &value ) const
 {
   double x_power = pow( fabs(location.x() ), two_over_alpha );
@@ -56,7 +53,7 @@
 }
 
 inline Vector const SuperEllipsoid::logarithmFunctionGradient(
-  Vector const &location,
+  Point const &location,
   double &value ) const
 {
   double x_power = pow( fabs(location.x() ), two_over_alpha );
@@ -73,64 +70,49 @@
     PreprocessContext const &,
     BBox &bbox ) const
 {
-  bbox.extendBySphere( center, radius );
+    bbox.extendByPoint( Point( -1., -1., -1. ) );
+    bbox.extendByPoint( Point( 1., 1., 1. ) );
 }
 
 void SuperEllipsoid::intersect(const RenderContext&, RayPacket& rays) const
 {
 
-  rays.normalizeDirections();
   rays.computeInverseDirections();
-  for ( int i = 0; i < rays.getSize(); ++i )
-  {
-    RayPacket::Element &e( rays.get( i ) );
-    Vector offset_center = e.ray.origin() - center;
+  for( int i = 0; i < rays.getSize(); ++i ) {
+    RayPacket::Element& e( rays.get(i) );
 
     // First check if the ray hits the bounding box and whether it could
     // remotely produce a hit of interest.
     // TODO: Maybe factor this out into a common routine?
-    double tnear, tfar, t1, t2;
-    t1 = ( -radius - offset_center.x() ) * e.inverseDirection.x();
-    t2 = ( radius - offset_center.x() ) * e.inverseDirection.x();
-    if( t1 > t2 ) {
-      double temp = t1;
-      t1 = t2;
-      t2 = temp;
-    }
-    tnear = t1;
-    tfar = t2;
-    t1 = ( -radius - offset_center.y() ) * e.inverseDirection.y();
-    t2 = ( radius - offset_center.y() ) * e.inverseDirection.y();
-    if( t1 > t2 ) {
-      double temp = t1;
-      t1 = t2;
-      t2 = temp;
-    }
-    using SCIRun::Max;
-    using SCIRun::Min;
-    tnear = Max( t1, tnear );
-    tfar = Min( t2, tfar );
-    t1 = ( -radius - offset_center.z() ) * e.inverseDirection.z();
-    t2 = ( radius - offset_center.z() ) * e.inverseDirection.z();
-    if( t1 > t2 ) {
-      double temp = t1;
-      t1 = t2;
-      t2 = temp;
-    }
-    tnear = Max( Max( t1, tnear ), T_EPSILON );
-    tfar = Min( Min( t2, tfar ), e.hitInfo.minT() );
-
-    if ( tfar < T_EPSILON || tnear >= e.hitInfo.minT() )
+    double t1 = ( -1.01 - e.ray.origin().x() ) * e.inverseDirection.x();
+    double t2 = ( 1.01 - e.ray.origin().x() ) * e.inverseDirection.x();
+    if ( t1 > t2 )
+      SWAP( t1, t2 );
+    double tnear = t1;
+    double tfar = t2;
+    t1 = ( -1.01 - e.ray.origin().y() ) * e.inverseDirection.y();
+    t2 = ( 1.01 - e.ray.origin().y() ) * e.inverseDirection.y();
+    if ( t1 > t2 )
+      SWAP( t1, t2 );
+    tnear = SCIRun::Max( t1, tnear );
+    tfar = SCIRun::Min( t2, tfar );
+    t1 = ( -1.01 - e.ray.origin().z() ) * e.inverseDirection.z();
+    t2 = ( 1.01 - e.ray.origin().z() ) * e.inverseDirection.z();
+    if ( t1 > t2 )
+      SWAP( t1, t2 );
+    tnear = SCIRun::Max( SCIRun::Max( t1, tnear ), T_EPSILON );
+    tfar = SCIRun::Min( SCIRun::Min( t2, tfar ), e.hitInfo.minT() );
+    if ( tnear > tfar || tfar <= T_EPSILON || tnear >= e.hitInfo.minT() )
       continue;
 
     // A few preliminary early exit tests...
     double near_value, far_value;
     double near_deriv = Dot( functionGradient(
-                                 ( offset_center + tnear * e.ray.direction() 
) * inv_radius,
+                                 ( e.ray.origin() + tnear * 
e.ray.direction() ),
                                  near_value ),
                              e.ray.direction() );
     double far_deriv = Dot( functionGradient(
-                                ( offset_center + tfar * e.ray.direction() ) 
* inv_radius,
+                                ( e.ray.origin() + tfar * e.ray.direction() 
),
                                 far_value ),
                             e.ray.direction() );
     if ( ( near_value < T_EPSILON && far_value < T_EPSILON ) ||
@@ -143,32 +125,28 @@
       double a_bracket = tnear;
       double b_bracket = tfar;
       double left = GOLDENMEAN * a_bracket + ( 1. - GOLDENMEAN ) * b_bracket;
-      double left_value = functionValue(
-        ( offset_center + left * e.ray.direction() ) * inv_radius );
+      double left_value = functionValue( e.ray.origin() + left * 
e.ray.direction() );
       double right = ( 1. - GOLDENMEAN ) * a_bracket + GOLDENMEAN * 
b_bracket;
-      double right_value = functionValue(
-        ( offset_center + right * e.ray.direction() ) * inv_radius );
+      double right_value = functionValue( e.ray.origin() + right * 
e.ray.direction() );
       while( fabs( b_bracket - a_bracket ) > BRACKETWIDTH ) {
         if ( left_value < right_value ) {
           b_bracket = right;
           right = left;
           right_value = left_value;
           left = GOLDENMEAN * a_bracket + ( 1. - GOLDENMEAN ) * b_bracket;
-          left_value = functionValue(
-            ( offset_center + left * e.ray.direction() ) * inv_radius );
+          left_value = functionValue( e.ray.origin() + left * 
e.ray.direction() );
         } else {
           a_bracket = left;
           left = right;
           left_value = right_value;
           right = ( 1. - GOLDENMEAN ) * a_bracket + GOLDENMEAN * b_bracket;
-          right_value = functionValue(
-            ( offset_center + right * e.ray.direction() ) * inv_radius );
+          right_value = functionValue( e.ray.origin() + right * 
e.ray.direction() );
         }
       }
 
       // If our minimum is positive, we missed the superquadric - it
       // couldn't have crossed zero.
-      if ( right_value >= T_EPSILON )
+      if ( right_value >= -T_EPSILON )
         continue;
 
       // Narrow the range with the location of the minimum found
@@ -191,7 +169,7 @@
     double troot = ( tnear + tfar ) * 0.5;
     double root_value;
     double root_deriv = Dot( logarithmFunctionGradient(
-                                 ( offset_center + troot * e.ray.direction() 
) * inv_radius,
+                                 e.ray.origin() + troot * e.ray.direction(),
                                  root_value ),
                              e.ray.direction() );
     int iterations = 0;
@@ -205,11 +183,11 @@
         tnear = troot;
         near_value = root_value;
       }
-      troot = troot - root_value / root_deriv;
+      troot -= root_value / root_deriv;
       if ( troot <= tnear || troot >= tfar )
         troot = ( tnear + tfar ) * 0.5;
       root_deriv = Dot( logarithmFunctionGradient(
-                            ( offset_center + troot * e.ray.direction() ) * 
inv_radius,
+                            e.ray.origin() + troot * e.ray.direction(),
                             root_value ),
                         e.ray.direction() );
     }
@@ -225,14 +203,14 @@
   RenderContext const &,
   RayPacket &rays ) const
 {
-  // FIXME: Calculate proper normal once intersection code works - sphere
-  // normal close enough for now...
   double ignored;
   rays.computeHitPositions();
   for( int i = 0; i < rays.getSize(); i++ ) {
     RayPacket::Element &e( rays.get( i ) );
-    e.normal = functionGradient( ( e.hitPosition - center ) * inv_radius, 
ignored );
+    e.normal = functionGradient( e.hitPosition, ignored );
+    e.normal.normalize();
   }
+  rays.setFlag(RayPacket::HaveNormals | RayPacket::HaveUnitNormals);
 }
 
 void SuperEllipsoid::computeTexCoords2(
@@ -242,13 +220,12 @@
   rays.computeHitPositions();
   for( int i = 0; i < rays.getSize(); i++ ) {
     RayPacket::Element &e( rays.get( i ) );
-    Vector n = ( e.hitPosition - center ) * inv_radius;
-    double angle = Clamp( n.z(), -1.0, 1.0 );
+    double angle = Clamp( e.hitPosition.z(), -1.0, 1.0 );
     double theta = acos( angle );
-    double phi = atan2( n.x(), n.y() );
+    double phi = atan2( e.hitPosition.x(), e.hitPosition.y() );
     e.texCoords = Point( ( phi + M_PI ) * ( 1. / ( 2 * M_PI ) ),
                          theta * ( 1. / M_PI ),
-                         0 );
+                         0. );
   }
   rays.setFlag( RayPacket::HaveTexture2 | RayPacket::HaveTexture3 );
 }
@@ -260,13 +237,12 @@
   rays.computeHitPositions();
   for( int i = 0; i < rays.getSize(); i++ ) {
     RayPacket::Element &e = rays.get( i );
-    Vector n = ( e.hitPosition - center ) * inv_radius;
-    double angle = Clamp( n.z(), -1.0, 1.0 );
+    double angle = Clamp( e.hitPosition.z(), -1.0, 1.0 );
     double theta = acos( angle );
-    double phi = atan2( n.x(), n.y() );
+    double phi = atan2( e.hitPosition.x(), e.hitPosition.y() );
     e.texCoords = Point( ( phi + M_PI ) * ( 1. / ( 2 * M_PI ) ),
                          theta * ( 1. / M_PI ),
-                         0 );
+                         0. );
   }
   rays.setFlag( RayPacket::HaveTexture2 | RayPacket::HaveTexture3 );
 }

Modified: trunk/Model/Primitives/SuperEllipsoid.h
==============================================================================
--- trunk/Model/Primitives/SuperEllipsoid.h     (original)
+++ trunk/Model/Primitives/SuperEllipsoid.h     Tue Jul 12 12:22:38 2005
@@ -9,7 +9,7 @@
 namespace Manta {
   class SuperEllipsoid : public PrimitiveCommon, public TexCoordMapper {
   public:
-    SuperEllipsoid(Material* material, const Point& center, double radius, 
double alpha, double beta);
+    SuperEllipsoid(Material* material, double alpha, double beta);
     virtual ~SuperEllipsoid();
 
     virtual void computeBounds(const PreprocessContext& context,
@@ -21,12 +21,9 @@
     virtual void computeTexCoords3(const RenderContext& context,
                                    RayPacket& rays) const;
   private:
-    double functionValue(const Vector& location) const;
-    Vector const functionGradient(const Vector& location, double& value ) 
const;
-    Vector const logarithmFunctionGradient(const Vector& location, double& 
value ) const;
-    Point center;
-    double radius;
-    double inv_radius;
+    double functionValue(const Point& location) const;
+    Vector const functionGradient(const Point& location, double& value ) 
const;
+    Vector const logarithmFunctionGradient(const Point& location, double& 
value ) const;
     double two_over_alpha;
     double two_over_beta;
     double alpha_over_beta;

Modified: trunk/scenes/primtest.cc
==============================================================================
--- trunk/scenes/primtest.cc    (original)
+++ trunk/scenes/primtest.cc    Tue Jul 12 12:22:38 2005
@@ -174,8 +174,8 @@
       for(int j=0;j<numy;j++){
         int idx = j*numx+i;
         double radius = (idx+1)/static_cast<double>(numx*numy)*scale/max;
-        Point p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+        Point p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 0.)*scale*2,
+                (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 0.)*scale*2,
                 0);
   Primitive* prim = new Sphere( matl, p, radius );
   if ( mapr )
@@ -183,28 +183,12 @@
         group->add( prim );
       }
     }
-  } else if (primtype == "simplesuperellipsoid"){
-    for(int i=0;i<numx;i++){
-      for(int j=0;j<numy;j++){
-        int idx = j*numx+i;
-        double radius = (idx+1)/static_cast<double>(numx*numy)*scale/max;
-        double alpha = (i+1)/static_cast<double>(numx)*2.;
-        double beta = (j+1)/static_cast<double>(numy)*2.;
-        Point p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                (j/static_cast<double>(numy-1) - 0.5)*scale*2,
-                0);
-        Primitive* prim = new SuperEllipsoid( matl, p, radius, alpha, beta );
-        if ( mapr )
-            prim->setTexCoordMapper( mapr );
-        group->add( prim );
-      }
-    }
   } else if(primtype == "simplebox"){
     Vector p2(scale/max, scale/max, scale/max);
     for(int i=0;i<numx;i++){
       for(int j=0;j<numy;j++){
-        Point p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+        Point p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 0.)*scale*2,
+                (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 0.)*scale*2,
                 0);
   Primitive* prim = new Cube( matl, p-p2, p2.x()*1.156, p2.y()*1.156, 
p2.z()*1.156 );
   if ( mapr )
@@ -253,6 +237,11 @@
     if ( mapr )
         prim->setTexCoordMapper( mapr );
     spinprim = prim;
+  } else if (primtype == "superellipsoid"){
+    Primitive* prim = new SuperEllipsoid(matl, 0.5, 1.5);
+    if ( mapr )
+        prim->setTexCoordMapper( mapr );
+    spinprim = prim;
   } else if(primtype == "ply"){
     AffineTransform t;
     t.initWithScale(Vector(scale/max*30, scale/max*30, scale/max*30));
@@ -273,11 +262,11 @@
     if(arraytype == "spin"){
       for(int i=0;i<numx;i++){
         for(int j=0;j<numy;j++){
-          Vector p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                   (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+          Vector p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 
0.)*scale*2,
+                   (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 
0.)*scale*2,
                    0);
-          double a1 = i/static_cast<double>(numx-1)*M_PI*2;
-          double a2 = j/static_cast<double>(numy-1)*M_PI*2;
+          double a1 = i/static_cast<double>(numx)*M_PI*2;
+          double a2 = j/static_cast<double>(numy)*M_PI*2;
           AffineTransform t;
           t.initWithIdentity();
           t.rotate(Vector(0,1,0), a1);
@@ -289,8 +278,8 @@
     } else if(arraytype == "shift"){
       for(int i=0;i<numx;i++){
         for(int j=0;j<numy;j++){
-          Vector p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                   (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+          Vector p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 
0.)*scale*2,
+                   (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 
0.)*scale*2,
                    0);
           group->add(new InstanceT(spinprim, p));
         }
@@ -298,8 +287,8 @@
     } else if(arraytype == "scale"){
       for(int i=0;i<numx;i++){
         for(int j=0;j<numy;j++){
-          Vector p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                   (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+          Vector p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 
0.)*scale*2,
+                   (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 
0.)*scale*2,
                    0);
           int idx = j*numx+i;
           double scale = (idx+1)/static_cast<double>(numx*numy);
@@ -309,8 +298,8 @@
     } else if(arraytype == "nuscale"){
       for(int i=0;i<numx;i++){
         for(int j=0;j<numy;j++){
-          Vector p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                   (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+          Vector p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 
0.)*scale*2,
+                   (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 
0.)*scale*2,
                    0);
           double xscale = (i+1)/static_cast<double>(numx);
           double yscale = (j+1)/static_cast<double>(numy);
@@ -320,11 +309,11 @@
     } else if(arraytype == "spinscale"){
       for(int i=0;i<numx;i++){
         for(int j=0;j<numy;j++){
-          Vector p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                   (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+          Vector p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 
0.)*scale*2,
+                   (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 
0.)*scale*2,
                    0);
-          double a1 = i/static_cast<double>(numx-1)*M_PI*2;
-          double a2 = j/static_cast<double>(numy-1)*M_PI*2;
+          double a1 = i/static_cast<double>(numx)*M_PI*2;
+          double a2 = j/static_cast<double>(numy)*M_PI*2;
           int idx = j*numx+i;
           double scale = (idx+1)/static_cast<double>(numx*numy);
           AffineTransform t;
@@ -339,11 +328,11 @@
     } else if(arraytype == "spinscale2"){
       for(int i=0;i<numx;i++){
         for(int j=0;j<numy;j++){
-          Vector p((i/static_cast<double>(numx-1) - 0.5)*scale*2,
-                   (j/static_cast<double>(numy-1) - 0.5)*scale*2,
+          Vector p((numx>1 ? i/static_cast<double>(numx-1) - 0.5 : 
0.)*scale*2,
+                   (numy>1 ? j/static_cast<double>(numy-1) - 0.5 : 
0.)*scale*2,
                    0);
-          double a1 = i/static_cast<double>(numx-1)*M_PI*2;
-          double a2 = j/static_cast<double>(numy-1)*M_PI*2;
+          double a1 = i/static_cast<double>(numx)*M_PI*2;
+          double a2 = j/static_cast<double>(numy)*M_PI*2;
           int idx = j*numx+i;
           double scale = (idx+1)/static_cast<double>(numx*numy);
           AffineTransform t;

---

• [MANTA] r429 - in trunk: Model/Primitives scenes, aek, 07/12/2005

---

Archive powered by MHonArc 2.6.16.