Trilinear Interpolation on Voxels at specific angle

Question

I'm currently attempting to implement this algorithm for volume rendering in Python, and am conceptually confused about their method of generating the LH histogram (see section 3.1, page 4).

I have a 3D stack of DICOM images, and calculated its gradient magnitude and the 2 corresponding azimuth and elevation angles with it (which I found out about here), as well as finding the second derivative.

Now, the algorithm is asking me to iterate through a set of voxels, and "track a path by integrating the gradient field in both directions...using the second order Runge-Kutta method with an integration step of one voxel".

What I don't understand is how to use the 2 angles I calculated to integrate the gradient field in said direction. I understand that you can use trilinear interpolation to get intermediate voxel values, but I don't understand how to get the voxel coordinates I want using the angles I have.

In other words, I start at a given voxel position, and want to take a 1 voxel step in the direction of the 2 angles calculated for that voxel (one in the x-y direction, the other in the z-direction). How would I take this step at these 2 angles and retrieve the new (x, y, z) voxel coordinates?

Apologies in advance, as I have a very basic background in Calc II/III, so vector fields/visualization of 3D spaces is still a little rough for me.

Creating 3D stack of DICOM images:

def collect_data(data_path):
    print "collecting data"
    files = []  # create an empty list
    for dirName, subdirList, fileList in os.walk(data_path):
        for filename in fileList:
            if ".dcm" in filename:
                files.append(os.path.join(dirName,filename))
    # Get reference file
    ref = dicom.read_file(files[0])

    # Load dimensions based on the number of rows, columns, and slices (along the Z axis)
    pixel_dims = (int(ref.Rows), int(ref.Columns), len(files))

    # Load spacing values (in mm)
    pixel_spacings = (float(ref.PixelSpacing[0]), float(ref.PixelSpacing[1]), float(ref.SliceThickness))

    x = np.arange(0.0, (pixel_dims[0]+1)*pixel_spacings[0], pixel_spacings[0])
    y = np.arange(0.0, (pixel_dims[1]+1)*pixel_spacings[1], pixel_spacings[1])
    z = np.arange(0.0, (pixel_dims[2]+1)*pixel_spacings[2], pixel_spacings[2])

    # Row and column directional cosines
    orientation = ref.ImageOrientationPatient

    # This will become the intensity values
    dcm = np.zeros(pixel_dims, dtype=ref.pixel_array.dtype)

    origins = []

    # loop through all the DICOM files
    for filename in files:
        # read the file
        ds = dicom.read_file(filename)
        #get pixel spacing and origin information
        origins.append(ds.ImagePositionPatient) #[0,0,0] coordinates in real 3D space (in mm)

        # store the raw image data
        dcm[:, :, files.index(filename)] = ds.pixel_array
    return dcm, origins, pixel_spacings, orientation

Calculating gradient magnitude:

def calculate_gradient_magnitude(dcm):
    print "calculating gradient magnitude"
    gradient_magnitude = []
    gradient_direction = []

    gradx = np.zeros(dcm.shape)
    sobel(dcm,0,gradx)
    grady = np.zeros(dcm.shape)
    sobel(dcm,1,grady)
    gradz = np.zeros(dcm.shape)
    sobel(dcm,2,gradz)

    gradient = np.sqrt(gradx**2 + grady**2 + gradz**2)

    azimuthal = np.arctan2(grady, gradx)
    elevation = np.arctan(gradz,gradient)

    azimuthal = np.degrees(azimuthal)
    elevation = np.degrees(elevation)

    return gradient, azimuthal, elevation

Converting to patient coordinate system to get actual voxel position:

def get_patient_position(dcm, origins, pixel_spacing, orientation):
    """
        Image Space --> Anatomical (Patient) Space is an affine transformation
        using the Image Orientation (Patient), Image Position (Patient), and
        Pixel Spacing properties from the DICOM header
    """
    print "getting patient coordinates"

    world_coordinates = np.empty((dcm.shape[0], dcm.shape[1],dcm.shape[2], 3))
    affine_matrix = np.zeros((4,4), dtype=np.float32)

    rows = dcm.shape[0]
    cols = dcm.shape[1]
    num_slices = dcm.shape[2]

    image_orientation_x = np.array([ orientation[0], orientation[1], orientation[2]  ]).reshape(3,1)
    image_orientation_y = np.array([ orientation[3], orientation[4], orientation[5]  ]).reshape(3,1)
    pixel_spacing_x = pixel_spacing[0]

    # Construct affine matrix
    # Method from:
    # http://nipy.org/nibabel/dicom/dicom_orientation.html
    T_1 = origins[0]
    T_n = origins[num_slices-1]


    affine_matrix[0,0] = image_orientation_y[0] * pixel_spacing[0]
    affine_matrix[0,1] = image_orientation_x[0] * pixel_spacing[1]
    affine_matrix[0,3] = T_1[0]
    affine_matrix[1,0] = image_orientation_y[1] * pixel_spacing[0]
    affine_matrix[1,1] = image_orientation_x[1] * pixel_spacing[1]
    affine_matrix[1,3] = T_1[1]
    affine_matrix[2,0] = image_orientation_y[2] * pixel_spacing[0]
    affine_matrix[2,1] = image_orientation_x[2] * pixel_spacing[1]
    affine_matrix[2,3] = T_1[2]
    affine_matrix[3,3] = 1

    k1 = (T_1[0] - T_n[0])/ (1 - num_slices)
    k2 = (T_1[1] - T_n[1])/ (1 - num_slices)
    k3 = (T_1[2] - T_n[2])/ (1 - num_slices)

    affine_matrix[:3, 2] = np.array([k1,k2,k3])

    for z in range(num_slices):
        for r in range(rows):
            for c in range(cols):
                vector = np.array([r, c, 0, 1]).reshape((4,1))
                result = np.matmul(affine_matrix, vector)
                result = np.delete(result, 3, axis=0)
                result = np.transpose(result)
                world_coordinates[r,c,z] = result
        # print "Finished slice ", str(z)
    # np.save('./data/saved/world_coordinates_3d.npy', str(world_coordinates))
    return world_coordinates

Now I'm at the point where I want to write this function:

def create_lh_histogram(patient_positions, dcm, magnitude, azimuthal, elevation):
    print "constructing LH histogram"
    # Get 2nd derivative
    second_derivative = gaussian_filter(magnitude, sigma=1, order=1)

    # Determine if voxels lie on boundary or not (thresholding)

    # Still have to code out: let's say the thresholded voxels are in 
    # a numpy array called voxels

    #Iterate through all thresholded voxels and integrate gradient field in
    # both directions using 2nd-order Runge-Kutta
    vox_it = voxels.nditer(voxels, flags=['multi_index'])
    while not vox_it.finished:
        # ???