brainglobe / brainglobe.github.io

brainglobe.info website
https://brainglobe.info/
Creative Commons Attribution 4.0 International
10 stars 7 forks source link

[Feature] Document cell matching #187

Open alessandrofelder opened 3 months ago

alessandrofelder commented 3 months ago

Is your feature request related to a problem? Please describe. @matham added functionality to compare two sets of cells in https://github.com/brainglobe/brainglobe-utils/pull/74 but this is not documented.

Describe the solution you'd like A simple script showing how to use the top-level cell matching functionality.

Describe alternatives you've considered

-\

Additional context

-\

matham commented 3 months ago

For reference, this is the scripts I used:

Loading:

import pickle
from brainglobe_utils.IO.cells import get_cells_xml
from brainglobe_utils.cells.cells import Cell, match_cells
import numpy as np
import matplotlib.pyplot as plt

cells1 = get_cells_xml("cells_main.xml")
cells2 = get_cells_xml("cells_cuda.xml")

Optimization

missing_c1, good_matches, missing_c2 = match_cells(cells1, cells2, threshold=10)

with open(r"cells_cmp.pkl", "wb") as fh:
    pickle.dump((cells1, cells2, missing_c1, good_matches, missing_c2), fh)

Analysis:

with open(r"C:\Users\Matthew Einhorn\Downloads\cells_cmp.pkl", "rb") as fh:
    cells1, cells2, cells1_missing, matches, cells2_missing = pickle.load(fh)

print(len(cells1_missing), len(matches), len(cells2_missing))
matches = np.array(matches)
cells1_missing = np.array(cells1_missing)

cells1 = to_numpy_pos(cells1)
cells2 = to_numpy_pos(cells2)
a = cells1[matches[:, 0], :]
b = cells2[matches[:, 1], :]
dist = np.sqrt(np.sum(np.square(a - b), axis=1))
plt.hist(dist, bins=10)
plt.xlabel("Distance")
plt.ylabel("Count")
plt.title("Distribution of distances between matching cell pairs")
plt.show()

for i, name in enumerate(["X", "Y", "Z"]):
    plt.hist(cells1[cells1_missing, i], bins=1000)
    plt.xlabel(f"{name}-axis")
    plt.ylabel("Count")
    plt.title("Distribution of cells with no match")
    plt.show()