amanatidesWoo3D.h

Go to the documentation of this file.

#ifndef _DLR_NUMERIC_AMANATIDESWOO3D_H_
#define _DLR_NUMERIC_AMANATIDESWOO3D_H_

#include <iostream>
#include <limits>
#include <dlrCommon/exception.h>
#include <dlrNumeric/vector3D.h>
#include <dlrNumeric/transform3D.h>
#include <dlrNumeric/amanatidesWoo3DIterator.h>

namespace dlr {

  namespace numeric {
    
    template <class ARRAY3D>
    class AmanatidesWoo3D {
    public:
      /* ================== Public typedefs ================== */

      typedef AmanatidesWoo3DIterator<ARRAY3D> iterator;

      // Iteration over const arrays will be done using a different
      // class.
      // typedef AmanatidesWoo3DConstIterator<ARRAY3D> const_iterator;
    
      /* ================== Public methods ================== */

      AmanatidesWoo3D(ARRAY3D& data,
                      const Transform3D& voxelTworld,
                      const Vector3D& rayOrigin,
                      const Vector3D& rayDirection,
                      bool downstreamOnly=true);

    
      AmanatidesWoo3D(const AmanatidesWoo3D& source);

      ~AmanatidesWoo3D();

      iterator
      begin();

      iterator
      end();

      ARRAY3D&
      getData() {return m_data;}

      inline bool
      validIntersection();
    
      AmanatidesWoo3D&
      operator=(const AmanatidesWoo3D& source);
  
    private:

      std::pair<double, double>
      findEntryAndExitPoints(const Vector3D& rayOriginVoxel,
                             const Vector3D& rayDirectionVoxel,
                             const ARRAY3D& m_data);

      double
      findIntersection(const Vector3D& rayOrigin,
                       const Vector3D& rayDirection,
                       const Vector3D& bVector,
                       double cConstant,
                       double defaultValue);

      ARRAY3D& m_data;

      int m_initialU;

      int m_initialV;

      int m_initialW;
    
      int m_stepU;

      int m_stepV;
    
      int m_stepW;
    
      double m_tDeltaU;

      double m_tDeltaV;
    
      double m_tDeltaW;
    
      double m_tMaxU;

      double m_tMaxV;

      double m_tMaxW;

      double m_tStart;

      bool m_validIntersection;
    };

  } // namespace numeric

} // namespace dlr


/* ======= Declarations to maintain compatibility with legacy code. ======= */

namespace dlr {

  using numeric::AmanatidesWoo3D;

} // namespace dlr


/* ========================================================= */
/* Implementation follows                                    */
/* ========================================================= */

namespace dlr {

  namespace numeric {
    
