<jupyter_start><jupyter_text>Hollender 3D image analysis Author: <NAME> Last updated: 2021-04-16 Goals:* Read in .ply files* Process images to get baseline skeletons that all subgroups will work with* Measure cool phenotypes Link to cookbook: http://www.open3d.org/docs/release/tutorial/geometry/pointcloud.html  Import modules<jupyter_code>import numpy as np
import open3d as o3d
from matplotlib import pyplot as plt
import skimage.morphology
import sknw
import scipy.ndimage<jupyter_output><empty_output><jupyter_text>Write functions<jupyter_code>#Write function to take pcd object and plane formula then subset points
def remove_plane(pcdObject, modelArray, upr, lwr):
    #Write list for output
    goodIndices = []
    #Convert pcd object to numpy array, save points; normals; and colors
    pcdPoints = np.asarray(pcdObject.points)
    pcdNormals = np.asarray(pcdObject.normals)
    pcdColors = np.asarray(pcdObject.colors)
    #Extract coefficients for plane equation
    a = plane_model[0]
    b = plane_model[1]
    c = plane_model[2]
    d = plane_model[3]
    #evaluate plane equation for each point
    for i in range(len(pcdPoints)):
        x = pcdPoints[i][0]
        y = pcdPoints[i][1]
        z = pcdPoints[i][2]
        if(a*x + b*y + c*z + d < upr and a*x + b*y + c*z + d > lwr):
            goodIndices.append(i)
    #subset array, convert array back to cloud
    pcdNew = o3d.geometry.PointCloud()
    pcdNew.points = o3d.utility.Vector3dVector(pcdPoints[goodIndices])
    pcdNew.normals = o3d.utility.Vector3dVector(pcdNormals[goodIndices])
    pcdNew.colors = o3d.utility.Vector3dVector(pcdColors[goodIndices])
    return pcdNew


######################################################################
def remove_xzplane(pcdObject, lwr, upr):
    #Write list for output
    goodIndices = []
    #Convert pcd object to numpy array, save points; normals; and colors
    pcdPoints = np.asarray(pcdObject.points)
    pcdNormals = np.asarray(pcdObject.normals)
    pcdColors = np.asarray(pcdObject.colors)
    #evaluate plane equation for each point
    for i in range(len(pcdPoints)):
        x = abs(pcdPoints[i][0])
        if(x < upr and x > lwr):
            goodIndices.append(i)
    #subset array, convert array back to cloud
    pcdNew = o3d.geometry.PointCloud()
    pcdNew.points = o3d.utility.Vector3dVector(pcdPoints[goodIndices])
    pcdNew.normals = o3d.utility.Vector3dVector(pcdNormals[goodIndices])
    pcdNew.colors = o3d.utility.Vector3dVector(pcdColors[goodIndices])
    return pcdNew

######################################################################
#Write function to take pcd object and plane formula then subset points
def remove_xyplane(pcdObject, lwr, upr):
    #Write list for output
    goodIndices = []
    #Convert pcd object to numpy array, save points; normals; and colors
    pcdPoints = np.asarray(pcdObject.points)
    pcdNormals = np.asarray(pcdObject.normals)
    pcdColors = np.asarray(pcdObject.colors)
    #evaluate plane equation for each point
    for i in range(len(pcdPoints)):
        y = pcdPoints[i][1]
        if(y < upr and y > lwr):
            goodIndices.append(i)
    #subset array, convert array back to cloud
    pcdNew = o3d.geometry.PointCloud()
    pcdNew.points = o3d.utility.Vector3dVector(pcdPoints[goodIndices])
    pcdNew.normals = o3d.utility.Vector3dVector(pcdNormals[goodIndices])
    pcdNew.colors = o3d.utility.Vector3dVector(pcdColors[goodIndices])
    return pcdNew

