Merge branch 'new_things' into 'v3'

Scanning and recentring

See merge request openflexure/openflexure-microscope-server!162

src/openflexure_microscope_server/things/auto_recentre_stage.py

@@ -0,0 +1,193 @@
+import numpy as np
+import logging
+import matplotlib.pyplot as plt
+import json
+import numpy as np
+import time
+
+from labthings_fastapi.thing import Thing
+from labthings_fastapi.dependencies.thing import direct_thing_client_dependency
+from labthings_fastapi.decorators import thing_action
+from labthings_sangaboard import SangaboardThing
+from labthings_picamera2.thing import StreamingPiCamera2
+from openflexure_microscope_server.things.autofocus import AutofocusThing
+from openflexure_microscope_server.things.camera_stage_mapping import CameraStageMapper
+
+StageDep = direct_thing_client_dependency(SangaboardThing, "/stage/")
+CamDep = direct_thing_client_dependency(StreamingPiCamera2, "/camera/")
+CSMDep = direct_thing_client_dependency(CameraStageMapper, "/camera_stage_mapping/")
+AutofocusDep = direct_thing_client_dependency(AutofocusThing, "/autofocus/")
+
+def _gui_description():
+    return {
+        "icon": "center_focus_weak",
+        "forms": [
+            {
+                "name": "Run recentre",
+                "route": "/recentre",
+                "submitLabel": "Run recentre",
+                "isTask": True,
+                "schema": [
+                    {
+                        "fieldType": "htmlBlock",
+                        "name": "status_display",
+                        "label": "current_status",
+                        "content": (
+                            "Find the centre of the range of motion of the microscope, based on the position of the highest point."
+                            "<br />"
+                            "<br />"
+                            "Requires a flat sample with fairly dense features to focus on."
+                        )
+                    },
+                ]
+            }
+        ]
+    }
+
+def turningpoints(lst):
+    dx = np.diff(lst)
+    return (dx[1:] * dx[:-1] < 0)
+
+def unpack_autofocus(scan_data):
+    """Extract z, sharpness data from a move_and_measure call
+    """
+    scan_data = dict(scan_data)
+    jpeg_times = scan_data['jpeg_times']
+    jpeg_sizes = scan_data['jpeg_sizes']
+    jpeg_sizes_MB = [x / 10**3 for x in jpeg_sizes]
+    stage_times = scan_data['stage_times']
+    stage_positions = scan_data['stage_positions']
+    stage_height = [pos['z'] for pos in stage_positions]
+
+    jpeg_heights = np.interp(jpeg_times,stage_times,stage_height)
+
+    turning = np.where(turningpoints(jpeg_heights))[0]+1
+
+    return jpeg_heights[turning[0]:turning[1]], jpeg_sizes_MB[turning[0]:turning[1]]
+
+
+class RecentringThing(Thing):
+    @thing_action
+    def looping_autofocus(self, autofocus: AutofocusDep, stage: StageDep, dz = 2000):
+        repeat = True
+        attempts = 0
+        while repeat and attempts < 10:
+            height_min = stage.position['z'] - dz / 2
+            height_max = stage.position['z'] + dz / 2
+            data = autofocus.fast_autofocus(dz = dz)
+            heights, _ = unpack_autofocus(data)
+            print(min(heights), max(heights), stage.position['z'])
+            time.sleep(0.3)
+            # TODO: max heights seems badly wrong! Something about turning?
+            if stage.position['z'] - height_min < dz / 5 or height_max - stage.position['z'] < dz / 5:
+                print(heights)
+                attempts += 1
+            else:
+                repeat = False
+
+    @thing_action
+    def recentre(self, autofocus: AutofocusDep, stage: StageDep, max_steps = 15, lateral_distance = 5000):
+            """ Recentre the OpenFlexure stage, based on the
+            location of the xy turning point
+            
+            Autofocuses at multiple points around the sample to 
+            find the overall maximum (or minimum) height, which
+            corresponds to the centre of the stage.
+
+            max_steps: The maximum number of moves in x or y before
+            aborting due to a poorly positioned stage or hard to focus
+            sample
+
+            lateral_distance: The xy distance between areas to check.
+            Below 3000 becomes less reliable, as focus shouldn't shift
+            much between these sites, making the procedure more sensitive
+            to noise or a failed autofocus.
+            """
+
+            max_steps = 20
+            dx = lateral_distance
+
+            centre = list(stage.position.values())
+
+            # A list of all the positions we've focused
+            focused_pos = [[],[]]
+
+            self.looping_autofocus(autofocus, stage)
+
+            for direction in [0, 1]:
+                # Start off with the current position, and moving in the positive direction
+                focused_pos[direction] = [list(stage.position.values())]
+                moves = +1
+                
+                print(centre)
+                
+                stage.move_absolute(x = centre[0], y = centre[1], z = centre[2])
+                steps = 0
+                all_heights = []
+
+                # We'll run this for x, then y
+                while True:
+                    # If we're moving in the positive direction, we want the highest point
+                    # Otherwise, we want the lowest
+                    if moves > 0:
+                        starting_point = np.max(np.array(focused_pos[direction])[:,direction])
+                    else:
+                        starting_point = np.min(np.array(focused_pos[direction])[:,direction])
+                    
+                    # Next location is an extra move in the direction we want
+                    destination = centre
+                    destination[direction] = starting_point + moves * dx
+                    stage.move_absolute(x = int(destination[0]), y = int(destination[1]), z = destination[2])
+                    self.looping_autofocus(autofocus, stage)
+                    position = list(stage.position.values())            
+                    focused_pos[direction].append(position)
+
+                    logging.info(focused_pos)
+
+                    steps += 1
+                    if steps > max_steps:
+                        logging.warning("Couldn't find a suitable position. Roughly centre the stage and check your sample is suitable for autofocus")
+                        break
+
+                    if len(focused_pos[direction]) > 4:
+                        
+                        all_heights = [x[2] for x in focused_pos[direction]]
+                        direction_index = [x[direction] for x in focused_pos[direction]]
+                        
+                        sorted_all_heights = [x for y, x in sorted(zip(direction_index, all_heights))]
+                        
+                        sorted_lateral = sorted(direction_index)
+                        quad_fit = np.polyfit(sorted_lateral, sorted_all_heights, 2)
+                        quad_fit_func = np.poly1d(quad_fit)
+
+                        turning = quad_fit_func.deriv()
+
+                        turning_loc = -turning[0]/(turning[1])
+
+                        logging.warning(sorted_all_heights)
+                        if np.argmax(sorted_all_heights) != 0 and np.argmax(sorted_all_heights) != len(all_heights) - 1:
+                            logging.info(f'Breaking because the highest point is at {np.argmax(sorted_all_heights)} in the list')
+                            # plt.plot(sorted_lateral, sorted_all_heights,'.')
+                            # plt.plot(sorted_lateral, quad_fit_func(sorted_lateral))
+                            # plt.show()
+                            break
+                        else:
+                            if turning_loc < np.min(sorted_lateral):
+                                moves = -1
+                            elif turning_loc > np.max(sorted_lateral):
+                                moves = 1
+                            else:
+                                # plt.plot(sorted_lateral, sorted_all_heights,'.')
+                                # plt.plot(sorted_lateral, quad_fit_func(sorted_lateral))
+                                # plt.show()
+                                pass
+
+                
+                # Centre value is replaced by the maximum value recorded in that axis
+                centre[direction] = focused_pos[direction][np.argmax(all_heights)][direction]
+                stage.move_absolute(x = centre[0], y = centre[1], z = centre[2])
+                self.looping_autofocus(autofocus, stage)
+
+            logging.info(f"Centre of ROM is at {centre, stage.position['z']} \n")
+
+            return focused_pos

