Manta Interactive Ray Tracer Development Mailing List

Text archives Help


[MANTA] r1193 - in trunk: Model/Groups SwigInterface


Chronological Thread 
  • From: abe@sci.utah.edu
  • To: manta@sci.utah.edu
  • Subject: [MANTA] r1193 - in trunk: Model/Groups SwigInterface
  • Date: Thu, 28 Sep 2006 13:11:44 -0600 (MDT)

Author: abe
Date: Thu Sep 28 13:11:43 2006
New Revision: 1193

Added:
   trunk/Model/Groups/TriangleIW.cc   (contents, props changed)
   trunk/Model/Groups/TriangleIW.h   (contents, props changed)
Modified:
   trunk/Model/Groups/Build_KDTreeDyn.cc
   trunk/Model/Groups/Build_KDTreeDyn.h
   trunk/Model/Groups/CMakeLists.txt
   trunk/Model/Groups/DynBVH.cc
   trunk/Model/Groups/KDTreeDyn.cc
   trunk/Model/Groups/KDTreeDyn.h
   trunk/SwigInterface/KDTreeDynFrame.py
   trunk/SwigInterface/kdtreedyn.py
   trunk/SwigInterface/wxManta.py
Log:

python ../SwigInterface/kdtreedyn.py -f /local/abe/models/defender-ext.obj


Added kdtreedyn invocation script. Provides KDTree explorer in new menu.
M    SwigInterface/kdtreedyn.py
M    SwigInterface/KDTreeDynFrame.py

Modified wxManta.py to facilitate adding menu items from external scripts.
M    SwigInterface/wxManta.py

Added O(nlogn) KDTree build. (Will shortly rename file to reflect class name.)
M    Model/Groups/Build_KDTreeDyn.cc
M    Model/Groups/CMakeLists.txt
M    Model/Groups/Build_KDTreeDyn.h
M    Model/Groups/KDTreeDyn.cc
M    Model/Groups/KDTreeDyn.h

Removed dbl_max constants which were causing compilation warnings. Replaced 
with numeric_limits
M    Model/Groups/DynBVH.cc

Moved TriangleIW into it's own file.
AM   Model/Groups/TriangleIW.cc
AM   Model/Groups/TriangleIW.h

Modified: trunk/Model/Groups/Build_KDTreeDyn.cc
==============================================================================
--- trunk/Model/Groups/Build_KDTreeDyn.cc       (original)
+++ trunk/Model/Groups/Build_KDTreeDyn.cc       Thu Sep 28 13:11:43 2006
@@ -31,9 +31,9 @@
 #include <SCIRun/Core/Math/MinMax.h>
 #include <SCIRun/Core/Thread/Time.h>
 
-
 #include <limits>
 #include <iostream>
+#include <iomanip>
 
 using namespace std;
 using namespace SCIRun;
@@ -41,235 +41,544 @@
 
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
-// FIND SPLIT APPROX CC  FIND SPLIT APPROX CC  FIND SPLIT APPROX CC  FIND 
SPLIT
+//
+// This file contains a C/C++ non-simd, non-parallel implementation of
+// Ingo Wald's O(N log N) KDTree build.
+//
+// Abe Stephens.
+//
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
 
-enum { MaxCandidates = 8 };
+// Flags to indicate if area of related triangle should be
+// added/subtracted to left or right sum.
+enum {
+  // Axis of this candidate.
+  AxisMask   = 0x3,
+  AxisX      = 0x0,
+  
+  // LeftBound indicates that the candidate is the minmum bounding 
coordinate,
+  // it's flag is greater than RightBound, because when sorting candidates, 
right
+  // candidates must be operated on first.
+  BoundMask  = 0xC,
+  BoundShift = 2,
+  RightBound = 0x4,
+  Planar     = 0x8,
+  LeftBound  = 0xC
+};
 