    template <class ARRAY3D>
    AmanatidesWoo3D<ARRAY3D>::
    AmanatidesWoo3D(ARRAY3D& data,
                    const Transform3D& voxelTworld,
                    const Vector3D& rayOrigin,
                    const Vector3D& rayDirection,
                    bool downstreamOnly)
      : m_data(data),
        m_initialU(-1),  // Initialize to illegal values.  These will be 
        m_initialV(-1),  // replaced with legal ones if there's a valid
        m_initialW(-1),  // intersection between the ray and the voxel
        // array.
        m_stepU(),
        m_stepV(),
        m_stepW(),
        m_tDeltaU(),
        m_tDeltaV(),
        m_tDeltaW(),
        m_tMaxU(),
        m_tMaxV(),
        m_tMaxW(),
        m_tStart(),
        m_validIntersection(true)    
    {
      // First convert everything into voxel coordinates.
      Vector3D rayOriginVoxel = voxelTworld * rayOrigin;
      Vector3D rayDirectionVoxel =
        (voxelTworld * (rayOrigin + rayDirection)) - rayOriginVoxel;

      // Now find points entry and exit from the CT volume.  These are
      // expressed in as parameter values tEntry and tExit such that
      // (rayOriginVoxel + tEntry * rayDirectionVoxel) is the entry
      // point and (rayOriginVoxel + tExit * rayDirectionVoxel) is the
      // exit point.
      std::pair<double, double> tEntry_tExit =
        this->findEntryAndExitPoints(rayOriginVoxel, rayDirectionVoxel, m_data);

      // Sometimes we want to disallow any voxels which are "behind" the
      // ray origin.  That is, sometimes we want to only trace those
      // parts of the ray which correspond to positive values of the
      // parameter t.
      if(downstreamOnly && (tEntry_tExit.first < 0.0)) {
        tEntry_tExit.first = 0.0;
      }

      // Sometimes the there's no intersection with the voxel array.  In
      // this case, tExit will be less than tEntry, and we don't bother
      // doing any more calculation.
      if(tEntry_tExit.second <= tEntry_tExit.first) {
        m_validIntersection = false;
        return;
      }

      // Now make the point of entry explicit.
      Vector3D entryPoint =
        rayOriginVoxel + tEntry_tExit.first * rayDirectionVoxel;

      // Correct for rounding error which sometimes makes entryPoint
      // have values like -1.0e-14.
      if(entryPoint.x() < 0.0) {
        entryPoint.x() = 0.0;
      }
      if(entryPoint.y() < 0.0) {
        entryPoint.y() = 0.0;
      }
      if(entryPoint.z() < 0.0) {
        entryPoint.z() = 0.0;
      }
    
      // Finally, assign the variables described in the Amanatides' and
      // Woo's paper.

      // Since we've already converted to voxel coords, finding the
      // first voxel on the path is trivial.
      m_initialU = static_cast<int>(entryPoint.x());
      m_initialV = static_cast<int>(entryPoint.y());
      m_initialW = static_cast<int>(entryPoint.z());

      // There's an similar case to the one handled by the rounding
      // error test above.  Member function findEntryAndExitPoints()
      // uses m_data.shape to define the maximum U, V, and W
      // coordinates, respectively.  This means that it's very possible
      // for m_initialU to be equal to m_data.columns(), which will
      // cause an out-of-bounds index. Correct for this now.
      if(m_initialU == static_cast<int>(m_data.shape()[2])) {
        --m_initialU;
      }
      if(m_initialV == static_cast<int>(m_data.shape()[1])) {
        --m_initialV;
      }
      if(m_initialW == static_cast<int>(m_data.shape()[0])) {
        --m_initialW;
      }

      // Sanity check.
      if((m_initialU >= static_cast<int>(m_data.shape()[2]))
         || (m_initialV >= static_cast<int>(m_data.shape()[1]))
         || (m_initialW >= static_cast<int>(m_data.shape()[0]))) {
        DLR_THROW3(LogicException,
                   "AmanatidesWoo3D::AmanatidesWoo3D(ARRAY3D&, ...)",
                   "Illegal value for m_initialU, m_initialV, or m_initialW.");
      }
    
      // m_tStart is just the same as tEntry.
      m_tStart = tEntry_tExit.first;
    
      // The remaining variables depend on whether U & V will be
      // increasing or decreasing as we travel along the ray, so we need
      // if clauses.  Please see the declaration for documentation on
      // what each of these member variables means.
      if(rayDirectionVoxel.x() > 0.0) {
        m_stepU = 1;
        m_tDeltaU = 1.0 / rayDirectionVoxel.x();
        m_tMaxU = m_tStart + (((m_initialU + 1) - entryPoint.x())
                              / rayDirectionVoxel.x());
      } else if(rayDirectionVoxel.x() < 0.0) {
        m_stepU = -1;
        m_tDeltaU = -(1.0 / rayDirectionVoxel.x());
        m_tMaxU = m_tStart + ((m_initialU - entryPoint.x())
                              / rayDirectionVoxel.x());
      } else { // rayDirectionVoxel.x() == 0.0;
        m_stepU = 0;
        m_tDeltaU = std::numeric_limits<double>::max();
        m_tMaxU = std::numeric_limits<double>::max();
      }
      if(rayDirectionVoxel.y() > 0.0) {
        m_stepV = 1;
        m_tDeltaV = 1.0 / rayDirectionVoxel.y();
        m_tMaxV = m_tStart + (((m_initialV + 1) - entryPoint.y())
                              / rayDirectionVoxel.y());
      } else if(rayDirectionVoxel.y() < 0.0) {
        m_stepV = -1;
        m_tDeltaV = -(1.0 / rayDirectionVoxel.y());
        m_tMaxV = m_tStart + ((m_initialV - entryPoint.y())
                              / rayDirectionVoxel.y());
      } else { // rayDirectionVoxel.y() == 0.0;
        m_stepV = 0;
        m_tDeltaV = std::numeric_limits<double>::max();
        m_tMaxV = std::numeric_limits<double>::max();
      }
      if(rayDirectionVoxel.z() > 0.0) {
        m_stepW = 1;
        m_tDeltaW = 1.0 / rayDirectionVoxel.z();
        m_tMaxW = m_tStart + (((m_initialW + 1) - entryPoint.z())
                              / rayDirectionVoxel.z());
      } else if(rayDirectionVoxel.z() < 0.0) {
        m_stepW = -1;
        m_tDeltaW = -(1.0 / rayDirectionVoxel.z());
        m_tMaxW = m_tStart + ((m_initialW - entryPoint.z())
                              / rayDirectionVoxel.z());
      } else { // rayDirectionVoxel.z() == 0.0;
        m_stepW = 0;
        m_tDeltaW = std::numeric_limits<double>::max();
        m_tMaxW = std::numeric_limits<double>::max();
      }
    }
    
