amanatidesWoo2D.h

Go to the documentation of this file.

#ifndef _DLR_NUMERIC_AMANATIDESWOO2D_H_
#define _DLR_NUMERIC_AMANATIDESWOO2D_H_

#include <iostream>
#include <limits>
#include <dlrCommon/exception.h>
#include <dlrNumeric/vector2D.h>
#include <dlrNumeric/transform2D.h>
#include <dlrNumeric/utilities.h>
#include <dlrNumeric/amanatidesWoo2DIterator.h>

namespace dlr {

  namespace numeric {
    
    template <class ARRAY2D>
    class AmanatidesWoo2D {
    public:
      /* ================== Public typedefs ================== */

      typedef AmanatidesWoo2DIterator<ARRAY2D> iterator;

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

      AmanatidesWoo2D(ARRAY2D& data,
                      const Transform2D& pixelTworld,
                      const Vector2D& rayOrigin,
                      const Vector2D& rayDirection,
                      bool downstreamOnly=true);

    
      AmanatidesWoo2D(const AmanatidesWoo2D& source);

      ~AmanatidesWoo2D();

      iterator
      begin();

      iterator
      end();

      ARRAY2D&
      getData() {return m_data;}
    
      inline bool
      validIntersection();

      AmanatidesWoo2D&
      operator=(const AmanatidesWoo2D& source);
  
    private:

      std::pair<double, double>
      findEntryAndExitPoints(const Vector2D& rayOriginPixel,
                             const Vector2D& rayDirectionPixel,
                             const ARRAY2D& data);

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

      ARRAY2D& m_data;

      int m_initialU;

      int m_initialV;

      int m_stepU;

      int m_stepV;
    
      double m_tDeltaU;

      double m_tDeltaV;
    
      double m_tMaxU;

      double m_tMaxV;

      double m_tStart;

      bool m_validIntersection;
    };

  } // namespace numeric

} // namespace dlr


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

namespace dlr {

  using numeric::AmanatidesWoo2D;

} // namespace dlr


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

namespace dlr {

  namespace numeric {
    
    template <class ARRAY2D>
    AmanatidesWoo2D<ARRAY2D>::
    AmanatidesWoo2D(ARRAY2D& data,
                    const Transform2D& pixelTworld,
                    const Vector2D& rayOrigin,
                    const Vector2D& 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
        // intersection between the ray and the pixel
        // array.
        m_stepU(),
        m_stepV(),
        m_tDeltaU(),
        m_tDeltaV(),
        m_tMaxU(),
        m_tMaxV(),
        m_tStart(),
        m_validIntersection(true)    
    {
      // First convert everything into pixel coordinates.
      Vector2D rayOriginPixel = pixelTworld * rayOrigin;
      Vector2D rayDirectionPixel =
        (pixelTworld * (rayOrigin + rayDirection)) - rayOriginPixel;

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

      // Sometimes we want to disallow any pixels 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 pixel 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.
      Vector2D entryPoint =
        rayOriginPixel + tEntry_tExit.first * rayDirectionPixel;

      // 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;
      }
    
      // Finally, assign the variables described in the Amanatides' and
      // Woo's paper.

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

      // There's an similar case to the one handled by the rounding
      // error test above.  Member function findEntryAndExitPoints()
      // uses m_data.columns() and m_data.rows() to define the maximum U
      // and V 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.columns())) {
        --m_initialU;
      }
      if(m_initialV == static_cast<int>(m_data.rows())) {
        --m_initialV;
      }

      // Sanity check.
      if((m_initialU >= static_cast<int>(m_data.columns()))
         || (m_initialV >= static_cast<int>(m_data.rows()))) {
        DLR_THROW3(LogicException,
                   "AmanatidesWoo2D::AmanatidesWoo2D(ARRAY2D&, ...)",
                   "Illegal value for m_initialU or m_initialV.");
      }
    
