From c6e68fc7615014e4780656cc1a8df478e5b40007 Mon Sep 17 00:00:00 2001 From: Johannes Elferich Date: Thu, 2 Nov 2023 09:45:31 -0400 Subject: [PATCH] Changes during datacollection 20231030 --- notest.txt | 9 ++ pyproject.toml | 3 +- src/decolace/acquisition/acquisition_area.py | 41 +++++- src/decolace/acquisition/cli_acquisition.py | 130 ++++++++++++++++++- src/decolace/acquisition/grid.py | 5 +- src/decolace/acquisition/session.py | 101 ++++++++++---- 6 files changed, 250 insertions(+), 39 deletions(-) create mode 100644 notest.txt diff --git a/notest.txt b/notest.txt new file mode 100644 index 0000000..2518bc1 --- /dev/null +++ b/notest.txt @@ -0,0 +1,9 @@ +- Better alignemnt of beam-image-beam_shift + - x : 5.6 um half gone + - y: 5.6 half gone + - Now using Kangkans FFI alignements + - So much better!!!!! +- Fix napari text bug - Done +- Improve serpentine +- Use previous BIS-calibration more +- Detect over-focus better \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index f45cb75..da21aa6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,7 +36,8 @@ dependencies = [ "mdocfile", "mrcfile", "contrasttransferfunction", - "pydantic" + "pydantic", + "tsp-solver2" ] [project.scripts] diff --git a/src/decolace/acquisition/acquisition_area.py b/src/decolace/acquisition/acquisition_area.py index edef9a0..aca6ce3 100644 --- a/src/decolace/acquisition/acquisition_area.py +++ b/src/decolace/acquisition/acquisition_area.py @@ -7,7 +7,7 @@ import numpy as np from sklearn.linear_model import HuberRegressor from timeit import default_timer as timer -from pydantic import BaseModel +from pydantic.v1 import BaseModel from typing import Optional from .serialem_helper import connect_sem @@ -128,7 +128,7 @@ class AcquisitionAreaSingleState(BaseModel): count_threshold_for_beamshift: float = 1500 count_threshold_for_ctf: float = 1500 - ctf_cc_threshold: float = 100 + ctf_cc_threshold: float = 70 ctf_step_when_unreliable: float = -0.03 ctf_max_step_initially: float = 2.0 ctf_max_step: float = 0.5 @@ -147,6 +147,7 @@ class AcquisitionAreaSingleState(BaseModel): positions_still_to_fasttrack: int = 0 fasttrack: bool = False + aborted: bool = False class Config: arbitrary_types_allowed = True @@ -182,6 +183,7 @@ def initialize_from_napari( self, map_navigator_id: int, center_coordinate, corner_coordinates, affine= None ): serialem = connect_sem() + serialem.SuspendNavRedraw() serialem.LoadOtherMap(map_navigator_id, "A") center_item = int( @@ -243,6 +245,21 @@ def calculate_acquisition_positions_from_napari(self, beam_radius, add_overlap=0 self.state.positions_acquired = np.zeros( (self.state.acquisition_positions.shape[0]), dtype="bool" ) + + def optimize_path(self): + from tsp_solver.greedy import solve_tsp + + # Create square distance matrix of self.state.acquisition_positions + dist_matrix = np.zeros((len(self.state.acquisition_positions), len(self.state.acquisition_positions))) + for i in range(len(self.state.acquisition_positions)): + for j in range(len(self.state.acquisition_positions)): + dist_matrix[i][j] = np.linalg.norm(self.state.acquisition_positions[i] - self.state.acquisition_positions[j]) + + # Solve TSP to find optimal path + path = solve_tsp(dist_matrix,endpoints=(0,len(self.state.acquisition_positions)-1)) + self.state.acquisition_positions = self.state.acquisition_positions[path] + + def predict_beamshift(self, specimen_coordinates, selection=-100): if self.state.beamshift_calibration_measurements is None: @@ -267,12 +284,16 @@ def acquire( serialem = connect_sem() last_bs_correction=0.0 num_bad_predictions = 0 + num_max_correction = 0 + offset_before_max_correction = 0 if self.state.acquisition_positions is None or self.state.positions_acquired is None: raise ValueError("No acquisition positions defined") for index in range(len(self.state.acquisition_positions)): report = {} if self.state.positions_acquired[index]: continue + if self.state.aborted: + return if index == 0: if initial_beamshift is not None and self.state.beamshift_correction: serialem.SetBeamShift(initial_beamshift[0], initial_beamshift[1]) @@ -281,6 +302,9 @@ def acquire( if index % 10 == 0: self.write_to_disk() + if established_lock: + serialem.SuppressReports() + serialem.DeferLogUpdates() report["position"] = index current_speciment_shift = np.array(serialem.ReportSpecimenShift()) diff = ( @@ -379,11 +403,24 @@ def acquire( self.state.desired_defocus = self.state.low_defocus + fraction_of_gradient * (self.state.high_defocus-self.state.low_defocus) offset = self.state.desired_defocus - measured_defocus + if abs(offset) > self.state.ctf_max_step and established_lock: offset = self.state.ctf_max_step * np.sign(offset) + if num_max_correction == 0: + offset_before_max_correction = abs(offset) + num_max_correction += 1 + else: if abs(offset) > self.state.ctf_max_step_initially: offset = self.state.ctf_max_step_initially * np.sign(offset) + if num_max_correction == 0: + offset_before_max_correction = abs(offset) + num_max_correction += 1 + else: + num_max_correction = 0 + if num_max_correction > 3 and abs(offset) > offset_before_max_correction: + report["potential_overfocus"] = True + if abs(offset) < 0.1 and not established_lock and last_bs_correction < 0.06: established_lock = True if abs(offset) > 0.001: diff --git a/src/decolace/acquisition/cli_acquisition.py b/src/decolace/acquisition/cli_acquisition.py index 8ec12e0..10f59dd 100644 --- a/src/decolace/acquisition/cli_acquisition.py +++ b/src/decolace/acquisition/cli_acquisition.py @@ -103,14 +103,24 @@ def prepare_beam_vacuum( typer.echo(f"Prepared beam vacuum for session {session_o.name}") @app.command() -def prepare_beam_cross( +def prepare_beam_cross_euc( 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(name, directory) - session_o.prepare_beam_cross() + session_o.prepare_beam_cross_euc() session_o.write_to_disk() - typer.echo(f"Prepared beam vacuum for session {session_o.name}") + typer.echo(f"Prepared beam cross 1 for session {session_o.name}") + +@app.command() +def prepare_beam_cross_final( + 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(name, directory) + session_o.prepare_beam_cross_final() + session_o.write_to_disk() + typer.echo(f"Prepared beam cross 2 for session {session_o.name}") @app.command() def print_session_state( @@ -120,6 +130,15 @@ def print_session_state( session_o = load_session(name, directory) print(session_o.state) +@app.command() +def print_aa_state( + 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(name, directory) + for aa in session_o.active_grid.acquisition_areas: + print(aa.state) + @app.command() def set_session_state( name: str = typer.Option(None, help="Name of the session"), @@ -207,6 +226,84 @@ def new_map( session_o.write_to_disk() typer.echo(f"Created new map for grid {session_o.active_grid.name}") +@app.command() +def update_cc( + 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) + session_o.state.min_defocus_for_ffsearch = -10.0 + session_o.state.max_defocus_for_ffsearch = 10.0 + #for aa in session_o.active_grid.acquisition_areas: + # aa.state.ctf_cc_threshold = 70 + # aa.write_to_disk() + session_o.write_to_disk() + +@app.command() +def optimize_path( + 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) + + for aa in session_o.active_grid.acquisition_areas: + aa.optimize_path() + serialem = connect_sem() + 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)) + + # Now test if there are already areas set up + existing_acquisiton_positions = False + if len(session_o.active_grid.state.acquisition_areas) > 0: + if sum([len(aa.state.acquisition_positions) for aa in session_o.active_grid.acquisition_areas]) > 0: + existing_acquisiton_positions = True + + if existing_acquisiton_positions: + affine = session_o.active_grid.draw_acquisition_positions_into_napari(viewer, map_navids, session_o.state.beam_radius) + + @magicgui(shapes={'label': 'Setup areas'}) + def my_widget(shapes: Shapes, use_square_beam: bool = False, start_from_bottom: bool = False): + 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(),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(beam_radius=session_o.state.beam_radius, use_square_beam=use_square_beam, start_from_bottom=start_from_bottom) + aa.write_to_disk() + session_o.active_grid.state.acquisition_areas.append([aa.name,aa.directory]) + session_o.active_grid.acquisition_areas.append(aa) + session_o.active_grid.write_to_disk() + session_o.active_grid.save_navigator() + session_o.active_grid.draw_acquisition_positions_into_napari(viewer, map_navids, session_o.state.beam_radius,use_square_beam=use_square_beam) + + viewer.window.add_dock_widget(my_widget) + napari.run() + print("Done") + typer.Exit() + @app.command() def add_acquisition_area(): @@ -435,14 +532,29 @@ def performance(): if i in [0,1,2,10,11,12,100,101,102,1000,1001,1002,9997,9998,9999]: print(f"Iteration {i}: {t2[0] - t1[0]:.4f} s") +def handle_problem(serialem, session_o: session, aa: AcquisitionAreaSingle, message: str): + from .session import ProblemPolicy + typer.echo(message) + if session_o.state.problem_policy == ProblemPolicy.CONTINUE: + return + elif session_o.state.problem_policy == ProblemPolicy.PAUSE: + session_o.active_grid.state.stepwise = True + elif session_o.state.problem_policy == ProblemPolicy.NEXT: + aa.state.aborted = True + else: + serialem.SetColumnOrGunValve(0) + typer.Exit() + + @app.command() def acquire( session_name: str = typer.Option(None, help="Name of the session"), directory: str = typer.Option(None, help="Directory to save session in"), stepwise: bool = typer.Option(False, help="Acquire stepwise"), - defocus_offset: float = typer.Option(6.0) + defocus_offset: float = typer.Option(2.0) ): from rich.prompt import Confirm + serialem = connect_sem() killer = GracefulKiller() session_o = load_session(session_name, directory) session_o.active_grid.state.stepwise = stepwise @@ -468,7 +580,7 @@ def acquire( padding=1, ) ) - progress.start_task(grid_task) + #progress.start_task(grid_task) def progress_callback(grid, report, acquisition_area, type=None): global numbad @@ -515,9 +627,14 @@ def progress_callback(grid, report, acquisition_area, type=None): log_string += f" {numbad} failed corrections" progress.console.log(log_string) - if killer.kill_now or numbad > 20: + if "potential_overfocus" in report and report["potential_overfocus"]: + handle_problem(serialem, session_o, acquisition_area, "Potential Overfocus detected") + if numbad > 20: + handle_problem(serialem,session_o,acquisition_area,"More than 20 bad predictions") + if killer.kill_now: grid.state.stepwise = True if grid.state.stepwise: + serialem.Exit(1) with PauseProgress(progress): cont = Confirm.ask("Continue?") if not cont: @@ -528,6 +645,7 @@ def progress_callback(grid, report, acquisition_area, type=None): with PauseProgress(progress): continous = Confirm.ask("Continous:") if continous: + killer.kill_now = False grid.state.stepwise = False else: print("Aborting!") diff --git a/src/decolace/acquisition/grid.py b/src/decolace/acquisition/grid.py index 04b506f..461741d 100644 --- a/src/decolace/acquisition/grid.py +++ b/src/decolace/acquisition/grid.py @@ -3,7 +3,7 @@ import time from functools import partial from pathlib import Path -from pydantic import BaseModel +from pydantic.v1 import BaseModel from typing import List, Optional from .serialem_helper import connect_sem @@ -138,7 +138,6 @@ def initialize_acquisition_areas(self, navigator_ids): def start_acquisition(self, initial_defocus, progress_callback=None): serialem = connect_sem() - for index, aa in enumerate(self.acquisition_areas): serialem.CloseLogOpenNew(1) #serialem.Exit(1) @@ -185,7 +184,7 @@ def draw_acquisition_positions_into_napari(self, viewer, map_navids, beam_radius 'string': '{order}', 'size': 10, 'color': 'white', - 'translation': np.array([0, 0]), + 'translation': np.array([0, 0, 0]), } order = [[] for a in map_navids] positions = [[] for a in map_navids] diff --git a/src/decolace/acquisition/session.py b/src/decolace/acquisition/session.py index ebc6e70..1e3fc0c 100644 --- a/src/decolace/acquisition/session.py +++ b/src/decolace/acquisition/session.py @@ -2,15 +2,20 @@ import os import time from pathlib import Path, PurePath -from pydantic import BaseModel +from pydantic.v1 import BaseModel from typing import List, Optional from .serialem_helper import connect_sem import numpy as np from rich import print - +from enum import Enum from .grid import grid +class ProblemPolicy(Enum): + CONTINUE = 'continue' + PAUSE = 'pause' + ABORT = 'abort' + NEXT = 'next' class SessionState(BaseModel): grids: List[List] = [] microscope_settings: dict = {} @@ -21,6 +26,11 @@ class SessionState(BaseModel): max_defocus_for_ffsearch: Optional[float] = 80 fringe_free_focus_cross_grating: Optional[float] = None dose_rate_e_per_pix_s: Optional[float] = None + unbinned_pixel_size_A: Optional[float] = None + cross_euc_Z_height: Optional[float] = None + cross_ff_Z_height: Optional[float] = None + euc_to_ff_offset: Optional[float] = None + problem_policy: ProblemPolicy = 'next' class Config: arbitrary_types_allowed = True @@ -103,9 +113,10 @@ def add_current_setting(self): self.state.microscope_settings[mode] = settings - def prepare_beam_cross(self): + def prepare_beam_cross_euc(self): serialem = connect_sem() - from contrasttransferfunction.spectrumhelpers import radial_average + serialem.Eucentricity(3) + self.state.cross_euc_Z_height = serialem.ReportStageXYZ()[2] serialem.SetBinning("R",4) serialem.SetExposure("R",1) serialem.SetDoseFracParams("R",0,0,0) @@ -123,35 +134,71 @@ def prepare_beam_cross(self): if center_tries > 10: print(f"Error could not center beam") return - def calculate_fringe_free_score(defocus): - serialem.SetDefocus(defocus) - serialem.Record() - serialem.CenterBeamFromImage(0, 1.0) + + print("Please adjust the focus until there are no beam fringes and then run prepare-beam-cross-final") + + def prepare_beam_cross_final(self): + from contrasttransferfunction import CtfFit + + serialem = connect_sem() + self.state.fringe_free_focus_cross_grating = serialem.ReportDefocus() + serialem.SetBinning("R",4) + serialem.SetExposure("R",1) + serialem.SetDoseFracParams("R",0,0,0) + + center_tries = 0 + adjustment = 99 + while adjustment > 0.01 * self.state.beam_radius: + beam_shift_before_centering = np.array(serialem.ReportBeamShift()) serialem.Record() serialem.CenterBeamFromImage(0, 1.0) - serialem.Record() + beam_shift_after_centering = np.array(serialem.ReportBeamShift()) + adjustment = np.linalg.norm(beam_shift_before_centering-beam_shift_after_centering) + print(f"Adjusting by {adjustment}") + center_tries+=1 + if center_tries > 10: + print(f"Error could not center beam") + return - beam_image = np.asarray(serialem.bufferImage("A")) - # Crop out largest possbile square from the center of slice - # Get the dimensions of the array - rows, cols = beam_image.shape + #serialem.MoveStage(0,0, -50.0) + wanted_defocus = -1.0 + + serialem.Record() + serialem.FFT("A") + powerspectrum = np.asarray(serialem.bufferImage("AF")) + fit_result = CtfFit.fit_1d( + powerspectrum, + pixel_size_angstrom=self.state.unbinned_pixel_size_A * 4, + voltage_kv=300.0, + spherical_aberration_mm=2.7, + amplitude_contrast=0.07) + measured_defocus = fit_result.ctf.defocus1_angstroms / -10000 + print(measured_defocus) + defocus_error = wanted_defocus - measured_defocus + num_tries=0 + while abs(defocus_error) > 0.1 and num_tries < 10: + num_tries += 1 + serialem.MoveStage(0,0,defocus_error) + serialem.Record() + serialem.FFT("A") + powerspectrum = np.asarray(serialem.bufferImage("AF")) + fit_result = CtfFit.fit_1d( + powerspectrum, + pixel_size_angstrom=self.state.unbinned_pixel_size_A * 4, + voltage_kv=300.0, + spherical_aberration_mm=2.7, + amplitude_contrast=0.07) + measured_defocus = fit_result.ctf.defocus1_angstroms / -10000 + print(measured_defocus) + defocus_error = wanted_defocus - measured_defocus + + self.state.cross_ff_Z_height = serialem.ReportStageXYZ()[2] + self.state.euc_to_ff_offset = self.state.cross_ff_Z_height - self.state.cross_euc_Z_height - # Calculate the size of the center square - size = min(rows, cols) + print(f"Measured stage Z offset between eucentric and fring-free position is {self.state.euc_to_ff_offset} um") - # Calculate the starting indices for the slice - start_row = (rows - size) // 2 - start_col = (cols - size) // 2 - # Slice the array to get the center square - center = beam_image[start_row:start_row+size, start_col:start_col+size] - ra = radial_average(center) - - minimal_slope = np.diff(ra).min() - return minimal_slope - for i in np.arange(-5.0,5.0,0.5): - print(calculate_fringe_free_score(self.state.fringe_free_focus_vacuum + i)) def prepare_beam_vacuum(self, coverage=0.9): serialem = connect_sem() @@ -164,7 +211,7 @@ def prepare_beam_vacuum(self, coverage=0.9): s_dim = min(x_dim,y_dim) s_dim_um = s_dim * pixel_size * 0.001 print(f"Pixel size is {pixel_size} at current magnification. Smallest camera dimension is {s_dim} pixels and {s_dim_um} um.") - + self.state.unbinned_pixel_size_A = pixel_size * 10 #Set binning to 4 and prevent saving of frames and alignment serialem.SetBinning("R",4) serialem.SetExposure("R",1)