Blame | Last modification | View Log | RSS feed
/* ----------------------------------------------------------------------- *
* This file is part of GEL, http://www.imm.dtu.dk/GEL
* Copyright (C) the authors and DTU Informatics
* For license and list of authors, see ../../doc/intro.pdf
* ----------------------------------------------------------------------- */
/**
* @file ThreeDDDA.h
* @brief a class for traversing cells in a regular grid pierced by a ray.
*/
#ifndef __GEOMETRY_THREEDDDA_H
#define __GEOMETRY_THREEDDDA_H
#ifdef _MSC_VER
typedef __int64 Integer64Bits;
#else
typedef int64_t Integer64Bits;
#endif
#include "../CGLA/Vec3f.h"
#include "../CGLA/Vec3i.h"
namespace Geometry
{
/** \brief A ThreeDDDA is a class for traversing a grid of cells.
We wish to
enumerate all cells pierced by the ray going from a point p0 to a
point p1. It is dangerous to use floating point arithmetic since rounding
errors may accumulate entailing that the 3ddda misses the cell containing
the end point of the ray.
Daniel Cohen-Or devised a technique for exact ray traversal using only
integer arithmetic. Using integer arithmetic, the computations are
exact, but the ray must start and terminate exactly at grid points.
I.e. the ray must start and terminate in cell corners.
To overcome this issue, I have implemented a scheme where the ray
starts and terminates on a grid that is finer than the cell grid.
This means that we have fast and exact computations, and the ray can
travel between two almost arbitrary points.
A central problem with this scheme was that - in order to simplify
things during the iterative traversal - the ray direction vector
is always mirrored into the first octant (+x +y +z). If the fine
grid used for the ray start and end points is a superset of the
grid points of the cell grid, this mirroring is almost impossible
to implement correctly. This is due to the fact that the first
cell on the right side of the y axis is called 0,y whereas the
first cell on the left is called -1,y. This lack of symmetry makes
things very difficult, but by ensuring that the endpoints of a ray
can never lie on an a corner, edge, or face of the cell grid, we
can make the problems go away.
Hence, in this implementation, the fine grid is obtained by dividing
the cell into NxNxN subcells, but the grid points of the fine grid lie
at the CENTRES of these subcells.
*/
class ThreeDDDA
{
/** Resolution of the grid. This number indicates how many subcells
(along each axis) a cell is divided into. */
const int PRECISION;
/// Sign indicates which octant contains the direction vector.
const CGLA::Vec3i sgn;
/// Fine grid position where the ray begins
const CGLA::Vec3i p0_int;
// Fine grid position where the ray ends.
const CGLA::Vec3i p1_int;
/// The direction of the ray
const CGLA::Vec3i dir;
/// The cell containing the initial point on the ray.
const CGLA::Vec3i first_cell;
/// The cell containing the last point on the ray.
const CGLA::Vec3i last_cell;
/** Number of steps. We can compute the number of steps that the
ray will take in advance. */
const int no_steps;
/** The direction vector, premultiplied by the step length.
This constant is used to update the equations governing the
traversal. */
const CGLA::Vec3i dir_pre_mul;
/** Discriminators. These values are delta[i]*dir[j] where delta[i]
is the distance from the origin of the ray to the current position
along the i axis multiplied by two). dir[j] is the direction vector
coordinate j. */
Integer64Bits disc[3][3];
/** The current cell containing the ray. This value along with the
discriminators mentioned above represent the state of a ThreeDDDA. */
CGLA::Vec3i cell;
/** Return the position on the fine grid of a vector v.
We simply multiply by the precision and round down. Note that
it is necessary to use floor since simply clamping the value
will always round toward zero.
*/
const CGLA::Vec3i grid_pos(const CGLA::Vec3f& v) const;
/// Compute the cell containing a given fine grid position.
CGLA::Vec3i get_cell(const CGLA::Vec3i& v) const;
/// Mirror a fine grid position
CGLA::Vec3i mirror(const CGLA::Vec3i& v) const;
public:
/** Create a 3DDDA based on the two end-points of the ray. An optional
last argument indicates the resolution of the grid. */
ThreeDDDA(const CGLA::Vec3f& p0, const CGLA::Vec3f& p1,
int _PRECISION = 0x100);
/// Return the total number of steps the ray will traverse.
int get_no_steps() const
{
return no_steps;
}
/// Yields the current cell containing the ray.
const CGLA::Vec3i& get_current_cell() const
{
return cell;
}
/// Get the initial cell containing the ray.
const CGLA::Vec3i& get_first_cell() const
{
return first_cell;
}
/// Get the last cell containing the ray.
const CGLA::Vec3i& get_last_cell() const
{
return last_cell;
}
/// Get the initial point containing the ray.
const CGLA::Vec3f get_first_point() const
{
return (CGLA::Vec3f(p0_int)+CGLA::Vec3f(0.5))/PRECISION;
}
/// Get the last cell containing the ray.
const CGLA::Vec3f get_last_point() const
{
return (CGLA::Vec3f(p1_int)+CGLA::Vec3f(0.5))/PRECISION;
}
/// Returns true if we have reached the last cell.
bool at_end() const
{
return cell == last_cell;
}
/// Take a step along the ray.
void operator++()
{
int step_dir;
if(disc[1][0] > disc[0][1])
if(disc[2][0] > disc[0][2])
step_dir = 0;
else
step_dir = 2;
else
if(disc[2][1] > disc[1][2])
step_dir = 1;
else
step_dir = 2;
const int k1 = (step_dir + 1) % 3;
const int k2 = (step_dir + 2) % 3;
cell[step_dir] += sgn[step_dir];
disc[step_dir][k1] += dir_pre_mul[k1];
disc[step_dir][k2] += dir_pre_mul[k2];
}
/// Get the initial point containing the ray.
const CGLA::Vec3f step();
};
}
#endif