Tetrahexacontree Voxel Volume Raymarching

variables

    tree bounding box (bounds)
    tree root node index
    nodes -- array of tree nodes
    volumes -- array of tree voxel volumes:
    far -- far clipping plane distance
    near -- near clipping plane distance
    origin -- ray origin
    direction -- ray direction
    t -- t distance along ray (starts at near)
    ray location --  origin + (direction * t)
    width -- width of child node bounds
    bounds0 -- bounds of current node
    bounds1 -- bounds of voxel volume
    stack -- stack of bounding boxes (bounds) and tree node indices
    frame -- current bounding box (bounds) and tree node index
    current node -- tree.nodes[index] or tree.nodes[frame.index]
    current volume -- tree.volumes[index]

algorithm

  * if tree bounds does not contain ray
    * updated ray with intersection of ray with tree bounds
      * if fails
        * return far
  * initialize width with tree bounds width divided by four
  * initialize stack
  * push tree bounds and tree node root index on to stack
  * set frame equal to value from top of stack
  * while true
    * set bounds0 to bounds of current node based on width and ray location -- use ray location aligned by width to calculate min and max of bounds
    * calculate octlet of child node based on frame bounds, ray location and width
    * if current node's child (nodes[frame.index].chilren[octlet]) is that of a valid child node (not null sentinal)
      * adjust frame bounds based on ray location
      * set frame index to index of child node
      * if width equal to four times that of a voxel volume width
        * while true
          * calculate bounds (bounds1) of voxel volume based on width divided by four
          * calculate octlet of voxel volume based on frame bounds, ray location and width divided by four
          * if current node's child (nodes[frame.index].chilren[octlet]) is that of a valid child node (not null sentinal)
            * get volume based on frame index and octlet (volumes_[nodes_[frame.index].child(octlet)]
            * initialize digital differential (dda) analizer from ray's origin, direction and t
            * while true -- step through volume's voxel using dda
              * convert dda's location to volume coordinates
              * if non-air voxel hit in volume at converted coordinates
                * return t and hit location
              * increment dda
              * if dda t greater than far, return far
              * if location of dda beyond bounds of volume (bounds1)
                * break from loop
            * set ray's t to that of dda's t
            * if ray location outside frame bounds
              * if ray's t greater than far or ray location outside of tree bounds
                * return far
              * break from loop
          * else
            * calculate intersection of next potential volume using volume's bounds (bounds1) and update ray location
            * if intersection fails, ray t greater than far or ray location outside of tree bounds
              * return far
            * if ray location outside of frame bounds, break from loop
        * if stack size == 1
          * set frame to value at top of stack (frame = stack.top())
        * else
          * pop from stack into frame until frame has ray location multiplying width by four for each pop
      * else
        * push frame onto stack
        * divide width by four
    * else
      * calculate intersection of next potential node using current node's bounds (bounds0) and update ray location
      * if intersection fails, ray's t greater than far or ray location outside of tree bounds
        * return far
      * if stack size == 1
        * set frame to value at top of stack
      * else
        * if ray location outside of frame bounds
          * pop from stack into frame until frame has ray location multiplying width by four for each pop