###Code to show which points are removed as outliers
#got code from here: http://www.open3d.org/docs/release/tutorial/geometry/pointcloud_outlier_removal.html#Statistical-outlier-removal
def display_inlier_outlier(cloud, ind):
    inlier_cloud = cloud.select_by_index(ind)
    outlier_cloud = cloud.select_by_index(ind, invert=True)

    print("Showing outliers (red) and inliers (gray): ")
    outlier_cloud.paint_uniform_color([1, 0, 0])
    inlier_cloud.paint_uniform_color([0.8, 0.8, 0.8])
    o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud])

######################################################################
#Input: open3d point cloud object, and a whole number of intervals to partition each dimmension into
#Output: a matrix of true/false values marking which voxels include at least one point
def pcd2bimg(pcdObject, resolution):
    #Extract coordinates from point cloud
    Coord = np.asarray(pcdObject.points)
    xCoord = [Coord[i][0] for i in range(len(Coord))]
    yCoord = [Coord[i][1] for i in range(len(Coord))]
    zCoord = [Coord[i][2] for i in range(len(Coord))]

    #Get vertices of overall bounding box
    xmin = min(xCoord)
    xmax = max(xCoord)
    xrange = xmax - xmin

    ymin = min(yCoord)
    ymax = max(yCoord)
    yrange = ymax - ymin

    zmin = min(zCoord)
    zmax = max(zCoord)
    zrange = zmax - zmin

    #calculate side lengths of voxels (i.e. the l's). Higher number of voxels (i.e. the n's) means higher resolution
    xl = xrange/resolution
    yl = yrange/resolution
    zl = zrange/resolution

    #Create empty 3d numpty array to store result
    bimg = np.zeros((resolution, resolution, resolution))

    #Loop over individual voxels, label true if they contain points
    for k in range(resolution):
        zvox1 = zmin + zl*k #determine bounds of voxel
        zvox2 = zmin + zl*(k+1)
        zvoxmin = min(zvox1, zvox2)
        zvoxmax = max(zvox1, zvox2)
        for j in range(resolution):
            yvox1 = ymin + yl*j
            yvox2 = ymin + yl*(j+1)
            yvoxmin = min(yvox1, yvox2)
            yvoxmax = max(yvox1, yvox2)
            for i in range(resolution):
                xvox1 = xmin + xl*i
                xvox2 = xmin + xl*(i+1)
                xvoxmin = min(xvox1, xvox2)
                xvoxmax = max(xvox1, xvox2)
                #Loop over the point cloud coordinates, if a coordinate falls within bounds of voxel, label entire voxel with a 1
                for row in Coord:
                    if xvoxmin < row[0] < xvoxmax and yvoxmin < row[1] < yvoxmax and zvoxmin < row[2] < zvoxmax:
                        bimg[i,j,k] = 1
                        break

    #Convert matrix to boolean type
    bimgBool = np.array(bimg, dtype=bool)
    
    return bimgBool<jupyter_output><empty_output><jupyter_text>Read in point cloud data<jupyter_code># Read .ply file
input_file = "./Hollender_arabidopsis/arabidopsis lab test 4 (col).ply"
pcd = o3d.io.read_point_cloud(input_file) # Read the point cloud

# Visualize the point cloud within open3d
o3d.visualization.draw_geometries([pcd])<jupyter_output><empty_output><jupyter_text>Isolate plant in point cloud  Downsample<jupyter_code>#Downsampling
downpcd = pcd.voxel_down_sample(voxel_size=0.015)
o3d.visualization.draw_geometries([downpcd])<jupyter_output><empty_output><jupyter_text>Remove plane<jupyter_code>#Identify plane in downsampled figure
plane_model, inliers = downpcd.segment_plane(distance_threshold=0.01,
                                         ransac_n=3,
                                         num_iterations=1000)
[a, b, c, d] = plane_model
print(f"Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0")

inlier_cloud = downpcd.select_by_index(inliers)
inlier_cloud.paint_uniform_color([1.0, 0, 0])
outlier_cloud = downpcd.select_by_index(inliers, invert=True)
o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud])

#Remove plane
noPlanePcd = remove_plane(downpcd, plane_model, 0.1, 0.005)
o3d.visualization.draw_geometries([noPlanePcd])<jupyter_output>Plane equation: 0.01x + -0.05y + 1.00z + 0.55 = 0<jupyter_text>Crop along the other axes to isolate the plant<jupyter_code>noXZplanePcd = remove_xzplane(noPlanePcd, 0, 0.15)
o3d.visualization.draw_geometries([noXZplanePcd])
noXYplanePcd = remove_xyplane(noXZplanePcd, -0.02, 0.21)
o3d.visualization.draw_geometries([noXYplanePcd])<jupyter_output><empty_output><jupyter_text>Finally, Use DBSCAN clustering to find largest point cloud (i.e. the plant)<jupyter_code>with o3d.utility.VerbosityContextManager(o3d.utility.VerbosityLevel.Debug) as cm:
    labels = np.array(noXYplanePcd.cluster_dbscan(eps=0.1, min_points=30, print_progress=True))

max_label = labels.max()
print(f"point cloud has {max_label + 1} clusters")
colors = plt.get_cmap("tab20")(labels / (max_label if max_label > 0 else 1))
colors[labels < 0] = 0
noXYplanePcd.colors = o3d.utility.Vector3dVector(colors[:, :3])
o3d.visualization.draw_geometries([noXYplanePcd])
finalPoints = np.asarray(noXYplanePcd.points)[np.where(labels==0)]
finalPcd = o3d.geometry.PointCloud()
finalPcd.points = o3d.utility.Vector3dVector(finalPoints)
o3d.visualization.draw_geometries([finalPcd])<jupyter_output><empty_output><jupyter_text>Filter outlier points<jupyter_code>pcdNoOut, ind = finalPcd.remove_statistical_outlier(nb_neighbors=20,std_ratio=1.25)
display_inlier_outlier(finalPcd, ind)
o3d.visualization.draw_geometries([pcdNoOut])<jupyter_output><empty_output><jupyter_text>Construct graph of plant point cloud  Parition points into voxels to be skeletonized later<jupyter_code>bimg = pcd2bimg(pcdNoOut, 200)<jupyter_output><empty_output><jupyter_text>Skeletonize image, visualize skeleton<jupyter_code>bimgSkel = skimage.morphology.skeletonize_3d(bimg).astype(np.uint16)
#bimgErod = skimage.morphology.binary_erosion(bimg)
bimgErod = scipy.ndimage.binary_erosion(bimg, iterations = 2)

graph = sknw.build_sknw(bimgSkel, multi=False)

#Extract location of individual nodes, convert back to point cloud object, draw edges of graph
lines = [[s,e] for (s,e) in graph.edges()]
points = [[graph.nodes[i]['pts'][0][0],graph.nodes[i]['pts'][0][1],graph.nodes[i]['pts'][0][2]] for i in range(len(graph.nodes()))]

colors = [[1, 0, 0] for i in range(len(lines))] #All lines will be colored red
line_set = o3d.geometry.LineSet(
    points=o3d.utility.Vector3dVector(points),
    lines=o3d.utility.Vector2iVector(lines),
)
line_set.colors = o3d.utility.Vector3dVector(colors)
o3d.visualization.draw_geometries([line_set])

#draw just point cloud
bimgPcd = o3d.geometry.PointCloud()
bimgPcd.points = o3d.utility.Vector3dVector(points)
o3d.visualization.draw_geometries([bimgPcd])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12550447878728221, 'fraction_numerical': 0.01978541194999508, 'mean_word_length': 4.204918032786885, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '150047', 'license': 'MIT', 'n_tokens_mistral': 3750, 'n_tokens_neox': 3593, 'n_words': 790, 'path': 's00_HollenderImageAnalysis_v8.ipynb', 'repo_name': 'thompsonmaizelab/CSS844-module3-Hollender', 'size': '10166'}