Integrated Analysis of Cell Shape and Movement in Moving Frame

# Add necessary packages

pkg"add  GaussianProcesses Conda"

# Add python package igl from conda-forge as one of MovingFrame dependency

using Conda

Conda.add("igl";channel="conda-forge")

# Add MovingFrame package from github repository

pkg"add https://github.com/yusri-dh/MovingFrame.jl"

pkg"precompile"

tar -xzvf 
cell_shape_final.tar.gz

using MovingFrame;

using GeometryBasics

using FileIO

using LinearAlgebra;

using GaussianProcesses;

import Plots;

const mva = MovingFrame

const plt = Plots

# Create MovingCells Type

struct MovingCells

    coordinates::Matrix{Float64}

    time::Vector{Float64}

    smooth_n_points::Int64

    smooth_coordinates::Matrix{Float64}

    smooth_time::Vector{Float64}

    T::Matrix{Float64}

    N::Matrix{Float64}

    B::Matrix{Float64}

    curvature::Vector{Float64}

    torsion::Vector{Float64}

    speed::Vector{Float64}

    original_shapes::Vector{GeometryBasics.Mesh}

    reoriented_shapes::Vector{GeometryBasics.Mesh}

    Spharm_Cx::Matrix{Float64}

    Spharm_Cy::Matrix{Float64}

    Spharm_Cz::Matrix{Float64}

end

const base_dir = pwd();

# You need to extract the cell_shape_final.zip first

const input_data_dir = joinpath(base_dir, "cell_shape_final/")

const input_files = readdir(input_data_dir);

const n_points = length(input_files)

# Calculating the cell trajectory

coordinate = Matrix{Float64}(undef, n_points, 3)

time = collect(1.0:n_points)

shapes = Vector{GeometryBasics.Mesh}(undef, n_points)

for i = 1:length(input_files)

    shape = load(joinpath(input_data_dir, input_files[i]))

    shapes[i] = shape

    verts = decompose(Point3f0, shape)

    # The trajectory is the mass centre of the cell

    coordinate[i, :] .= reduce(+, verts) / length(verts)

end

# Normalizing volume and translation

shapes = (volume_normalizing ∘ mesh_centering).(shapes)

new_coordinate = coordinate .- coordinate[1, :]';

exp_lengthscale = 3.4 # lengthscale = e^3.4

exp_variance= 0.0 #variance = e^0.0

# Gaussian Processs using Zero mean function and 

# Squared exponential kernel (note that hyperparameters are on the log scale)

smoother =  GaussianProcessSmoother(MeanZero(), SE(exp_lengthscale, exp_variance))

smooth_coordinates = smoothing(new_coordinate, smoother, step=1.)

smooth_time = range(1.0, stop=time[end], length=size(smooth_coordinates, 1)) |> collect

T, N, B, curvature, torsion, speed = parallel_transport(smooth_coordinates);

# Reoriented the cell shape to the Moving Frame basis

reoriented_shape = [

    mesh_reorientation(shapes[i], T[i, :], N[i, :], B[i, :])

    for i in eachindex(shapes)

];

# Spherical parameterization

spheres = MovingFrame.spherization.(reoriented_shape, n_step=100);

# Calculating spherical harmonics coefficients

spharm_coefs = [

    MovingFrame.spharm_descriptor(

        reoriented_shape[i],

        GeometryBasics.coordinates(spheres[i]),

    ) for i in eachindex(reoriented_shape)

]

Cx = MovingFrame.array_of_array_to_matrix([coef[:Cx] for coef in spharm_coefs])

Cy = MovingFrame.array_of_array_to_matrix([coef[:Cy] for coef in spharm_coefs])

Cz = MovingFrame.array_of_array_to_matrix([coef[:Cz] for coef in spharm_coefs])

xy_ecc = Cx[:, 4] ./ Cy[:, 2]

xz_ecc = Cx[:, 4] ./ Cz[:, 3]

yz_ecc = Cy[:, 2] ./ Cz[:, 3]

plt.plot(smooth_coordinates[:,1],smooth_coordinates[:,2],smooth_coordinates[:,3],label ="SMOOTH_TRAJECTORY")

plt.scatter!(new_coordinate[:,1],new_coordinate[:,2],new_coordinate[:,3],label ="OBSERVATION",markersize=1.0)

plt.plot(smooth_time,[xy_ecc xz_ecc],label=["XY Eccentricity" "XZ Eccentricity"])

begin

	p1 = plt.plot(smooth_time,speed,ylabel = "Speed")

	p2 = plt.plot(smooth_time,curvature,ylabel = "Curvature",)

	p3 = plt.plot(smooth_time,torsion,ylabel="Torsion",)

	plt.plot(p1,p2,p3,legend=false)

end