Skip to content

Commit

Permalink
Minor tweaks to processing
Browse files Browse the repository at this point in the history
  • Loading branch information
jojoelfe committed Oct 10, 2023
1 parent d5b231b commit 335ca73
Show file tree
Hide file tree
Showing 5 changed files with 251 additions and 38 deletions.
105 changes: 105 additions & 0 deletions src/decolace/acquisition/cli_acquisition_test.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,105 @@
import glob
import os
from enum import Enum
from pathlib import Path

import typer
from rich.panel import Panel
from rich.progress import Progress
import numpy as np
import shapely
import signal

class GracefulKiller:
kill_now = False
def __init__(self):
signal.signal(signal.SIGINT, self.exit_gracefully)
signal.signal(signal.SIGTERM, self.exit_gracefully)

def exit_gracefully(self, *args):
print("Recieved SIGTERM")
self.kill_now = True



app = typer.Typer()


def get_image_fn(dummy):
import serialem as sem
sem.ConnectToSEM(48888,"192.168.122.149")
image = np.asarray(sem.bufferImage("P")).copy()
return image


@app.command()
def get_image():
import napari
from napari.layers import Image
from magicgui import magicgui
from multiprocessing import Pool

@magicgui(call_button="Get image")
def my_widget() -> Image:
pool = Pool(processes = 1)
image = pool.map(get_image_fn,[0])[0]
return Image(image)

viewer = napari.Viewer()
viewer.window.add_dock_widget(my_widget)
napari.run()
print("Done")

@app.command()
def setup_areas(
session_name: str = typer.Option(None, help="Name of the session"),
directory: str = typer.Option(None, help="Directory to save session in"),
):
session_o = load_session(session_name, directory)

num_items = serialem.ReportNumTableItems()
maps = []
map_coordinates = []
map_navids = []
for i in range(1,int(num_items)+1):
nav_item_info = serialem.ReportOtherItem(i)
nav_id = int(nav_item_info[0])
nav_note = serialem.GetVariable("navNote")

if nav_note == "decolace_acquisition_map":
serialem.LoadOtherMap(i,"A")
image = np.asarray(serialem.bufferImage("A")).copy()
maps.append(image)
map_navids.append(nav_id)
map_coordinates.append(nav_item_info[1:3])
import napari
from napari.layers import Shapes
from magicgui import magicgui
viewer = napari.view_image(np.array(maps))

@magicgui(shapes={'label': 'Setup areas'})
def my_widget(shapes: Shapes):
points = []
areas = shapes.data
for area in areas:
map_id = area[0,0]
if np.sum(area[:,0] - map_id) != 0:
raise("Error: Map ID is not the same for all points in the polygon")
name = f"area{len(session_o.active_grid.state['acquisition_areas'])+2}"
polygon = shapely.geometry.Polygon(area[:,1:3])
aa = AcquisitionAreaSingle(name,Path(session_o.active_grid.directory,name).as_posix(),beam_radius=session_o.state["beam_radius"],tilt=session_o.active_grid.state["tilt"])
aa.initialize_from_napari(map_navids[int(map_id)], [polygon.centroid.y, polygon.centroid.x], area[:,1:3])
aa.calculate_acquisition_positions_from_napari()
aa.write_to_disk()
session_o.active_grid.state["acquisition_areas"].append([aa.name,aa.directory])
session_o.active_grid.write_to_disk()
session_o.active_grid.save_navigator()



viewer.window.add_dock_widget(my_widget)
napari.run()
print("Done")
typer.Exit()