      // 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(rayDirectionPixel.x() > 0.0) {
        m_stepU = 1;
        m_tDeltaU = 1.0 / rayDirectionPixel.x();
        m_tMaxU = m_tStart + (((m_initialU + 1) - entryPoint.x())
                              / rayDirectionPixel.x());
      } else if(rayDirectionPixel.x() < 0.0) {
        m_stepU = -1;
        m_tDeltaU = -(1.0 / rayDirectionPixel.x());
        m_tMaxU = m_tStart + ((m_initialU - entryPoint.x())
                              / rayDirectionPixel.x());
      } else { // rayDirectionPixel.x() == 0.0;
        m_stepU = 0;
        m_tDeltaU = std::numeric_limits<double>::max();
        m_tMaxU = std::numeric_limits<double>::max();
      }
      if(rayDirectionPixel.y() > 0.0) {
        m_stepV = 1;
        m_tDeltaV = 1.0 / rayDirectionPixel.y();
        m_tMaxV = m_tStart + (((m_initialV + 1) - entryPoint.y())
                              / rayDirectionPixel.y());
      } else if(rayDirectionPixel.y() < 0.0) {
        m_stepV = -1;
        m_tDeltaV = -(1.0 / rayDirectionPixel.y());
        m_tMaxV = m_tStart + ((m_initialV - entryPoint.y())
                              / rayDirectionPixel.y());
      } else { // rayDirectionPixel.y() == 0.0;
        m_stepV = 0;
        m_tDeltaV = std::numeric_limits<double>::max();
        m_tMaxV = std::numeric_limits<double>::max();
      }
    }


    // The copy constructor deep copies its argument.
    template <class ARRAY2D>
    AmanatidesWoo2D<ARRAY2D>::
    AmanatidesWoo2D(const AmanatidesWoo2D& source)
      : m_data(source.m_data),
        m_initialU(source.m_initialU),
        m_initialV(source.m_initialV),
        m_stepU(source.m_stepU),
        m_stepV(source.m_stepV),
        m_tDeltaU(source.m_tDeltaU),
        m_tDeltaV(source.m_tDeltaV),
        m_tMaxU(source.m_tMaxU),
        m_tMaxV(source.m_tMaxV),
        m_tStart(source.m_tStart),
        m_validIntersection(source.m_validIntersection)
    {
      // Empty
    }

    template <class ARRAY2D>
    AmanatidesWoo2D<ARRAY2D>::
    ~AmanatidesWoo2D() {};

    template <class ARRAY2D>
    typename AmanatidesWoo2D<ARRAY2D>::iterator
    AmanatidesWoo2D<ARRAY2D>::
    begin()
    {
      return AmanatidesWoo2DIterator<ARRAY2D>(
        m_data, m_initialU, m_initialV, m_stepU, m_stepV, m_tMaxU, m_tMaxV,
        m_tDeltaU, m_tDeltaV, m_tStart);
    }

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

    template <class ARRAY2D>
    inline bool
    AmanatidesWoo2D<ARRAY2D>::
    validIntersection()
    {
      return m_validIntersection;
    }

    template <class ARRAY2D>
    AmanatidesWoo2D<ARRAY2D>&
    AmanatidesWoo2D<ARRAY2D>::
    operator=(const AmanatidesWoo2D& source)
    {
      m_data = source.m_data;
      m_initialU = source.m_initialU;
      m_initialV = source.m_initialV;
      m_stepU = source.m_stepU;
      m_stepV = source.m_stepV;
      m_tDeltaU = source.m_tDeltaU;
      m_tDeltaV = source.m_tDeltaV;
      m_tMaxU = source.m_tMaxU;
      m_tMaxV = source.m_tMaxV;
      m_tStart = source.m_tStart;
      m_validIntersection = source.m_validIntersection;
      return *this;
    }
    
    template <class ARRAY2D>
    std::pair<double, double>
    AmanatidesWoo2D<ARRAY2D>::
    findEntryAndExitPoints(const Vector2D& rayOriginPixel,
                           const Vector2D& rayDirectionPixel,
                           const ARRAY2D& data)
    {
      // First find intersection with each boundary line.

      // ... Find the intersection with the line U = 0, or else a really
      // small number if rayDirection is parallel to the rows of the
      // array.
      double tIntersectU0 = findIntersection(
        rayOriginPixel, rayDirectionPixel, Vector2D(1.0, 0.0), 0.0,
        std::numeric_limits<double>::min());

      // ... Find the intersection with the line U = data.columns(),
      // or else a really big number if rayDirection is parallel to the
      // rows of the array.
      double tIntersectU1 = findIntersection(
        rayOriginPixel, rayDirectionPixel, Vector2D(1.0, 0.0), data.columns(),
        std::numeric_limits<double>::max());

      // ... Find the intersection with the line V = 0, or else a really
      // small number if rayDirection is parallel to the columns of the
      // array.
      double tIntersectV0 = findIntersection(
        rayOriginPixel, rayDirectionPixel, Vector2D(0.0, 1.0), 0.0,
        std::numeric_limits<double>::min());
    
      // ... Find the intersection with the line V = data.rows(), or
      // else a really big number if rayDirection is parallel to the
      // columns of the array.
      double tIntersectV1 = findIntersection(
        rayOriginPixel, rayDirectionPixel, Vector2D(0.0, 1.0), data.rows(),
        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 tMaxU = std::max(tIntersectU0, tIntersectU1);
      double tMaxV = std::max(tIntersectV0, tIntersectV1);

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

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

    template <class ARRAY2D>
    double
    AmanatidesWoo2D<ARRAY2D>::
    findIntersection(const Vector2D& rayOrigin,
                     const Vector2D& rayDirection,
                     const Vector2D& 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 + t * rayDirection.
      // Substituting, we have:
      //   dot(rayOrigin + t * rayDirection, bVector) = cConstant.
      // Which we rewrite:
      //   dot(rayOrigin, bVector) + t * dot(rayDirection, bVector)
      //     = cConstant.
      // Solving for t, we have:
      //   t = (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_AMANATIDESWOO2D_H_