    template <class ARRAY3D>
    AmanatidesWoo3D<ARRAY3D>::
    AmanatidesWoo3D(const AmanatidesWoo3D& source)
      : m_data(source.m_data),
        m_initialU(source.m_initialU),
        m_initialV(source.m_initialV),
        m_initialW(source.m_initialW),
        m_stepU(source.m_stepU),
        m_stepV(source.m_stepV),
        m_stepW(source.m_stepW),
        m_tDeltaU(source.m_tDeltaU),
        m_tDeltaV(source.m_tDeltaV),
        m_tDeltaW(source.m_tDeltaW),
        m_tMaxU(source.m_tMaxU),
        m_tMaxV(source.m_tMaxV),
        m_tMaxW(source.m_tMaxW),
        m_tStart(source.m_tStart),
        m_validIntersection(source.m_validIntersection)
    {
      // Empty
    }

    template <class ARRAY3D>
    AmanatidesWoo3D<ARRAY3D>::
    ~AmanatidesWoo3D() {};

    template <class ARRAY3D>
    typename AmanatidesWoo3D<ARRAY3D>::iterator
    AmanatidesWoo3D<ARRAY3D>::
    begin()
    {
      return AmanatidesWoo3DIterator<ARRAY3D>(
        m_data, m_initialU, m_initialV, m_initialW, m_stepU, m_stepV, m_stepW,
        m_tMaxU, m_tMaxV, m_tMaxW, m_tDeltaU, m_tDeltaV, m_tDeltaW, m_tStart);
    }

    template <class ARRAY3D>
    typename AmanatidesWoo3D<ARRAY3D>::iterator
    AmanatidesWoo3D<ARRAY3D>::
    end()
    {
      // AmanatidesWoo3DIterator<ARRAY3D>::operator==(...) considers all
      // iterators with illegal voxel coordinates to be equal, so all we
      // have to do here is return an iterator which references an
      // illegal voxel.
      return AmanatidesWoo3DIterator<ARRAY3D>(
        m_data, -1, -1, -1, m_stepU, m_stepV, m_stepW,
        m_tMaxU, m_tMaxV, m_tMaxW, m_tDeltaU, m_tDeltaV, m_tDeltaW, m_tStart);
    }

    template <class ARRAY3D>
    inline bool
    AmanatidesWoo3D<ARRAY3D>::
    validIntersection()
    {
      return m_validIntersection;
    }

    template <class ARRAY3D>
    AmanatidesWoo3D<ARRAY3D>&
    AmanatidesWoo3D<ARRAY3D>::
    operator=(const AmanatidesWoo3D& source)
    {
      m_data = source.m_data;
      m_initialU = source.m_initialU;
      m_initialV = source.m_initialV;
      m_initialW = source.m_initialW;
      m_stepU = source.m_stepU;
      m_stepV = source.m_stepV;
      m_stepW = source.m_stepW;
      m_tDeltaU = source.m_tDeltaU;
      m_tDeltaV = source.m_tDeltaV;
      m_tDeltaW = source.m_tDeltaW;
      m_tMaxU = source.m_tMaxU;
      m_tMaxV = source.m_tMaxV;
      m_tMaxW = source.m_tMaxW;
      m_tStart = source.m_tStart;
      m_validIntersection = source.m_validIntersection;
      return *this;
    }
    
