Better montage processing

src/decolace/processing/cli_edittags.py

@@ -0,0 +1,40 @@
+import typer
+app = typer.Typer()
+
+
+@app.command()
+def edit_experimental(
+    ctx: typer.Context,
+):
+    try:
+        from fastapi import FastAPI
+        from starlette.staticfiles import StaticFiles        
+        import uvicorn
+        import os 
+        from pathlib import Path
+        from decolace.processing.project_managment import AcquisitionAreaPreProcessing, process_experimental_conditions
+        from pydantic import BaseModel
+        class AcquisitionAreaExperimentalConditions(BaseModel):
+            area_name: str
+            experimental_conditions: dict
+        current_directory = Path(__file__).parent.absolute()
+        api = FastAPI()
+        @api.get("/aa")
+        async def get_aa() -> list[AcquisitionAreaExperimentalConditions]:
+            experiment_consitions_dicts = process_experimental_conditions(ctx.obj.acquisition_areas)
+            return [AcquisitionAreaExperimentalConditions(area_name=aa.area_name, experimental_conditions=experiment_consitions_dicts[i]) for i, aa in enumerate(ctx.obj.acquisition_areas)]
+        api.mount("/", StaticFiles(directory=current_directory/'web_components/acquisition_area_spreadsheet/',html=True), name="static")
+
+
+        # serve the ctx.aa object using FastAPI
+
+
+
+
+
+        uvicorn.run(api, host="127.0.0.1", port=8000, log_level="info")
+    except ModuleNotFoundError as e:
+        print(e)
+        typer.echo("FastAPI not installed. Please install using 'pip install fastapi'.")
+
+

src/decolace/processing/cli_match_template.py

