diff --git a/PYME/Acquire/Protocols/htsms-calibration.py b/PYME/Acquire/Protocols/htsms-cal-psf.py similarity index 64% rename from PYME/Acquire/Protocols/htsms-calibration.py rename to PYME/Acquire/Protocols/htsms-cal-psf.py index cd82615e5..624648977 100644 --- a/PYME/Acquire/Protocols/htsms-calibration.py +++ b/PYME/Acquire/Protocols/htsms-cal-psf.py @@ -3,8 +3,17 @@ # T(when, what, *args) creates a new task. "when" is the frame number, "what" is a function to # be called, and *args are any additional arguments. taskList = [ - T(-1, scope.l642.TurnOn), - T(-1, scope.l560.TurnOn), + T(-1, scope.state.update, { + 'Lasers.MPB560.On': True, + 'Lasers.MPB560.Power': 0.0, + 'Lasers.MPB642.On': True, + 'Lasers.MPB642.Power': 0.0, + 'Lasers.OBIS405.On': False, + 'Lasers.OBIS488.On': False, + 'Multiview.ActiveViews': [0, 1, 2, 3], + 'Multiview.ROISize': [256, 256], + 'Camera.IntegrationTime': 0.1, + }), T(0, scope.focus_lock.DisableLock), T(maxint, scope.turnAllLasersOff), # T(maxint, scope.focus_lock.EnableLock), @@ -19,5 +28,5 @@ # must be defined for protocol to be discovered PROTOCOL = TaskListProtocol(taskList, metaData, preflight, filename=__file__) -PROTOCOL_STACK = ZStackTaskListProtocol(taskList, 1, 3, metaData, preflight, slice_order='triangle', +PROTOCOL_STACK = ZStackTaskListProtocol(taskList, 1, 3, metaData, preflight, slice_order='saw', require_camera_restart=True, filename=__file__) \ No newline at end of file diff --git a/PYME/Acquire/Protocols/htsms-cal-registration.py b/PYME/Acquire/Protocols/htsms-cal-registration.py new file mode 100644 index 000000000..3b544b997 --- /dev/null +++ b/PYME/Acquire/Protocols/htsms-cal-registration.py @@ -0,0 +1,52 @@ +from PYME.Acquire.protocol import * +from PYME.Acquire.Utils.pointScanner import PointScanner +import numpy as np +import logging + +logger = logging.getLogger(__name__) + +class Scanner(PointScanner): + def __init__(self, steps=10): + self._steps = steps + + def setup(self): + fs = np.array((scope.cam.size_x, scope.cam.size_y)) + fov_size = 2 * np.mean(fs * np.array(scope.GetPixelSize())) + + PointScanner.__init__(self, scope, self._steps, + fov_size/self._steps, 1, 0, + False, True, trigger=True, + stop_on_complete=True, + return_to_start=True) + self.on_stop.connect(scope.spoolController.StopSpooling) + self.genCoords() + +scanner = Scanner() + +taskList = [ + T(-1, scope.state.update, { + 'Lasers.MPB560.On': True, + 'Lasers.MPB560.Power': 0.0, + 'Lasers.MPB642.On': True, + 'Lasers.MPB642.Power': 0.0, + 'Lasers.OBIS405.On': False, + 'Lasers.OBIS488.On': False, + 'Multiview.ActiveViews': [0, 1, 2, 3], + 'Multiview.ROISize': [256, 256], + 'Camera.IntegrationTime': 0.1, + }), + T(-1, scanner.setup), + T(0, scanner.start), + T(maxint, scope.turnAllLasersOff), +] + +metaData = [ + ('Protocol.DataStartsAt', 0), +] + +preflight = [] # no preflight checks + +# must be defined for protocol to be discovered +PROTOCOL = TaskListProtocol(taskList, metaData, preflight, filename=__file__) +PROTOCOL_STACK = ZStackTaskListProtocol(taskList, 1, 3, metaData, preflight, slice_order='saw', + require_camera_restart=True, filename=__file__) diff --git a/PYME/Acquire/SpoolController.py b/PYME/Acquire/SpoolController.py index ae1dd27c9..666414bbe 100644 --- a/PYME/Acquire/SpoolController.py +++ b/PYME/Acquire/SpoolController.py @@ -414,7 +414,8 @@ def StartSpooling(self, fn=None, stack=None, compLevel=None, if stack: protocol = protocol_z protocol.dwellTime = z_dwell - #print(protocol) + stack_settings = stack if hasattr(stack, 'keys') else None + protocol.set_stack_positions(stack_settings) else: protocol = protocol diff --git a/PYME/Acquire/protocol.py b/PYME/Acquire/protocol.py index 056d6fc77..30bd41f5a 100644 --- a/PYME/Acquire/protocol.py +++ b/PYME/Acquire/protocol.py @@ -37,6 +37,7 @@ def __init__(self, filename=None): # The filename parameter exists to allow setting the filename in protocols which are not instantiated through the spool controller, and # requires passing __name__ to the constructor in the protocol itself. Both solutions are a bit gross, and may be revisited in the future. self.filename = filename + self.dwellTime = None def Init(self, spooler): pass @@ -185,6 +186,7 @@ def __init__(self, taskList, startFrame, dwellTime, metadataEntries=[], prefligh self.startFrame = startFrame self.dwellTime = dwellTime self.randomise = randomise + self.jitter = False if self.randomise: logger.warning("Use slice_order='random' instead of randomise=True") self.slice_order = 'random' @@ -192,11 +194,14 @@ def __init__(self, taskList, startFrame, dwellTime, metadataEntries=[], prefligh self.slice_order = slice_order self.require_camera_restart = require_camera_restart - - def Init(self, spooler): - self.zPoss = np.arange(scope.stackSettings.GetStartPos(), - scope.stackSettings.GetEndPos() + .95 * scope.stackSettings.GetStepSize(), - scope.stackSettings.GetStepSize() * scope.stackSettings.GetDirection()) + + def set_stack_positions(self, stack_mdh=None): + stack_mdh = dict() if stack_mdh is None else stack_mdh + start = stack_mdh.get('StackSettings.StartPos', scope.stackSettings.GetStartPos()) + stop = stack_mdh.get('StackSettings.EndPos', scope.stackSettings.GetEndPos()) + direction = 1 if stop > start else -1 + step = direction * stack_mdh.get('StackSettings.StepSize', scope.stackSettings.GetStepSize()) + self.zPoss = np.arange(start, stop + .95 * step, step) if self.slice_order != 'saw': if self.slice_order == 'random': @@ -208,9 +213,10 @@ def Init(self, spooler): else: # even self.zPoss = np.concatenate([self.zPoss[::2], self.zPoss[-1::-2]]) + + self.piezoName = 'Positioning.%s' % stack_mdh.get('StackSettings.ScanPiezo', scope.stackSettings.GetScanChannel()) - - self.piezoName = 'Positioning.%s' % scope.stackSettings.GetScanChannel() + def Init(self, spooler): self.startPos = scope.state[self.piezoName + '_target'] #FIXME - _target positions shouldn't be part of scope state self.pos = 0 @@ -225,13 +231,18 @@ def Init(self, spooler): TaskListProtocol.Init(self,spooler) def OnFrame(self, frameNum): + curr_cycle, self.jitter_val = 0, 0 if frameNum > self.startFrame: fn = int(floor((frameNum - self.startFrame)/self.dwellTime) % len(self.zPoss)) + cycle = int(floor((frameNum - self.startFrame)/self.dwellTime) // len(self.zPoss)) + if cycle != curr_cycle and self.jitter: + self.jitter_val = (np.random.rand(1) - 0.5) * (self.zPoss[1] - self.zPoss[0]) + curr_cycle = cycle if not fn == self.pos: self.pos = fn #self.piezo.MoveTo(self.piezoChan, self.zPoss[self.pos]) scope.state.setItem(self.piezoName, self.zPoss[self.pos], stopCamera=self.require_camera_restart) - eventLog.logEvent('ProtocolFocus', '%d, %3.3f' % (frameNum, self.zPoss[self.pos])) + eventLog.logEvent('ProtocolFocus', '%d, %3.3f' % (frameNum, self.zPoss[self.pos] + self.jitter_val)) TaskListProtocol.OnFrame(self, frameNum) diff --git a/PYME/recipes/acquisition.py b/PYME/recipes/acquisition.py index eca9f074b..2ae637b87 100644 --- a/PYME/recipes/acquisition.py +++ b/PYME/recipes/acquisition.py @@ -5,6 +5,9 @@ import json import numpy as np import time +import logging + +logger = logging.getLogger(__name__) def format_values(d, context): for k in d.keys(): @@ -155,3 +158,124 @@ def save(self, namespace, context={}): } }]), headers={'Content-Type': 'application/json'}) + + +@register_module('QueueZAcquisition') +class QueueZAcquisition(OutputModule): + """ + Queue move-to and start-spool acquisition tasks for each input position + using the ActionServer web wrapper queue_action endpoint. + + Parameters + ---------- + input_positions : Input + PYME.IO.tabular containing 'x_um' and 'y_um' coordinates in units of + micrometers (preferred) *or* 'x' and 'y' coordinates in units of + nanometers. + action_server_url : CStr + URL of the microscope-side action server process. + spool_settings : DictStrAny + settings to be passed to `PYME.Acquire.SpoolController.StartSpooling` as + key-word arguments. Ones that make sense in the context of this recipe + module include: + max_frames : int + number of frames to spool per series + method : str + 'File', 'Queue' (py2 only), or 'Cluster'. + hdf_compression_level: int + pytables compression level, valid for file/queue methods only. + z_stepped : bool + flag to toggle z stepping or standard protocol during spool + z_dwell : int + number of frames to acquire at each z position for z-stepped + spools + cluster_h5 : bool + spool to single h5 file on cluster (True) or pzf files (False). + Only relevant for `Cluster` method. + pzf_compression_settings : dict + see PYME.Acquire.HTTPSpooler + protocol_name : str + filename of the acquisition protocol to follow while spooling + subdirectory : str + firectory within current SpoolController set directory to spool + a given series. The directory will be created if it doesn't + already exist. + timeout : Float + time in seconds after which the acquisition tasks associated with these + positions will be ignored/unqueued from the action manager. + nice: Int + priority at which acquisition tasks should execute (default=10) + max_duration : float + A generous estimate, in seconds, of how long the task might take, after + which the lasers will be automatically turned off and the action queue + paused. + """ + input_positions = Input('input') + action_server_url = CStr('http://127.0.0.1:9393') + spool_settings = DictStrAny() + timeout = Float(np.finfo(float).max) + max_duration = Float(np.finfo(float).max) + nice = Int(10) + + def save(self, namespace, context={}): + """ + Parameters + ---------- + namespace : dict + The recipe namespace + context : dict + Information about the source file to allow pattern substitution to + generate the output name. At least 'basedir' (which is the fully + resolved directory name in which the input file resides) and + 'file_stub' (which is the filename without any extension) should be + resolved. + + Notes + ----- + str spool_settings values can context-substitute templated parameters, + e.g. spool_settings = {'subdirectory': '{file_stub}', + 'extra_metadata: {'Samples.Well': '{file_stub}'}} + """ + # substitute spool settings + spool_settings = self.spool_settings.copy() + format_values(spool_settings, context) + pos = namespace[self.input_positions] + stack_settings = dict() + for k in pos.mdh.keys(): + if k.startswith('StackSettings'): + stack_settings[k] = pos.mdh[k] + logger.debug('stack settings: %s' % str(stack_settings)) + spool_settings['stack'] = stack_settings + + if 'Sample.Well' in pos.mdh.keys(): + if 'extra_metadata' not in spool_settings.keys(): + spool_settings['extra_metadata'] = {} + spool_settings['extra_metadata']['Sample.Well'] = pos.mdh['Sample.Well'] + spool_settings['subdirectory'] = pos.mdh['Sample.Well'] + + try: # get positions in units of micrometers + positions = np.stack((pos['x_um'], + pos['y_um']), + axis=1) # (N, 2), [um] + except KeyError: # assume x and y are in nanometers + positions = np.stack((pos['x'], + pos['y']), + axis=1) / 1e3 # (N, 2), [nm] -> [um] + + dest = self.action_server_url + '/queue_actions' + session = requests.Session() + actions = list() + for ri in range(positions.shape[0]): + actions.append({ + 'CentreROIOn':{ + 'x': positions[ri, 0], 'y': positions[ri, 1], + 'then': { + 'SpoolSeries' : spool_settings + } + } + }) + session.post(dest + '?timeout=%f&nice=%d&max_duration=%f' % (self.timeout, + self.nice, + self.max_duration), + data=json.dumps(actions), + headers={'Content-Type': 'application/json'}) diff --git a/PYME/recipes/focus.py b/PYME/recipes/focus.py new file mode 100644 index 000000000..990c73891 --- /dev/null +++ b/PYME/recipes/focus.py @@ -0,0 +1,153 @@ + +from .base import register_module, ModuleBase, Filter +from .traits import Input, Output, Float +import numpy as np +from PYME.IO import tabular, MetaDataHandler +import logging + +logger = logging.getLogger(__name__) + + +class GaussFitter1D(object): + """ + 1D gaussian fitter for use with focus locks which either have line-cameras, or whose frames are summed alone one + direction to create a line profile, the peak position of which indicates the current focal position. + """ + def __init__(self, maxfev=200, min_amp=0, max_sigma=np.finfo(float).max): + """ + + Parameters + ---------- + maxfev: int + see scipy.optimize.leastsq argument by the same name + min_amp : float + minimum fit result amplitude which we are willing to accept as a + successful fit. + max_sigma : float + maximum fit result sigma which we are willing to accept as a + successful fit. + """ + self.maxfev = maxfev + self._min_amp = min_amp + self._max_sigma = max_sigma + + def _model_function(self, parameters, position): + """ + 1D gaussian + Parameters + ---------- + parameters : tuple + fit model parameters + distance : ndarray + 1D position array [pixel] + Returns + ------- + model : ndarray + """ + amplitude, center, sigma, b = parameters + return amplitude * np.exp(-((position - center) ** 2) / (2 * sigma ** 2)) + b + + def _error_function(self, parameters, position, data): + """ + """ + return data - self._model_function(parameters, position) + + def _calc_guess(self, position, data): + offset = data.min() + p95 = np.percentile(data, 95) + amplitude = p95 - offset + max_ind = np.argmin(np.abs(data - p95)) + fwhm = np.sum(data > offset + 0.5 * amplitude) + # amplitude, center, sigma, bx, b = parameters + return amplitude, position[max_ind], fwhm / 2.355, offset + + def fit(self, position, data): + from scipy import optimize + + guess = self._calc_guess(position, data) + + (res, cov_x, infodict, mesg, res_code) = optimize.leastsq(self._error_function, guess, args=(position, data), + full_output=True, maxfev=self.maxfev) + + success = res_code > 0 and res_code < 5 and res[0] > self._min_amp and res[2] < self._max_sigma + return tuple(res.astype('f')), success + + +@register_module('EstimateStackSettings') +class EstimateStackSettings(ModuleBase): + """ + + """ + input_stack = Input('input') + z_bottom = Float(49.0) + min_range = Float(3) + step_size = Float(0) + z_max = Float(58) + output = Output('with_stack_settings') + + def execute(self, namespace): + from scipy.ndimage import laplace + from PYME.Analysis.piezo_movement_correction import correct_target_positions + from PYME.recipes.processing import Threshold + + im = namespace[self.input_stack] + + z = correct_target_positions(np.arange(im.data.shape[2]), im.events, im.mdh) + if 'Multiview.ActiveViews' in im.mdh: + # dodge striping in the middle + from PYME.recipes.multiview import ExtractMultiviewChannel + lps = [] + for view in im.mdh['Multiview.ActiveViews']: + chan = ExtractMultiviewChannel(view_number=view).apply_simple(im) + lps.append(np.stack([laplace(chan.data[:,:,ind,0].squeeze()) for ind in range(chan.data.shape[2])], axis=2)) + lp = np.concatenate(lps, axis=0) + else: + lp = np.stack([laplace(im.data[:,:,ind,0].squeeze()) for ind in range(im.data.shape[2])], axis=2) + + + otsu = Threshold(method='otsu').apply_simple(im) + masked = otsu.data[:,:,:,0] * lp + metric = np.sum(masked ** 2, axis=(0, 1)) + + step_size = im.mdh['StackSettings.StepSize'] + next_stack = np.argmin(np.abs(z - z[z < z.min() + 0.5 * step_size][0])) + 1 + + fitter = GaussFitter1D() + res, success = fitter.fit(z[:next_stack], metric[:next_stack]) + if not success: + raise RuntimeError('Fit did not converge') + + if z.min() > self.z_bottom: + logger.error('stack does not go below bottom z, going to be some guesswork') + zz = np.linspace(min(z.min(), self.z_bottom), max(z.max(), self.z_max), 100) + fitted = fitter._model_function(res, zz) + # plt.plot(zz, fitted) + + target_lap = fitter._model_function(res, self.z_bottom) + top = zz[fitted > target_lap][-1] + if top - self.z_bottom < self.min_range: + logger.error('found focus smaller than minimum range, increasing') + top = self.z_bottom + self.min_range + + mdh = MetaDataHandler.DictMDHandler({ + 'StackSettings.StartPos': self.z_bottom, + 'StackSettings.EndPos': top + }) + logger.debug('StartPos %.3f, EndPos %.3f' % (self.z_bottom, top)) + if self.step_size != 0: + mdh['StackSettings.StepSize'] = self.step_size + for k in im.mdh.keys(): + if k.startswith('Sample'): + mdh[k] = im.mdh[k] + + vx, vy, _ = im.mdh.voxelsize_nm + + roi_position = tabular.DictSource({ + # set position to be center of the field of view imaged + 'x': np.asarray([1e3 * im.mdh['Positioning.x'] + vx * (im.mdh['Multiview.ROI0Origin'][0] + 0.5 * im.mdh['Multiview.ROISize'][0])]), + 'y': np.asarray([1e3 * im.mdh['Positioning.y'] + vy * (im.mdh['Multiview.ROI0Origin'][1] + 0.5 * im.mdh['Multiview.ROISize'][1])]) + }) + + roi_position.mdh = mdh + + namespace[self.output] = roi_position diff --git a/PYME/recipes/modules.py b/PYME/recipes/modules.py index 3baed6f6b..5646c1490 100644 --- a/PYME/recipes/modules.py +++ b/PYME/recipes/modules.py @@ -22,6 +22,7 @@ from . import acquisition from . import supertile from . import pointcloud +from . import focus try: from . import skfilters except: