Measure diameter of curved 3d tube

Hello everyone,

I am working with a 3d curved tubular structure. I have extracted the centerline at various points and now I am trying to calculate the diameter of the tube at those points. At the curved regions, I am trying to use a 3D Euler rotation to align the tube with the z-axis so I can compute the diameter.

I am struggling to rotate the tube around the point so that I can extract the cross-section and the transformed point in the rotated space.

My current approach is to calculate the displacement between the two closest points in the centerline. From there, I compute the Euler rotation angles and apply the rotation to the data and a transformation to the center line point.

It appears that my rotation of the data is correct, however I am unable to verify this because the point transformation does not appear correct.

Here is the code I am using (please note the point is in numpy format so I have to change the axis)

def rotation3d(image, theta_x, theta_y, theta_z, point):
    This function rotates an image across each of the x, y, z axes by theta_x, theta_y, and theta_z degrees
    :param image: An sitk MRI image
    :param theta_x: The amount of degrees the user wants the image rotated around the x axis
    :param theta_y: The amount of degrees the user wants the image rotated around the y axis
    :param theta_z: The amount of degrees the user wants the image rotated around the z axis
    :return: The rotated image
    point = (int(point[2]), int(point[1]), int(point[0]))    
    euler_transform = sitk.Euler3DTransform(image.TransformContinuousIndexToPhysicalPoint([(sz-1)/2.0 for sz in image.GetSize()]), theta_x, theta_y, theta_z, (0, 0, 0))

    euler_transform.SetRotation(theta_x, theta_y, theta_z)

    max_indexes = [sz-1 for sz in image.GetSize()]
    extreme_indexes = list(itertools.product(*(list(zip([0]*image.GetDimension(),max_indexes)))))
    extreme_points_transformed = [euler_transform.TransformPoint(image.TransformContinuousIndexToPhysicalPoint(p)) for p in extreme_indexes]
    output_min_coordinates = np.min(extreme_points_transformed, axis=0)
    output_max_coordinates = np.max(extreme_points_transformed, axis=0)

    output_spacing = image.GetSpacing()

    output_origin = output_min_coordinates
    output_size = [int(((omx-omn)/ospc)+0.5)  for ospc, omn, omx in zip(output_spacing, output_min_coordinates, output_max_coordinates)]
    output_direction = [1,0,0,0,1,0,0,0,1]
    output_pixeltype = image.GetPixelIDValue()
    output = sitk.Resample(image, 
    #resampled_image = resample(image, euler_transform)

    point = euler_transform.TransformPoint(point)
    point = (int(point[2]), int(point[1]), int(point[0]))

    return output, point

This code is adapted from SimpleITK Euler3D transform, problem with output size/resampling - #2 by dzenanz

Thank you and any help is appreciated!

@Stephen_Aylward might have some suggestions. If he converts his tube to ITK tube, will he get GetRadiusInObjectSpace for free?
@cams2b Take a look at Demo-ManipulateTubes.