Merge branch 'main' of https://github.com/jojoelfe/decolace

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

@@ -16,6 +16,8 @@ def run_matchtemplate(
     defocus_step: float = typer.Option(0.0, help="Defocus step for template matching"),
     defocus_range: float = typer.Option(0.0, help="Defocus range for template matching"),
     save_mip: bool = typer.Option(False, help="Save MIP of template"),
+    symmetry: str = typer.Option("C1", help="Symmetry of the template"),
+    run_on: str = typer.Option("all", help="Run on all or a subset of acquisition areas"),
 ):
     """Runs match template on all images in the acquisition areas"""
     from decolace.processing.project_managment import MatchTemplateRun
@@ -25,7 +27,7 @@ def run_matchtemplate(
     from pycistem.database import get_already_processed_images
     import pandas as pd
     logging.basicConfig(
-        level=logging.DEBUG,
+        level=logging.INFO,
         format="%(message)s",
         handlers=[
             RichHandler(),
@@ -45,7 +47,9 @@ def run_matchtemplate(
         "milano": 8,
         "sofia": 8,
         "manchester": 8,
-    }   
+    }  
+    if run_on != "all":
+        run_profile = {run_on: run_profile[run_on]}
 
     new_mtr = True
     if ctx.obj.match_template_job is None:
@@ -65,6 +69,7 @@ def run_matchtemplate(
         defocus_step = ctx.obj.match_template_job.defocus_step
         defocus_range = ctx.obj.match_template_job.defocus_range
         match_template_job_id = ctx.obj.match_template_job.run_id
+        symmetry = ctx.obj.match_template_job.symmetry
         typer.echo(f"Template match job id already exists, continuing job {ctx.obj.match_template_job.run_name}")
         typer.echo(f"template_filename={template.absolute().as_posix()}, angular_step={angular_step}, in_plane_angular_step={in_plane_angular_step} defous_step={defocus_step}, defocus_range={defocus_range}, decolace=True)")
 
@@ -78,6 +83,7 @@ def run_matchtemplate(
                 in_plane_angular_step=in_plane_angular_step,
                 defocus_step=defocus_step,
                 defocus_range=defocus_range,
+                symmetry=symmetry
             )
         )
         ctx.obj.project.write()
@@ -94,6 +100,7 @@ def run_matchtemplate(
             par.defocus_step = defocus_step
             par.defocus_search_range = defocus_range
             par.max_threads = 2
+            par.my_symmetry = symmetry
             if save_mip:
                 par.mip_output_file = par.scaled_mip_output_file.replace("_scaled_mip.mrc", "_mip.mrc")
 
@@ -102,7 +109,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")
@@ -117,6 +124,8 @@ def run_matchtemplate(
 @app.command()
 def create_tmpackage(
     ctx: typer.Context,
+    force_incomplete: bool = typer.Option(False, help="Force creation of tm package even if not all images have been processed"),
+    swap_phi_and_psi: bool = typer.Option(False, help="Swap phi and psi in tm package"),
 ):
     """Creates a Template-Matches Package for each acquisition area (i.e. a star file containing all matches)"""
     from pycistem.database import get_num_already_processed_images, write_match_template_to_starfile, insert_tmpackage_into_db, get_num_images
@@ -125,14 +134,16 @@ def create_tmpackage(
     for aa in ctx.obj.acquisition_areas:
         if get_num_already_processed_images(aa.cistem_project, ctx.obj.match_template_job.run_id) != get_num_images(aa.cistem_project):
             typer.echo(f"No images processed for {aa.area_name}")
-            continue
+            typer.echo(f"Processed {get_num_already_processed_images(aa.cistem_project, ctx.obj.match_template_job.run_id)} out of {get_num_images(aa.cistem_project)}")
+            if not force_incomplete:
+                continue
 
         typer.echo(f"Creating tm package for {aa.area_name}")
         output_filename = Path(aa.cistem_project).parent / "Assets" / "TemplateMatching" / f"{aa.area_name}_{ctx.obj.match_template_job.run_id}_tm_package.star"
         if output_filename.exists():
             typer.echo(f"Package already exists for {aa.area_name}")
             continue
-        write_match_template_to_starfile(aa.cistem_project, ctx.obj.match_template_job.run_id, output_filename)  
+        write_match_template_to_starfile(aa.cistem_project, ctx.obj.match_template_job.run_id, output_filename, switch_phi_psi=swap_phi_and_psi)  
         insert_tmpackage_into_db(aa.cistem_project, f"DeCOoLACE_{aa.area_name}_TMRun_{ctx.obj.match_template_job.run_id}", output_filename)
 
 @app.command()
@@ -153,7 +164,7 @@ def run_refinetemplate(
             input_starfile=tm_package_file.as_posix(),
             output_starfile=tm_package_file.with_suffix('').as_posix()+'_refined.star',
             input_template=Path(ctx.obj.match_template_job.template_path).as_posix(),
-            num_threads=10
+            num_threads=10,
         )
         if Path(par.output_starfile).exists():
             typer.echo(f"Refined tm package already exists for {aa.area_name}")
@@ -293,7 +304,7 @@ def join_matches(
         filtered_matches_starfile = Path(aa.cistem_project).parent / "Assets" / "TemplateMatching" / f"{aa.area_name}_{ctx.obj.match_template_job.run_id}_tm_package_filtered_{use_filtered}.star"
         if filtered_matches_starfile.exists():
             new_matches = starfile.read(filtered_matches_starfile)
-            if len(new_matches) > 0:
+            if len(new_matches) > 10:
                 matches.append(new_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"

src/decolace/processing/cli_montaging.py

@@ -1,58 +1,52 @@
 import typer
 from pathlib import Path
-from typing import List
+from typing import List, Optional
 
 app = typer.Typer()
 
 @app.command()
 def reset_montage(
-    project_main: Path = typer.Option(None, help="Path to wanted project file"),
+    ctx: typer.Context,
 ):
-    if project_main is None:
-        project_path = Path(glob.glob("*.decolace")[0])
-    project = ProcessingProject.read(project_path)
-
-    for aa in project.acquisition_areas:
+    for aa in ctx.obj.acquisition_areas:
         aa.montage_image = None
         aa.montage_star = None
         aa.initial_tile_star = None
         aa.refined_tile_star = None
-    project.write()
+    ctx.obj.project.write()
 
 @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"),
+    max_mean_density: Optional[float] = typer.Option(None, help="Maximum mean density of the tiles"),
+    cc_cutoff_as_fraction_of_median: float = typer.Option(0.5, help="Cutoff for the cross-correlation as a fraction of the median cross-correlation"),
 ):
     from rich.console import Console
     import starfile
     from numpy.linalg import LinAlgError
+    import pandas as pd
 
     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, prune_tiles
     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."
@@ -66,16 +60,21 @@ def run_montage(
 
             output_directory.mkdir(parents=True, exist_ok=True)
             output_path = output_directory / f"{aa.area_name}_tile_metadata.star"
+            if type(decolace_data) is not dict:
+                decolace_data = dict(decolace_data)
             tile_metadata_result = create_tile_metadata(cistem_data, decolace_data, output_path)
-            tile_metadata = tile_metadata_result["tiles"] 
+            tile_metadata = tile_metadata_result["tiles"]
+            tile_metadata, message = prune_tiles(tile_metadata, max_mean_density=max_mean_density)
+            console.log(message) 
             aa.initial_tile_star = output_path
             typer.echo(f"Wrote tile metadata to {output_path}.")
         else:
             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}.")
-
-        if aa.refined_tile_star is None:
+
+        # I should sort by acquisition time here
+        if aa.refined_tile_star is None or redo:
             estimated_distance_threshold = np.median(
                 tile_metadata["tile_x_size"] * tile_metadata["tile_pixel_size"]
             )
@@ -87,11 +86,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,cc_cutoff_as_fraction_of_media=cc_cutoff_as_fraction_of_median)
+                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
@@ -128,16 +139,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_preprocessing.py

@@ -77,7 +77,13 @@ def run_unblur(
 @app.command()
 def run_ctffind(
     ctx: typer.Context,
+    cmd_prefix: str = typer.Option("", help="Prefix of run command"),
+    cmd_suffix: str = typer.Option("", help="Suffix of run command"),
     num_cores: int = typer.Option(10, help="Number of cores to use"),
+    fit_nodes: bool = typer.Option(True, help="Fit nodes"),
+    fit_nodes_brute_force: bool = typer.Option(True, help="Fit nodes brute force"),
+    fit_nodes_lowres: float = typer.Option(30.0, help="Fit nodes lowres"),
+    fit_nodes_highres: float = typer.Option(4.0, help="Fit nodes highres"),
 ):
     """
     Run ctffind for each acquisition area
@@ -91,7 +97,12 @@ def run_ctffind(
         typer.echo(f"Running ctffind for {aa.area_name}")
         pars, image_info = ctffind.parameters_from_database(aa.cistem_project,decolace=True)
 
-        res = ctffind.run(pars,num_procs=num_cores)
+        for par in pars:
+            par.fit_nodes = fit_nodes
+            par.fit_nodes_1D_brute_force = fit_nodes_brute_force
+            par.fit_nodes_low_resolution_limit = fit_nodes_lowres
+            par.fit_nodes_high_resolution_limit = fit_nodes_highres
+        res = ctffind.run(pars,num_procs=num_cores,cmd_prefix=cmd_prefix,cmd_suffix=cmd_suffix)
 
         ctffind.write_results_to_database(aa.cistem_project,pars,res,image_info)
         aa.ctffind_run = True
@@ -102,7 +113,8 @@ def update_database(
     project_main: Path = typer.Option(None, help="Path to wanted project file")
 ):
     from pycistem.core import Project
-
+    import glob 
+    from decolace.processing.project_managment import ProcessingProject
     if project_main is None:
         project_path = Path(glob.glob("*.decolace")[0])
     project = ProcessingProject.read(project_path)