def trace
      # Set up instance varibles which depends on the initial conditions.
      # Delta positions determines whether the ray's travel is positive
      # or negative in the three directions.
      delta_x = @p1.x < @p2.x ? 1 : -1
      delta_y = @p1.y < @p2.y ? 1 : -1
      delta_z = @p1.z < @p2.z ? 1 : -1
      # Delta indices determines whether the ray's voxel space indices increments
      # in a positive or negative fashion as it travels through the voxel space.
      # Note that for rays travelling perpendicular on a particular axis, the
      # delta value of that axis will be undefined (zero).
      @delta_i = @p1.x == @p2.x ? 0 : delta_x
      @delta_j = @p1.y == @p2.y ? 0 : delta_y
      @delta_k = @p1.z == @p2.z ? 0 : delta_z
      # These variables describe how much the alpha fraction changes (in a
      # given axis) when we follow the ray through one voxel (across two planes).
      # This value is high if the ray is perpendicular on the axis (angle high, or close to 90 degrees),
      # and low if the angle is parallell with the axis (angle low, or close to 0 degrees).
      @delta_ax = @vs.delta_x / (@p2.x - @p1.x).abs
      @delta_ay = @vs.delta_y / (@p2.y - @p1.y).abs
      @delta_az = @vs.delta_z / (@p2.z - @p1.z).abs
      # Determines the ray length (from p1 to p2).
      @length = Math.sqrt((@p2.x-@p1.x)**2 + (@p2.y-@p1.y)**2 + (@p2.z-@p1.x)**2)
      # Perform the ray tracing:
      # Number of voxels: nx, ny, nz
      # Number of planes: nx+1, ny+1, nz+1
      # Voxel indices: 0..(nx-1), etc.
      # Plane indices: 0..nx, etc.
      intersection_with_voxel_space
      # Perform the ray trace only if the ray's intersection with the
      # voxel space is between the ray's start and end points:
      if @alpha_max > 0 && @alpha_min < 1
        min_and_max_indices
        number_of_planes
        axf, ayf, azf = first_intersection_point_in_voxel_space
        # If we didn't get any intersection points, there's no need to proceed with the ray trace:
        if [axf, ayf, azf].any?
          # Initialize the starting alpha values:
          initialize_alphas(axf, ayf, azf)
          @i, @j, @k = indices_first_intersection(axf, ayf, azf)
          # Initiate the ray tracing if we got valid starting coordinates:
          if @i && @j && @k
            # Calculate the first move:
            update
            # Next interactions:
            # To avoid float impresision issues we choose to round here before
            # doing the comparison. How many decimals should we choose to round to?
            # Perhaps there is a more prinipal way than the chosen solution?
            alpha_max_rounded = @alpha_max.round(8)
            while @ac.round(8) < alpha_max_rounded do
              update
            end
          end
        end
      end
    end
def trace
      # Set up instance varibles which depends on the initial conditions.
      # Delta positions determines whether the ray's travel is positive
      # or negative in the three directions.
      delta_x = @p1.x < @p2.x ? 1 : -1
      delta_y = @p1.y < @p2.y ? 1 : -1
      delta_z = @p1.z < @p2.z ? 1 : -1
      # Delta indices determines whether the ray's voxel space indices increments
      # in a positive or negative fashion as it travels through the voxel space.
      # Note that for rays travelling perpendicular on a particular axis, the
      # delta value of that axis will be undefined (zero).
      @delta_i = @p1.x == @p2.x ? 0 : delta_x
      @delta_j = @p1.y == @p2.y ? 0 : delta_y
      @delta_k = @p1.z == @p2.z ? 0 : delta_z
      # These variables describe how much the alpha fraction changes (in a
      # given axis) when we follow the ray through one voxel (across two planes).
      # This value is high if the ray is perpendicular on the axis (angle high, or close to 90 degrees),
      # and low if the angle is parallell with the axis (angle low, or close to 0 degrees).
      @delta_ax = @vs.delta_x / (@p2.x - @p1.x).abs
      @delta_ay = @vs.delta_y / (@p2.y - @p1.y).abs
      @delta_az = @vs.delta_z / (@p2.z - @p1.z).abs
      # Determines the ray length (from p1 to p2).
      @length = Math.sqrt((@p2.x-@p1.x)**2 + (@p2.y-@p1.y)**2 + (@p2.z-@p1.x)**2)
      # Perform the ray tracing:
      # Number of voxels: nx, ny, nz
      # Number of planes: nx+1, ny+1, nz+1
      # Voxel indices: 0..(nx-1), etc.
      # Plane indices: 0..nx, etc.
      intersection_with_voxel_space
      # Perform the ray trace only if the ray's intersection with the
      # voxel space is between the ray's start and end points:
      if @alpha_max > 0 && @alpha_min < 1
        min_and_max_indices
        number_of_planes
        axf, ayf, azf = first_intersection_point_in_voxel_space
        # If we didn't get any intersection points, there's no need to proceed with the ray trace:
        if [axf, ayf, azf].any?
          # Initialize the starting alpha values:
          initialize_alphas(axf, ayf, azf)
          @i, @j, @k = indices_first_intersection(axf, ayf, azf)
          # Initiate the ray tracing if we got valid starting coordinates:
          if @i && @j && @k
            # Calculate the first move:
            update
            # Next interactions:
            # To avoid float impresision issues we choose to round here before
            # doing the comparison. How many decimals should we choose to round to?
            # Perhaps there is a more prinipal way than the chosen solution?
            alpha_max_rounded = @alpha_max.round(8)
            while @ac.round(8) < alpha_max_rounded do
              update
            end
          end
        end
      end
    end
Performs ray tracing, where the ray's possible intersection of the
 associated voxel space is investigated for the ray's movement from
 its source coordinate to its target coordinate.

 The resulting density experienced by the ray through the voxel space
 is stored in the 'd' attribute. The indices of the voxel space
 intersected by the ray is stored in the 'indices' attribute (if the
 'path' option has been set).