Skip to content

Commit

Permalink
Changes from 10062023 session
Browse files Browse the repository at this point in the history
  • Loading branch information
@jojoelfe
jojoelfe committed Oct 9, 2023
1 parent 77e3c50 commit d5b231b
Show file tree
Hide file tree
Showing 4 changed files with 143 additions and 90 deletions.
25 changes: 19 additions & 6 deletions src/decolace/acquisition/acquisition_area.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@
from contrasttransferfunction import CtfFit


def _hexagonal_cover(polygon, radius):
def _hexagonal_cover(polygon, radius, y_direction=1):
"""
Compute hexagonal grid covering the input polygon using spheres of the given radius.
Expand Down Expand Up @@ -53,6 +53,7 @@ def _hexagonal_cover(polygon, radius):

# Loop over each hexagon in the grid and test if it intersects the input polygon
for j in range(-ny, ny):
j *= y_direction
y = miny + offsety + (j * radius * 3 / 2)
direction = int(((j % 2) - 0.5) * 2)
for i in range(-nx*direction, nx*direction, direction):
Expand Down Expand Up @@ -151,11 +152,12 @@ class Config:
arbitrary_types_allowed = True

class AcquisitionAreaSingle:
state: AcquisitionAreaSingleState = AcquisitionAreaSingleState()


def __init__(self, name, directory, defocus=-1.0, tilt=0.0):
self.name = name
self.directory = directory
self.state:AcquisitionAreaSingleState = AcquisitionAreaSingleState()
Path(directory).mkdir(parents=True, exist_ok=True)
self.frames_directory = os.path.join(directory, "frames")
Path(self.frames_directory).mkdir(parents=True, exist_ok=True)
Expand Down Expand Up @@ -220,7 +222,7 @@ def initialize_from_napari(
)
self.state.corner_positions_image = np.array(corner_coordinates)

def calculate_acquisition_positions_from_napari(self, beam_radius, add_overlap=0.05, use_square_beam=False):
def calculate_acquisition_positions_from_napari(self, beam_radius, add_overlap=0.05, use_square_beam=False, start_from_bottom=False):

polygon = Polygon(self.state.corner_positions_specimen)

Expand All @@ -230,19 +232,22 @@ def calculate_acquisition_positions_from_napari(self, beam_radius, add_overlap=0
polygon, beam_radius * (1 - add_overlap)
)
else:
direction=1
if start_from_bottom:
direction=-1
center = _hexagonal_cover(
polygon, beam_radius * (1 - add_overlap)
polygon, beam_radius * (1 - add_overlap), y_direction=direction
)

self.state.acquisition_positions = center
self.state.positions_acquired = np.zeros(
(self.state.acquisition_positions.shape[0]), dtype="bool"
)

def predict_beamshift(self, specimen_coordinates):
def predict_beamshift(self, specimen_coordinates, selection=-100):
if self.state.beamshift_calibration_measurements is None:
return None
data = np.array(self.state.beamshift_calibration_measurements[-100:])
data = np.array(self.state.beamshift_calibration_measurements[selection:])
if len(data) < self.state.beamshift_calibrations_required:
return None
reg_x = HuberRegressor().fit(data[:, [0, 1]], data[:, 3])
Expand All @@ -261,6 +266,7 @@ def acquire(
):
serialem = connect_sem()
last_bs_correction=0.0
num_bad_predictions = 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)):
Expand Down Expand Up @@ -332,6 +338,13 @@ def acquire(
beam_shift_after_centering[1],
]
)
num_bad_predictions = 0
else:
if 'using_beamshift_prediction' in report and report['using_beamshift_prediction']:
num_bad_predictions += 1
if num_bad_predictions > 3:
self.state.beamshift_calibration_measurements = []


if counts < self.state.count_threshold_for_ctf:
if progress_callback is not None:
Expand Down
182 changes: 110 additions & 72 deletions src/decolace/acquisition/cli_acquisition.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,23 @@ def exit_gracefully(self, *args):
print("Recieved SIGTERM")
self.kill_now = True

class PauseProgress:
def __init__(self, progress: Progress) -> None:
self._progress = progress

def _clear_line(self) -> None:
UP = "\x1b[1A"
CLEAR = "\x1b[2K"
for _ in self._progress.tasks:
print(UP + CLEAR + UP)

def __enter__(self):
self._progress.stop()
self._clear_line()
return self._progress

def __exit__(self, exc_type, exc_value, exc_traceback):
self._progress.start()

app = typer.Typer()

Expand Down Expand Up @@ -382,7 +398,7 @@ def setup_areas(
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):
def my_widget(shapes: Shapes, use_square_beam: bool = False, start_from_bottom: bool = False):
points = []
areas = shapes.data
for area in areas:
Expand All @@ -392,8 +408,8 @@ def my_widget(shapes: Shapes, use_square_beam: bool = False):
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],affine=affine)
aa.calculate_acquisition_positions_from_napari(beam_radius=session_o.state.beam_radius, use_square_beam=use_square_beam)
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)
Expand All @@ -406,7 +422,18 @@ def my_widget(shapes: Shapes, use_square_beam: bool = False):
print("Done")
typer.Exit()

numbad = 0

@app.command()
def performance():
serialem = connect_sem()
import time
for i in range(1200):
t1 = time.perf_counter(), time.process_time()
serialem.ReportBinning("R")
t2 = time.perf_counter(), time.process_time()
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")

@app.command()
def acquire(
Expand All @@ -415,6 +442,7 @@ def acquire(
stepwise: bool = typer.Option(False, help="Acquire stepwise"),
defocus_offset: float = typer.Option(6.0)
):
from rich.prompt import Confirm
killer = GracefulKiller()
session_o = load_session(session_name, directory)
session_o.active_grid.state.stepwise = stepwise
Expand All @@ -424,78 +452,88 @@ def acquire(
for aa in session_o.active_grid.acquisition_areas
]
)

progress = Progress(auto_refresh=False)
grid_task = progress.add_task("grid", total=total_shots)
aa_task = progress.add_task("area")

progress.console.print(
Panel(
"[bold blue]Starting acquisition of grid {session_o.active_grid.name}.",
padding=1,
)
alreaddy_acquired = sum(
[
aa.state.positions_acquired.sum()
for aa in session_o.active_grid.acquisition_areas
]
)

def progress_callback(grid, report, acquisition_area, type=None):
if type == "start_new_area":
progress.log(
f"Starting acquisition of area {acquisition_area.name} in grid {grid.name}"
)
progress.reset(
aa_task,
total=len(acquisition_area.state.acquisition_positions),
description=f"{acquisition_area.name}",
)
progress.start_task(aa_task)
elif type == "resume_area":
progress.log(
f"Resuming acquisition of area {acquisition_area.name} in grid {grid.name}"
with Progress(auto_refresh=False) as progress:
grid_task = progress.add_task("grid", total=total_shots, completed=alreaddy_acquired)

progress.console.print(
Panel(
"[bold blue]Starting acquisition of grid {session_o.active_grid.name}.",
padding=1,
)
if report is not None:
progress.update(aa_task, advance=1)
progress.update(grid_task, advance=1)
log_string = f"{report['position']}/{len(acquisition_area.state.acquisition_positions)} "
if "using_focus_prediction" in report:
log_string += "FP :green_heart: "
else:
log_string += "FP :x: "
if "using_beamshift_prediction" in report:
log_string += "BSP :green_heart: "
else:
log_string += "BSP :x: "
if "fasttracked" in report:
log_string += "FT :green_heart: "
else:
log_string += "FT :x: "

log_string += f"Counts {report['counts']:.1f} "

if "beamshift_correction" in report:
log_string += f"BSC: {report['beamshift_correction']:.4f}um "

if "measured_defocus" in report:
log_string += f"MF: {report['measured_defocus']:.2f}um {report['ctf_cc']:.3f}CC {report['ctf_res']:.2f}A "

if "defocus_adjusted_by" in report:
log_string += f"DA: {report['defocus_adjusted_by']:.3f}um"

progress.log(log_string)
if killer.kill_now:
grid.state.stepwise = True
if grid.state.stepwise:
cont = typer.confirm("Continue?")
if not cont:
save = typer.confirm("Save?")
if save:
acquisition_area.write_to_disk()
continous = typer.confirm("Continous:")
if continous:
grid.state.stepwise = False
else:
print("Aborting!")
raise typer.Abort()

session_o.active_grid.start_acquisition(initial_defocus=session_o.state.fringe_free_focus_cross_grating-defocus_offset,progress_callback=progress_callback)
)
progress.start_task(grid_task)

def progress_callback(grid, report, acquisition_area, type=None):
global numbad
if type == "start_new_area":
progress.console.log(
f"Starting acquisition of area {acquisition_area.name} in grid {grid.name}"
)

elif type == "resume_area":
progress.console.log(
f"Resuming acquisition of area {acquisition_area.name} in grid {grid.name}"
)
if report is not None:
progress.update(grid_task, advance=1,refresh=True)
numbad += 1
log_string = f"{report['position']}/{len(acquisition_area.state.acquisition_positions)} "
if "using_focus_prediction" in report:
log_string += "FP :green_heart: "
else:
log_string += "FP :x: "
if "using_beamshift_prediction" in report:
log_string += "BSP :green_heart: "
else:
log_string += "BSP :x: "
if "fasttracked" in report:
log_string += "FT :green_heart: "
else:
log_string += "FT :x: "

log_string += f"Counts {report['counts']:.1f} "

if "beamshift_correction" in report:
log_string += f"BSC: {report['beamshift_correction']:.4f}um "

if "measured_defocus" in report:
log_string += f"MF: {report['measured_defocus']:.2f}um {report['ctf_cc']:.3f}CC {report['ctf_res']:.2f}A "

if "defocus_adjusted_by" in report:
numbad = 0
log_string += f"DA: {report['defocus_adjusted_by']:.3f}um"
if "perf" in report:
log_string += report["perf"]
if numbad > 0:
log_string += f" {numbad} failed corrections"

progress.console.log(log_string)
if killer.kill_now or numbad > 20:
grid.state.stepwise = True
if grid.state.stepwise:
with PauseProgress(progress):
cont = Confirm.ask("Continue?")
if not cont:
with PauseProgress(progress):
save = Confirm.ask("Save?")
if save:
acquisition_area.write_to_disk()
with PauseProgress(progress):
continous = Confirm.ask("Continous:")
if continous:
grid.state.stepwise = False
else:
print("Aborting!")
raise typer.Abort()

session_o.active_grid.start_acquisition(initial_defocus=session_o.state.fringe_free_focus_cross_grating-defocus_offset,progress_callback=progress_callback)

@app.callback()
def main(
Expand Down
11 changes: 8 additions & 3 deletions src/decolace/acquisition/grid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,6 +124,7 @@ def nice_view(self):

def take_map(self):
serialem = connect_sem()
serialem.SetDefocus(-10)
serialem.View()
self.ensure_view_file_is_open()
serialem.Save()
Expand All @@ -136,7 +137,11 @@ 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)
if np.sum(aa.state.positions_acquired) == len(
aa.state.positions_acquired
):
Expand All @@ -160,7 +165,7 @@ def start_acquisition(self, initial_defocus, progress_callback=None):
serialem.GoToLowDoseArea("R")
initial_beamshift = None
if index > 0:
initial_beamshift = self.acquisition_areas[index-1].predict_beamshift(aa.state.acquisition_positions[0])
initial_beamshift = self.acquisition_areas[index-1].predict_beamshift(aa.state.acquisition_positions[0],selection=0)


serialem.ManageDewarsAndPumps(-1)
Expand Down Expand Up @@ -201,8 +206,8 @@ def draw_acquisition_positions_into_napari(self, viewer, map_navids, beam_radius

positions[map_index].extend(pos)
order[map_index].extend(np.array(range(len(aa.state.acquisition_positions))))
pos = np.concatenate(positions, axis = 0)
order = np.concatenate(order, axis = 0)
pos = np.concatenate([a for a in positions if len(a) > 0], axis = 0)
order = np.concatenate([a for a in order if len(a) > 0], axis = 0)

if use_square_beam:
symbol='square'
Expand Down
15 changes: 6 additions & 9 deletions src/decolace/acquisition/session.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,15 +150,8 @@ def calculate_fringe_free_score(defocus):
minimal_slope = np.diff(ra).min()
return minimal_slope

from scipy.optimize import minimize_scalar

res = minimize_scalar(calculate_fringe_free_score,
bounds=(self.state.fringe_free_focus_vacuum-5.0, self.state.fringe_free_focus_vacuum+5.0),
method='bounded',
options={"maxiter":10,"disp":True})
self.state.fringe_free_focus_cross_grating = res.x
serialem.SetDefocus(res.x)
serialem.Record()
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()
Expand Down Expand Up @@ -225,6 +218,10 @@ def calculate_fringe_free_score(defocus):
bounds=(self.state.min_defocus_for_ffsearch, self.state.max_defocus_for_ffsearch),
method='bounded',
options={"maxiter":10,"disp":True})
res = minimize_scalar(calculate_fringe_free_score,
bounds=(res.x-10.0, res.x+10.0),
method='bounded',
options={"maxiter":10,"disp":True})
self.state.fringe_free_focus_vacuum = res.x
serialem.SetDefocus(res.x)
serialem.Record()
Expand Down

0 comments on commit d5b231b

Please sign in to comment.