added python implementation of capsule projector

tools/projectSphere.py

@@ -135,7 +135,7 @@ def singleSphereToDetectorPixelRange(spherePositionMM, radiusMM, radiusMargin=0.
     return numpy.array([[minRC[0], maxRC[0]], [minRC[1], maxRC[1]]])
 
 
-def projectSphereMM(spheresPositionMM, radiiMM, ROIcentreMM=None, ROIradiusMM=None, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512,512], projector='C', transformationCentreMM=None, transformationMatrix=None, blur=None):
+def projectSphereMM(spheresPositionMM, radiiMM, ROIcentreMM=None, ROIradiusMM=None, sourceDetectorDistMM=100, pixelSizeMM=0.1, detectorResolution=[512,512], projector='C', transformationCentreMM=None, transformationMatrix=None, blur=None, displacementsMM=None):
     """
     This is the python wrapping function for the C++ projector, it gets projection geometry,
     list of particle positions and radii, projected in the X direction.
@@ -258,17 +258,42 @@ def projectSphereMM(spheresPositionMM, radiiMM, ROIcentreMM=None, ROIradiusMM=No
         if ROIcentreMM is None or ROIradiusMM is None:
             projectionXmm = numpy.zeros(detectorResolution, dtype=('<f4'))
 
-            # This C++ function fill in the passed projectionXmm array
-            projectSphereC3.project_func(numpy.array([sourceDetectorDistMM], dtype='<f4'),
-                                         radiiMM.astype('<f4'),
-                                         numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
-                                                         pixelSizeMM*detectorResolution[0]/2.,
-                                                         detectorResolution[0]).astype('<f4'),
-                                         numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
-                                                         pixelSizeMM*detectorResolution[1]/2.,
-                                                         detectorResolution[1]).astype('<f4'),
-                                         spheresPositionMM.astype('<f4'),
-                                         projectionXmm)
+            if displacementsMM is None:
+                # This C++ function fill in the passed projectionXmm array
+                projectSphereC3.project_func(numpy.array([sourceDetectorDistMM], dtype='<f4'),
+                                             radiiMM.astype('<f4'),
+                                             numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
+                                                             pixelSizeMM*detectorResolution[0]/2.,
+                                                             detectorResolution[0]).astype('<f4'),
+                                             numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
+                                                             pixelSizeMM*detectorResolution[1]/2.,
+                                                             detectorResolution[1]).astype('<f4'),
+                                             spheresPositionMM.astype('<f4'),
+                                             projectionXmm)
+            else:
+                # # This C++ function fill in the passed projectionXmm array
+                # projectSphereC3.project_capsule_func(numpy.array([sourceDetectorDistMM], dtype='<f4'),
+                #                              radiiMM.astype('<f4'),
+                #                              numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
+                #                                              pixelSizeMM*detectorResolution[0]/2.,
+                #                                              detectorResolution[0]).astype('<f4'),
+                #                              numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
+                #                                              pixelSizeMM*detectorResolution[1]/2.,
+                #                                              detectorResolution[1]).astype('<f4'),
+                #                              spheresPositionMM.astype('<f4'),
+                #                              displacementsMM.astype('<f4'),
+                #                              projectionXmm)
+                projectionXmm = project_capsule(numpy.array([sourceDetectorDistMM], dtype='<f4'),
+                                             radiiMM.astype('<f4'),
+                                             numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
+                                                             pixelSizeMM*detectorResolution[0]/2.,
+                                                             detectorResolution[0]).astype('<f4'),
+                                             numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
+                                                             pixelSizeMM*detectorResolution[1]/2.,
+                                                             detectorResolution[1]).astype('<f4'),
+                                             spheresPositionMM.astype('<f4'),
+                                             displacementsMM.astype('<f4'),
+                                             projectionXmm)
         elif ROIcentreMM is not None and ROIradiusMM is not None:
             # Make sure there's only one sphere:
             assert(len(spheresPositionMM.ravel())==3), "projectSphere.projectSphereMM(): in ROI mode I want only one sphere"
@@ -286,17 +311,31 @@ def projectSphereMM(spheresPositionMM, radiiMM, ROIcentreMM=None, ROIradiusMM=No
             # Define (smaller) projection array to fill in
             projectionXmm = numpy.zeros((limits[0,1]-limits[0,0], limits[1,1]-limits[1,0]), dtype=('<f4'))
 
-            # This C++ function fill in the passed projectionXmm array -- not limits after linspace
-            projectSphereC3.project_func(numpy.array([sourceDetectorDistMM], dtype='<f4'),
-                                         radiiMM.astype('<f4'),
-                                         numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
-                                                         pixelSizeMM*detectorResolution[0]/2.,
-                                                         detectorResolution[0]).astype('<f4')[limits[0,0]:limits[0,1]],
-                                         numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
-                                                         pixelSizeMM*detectorResolution[1]/2.,
-                                                         detectorResolution[1]).astype('<f4')[limits[1,0]:limits[1,1]],
-                                         spheresPositionMM.astype('<f4'),
-                                         projectionXmm)
+            if displacementsMM is None:
+                # This C++ function fill in the passed projectionXmm array -- not limits after linspace
+                projectSphereC3.project_func(numpy.array([sourceDetectorDistMM], dtype='<f4'),
+                                             radiiMM.astype('<f4'),
+                                             numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
+                                                             pixelSizeMM*detectorResolution[0]/2.,
+                                                             detectorResolution[0]).astype('<f4')[limits[0,0]:limits[0,1]],
+                                             numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
+                                                             pixelSizeMM*detectorResolution[1]/2.,
+                                                             detectorResolution[1]).astype('<f4')[limits[1,0]:limits[1,1]],
+                                             spheresPositionMM.astype('<f4'),
+                                             projectionXmm)
+            else:
+                 # This C++ function fill in the passed projectionXmm array
+                 projectSphereC3.project_capsule_func(numpy.array([sourceDetectorDistMM], dtype='<f4'),
+                                              radiiMM.astype('<f4'),
+                                              numpy.linspace(-pixelSizeMM*detectorResolution[0]/2.,
+                                                              pixelSizeMM*detectorResolution[0]/2.,
+                                                              detectorResolution[0]).astype('<f4'),
+                                              numpy.linspace(-pixelSizeMM*detectorResolution[1]/2.,
+                                                              pixelSizeMM*detectorResolution[1]/2.,
+                                                              detectorResolution[1]).astype('<f4'),
+                                              spheresPositionMM.astype('<f4'),
+                                              displacementsMM.astype('<f4'),
+                                              projectionXmm)
         else:
             print("projectSphere.projectSphereMM(): If you set ROIcentreMM you must also set ROIradiusMM")
 
@@ -430,3 +469,93 @@ def computeMotionKernel(posPrev, posPost, displacementThreshold=1):
         #kernelR = kernel/kernelSize
 
     return kernelR
+
+def project_capsule(SDD, radiiMM, widths, heights, positionsMM, displacementsMM, proj):
+    for sphere in range(len(positionsMM)):
+        axis_length = numpy.linalg.norm(displacementsMM[sphere])
+        V_sphere  = 4/3*numpy.pi*radiiMM[sphere]**3
+        V_capsule = V_sphere + numpy.pi*radiiMM[sphere]**2*axis_length
+        path_reduction_factor = V_sphere/V_capsule # a capsule has a lower attenuation because the same mass is smeared through a larger volume, so we need to lower the effective path length through the material
+
+        # define a local coordinate reference system with z along the direction of displacement, and x and y normal to that
+        z_unit_vector = displacementsMM[sphere] / axis_length
+        y_vector = numpy.cross(z_unit_vector, z_unit_vector*numpy.random.rand(1,3)) # dont really care which direction y is in as long as it is normal to z
+        x_vector = numpy.cross(y_vector, z_unit_vector) # and a third normal direction
+        A = numpy.vstack([x_vector/numpy.linalg.norm(x_vector), y_vector/numpy.linalg.norm(y_vector), z_unit_vector]) # rotation matrix to go from global to local coordinates
+        ray_source_rotated = numpy.dot(A,-positionsMM[sphere]) # apply the rotation matrix
+
+        for i,y in enumerate(widths):
+            for j,z in enumerate(heights):
+
+                ray_vector = numpy.array([SDD,y,z])
+                rotated_vector = numpy.dot(A,ray_vector)
+
+                path_length = _get_capsule_path_length(rotated_vector,
+                                                       ray_source_rotated,
+                                                       axis_length,
+                                                       radiiMM[sphere]
+                                                       )
+                proj[i,j] += path_length*path_reduction_factor
+
+    return proj
+
+def _get_capsule_path_length(ray_vector, ray_source, axis_length, radius):
+    path_length = 0
+
+    mag = numpy.linalg.norm(ray_vector)
+    ray_unit_vector = ray_vector/mag
+
+    # check for intersection with cylinder at ANY LOCAL z, i.e. does it cross the circle defined by the capsule. not a guarantee of collision with the capsule, just a collision with the infinite cylinder.
+    a = ray_unit_vector[0]**2 + ray_unit_vector[1]**2
+    b = 2*ray_unit_vector[0]*ray_source[0] + 2*ray_unit_vector[1]*ray_source[1]
+    c = ray_source[0]**2 + ray_source[1]**2 - radius**2
+
+    discriminant_squared = b**2 - 4*a*c
+    if discriminant_squared > 0: # there must be two solutions, check them both for collisions
+        l_0 = (-b + numpy.sqrt(discriminant_squared))/2/a
+        l_1 = (-b - numpy.sqrt(discriminant_squared))/2/a
+        l_2, x_0_good = _get_capsule_collision(l_0, ray_unit_vector, ray_source, radius, axis_length,  1)
+        l_3, x_1_good = _get_capsule_collision(l_1, ray_unit_vector, ray_source, radius, axis_length, -1)
+
+        if x_0_good and x_1_good: # if both have collisions, we have a ray that enters and leaves the capsule
+            path_length = l_2 - l_3
+
+    return path_length
+
+def _get_capsule_collision(l, u_hat, o, r, a, sign):
+    z = o[2] + u_hat[2]*l # z position at intercept with cylinder
+    if numpy.abs(z) < a/2: return [l, True] # inside the cylindrical part
+
+    if ( z > 0 ): z_offset = a/2 # pick one sphere to check
+    else: z_offset = -a/2
+    c = numpy.array([0,0,z_offset]) # centre of sphere
+    det = numpy.dot(u_hat, o - c)**2 - (numpy.linalg.norm(o - c)**2 - r**2)
+    if det > 0: # if the ray collides with the sphere
+        l = -numpy.dot(u_hat, o-c) + sign*numpy.sqrt(det) # pick the correct part of the sphere
+        return [l, True]
+    return [None, False] # nothing found
+
+if __name__ == '__main__':
+    # detector properties
+    ny = nz = 100
+    Ly = Lz = 10 # width and height of panel
+
+    SDD = 100 # distance from source to detector
+    w = numpy.linspace(-Ly/2,Ly/2,ny)
+    h = numpy.linspace(-Lz/2,Lz/2,nz)
+    nspheres = 10
+    radiiMM = (numpy.random.rand(nspheres) + 2)/3
+    pos = 8*numpy.random.rand(nspheres,3) - 4
+    pos[:,0] += 90
+    disp = numpy.random.rand(nspheres,3)
+    # disp = numpy.zeros([nspheres,3])
+
+    im = numpy.zeros([ny,nz])
+    project_capsule(SDD, radiiMM, w, h, pos, disp, im)
+
+    import matplotlib.pyplot as plt
+    plt.pcolormesh(w,h,im.T)
+    plt.xlabel('y')
+    plt.ylabel('z')
+    plt.colorbar()
+    plt.show()