<html><head><style type="text/css"><!-- DIV {margin:0px;} --></style></head><body><div style="font-family:arial, helvetica, sans-serif;font-size:10pt;color:#000000;">Hi All,<div><br></div><div>While tracking some lost particles I found a bug in MBOrientedBox::intersect_ray. It was due to subtracting a squared distance. More specifically, the fact that a squared distance becomes larger if the distance is greater than one, but smaller if the distance is less than one (assuming distance is always positive). With Jason's help I tried to make the tolerances used in the test more consistent. Below is a diff of the fix. The fix reduces the number of lost particles by ~50%.</div><div><br></div><div>Are there any changes/suggestions to this diff?</div><div><br></div><div>Brandon</div><div><br></div><div><div>Index: MBOrientedBox.hpp</div><div>===================================================================</div><div>--- MBOrientedBox.hpp (revision
3533)</div><div>+++ MBOrientedBox.hpp (working copy)</div><div>@@ -53,8 +53,8 @@</div><div> </div><div> inline double inner_radius() const; //!< radius of inscribed sphere</div><div> inline double outer_radius() const; //!< radius of circumscribed sphere</div><div>- inline double outer_radius_squared() const; //!< square of radius of circumsphere</div><div>- inline double inner_radius_squared() const; //!< square of radius if inscribed sphere</div><div>+ inline double outer_radius_squared(const double reps) const; //!< square of (radius+at least epsilon) of circumsphere</div><div>+ inline double inner_radius_squared(const double reps) const; //!< square of (radius-epsilon) of inscribed sphere</div><div> inline double volume() const; //!< volume of box</div><div> inline MBCartVect dimensions() const;
//!< number of dimensions for which box is not flat</div><div> inline double area() const; //!< largest side area</div><div>@@ -176,24 +176,30 @@</div><div> #endif</div><div> }</div><div> </div><div>-double MBOrientedBox::outer_radius_squared() const</div><div>+// Add at least epsilon to the radius, before squaring it.</div><div>+double MBOrientedBox::outer_radius_squared(const double reps) const</div><div> {</div><div> #if MB_ORIENTED_BOX_OUTER_RADIUS</div><div>- return radius * radius;</div><div>+ return (radius+reps)*(radius+reps);</div><div> #elif MB_ORIENTED_BOX_UNIT_VECTORS</div><div>- return length % length;</div><div>+ MBCartVect tmp(length[0]+reps,length[1]+reps,length[2]+reps);</div><div>+ return tmp % tmp;</div><div> #else</div><div>- const MBCartVect half_diag = axis[0] +
axis[1] + axis[2];</div><div>+ MBCartVect half_diag = axis[0] + axis[1] + axis[2];</div><div>+ half_diag += MBCartVect(reps,reps,reps);</div><div> return half_diag % half_diag;</div><div> #endif</div><div> }</div><div> </div><div>-double MBOrientedBox::inner_radius_squared() const</div><div>+// Subtract epsilon from the length of the shortest axis, before squaring it.</div><div>+double MBOrientedBox::inner_radius_squared(const double reps) const</div><div> {</div><div> #if MB_ORIENTED_BOX_UNIT_VECTORS</div><div>- return length[0] * length[0];</div><div>+ return (length[0]-reps) * (length[0]-reps);</div><div> #else</div><div>- return axis[0] % axis[0];</div><div>+ MBCartVect tmp = axis[0];</div><div>+ tmp -= MBCartVect(reps,reps,reps);</div><div>+ return (tmp % tmp);</div><div> #endif</div><div> }</div><div> </div></div><div><br></div><div><div>Index:
MBOrientedBox.cpp</div><div>===================================================================</div><div>--- MBOrientedBox.cpp (revision 3533)</div><div>+++ MBOrientedBox.cpp (working copy)</div><div>@@ -477,16 +477,27 @@</div><div> const MBCartVect cx = center - b;</div><div> double dist_s = cx % m;</div><div> double dist_sq = cx % cx - (dist_s*dist_s);</div><div>- double max_diagsq = outer_radius_squared();</div><div>+ double max_diagsq = outer_radius_squared(reps);</div><div> </div><div> // if greater than the longest diagonal, we don't hit</div><div>- if (dist_sq > max_diagsq+reps)</div><div>+ if (dist_sq > max_diagsq)</div><div> return false;</div><div>- if (len && dist_s - max_diagsq > *len)</div><div>- return false;</div><div>- </div><div>+</div><div>+ // If the closest
possible hit is farther than len, we don't want the hit.</div><div>+ // Problem: the previous method was wrong because max_diagsq will be greater</div><div>+ // then max_diag is max_diag>1 but less than max_diag if max_diag<1.</div><div>+ if (len) {</div><div>+ if (max_diagsq>1) {</div><div>+ if (dist_s - max_diagsq > *len)</div><div>+ return false;</div><div>+ } else {</div><div>+ if (dist_s - 1 > *len)</div><div>+ return false;</div><div>+ }</div><div>+ }</div><div>+</div><div> // if smaller than shortest diagonal, we do hit</div><div>- if (dist_sq < inner_radius_squared() - reps && dist_s >= 0.0)</div><div>+ if (dist_sq < inner_radius_squared(reps) && dist_s >= 0.0)</div><div>
return true;</div><div> </div><div> // get transpose of axes</div></div></div><br>
</body></html>