    template <class ARRAY3D>
    std::pair<double, double>
    AmanatidesWoo3D<ARRAY3D>::
    findEntryAndExitPoints(const Vector3D& rayOriginVoxel,
                           const Vector3D& rayDirectionVoxel,
                           const ARRAY3D& m_data)
    {
      // First find intersection with each boundary line.

      // ... Find the intersection with the plane U = 0, or else a
      // really small number if rayDirection is parallel to the U axis.
      double tIntersectU0 = findIntersection(
        rayOriginVoxel, rayDirectionVoxel, Vector3D(1.0, 0.0, 0.0),
        0.0, std::numeric_limits<double>::min());
    
      // ... Find the intersection with the plane U = m_data.shape()[2],
      // or else a really big number if rayDirection is parallel to the
      // U axis.
      double tIntersectU1 = findIntersection(
        rayOriginVoxel, rayDirectionVoxel, Vector3D(1.0, 0.0, 0.0),
        m_data.shape()[2], std::numeric_limits<double>::max());
    
      // ... Find the intersection with the plane V = 0, or else a
      // really small number if rayDirection is parallel to the V axis.
      double tIntersectV0 = findIntersection(
        rayOriginVoxel, rayDirectionVoxel, Vector3D(0.0, 1.0, 0.0),
        0.0, std::numeric_limits<double>::min());
    
      // ... Find the intersection with the plane V = m_data.shape()[1],
      // or else a really big number if rayDirection is parallel to the
      // V axis.
      double tIntersectV1 = findIntersection(
        rayOriginVoxel, rayDirectionVoxel, Vector3D(0.0, 1.0, 0.0),
        m_data.shape()[1], std::numeric_limits<double>::max());
    
      // ... Find the intersection with the plane W = 0, or else a
      // really small number if rayDirection is parallel to the W axis.
      double tIntersectW0 = findIntersection(
        rayOriginVoxel, rayDirectionVoxel, Vector3D(0.0, 0.0, 1.0),
        0.0, std::numeric_limits<double>::min());
    
      // ... Find the intersection with the plane W = m_data.shape()[0],
      // or else a really big number if rayDirection is parallel to the
      // W axis.
      double tIntersectW1 = findIntersection(
        rayOriginVoxel, rayDirectionVoxel, Vector3D(0.0, 0.0, 1.0),
        m_data.shape()[0], std::numeric_limits<double>::max());

      // Now find the closer and farther of each pair of intersections.
      double tMinU = std::min(tIntersectU0, tIntersectU1);
      double tMinV = std::min(tIntersectV0, tIntersectV1);
      double tMinW = std::min(tIntersectW0, tIntersectW1);
      double tMaxU = std::max(tIntersectU0, tIntersectU1);
      double tMaxV = std::max(tIntersectV0, tIntersectV1);
      double tMaxW = std::max(tIntersectW0, tIntersectW1);

      // Compute closest point which could possibly intersect with the volume.
      double tEntry = std::max(tMinU, std::max(tMinV, tMinW));
      // Compute farthest point which could possibly intersect with the volume.
      double tExit = std::min(tMaxU, std::min(tMaxV, tMaxW));

      // Return our findings.
      return std::make_pair(tEntry, tExit);
    }

    template <class ARRAY3D>
    double
    AmanatidesWoo3D<ARRAY3D>::
    findIntersection(const Vector3D& rayOrigin,
                     const Vector3D& rayDirection,
                     const Vector3D& bVector,
                     double cConstant,
                     double defaultValue)
    {
      // bVector and cConstant describe the desired plane:
      //   dot(x, bVector) = cConstant.
      // Now, x is constrained to lie on the specified line:
      //   x = rayOrigin + alpha * rayDirection.
      // Substituting, we have:
      //   dot(rayOrigin + alpha * rayDirection, bVector) = cConstant.
      // Which we rewrite:
      //   dot(rayOrigin, bVector) + alpha * dot(rayDirection, bVector)
      //     = cConstant.
      // Solving for alpha, we have:
      //   alpha = (cConstant - dot(rayOrigin, bVector))
      //            / dot(rayDirection, bVector).
      double denominator = dot(rayDirection, bVector);
      if(denominator == 0.0) {
        return defaultValue;
      }
      // else
      return (cConstant - dot(rayOrigin, bVector)) / denominator;
    }      

  } // namespace numeric

} // namespace dlr

#endif // #ifndef _DLR_NUMERIC_AMANATIDESWOO3D_H_