61 changes: 25 additions & 36 deletions src/decolace/analysis/cli_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,11 +51,6 @@ def cluster_by_distance_and_orientation(
import starfile
import numpy as np
from scipy.spatial import KDTree
from scipy.sparse.csgraph import connected_components
from healpy import pix2ang
from collections import defaultdict
from sympy import Plane, Point3D
from decolace.analysis.orientations import calculate_angle_between_vectors
import eulerangles

input_dir = Path(ctx.obj.project.project_path) / "Matches"
Expand All @@ -77,41 +72,35 @@ def cluster_by_distance_and_orientation(
pairs = tree.query_pairs(r=distance_threshold, output_type='ndarray')


data['NUMCLUSTERS'] = 0
data['CLUSTER_SCORE'] = 0
print(f"Distance {distance_threshold} has {len(pairs)} pairs")
# Get all pairs that ar aligned correctly
vectors1 = np.dot(rotation_matrices[pairs[:,0]], np.array([1,0,0]))
vectors2 = np.dot(rotation_matrices[pairs[:,1]], np.array([1,0,0]))
dotproducts = np.sum(vectors1 * vectors2, axis=1)
vectors = np.dot(rotation_matrices, np.array([0,0,1]))
#vectors1 = np.dot(rotation_matrices[pairs[:,0]], np.array([1,0,0]))
#vectors2 = np.dot(rotation_matrices[pairs[:,1]], np.array([1,0,0]))
dotproducts = np.sum(vectors[pairs[:,0]] * vectors[pairs[:,1]], axis=1)
angles = np.degrees(np.arccos(dotproducts)) # No normalization needed as long as rotated vector is a unit vector
potential_clusers = pairs[angles < orientation_threshold]
potential_clusters_midpoints = (coordinates[potential_clusers[:,0]] + coordinates[potential_clusers[:,1]]) / 2
potential_additional_candidates = tree.query_ball_point(potential_clusters_midpoints,r=distance_threshold)
return
for pair in pairs:
rotation, uv1, uv2 = calculate_angle_between_vectors(euler_angles[pair[0]], euler_angles[pair[1]],[1,0,0])

if rotation < orientation_threshold:
average_v = (uv1 + uv2) / np.linalg.norm(uv1 + uv2)
plane = Plane(Point3D((coordinates[pair[0]]+coordinates[pair[1]])/2), normal_vector=average_v)
potential_clusers.append((pair[0],pair[1],plane))
potential_clusters = pairs[angles < orientation_threshold]
print(f"There are {len(potential_clusters)} pairs with aligned orientation")
potential_clusters_midpoints = (coordinates[potential_clusters[:,0]] + coordinates[potential_clusters[:,1]]) / 2
potential_clusters_normal = vectors[potential_clusters[:,0]] + vectors[potential_clusters[:,1]]
potential_clusters_normal /= np.linalg.norm(potential_clusters_normal, axis=1)[:,None]
potential_additional_candidates = tree.query_ball_point(potential_clusters_midpoints,r=distance_threshold*2)

for i,potential_cluster in enumerate(potential_clusers):
pair_1, pair_2, plane = potential_cluster
res = tree.query_ball_point(plane.p1,r=distance_threshold)
add_points = []
for potential_additional_candidate in res:
if potential_additional_candidate in [pair_1, pair_2]:
continue
if plane.distance(Point3D(coordinates[potential_additional_candidate])).evalf() < 400 and ( calculate_angle_between_vectors(euler_angles[pair_1], euler_angles[potential_additional_candidate],[1,0,0])[0] < orientation_threshold or calculate_angle_between_vectors(euler_angles[pair_2], euler_angles[potential_additional_candidate],[1,0,0])[0] < orientation_threshold ):
add_points.append(potential_additional_candidate)
print(f"For pc {i} there are {len(add_points)} additional candidates")

if len(add_points) > 1:
data.loc[pair_1,'NUMCLUSTERS'] += 1
data.loc[pair_2,'NUMCLUSTERS'] += 1
for add_point in add_points:
data.loc[add_point, 'NUMCLUSTERS'] += 1
valid = 0
for i,potential_cluster in enumerate(potential_additional_candidates):
potential_cluster = np.array(potential_cluster)
dotproducts = np.sum(potential_clusters_normal[i] * vectors[potential_cluster], axis=1)
angles = np.degrees(np.arccos(dotproducts))
potential_cluster = potential_cluster[angles < orientation_threshold]

distance_to_plane = np.abs(np.sum((coordinates[potential_cluster]-potential_clusters_midpoints[i]) * potential_clusters_normal[i], axis=1))
potential_cluster = potential_cluster[distance_to_plane < 400]
if len(potential_cluster) > 3:

data.iloc[potential_cluster, data.columns.get_loc('CLUSTER_SCORE')] += len(potential_cluster)
print(data['CLUSTER_SCORE'].describe())
return
output_path = input_dir / f"{aa.area_name}_{ctx.obj.match_template_job.run_name}_{ctx.obj.match_template_job.run_id}_filtered_numclusters.star"
starfile.write(data, output_path, overwrite=True)

Expand Down
70 changes: 68 additions & 2 deletions src/decolace/processing/cli_match_template.py
Original file line number Diff line number Diff line change
Expand Up @@ -295,11 +295,20 @@ def join_matches(
new_matches = starfile.read(filtered_matches_starfile)
if len(new_matches) > 0:
matches.append(new_matches)
combined_matches = pd.concat(matches)
combined_matches = pd.concat(matches, ignore_index=True)
combined_matches_starfile = Path(ctx.obj.project.project_path) / "Matches" / f"combined_{ctx.obj.match_template_job.run_name}_{ctx.obj.match_template_job.run_id}_{name}.star"
if use_other_images != "":
combined_matches_starfile = combined_matches_starfile.with_name(combined_matches_starfile.stem + f"_{use_other_images}.star")
combined_matches['cisTEMOriginalImageFilename'].str.replace(".mrc", f"_{use_other_images}.mrc")
for i, row in combined_matches.iterrows():
if type(row['cisTEMOriginalImageFilename']) is float:
continue
new_filename = Path(row['cisTEMOriginalImageFilename'].strip("'").replace("_auto", f"_auto_{use_other_images}"))
if not new_filename.exists():
new_filename = str(new_filename).replace(".mrc", f"_{int(str(new_filename).split('_')[-3])-1}.mrc")
if not Path(new_filename).exists():
print(f"Ouch {new_filename} does not exist")
return
combined_matches.iloc[i,combined_matches.columns.get_loc('cisTEMOriginalImageFilename')] = "'"+str(new_filename)+"'"
if use_different_pixel_size is not None:
combined_matches['cisTEMPixelSize'] = use_different_pixel_size
starfile.write(combined_matches, combined_matches_starfile, overwrite=True)
Expand Down Expand Up @@ -340,3 +349,60 @@ def calculate_changes_during_refine_template(
refined_matches.loc[j, "y_change"] = original_match["cisTEMOriginalYPosition"] - refined_match["cisTEMOriginalYPosition"]
filtered_matches_starfile = Path(aa.cistem_project).parent / "Assets" / "TemplateMatching" / f"{aa.area_name}_{ctx.obj.match_template_job.run_id}_tm_package_filtered.star"
starfile.write(refined_matches, filtered_matches_starfile, overwrite=True)

@app.command()
def convert_coordinates(
starfile_path: Path = typer.Argument(..., help="Path to starfile"),
cistem_project: Path = typer.Argument(..., help="Path to cistem project"),
):
import starfile
import pandas as pd
import sqlite3
import mrcfile

data = starfile.read(starfile_path)

current_project = None
current_db = None

current_micrograph = None
current_micrograph_info = None
original_micrograph_info = None
current_micrograph_header = None
original_micrograph_header = None
for i, row in data.iterrows():
project = Path(row["cisTEMOriginalImageFilename"].strip("'")).parent.parent.parent
cc = "'"
project = project / (str(project.name) +'.db')
if project != current_project:
current_project = project
current_db = sqlite3.connect(current_project)
if current_micrograph != row["cisTEMOriginalImageFilename"].strip(cc):
current_micrograph = row["cisTEMOriginalImageFilename"].strip(cc)
current_micrograph_info = pd.read_sql(f'SELECT * FROM MOVIE_ALIGNMENT_LIST WHERE OUTPUT_FILE="{current_micrograph}"', current_db).iloc[-1]
with mrcfile.open(current_micrograph) as mrc:
current_micrograph_header = mrc.header
original_micrograph_info = pd.read_sql(f'SELECT * FROM MOVIE_ALIGNMENT_LIST WHERE ALIGNMENT_JOB_ID=1 AND MOVIE_ASSET_ID="{current_micrograph_info["MOVIE_ASSET_ID"]}"', current_db).iloc[-1]
with mrcfile.open(original_micrograph_info["OUTPUT_FILE"].strip(cc)) as mrc:
original_micrograph_header = mrc.header
#print(f' current: {current_micrograph_info["CROP_CENTER_X"]} {current_micrograph_info["CROP_CENTER_Y"]}')
#print(f' current: {current_micrograph_header["nx"]} {current_micrograph_header["ny"]}')
extraction_coordinate_original_x = row["cisTEMOriginalXPosition"] / 2.0 - original_micrograph_header["nx"] / 2.0
extraction_coordinate_original_y = row["cisTEMOriginalYPosition"] / 2.0 - original_micrograph_header["ny"] / 2.0
extraction_coordinate_original_unbinned_x = extraction_coordinate_original_x + original_micrograph_info["CROP_CENTER_X"]
extraction_coordinate_original_unbinned_y = extraction_coordinate_original_y + original_micrograph_info["CROP_CENTER_Y"]
extraction_coordinate_current_unbinned_x = extraction_coordinate_original_unbinned_x * 2.0
extraction_coordinate_current_unbinned_y = extraction_coordinate_original_unbinned_y * 2.0
extraction_coordinate_current_x = extraction_coordinate_current_unbinned_x - current_micrograph_info["CROP_CENTER_X"]
extraction_coordinate_current_y = extraction_coordinate_current_unbinned_y - current_micrograph_info["CROP_CENTER_Y"]
extraction_coordinate_A_x = extraction_coordinate_current_x + current_micrograph_header["nx"] / 2.0
extraction_coordinate_A_y = extraction_coordinate_current_y + current_micrograph_header["ny"] / 2.0
data.loc[i, "cisTEMOriginalXPosition"] = extraction_coordinate_A_x
data.loc[i, "cisTEMOriginalYPosition"] = extraction_coordinate_A_y
print(f' X-difference: {extraction_coordinate_A_x - row["cisTEMOriginalXPosition"]}, Y-difference: {extraction_coordinate_A_y - row["cisTEMOriginalYPosition"]}')
#print(f' original: {original_micrograph_info["CROP_CENTER_X"]} {original_micrograph_info["CROP_CENTER_Y"]}')
#print(f' original: {original_micrograph_header["nx"]} {original_micrograph_header["ny"]}')
output_filename = starfile_path.parent / (starfile_path.stem + "_converted.star")
starfile.write(data, output_filename, overwrite=True)


5 changes: 5 additions & 0 deletions src/decolace/processing/cli_processing.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
from decolace.processing.cli_preprocessing import app as preprocessing_app
from decolace.processing.cli_montaging import app as montaging_app
from decolace.processing.cli_match_template import app as match_template_app
from decolace.processing.cli_single_particle import app as single_particle_app

class OrderCommands(TyperGroup):
def list_commands(self, ctx: typer.Context):
Expand All @@ -35,6 +36,10 @@ def list_commands(self, ctx: typer.Context):
for command in match_template_app.registered_commands:
command.rich_help_panel="Template Matching Commands"
app.registered_commands += match_template_app.registered_commands
for command in single_particle_app.registered_commands:
command.rich_help_panel="Single Particle Commands"
app.registered_commands += single_particle_app.registered_commands


@app.callback()
def main(
Expand Down
48 changes: 48 additions & 0 deletions src/decolace/processing/cli_single_particle.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
import typer
import logging
from rich.logging import RichHandler
from rich import print
from pathlib import Path


app = typer.Typer()

@app.command()
def split_particles_into_optical_groups(
ctx: typer.Context,
tmpackage_star_file: Path,
parameters_star_file: Path,
split_by_session: bool = True,
beam_image_shift_threshold: float = 1.0
):
import starfile
import pandas as pd
import sqlite3
from sklearn.cluster import KMeans

particle_info = starfile.read(tmpackage_star_file)

aa_fd = pd.DataFrame([aa.model_dump() for aa in ctx.obj.acquisition_areas])


aa_fd["image_folder"] = aa_fd["cistem_project"].map(lambda x: str(Path(x).parent / "Assets" / "Images" ))
aa_fd["session"] = aa_fd["frames_folder"].map(lambda x: str(Path(x).parent.parent.parent.name))
particle_info["image_folder"] = particle_info["cisTEMOriginalImageFilename"].str.strip("'").map(lambda x: str(Path(x).parent))

particle_info["cisTEMOriginalImageFilename"] = particle_info["cisTEMOriginalImageFilename"].str.strip("'")
particle_info = particle_info.join(aa_fd.set_index("image_folder"), on="image_folder", rsuffix="_aa")
for i, cistem_project in enumerate(particle_info["cistem_project"].unique()):
db = sqlite3.connect(cistem_project)
microgaph_info = pd.read_sql_query("SELECT IMAGE_ASSETS.FILENAME, IMAGE_SHIFT_X, IMAGE_SHIFT_Y FROM IMAGE_ASSETS INNER JOIN MOVIE_ASSETS ON IMAGE_ASSETS.PARENT_MOVIE_ID=MOVIE_ASSETS.MOVIE_ASSET_ID INNER JOIN MOVIE_ASSETS_METADATA ON MOVIE_ASSETS.MOVIE_ASSET_ID=MOVIE_ASSETS_METADATA.MOVIE_ASSET_ID", db)
if i == 0:
particle_info = particle_info.merge(microgaph_info, right_on="FILENAME", left_on="cisTEMOriginalImageFilename",suffixes=(None,None),how="left")
particle_info.set_index("cisTEMOriginalImageFilename", inplace=True)
else:
particle_info.update(microgaph_info.set_index("FILENAME"))

particle_info["image_shift_label"] = particle_info.apply(lambda x: f"{x['session']}_{x['IMAGE_SHIFT_X'] // 1}_{x['IMAGE_SHIFT_Y'] // 1}", axis=1)
particle_info.reset_index(inplace=True)
print(particle_info["image_shift_label"].value_counts())



0 comments on commit 335ca73

Please sign in to comment.