src/openflexure_microscope_server/things/smart_scan.py

@@ -0,0 +1,315 @@
+import cv2
+import numpy as np
+import os
+import time
+from PIL import Image
+from scipy.stats import norm
+import pprint
+import logging
+from copy import deepcopy
+
+from labthings_fastapi.thing import Thing
+from labthings_fastapi.dependencies.thing import direct_thing_client_dependency
+from labthings_fastapi.decorators import thing_action
+from labthings_sangaboard import SangaboardThing
+from labthings_picamera2.thing import StreamingPiCamera2
+from openflexure_microscope_server.things.autofocus import AutofocusThing
+from openflexure_microscope_server.things.camera_stage_mapping import CameraStageMapper
+from openflexure_microscope_server.things.auto_recentre_stage import RecentringThing
+from openflexure_microscope_server.things.settings_manager import SettingsManager
+
+StageDep = direct_thing_client_dependency(SangaboardThing, "/stage/")
+CamDep = direct_thing_client_dependency(StreamingPiCamera2, "/camera/")
+CSMDep = direct_thing_client_dependency(CameraStageMapper, "/camera_stage_mapping/")
+AutofocusDep = direct_thing_client_dependency(AutofocusThing, "/autofocus/")
+RecentreStage = direct_thing_client_dependency(RecentringThing, "/auto_recentre_stage/")
+SettingsDep = direct_thing_client_dependency(SettingsManager, "/settings/")
+
+def check_dist(image, stats_list, multiple):
+    '''tests whether each pixel in image is within multiple stdevs (stats_list[1]) of the mean (stats_list[0]) for all three channels
+    returns a mask with True for pixels within that range (similar to stats list) and False otherwise
+    '''
+    return(
+        np.all(np.abs(image - stats_list[np.newaxis, np.newaxis, 0, :]) < stats_list[np.newaxis, np.newaxis, 1, :] * multiple, axis = 2)
+    )
+
+def closest(current, focused_path):
+    ''' Finds the index of the closest x-y position in a list from the current position,
+    with ties split by the later element in the list (most recently taken)
+    
+    must be float64 to deal with the huge numbers involved!'''
+    
+    current_pos = np.array(current[:2], dtype='float64')
+    path_pos = np.asarray(focused_path, dtype='float64').T[:2].T
+
+    dist_2 = np.sqrt(np.sum((path_pos - current_pos)**2, axis=1, dtype='float64'), dtype = 'float64')
+    min_dist = np.argmin(dist_2)
+    mask = np.where(dist_2 == dist_2[min_dist], 1, 0)
+    try:
+        closest = np.max(np.nonzero(mask))
+    except:
+        closest = 0
+    return closest 
+
+def unpack_autofocus(scan_data):
+    """Extract z, sharpness data from a move_and_measure call
+    
+    Data will start at `start_index`, i.e. `start_index` points are dropped
+    from the beginning of the array.
+    """
+    jpeg_times = scan_data['jpeg_times']
+    jpeg_sizes = scan_data['jpeg_sizes']
+    jpeg_sizes_MB = [x / 10**3 for x in jpeg_sizes]
+    stage_times = scan_data['stage_times']
+    stage_positions = scan_data['stage_positions']
+    stage_height = [pos[2] for pos in stage_positions]
+
+    jpeg_heights = np.interp(jpeg_times,stage_times,stage_height)
+
+    def turningpoints(lst):
+        dx = np.diff(lst)
+        return dx[1:] * dx[:-1] < 0
+
+    turning = np.where(turningpoints(jpeg_heights))[0]+1
+
+    return jpeg_heights[turning[0]:turning[1]], jpeg_sizes_MB[turning[0]:turning[1]]
+
+def limit_focus_change(prev_pos, prev_z, new_pos, new_z, limit):
+    # limit is the largest ratio of change in z to change in xy that's allowed
+
+    prev_xy = np.asarray(prev_pos, dtype = 'float64')
+    new_xy = np.asarray(new_pos, dtype = 'float64')
+
+    dist = np.sqrt(np.sum(((new_xy - prev_xy)/10**4)**2, dtype = 'float64'))
+
+    focus_change = abs(new_z - prev_z)/10**4
+    if dist == 0:
+        print(f'Not moved between {prev_pos} and {new_pos}')
+        movement_ratio = 0
+    else:
+        movement_ratio = np.divide(focus_change, dist, dtype='float64')
+
+    # print('Movement ratio is {0} in z per lateral step. The limit is {1}'.format(round(movement_ratio, 4), round(limit,4)))
+    # print(f'This is the distance between {prev_pos}, {prev_z} and {new_pos}, {new_z}')
+
+    if movement_ratio > limit:
+        return 'reject'
+    else:
+        return 'accept'
+
+def distance_to_site(current, next):
+    current = np.array(current, dtype='float64') 
+    next = np.array(next, dtype='float64')
+    if (next[1] - current[1])**2 + (next[0] - current[0])**2 < 0:
+        print(f'Negative distance between {next} and {current}')
+    return np.sqrt((next[1] - current[1])**2 + (next[0] - current[0])**2, dtype='float64')
+
+# def set_template(microscope, pos):
+#     microscope.move(pos)
+#     background = microscope.grab_image_array()
+#     background_LUV = cv2.cvtColor(background, cv2.COLOR_RGB2LUV)
+
+#     ch1 = (background_LUV.T[0]).flatten()
+#     ch2 = (background_LUV.T[1]).flatten()
+#     ch3 = (background_LUV.T[2]).flatten()
+
+#     points = np.array([np.asarray(ch1),np.asarray(ch2),np.asarray(ch3)]).T
+
+#     # we get the mean and standard deviation of values in each channel
+
+#     mu, std = np.apply_along_axis(norm.fit, 0, points)
+#     stats_list = np.vstack([mu, std])
+#     return stats_list
+
+def distance_to_site(current, next):
+    next = np.array(next, dtype='float64')
+    current = np.array(current, dtype='float64')
+    return np.sqrt( (next[1] - current[1])**2 + (next[0] - current[0])**2 )
+
+class SmartScanThing(Thing):
+    @thing_action
+    def set_template(self, settings: SettingsDep, cam: CamDep):
+        background = cam.grab_jpeg()
+        background = np.array(Image.open(background.open()))
+
+        # we're working in the LUV colourspace as it collect colours together in a human-intuitive way
+        background_LUV = cv2.cvtColor(background, cv2.COLOR_RGB2LUV)
+
+        ch1 = (background_LUV.T[0]).flatten()
+        ch2 = (background_LUV.T[1]).flatten()
+        ch3 = (background_LUV.T[2]).flatten()
+
+        points = np.array([np.asarray(ch1),np.asarray(ch2),np.asarray(ch3)]).T
+
+        # we get the mean and standard deviation of values in each channel
+
+        mu, std = np.apply_along_axis(norm.fit, 0, points)
+        stats_list = np.vstack([mu, std])
+
+        settings.update_external_metadata(data = {"background_stats": stats_list.tolist()})
+
+        current_keys = settings.external_metadata_in_state
+        logging.info(type(current_keys))
+
+        # if "1background_stats" not in current_keys:
+        current_keys.append("background_stats")
+        settings.external_metadata_in_state = deepcopy(current_keys)
+
+        logging.info(settings.external_metadata_in_state)
+        logging.info(settings.external_metadata)
+    
+    @thing_action
+    def sample_scan(self, autofocus: AutofocusDep, stage: StageDep, cam: CamDep, csm: CSMDep, settings: SettingsDep, recentre: RecentreStage):
+        # try:
+        stats_list = np.array(settings.external_metadata['background_stats'])
+        #     print(stats_list)
+        # except:
+        #     logging.warning("No template set")
+        #     raise Exception
+
+        previous_dir = 0
+
+        # if necessary, move to the starting point for your scan
+
+        starting_loc = list(stage.position.values())
+
+        recentre.looping_autofocus()
+
+        dx = 60
+        dy = 60
+
+        r = cam.grab_jpeg()
+        arr = np.array(Image.open(r.open()))
+
+        camera_stage_mapping_matrix = csm.image_to_stage_displacement_matrix
+
+        #TODO: CHECK HOW TO GET CALIBRATION IMAGE SIZE
+
+        dx = dx * arr.shape[1] / 100
+        dy = dy * arr.shape[0] / 100
+
+        dx = np.array(np.dot(np.array([0,dx]), camera_stage_mapping_matrix), dtype=int)[0]
+        dy = np.array(np.dot(np.array([dy,0]), camera_stage_mapping_matrix), dtype=int)[1]
+
+        dz = 3000
+
+        sample_coverage = 7
+
+        print(dx, dy)
+
+        # construct a 2D scan path
+        path = [[stage.position['x'], stage.position['y']]]
+
+        # a list of the sites images have been taken at, and sites with a successful autofocus
+        focused_path = []
+        true_path = []
+
+        i = 0
+
+        ids = []
+
+        start_time = time.strftime("%H_%M_%S-%d_%m_%Y")
+
+        folder_path = r'scans'
+
+        if not os.path.exists(folder_path):
+            os.makedirs(folder_path)
+        j = 0
+        while os.path.exists(os.path.join(folder_path, str(j))):
+            j += 1
+        folder_path = os.path.join(folder_path, str(j))
+        os.makedirs(folder_path)
+
+        # move to each x-y position. in z, move to the height of the closest x-y position that successfully focused
+        while len(path) > 0:
+            loc = [path[0][0], path[0][1], stage.position['z']]
+
+            path.remove(loc[:2])
+
+            stage.move_absolute(x = int(loc[0]), y = int(loc[1]), z = int(loc[2]))
+
+            if len(focused_path) > 1:
+                z_index = closest(loc, focused_path)
+                # print('Moving to {0}'.format([coords[0], coords[1], focused_path[z_index][2]]))
+                # print(focused_path)
+                stage.move_absolute(x = int(loc[0]), y = int(loc[1]), z = int(focused_path[z_index][2]))
+
+            # capture an image, convert it to LUV, then make lists of the 3 channels' values
+            img = cam.grab_jpeg()
+            img = np.array(Image.open(img.open()))
+
+            img2 = cv2.cvtColor(img, cv2.COLOR_RGB2LUV)
+            ch4 = (img2.T[0]).flatten()
+            ch5 = (img2.T[1]).flatten()
+            ch6 = (img2.T[2]).flatten()
+
+            # mask the image to only include pixels outside 5 stds of the mean of all three channels.
+            # get the percent of pixels that are in the mask (assumed sample)
+            img_mask = check_dist(img2, stats_list, 3)
+            background_coverage = round(100*np.count_nonzero(img_mask)/img_mask.size, 1)
+
+
+            # if more than 92% of the image is background, treat it as background and continue
+
+            if 100 - background_coverage < sample_coverage:
+                category = 'background'
+                # fig.set_facecolor("gray")
+            else:
+                # if not, it's sample. run an autofocus and use the updated height
+                new_pos = [[stage.position['x'] - dx, stage.position['y']], [stage.position['x'] + dx, stage.position['y']], [stage.position['x'], stage.position['y'] - dy], [stage.position['x'], stage.position['y'] + dy]]
+                for pos in new_pos:
+                    if pos not in [sublist[:2] for sublist in true_path] and pos not in path:
+                        path.append(pos)
+
+                category = 'sample'
+                # fig.set_facecolor("pink")
+
+                attempts = 0
+
+                while True:
+                    recentre.looping_autofocus(dz=3000)
+                    current_height = stage.position['z']
+
+                    # if there have been successful autofocuses in this scan, find the closest one in x-y
+                    # test if the change in z between them exceeds a ratio (indicating a failed autofocus)
+
+                    if len(focused_path) > 0:
+                        nearest_focused_site = focused_path[closest(loc, focused_path)]
+                        result = limit_focus_change(nearest_focused_site[0:2], nearest_focused_site[-1], loc[0:2], current_height, 0.3)
+
+                    # if there haven't been any previous autofocuses, we have to assume this one worked
+                    else:
+                        result = 'accept'
+
+                    # if the autofocus worked, take a new, more focused image. add the current position to the list of successful locations
+                    if result == 'accept':
+                        loc = list(stage.position.values())
+                        # img = microscope.grab_image_array()
+                        focused_path.append(loc)
+                        break
+                    else:
+                        if attempts >= 5:
+                            break
+                        else:
+                        # if the autofocus was rejected, we return to the height of the closest successful autofocus. not perfect, but better than wandering out of focus
+                            print('Resetting from a previous focus position. Options are {0}, we chose {1}'.format(focused_path, nearest_focused_site))
+                            stage.move_absolute(z = int(focused_path[z_index][2]))
+                            attempts += 1
+
+                img = cam.capture_jpeg()
+                img = np.array(Image.open(img.open()))
+                img = cv2.resize(img, (0,0), fx = 0.5, fy = 0.5)
+                img = Image.fromarray(img)
+                img.save(os.path.join(folder_path, f"rabbit_{loc[0]}_{loc[1]}.jpg"))
+
+            img_preview = cam.grab_jpeg()
+            img_preview = np.array(Image.open(img_preview.open()))
+            # add the current position to the list of all positions visited
+            true_path.append(loc)
+
+            path = sorted(path, key = lambda x: distance_to_site(loc[:2], x))
+
+            if len(true_path) > 750:
+                break
+        
+        stage.move_absolute(x = starting_loc[0], y = starting_loc[1], z = starting_loc[2])