-int find_split_approx_cc( const BoundsList &geometry,
-                          const IndexList  &indices,
+typedef unsigned int LocalIndex;
 
-                          const BBox &parent_bounds,
-                          const int dim,
-                           
-                          // Output arguments, 
-                          Real candidate[MaxCandidates],
-                          Real costs    [MaxCandidates],
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// INDEXED CANDIDATE
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
-                          Integer left_size,
-                          Integer right_size )
-{
+// Candidate consisting of a split location, associated index to a primitive,
+// and flag indicating which bound the candidate was generated by.
 
-  
/////////////////////////////////////////////////////////////////////////////
-  // Determine where the candidate split planes are.
-  const Vector extent = (parent_bounds[1] - parent_bounds[0]);
-  const Vector delta  = (parent_bounds[1] - parent_bounds[0]) * 
((Real)1.0)/((Real)MaxCandidates);
-
-  // Real candidates[MaxCandidates];
-  for (int i=0;i<MaxCandidates;++i) {
-    candidate[i] = parent_bounds[0][dim] + (delta[dim] * (Real)i);
+class IndexedCandidate {
+public:
+  IndexedCandidate( Real split_, LocalIndex index_, int flags_ ) :
+    split( split_ ), index( index_ ), flags( flags_ ) { }
+  IndexedCandidate() { };
+  
+  Real         split;
+  LocalIndex   index;
+  int          flags;
+
+  // Flag accessors.
+  int axis() const  { return flags & AxisMask; }
+  int bound() const { return (flags & BoundMask); }
+
+  bool operator< (const IndexedCandidate &c ) const {
+    return 
+      (split < c.split) || ((split == c.split) &&
+                            ((axis() < c.axis()) || ((axis() == c.axis()) &&
+                                                     (bound() < 
c.bound()))));
   }
+};
 
-  // Initialize counters.
-  size_t left_side [MaxCandidates];
-  size_t right_side[MaxCandidates];
-
-  for (int i=0;i<MaxCandidates;++i) {
-    left_side[i]  = 0;
-    right_side[i] = 0;
-  }
+// Candidate List types.
+typedef vector<IndexedCandidate>   IndexedCandidateList;
 
-  
/////////////////////////////////////////////////////////////////////////////
-  // Iterate over the indices
-  for (size_t index=0;index<indices.size();++index) {
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// SPLIT PLANE STRUCT
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
-    const BBox &bounds = geometry[ indices[index] ];
-    
-    // Iterate over the candidates.
-    for (int i=0;i<MaxCandidates;++i) {
-      
-      const Real &candidate_split = candidate[i];
+enum { LeftSide, RightSide }; // Planar flags.
 
-      // Check the left and right split candidates.
-      left_side[i]  += ((bounds.getMin()[dim] <= candidate_split));
-      right_side[i] += ((bounds.getMax()[dim] > candidate_split));
-    }
-  }
+class SplitPlane {
+public:
+  Real cost;
+  Real split;
+  Integer axis;
+  Integer planar; // Which side to count planar candidates on.
+};
 
-  
/////////////////////////////////////////////////////////////////////////////
-  // Compute cost of each split candidate.
-  Real min_cost = std::numeric_limits<Real>::max();
-  int  min_candidate;
 
-  const int dim0 = (dim+1)%3;
-  const int dim1 = (dim+2)%3;
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// Build frame.
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+struct Frame {
+  Frame( const BBox &bounds_,  NodeOffset node_, Integer depth_ ) :
+    bounds( bounds_ ),
+    node( node_ ),
+    depth( depth_ ) { };
 
-  const Real ends_area = extent[dim0] * extent[dim1] * 2;
-  const Real cross_section  = (extent[dim0] + extent[dim1]) * 2;
+  // candidates index's map into the local index_list which maps into the
+  // global face list.
   
-  for (int i=0;i<MaxCandidates;++i) {
+  IndexedCandidateList  candidates;  // Array of candidates.
+  IndexList             index_list;  // Array of global face indices.
+  
+  BBox       bounds;            // Node bounds.
+  NodeOffset node;              // Index offset of allocated node.
+  Integer    depth;             // Depth.
+
+  // (Rely on "index_list.size()" instead.)
+  // Integer    size;              // Total number of bounding boxes (not # 
candidates)
+};
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// Geometry Classification
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+enum { Both, RightOnly, LeftOnly };
+class Classification {
+public:
+  unsigned int side; // LeftOnly, RightOnly, Both.
+  LocalIndex left;   // Mapping for this index in the left child.
+  LocalIndex right;  // Mapping for this index in the right child.
+};
+
+typedef vector<Classification> BoundsClassificationList;
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// Helper Function Prototypes.
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+void find_bbox_candidates( const BoundsList &bounds,
+                           IndexedCandidateList &candidates );
+
+void find_sah_split_all_axes( const int &total_triangles,
+                              const IndexedCandidateList &candidates,
+                              const BBox &node_bounds,
+                              SplitPlane &result );
+
+void classify_candidates( const SplitPlane &result,
+                          const IndexedCandidateList &candidates,            
              
+                          BoundsClassificationList &class_list/*,
+                                                                Integer 
&left_side,
+                                                                Integer 
&right_side */);
 
-    const Real left_area  = (candidate[i] - parent_bounds[0][dim]) * 
cross_section + ends_area;
-    const Real right_area = (parent_bounds[1][dim] - candidate[i]) * 
cross_section + ends_area;
+void splice_indices( const IndexList &source,
+                     BoundsClassificationList &class_map, // Updates left 
and right child index mapping.
+                     IndexList &left_index_list,
+                     IndexList &right_index_list );
 
-    costs[i] =
-      (Real)(left_side[i])*(left_area) +
-      (Real)(right_side[i])*(right_area);
+void splice_candidates( const IndexedCandidateList &source,
+                        const BoundsClassificationList &class_map,
+                        IndexedCandidateList &left_candidates,
+                        IndexedCandidateList &right_candidates );
 
-    if (min_cost > costs[i]) {
-      min_cost = costs[i];
-      min_candidate = i;
+
+
+void debug_print_candidates( const IndexedCandidateList &candidates );
+
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// Find Candidate functions.
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+// Uses planes of each bounding box as split candidates.
+void find_bbox_candidates( // Assume bounds is in the same order as the 
frame index_list.
+                           const BoundsList &bounds,
+                           IndexedCandidateList &candidates )
+{
+  ASSERT(candidates.size() >= (bounds.size()*6));
+
+  int total_planar = 0;
+  int total_leftright = 0;
+  
+  // Find candidates.
+  int count = 0;
+  for (int i=0;i<bounds.size();++i) {
+
+    const BBox &box = bounds[i];
+    
+    // Iterate over dimensions.
+    for (int d=0;d<3;++d) {
+
+      const unsigned int axis_flag = AxisX+d;
+
+      // Check to see if the box is planar in this dimension.
+      if ((box[1][d] - box[0][d]) < T_EPSILON) {
+
+        // Add a planar candidate.
+        candidates[count++] = IndexedCandidate( box[0][d], i, Planar | 
axis_flag );
+        total_planar++;
+      }
+      else {
+        
+        // Right candidate.
+        candidates[count++] = IndexedCandidate( box[1][d], i, RightBound | 
axis_flag );
+
+        // Left candidate.
+        candidates[count++] = IndexedCandidate( box[0][d], i, LeftBound | 
axis_flag );
+
+        total_leftright++;
+      }
     }
   }
 
-  left_size = left_side[min_candidate];
-  right_size = right_side[min_candidate];
+  // Reduce candidate array size.
+  candidates.resize( count );  
+
+  // Sort candidates.
+  sort( candidates.begin(), candidates.end() );
+
+  std::cerr << "Total planar candidates: " << total_planar << std::endl;
+  std::cerr << "Total left/right (same): " << total_leftright << std::endl;
   
-  return min_candidate;
 }
 
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
-// SAH SPLIT  SAH SPLIT  SAH SPLIT  SAH SPLIT  SAH SPLIT  SAH SPLIT  SAH 
SPLIT 
+// FULL SAH SPLIT  FULL SAH SPLIT  FULL SAH SPLIT  FULL SAH SPLIT  FULL SAH 
SPL
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
 
-#if 0
+void find_sah_split_all_axes( const int &total_triangles,
+                              const IndexedCandidateList &candidates,
+                              const BBox &node_bounds,
+                              SplitPlane &result
+                              )
+{
+  // Candidate list must be sorted and contain candidates from all axes.
 
-// It would be better to sort split candidates:
-class Candidate {
-public:
-  Candidate( Real split_, size_t index_, int flags_ ) :
-    split( split_ ), index( index_ ), flags( flags_ ) { };
+  
/////////////////////////////////////////////////////////////////////////////
+  // Precompute extents and cross section area for each dimension.
+  Vector ends_area;
+  Vector cross_section;
+  
+  const Vector extent = node_bounds[1] - node_bounds[0];  
+  for (int axis=0;axis<3;++axis) {
+    const int dim0 = (axis+1)%3;
+    const int dim1 = (axis+2)%3;
+    ends_area[axis]      = extent[dim0] * extent[dim1] * 2;
+    cross_section[axis]  = (extent[dim0] + extent[dim1]) * 2;
+  }
   
-  Real   split;
-  unsigned int index;
+  
/////////////////////////////////////////////////////////////////////////////
+  // Initialize counters for each dimension.
+  Integer N_left[]   = { 0, 0, 0 };
+  Integer N_planar[] = { 0, 0, 0 };
+  Integer N_right[]  = { total_triangles, total_triangles, total_triangles };
 
-  // Flags to indicate if area of related triangle should be
-  // added/subtracted to left or right sum.
-  enum {LeftBound,Middle,RightBound};
-  int    flags;
+  // Initialize min cost.
+  result.cost = std::numeric_limits<Real>::max();
+  
+  
/////////////////////////////////////////////////////////////////////////////
+  // Iterate over all candidates.
+  for (int i=0;i<candidates.size();++i) {
 
-  inline bool operator < (const Candidate &c ) const { return split < 
c.split; }
-};
+    // Determine the split and axis of the candidate.
+    const Integer c_axis  = candidates[i].axis();
+    const Real    c_split = candidates[i].split;
+
+    
///////////////////////////////////////////////////////////////////////////
+    // Find runs of candidates with the same split on a the same axis.
+    Integer run_right = 0;
+    while (
+           i<candidates.size() &&
+           candidates[i].bound() == RightBound && 
+           candidates[i].axis()  == c_axis &&
+           candidates[i].split   == c_split           
+           ) {
+      run_right++;
+      i++;
+    }    
+    
+    Integer run_planar = 0;
+    while (
+           i<candidates.size() &&
+           candidates[i].bound() == Planar && 
+           candidates[i].axis()  == c_axis &&
+           candidates[i].split   == c_split           
+           ) {
+      run_planar++;
+      i++;
+    }        
+
+    Integer run_left = 0;
+    while (
+           i<candidates.size() &&
+           candidates[i].bound() == LeftBound &&
+           candidates[i].axis()  == c_axis &&
+           candidates[i].split   == c_split            
+           ) {
+      run_left++;
+      i++;
+    }
 
+    
///////////////////////////////////////////////////////////////////////////
+    // Pre-Increment Right Side and Planar count.
+    
+    // Update the number of bounds exactly on the split.
+    N_planar[c_axis] = run_planar;
 
-void sah_find_candidates( TriangleList &geometry, IndexList axis[3] ) {
+    // Count any bounds ending on the split as not being on the right side.
+    N_right[c_axis] -= run_planar; 
+    N_right[c_axis] -= run_right;  
+
+    
///////////////////////////////////////////////////////////////////////////
+    // Apply Surface Area Heuristic.
+
+    // Compute child node areas.
+    const Real left_area  =
+      (c_split - node_bounds[0][c_axis]) * cross_section[c_axis] + 
ends_area[c_axis];
+    const Real right_area =
+      (node_bounds[1][c_axis] - c_split) * cross_section[c_axis] + 
ends_area[c_axis];
+    
+    // Compute SAH cost twice, once with the planar triangles on the left 
and once
+    // with the planar triangles on the right.
+    const Real left_planar_cost =
+      left_area *  ((Real)N_left[c_axis] + (Real)N_planar[c_axis]) +
+      right_area * (Real)N_right[c_axis];
+
+    const Real right_planar_cost =
+      right_area * ((Real)N_right[c_axis] + (Real)N_planar[c_axis]) +
+      left_area *  (Real)N_left[c_axis];
+    
+    // Update min cost.
+    if (left_planar_cost < result.cost) {
+      result.cost  = left_planar_cost;
+      result.axis  = c_axis;
+      result.split = c_split;     
+      result.planar = LeftSide;
+    }
 
-  // Find candidates.
-  // for (int i=0;i<)
+    if (right_planar_cost < result.cost) {
+      result.cost  = right_planar_cost;
+      result.axis  = c_axis;
+      result.split = c_split;      
+      result.planar = RightSide;
+    }
 
-  // Sort candidates.
-  start_time = Time::currentSeconds();
-  for (int i=0;i<3;++i) {
-    sort(axis[i].begin(),axis[i].end());
+    
///////////////////////////////////////////////////////////////////////////
+    // Post-increment Left Side.
+    
+    // Count any bounds starting on the split as left side.
+    N_left[c_axis] += run_left;
+    N_left[c_axis] += run_planar;
   }
-  end_time = Time::currentSeconds();
-
-  // Eliminate duplicate candidates.
-  
 
-  std::cerr << "sah sort time " << (end_time-start_time) << std::endl;
 }
 
-void sah_split( const BBox &node_bounds, int dim ) {
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// Classify Candidates
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
 
-  // Precompute node surface area constants.
-  const Vector extent = (node_bounds[1] - node_bounds[0]);
-  const int dim0 = (dim+1)%3;
-  const int dim1 = (dim+2)%3;
+// This function produces a list of triangle classifications and counts the
+// total amount of _geometry_ (not candidates) on each side.
+void classify_candidates( const SplitPlane &result,
+                          const IndexedCandidateList &candidates,
+                          
+                          BoundsClassificationList &class_list/*,
+                                                                Integer 
&left_side,
+                                                                Integer 
&right_side */)
+{
   
-  const Real ends_area = extent[dim0] * extent[dim1] * 2;
-  const Real cross_section = 2*(extent[dim0] + extent[dim1]);
+  // class_list size must be equal to number of input bounds indexed by
+  // candidate list.
 
-  // 
+  
/////////////////////////////////////////////////////////////////////////////
+  // Classify all as both.
+  for (int i=0;i<class_list.size();++i) {
+    class_list[i].side = Both;
+  }
 
-  
-  for (int index=0;index<triangle_array.size();index+=delta) {
+  // Reset counters.
+  int both = class_list.size(); // All _trianges_ classified as both.
+  Integer left_side = 0;
+  Integer right_side = 0; 
 
-    const Triangle & triangle = triangle_array[index];
+  
/////////////////////////////////////////////////////////////////////////////
+  // Iterate over candidates.
+  for (int i=0;i<candidates.size();++i) {
     
-    // Iterate over bounding box faces.
-      
-      Real candidate_split = triangle_array[index].bounds[0][dim];
+    // Only classify triangles using candidates on the split axis.
+    if (candidates[i].axis() == result.axis) {
 
-      Real left_side = 0;
-      Real right_side = 0;
+      const unsigned int &index = candidates[i].index;
       
-      // Iterate over all triangles and sum up surface area.
-      for (int i=0;i<triangle_array.size();++i) {
-        if ((triangle_array[i].bounds[0][dim] <= candidate_split)) {
-          left_side += triangle_array[i].area;
-        }
-        if ((triangle_array[i].bounds[1][dim] > candidate_split)) {
-          right_side += triangle_array[i].area;
-        }
+      // If the right bound is less than the split, then the triangle is
+      // only on the left side.
+      if (candidates[i].bound() == RightBound &&
+          candidates[i].split <= result.split) {
+        
+        class_list[index].side = LeftOnly;
+
+        // Update count.
+        left_side ++;
+        both --;
       }
 
-      // Compute the node surface area      
-      const Real left_area  = (candidate_split - model_bounds[0][dim]) * 
cross_section + ends_area;
-      const Real right_area = (model_bounds[1][dim] - candidate_split) * 
cross_section + ends_area;
-      
-      const Real cost =
-        (Real)(left_side)*(left_area) +
-        (Real)(right_side)*(right_area);
-
-      // Output the candidate.
-      file << candidate_split << "\t"
-           << cost
-           << std::endl;
+      // If the left bound is greater than the split, then the triangle
+      // is only on the right side.
+      else if ((candidates[i].bound() == LeftBound) &&
+               (candidates[i].split >= result.split)) {
+        
+        class_list[index].side = RightOnly;
+
+        right_side ++;
+        both --;        
+      }
+
+      else if (candidates[i].bound() == Planar) {
+
+        // If the planar candidate is less than the split plane, or if it's 
on
+        // the split plane and we are counting the Planar triangles on the 
left
+        // side.
+        if ((candidates[i].split < result.split) || // OR
+            ((candidates[i].split == result.split) && (result.planar == 
LeftSide))) {
+          
+          class_list[index].side = LeftOnly;
+
+          left_side ++;
+          both --;
+        }
+
+        // If the planar candidate is greater than the split plane, or if 
it's on
+        // the split and we are counting planar triangles on the right side.
+        else if ((candidates[i].split > result.split) || // OR
+                 ((candidates[i].split == result.split) && (result.planar == 
RightSide))) {
+
+          class_list[index].side = RightOnly;
 
+          right_side ++;
+          both --;
+        }
+      }      
+    }
+  }
+
+  // Increment each side by the both count.
+  if (both > 0) {
+    right_side += both;
+    left_side += both;
   }
-  file.close();  
 }
 
-#endif 
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+// SPLICE INDICES  SPLICE INDICES  SPLICE INDICES  SPLICE INDICES  SPLICE 
INDIC
+///////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+
+void splice_indices( const IndexList &source,
+                     BoundsClassificationList &class_map, // Updates left 
and right child index mapping.
+                     IndexList &left_index_list,
+                     IndexList &right_index_list )
+{
+
+  // Purge both lists.
+  left_index_list.clear();
+  right_index_list.clear();
+
+  for (int i=0;i<source.size();++i) {
+
+    // Add the triangle indices to either side (or both)
+    // and record the new position in the child frame index list.
+    
+    switch (class_map[i].side) { 
+    case LeftOnly:
+      class_map[i].left = left_index_list.size();
+      left_index_list.push_back( source[i] );
+      break;
+    case RightOnly:
+      class_map[i].right = right_index_list.size();
+      right_index_list.push_back( source[i] );
+      break;
+    case Both:
+      class_map[i].left = left_index_list.size();
+      left_index_list.push_back( source[i] );
+
+      class_map[i].right = right_index_list.size();
+      right_index_list.push_back( source[i] );
+      break;
+    };
+  }  
+  
+}
 
 
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
-// BUCKET SORT  BUCKET SORT  BUCKET SORT  BUCKET SORT  BUCKET SORT  BUCKET 
SORT
+// SPLICE CANDIDATES  SPLICE CANDIDATES  SPLICE CANDIDATES  SPLICE 
CANDIDATES  
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
 
-void bucket_sort_cc( const BoundsList  &geometry,
-                     
-                     const IndexList &source_indices,
-                     IndexList &left_indices,
-                     IndexList &right_indices,
-                     
-                     const Real split, const int dim )
+// Splice the source candidate list into a left list and a right list.
+// The sift is based on the classification of each indexed triangle by
+// each candidate and therefore is independent of candidate axis.
+
+void splice_candidates( const IndexedCandidateList &source,
+                        const BoundsClassificationList &class_map,
+                        IndexedCandidateList &left_candidates,
+                        IndexedCandidateList &right_candidates )
 {
-  const int source_size = source_indices.size();
-  
-  left_indices.clear();
-  right_indices.clear();
 
-  
-  for (int i=0;i<source_size;++i) {
+  // Purge both lists.
+  left_candidates.clear();
+  right_candidates.clear();
 
-    const Integer index = source_indices[i];
-    const BBox &bounds = geometry[index];
+  for (int i=0;i<source.size();++i) {
 
-    if (bounds[0][dim] < split) {
-      left_indices.push_back( index );
-    }
+    const IndexedCandidate &candidate = source[i];
+    const Classification &c = class_map[candidate.index];
     
-    if (bounds[1][dim] > split)
-      right_indices.push_back( index );
+    switch (c.side) { 
+    case LeftOnly:
+      left_candidates.push_back( candidate ); // Push the candidate on one 
side.
+      left_candidates.back().index = c.left;  // Update the candidate index 
mapping.     
+      break;
+    case RightOnly:
+      right_candidates.push_back( candidate );
+      right_candidates.back().index = c.right;
+      break;
+    case Both:
+      left_candidates.push_back( candidate );
+      left_candidates.back().index = c.left;
+      right_candidates.push_back( candidate );
+      right_candidates.back().index = c.right;
+      break;
+    };
   }
-
 }
 
+
+    
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
-// BUILD APPROX CC  BUILD APPROX CC  BUILD APPROX CC  BUILD APPROX CC  BUILD 
AP
+// FULL SAH BUILD  FULL SAH BUILD  FULL SAH BUILD  FULL SAH BUILD  FULL SAH 
BUI
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
 
-
-void Approx_cc::build( KDTreeDyn *kdtree  ) throw (SCIRun::Exception&) {
+void Build_OnLogn_cc::build( KDTreeDyn *kdtree  ) throw (SCIRun::Exception&) 
{
 
   Real start_time = Time::currentSeconds();
   
@@ -280,38 +589,39 @@
   Integer total_nodes = 0;
   Integer total_leaves = 0;
   Integer total_indices = 0;
-  Integer heuristic_trigger = 0;
+  Integer cost_trigger = 0;
 
   
/////////////////////////////////////////////////////////////////////////////
-  // Create a build_stack.
+  
/////////////////////////////////////////////////////////////////////////////
+  // Initialize Build Stack
+  
/////////////////////////////////////////////////////////////////////////////
+  
/////////////////////////////////////////////////////////////////////////////
+
   vector<Frame *> build_stack;
   build_stack.reserve( 128 );
   
   
/////////////////////////////////////////////////////////////////////////////
   // Add the root node to the build stack.
-  {
-    // Allocate the root node.
-    NodeOffset root_offset = kdtree->allocateNodePair();
-
-    // Allocate the root stack frame.
-    Frame *root_frame = new Frame( kdtree->bbox,
-                                   root_offset,
-                                   0 );
-
-    // Create initial index list in the root frame.
-    root_frame->array.resize( bounds_list.size() );
-
-    for (int i=0;i<bounds_list.size();++i)
-      root_frame->array[i] = i;
-    
-    // Push the root frame on the stack.
-    build_stack.push_back( root_frame );
+  
+  // Allocate the root stack frame.
+  const NodeOffset root_offset = kdtree->allocateNodePair();
+  Frame *root_frame = new Frame( kdtree->bbox,
+                                         root_offset,
+                                         0 );  
+
+  // Find initial candidates.
+  root_frame->candidates.resize( bounds_list.size()*6 );
+  find_bbox_candidates( bounds_list, root_frame->candidates );
+
+  // Copy inital indices into root frame's index_list.
+  root_frame->index_list.resize( bounds_list.size() );
+  for (int i=0;i<bounds_list.size();++i) {
+    root_frame->index_list[i] = i;
   }
-
+  
   
/////////////////////////////////////////////////////////////////////////////
-  // Initialize candidate arrays.
-  Real candidate[MaxCandidates];
-  Real costs    [MaxCandidates];
+  // Push the root frame on the stack.
+  build_stack.push_back( root_frame );
 
   
/////////////////////////////////////////////////////////////////////////////
   
/////////////////////////////////////////////////////////////////////////////
@@ -330,31 +640,33 @@
     while (1) {
 
       
///////////////////////////////////////////////////////////////////////////
-      // Check to see if this node is a leaf.
-      if ((frame->array.size() < leaf_threshold) || // Condition: Leaf 
threshold.
+      
///////////////////////////////////////////////////////////////////////////
+      // Create a Leaf.
+      
///////////////////////////////////////////////////////////////////////////
+      
///////////////////////////////////////////////////////////////////////////
+      if ((frame->candidates.size() < leaf_threshold) || // Condition: Leaf 
threshold.
           (frame->depth >= depth_threshold)      || // Condition: Depth 
threshold.
           (break_out)
           ) {
           
         
-        // How to check for heuristic performance?
-        
         
///////////////////////////////////////////////////////////////////////
         // Set leaf fields.
-        KDTreeType::Node &leaf = kdtree->getNode(frame->node);
+        KDTreeDyn::Node &leaf = kdtree->getNode(frame->node);
 
         // Copy the index list.
-        leaf.listBegin = kdtree->allocateLeafList( frame->array );
-        leaf.listLen   = frame->array.size();
+        leaf.listBegin = kdtree->allocateLeafList( frame->index_list );
+        leaf.listLen   = frame->index_list.size();
         leaf.type      = KDTreeType::LeafType;
 
         // DEBUG
         leaf.size      = leaf.listLen;
 
         // Update stats
-        total_indices += frame->array.size();
-        max_leaf_size = SCIRun::Max( max_leaf_size, frame->array.size() );
+        total_indices += frame->index_list.size();
+        max_leaf_size = SCIRun::Max( max_leaf_size, frame->index_list.size() 
);
         ++total_leaves;
+        ++total_nodes;
 
         // Delete the frame.
         delete frame;
@@ -362,63 +674,46 @@
         break;
       }
       else {
-        
-        
///////////////////////////////////////////////////////////////////////////   
 
-        // Find best split.
-        Real min_split;
-        int min_dim;
-        Real min_cost = std::numeric_limits<Real>::max();
 
-        Integer min_left_size;
-        Integer min_right_size;
-        
-        for (int d=0;d<3;++d) {
+        
///////////////////////////////////////////////////////////////////////////
+        
///////////////////////////////////////////////////////////////////////////
+        // Find the best split.
+        
///////////////////////////////////////////////////////////////////////////
+        
///////////////////////////////////////////////////////////////////////////   
 
 
-          Integer left_size = 0;
-          Integer right_size = 0;
-          
-          const int initial = find_split_approx_cc( bounds_list,
-                                                    frame->array,
-                                                    frame->bounds,
-                                                    d,
-                                                    candidate,
-                                                    costs,
-                                                    left_size,
-                                                    right_size );
-
-          // *** Refinement step goes here ***
-
-          // Check for min cost.
-          if (costs[initial] < min_cost) {
-            min_cost  = costs[initial];
-            min_split = candidate[initial];
-            min_dim   = d;
-
-            min_left_size = left_size;
-            min_right_size = right_size;
-          }
-        }
+        // Invoke n log n SAH.
+        SplitPlane plane;
+        find_sah_split_all_axes( frame->index_list.size(),
+                                 frame->candidates,
+                                 frame->bounds,
+                                 plane );
 
+        
///////////////////////////////////////////////////////////////////////
         // Check to see if the heuristic is improving anything.
-        if (min_left_size == min_right_size == frame->array.size()) {
+        const Real parent_cost = frame->bounds.computeArea() * 
(Real)frame->index_list.size();
+        
+        // Compare split cost with that of not splitting.
+        if (plane.cost > parent_cost*cost_threshold) {
 
-          heuristic_trigger++;
+          cost_trigger++;
           
+          // Make a leaf instead.
           break_out = true;
-          continue;          
+          continue;
         }
         
-        const Real split = min_split;
-        const int dim = min_dim;
-        
+        
///////////////////////////////////////////////////////////////////////
         
///////////////////////////////////////////////////////////////////////
         // Allocate two children.
+        
///////////////////////////////////////////////////////////////////////
+        
///////////////////////////////////////////////////////////////////////
+        
         NodeOffset left_child_offset = kdtree->allocateNodePair();
 
         // Update parent pointer.
         KDTreeType::Node &parent_node = kdtree->getNode( frame->node );
         
-        parent_node.size  = frame->array.size();
+        parent_node.size  = frame->index_list.size();
         parent_node.leftOffset = left_child_offset;
 
         // Allocate a new Stack frame for the right child.
@@ -426,28 +721,44 @@
         Frame *right_frame = new Frame( frame->bounds, left_child_offset+1, 
frame->depth+1 );
 
         
///////////////////////////////////////////////////////////////////////
-        // Use the existing frame for the left child and bucket sort.
+        // Splice the candidate list and index list into the children.
 
-        // frame->array.reserve( min_left_size );
-        // right_frame->array.reserve( min_right_size );
-        
-        bucket_sort_cc( bounds_list,
-                        frame->array,
-                        left_frame->array,
-                        right_frame->array,
-                        split, dim );
+        // class_map[i] contains the classification of triangle
+        // index frame->index_list[i].
+
+        BoundsClassificationList class_map( frame->index_list.size() );
+
+        // Classify each triangle.
+        classify_candidates( plane,
+                             frame->candidates,                             
+                             class_map );
+
+        // Splice indices. Add mapping from parent index list to child index 
list(s)
+        splice_indices( frame->index_list,
+                        class_map,
+                        left_frame->index_list,
+                        right_frame->index_list );
+        
+        // Splice candidates.
+        splice_candidates( frame->candidates,
+                           class_map,
+                           left_frame->candidates,
+                           right_frame->candidates );
+
+        
///////////////////////////////////////////////////////////////////////
         
         // Set the node split and axis fields.
-        parent_node.split = split;
-        parent_node.type  = dim;
+        parent_node.split = plane.split;
+        parent_node.type  = plane.axis;
         
         // Update the bounds of both child frames.
-        left_frame->bounds[1][dim] = split;
-        right_frame->bounds[0][dim] = split;
+        left_frame->bounds[1][plane.axis]  = plane.split;
+        right_frame->bounds[0][plane.axis] = plane.split;
 
         // Push the right frame on the stack.
         build_stack.push_back( right_frame );
 
+        // Delete the old build frame.
         Frame *old = frame;
         frame = left_frame;
         delete old;
@@ -465,22 +776,64 @@
   Real end_time = Time::currentSeconds();
 
   std::cerr << "Build information: " << std::endl
-            << "max_depth     " << max_depth     << " (threshold " << 
depth_threshold << ")" << std::endl
-            << "max_leaf_size " << max_leaf_size << " (threshold " << 
leaf_threshold  << ")" << std::endl
-            << "avg_leaf_size " << 
((Real)total_indices)/((Real)total_leaves) << std::endl
-            << "total_nodes   " << total_nodes   << std::endl
-            << "total_leaves  " << total_leaves  << std::endl
-            << "total_indices " << total_indices << std::endl
-            << "total faces   " << kdtree->triangles.size() << std::endl
-            << "heuristic     " << heuristic_trigger << std::endl
-            << "build time    " << (end_time-start_time) << "s" << std::endl;
+            << "max_depth      " << max_depth     << " (threshold " << 
depth_threshold << ")" << std::endl
+            << "max_leaf_size  " << max_leaf_size << " (threshold " << 
leaf_threshold  << ")" << std::endl
+            << "avg_leaf_size  " << 
((Real)total_indices)/((Real)total_leaves) << std::endl
+            << "total_nodes    " << total_nodes   << std::endl
+            << "total_leaves   " << total_leaves  << std::endl
+            << "total_indices  " << total_indices << std::endl
+            << "total faces    " << kdtree->triangles.size() << std::endl
+            << "cost_threshold " << cost_trigger << std::endl
+            << "build time     " << (end_time-start_time) << "s" << 
std::endl;
 }
 
+
+
+
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
-// TEMPLATES  TEMPLATES  TEMPLATES  TEMPLATES  TEMPLATES  TEMPLATES  
TEMPLATES 
+// DEBUG PRINT CANDIDATES  DEBUG PRINT CANDIDATES  DEBUG PRINT CANDIDATES  
DEBU
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
 
+void debug_print_candidates( const IndexedCandidateList &candidates ) {
 
+  std::cerr << "Total candidates: " << candidates.size() << std::endl;
+  
+  // Iterate over all of the candidates.
+  for (int i=0;i<candidates.size();++i) {
 
+    if (i>0 && (candidates[i].split != candidates[i-1].split)) {
+      cerr << "--" << std::endl;
+    }
+    
+    cerr << setiosflags(ios::fixed|ios::showpoint)
+         << setprecision(2)
+         << candidates[i].split << "\t";
+
+    unsigned int axis = candidates[i].flags&AxisMask;
+    axis = axis < 3 ? axis : 3;
+    static const char axis_char[] = { 'X', 'Y', 'Z', 'E' };
+    
+    cerr << axis << "\t";
+    
+    switch (candidates[i].flags&BoundMask) {
+    case RightBound:
+      cerr << "R";
+      break;
+    case Planar:
+      cerr << "P";
+      break;
+    case LeftBound:
+      cerr << "L";
+      break;
+    default:
+      cerr << "E";
+      break;
+    }
+
+    cerr << "\t" << candidates[i].index;
+
+    cerr << endl;
+  }
+}

Modified: trunk/Model/Groups/Build_KDTreeDyn.h
==============================================================================
--- trunk/Model/Groups/Build_KDTreeDyn.h        (original)
+++ trunk/Model/Groups/Build_KDTreeDyn.h        Thu Sep 28 13:11:43 2006
@@ -28,6 +28,7 @@
 
 #include <MantaTypes.h>
 #include <Core/Geometry/BBox.h>
+
 #include <Model/Groups/KDTreeDyn.h>
 #include <SCIRun/Core/Exceptions/Exception.h>
 
@@ -35,32 +36,22 @@
 namespace Manta {
   
   // Class containing build method.
-  class Approx_cc {
+  class Build_OnLogn_cc {
   public:
     typedef KDTreeDyn KDTreeType;
     
-
     enum { MaxCandidates = 8 };
     
     // Thresholds.
     int leaf_threshold;
     int depth_threshold;
+    Real cost_threshold;
 
-    // Build frame.
-    struct Frame {
-      Frame( const BBox &bounds_,  NodeOffset node_, Integer depth_ ) :
-        bounds( bounds_ ), node( node_ ), depth( depth_ ) { };
-      
-      IndexList  array;
-      BBox       bounds;
-      NodeOffset node;
-      Integer    depth;
-    };
-    
     // Constructors.
-    Approx_cc( int leaf_threshold_, int depth_threshold_ ) :
+    Build_OnLogn_cc( int leaf_threshold_, int depth_threshold_, Real 
cost_threshold_ ) :
       leaf_threshold( leaf_threshold_ ),
-      depth_threshold( depth_threshold_ ) { }
+      depth_threshold( depth_threshold_ ),
+      cost_threshold( cost_threshold_ ) { }
     
     // Build methods.
     void build( KDTreeType *kdtree ) throw(SCIRun::Exception&);

Modified: trunk/Model/Groups/CMakeLists.txt
==============================================================================
--- trunk/Model/Groups/CMakeLists.txt   (original)
+++ trunk/Model/Groups/CMakeLists.txt   Thu Sep 28 13:11:43 2006
@@ -36,6 +36,8 @@
      Groups/RealisticBvh.h
      Groups/TransparentKDTree.cc
      Groups/TransparentKDTree.h
+     Groups/TriangleIW.cc
+     Groups/TriangleIW.h
      Groups/VolumeGrid.h
 )
 

Modified: trunk/Model/Groups/DynBVH.cc
==============================================================================
--- trunk/Model/Groups/DynBVH.cc        (original)
+++ trunk/Model/Groups/DynBVH.cc        Thu Sep 28 13:11:43 2006
@@ -2,6 +2,8 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <float.h>
+#include <limits>
+
 using namespace Manta;
 
 // these constants control the SAH cost model
@@ -18,10 +20,10 @@
     IAData ia_data;
     for (int axis = 0; axis < 3; axis++ )
     {
-        ia_data.min_rcp[axis]     =  DBL_MAX;
-        ia_data.max_rcp[axis]     = -DBL_MAX;
-        ia_data.min_org_rcp[axis] =  DBL_MAX;
-        ia_data.max_org_rcp[axis] = -DBL_MAX;
+        ia_data.min_rcp[axis]     =  std::numeric_limits<float>::max();
+        ia_data.max_rcp[axis]     = -std::numeric_limits<float>::max();
+        ia_data.min_org_rcp[axis] =  std::numeric_limits<float>::max();
+        ia_data.max_org_rcp[axis] = -std::numeric_limits<float>::max();
     }
 
     // TODO(boulos): provide an SSE version

Modified: trunk/Model/Groups/KDTreeDyn.cc
==============================================================================
--- trunk/Model/Groups/KDTreeDyn.cc     (original)
+++ trunk/Model/Groups/KDTreeDyn.cc     Thu Sep 28 13:11:43 2006
@@ -47,67 +47,6 @@
 
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////
-// TRIANGLE IW  TRIANGLE IW  TRIANGLE IW  TRIANGLE IW  TRIANGLE IW  TRIANGLE 
IW
-///////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////////////////////
-
-
-
-
-using SCIRun::Abs;
-
-TriangleIW::TriangleIW( const Vector& _p1,
-                        const Vector& _p2,
-                        const Vector& _p3 )
-{
-  Vector normal = Cross(_p2-_p1, _p3-_p1);
-  normal.normalize();
-
-  unsigned int n, u, v;
-
-  n = 0;
-  if ( Abs(normal[1]) > Abs(normal[n])) n = 1;
-  if ( Abs(normal[2]) > Abs(normal[n])) n = 2;
-
-  switch ( n ) {
-  case 0: u = 1; v = 2; break;
-  case 1: u = 2; v = 0; break;
-  default: u = 0; v = 1; break;
-  };
-
-  // copy to struct
-  k = n;
-  n_u = normal[u] / normal[k];
-  n_v = normal[v] / normal[k];
-  n_d = _p1[u] * n_u + _p1[v] * n_v + _p1[k];
-
-  float s;
-
-  c_nu = + _p2[v] - _p1[v];
-  c_nv = - _p2[u] + _p1[u];
-  c_d  = - (_p1[u] * c_nu + _p1[v] * c_nv);
-
-  s = 1.f / ( _p3[u] * c_nu + _p3[v] * c_nv + c_d );
-
-  c_nu *= s;
-  c_nv *= s;
-  c_d  *= s;
-
-  b_nu = + _p3[v] - _p1[v];
-  b_nv = - _p3[u] + _p1[u];
-  b_d  = - (_p1[u] * b_nu + _p1[v] * b_nv);
-
-  s = 1.f / ( _p2[u] * b_nu + _p2[v] * b_nv + b_d );
-
-  b_nu *= s;
-  b_nv *= s;
-  b_d  *= s;
-
-  n_k = normal[k];
-}
-
-///////////////////////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////////////////////
 // KDTREE DYN  KDTREE DYN  KDTREE DYN  KDTREE DYN  KDTREE DYN  KDTREE DYN  
KDTR
 
///////////////////////////////////////////////////////////////////////////////
 
///////////////////////////////////////////////////////////////////////////////

Modified: trunk/Model/Groups/KDTreeDyn.h
==============================================================================
--- trunk/Model/Groups/KDTreeDyn.h      (original)
+++ trunk/Model/Groups/KDTreeDyn.h      Thu Sep 28 13:11:43 2006
@@ -29,6 +29,7 @@
 #ifndef Manta_Groups_KDTreeDyn_h
 #define Manta_Groups_KDTreeDyn_h
 
+#include <Model/Groups/TriangleIW.h>
 #include <Model/Primitives/PrimitiveCommon.h>
 #include <Core/Geometry/Vector.h>
 #include <Core/Geometry/Ray.h>
@@ -42,85 +43,6 @@
 #include <vector>
 
 namespace Manta {
-
-  
/////////////////////////////////////////////////////////////////////////////
-  
/////////////////////////////////////////////////////////////////////////////
-  // TRIANGLE IW  TRIANGLE IW  TRIANGLE IW  TRIANGLE IW  TRIANGLE IW  
TRIANGLE 
-  
/////////////////////////////////////////////////////////////////////////////
-  
///////////////////////////////////////////////////////////////////////////// 
 
-  
-  class TriangleIW {
-  public:
-    
-    float n_u;
-    float n_v;
-    float n_d;
-    unsigned int k;
-    
-    float b_nu;
-    float b_nv;
-    float b_d;
-    float n_k;
-    
-    float c_nu;
-    float c_nv;
-    float c_d;
-
-    TriangleIW() {}
-    TriangleIW(const Vector& _p1, const Vector& _p2, const Vector& _p3);
-    
-    
///////////////////////////////////////////////////////////////////////////
-    // Single ray intersection.
-    // Returns true if the ray intersects the triangle, returns false 
otherwise.
-    // Mangles output parameter t in either case.
-    
///////////////////////////////////////////////////////////////////////////
-    inline
-    bool intersect_single_ray( // Input.
-                               const Vector &origin,
-                               const Vector &direction,
-                               const float &minT,
-
-                               // Output.
-                               float &t ) const {
-
-      int axis = k;
-      const int ku = (axis==2)?0:axis+1;
-      const int kv = (axis==0)?2:axis-1;
-            
-      const float org_k  = origin[axis];
-      const float org_ku = origin[ku];
-      const float org_kv = origin[kv];
-      const float f0     = n_d - (org_k + n_u * org_ku + n_v * org_kv);
-      
-      const float dir_k  = direction[axis];
-      const float dir_ku = direction[ku];
-      const float dir_kv = direction[kv];
-      
-      const float nd0 = n_u * dir_ku + n_v * dir_kv + dir_k;
-      const float nd  = 1.f/nd0;
-      
-      t = f0 * nd;
-
-      // plane test
-      if ( t < T_EPSILON || t > minT )
-        return false;
-      
-      const float hu = org_ku + t*dir_ku;
-      const float hv = org_kv + t*dir_kv;
-      const float lambda = b_d + hu*b_nu + hv * b_nv;
-      
-      // barycentric test
-      if ( lambda < 0.f )
-        return false;
-      
-      const float mue = c_d + hu * c_nu + hv * c_nv;
-      if ( mue < 0.f || mue + lambda > 1.f )
-        return false;
-
-      return true;
-    }
-  };  
-  
   
   
/////////////////////////////////////////////////////////////////////////////
   
/////////////////////////////////////////////////////////////////////////////
@@ -128,7 +50,6 @@
   
/////////////////////////////////////////////////////////////////////////////
   
/////////////////////////////////////////////////////////////////////////////
 
-  
   using std::vector;
   
   typedef unsigned int Integer;

Added: trunk/Model/Groups/TriangleIW.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/TriangleIW.cc    Thu Sep 28 13:11:43 2006
@@ -0,0 +1,89 @@
+/*
+  For more information, please see: http://software.sci.utah.edu
+
+  The MIT License
+
+  Copyright (c) 2005-2006
+  Scientific Computing and Imaging Institute, University of Utah
+
+  License for the specific language governing rights and limitations under
+  Permission is hereby granted, free of charge, to any person obtaining a
+  copy of this software and associated documentation files (the "Software"),
+  to deal in the Software without restriction, including without limitation
+  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+  and/or sell copies of the Software, and to permit persons to whom the
+  Software is furnished to do so, subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be included
+  in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+  DEALINGS IN THE SOFTWARE.
+*/
+
+#include <Model/Groups/TriangleIW.h>
+
+#include <Core/Geometry/Vector.h>
+#include <SCIRun/Core/Math/MiscMath.h>
+
+using namespace Manta;
+
+using SCIRun::Abs;
+
+TriangleIW::TriangleIW( const Vector& _p1,
+                        const Vector& _p2,
+                        const Vector& _p3 )
+{
+  Vector normal = Cross(_p2-_p1, _p3-_p1);
+  normal.normalize();
+
+  unsigned int n, u, v;
+
+  n = 0;
+  if ( Abs(normal[1]) > Abs(normal[n])) n = 1;
+  if ( Abs(normal[2]) > Abs(normal[n])) n = 2;
+
+  switch ( n ) {
+  case 0: u = 1; v = 2; break;
+  case 1: u = 2; v = 0; break;
+  default: u = 0; v = 1; break;
+  };
+
+  // copy to struct
+  k = n;
+  n_u = normal[u] / normal[k];
+  n_v = normal[v] / normal[k];
+  n_d = _p1[u] * n_u + _p1[v] * n_v + _p1[k];
+
+  float s;
+
+  c_nu = + _p2[v] - _p1[v];
+  c_nv = - _p2[u] + _p1[u];
+  c_d  = - (_p1[u] * c_nu + _p1[v] * c_nv);
+
+  s = 1.f / ( _p3[u] * c_nu + _p3[v] * c_nv + c_d );
+
+  c_nu *= s;
+  c_nv *= s;
+  c_d  *= s;
+
+  b_nu = + _p3[v] - _p1[v];
+  b_nv = - _p3[u] + _p1[u];
+  b_d  = - (_p1[u] * b_nu + _p1[v] * b_nv);
+
+  s = 1.f / ( _p2[u] * b_nu + _p2[v] * b_nv + b_d );
+
+  b_nu *= s;
+  b_nv *= s;
+  b_d  *= s;
+
+  n_k = normal[k];
+}
+
+
+

Added: trunk/Model/Groups/TriangleIW.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/TriangleIW.h     Thu Sep 28 13:11:43 2006
@@ -0,0 +1,144 @@
+
+#ifndef MODEL_GROUPS_TRIANGLEIW__H
+#define MODEL_GROUPS_TRIANGLEIW__H
+
+/*
+  For more information, please see: http://software.sci.utah.edu
+
+  The MIT License
+
+  Copyright (c) 2005-2006
+  Scientific Computing and Imaging Institute, University of Utah
+
+  License for the specific language governing rights and limitations under
+  Permission is hereby granted, free of charge, to any person obtaining a
+  copy of this software and associated documentation files (the "Software"),
+  to deal in the Software without restriction, including without limitation
+  the rights to use, copy, modify, merge, publish, distribute, sublicense,
+  and/or sell copies of the Software, and to permit persons to whom the
+  Software is furnished to do so, subject to the following conditions:
+
+  The above copyright notice and this permission notice shall be included
+  in all copies or substantial portions of the Software.
+
+  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+  DEALINGS IN THE SOFTWARE.
+*/
+
+#include <Core/Geometry/Vector.h>
+#include <Interface/Parameters.h>
+
+#include <Interface/Context.h>
+#include <Interface/RayPacket.h>
+#include <Core/Geometry/BBox.h>
+
+namespace Manta {
+  
+  class TriangleIW {
+  public:
+    
+    float n_u;
+    float n_v;
+    float n_d;
+    unsigned int k;
+    
+    float b_nu;
+    float b_nv;
+    float b_d;
+    float n_k;
+    
+    float c_nu;
+    float c_nv;
+    float c_d;
+
+    TriangleIW() {}
+    TriangleIW(const Vector& _p1, const Vector& _p2, const Vector& _p3);
+    
+    
///////////////////////////////////////////////////////////////////////////
+    // Single ray intersection.
+    // Returns true if the ray intersects the triangle, returns false 
otherwise.
+    // Mangles output parameter t in either case.
+    
///////////////////////////////////////////////////////////////////////////
+    inline
+    bool intersect_single_ray( // Input.
+                               const Vector &origin,
+                               const Vector &direction,
+                               const float &minT,
+
+                               // Output.
+                               float &t ) const {
+
+      int axis = k;
+      const int ku = (axis==2)?0:axis+1;
+      const int kv = (axis==0)?2:axis-1;
+            
+      const float org_k  = origin[axis];
+      const float org_ku = origin[ku];
+      const float org_kv = origin[kv];
+      const float f0     = n_d - (org_k + n_u * org_ku + n_v * org_kv);
+      
+      const float dir_k  = direction[axis];
+      const float dir_ku = direction[ku];
+      const float dir_kv = direction[kv];
+      
+      const float nd0 = n_u * dir_ku + n_v * dir_kv + dir_k;
+      const float nd  = 1.f/nd0;
+      
+      t = f0 * nd;
+
+      // plane test
+      if ( t < T_EPSILON || t > minT )
+        return false;
+      
+      const float hu = org_ku + t*dir_ku;
+      const float hv = org_kv + t*dir_kv;
+      const float lambda = b_d + hu*b_nu + hv * b_nv;
+      
+      // barycentric test
+      if ( lambda < 0.f )
+        return false;
+      
+      const float mue = c_d + hu * c_nu + hv * c_nv;
+      if ( mue < 0.f || mue + lambda > 1.f )
+        return false;
+
+      return true;
+    }
+  };
+
+  class TriangleIW_Object : public Object,  public TriangleIW {
+  public:
+    // Constructors.
+    TriangleIW_Object() { };
+    TriangleIW_Object(const Vector& _p1, const Vector& _p2, const Vector& 
_p3)
+      : TriangleIW( _p1, _p2, _p3 )
+    {
+      bounds.extendByPoint( _p1 );
+      bounds.extendByPoint( _p2 );
+      bounds.extendByPoint( _p3 );
+    }
+    
+    
///////////////////////////////////////////////////////////////////////////
+    // Object Interface.
+    
///////////////////////////////////////////////////////////////////////////
+    virtual void preprocess   (const PreprocessContext& context) { };
+    virtual void computeBounds(const PreprocessContext&, BBox& bbox) const { 
bbox.extendByBox( bounds ); }
+    
+    virtual void intersect    (const RenderContext&, RayPacket& rays) const {
+      
+      // Currently unimplemented. intersect_single_ray is called explicitly
+      // from acceleration structure code.
+    }
+
+  private:
+    BBox bounds;
+  };
+};
+
+#endif
+

Modified: trunk/SwigInterface/KDTreeDynFrame.py
==============================================================================
--- trunk/SwigInterface/KDTreeDynFrame.py       (original)
+++ trunk/SwigInterface/KDTreeDynFrame.py       Thu Sep 28 13:11:43 2006
@@ -62,7 +62,7 @@
         self.engine.addTransaction( "update plane",
                                     manta_new( createMantaTransaction( 
self.UpdateFrame,
                                                                         () 
)))        
-
+        
     def UpdateFrame(self):
         
         # Get the current root node.
@@ -135,6 +135,9 @@
     ## Event Callbacks
     
###########################################################################
     
###########################################################################
+
+    def OnShow(self, event):
+        self.Show()
 
     def OnResetRoot(self, event):
 

Modified: trunk/SwigInterface/kdtreedyn.py
==============================================================================
--- trunk/SwigInterface/kdtreedyn.py    (original)
+++ trunk/SwigInterface/kdtreedyn.py    Thu Sep 28 13:11:43 2006
@@ -15,8 +15,9 @@
 
 # Remember to declare these are "global" when defining them.
 file_name = ""
-leaf_threshold = 4
-depth_threshold = 16
+leaf_threshold       = 8
+depth_threshold      = 64
+cost_threshold       = 1.0
 
 
###############################################################################
 
###############################################################################
@@ -26,7 +27,7 @@
 
 def initialize_scene( frame, engine ):
 
-    global leaf_threshold, depth_threshold
+
 
     
###########################################################################
     # Create basic scene with background.
@@ -34,6 +35,8 @@
     scene.getRenderParameters().maxDepth = 5
     
scene.setBackground(manta_new(ConstantBackground(ColorDB.getNamedColor("FloralWhite"))))
 
+    engine.selectShadowAlgorithm( "noshadows" )
+
     # Create lights.
     lights = manta_new(LightSet())
     lights.add(manta_new(PointLight(Vector(0,5,8), 
Color(RGBColor(.6,.1,.1)))))
@@ -56,7 +59,10 @@
 
     
###########################################################################
     # Build the kdtree.
-    builder = Approx_cc( leaf_threshold, depth_threshold )
+    global leaf_threshold, depth_threshold, cost_threshold
+    builder = Build_OnLogn_cc( leaf_threshold,
+                               depth_threshold,
+                               cost_threshold )
 
     # Add the kdtree to the scene.
     world = manta_new( Group() )
@@ -79,11 +85,14 @@
     engine.setScene( scene )
 
     # Create the explorer.
-    # global explorer
+    global explorer
     explorer = KDTreeDynFrame( frame, frame.engine, kdtree )
-    explorer.Show()
-
 
+    # Add the explorer to the main frame menu.
+    kd_menu = wx.Menu()
+    frame.Bind(wx.EVT_MENU, explorer.OnShow,
+               kd_menu.Append(wx.NewId(), "Explorer"))
+    frame.menuBar.Append( kd_menu, "&KDTreeDyn" );
 
 
 
###############################################################################
@@ -98,11 +107,15 @@
     print "-n --np=<threads>"
     print "-f --file=<file name>"
     print "   --depth_threshold=<n>"
+    print "      Create a leaf if depth is equal to n"
     print "   --leaf_threshold=<n>"
+    print "      Create a leaf if the total number of faces is less than n."
+    print "   --cost_threshold=<f>"
+    print "      Make a leaf if the (cost of splitting)*f > (leaf cost)."
 
 def main():
 
-    global file_name, leaf_threshold, depth_threshold
+    global file_name, leaf_threshold, depth_threshold, refinement_threshold, 
cost_threshold
 
     # Default options.
     num_workers = 1
@@ -115,28 +128,30 @@
                                    ["np=",
                                     "file=",
                                     "leaf_threshold=",
-                                    "depth_threshold="] )
+                                    "depth_threshold=",
+                                    "cost_threshold=" ] )
 
     except getopt.GetoptError:
         usage()
         sys.exit(2)
         
     for o, a in opts:
-        if o in ("-n", "--np"):
-            try:
+        try:
+            if o in ("-n", "--np"):
                 num_workers = int(a)
-            except ValueError:
-                usage()
-                sys.exit(2)
-                
-        if o in ("-f", "--file"):
-            file_name = a;
-
-        if o in ("--leaf_threshold"):
-            leaf_threshold = int(a)
+            elif o in ("-f", "--file"):
+                file_name = a;      
+            elif o in ("--leaf_threshold"):
+                leaf_threshold = int(a)
+            elif o in ("--depth_threshold"):
+                depth_threshold = int(a)
+            elif o in ("--cost_threshold"):
+                cost_threshold = float(a)
+        except ValueError:
+            print "Error parsing " + a + " as a number."
+            usage()
+            sys.exit(2)
 
-        if o in ("--depth_threshold"):
-            depth_threshold = int(a)
 
         # Add additional command line args here.
         

Modified: trunk/SwigInterface/wxManta.py
==============================================================================
--- trunk/SwigInterface/wxManta.py      (original)
+++ trunk/SwigInterface/wxManta.py      Thu Sep 28 13:11:43 2006
@@ -632,8 +632,6 @@
     
###########################################################################
     def OnLeftDClick( self, event ):
 
-        print "Double click"
-
         # Compute pixel coordinates.
         mouse_pos = event.GetPosition()
         window_size = self.canvas.GetSize()
@@ -651,18 +649,24 @@
         data = RayPacketData()
         rays = RayPacket( data, RayPacket.UnknownShape, 0, 1, 0, 0 )
 
-        # Shoot one ray.
+        # Pass a color to store the result.
         color = Color(RGBColor())
 
+        # Shoot a single ray.
         self.engine.shootOneRay( color, rays, xpos, ypos, channel_index )
 
         if (rays.wasHit(0)):
-            position = rays.getHitPosition(0)
-
-            print "Location = (%s, %s, %s)" % (position.x(), position.y(), 
position.z())
-        else:
-            print "Missed"
+            lookat = rays.getHitPosition(0)
 
+            # print "Location = (%s, %s, %s)" % (lookat.x(), lookat.y(), 
lookat.z())
+            
+            # Change the camera look at position.
+            camera = self.engine.getCamera( self.channelID )
+            camera.reset( Vector(camera.getPosition()),
+                          Vector(camera.getUp()),
+                          Vector(lookat)
+                          
+                          )
         
     
###########################################################################
     ## mouseRotate




  • [MANTA] r1193 - in trunk: Model/Groups SwigInterface, abe, 09/28/2006

Archive powered by MHonArc 2.6.16.

Top of page