@@ -102,7 +102,7 @@ def run_matchtemplate(
         all_image_info.append(image_info)
         typer.echo(f"{len(image_info)} tiles out of {orig_num} still to process")
         if len(image_info) == 0:
-            typer.echo(f"All tiles already processed")
+            typer.echo("All tiles already processed")
             continue
     all_image_info = pd.concat(all_image_info)
     typer.echo(f"Total of {len(all_image_info)} tiles to process")

src/decolace/processing/cli_montaging.py

@@ -21,41 +21,38 @@ def reset_montage(
 
 @app.command()
 def run_montage(
-    project_main: Path = typer.Option(None, help="Path to wanted project file"),
-    acquisition_areas: List[int] = typer.Option(None, help="List of acquisition areas to run montage for"),
-    iterations: int = typer.Option(2, help="Number of iterations"),
+    ctx: typer.Context,
+    iterations: int = typer.Option(3, help="Number of iterations"),
     num_procs: int = typer.Option(10, help="Number of processors to use"),
     bin_to: float = typer.Option(10.0, help="Pixel size to bin to"),
-    erode_mask_shifts: int = typer.Option(0, help="Erode mask by this many pixels before calculating cross-correlation"),
-    erode_mask_montage: int = typer.Option(0, help="Erode mask by this many pixels before assembling montage"), 
+    erode_mask_shifts: int = typer.Option(50, help="Erode mask by this many pixels before calculating cross-correlation"),
+    erode_mask_montage: int = typer.Option(50, help="Erode mask by this many pixels before assembling montage"), 
     filter_cutoff_frequency_ratio: float = typer.Option(0.02, help="Cutoff ratio of the butterworth filter"), 
-    ilter_order: float = typer.Option(4.0, help="Order of the butterworth filter"), 
+    filter_order: float = typer.Option(4.0, help="Order of the butterworth filter"), 
     mask_size_cutoff: int = typer.Option(100, help="If the mask size is smaller than this, the pair will be discarded"),
-    overlap_ratio: float = typer.Option(0.1, help="Overlap ratio parameter for masked crosscorrelation"),
+    overlap_ratio: float = typer.Option(0.2, help="Overlap ratio parameter for masked crosscorrelation"),
+    redo: bool = typer.Option(False, help="Redo the montage even if it already exists"),
+    redo_montage: bool = typer.Option(False, help="Redo only the creatin of the montage even if it already exists"),
 ):
     from rich.console import Console
     import starfile
     from numpy.linalg import LinAlgError
+    import pandas as pd
     import glob
     from decolace.processing.project_managment import ProcessingProject
 
 
     console = Console()
 
-    from decolace.processing.decolace_processing import read_data_from_cistem, read_decolace_data, create_tile_metadata, find_tile_pairs, calculate_shifts, calculate_refined_image_shifts, calculate_refined_intensity, create_montage_metadata, create_montage
+    from decolace.processing.decolace_processing import read_data_from_cistem, read_decolace_data, create_tile_metadata, find_tile_pairs, calculate_shifts, calculate_refined_image_shifts, calculate_refined_intensity, create_montage_metadata, create_montage, prune_bad_shifts
     import numpy as np
-    if project_main is None:
-        project_path = Path(glob.glob("*.decolace")[0])
-    project = ProcessingProject.read(project_path)
-    print(acquisition_areas)
-    for i, aa in enumerate(project.acquisition_areas):
-        if aa.montage_image is not None:
-            continue
-        if len(acquisition_areas) > 0 and i not in acquisition_areas:
+
+    for i, aa in enumerate(ctx.obj.acquisition_areas):
+        if aa.montage_image is not None and not redo and not redo_montage:
             continue
-        output_directory = project.project_path / "Montages" / aa.area_name
+        output_directory = ctx.obj.project.project_path / "Montages" / aa.area_name
         typer.echo(f"Running montage for {aa.area_name}")
-        if aa.initial_tile_star is None:
+        if aa.initial_tile_star is None or redo:
             cistem_data = read_data_from_cistem(aa.cistem_project)
             typer.echo(
                 f"Read data about {len(cistem_data)} tiles."
@@ -79,8 +76,10 @@ def run_montage(
             tile_metadata_result = starfile.read(aa.initial_tile_star,always_dict=True)
             tile_metadata = tile_metadata_result["tiles"]
             typer.echo(f"Read tile metadata from {aa.initial_tile_star}.")
+
+        # I should sort by acquisition time here
 
-        if aa.refined_tile_star is None:
+        if aa.refined_tile_star is None or redo:
             estimated_distance_threshold = np.median(
                 tile_metadata["tile_x_size"] * tile_metadata["tile_pixel_size"]
             )
@@ -92,11 +91,23 @@ def run_montage(
                 console.log(
                     f"In {len(tile_metadata)} tiles {len(tile_pairs)} tile pairs were found using distance threshold of {estimated_distance_threshold:.2f} A."
                 )
+                # I should sort by image_order here
                 # Calculate shifts of tile pairs
-                shifts = calculate_shifts(tile_pairs, num_procs, erode_mask=erode_mask_shifts, filter_cutoff_frequency_ratio=filter_cutoff_frequency_ratio, filter_order=ilter_order, mask_size_cutoff=mask_size_cutoff, overlap_ratio=overlap_ratio)
+                shifts = calculate_shifts(tile_pairs, num_procs, erode_mask=erode_mask_shifts, filter_cutoff_frequency_ratio=filter_cutoff_frequency_ratio, filter_order=filter_order, mask_size_cutoff=mask_size_cutoff, overlap_ratio=overlap_ratio)
                 console.log(
                     f"Shifts were adjusted. Mean: {np.mean(shifts['add_shift']):.2f} A, Median: {np.median(shifts['add_shift']):.2f} A, Std: {np.std(shifts['add_shift']):.2f} A. Min: {np.min(shifts['add_shift']):.2f} A, Max: {np.max(shifts['add_shift']):.2f} A."
                 )
+                # TODO: prune bad shifts and then tiles with no shifts
+                shifts, message = prune_bad_shifts(shifts)
+                shifts = shifts.copy()
+                console.log(message)
+
+                # Prune out tiles without any shifts
+                init_len = len(tile_metadata)
+                tile_metadata["filename_index"] = tile_metadata["tile_filename"]
+                tile_metadata.set_index("filename_index", inplace=True)
+                tile_metadata = tile_metadata.loc[pd.concat([shifts["image_1"],shifts["image_2"]]).unique()].copy()
+                console.log(f"Pruned out {init_len - len(tile_metadata)} tiles without any shifts.")
                 output_path = output_directory / f"{aa.area_name}_cc_shifts_{iteration}.star"
                 starfile.write(shifts, output_path, overwrite=True)
                 # Optimize tile positions
@@ -133,16 +144,26 @@ def run_montage(
         )
         aa.montage_star = output_path_metadata
         aa.montage_image = output_path_montage
-        project.write()
+        ctx.obj.project.write()
+
+
+from enum import Enum
+class KnotsOptions(str, Enum):
+    uniform = "uniform"
+    quantile = "quantile"
 
 @app.command()
 def filter_montage(
     ctx: typer.Context,
+    n_knots: int = typer.Option(5, help="Number of knots for spline"),
+    knots: KnotsOptions = typer.Option(KnotsOptions.uniform, help="Knots for spline"),
 ):
     for aa in ctx.obj.acquisition_areas:
         if aa.montage_image is None:
             continue
+        typer.echo(f"Filtering montage for {aa.area_name}")
         output_directory = ctx.obj.project.project_path / "Montages" / aa.area_name
         output_path = output_directory / f"{aa.area_name}_montage_filtered.mrc"
-        from decolace.processing.montage_filtering import highpassfilter_montage
-        highpassfilter_montage(aa.montage_image, output_path)
+        from decolace.processing.montage_filtering import subtract_linear_background_model
+        #highpassfilter_montage(aa.montage_image, output_path)
+        subtract_linear_background_model(aa.montage_image, output_path, n_knots=n_knots, knots=knots)

src/decolace/processing/cli_processing.py

@@ -7,13 +7,15 @@
 from rich.logging import RichHandler
 from typer.core import TyperGroup
 
-from decolace.processing.project_managment import ProcessingProject
+from decolace.processing.project_managment import ProcessingProject, DLContext, DLGlobals
 
 from decolace.processing.cli_project_managment import app as project_managment_app
 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
+from decolace.processing.cli_edittags import app as edittags_app
+from decolace.processing.cli_visualization import app as visualization_app
 
 class OrderCommands(TyperGroup):
   def list_commands(self, ctx: typer.Context):
@@ -39,11 +41,17 @@ def list_commands(self, ctx: typer.Context):
 for command in single_particle_app.registered_commands:
     command.rich_help_panel="Single Particle Commands"
 app.registered_commands += single_particle_app.registered_commands
+for command in edittags_app.registered_commands:
+    command.rich_help_panel="Edit Tags Commands"
+app.registered_commands += edittags_app.registered_commands
+for command in visualization_app.registered_commands:
+    command.rich_help_panel="Visualization Commands"
+app.registered_commands += visualization_app.registered_commands
 
 
 @app.callback()
 def main(
-    ctx: typer.Context,
+    ctx: DLContext,
     project: Path = typer.Option(None, help="Path to wanted project file",rich_help_panel="Expert Options"),
     acquisition_area_name: List[str] = typer.Option(None, help="List of acquisition areas names to process"),
     select_condition: str = typer.Option(None, help="Condition to select acquisition areas"),
@@ -72,4 +80,4 @@ def main(
         mtr = project_object.match_template_runs[array_position]
     else:
         mtr = None
-    ctx.obj = SimpleNamespace(project = project_object, acquisition_areas = aas_to_process, match_template_job = mtr, cistem_path = cistem_path)
+    ctx.obj = DLGlobals(project = project_object, acquisition_areas = aas_to_process, match_template_job = mtr, cistem_path = cistem_path)

src/decolace/processing/cli_project_managment.py

@@ -102,7 +102,6 @@ def status(
 ):
     from pycistem.database import get_num_movies, get_num_images, get_num_already_processed_images, get_num_matches
     from rich.table import Table
-    from rich.live import Live
     from rich import print
     from glob import glob
     import json
@@ -122,68 +121,78 @@ def status(
     table.add_column("Ctffind")
     table.add_column("Montage")
 
-    with Live(table,refresh_per_second=4,vertical_overflow="visible") as live:
 
-        for i, aa in enumerate(ctx.obj.acquisition_areas):
-            # Get size and modification time of the file aa.cistem_project
-            if aa.cistem_project is None:
-                table.add_row(
-                       str(i),
-                       aa.area_name,
-                       ":x:",
-                       ":x:",
-                       ":x:",
-                       ":x:")
-                continue
-            cistem_project_path = Path(aa.cistem_project)
-            if not cistem_project_path.exists():
-                table.add_row(
+    for i, aa in enumerate(ctx.obj.acquisition_areas):
+        # Get size and modification time of the file aa.cistem_project
+        if aa.cistem_project is None:
+            table.add_row(
                     str(i),
                     aa.area_name,
                     ":x:",
                     ":x:",
                     ":x:",
                     ":x:")
-                continue
-            cistem_project_size = str(cistem_project_path.stat().st_size)
-            cistem_project_mtime = str(cistem_project_path.stat().st_mtime)
-
-            if f"{aa.area_name}_preprocessing" in status_cache and status_cache[f"{aa.area_name}_preprocessing"][0] == cistem_project_size and status_cache[f"{aa.area_name}_preprocessing"][1] == cistem_project_mtime:
-                num_movies = status_cache[f"{aa.area_name}_preprocessing"][2]
-                num_images = status_cache[f"{aa.area_name}_preprocessing"][3]
-                ctffind_run = status_cache[f"{aa.area_name}_preprocessing"][4]
-                montage_run = status_cache[f"{aa.area_name}_preprocessing"][5]
-            else:
-                num_movies = f"✓ {get_num_movies(aa.cistem_project)}" if aa.cistem_project is not None else ":x:"
-                num_images = f"✓ {get_num_images(aa.cistem_project)}" if aa.unblur_run else ":x:"
-                ctffind_run = "✓" if aa.ctffind_run else ":x:"
-                montage_run = "✓" if aa.montage_image is not None else ":x:"
-                status_cache[f"{aa.area_name}_preprocessing"] = [cistem_project_size, cistem_project_mtime, num_movies, num_images, ctffind_run, montage_run]
-                json.dump(status_cache, open(status_chache_file, "w"))
-
+            continue
+        cistem_project_path = Path(aa.cistem_project)
+        if not cistem_project_path.exists():
             table.add_row(
                 str(i),
                 aa.area_name,
-                num_movies,
-                num_images,
-                ctffind_run,
-                montage_run
-            )
-        live.update(table)
+                ":x:",
+                ":x:",
+                ":x:",
+                ":x:")
+            continue
+        cistem_project_size = str(cistem_project_path.stat().st_size)
+        cistem_project_mtime = str(cistem_project_path.stat().st_mtime)
+
+        if f"{aa.area_name}_preprocessing" in status_cache and status_cache[f"{aa.area_name}_preprocessing"][0] == cistem_project_size and status_cache[f"{aa.area_name}_preprocessing"][1] == cistem_project_mtime:
+            num_movies = status_cache[f"{aa.area_name}_preprocessing"][2]
+            num_images = status_cache[f"{aa.area_name}_preprocessing"][3]
+            ctffind_run = status_cache[f"{aa.area_name}_preprocessing"][4]
+            montage_run = status_cache[f"{aa.area_name}_preprocessing"][5]
+        else:
+            num_movies = f"✓ {get_num_movies(aa.cistem_project)}" if aa.cistem_project is not None else ":x:"
+            num_images = f"✓ {get_num_images(aa.cistem_project)}" if aa.unblur_run else ":x:"
+            ctffind_run = "✓" if aa.ctffind_run else ":x:"
+            montage_run = "✓" if aa.montage_image is not None else ":x:"
+            status_cache[f"{aa.area_name}_preprocessing"] = [cistem_project_size, cistem_project_mtime, num_movies, num_images, ctffind_run, montage_run]
+            json.dump(status_cache, open(status_chache_file, "w"))
+
+        table.add_row(
+            str(i),
+            aa.area_name,
+            num_movies,
+            num_images,
+            ctffind_run,
+            montage_run
+        )
+    print(table)
 
     match_template_table = Table(title="Match Template Status")
     match_template_table.add_column("AA")
 
     for mtm in ctx.obj.project.match_template_runs:
         match_template_table.add_column(f"{mtm.run_id}: {mtm.run_name}")
-    with Live(match_template_table,refresh_per_second=4,vertical_overflow="visible") as live:
 
-        for aa in ctx.obj.acquisition_areas:
+    for aa in ctx.obj.acquisition_areas:
+        if aa.cistem_project is None:
+            continue
+        cistem_project_path = Path(aa.cistem_project)
+        if not cistem_project_path.exists():
+            continue
+        cistem_project_size = str(cistem_project_path.stat().st_size)
+        cistem_project_mtime = str(cistem_project_path.stat().st_mtime)
+
+        if f"{aa.area_name}_match_template" in status_cache and status_cache[f"{aa.area_name}_match_template"][0] == cistem_project_size and status_cache[f"{aa.area_name}_match_template"][1] == cistem_project_mtime:
+            mtm_status = status_cache[f"{aa.area_name}_match_template"][2]
+        else:
             mtm_status = []
             for mtm in ctx.obj.project.match_template_runs:        
                 mtm_status.append(f"{get_num_already_processed_images(aa.cistem_project, mtm.run_id)}/{get_num_images(aa.cistem_project)} {get_num_matches(aa.cistem_project, mtm.run_id)} Matches")
-
-            match_template_table.add_row(
-                aa.area_name, *mtm_status)
-
-        live.update(match_template_table)
+            status_cache[f"{aa.area_name}_match_template"] = [cistem_project_size, cistem_project_mtime, mtm_status]
+            json.dump(status_cache, open(status_chache_file, "w"))   
+        match_template_table.add_row(
+            aa.area_name, *mtm_status)
+
+    print(match_template_table)

src/decolace/processing/cli_visualization.py

@@ -0,0 +1,21 @@
+import typer
+from pathlib import Path
+from typing import List, Optional
+from decolace.processing.project_managment import DLContext
+
+app = typer.Typer()
+
+@app.command()
+def render_matches_via_blender(
+    ctx: DLContext,
+    filter_set_name: str = typer.Argument(..., help="Name of filter set"),
+):
+    """
+    Render matches via blender
+    """
+    from decolace.processing.match_visualization import render_aa
+
+    for aa in ctx.obj.project.acquisition_areas:
+        filename = ctx.obj.project.project_path / "visualization" / "matches" / f"{aa.area_name}_{ctx.obj.match_template_job.run_id}_tm_package_filtered_{filter_set_name}.png"
+        filename.parent.mkdir(parents=True, exist_ok=True)
+        render_aa(ctx.obj.project, aa, ctx.obj.match_template_job, filter_set_name, filename)