diff --git a/.gitignore b/.gitignore
index 563e733a..2e927d4e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -33,6 +33,7 @@ nosetests.xml
.mr.developer.cfg
.project
.pydevproject
+.idea
*.egg-info/
*.iml
@@ -45,4 +46,4 @@ build/
dist/
frozen_version.py
plugins/
-src/
\ No newline at end of file
+src/
diff --git a/README.md b/README.md
index d5c83ff6..5e4c4850 100644
--- a/README.md
+++ b/README.md
@@ -13,6 +13,8 @@ Please see help here: https://github.com/CellProfiler/CellProfiler/blob/master/c
cd PLUGIN_DIRECTORY
git clone https://github.com/CellProfiler/CellProfiler-plugins.git
```
+
+ Alternatively download zip and manually extract to PLUGIN_DIRECTORY.
1. Install required dependencies:
```
diff --git a/cellstar/__init__.py b/cellstar/__init__.py
new file mode 100644
index 00000000..e06d1359
--- /dev/null
+++ b/cellstar/__init__.py
@@ -0,0 +1,8 @@
+# -*- coding: utf-8 -*-
+"""
+CellStar package providing CellStar algorithm for segmentation of yeast cells in brightfield imagery.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+__author__ = 'Adam Kaczmarek, Filip Mróz, Szymon Stoma'
+__all__ = ["segmentation"]
diff --git a/cellstar/core/__init__.py b/cellstar/core/__init__.py
new file mode 100644
index 00000000..c4540384
--- /dev/null
+++ b/cellstar/core/__init__.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""
+Core package including main components used in CellStar segmentation.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+__all__ = ["image_repo", "point", "seed", "seeder", "snake", "snake_filter"]
diff --git a/cellstar/core/config.py b/cellstar/core/config.py
new file mode 100644
index 00000000..73e89bbb
--- /dev/null
+++ b/cellstar/core/config.py
@@ -0,0 +1,84 @@
+# -*- coding: utf-8 -*-
+"""
+Config is a storage for default CellStar configuration.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+
+def default_config():
+ return {
+ 'segmentation': {
+ 'foreground': {
+ 'FillHolesWithAreaSmallerThan': 2.26,
+ 'MaskDilation': 0.136,
+ 'MaskMinRadius': 0.34,
+ 'MaskThreshold': 0.03,
+ 'pickyDetection': False,
+ 'blur': 1,
+ 'MinCellClusterArea': 0.85
+ },
+ 'avgCellDiameter': 35,
+ 'background': {
+ 'blurSteps': 50,
+ 'computeByBlurring': 0.5,
+ 'blur': 0.3
+ },
+ 'ranking': {
+ 'avgInnerBrightnessWeight': 10,
+ 'avgBorderBrightnessWeight': 300,
+ 'stickingWeight': 60,
+ 'shift': 0.68,
+ 'maxInnerBrightnessWeight': 10,
+ 'logAreaBonus': 18,
+ 'maxRank': 100,
+ 'avgInnerDarknessWeight': 0
+ },
+ 'minArea': 0.07,
+ 'cellBorder': {
+ 'medianFilter': 0.1
+ },
+ 'maxFreeBorder': 0.4,
+ 'cellContent': {
+ 'MaskThreshold': 0.0,
+ 'medianFilter': 0.17,
+ 'blur': 0.6
+ },
+ 'steps': 2,
+ 'maxArea': 2.83,
+ 'stars': {
+ 'cumBrightnessWeight': 304.45,
+ 'maxSize': 1.67,
+ 'gradientWeight': 15.482,
+ 'sizeWeight': [189.4082],
+ 'brightnessWeight': 0.0442,
+ 'step': 0.0335,
+ 'points': 28,
+ 'borderThickness': 0.1,
+ 'unstick': 0.3,
+ 'backgroundWeight': 0.0,
+ 'smoothness': 7.0,
+ 'gradientBlur': 0.0
+ },
+ 'minAvgInnerDarkness': 0.1,
+ 'maxOverlap': 0.3,
+ 'seeding': {
+ 'from': {
+ 'cellContentRandom': 0,
+ 'cellBorderRemovingCurrSegmentsRandom': 0,
+ 'cellContentRemovingCurrSegments': 1,
+ 'snakesCentroids': 0,
+ 'cellContent': 0,
+ 'cellContentRemovingCurrSegmentsRandom': 0,
+ 'cellBorderRemovingCurrSegments': 0,
+ 'cellBorder': 1,
+ 'snakesCentroidsRandom': 0,
+ 'cellBorderRandom': 0
+ },
+ 'ContentBlur': 2,
+ 'randomDiskRadius': 0.33,
+ 'minDistance': 0.27,
+ 'BorderBlur': 2
+ }
+ }
+ }
diff --git a/cellstar/core/image_repo.py b/cellstar/core/image_repo.py
new file mode 100644
index 00000000..fe74435e
--- /dev/null
+++ b/cellstar/core/image_repo.py
@@ -0,0 +1,318 @@
+# -*- coding: utf-8 -*-
+"""
+Images repository calculates and stores intermediate images used in segmentation process.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import math
+
+from cellstar.utils.image_util import *
+
+
+class ImageRepo(object):
+ """
+ After initialization with input image it calculates on demand intermediate images such as foreground,
+ background mask, and others.
+ """
+
+ @property
+ def background(self):
+ if self._background is None:
+ self.calculate_background()
+
+ return self._background
+
+ @background.setter
+ def background(self, new_background):
+ self._background = new_background
+
+ @property
+ def brighter(self):
+ if self._brighter is None:
+ self.calculate_brighter()
+
+ return self._brighter
+
+ @property
+ def brighter_original(self):
+ if self._brighter is None:
+ self.calculate_brighter_original()
+
+ return self._brighter_original
+
+ @property
+ def darker(self):
+ if self._darker is None:
+ self.calculate_darker()
+
+ return self._darker
+
+ @property
+ def darker_original(self):
+ if self._darker_original is None:
+ self.calculate_darker_original()
+
+ return self._darker_original
+
+ @property
+ def image_back_difference(self):
+ if self._clean_original is None:
+ self.calculate_clean_original()
+
+ return self._clean_original
+
+ @property
+ def image_back_difference_blurred(self):
+ if self._clean is None:
+ self.calculate_clean()
+
+ return self._clean
+
+ @property
+ def foreground_mask(self):
+ if self._foreground_mask is None:
+ self.calculate_forebackground_masks()
+
+ return self._foreground_mask
+
+ @property
+ def background_mask(self):
+ if self._background_mask is None:
+ self.calculate_forebackground_masks()
+
+ return self._background_mask
+
+ @property
+ def cell_content_mask(self):
+ if self._cell_content_mask is None:
+ self.calculate_cell_border_content_mask()
+
+ return self._cell_content_mask
+
+ @property
+ def cell_border_mask(self):
+ if self._cell_border_mask is None:
+ self.calculate_cell_border_content_mask()
+
+ return self._cell_border_mask
+
+ @property
+ def segmentation(self):
+ if self._segmentation is None:
+ self.init_segmentation()
+
+ return self._segmentation
+
+ @property
+ def mask(self):
+ if self._mask is None:
+ self._mask = np.ones(self._image_original.shape, dtype=bool)
+
+ return self._mask
+
+ def apply_mask(self, ignore_mask):
+ """
+ @param ignore_mask: binary mask of places which are to be ignored
+ @type ignore_mask: np.ndarray
+ """
+ use_mask = np.logical_not(ignore_mask)
+ masked_image = self._image_original * use_mask
+ filler_value = np.median(masked_image)
+ self.image = masked_image + filler_value * ignore_mask
+ self._mask = use_mask
+
+ def __init__(self, image, parameters):
+ """
+ @param image: 0-1 float array
+ @type image: np.ndarray
+ @type parameters: dict
+ """
+ self.parameters = parameters
+
+ # float arrays
+ self.image = image
+ self._image_original = image
+ self._background = None
+ self._mask = None
+ self._brighter_original = None
+ self._darker_original = None
+ self._clean_original = None
+ self._brighter = None
+ self._darker = None
+ self._clean = None
+
+ # binary masks
+ self._foreground_mask = None
+ self._background_mask = None
+ self._cell_border_mask = None
+ self._cell_content_mask = None
+
+ # cache
+ self._blurred = [] # (image,blur,image_blurred)
+
+ # segmentation labels
+ self._segmentation = None
+
+ def init_segmentation(self):
+ self._segmentation = np.zeros(self.image.shape[:2], int)
+
+ def calculate_background(self, background_mask=None):
+ """
+ Fills in background pixels based on foreground pixels using blur.
+ @param background_mask: pixels that belong to the background, if not given background is calculated from image
+ """
+
+ # No background mask is provided so calculate one based on edges in the image.
+ if background_mask is None:
+ smoothed = image_smooth(self.image, int(self.parameters["segmentation"]["background"]["computeByBlurring"]
+ * self.parameters["segmentation"]["avgCellDiameter"]))
+ foreground_mask = image_normalize(abs(self.image - smoothed)) \
+ > self.parameters["segmentation"]["foreground"]["MaskThreshold"]
+
+ foreground_mask = \
+ fill_foreground_holes(foreground_mask,
+ self.parameters["segmentation"]["foreground"]["MaskDilation"]
+ * self.parameters["segmentation"]["avgCellDiameter"],
+ self.parameters["segmentation"]["foreground"]["FillHolesWithAreaSmallerThan"]
+ * self.parameters["segmentation"]["avgCellDiameter"] ** 2 * math.pi / 4,
+ self.parameters["segmentation"]["foreground"]["MinCellClusterArea"]
+ * self.parameters["segmentation"]["avgCellDiameter"] ** 2 * math.pi / 4,
+ self.parameters["segmentation"]["foreground"]["MaskMinRadius"]
+ * self.parameters["segmentation"]["avgCellDiameter"]
+ )
+ background_mask = np.logical_not(foreground_mask)
+
+ if self._mask is not None:
+ background_mask |= np.logical_not(self._mask)
+
+ filler_value = np.median(self.image)
+ background = self.image.astype(float)
+ foreground_mask = np.logical_not(background_mask)
+
+ if background_mask.any():
+ filler_value = np.median(background[background_mask])
+
+ background = background * background_mask + filler_value * foreground_mask
+
+ # Spread foreground to background pixels
+ smooth_radius = round(self.parameters["segmentation"]["background"]["blur"]
+ * self.parameters["segmentation"]["avgCellDiameter"])
+ steps = self.parameters["segmentation"]["background"]["blurSteps"]
+
+ for i in xrange(steps):
+ background = image_smooth(background, 1 + round(smooth_radius * ((steps - i) / steps) ** 2), i != 0)
+ background = background * foreground_mask + self.image * background_mask
+
+ self._background = background
+ self._foreground_mask = foreground_mask
+
+ def calculate_clean_original(self):
+ self._clean_original = image_normalize(self.image - self.background)
+
+ def calculate_brighter_original(self):
+ self._brighter_original = self.image - self.background
+ self._brighter_original = image_normalize(np.maximum(self._brighter_original, 0))
+
+ def calculate_darker_original(self):
+ self._darker_original = self.background - self.image
+ self._darker_original = image_normalize(np.maximum(self._darker_original, 0))
+
+ def calculate_forebackground_masks(self):
+ """
+ Determines foreground which is the part of the image different substantialy from provided background.
+ Holes in initial foreground are filled using morphological operations
+ """
+
+ # Calculate foreground based on blurred brighter/darker.
+ darker_blurred = self.get_blurred(self.darker_original, self.parameters["segmentation"]["foreground"]["blur"])
+ brighter_blurred = self.get_blurred(self.brighter_original,
+ self.parameters["segmentation"]["foreground"]["blur"])
+ self._foreground_mask = \
+ (darker_blurred + brighter_blurred) > self.parameters["segmentation"]["foreground"]["MaskThreshold"]
+
+ if self.parameters["segmentation"]["foreground"]["pickyDetection"]:
+ smooth_coefficient = round(self.parameters["segmentation"]["background"]["computeByBlurring"]
+ * self.parameters["segmentation"]["avgCellDiameter"])
+
+ temp_blurred = image_smooth(self.image, smooth_coefficient)
+ temp_fg_mask = image_normalize(np.abs((self.image - temp_blurred))) > \
+ self.parameters["segmentation"]["foreground"]["MaskThreshold"]
+
+ self._foreground_mask = np.logical_and(self.foreground_mask, temp_fg_mask)
+
+ self._foreground_mask = \
+ fill_foreground_holes(self.foreground_mask,
+ self.parameters["segmentation"]["foreground"]["MaskDilation"]
+ * self.parameters["segmentation"]["avgCellDiameter"],
+ self.parameters["segmentation"]["foreground"]["FillHolesWithAreaSmallerThan"]
+ * self.parameters["segmentation"]["avgCellDiameter"] ** 2 * math.pi / 4,
+ self.parameters["segmentation"]["foreground"]["MinCellClusterArea"]
+ * self.parameters["segmentation"]["avgCellDiameter"] ** 2 * math.pi / 4,
+ self.parameters["segmentation"]["foreground"]["MaskMinRadius"]
+ * self.parameters["segmentation"]["avgCellDiameter"]
+ )
+
+ if self._mask is not None:
+ self._foreground_mask &= self._mask
+
+ self._background_mask = np.logical_not(self.foreground_mask)
+
+ def calculate_clean(self):
+ image_diff_med_filter_size = round(self.parameters["segmentation"]["stars"]["gradientBlur"]
+ * self.parameters["segmentation"]["avgCellDiameter"])
+ self._clean = image_smooth(self.image_back_difference, image_diff_med_filter_size)
+
+ clean_mean = self._clean_original.mean()
+ masked = self._clean * self.foreground_mask
+ mean_negative_masked = clean_mean * self.background_mask
+ self._clean = masked + mean_negative_masked
+
+ def calculate_brighter(self):
+ brighter_med_filter_size = np.round(self.parameters["segmentation"]["cellBorder"]["medianFilter"]
+ * self.parameters["segmentation"]["avgCellDiameter"])
+ self._brighter = image_median_filter(self.brighter_original, brighter_med_filter_size)
+ self._brighter *= self.foreground_mask
+
+ def calculate_darker(self):
+ darker_med_filter_size = round(self.parameters["segmentation"]["cellContent"]["medianFilter"]
+ * self.parameters["segmentation"]["avgCellDiameter"])
+ self._darker = image_median_filter(self.darker_original, darker_med_filter_size)
+ self._darker *= self.foreground_mask
+
+ def calculate_cell_border_content_mask(self):
+ """
+ Splits foreground pixels into cell content and cell borders.
+ """
+ blur_iterations = int(math.ceil(self.parameters["segmentation"]["cellContent"]["blur"]
+ * self.parameters["segmentation"]["avgCellDiameter"] - 0.5))
+ darker_blurred = self.get_blurred(self.darker, blur_iterations)
+
+ if self.parameters["segmentation"]["cellContent"]["MaskThreshold"] == 0:
+ eroded_foreground = image_erode(self.foreground_mask,
+ self.parameters["segmentation"]["foreground"]["MaskDilation"]
+ * self.parameters["segmentation"]["avgCellDiameter"])
+
+ if not eroded_foreground.any():
+ cell_content_mask_threshold = 0.0
+ else:
+ cell_content_mask_threshold = np.median(self._darker[eroded_foreground])
+ if np.max(self._darker[eroded_foreground]) == cell_content_mask_threshold: # there is only one value
+ cell_content_mask_threshold = 0.0 # take all
+ else:
+ cell_content_mask_threshold = self.parameters["segmentation"]["cellContent"]["MaskThreshold"]
+
+ self._cell_content_mask = (self.brighter == 0) & self.foreground_mask & \
+ (darker_blurred > cell_content_mask_threshold)
+
+ self._cell_border_mask = self.foreground_mask & np.logical_not(self.cell_content_mask)
+
+ def get_blurred(self, image, blur_param, cache_result=False):
+ cache = [c for a, b, c in self._blurred if a is image and b == blur_param]
+ if cache:
+ return cache[0]
+ else:
+ blurred = image_blur(image, blur_param)
+ if cache_result:
+ self._blurred.append((image, blur_param, blurred))
+ return blurred
diff --git a/cellstar/core/parallel/__init__.py b/cellstar/core/parallel/__init__.py
new file mode 100644
index 00000000..f97ef5af
--- /dev/null
+++ b/cellstar/core/parallel/__init__.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""
+Parallel is an attempt to use multiprocessing in snake growing thus improving the speed of the segmentation.
+Currently not tested and on hold.
+Date: 2015-2016
+Website: http://cellstar-algorithm.org/
+"""
diff --git a/cellstar/core/parallel/snake_grow.py b/cellstar/core/parallel/snake_grow.py
new file mode 100644
index 00000000..aa123a53
--- /dev/null
+++ b/cellstar/core/parallel/snake_grow.py
@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+"""
+Snake grow is an attempt to run time expensive snake in parallel using multiprocessing package.
+Currently not tested and on hold. Tried only once for parameter fitting.
+Date: 2015-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import ctypes
+from copy import copy
+from multiprocessing import Pool, Array, Manager
+
+import numpy as np
+
+from cellstar.core.image_repo import ImageRepo
+from cellstar.core.polar_transform import PolarTransform
+from cellstar.core.snake import Snake
+from cellstar.parameter_fitting.pf_snake import PFSnake
+
+
+def conv_single_image(image):
+ shared_array_base = Array(ctypes.c_double, image.size)
+ shared_array = np.ctypeslib.as_array(shared_array_base.get_obj())
+ shared_array = shared_array.reshape(image.shape)
+ shared_array[:] = image
+
+ return shared_array
+
+
+def conv_image_repo(images):
+ return map(conv_single_image, [
+ images.foreground_mask,
+ images.brighter,
+ images.darker,
+ images.image_back_difference_blurred,
+ images.image_back_difference,
+ images.cell_content_mask,
+ images.cell_border_mask,
+ images.background,
+ images.background_mask
+ ])
+
+
+def grow_single_snake(frame, images, parameters, seed):
+ #
+ #
+ # RECONSTRUCT INPUT
+ #
+ #
+
+ ir = ImageRepo(frame, parameters)
+ ir._foreground_mask, ir._brighter, ir._darker, ir._clean, ir._clean_original, \
+ ir._cell_content_mask, ir._cell_border_mask, ir._background, ir._background_mask = images
+
+ #
+ #
+ # CREATE AND GROW SNAKE
+ #
+ #
+
+ polar_transform = PolarTransform.instance(parameters["segmentation"]["avgCellDiameter"],
+ parameters["segmentation"]["stars"]["points"],
+ parameters["segmentation"]["stars"]["step"],
+ parameters["segmentation"]["stars"]["maxSize"])
+
+ s = Snake.create_from_seed(parameters, seed, parameters["segmentation"]["stars"]["points"], ir)
+
+ size_weight_list = parameters["segmentation"]["stars"]["sizeWeight"]
+ snakes_to_grow = [(copy(s), w) for w in size_weight_list]
+
+ for snake, weight in snakes_to_grow:
+ snake.grow(size_weight=weight, polar_transform=polar_transform)
+ snake.evaluate(polar_transform)
+
+ best_snake = sorted(snakes_to_grow, key=lambda (sn, _): sn.rank)[0][0]
+
+ pf_s = PFSnake(None, None, None, best_snake=best_snake)
+ pf_s.best_snake = best_snake
+
+ return pf_s
+
+
+def grow_fun((seed, frame, images, parameters)):
+ return grow_single_snake(frame, images, parameters, seed)
+
+
+def add_snake(snakes, snake):
+ snakes.append(snake)
+
+
+def mp_snake_grow(images, parameters, seeds):
+ snakes = []
+ manager = Manager()
+ shared_parameters = manager.dict(parameters)
+ shared_images = conv_image_repo(images)
+ shared_frame = conv_single_image(images.image)
+
+ snakes = []
+
+ pool = Pool(processes=8)
+ snakes = pool.map_async(grow_fun, [(seed, shared_frame, shared_images, shared_parameters) for seed in seeds]).get()
+ pool.close()
+ pool.join()
+
+ return snakes
+
diff --git a/cellstar/core/point.py b/cellstar/core/point.py
new file mode 100644
index 00000000..2142819f
--- /dev/null
+++ b/cellstar/core/point.py
@@ -0,0 +1,49 @@
+# -*- coding: utf-8 -*-
+"""
+Integer point wrapper.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import math
+
+
+class Point(object):
+ """
+ @ivar x: x coordinate of point
+ @ivar y: y coordinate of point
+ """
+
+ def __init__(self, x, y):
+ """
+ @type x: int
+ @type y: int
+ """
+ self.x = x
+ self.y = y
+
+ def __repr__(self):
+ return "Seed(x={0},y={1})".format(self.x, self.y)
+
+ def __eq__(self, other):
+ return self.x == other.x and self.y == other.y
+
+ def polar_coords(self, origin):
+ """
+ @type origin : Point
+ @return: radius, angle
+ @rtype: (int,int)
+ """
+ r_x = self.x - origin.x
+ r_y = self.y - origin.y
+
+ radius = math.sqrt(r_x ** 2 + r_y ** 2)
+ angle = math.atan2(r_y, r_x)
+
+ return radius, angle
+
+ def as_xy(self):
+ return self.x, self.y
+
+ def euclidean_distance_to(self, other_point):
+ return math.sqrt((self.x - other_point.x) ** 2 + (self.y - other_point.y) ** 2)
diff --git a/cellstar/core/polar_transform.py b/cellstar/core/polar_transform.py
new file mode 100644
index 00000000..73a77f0f
--- /dev/null
+++ b/cellstar/core/polar_transform.py
@@ -0,0 +1,171 @@
+# -*- coding: utf-8 -*-
+"""
+Polar transform provides tools for fast and convenient operation in polar domain.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import math
+import threading
+
+import numpy as np
+
+from cellstar.utils.calc_util import to_int
+from cellstar.utils.calc_util import sub2ind
+from cellstar.utils.image_util import image_dilate_with_element, get_circle_kernel
+
+
+class PolarTransform(object):
+ """
+ Object wrapping polar transform calculations and cached properties
+
+ @type N: int
+ @ivar N: number of points for polar transform calculation
+
+ @type distance: float
+ @ivar distance: maximal distance from cell center to cell border
+
+ @type step: float
+ @ivar step: length of step for active contour along its axis in pixels
+
+ @type steps: int
+ @ivar steps: number of steps considered for active contour along single axis
+
+ @type max_r: int
+ @ivar max_r: maximal radius of active contour
+
+ @type R: numpy.array
+ @ivar R: consecutive radii values for single axis of active contour
+
+ @type center: int
+ @ivar center: Polar transform center coordinate
+
+ @type edge: int
+ @ivar edge: dimension of polar transform
+
+ @type halfedge: int
+ @ivar halfedge: half of polar transform dimension
+
+ @type t: numpy.array
+ @ivar t: angles (in radians) of consecutive rays casted from 'center'
+
+ @type x: numpy.array
+ @ivar x: cartesian x-coordinates of points in polar coordinates system
+ coordinates ordered by radius of polar points --> x[r,a] = P(r,a).x
+
+ @type y: numpy.array
+ @ivar y: cartesian y-coordinates of points in polar coordinates system
+ coordinates ordered by radius of polar points --> y[r,a] = P(r,a).y
+
+ @type dot_voronoi: numpy.array
+ @ivar dot_voronoi - voronoi - "gravity field" of contour points
+ dot_voronoi[x,y] = id(closest_contour_point(x,y))
+
+ @type to_polar: dict
+ @ivar to_polar - dictionary of lists
+ for each point:
+ to_polar[index(P(R,a)] - list of point id in voronoi of contour points {P(r,a)| 0 < r < R}
+ to_polar[index(P(R,a)] = [gravity_field(dot_voronoi, p) for p in {P(r,a) | 0 < r < R}]
+ to_polar[index(P(R,a)] =
+ [index(x,y) for x,y in range((0,0),(edge,edge)) if dot_voronoi[x,y] == gravity_index(p) for p in {P(r,a) | 0 < r < R}]
+ """
+
+ __singleton_lock = threading.Lock()
+ __singleton_instances = {}
+
+ @classmethod
+ def instance(cls, avg_cell_diameter, points, step, max_size):
+ init_params = avg_cell_diameter, points, step, max_size
+ if not cls.__singleton_instances.get(init_params, False):
+ with cls.__singleton_lock:
+ if not cls.__singleton_instances.get(init_params, False):
+ cls.__singleton_instances[init_params] = cls(avg_cell_diameter, points, step, max_size)
+ return cls.__singleton_instances.get(init_params, None)
+
+ def __init__(self, avg_cell_diameter, points_number, step, max_size):
+ self.N = points_number
+ self.distance = max_size * avg_cell_diameter
+ self.step = max(step * avg_cell_diameter, 0.2)
+ self.steps = 1 + int(round((self.distance + 2) / self.step))
+ self.max_r = min(1 + int(round(self.distance / self.step)), self.steps - 1)
+
+ self.R = None
+ self.center = None
+ self.edge = None
+ self.half_edge = None
+ self.x = None
+ self.y = None
+
+ self.dot_voronoi = None
+ self.to_polar = {}
+
+ self._calculate_polar_transform()
+
+ def _calculate_polar_transform(self):
+ self.R = np.arange(1, self.steps + 1).reshape((1, self.steps)).transpose() * self.step
+
+ # rays angle from cell center
+ self.t = np.linspace(0, 2 * math.pi, self.N + 1)
+ self.t = self.t[:-1]
+
+ # sinus and cosinus of rays angle repeated steps-times
+ # function value for angles and every radius (for a given angle same for every radius)
+ sin_t = np.kron(np.ones((len(self.R), 1)), np.sin(self.t))
+ cos_t = np.kron(np.ones((len(self.R), 1)), np.cos(self.t))
+
+ # N-times repeated vector of subsequent radiuses
+ RR = np.kron(np.ones((1, len(self.t))), self.R)
+
+ # From polar definition:
+ # x - matrix of xs for angle alpha and radius R
+ # y - matrix of ys for angle alpha and radius R
+ self.x = RR * cos_t
+ self.y = RR * sin_t
+
+ self.half_edge = math.ceil(self.R[-1] + 2)
+ self.center = to_int(self.half_edge + 1)
+ self.edge = to_int(self.center + self.half_edge)
+
+ # clear black image [edge x edge]
+ self.dot_voronoi = np.zeros((self.edge, self.edge), dtype=int)
+ px = self.center + self.x
+ py = self.center + self.y
+
+ # create list of coordinates (x,y) on the checked contour
+ index = np.column_stack(((py - .5).astype(int).T.flat, (px - .5).astype(int).T.flat))
+
+ # create list of subsequent id for above points
+ cont = np.arange(1, px.size + 1)
+
+ # mark on 'dot_voronoi' every point using unique id
+ self.dot_voronoi[tuple(index.T)] = cont
+
+ # in every iteration smooth 'dot_voronoi' marking gravity field of given points
+ for i in range(0, int(self.center)):
+ ndv = image_dilate_with_element(self.dot_voronoi, 3)
+ mask = np.logical_and((self.dot_voronoi == 0), (ndv != 0))
+ mask = mask.nonzero()
+ self.dot_voronoi[mask] = ndv[mask]
+
+ # apply circle mask on 'dot_voronoi'
+ circ_mask = get_circle_kernel(self.half_edge)
+ self.dot_voronoi[np.logical_not(circ_mask)] = 0
+ self.dot_voronoi[self.center - 1, self.center - 1] = 0
+
+ # for every angle
+ for a in range(self.t.size):
+ # create new clear mask
+ mask = np.zeros((self.edge, self.edge), dtype=bool)
+ # for point
+ for r in range(self.R.size):
+ # find index of point P(r,a)
+ idx = sub2ind(px.shape[0], (r, a))
+ val = idx + 1
+ # find places which belong to that index in 'dot_voronoi'
+ indices = np.array(zip(*np.nonzero(self.dot_voronoi == val)))
+ # set mask to 1 in above places
+ if len(indices) > 0:
+ mask[tuple(indices.T)] = 1
+
+ # to_polar[idx] is a list of coordinates (x,y) from points on the mask
+ self.to_polar[idx] = map(lambda pair: pair, zip(*np.nonzero(mask)))
diff --git a/cellstar/core/seed.py b/cellstar/core/seed.py
new file mode 100644
index 00000000..d2f4a299
--- /dev/null
+++ b/cellstar/core/seed.py
@@ -0,0 +1,25 @@
+# -*- coding: utf-8 -*-
+"""
+Contour seed.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+from cellstar.core.point import Point
+
+
+class Seed(Point):
+ """
+ Object of a seed.
+ @ivar x: x coordinate of seed
+ @ivar y: y coordinate of seed
+ @ivar origin: where seed comes from ('content' or 'background' or 'snakes')
+ """
+
+ def __init__(self, x, y, origin):
+ """
+ @type x: int
+ @type y: int
+ """
+ super(Seed, self).__init__(x, y)
+ self.origin = origin
diff --git a/cellstar/core/seeder.py b/cellstar/core/seeder.py
new file mode 100644
index 00000000..af27a85a
--- /dev/null
+++ b/cellstar/core/seeder.py
@@ -0,0 +1,320 @@
+# -*- coding: utf-8 -*-
+"""
+Seeder is responsible for finding places where might be cells.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+import random
+
+from cellstar.core.seed import Seed
+from cellstar.utils.calc_util import *
+from cellstar.utils.image_util import *
+
+
+class Seeder(object):
+ """
+ Find places for seeds based on images and previously found snakes.
+ """
+
+ def __init__(self, images, parameters):
+ """
+ @type images: core.image_repo.ImageRepo
+ @type parameters: dict
+ """
+ self.images = images
+ np.random.seed(abs(hash(abs(np.sum(images.image)))))
+ random.seed(abs(np.sum(images.image)))
+ self.parameters = parameters
+ self.cluster_min_distance = self.parameters["segmentation"]["seeding"]["minDistance"] \
+ * self.parameters["segmentation"]["avgCellDiameter"]
+ self.random_radius = self.parameters["segmentation"]["seeding"]["randomDiskRadius"] \
+ * self.parameters["segmentation"]["avgCellDiameter"]
+
+ def cluster_seeds(self, seeds):
+ """
+ @type seeds: list[Seed]
+ """
+ origin = 'cluster'
+ if len(seeds) > 0:
+ origin = seeds[0].origin
+ seeds = np.array(map(lambda s: (s.x, s.y), seeds))
+ my_inf = \
+ np.inf if seeds.dtype.kind == 'f' else np.iinfo(seeds.dtype).max
+ while len(seeds) > 1:
+ #
+ # find p1 and p2 closest to each other - p1 is closest to p2
+ # and p2 is closest to p1
+ #
+ p1 = np.arange(seeds.shape[0])
+ d = seeds[:, np.newaxis, :] - seeds[np.newaxis, :, :]
+ d2 = np.sum(d * d, 2)
+ d2[p1, p1] = my_inf
+ p2 = np.argmin(d2, 0)
+ #
+ # Eliminate p1 / p2 if p1 is closest to p2 and p2 is closest to
+ # someone else
+ #
+ good = p1 == p2[p2]
+ p1, p2 = p1[good], p2[good]
+ #
+ # Eliminate p1 / p2 if p1 > p2 (get rid of (2, 1) since there is
+ # a (1, 2)
+ #
+ good = p1 < p2
+ p1, p2 = p1[good], p2[good]
+ #
+ # Eliminate p1 / p2 if < self.cluster_min_distance
+ #
+ good = d2[p1, p2] < self.cluster_min_distance
+ p1, p2 = p1[good], p2[good]
+ if len(p1) == 0:
+ break
+ #
+ # coalesce
+ #
+ new_seeds = (seeds[p1, :] + seeds[p2, :]) / 2.0
+ to_keep = np.ones(seeds.shape[0], bool)
+ to_keep[p1] = False
+ to_keep[p2] = False
+ seeds = np.vstack((seeds[to_keep, :], new_seeds))
+
+ seed_points = point_list_as_seeds(seeds, origin)
+ return seed_points
+
+ @staticmethod
+ def rand_seeds(max_random_radius, times, seeds, min_random_radius=0):
+ """
+ Generate (len(seeds) * times) random seeds within (min_random_radius:max_random_radius).
+ @type seeds: list[Seed]
+ """
+
+ seeds_number = len(seeds)
+ new_seeds_number = int(seeds_number * times)
+
+ shuffled_seeds = list(seeds)
+ np.random.shuffle(shuffled_seeds)
+ px = multiply_list([s.x for s in shuffled_seeds], times)
+ py = multiply_list([s.y for s in shuffled_seeds], times)
+
+ random_angle = np.random.random(new_seeds_number) * 2 * math.pi
+ random_radius = np.random.random(new_seeds_number) * (max_random_radius - min_random_radius) + min_random_radius
+
+ rpx = px + random_radius * np.cos(random_angle)
+ rpy = py + random_radius * np.sin(random_angle)
+
+ return [Seed(x, y, 'rand') for x, y in zip(rpx, rpy)]
+
+ def find_seeds_from_border_or_content(self, image, foreground_mask, segments, mode):
+ """
+ Finds seeds from given image.
+ @param image: image (border or content) from which seeds are being extracted
+ @type image: np.ndarray
+ @param foreground_mask: binary foreground mask
+ @type foreground_mask: np.ndarray
+ @param segment: array with segments marked
+ @type segments: np.ndarray
+ @param mode: 'border' or 'content' determining if look for maxima or minima
+ """
+
+ im_name = ''
+ origin = ''
+
+ if mode == 'border':
+ blur = self.parameters["segmentation"]["seeding"]["BorderBlur"] \
+ * self.parameters["segmentation"]["avgCellDiameter"]
+ im_name = 'border' + im_name
+ image = set_image_border(image, 1)
+ excl_value = 1
+ else:
+ blur = self.parameters["segmentation"]["seeding"]["ContentBlur"] \
+ * self.parameters["segmentation"]["avgCellDiameter"]
+ im_name = 'content' + im_name
+ excl_value = 0
+
+ excluding = segments is not None
+ if excluding:
+ image = exclude_segments(image, segments, excl_value)
+ origin = ' no segments'
+ im_name += ' excluding current segments'
+
+ origin = im_name + origin
+
+ blurred = self.images.get_blurred(image, blur, cache_result=not excluding)
+
+ if mode == 'border':
+ blurred = 1 - blurred
+
+ maxima = find_maxima(blurred) * foreground_mask
+
+ maxima_coords = sp.nonzero(maxima)
+ seeds = zip(maxima_coords[1], maxima_coords[0])
+
+ seed_points = point_list_as_seeds(seeds, origin)
+ seed_points = self.cluster_seeds(seed_points)
+
+ return seed_points
+
+ def find_seeds_from_snakes(self, snakes):
+ """
+ Finds seeds from snakes centroids
+ @param snakes: Grown snakes from previous frame
+ @type snakes: list[Snake]
+ """
+ return point_list_as_seeds([snake.centroid for snake in snakes], 'snake centroid')
+
+ def find_seeds(self, snakes, all_seeds, exclude_current_segments):
+ seeds = []
+
+ # First iteration of segmentation.
+ if not exclude_current_segments:
+ if self.parameters["segmentation"]["seeding"]["from"]["cellBorder"]:
+ new_seeds = self.find_seeds_from_border_or_content(self.images.brighter, self.images.foreground_mask,
+ None, 'border')
+ new_seeds += self.rand_seeds(
+ self.random_radius,
+ self.parameters["segmentation"]["seeding"]["from"]["cellBorderRandom"],
+ new_seeds
+ )
+
+ seeds += new_seeds
+
+ if self.parameters["segmentation"]["seeding"]["from"]["cellContent"]:
+ new_seeds = self.find_seeds_from_border_or_content(self.images.darker, self.images.foreground_mask,
+ None, 'content')
+ new_seeds += self.rand_seeds(
+ self.random_radius,
+ self.parameters["segmentation"]["seeding"]["from"]["cellContentRandom"],
+ new_seeds
+ )
+
+ seeds += new_seeds
+
+ elif len(snakes) > 0:
+ if self.parameters["segmentation"]["seeding"]["from"]["cellBorderRemovingCurrSegments"]:
+ new_seeds = self.find_seeds_from_border_or_content(self.images.brighter, self.images.foreground_mask,
+ self.images.segmentation, 'border')
+ new_seeds += self.rand_seeds(
+ self.random_radius,
+ self.parameters["segmentation"]["seeding"]["from"]["cellBorderRemovingCurrSegmentsRandom"],
+ new_seeds
+ )
+
+ seeds += new_seeds
+
+ if self.parameters["segmentation"]["seeding"]["from"]["cellContentRemovingCurrSegments"]:
+ new_seeds = self.find_seeds_from_border_or_content(self.images.darker, self.images.foreground_mask,
+ self.images.segmentation, 'content')
+ new_seeds += self.rand_seeds(
+ self.random_radius,
+ self.parameters["segmentation"]["seeding"]["from"]["cellContentRemovingCurrSegmentsRandom"],
+ new_seeds
+ )
+
+ seeds += new_seeds
+
+ # Seeds from existing contours
+ if self.parameters["segmentation"]["seeding"]["from"]["snakesCentroids"]:
+ new_seeds = self.find_seeds_from_snakes(snakes)
+ new_seeds += self.rand_seeds(
+ self.random_radius,
+ self.parameters["segmentation"]["seeding"]["from"]["snakesCentroidsRandom"],
+ new_seeds
+ )
+ seeds += new_seeds
+
+ seeds = self.filter_seeds(self.images.image.shape, seeds, all_seeds)
+
+ return seeds
+
+ def filter_seeds(self, image_shape, seeds, all_seeds):
+ """
+ Ignore seeds that are close image borders or already processed seeds.
+ @param image_shape: (y,x) image size
+ @type seeds: [cellstar.core.seed.Seed]
+ @type all_seeds: [cellstar.core.seed.Seed]
+ @return:
+ """
+ distance = self.parameters["segmentation"]["stars"]["step"] * self.parameters["segmentation"]["avgCellDiameter"]
+ distance = float(max(distance, 0.5)) # not less than half of pixel length
+ ok_seeds = np.array([False for _ in seeds])
+
+ grid_size = int(round(max(image_shape) * 1.1 / distance))
+ im_y, im_x = image_shape
+
+ # Create grid
+ seeds_grid = SeedGrid(grid_size)
+
+ # Fill grid with previous seeds
+ for seed in all_seeds:
+ x = int(round(seed.x / distance))
+ y = int(round(seed.y / distance))
+ for xx in [x - 1, x, x + 1]:
+ for yy in [y - 1, y, y + 1]:
+ seeds_grid.add_seed(xx, yy, seed)
+
+ # Fill grid with current seeds
+ for i in xrange(len(seeds)):
+ seed = seeds[i]
+ x = max(0, int(round(seed.x / distance)))
+ y = max(0, int(round(seed.y / distance)))
+
+ # Validate seed if it lies inside grid
+ if seeds_grid.inside(x, y):
+ ok_seeds[i] = seed_is_new(seed, seeds_grid.get(x, y), distance)
+
+ for xx in [x - 1, x, x + 1]:
+ for yy in [y - 1, y, y + 1]:
+ seeds_grid.add_seed(xx, yy, seed)
+
+ # Remove seeds in image borders
+ seeds_x_ok = np.array([im_x - 0.5 > seed.x > 0.5 for seed in seeds])
+ seeds_y_ok = np.array([im_y - 0.5 > seed.y > 0.5 for seed in seeds])
+
+ # Filtered seeds boolean vector
+ ok_seeds = np.logical_and(ok_seeds, np.logical_and(seeds_y_ok, seeds_x_ok))
+
+ return [seeds[i] for i in range(len(seeds)) if ok_seeds[i]]
+
+
+def seed_is_new(seed, current_seeds, distance):
+ return all([euclidean_norm(seed.as_xy(), cs.as_xy()) >= distance for cs in current_seeds])
+
+
+def point_list_as_seeds(pointsYX, origin):
+ return [Seed(point[0], point[1], origin) for point in pointsYX]
+
+
+class SeedGrid(object):
+ def __init__(self, size=10):
+ self._size_x = size
+ self._size_y = size
+ self._grid = []
+ self._init_grid()
+
+ def _init_grid(self):
+ self._grid = np.empty((self._size_x, self._size_y), dtype=np.object)
+ self._grid.fill([])
+
+ def _resize(self, to_x, to_y):
+ self._size_y = max(self._size_y, to_y)
+ self._size_x = max(self._size_x, to_x)
+
+ self._grid = np.resize(self._grid, (self._size_x, self._size_y))
+
+ def get(self, x, y):
+ return self._grid[x][y]
+
+ def size(self):
+ return self._size_x, self._size_y
+
+ def add_seed(self, x, y, seed):
+ if x < 0 or y < 0:
+ return
+
+ if x >= self._size_x or y >= self._size_y:
+ self._resize(x + 1, y + 1)
+
+ self._grid[x][y] = self._grid[x][y] + [seed]
+
+ def inside(self, x, y):
+ return x < self._size_x and y < self._size_y
diff --git a/cellstar/core/snake.py b/cellstar/core/snake.py
new file mode 100644
index 00000000..92601324
--- /dev/null
+++ b/cellstar/core/snake.py
@@ -0,0 +1,380 @@
+# -*- coding: utf-8 -*-
+"""
+Snake object responsible for growing best contours from given seed.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import math
+
+import numpy as np
+
+from cellstar.core.point import Point
+from cellstar.utils import calc_util, image_util
+from cellstar.utils.debug_util import *
+from cellstar.utils.index import Index
+
+
+class Snake(object):
+ """
+ Contour representation of snake along with all its properties.
+ @ivar in_polygon : cropped array representation of contour
+ @ivar in_polygon_yx : position of in_polygon
+ @ivar rank: ranking of the contour (the smaller the better)
+ """
+ epsilon = 1e-10
+ max_rank = 10000
+
+ @property
+ def xs(self):
+ return [p.x for p in self.points]
+
+ @property
+ def ys(self):
+ return [p.y for p in self.points]
+
+ @property
+ def centroid(self):
+ return self.centroid_x, self.centroid_y
+
+ @property
+ def in_polygon_slice(self):
+ if self.in_polygon_yx is None or self.in_polygon is None:
+ return None
+
+ yx = self.in_polygon_yx
+ local = self.in_polygon
+ return slice(yx[0], yx[0] + local.shape[0]), slice(yx[1], yx[1] + local.shape[1])
+
+ @classmethod
+ def create_from_seed(cls, parameters, seed, points_num, images):
+ points = [seed] * points_num
+ return cls(parameters, seed, points, images)
+
+ def __init__(self, parameters, seed, points=None, images=None):
+ """
+ @type seed: Seed
+ @type images: core.image_repo.ImageRepo
+ @type points: list[Point]
+ @type parameters: dict
+ """
+ self.seed = seed
+ self.parameters = parameters
+ self.images = images
+
+ # Snake grow.
+ self.points = points
+ self.original_edgepoints = None
+ self.final_edgepoints = None
+ self.polar_coordinate_boundary = None
+
+ # Snake evaluation properties.
+ self.in_polygon = None
+ self.in_polygon_yx = None
+ self.rank = None
+ self.area = 0.0
+ self.avg_out_border_brightness = 0.0
+ self.max_out_border_brightness = 0.0
+ self.avg_in_border_brightness = 0.0
+ self.avg_inner_brightness = 0.0
+ self.max_inner_brightness = 0.0
+ self.avg_inner_darkness = 0.0
+ self.centroid_x = 0.0
+ self.centroid_y = 0.0
+ self.max_contiguous_free_border = 0
+ self.free_border_entropy = 0.0
+ self.properties_vector_cached = {}
+
+ def grow(self, size_weight, polar_transform):
+ """
+ Grow the snake from seed.
+ @type polar_transform: cellstar.core.vectorized.polar_transform.PolarTransform
+ @type size_weight: float
+ """
+
+ #
+ # Initialization
+ #
+
+ self.centroid_x = self.seed.x
+ self.centroid_y = self.seed.y
+
+ points_number = polar_transform.N
+ step = polar_transform.step
+ unstick = self.parameters["segmentation"]["stars"]["unstick"]
+ avg_cell_diameter = self.parameters["segmentation"]["avgCellDiameter"]
+ smoothness = self.parameters["segmentation"]["stars"]["smoothness"]
+ gradient_weight = self.parameters["segmentation"]["stars"]["gradientWeight"]
+ brightness_weight = self.parameters["segmentation"]["stars"]["brightnessWeight"]
+ size_weight = float(size_weight) / avg_cell_diameter
+ cum_brightness_weight = self.parameters["segmentation"]["stars"]["cumBrightnessWeight"] / avg_cell_diameter
+ background_weight = self.parameters["segmentation"]["stars"]["backgroundWeight"] / avg_cell_diameter
+ border_thickness_steps = \
+ 1 + math.floor(float(self.parameters["segmentation"]["stars"]["borderThickness"]) / float(step))
+
+ im = self.images.image_back_difference_blurred
+ imb = self.images.brighter
+ imfg = self.images.foreground_mask
+
+ steps = polar_transform.steps
+ max_r = polar_transform.max_r
+ R = polar_transform.R.flat
+ t = polar_transform.t
+ max_diff = (abs(smoothness) * np.arange(1, steps + 1) / points_number + 0.5).astype(int)
+
+ px = float(self.centroid_x) + polar_transform.x
+ px = np.maximum(px, 0)
+ px = np.minimum(px, im.shape[1] - 1)
+
+ py = float(self.centroid_y) + polar_transform.y
+ py = np.maximum(py, 0)
+ py = np.minimum(py, im.shape[0] - 1)
+
+ #
+ # Contour quality function calculation and finding the best candidates.
+ #
+ #
+ # index is ordered first by angle then by radius
+ # for angle:
+ # for radius:
+
+ index = Index.create(px.round(), py.round()).reshape(
+ (polar_transform.x.shape[0], polar_transform.x.shape[1], 2))
+
+ #
+ # pre_f - part of function quality for interior
+ # f_tot - final quality array
+ #
+
+ numpy_index = Index.to_numpy(index)
+ pre_f = (cum_brightness_weight * imb[numpy_index] \
+ + background_weight * (1 - imfg[numpy_index])) * step
+ f = np.cumsum(pre_f, axis=0)
+
+ im_diff = calc_util.get_gradient(im, index, border_thickness_steps)
+
+ f_tot = f \
+ - size_weight * np.kron(np.log(R), np.ones((t.size, 1))).T \
+ - gradient_weight * im_diff \
+ - brightness_weight * im[numpy_index]
+
+ f_tot = f_tot.T
+ # f_tot = image_util.image_smooth(f_tot, 1)
+
+ # Scale entire array to 0-1 then scale individual angles.
+ f_tot = (f_tot - f_tot.min()) / (f_tot.max() - f_tot.min() + self.epsilon)
+ f_tot /= (np.kron(np.ones((f_tot.shape[1], 1)), f_tot.max(axis=1)).T + self.epsilon)
+
+ # Find initial contour
+ _, best_radius = np.where(f_tot == np.kron(np.ones((f_tot.shape[1], 1)), f_tot.min(axis=1)).T)
+
+ #
+ # Contour smoothing
+ #
+
+ smoothed_radius, radius_bounds = self.smooth_contour(best_radius, max_diff, points_number, f_tot)
+
+ self.original_edgepoints = (smoothed_radius != radius_bounds) | \
+ ((smoothed_radius == best_radius) & (smoothed_radius < max_r - 2))
+ smoothed_radius = np.minimum(smoothed_radius, max_r)
+
+ # Determine final edgepoint that will be used to construct contour points
+ self.final_edgepoints = \
+ calc_util.unstick_contour(self.original_edgepoints, unstick)
+
+ # Interpolate points where no reliable points had been found.
+ calc_util.interpolate_radiuses(self.final_edgepoints, points_number, smoothed_radius)
+
+ #
+ # Create contour points list
+ #
+
+ final_radius = np.minimum(np.maximum(smoothed_radius, 1), max_r - 1)
+ #final_radius = np.minimum(np.maximum(np.round(smoothed_radius + 1), 1), max_r - 1)
+
+ px = self.seed.x + step * final_radius * np.cos(t.T)
+ py = self.seed.y + step * final_radius * np.sin(t.T)
+
+ self.polar_coordinate_boundary = final_radius
+ self.points = [Point(x, y) for x, y in zip(px, py)]
+
+ def smooth_contour(self, radius, max_diff, points_number, f_tot):
+ """
+ Smoothing contour using greedy length cut. Rotating from min radius clockwise and anti.
+ @type radius: np.ndarray
+ @param max_diff: max change of ray length per iter.
+ @type max_diff np.ndarray
+ @type points_number int
+ @param f_tot: quality function array
+ @type f_tot: np.ndarray
+
+ @rtype (np.ndarray, np.ndarray)
+ @return (smoothed_radius, used_radius_bounds)
+ """
+ min_angle = radius.argmin()
+ istart = min_angle
+
+ xmins2 = np.copy(radius)
+ xmaxs = np.copy(radius)
+
+ changed = True
+
+ def cut_rotate(start, step):
+ last_ok = False
+ iteration = 0
+ any_cut = False
+ while not (iteration >= points_number and last_ok):
+ current = (start + iteration * step) % points_number
+ previous = (current - 1 * step) % points_number
+
+ if xmins2[current] - xmins2[previous] > max_diff[xmins2[previous]]:
+ xmaxs[current] = xmins2[previous] + max_diff[xmins2[previous]]
+ f_tot_slice = f_tot[current, :xmaxs[current] + 1]
+ xmins2[current] = f_tot_slice.argmin()
+ if step < 0:
+ xmins2[current] = max(xmins2[current], xmins2[previous] - max_diff[xmins2[previous]])
+ last_ok = False
+ any_cut = True
+ else:
+ last_ok = True
+
+ iteration += 1
+ return (start + iteration * step) % points_number, any_cut
+
+ current_position = istart
+ while changed:
+ changed = False
+
+ current_position, any_cut = cut_rotate(current_position, 1)
+ changed = changed or any_cut
+
+ current_position, any_cut = cut_rotate(current_position, -1)
+ changed = changed or any_cut
+
+ return xmins2, xmaxs
+
+ def evaluate(self, polar_transform):
+ """
+ Analyse contour and calculate all it properties and ranking.
+ @type polar_transform: cellstar.core.vectorized.polar_transform.PolarTransform
+ """
+
+ # Potentially prevent unnecessary calculations
+ if self.rank is not None:
+ return
+
+ self.properties_vector_cached = {}
+
+ avg_cell_diameter = self.parameters["segmentation"]["avgCellDiameter"]
+ min_border_r = 0.055 * avg_cell_diameter
+ max_border_r = 0.1 * avg_cell_diameter
+ ranking_params = self.parameters["segmentation"]["ranking"]
+
+ image_bounds = calc_util.get_cartesian_bounds(self.polar_coordinate_boundary, self.seed.x, self.seed.y,
+ polar_transform)
+ dilated_bounds = image_util.extend_slices(image_bounds, int(max_border_r * 2))
+
+ original_clean = self.images.image_back_difference[dilated_bounds]
+ brighter = self.images.brighter[dilated_bounds]
+ cell_content_mask = self.images.cell_content_mask[dilated_bounds]
+ origin_in_slice = calc_util.inslice_point((self.seed.y, self.seed.x), dilated_bounds)
+ origin_in_slice = Point(x=origin_in_slice[1], y=origin_in_slice[0])
+
+ segment, self.in_polygon, in_polygon_local_yx = calc_util.star_in_polygon(brighter.shape,
+ self.polar_coordinate_boundary,
+ origin_in_slice.x, origin_in_slice.y,
+ polar_transform)
+ self.in_polygon_yx = calc_util.unslice_point(in_polygon_local_yx, dilated_bounds)
+
+ self.area = np.count_nonzero(self.in_polygon) + self.epsilon
+ approx_radius = math.sqrt(self.area / math.pi)
+ border_radius = max(min(approx_radius, max_border_r), min_border_r)
+
+ dilation = round(border_radius / polar_transform.step)
+
+ dilated_boundary = np.minimum(self.polar_coordinate_boundary + dilation, len(polar_transform.R) - 1)
+ eroded_boundary = np.maximum(self.polar_coordinate_boundary - dilation, 1)
+
+ dilated, _, _ = calc_util.star_in_polygon(brighter.shape, dilated_boundary,
+ origin_in_slice.x, origin_in_slice.y, polar_transform)
+
+ eroded, _, _ = calc_util.star_in_polygon(brighter.shape, eroded_boundary,
+ origin_in_slice.x, origin_in_slice.y, polar_transform)
+
+ out_border = dilated ^ segment
+ in_border = segment ^ eroded
+
+ out_border_area = np.count_nonzero(out_border) + self.epsilon
+ in_border_area = np.count_nonzero(in_border) + self.epsilon
+
+ # Calculate outer border brightness
+ tmp = original_clean[out_border]
+ self.avg_out_border_brightness = tmp.sum() / out_border_area
+ self.max_out_border_brightness = tmp.max() if tmp != [] else 0
+
+ # Calculate inner border brightness
+ tmp = original_clean[in_border]
+ self.avg_in_border_brightness = tmp.sum() / in_border_area
+
+ # Calculate inner brightness
+ tmp = brighter[segment]
+ self.avg_inner_brightness = tmp.sum() / self.area
+ self.max_inner_brightness = tmp.max() if tmp != [] else 0
+
+ # Calculate inner darkness
+ tmp = cell_content_mask[segment]
+ self.avg_inner_darkness = tmp.sum() / self.area
+
+ # Calculate snake centroid
+ if self.area > 2 * self.epsilon:
+ area2 = float(self.area - self.epsilon)
+ self.centroid_x = (self.in_polygon.sum(0) * np.arange(1, self.in_polygon.shape[1] + 1)).sum() / area2
+ self.centroid_y = (self.in_polygon.sum(1) * np.arange(1, self.in_polygon.shape[0] + 1)).sum() / area2
+ self.centroid_x += self.in_polygon_yx[1]
+ self.centroid_y += self.in_polygon_yx[0]
+ else:
+ self.centroid_x = self.xs[0]
+ self.centroid_y = self.ys[0]
+
+ # Calculate free border fragments - fragments of border, which endpoints have been discarded
+ fb, _, _ = calc_util.loop_connected_components(np.logical_not(self.final_edgepoints))
+
+ # Calculate free border entropy
+ fb_entropy = 0
+ if min(fb.shape) != 0:
+ fb_entropy = float(np.sum(fb * fb.T)) / len(self.xs) ** 2
+
+ # The longest continuous fragment of free border
+ self.max_contiguous_free_border = fb.max() if fb.size > 0 else 0
+
+ self.free_border_entropy = fb_entropy
+
+ if all(self.polar_coordinate_boundary <= 1): # minimal step is 1 (see snake_grow)
+ self.rank = Snake.max_rank
+ else:
+ self.rank = self.star_rank(ranking_params, avg_cell_diameter)
+
+ def star_rank(self, ranking_params, avg_cell_diameter):
+ return np.dot(self.ranking_parameters_vector(ranking_params), self.properties_vector(avg_cell_diameter))
+
+ def ranking_parameters_vector(self, ranking_params):
+ return np.array([
+ float(ranking_params["maxInnerBrightnessWeight"]),
+ float(ranking_params["avgInnerBrightnessWeight"]),
+ float(ranking_params["avgBorderBrightnessWeight"]),
+ float(ranking_params["avgInnerDarknessWeight"]),
+ float(ranking_params["logAreaBonus"]),
+ float(ranking_params["stickingWeight"]),
+ ])
+
+ def properties_vector(self, avg_cell_diameter):
+ if avg_cell_diameter not in self.properties_vector_cached:
+ self.properties_vector_cached[avg_cell_diameter] = np.array([
+ self.max_inner_brightness,
+ self.avg_inner_brightness,
+ self.avg_in_border_brightness - self.avg_out_border_brightness,
+ -self.avg_inner_darkness,
+ -math.log(self.area ** (1.0 / avg_cell_diameter)),
+ self.free_border_entropy
+ ])
+ return self.properties_vector_cached[avg_cell_diameter]
diff --git a/cellstar/core/snake_filter.py b/cellstar/core/snake_filter.py
new file mode 100644
index 00000000..b12267de
--- /dev/null
+++ b/cellstar/core/snake_filter.py
@@ -0,0 +1,118 @@
+# -*- coding: utf-8 -*-
+"""
+SnakeFilter is responsible for ranking and filtering out contours that as incorrect or overlap better ones.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import logging
+import math
+
+import numpy as np
+
+from cellstar.core.snake import Snake
+
+
+class SnakeFilter(object):
+ """
+ Order snakes based on their ranking and checks them for violating other constraints.
+ Discard snakes that overlap with already approved ones.
+ """
+
+ def __init__(self, images, parameters):
+ """
+ @type parameters: dict
+ @type images: core.image_repo.ImageRepo
+ """
+ self.parameters = parameters
+ self.images = images
+
+ def is_single_snake_discarded(self, snake):
+ """
+ @type snake: Snake
+ @rtype: bool
+ """
+ if snake.rank >= Snake.max_rank:
+ return True
+
+ if snake.avg_inner_darkness < self.parameters["segmentation"]["minAvgInnerDarkness"]:
+ return True
+
+ max_area = self.parameters["segmentation"]["maxArea"] * self.parameters["segmentation"]["avgCellDiameter"]**2 * math.pi / 4
+ if snake.area > max_area:
+ return True
+
+ min_area = self.parameters["segmentation"]["minArea"] * self.parameters["segmentation"]["avgCellDiameter"]**2 * math.pi / 4
+ if snake.area < min_area:
+ return True
+
+ max_free_border = self.parameters["segmentation"]["stars"]["points"] * self.parameters["segmentation"]["maxFreeBorder"]
+ if snake.max_contiguous_free_border > max_free_border:
+ return True
+
+ return False
+
+ def filter(self, snakes):
+ """
+ @type snakes: list[Snake]
+ @rtype: list[Snake]
+ """
+ logging.basicConfig(format='%(asctime)-15s %(message)s', level=logging.DEBUG)
+ logger = logging.getLogger(__package__)
+ log_message = "Discarding snake {0} for {1}: {2}"
+
+ original = self.images.image
+ filtered_snakes = []
+ segments = np.zeros(original.shape, dtype=int)
+
+ # do not allow cells on masked areas
+ segments[self.images.mask == 0] = -1
+
+ if len(snakes) > 0:
+ snakes_sorted = sorted(enumerate(snakes), key=lambda x: x[1].rank)
+ current_accepted_snake_index = 1
+ for i in xrange(len(snakes_sorted)):
+ curr_snake = snakes_sorted[i][1]
+ snake_index = snakes_sorted[i][0]
+
+ if curr_snake.rank >= Snake.max_rank:
+ logger.debug(log_message.format(snake_index, 'too high rank', curr_snake.rank))
+ break
+
+ local_snake = curr_snake.in_polygon
+ sxy = curr_snake.in_polygon_slice
+ local_segments = segments[sxy]
+
+ overlap_area = np.count_nonzero(np.logical_and(local_segments, local_snake))
+ overlap = float(overlap_area) / curr_snake.area
+
+ if overlap > self.parameters["segmentation"]["maxOverlap"]:
+ logger.debug(log_message.format(snake_index, 'too much overlapping', overlap))
+ else:
+ vacant_snake = np.logical_and(local_snake, local_segments == 0)
+ vacant_cell_content = vacant_snake[self.images.cell_content_mask[sxy]]
+ curr_snake.area = np.count_nonzero(vacant_snake) + Snake.epsilon
+ avg_inner_darkness = float(np.count_nonzero(vacant_cell_content)) / float(curr_snake.area)
+ if avg_inner_darkness < self.parameters["segmentation"]["minAvgInnerDarkness"]:
+ logger.debug(log_message.format(snake_index, 'too low inner darkness', '...'))
+ else:
+ if curr_snake.area > (self.parameters["segmentation"]["maxArea"] * self.parameters["segmentation"]["avgCellDiameter"]**2 * math.pi / 4):
+ logger.debug(log_message.format(snake_index, 'too big area', str(curr_snake.area)))
+ else:
+ if curr_snake.area < (self.parameters["segmentation"]["minArea"] * self.parameters["segmentation"]["avgCellDiameter"]**2 * math.pi / 4):
+ logger.debug(log_message.format(snake_index, 'too small area:', str(curr_snake.area)))
+ else:
+ max_free_border = self.parameters["segmentation"]["stars"]["points"] * self.parameters["segmentation"]["maxFreeBorder"]
+ if curr_snake.max_contiguous_free_border > max_free_border:
+ logger.debug(log_message.format(snake_index,
+ 'too long contiguous free border',
+ str(curr_snake.max_contiguous_free_border) +
+ ' over ' + str(max_free_border)))
+ else:
+ local_segments[[vacant_snake]] = current_accepted_snake_index
+ filtered_snakes.append(curr_snake)
+ current_accepted_snake_index += 1
+
+ segments *= self.images.mask # clear mask
+ self.images._segmentation = segments
+ return filtered_snakes
diff --git a/cellstar/parameter_fitting/__init__.py b/cellstar/parameter_fitting/__init__.py
new file mode 100644
index 00000000..ec8deddf
--- /dev/null
+++ b/cellstar/parameter_fitting/__init__.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""
+Parameter fitting package includes all components for automated parameter fitting based on a provided ground truth/
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+__all__ = ["pf_process", "pf_snake"]
\ No newline at end of file
diff --git a/cellstar/parameter_fitting/pf_auto_params.py b/cellstar/parameter_fitting/pf_auto_params.py
new file mode 100644
index 00000000..56a2e337
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_auto_params.py
@@ -0,0 +1,148 @@
+# -*- coding: utf-8 -*-
+"""
+Pf auto params contains parameters that are optimised as well as encode / decode functions.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import numpy as np
+
+parameters_range = {#"borderThickness": (0.001, 1.0),
+ "brightnessWeight": (-0.4, 0.4),
+ "cumBrightnessWeight": (0, 500),
+ "gradientWeight": (-30, 30),
+ "sizeWeight": (10, 300),
+ #"smoothness": (3, 8)
+}
+
+rank_parameters_range = {"avgBorderBrightnessWeight": (0, 600),
+ "avgInnerBrightnessWeight": (-100, 100),
+ "avgInnerDarknessWeight": (-100, 100),
+ "logAreaBonus": (5, 50),
+ "maxInnerBrightnessWeight": (-10, 50),
+ # "maxRank": (5, 300),
+ # "stickingWeight": (0, 120) # this is set to 60 rest of parameters should adapt to it
+}
+
+
+class OptimisationBounds(object):
+ def __init__(self, size=None, xmax=1, xmin=0):
+ self.xmax = xmax
+ self.xmin = xmin
+ if size is not None:
+ self.xmax = [xmax] * size
+ self.xmin = [xmin] * size
+
+ @staticmethod
+ def from_ranges(ranges_dict):
+ bounds = OptimisationBounds()
+ bounds.xmin = []
+ bounds.xmax = []
+ # bound only two parameters
+ for k, v in list(sorted(ranges_dict.iteritems())):
+ if k == "borderThickness":
+ bounds.xmin.append(0.001)
+ bounds.xmax.append(2)
+ elif k == "smoothness":
+ bounds.xmin.append(4.0)
+ bounds.xmax.append(10.0)
+ else:
+ bounds.xmin.append(-1000000)
+ bounds.xmax.append(1000000)
+
+ # bounds.xmin, bounds.xmax = zip(*zip(*list(sorted(ranges_dict.iteritems())))[1])
+ return bounds
+
+ def __call__(self, **kwargs):
+ x = kwargs["x_new"]
+ tmax = bool(np.all(x <= self.xmax))
+ tmin = bool(np.all(x >= self.xmin))
+ return tmax and tmin
+
+
+ContourBounds = OptimisationBounds.from_ranges(parameters_range)
+RankBounds = OptimisationBounds(size=len(rank_parameters_range))
+
+
+#
+#
+# PARAMETERS ENCODE DECODE
+#
+#
+
+def pf_parameters_encode(parameters):
+ """
+ brightnessWeight: 0.0442 +brightness on cell edges
+ cumBrightnessWeight: 304.45 -brightness in the cell center
+ gradientWeight: 15.482 +gradent on the cell edges
+ sizeWeight: 189.4082 (if list -> avg. will be comp.) +big cells
+ @param parameters: dictionary segmentation.stars
+ """
+ parameters = parameters["segmentation"]["stars"]
+ point = []
+ for name, (_, _) in sorted(parameters_range.iteritems()):
+ val = parameters[name]
+ if name == "sizeWeight":
+ if not isinstance(val, float):
+ val = np.mean(val)
+ point.append(val) # no scaling
+ return point
+
+
+def pf_parameters_decode(param_vector, org_size_weights_list):
+ """
+ sizeWeight is one number (mean of the future list)
+ @type param_vector: numpy.ndarray
+ @return:
+ """
+ parameters = {}
+ for (name, (_, _)), val in zip(sorted(parameters_range.iteritems()), param_vector):
+ if name == "sizeWeight":
+ val = list(np.array(org_size_weights_list) * (val / np.mean(org_size_weights_list)))
+ elif name == "borderThickness":
+ val = min(max(0.001, val), 3)
+ parameters[name] = val
+
+ # set from default
+ parameters["borderThickness"] = 0.1
+ parameters["smoothness"] = 6
+ return parameters
+
+
+def pf_rank_parameters_encode(parameters, complete_params_given=True):
+ """
+ avgBorderBrightnessWeight: 300
+ avgInnerBrightnessWeight: 10
+ avgInnerDarknessWeight: 0
+ logAreaBonus: 18
+ maxInnerBrightnessWeight: 10
+ @param parameters: dictionary all ranking params or a complete
+ @param complete_params_given: is parameters a complete dictionary
+ """
+ if complete_params_given:
+ parameters = parameters["segmentation"]["ranking"]
+
+ point = []
+ for name, (vmin, vmax) in sorted(rank_parameters_range.iteritems()):
+ val = parameters[name]
+ if vmax - vmin == 0:
+ point.append(0)
+ else:
+ point.append((val - vmin) / float(vmax - vmin)) # scaling to [0,1]
+ return point
+
+
+def pf_rank_parameters_decode(param_vector):
+ """
+ @type param_vector: numpy.ndarray
+ @return: only ranking parameters as a dict
+ """
+ parameters = {}
+ for (name, (vmin, vmax)), val in zip(sorted(rank_parameters_range.iteritems()), param_vector):
+ rescaled = vmin + val * (vmax - vmin)
+ parameters[name] = rescaled
+
+ # set from default
+ parameters["stickingWeight"] = 60
+
+ return parameters
diff --git a/cellstar/parameter_fitting/pf_mutator.py b/cellstar/parameter_fitting/pf_mutator.py
new file mode 100644
index 00000000..4a9a3563
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_mutator.py
@@ -0,0 +1,80 @@
+# -*- coding: utf-8 -*-
+"""
+Mutator can be used to change (ie mutate) existing snakes to provide higher variability in ground_truth pool.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import copy
+import random
+
+import numpy as np
+
+from cellstar.core.point import *
+from cellstar.utils.calc_util import polar_to_cartesian
+
+
+def add_mutations(gt_and_grown, avg_cell_diameter):
+ mutants = []
+ mutation_radiuses = 0.2 * avg_cell_diameter
+ for (gt, grown) in gt_and_grown:
+ mutants += [
+ (gt, grown.create_mutation(mutation_radiuses * 2, random_poly=True)),
+ (gt, grown.create_mutation(-mutation_radiuses * 2, random_poly=True)),
+ (gt, grown.create_mutation(mutation_radiuses)), (gt, grown.create_mutation(-mutation_radiuses)),
+ ]
+ return gt_and_grown + mutants
+
+
+def create_mutant_from_change(org_snake, polar_transform, boundary_change):
+ mutant_snake = copy.copy(org_snake)
+ # zero rank so it recalculates
+ mutant_snake.rank = None
+
+ # constrain change
+ new_boundary = mutant_snake.polar_coordinate_boundary + boundary_change
+ while (new_boundary <= 3).all() and abs(boundary_change).max() > 3:
+ new_boundary = np.maximum(np.minimum(
+ mutant_snake.polar_coordinate_boundary + boundary_change,
+ len(polar_transform.R) - 1), 3)
+ boundary_change /= 1.3
+
+ px, py = polar_to_cartesian(new_boundary, mutant_snake.seed.x, mutant_snake.seed.y, polar_transform)
+
+ mutant_snake.polar_coordinate_boundary = new_boundary
+ mutant_snake.points = [Point(x, y) for x, y in zip(px, py)]
+
+ # need to update self.final_edgepoints to calculate properties (for now we ignore this property)
+ mutant_snake.evaluate(polar_transform)
+
+ return mutant_snake
+
+
+def create_poly_mutation(org_snake, polar_transform, max_diff):
+ # change to pixels
+ length = org_snake.polar_coordinate_boundary.size
+ max_diff /= polar_transform.step
+
+ def polynomial(x1, x2):
+ def eval(x):
+ return x * (x - length) * (x - x1) * (x * 0.4 - x2)
+
+ return eval
+
+ poly = polynomial(random.uniform(0.001, length), random.uniform(0.001, length))
+ boundary_change = np.array([poly(x) for x in range(length)])
+
+ M = abs(boundary_change).max()
+ boundary_change = boundary_change / M * max_diff
+
+ mutant_snake = create_mutant_from_change(org_snake, polar_transform, boundary_change)
+ return mutant_snake
+
+
+def create_mutation(org_snake, polar_transform, dilation):
+ # change to pixels
+ dilation /= polar_transform.step
+ boundary_change = np.array([dilation for _ in range(org_snake.polar_coordinate_boundary.size)])
+
+ mutant_snake = create_mutant_from_change(org_snake, polar_transform, boundary_change)
+ return mutant_snake
diff --git a/cellstar/parameter_fitting/pf_process.py b/cellstar/parameter_fitting/pf_process.py
new file mode 100644
index 00000000..a1a71fdd
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_process.py
@@ -0,0 +1,364 @@
+# -*- coding: utf-8 -*-
+"""
+Pf process is the core of contour parameters fitting.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import copy
+import random
+import sys
+import time
+from multiprocessing import Process, Queue
+
+import numpy as np
+import scipy.optimize as opt
+from scipy.linalg import norm
+
+random.seed(1)
+np.random.seed(1)
+
+import logging
+
+logger = logging.getLogger(__name__)
+
+from cellstar.utils.params_util import *
+from cellstar.core.seed import Seed
+from cellstar.core.image_repo import ImageRepo
+from cellstar.parameter_fitting.pf_snake import PFSnake, GTSnake
+from cellstar.core.seeder import Seeder
+from cellstar.utils.debug_util import explore_cellstar
+from cellstar.parameter_fitting.pf_auto_params import pf_parameters_encode, pf_parameters_decode
+
+try:
+ from cellprofiler.preferences import get_max_workers
+except:
+ get_max_workers = lambda: 1
+
+min_number_of_chosen_seeds = 6
+
+MAX_NUMBER_OF_CHOSEN_SNAKES_NORMAL = 20
+MAX_NUMBER_OF_CHOSEN_SNAKES_SUPERFIT = 40
+max_number_of_chosen_snakes = 20
+
+SEARCH_LENGTH_NORMAL = 100
+SEARCH_LENGTH_SUPERFIT = 400
+
+#
+#
+# PROGRESS CALLBACKS
+#
+#
+
+ESTIMATED_CALCULATIONS_NUMBER_NORMAL = 3000.0
+ESTIMATED_CALCULATIONS_NUMBER_SUPERFIT = ESTIMATED_CALCULATIONS_NUMBER_NORMAL \
+ * SEARCH_LENGTH_SUPERFIT / SEARCH_LENGTH_NORMAL * 0.6
+estimated_calculations_number = ESTIMATED_CALCULATIONS_NUMBER_NORMAL
+
+callback_progress = None
+
+
+def show_progress(current_distance, calculation):
+ if calculation % 100 == 0:
+ logger.debug("Current distance: %f, Best: %f, Calc %d" % (current_distance, best_so_far, calculation))
+ if callback_progress is not None and calculation % (estimated_calculations_number / 50) == 0:
+ callback_progress(float(calculation) / estimated_calculations_number)
+
+
+#
+#
+# COST FUNCTION AND FITNESS
+#
+#
+
+best_so_far = 1
+calculations = 0
+best_3 = []
+
+
+def keep_3_best(partial_parameters, distance):
+ global best_3
+ best_3.append((distance, partial_parameters))
+ best_3.sort(key=lambda x: x[0])
+ best_3 = best_3[:3]
+ if best_3[0][0] == best_3[-1][0]:
+ best_3 = [best_3[0]]
+
+
+def distance_norm(fitnesses):
+ global calculations, best_so_far
+ # Mean-Squared Error
+ distance = norm((np.ones(fitnesses.shape) - fitnesses)) / np.sqrt(fitnesses.size)
+ best_so_far = min(best_so_far, distance)
+ calculations += 1
+
+ show_progress(distance, calculations)
+ return distance
+
+
+def grow_single_seed(seed, images, init_params, pf_param_vector):
+ pfsnake = PFSnake(seed, images, init_params)
+ return pfsnake.grow(pf_parameters_decode(pf_param_vector, pfsnake.orig_size_weight_list))
+
+
+def snakes_fitness(gt_snake_seed_pairs, images, parameters, pf_param_vector, debug=False):
+ gt_snake_grown_seed_pairs = [(gt_snake, grow_single_seed(seed, images, parameters, pf_param_vector)) for
+ gt_snake, seed in gt_snake_seed_pairs]
+
+ return np.array([pf_s.multi_fitness(gt_snake) for gt_snake, pf_s in gt_snake_grown_seed_pairs])
+
+
+#
+#
+# PREPARE DATA
+#
+#
+
+
+def get_gt_snake_seeds(gt_snake, number, max_radius, min_radius=0):
+ """
+ Create random seeds inside gt snake.
+ @type gt_snake: GTSnake
+ """
+ seed = Seed(gt_snake.centroid_x, gt_snake.centroid_y, "optimize_star_parameters")
+ seeds = [seed]
+ left = number
+ while left > 0:
+ random_seeds = Seeder.rand_seeds(max_radius, left, [seed], min_random_radius=min_radius)
+ inside_seeds = [s for s in random_seeds if gt_snake.is_inside(s.x, s.y)]
+ seeds += inside_seeds
+ left = number - (len(seeds) - 1)
+ min_radius /= 1.1 # make sure that it finish
+
+ return seeds
+
+
+def get_size_weight_list(params):
+ size_weight = params["segmentation"]["stars"]["sizeWeight"]
+ if isinstance(size_weight, float):
+ size_weight = [size_weight]
+ return size_weight
+
+
+def prepare_snake_seed_pairs(gt_snakes, initial_parameters):
+ radius = initial_parameters["segmentation"]["seeding"]["randomDiskRadius"] * initial_parameters["segmentation"][
+ "avgCellDiameter"]
+ for gt_snake in gt_snakes:
+ gt_snake.set_erosion(4)
+ gt_snake_seed_pairs = [(gt_snake, seed) for gt_snake in gt_snakes for seed in
+ get_gt_snake_seeds(gt_snake, number=3, max_radius=radius, min_radius=2 * radius / 3.0)]
+ random.shuffle(gt_snake_seed_pairs)
+ return gt_snake_seed_pairs
+
+
+def pf_get_distances(gt_snakes, images, initial_parameters, callback=None):
+ gt_snake_seed_pairs = prepare_snake_seed_pairs(gt_snakes, initial_parameters)
+ pick_seed_pairs = max(min_number_of_chosen_seeds, max_number_of_chosen_snakes /
+ len(initial_parameters["segmentation"]["stars"]["sizeWeight"]))
+ chosen_gt_snake_seed_pairs = gt_snake_seed_pairs[:pick_seed_pairs]
+
+ explore_cellstar(image=images.image, images=images, params=initial_parameters,
+ seeds=[sp[1] for sp in gt_snake_seed_pairs],
+ snakes=[])
+
+ def create_distance_function(pairs_to_use):
+ def distance(partial_parameters, debug=False):
+ # random.shuffle(pairs_to_use)
+ # randoms_pair = pairs_to_use[:pick_seed_pairs]
+ current_distance = distance_norm(
+ snakes_fitness(pairs_to_use, images, initial_parameters, partial_parameters, debug=debug)
+ )
+
+ # keep 3 best results
+ keep_3_best(partial_parameters, current_distance)
+
+ if callback is not None:
+ callback(partial_parameters, current_distance)
+
+ return current_distance
+
+ return distance
+
+ return create_distance_function(gt_snake_seed_pairs), create_distance_function(chosen_gt_snake_seed_pairs)
+
+
+#
+#
+# OPTIMIZATION
+#
+#
+
+def run(image, gt_snakes, precision, avg_cell_diameter, method='brute', initial_params=None, background_image=None,
+ ignore_mask=None):
+ global best_3, calculations
+ """
+ :param image: input image
+ :param gt_snakes: gt snakes label image
+ :param precision: if initial_params is None then it is used to calculate parameters
+ :param avg_cell_diameter: if initial_params is None then it is used to calculate parameters
+ :param method: optimization engine
+ :param initial_params: overrides precision and avg_cell_diameter
+ :return:
+ """
+ logger.info("Parameter fitting started...")
+ if initial_params is None:
+ params = default_parameters(segmentation_precision=precision, avg_cell_diameter=avg_cell_diameter)
+ else:
+ params = copy.deepcopy(initial_params)
+ images = ImageRepo(image, params)
+ images.background = background_image
+ if ignore_mask is not None:
+ images.apply_mask(ignore_mask)
+
+ start = time.clock()
+ best_3 = []
+ calculations = 0
+ best_arg, best_score = optimize(method, gt_snakes, images, params, precision, avg_cell_diameter)
+
+ best_params = pf_parameters_decode(best_arg, get_size_weight_list(params))
+
+ stop = time.clock()
+ logger.debug("Best: \n" + "\n".join([k + ": " + str(v) for k, v in sorted(best_params.iteritems())]))
+ logger.debug("Time: %d" % (stop - start))
+ logger.info("Parameter fitting finished with best score %f" % best_score)
+ return PFSnake.merge_parameters(params, best_params), best_arg, best_score
+
+
+def optimize(method_name, gt_snakes, images, params, precision, avg_cell_diameter):
+ global max_number_of_chosen_snakes, estimated_calculations_number
+
+ search_length = SEARCH_LENGTH_NORMAL
+ max_number_of_chosen_snakes = MAX_NUMBER_OF_CHOSEN_SNAKES_NORMAL
+ estimated_calculations_number = ESTIMATED_CALCULATIONS_NUMBER_NORMAL
+ if method_name == 'superfit':
+ # changes time limits to much longer
+ search_length = SEARCH_LENGTH_SUPERFIT
+ max_number_of_chosen_snakes = MAX_NUMBER_OF_CHOSEN_SNAKES_SUPERFIT
+ estimated_calculations_number = ESTIMATED_CALCULATIONS_NUMBER_SUPERFIT
+ method_name = 'brutemax3basin'
+ pass
+
+ encoded_params = pf_parameters_encode(params)
+ complete_distance, fast_distance = pf_get_distances(gt_snakes, images, params)
+ initial_distance = fast_distance(encoded_params)
+ initial_complete_distance = complete_distance(encoded_params)
+ logger.debug("Initial parameters complete-distance is (%f)." % (initial_complete_distance))
+ logger.debug("Initial parameters distance is (%f)." % (initial_distance))
+ logger.debug("Initial parameters are %s." % params)
+ if method_name.startswith("mp") and getattr(sys, "frozen", False) and sys.platform == 'win32':
+ # multiprocessing do not work then
+ method_name = "brutemaxbasin"
+ if method_name == "mp":
+ best_params_encoded, distance = multiproc_multitype_fitness(images.image, gt_snakes, precision,
+ avg_cell_diameter, "brutemaxbasin", params)
+ elif method_name == "mp_superfit":
+ best_params_encoded, distance = multiproc_multitype_fitness(images.image, gt_snakes, precision,
+ avg_cell_diameter, "superfit", params)
+ else:
+ if method_name == 'brute':
+ best_params_encoded, distance = optimize_brute(encoded_params, fast_distance)
+ elif method_name == 'brutemaxbasin':
+ best_params_encoded, distance = optimize_brute(encoded_params, fast_distance)
+ logger.debug("Best grid parameters distance is (%f)." % distance)
+ best_params_encoded, distance = optimize_basinhopping(best_params_encoded, fast_distance, time_percent=search_length)
+ elif method_name == 'brutemax3basin':
+ _, _ = optimize_brute(encoded_params, fast_distance)
+ logger.debug("Best grid parameters distance are %s." % str(zip(*best_3)[0]))
+ logger.debug("3-best grid parameters are %s." % str(zip(*best_3)[1]))
+
+ best_basins = []
+ for candidate in list(best_3):
+ best_basins.append(optimize_basinhopping(candidate[1], fast_distance, time_percent=search_length/3))
+ best_basins.sort(key=lambda x: x[1])
+
+ best_params_encoded, distance = best_basins[0]
+ elif method_name == 'basin':
+ best_params_encoded, distance = optimize_basinhopping(encoded_params, fast_distance)
+
+ complete_distance = complete_distance(best_params_encoded)
+ logger.debug("Final parameters complete-distance is (%f)." % (complete_distance))
+ if initial_complete_distance <= complete_distance:
+ logger.debug("Initial parameters (%f) are not worse than the best found (%f)." % (
+ initial_complete_distance, complete_distance))
+ return encoded_params, initial_complete_distance
+ else:
+ return best_params_encoded, complete_distance
+
+
+def optimize_brute(params_to_optimize, distance_function):
+ lower_bound = params_to_optimize - np.maximum(np.abs(params_to_optimize), 0.1)
+ upper_bound = params_to_optimize + np.maximum(np.abs(params_to_optimize), 0.1)
+
+ # introduce random shift (0,grid step) # max 20%
+ number_of_steps = 5
+ step = (upper_bound - lower_bound) / float(number_of_steps)
+ random_shift = np.array([random.random() * 2 / 10 for _ in range(len(lower_bound))])
+ lower_bound += random_shift * step
+ upper_bound += random_shift * step
+
+ logger.debug("Search range: " + str(zip(lower_bound, upper_bound)))
+ result = opt.brute(distance_function, zip(lower_bound, upper_bound), Ns=number_of_steps, disp=True, finish=None,
+ full_output=True)
+ logger.debug("Opt finished:" + str(result[:2]))
+ return result[0], result[1]
+
+
+def optimize_basinhopping(params_to_optimize, distance_function, time_percent=100):
+ minimizer_kwargs = {"method": "COBYLA"}
+ # bounds = ContourBounds
+ # minimizer_kwargs = {"method": "L-BFGS-B", "bounds" : zip(bounds.xmin,bounds.xmax)}
+ bounds = None
+ result = opt.basinhopping(distance_function, params_to_optimize, accept_test=bounds,
+ minimizer_kwargs=minimizer_kwargs, niter=35 * time_percent / 100)
+ logger.debug("Opt finished: " + str(result))
+ return result.x, result.fun
+
+
+#
+#
+# MULTIPROCESSING - MULTIPLE METHODS
+#
+#
+
+def general_multiproc_fitting(run_wrapper, *args):
+ result_queue = Queue()
+ update_queue = Queue()
+ workers_num = get_max_workers()
+
+ optimizers = [
+ Process(target=run_wrapper, args=(result_queue, update_queue) + args)
+ for _ in range(workers_num)]
+
+ for optimizer in optimizers:
+ optimizer.start()
+
+ optimizers_left = workers_num
+ results = []
+ while optimizers_left > 0:
+ time.sleep(0.1)
+ if not update_queue.empty() and callback_progress is not None:
+ callback_progress(update_queue.get())
+
+ if not result_queue.empty():
+ results.append(result_queue.get())
+ optimizers_left -= 1
+
+ for optimizer in optimizers:
+ optimizer.join()
+
+ sorted_results = sorted(results, key=lambda x: x[2])
+ logger.debug(str(sorted_results[0]))
+ return sorted_results[0][1], sorted_results[0][2]
+
+
+def run_wrapper(queue, update_queue, *args):
+ global callback_progress
+ random.seed() # reseed with random
+ callback_progress = lambda p: update_queue.put(p)
+ result = run(*args)
+ queue.put(result)
+
+
+def multiproc_multitype_fitness(image, gt_snakes, precision, avg_cell_diameter, method, init_params=None):
+ return general_multiproc_fitting(run_wrapper, image, gt_snakes, precision, avg_cell_diameter, method,
+ init_params)
diff --git a/cellstar/parameter_fitting/pf_rank_process.py b/cellstar/parameter_fitting/pf_rank_process.py
new file mode 100644
index 00000000..4e6caec5
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_rank_process.py
@@ -0,0 +1,326 @@
+# -*- coding: utf-8 -*-
+"""
+Pf rank process is the core of ranking parameters fitting.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import operator as op
+import random
+import time
+
+import numpy as np
+import scipy.optimize as opt
+from scipy.linalg import norm
+
+random.seed(1)
+np.random.seed(1)
+
+import logging
+
+logger = logging.getLogger(__name__)
+
+from cellstar.utils.params_util import *
+from cellstar.core.image_repo import ImageRepo
+from cellstar.core.snake_filter import SnakeFilter
+from cellstar.utils.debug_util import explore_cellstar
+from cellstar.parameter_fitting.pf_process import get_gt_snake_seeds, grow_single_seed, \
+ general_multiproc_fitting
+from cellstar.parameter_fitting.pf_rank_snake import PFRankSnake
+from cellstar.parameter_fitting.pf_auto_params import pf_parameters_encode, pf_rank_parameters_encode, \
+ pf_rank_parameters_decode, RankBounds
+from cellstar.parameter_fitting.pf_mutator import *
+
+#
+#
+# PROGRESS CALLBACKS
+#
+#
+ESTIMATED_CALCULATIONS_NUMBER = 20000.0
+callback_progress = None
+
+
+def show_progress(current_distance, calculation):
+ if calculations % 100 == 0:
+ logger.debug("Rank current: %f, Best: %f, Calc %d" % (current_distance, best_so_far, calculation))
+ if callback_progress is not None and calculation % (ESTIMATED_CALCULATIONS_NUMBER / 50) == 0:
+ callback_progress(float(calculation) / ESTIMATED_CALCULATIONS_NUMBER)
+
+#
+#
+# COST FUNCTION AND FITNESS
+#
+#
+
+best_so_far = 1000000000
+calculations = 0
+
+
+def maximal_distance(n):
+ l = n - 1
+ return l * (l + 1) * (l + 2) / 6.0 * 2
+
+
+def distance_smooth_norm(expected, result):
+ """
+ Calculates 2-norm from difference in fitness between expected and given snakes
+ @param expected: array of expected fitness
+ @param result: array of given fitness
+ @return:
+ """
+ global best_so_far, calculations
+ n = result.size
+ differences = abs(expected - result) ** 4 * np.arange(n * 2, 0, -2)
+ distance = norm(differences) / np.sqrt(n)
+
+ best_so_far = min(best_so_far, distance)
+ calculations += 1
+
+ show_progress(distance, calculations)
+ return distance
+
+
+def distance_norm_list(expected, result):
+ """
+ Calculates number of derangments between two sequences
+ @param expected: expected order
+ @param result: given order
+ @return:
+ """
+ global best_so_far, calculations
+ length = len(expected)
+ exp_position = dict([(obj, i) for (i, obj) in enumerate(expected)])
+ given_position = dict([(obj, i) for (i, obj) in enumerate(result)])
+ positions = enumerate(result)
+ distance = sum(
+ [abs(exp_position[obj] - i) ** 2 / (exp_position[obj] + 1) ** 2 for (i, obj) in positions]) / maximal_distance(
+ length) # scaling to [0,1]
+
+ best_so_far = min(best_so_far, distance)
+ calculations += 1
+
+ show_progress(distance, calculations)
+ return distance
+
+
+def calc_ranking(rank_snakes, pf_param_vector):
+ rank_params_decoded = pf_rank_parameters_decode(pf_param_vector)
+ fitness_order = sorted(rank_snakes, key=lambda x: -x.fitness)
+ ranking_order = sorted(rank_snakes, key=lambda x: x.calculate_ranking(rank_params_decoded))
+ return distance_norm_list(fitness_order, ranking_order)
+
+
+def calc_smooth_ranking(rank_snakes, pf_param_vector):
+ rank_params_decoded = pf_rank_parameters_decode(pf_param_vector)
+ fitness_order = np.array([r.fitness for r in sorted(rank_snakes, key=lambda x: -x.fitness)])
+ ranking_order = np.array(
+ [r.fitness for r in sorted(rank_snakes, key=lambda x: x.calculate_ranking(rank_params_decoded))])
+ return distance_smooth_norm(fitness_order, ranking_order)
+
+
+def pf_rank_get_ranking(rank_snakes, initial_parameters):
+ fitness = lambda partial_parameters, debug=False: \
+ calc_smooth_ranking(
+ rank_snakes,
+ partial_parameters
+ )
+
+ return fitness
+
+
+def filter_snakes_as_singles(parameters, images, snakes):
+ """
+ @type snakes: list[(GTSnake, PFRankSnake)]
+ """
+ filterer = SnakeFilter(images, parameters)
+ proper_snakes = [(gt, snake) for gt, snake in snakes if not filterer.is_single_snake_discarded(snake.grown_snake)]
+ logger.debug("Filtering left %d out of %d rank snakes" % (len(snakes) - len(proper_snakes), len(snakes)))
+ return proper_snakes
+
+
+#
+#
+# OPTIMIZATION
+#
+#
+
+def run_multiprocess(image, gt_snakes, precision=None, avg_cell_diameter=None, method='brute', initial_params=None,
+ background_image=None, ignore_mask=None):
+ """
+ :param gt_snakes: gt snakes label image
+ :param precision: if initial_params is None then it is used to calculate parameters
+ :param avg_cell_diameter: if initial_params is None then it is used to calculate parameters
+ :param method: optimization engine
+ :param initial_params: overrides precision and avg_cell_diameter
+ :return:
+ """
+ logger.info("Ranking parameter fitting (mp) started...")
+
+ if initial_params is None:
+ params = default_parameters(segmentation_precision=precision, avg_cell_diameter=avg_cell_diameter)
+ else:
+ params = copy.deepcopy(initial_params)
+ avg_cell_diameter = params["segmentation"]["avgCellDiameter"]
+
+ start = time.clock()
+ best_params, distance = multiproc_optimize((image, background_image, ignore_mask), gt_snakes, method, params)
+ best_params_full = PFRankSnake.merge_rank_parameters(params, best_params)
+ stop = time.clock()
+
+ logger.debug("Best: \n" + "\n".join([k + ": " + str(v) for k, v in sorted(best_params.iteritems())]))
+ logger.debug("Time: %d" % (stop - start))
+ logger.info("Ranking parameter fitting (mp) finished with best score %f" % distance)
+ return best_params_full, best_params, distance
+
+
+def run_singleprocess(image, gt_snakes, precision=None, avg_cell_diameter=None, method='brute', initial_params=None,
+ background_image=None, ignore_mask=None):
+ """
+ :param gt_snakes: gt snakes label image
+ :param precision: if initial_params is None then it is used to calculate parameters
+ :param avg_cell_diameter: if initial_params is None then it is used to calculate parameters
+ :param method: optimization engine
+ :param initial_params: overrides precision and avg_cell_diameter
+ :return:
+ """
+ global calculations
+ logger.info("Ranking parameter fitting started...")
+
+ if initial_params is None:
+ params = default_parameters(segmentation_precision=precision, avg_cell_diameter=avg_cell_diameter)
+ else:
+ params = copy.deepcopy(initial_params)
+ avg_cell_diameter = params["segmentation"]["avgCellDiameter"]
+
+ start = time.clock()
+
+ images = ImageRepo(image, params)
+ images.background = background_image
+ if ignore_mask is not None:
+ images.apply_mask(ignore_mask)
+
+ # prepare seed and grow snakes
+ encoded_star_params = pf_parameters_encode(params)
+ radius = params["segmentation"]["seeding"]["randomDiskRadius"] * params["segmentation"]["avgCellDiameter"]
+ radius_big = params["segmentation"]["avgCellDiameter"] * 1.5
+ gt_snake_seed_pairs = [(gt_snake, seed) for gt_snake in gt_snakes for seed in
+ get_gt_snake_seeds(gt_snake, max_radius=radius, number=8, min_radius=2 * radius / 3.0)
+ + get_gt_snake_seeds(gt_snake, max_radius=radius, number=8, min_radius=4 * radius / 5.0)
+ + get_gt_snake_seeds(gt_snake, max_radius=radius_big, number=8,
+ min_radius=3 * radius_big / 4.0)
+ ]
+
+ gt_snake_grown_seed_pairs = \
+ [(gt_snake, grow_single_seed(seed, images, params, encoded_star_params)) for gt_snake, seed in
+ gt_snake_seed_pairs]
+
+ gt_snake_grown_seed_pairs_all = reduce(op.add,
+ [PFRankSnake.create_all(gt, grown, params) for (gt, grown) in
+ gt_snake_grown_seed_pairs])
+
+ # gt_snake_grown_seed_pairs_filtered = filter_snakes_as_singles(params, images, gt_snake_grown_seed_pairs_all)
+ gt_snake_grown_seed_pairs_filtered = gt_snake_grown_seed_pairs_all
+
+ # gts_snakes_with_mutations = add_mutations(gt_snake_grown_seed_pairs_all, avg_cell_diameter)
+ gts_snakes_with_mutations = gt_snake_grown_seed_pairs_filtered
+ ranked_snakes = zip(*gts_snakes_with_mutations)[1]
+
+ explore_cellstar(image=images.image, images=images, params=params,
+ seeds=[sp[1].grown_snake.seed for sp in gts_snakes_with_mutations],
+ snakes=[sp[1].grown_snake for sp in gts_snakes_with_mutations])
+
+ calculations = 0
+ best_params_encoded, distance = optimize(
+ method,
+ pf_rank_parameters_encode(params),
+ pf_rank_get_ranking(ranked_snakes, params)
+ )
+
+ stop = time.clock()
+
+ best_params_org = pf_rank_parameters_decode(best_params_encoded)
+ best_params_full = PFRankSnake.merge_rank_parameters(params, best_params_org)
+
+ explore_cellstar(image=images.image, images=images, params=best_params_full,
+ seeds=[sp[1].grown_snake.seed for sp in gts_snakes_with_mutations],
+ snakes=[sp[1].grown_snake for sp in gts_snakes_with_mutations])
+
+ logger.debug("Best: \n" + "\n".join([k + ": " + str(v) for k, v in sorted(best_params_org.iteritems())]))
+ logger.debug("Time: %d" % (stop - start))
+ logger.info("Ranking parameter fitting finished with best score %f" % distance)
+ return best_params_full, best_params_org, distance
+
+
+#
+#
+# OPTIMISATION METHODS
+#
+#
+
+def optimize(method_name, encoded_params, distance_function):
+ initial_distance = distance_function(encoded_params)
+ logger.debug("Initial parameters distance is (%f)." % initial_distance)
+ if method_name == 'brute':
+ best_params_encoded, distance = optimize_brute(encoded_params, distance_function)
+ elif method_name == 'brutemaxbasin' or method_name == 'superfit':
+ best_params_encoded, distance = optimize_brute(encoded_params, distance_function)
+ logger.debug("Best grid parameters distance is (%f)." % distance)
+ best_params_encoded, distance = optimize_basinhopping(best_params_encoded, distance_function)
+ else:
+ raise Exception("No such optimization method.")
+
+ if initial_distance <= distance:
+ logger.debug("Initial parameters (%f) are not worse than the best found (%f)." % (initial_distance, distance))
+ return encoded_params, initial_distance
+ else:
+ return best_params_encoded, distance
+
+
+def optimize_brute(params_to_optimize, distance_function):
+ lower_bound = np.zeros(len(params_to_optimize), dtype=float)
+ upper_bound = np.ones(len(params_to_optimize), dtype=float)
+
+ # introduce random shift (0,grid step) # max 10%
+ number_of_steps = 6
+ step = (upper_bound - lower_bound) / float(number_of_steps)
+ random_shift = np.array([random.random() * 1 / 10 for _ in range(len(lower_bound))], dtype=float)
+ lower_bound += random_shift * step
+ upper_bound += random_shift * step
+
+ start = time.clock()
+ result = opt.brute(distance_function, zip(lower_bound, upper_bound), finish=None, Ns=number_of_steps, disp=True,
+ full_output=True)
+ elapsed = time.clock() - start
+
+ logger.debug("Opt finished: " + str(result[:2]) + " Elapsed[s]: " + str(elapsed))
+
+ return result[0], result[1]
+
+
+def optimize_basinhopping(params_to_optimize, distance_function):
+ bounds = RankBounds
+ # minimizer_kwargs = {"method": "COBYLA", bounds=bounds}
+ minimizer_kwargs = {"method": "L-BFGS-B", "bounds": zip(bounds.xmin, bounds.xmax)}
+ result = opt.basinhopping(distance_function, params_to_optimize, accept_test=bounds,
+ minimizer_kwargs=minimizer_kwargs, niter=170)
+ logger.debug("Opt finished: " + str(result))
+ return result.x, result.fun
+
+
+#
+#
+# MULTIPROCESSING METHODS
+#
+#
+
+def run_wrapper(queue, update_queue, images, gt_snakes, method, params):
+ global callback_progress
+ random.seed() # reseed with random
+ callback_progress = lambda p: update_queue.put(p)
+ result = run_singleprocess(images[0], gt_snakes, method=method, initial_params=params, background_image=images[1],
+ ignore_mask=images[2])
+ queue.put(result)
+
+
+def multiproc_optimize(images, gt_snakes, method='brute', initial_params=None):
+ return general_multiproc_fitting(run_wrapper, images, gt_snakes, method, initial_params)
diff --git a/cellstar/parameter_fitting/pf_rank_snake.py b/cellstar/parameter_fitting/pf_rank_snake.py
new file mode 100644
index 00000000..ca418642
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_rank_snake.py
@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+PFRankSnake represents one ground_truth contour for ranking parameters fitting.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import copy
+import random
+
+random.seed(1) # make it deterministic
+
+from cellstar.core.polar_transform import PolarTransform
+from cellstar.parameter_fitting.pf_snake import PFSnake
+import pf_mutator
+
+
+class PFRankSnake(object):
+ def __init__(self, gt_snake, grown_snake, avg_cell_diameter, params):
+ self.gt_snake = gt_snake
+ self.grown_snake = grown_snake
+ self.avg_cell_diameter = avg_cell_diameter
+ self.initial_parameters = params
+ self.fitness = PFSnake.fitness_with_gt(grown_snake, gt_snake)
+ self.rank_vector = grown_snake.properties_vector(avg_cell_diameter)
+ self.polar_transform = PolarTransform.instance(params["segmentation"]["avgCellDiameter"],
+ params["segmentation"]["stars"]["points"],
+ params["segmentation"]["stars"]["step"],
+ params["segmentation"]["stars"]["maxSize"])
+
+ @staticmethod
+ def create_all(gt_snake, grown_pf_snake, params):
+ return [(gt_snake, PFRankSnake(gt_snake, snake, grown_pf_snake.avg_cell_diameter, params)) for snake in
+ grown_pf_snake.snakes]
+
+ def create_mutation(self, dilation, random_poly=False):
+ if random_poly:
+ mutant = pf_mutator.create_poly_mutation(self.grown_snake, self.polar_transform, dilation)
+ else:
+ mutant = pf_mutator.create_mutation(self.grown_snake, self.polar_transform, dilation)
+ return PFRankSnake(self.gt_snake, mutant, self.avg_cell_diameter, self.initial_parameters)
+
+ @staticmethod
+ def merge_rank_parameters(initial_parameters, new_params):
+ params = copy.deepcopy(initial_parameters)
+ for k, v in new_params.iteritems():
+ params["segmentation"]["ranking"][k] = v
+
+ return params
+
+ def merge_parameters_with_me(self, new_params):
+ return PFRankSnake.merge_rank_parameters(self.initial_parameters, new_params)
+
+ def calculate_ranking(self, ranking_params):
+ return self.grown_snake.star_rank(ranking_params, self.avg_cell_diameter)
diff --git a/cellstar/parameter_fitting/pf_runner.py b/cellstar/parameter_fitting/pf_runner.py
new file mode 100644
index 00000000..ebb73097
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_runner.py
@@ -0,0 +1,98 @@
+# -*- coding: utf-8 -*-
+"""
+Entry point for running fitting process both parameter sets: contour and ranking.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+import logging
+
+import sys
+import numpy as np
+import scipy as sp
+
+import cellstar.parameter_fitting.pf_process as pf_process
+import cellstar.parameter_fitting.pf_rank_process as pf_rank
+from cellstar.parameter_fitting.pf_snake import GTSnake
+
+try:
+ from cellprofiler.preferences import get_max_workers
+except:
+ get_max_workers = lambda: 1
+
+logger = logging.getLogger(__name__)
+
+
+def single_mask_to_snake(bool_mask, seed=None):
+ return GTSnake(bool_mask, seed)
+
+
+def gt_label_to_snakes(components):
+ num_components = components.max()
+ return [single_mask_to_snake(components == label) for label in range(1, num_components + 1)]
+
+
+def image_to_label(image):
+ values = np.unique(image)
+ if len(values) == 2: # it is a mask
+ components, num_components = sp.ndimage.label(image, np.ones((3, 3)))
+ return components
+ else: # remap labels to [1..] values
+ curr = 1
+ label_image = image.copy()
+ for v in values[1:]: # zero is ignored
+ label_image[image == v] = curr
+ curr += 1
+ return label_image
+
+
+def run_pf(input_image, background_image, ignore_mask_image, gt_label, parameters, precision, avg_cell_diameter,
+ callback_progress=None):
+ """
+ :param input_image:
+ :param gt_label:
+ :param parameters:
+ :return: Best complete parameters settings, best distance
+ """
+
+ gt_mask = image_to_label(gt_label)
+ pf_process.callback_progress = callback_progress
+
+ gt_snakes = gt_label_to_snakes(gt_mask)
+ if get_max_workers() > 1:
+ best_complete_params, _, best_score = pf_process.run(input_image, gt_snakes, precision=precision,
+ avg_cell_diameter=avg_cell_diameter, initial_params=parameters,
+ method='mp', background_image=background_image,
+ ignore_mask=ignore_mask_image)
+ else:
+ best_complete_params, _, best_score = pf_process.run(input_image, gt_snakes, precision=precision,
+ avg_cell_diameter=avg_cell_diameter, initial_params=parameters,
+ method='brutemaxbasin', background_image=background_image,
+ ignore_mask=ignore_mask_image)
+
+ return best_complete_params, best_score
+
+
+def run_rank_pf(input_image, background_image, ignore_mask_image, gt_mask, parameters, callback_progress=None):
+ """
+ :return: Best complete parameters settings, best distance
+ """
+
+ gt_mask = image_to_label(gt_mask)
+ pf_rank.callback_progress = callback_progress
+
+ gt_snakes = gt_label_to_snakes(gt_mask)
+ if get_max_workers() > 1 and not (getattr(sys, "frozen", False) and sys.platform == 'win32'):
+ # multiprocessing do not work if frozen on win32
+ best_complete_params, _, best_score = pf_rank.run_multiprocess(input_image, gt_snakes,
+ initial_params=parameters,
+ method='brutemaxbasin',
+ background_image=background_image,
+ ignore_mask=ignore_mask_image)
+ else:
+ best_complete_params, _, best_score = pf_rank.run_singleprocess(input_image, gt_snakes,
+ initial_params=parameters,
+ method='brutemaxbasin',
+ background_image=background_image,
+ ignore_mask=ignore_mask_image)
+
+ return best_complete_params, best_score
\ No newline at end of file
diff --git a/cellstar/parameter_fitting/pf_snake.py b/cellstar/parameter_fitting/pf_snake.py
new file mode 100644
index 00000000..2756e88e
--- /dev/null
+++ b/cellstar/parameter_fitting/pf_snake.py
@@ -0,0 +1,149 @@
+# -*- coding: utf-8 -*-
+"""
+PFSnake represents one grown from a seed contour within a ground truth contour used in contour parameters fitting.
+GTSnake represents one ground truth contour.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import copy
+import random
+
+random.seed(1) # make it deterministic
+import numpy as np
+import scipy.ndimage.morphology as morph
+import scipy.ndimage.measurements as measure
+
+from cellstar.utils.calc_util import to_int
+from cellstar.core.seed import Seed
+from cellstar.core.snake import Snake
+from cellstar.core.polar_transform import PolarTransform
+
+
+class PFSnake(object):
+ def __init__(self, seed, image_repo, params, best_snake=None):
+ if seed is not None:
+ self.fit = 0.0
+ self.seed = seed
+ self.snakes = []
+ self.images = image_repo
+ self.initial_parameters = params
+ self.point_number = params["segmentation"]["stars"]["points"]
+ self.orig_size_weight_list = params["segmentation"]["stars"]["sizeWeight"]
+
+ if isinstance(self.orig_size_weight_list, float):
+ self.orig_size_weight_list = [self.orig_size_weight_list]
+ self.avg_cell_diameter = params["segmentation"]["avgCellDiameter"]
+ self.step = params["segmentation"]["stars"]["step"]
+ self.max_size = params["segmentation"]["stars"]["maxSize"]
+ self.polar_transform = PolarTransform.instance(params["segmentation"]["avgCellDiameter"],
+ params["segmentation"]["stars"]["points"],
+ params["segmentation"]["stars"]["step"],
+ params["segmentation"]["stars"]["maxSize"])
+
+ self.best_snake = best_snake
+
+ @staticmethod
+ def merge_parameters(initial_parameters, new_params):
+ params = copy.deepcopy(initial_parameters)
+ for k, v in new_params.iteritems():
+ params["segmentation"]["stars"][k] = v
+
+ return params
+
+ def merge_parameters_with_me(self, new_params):
+ return PFSnake.merge_parameters(self.initial_parameters, new_params)
+
+ def grow(self, supplementary_parameters=None):
+ if supplementary_parameters is None:
+ new_parameters = copy.deepcopy(self.initial_parameters)
+ else:
+ new_parameters = self.merge_parameters_with_me(supplementary_parameters)
+
+ s = Snake.create_from_seed(new_parameters, self.seed, self.point_number, self.images)
+
+ size_weight_list = new_parameters["segmentation"]["stars"]["sizeWeight"]
+ snakes_to_grow = [(copy.copy(s), w) for w in size_weight_list]
+
+ for snake, weight in snakes_to_grow:
+ snake.grow(size_weight=weight, polar_transform=self.polar_transform)
+ snake.evaluate(self.polar_transform)
+
+ self.snakes = [grown_snake for grown_snake, _ in snakes_to_grow]
+ self.best_snake = sorted(snakes_to_grow, key=lambda (sn, _): sn.rank)[0][0]
+
+ return self
+
+ def extract_total_mask(self, tot_shape):
+ return PFSnake.extract_total_mask_of_snake(self.best_snake, tot_shape)
+
+ @staticmethod
+ def extract_total_mask_of_snake(snake, tot_shape):
+ yx = snake.in_polygon_yx
+ in_polygon = snake.in_polygon
+ in_polygon_bounds = np.array([yx, np.array(yx) + in_polygon.shape]).flatten()
+
+ mask = np.zeros(tot_shape, dtype=bool)
+ mask[in_polygon_bounds[0]:in_polygon_bounds[2], in_polygon_bounds[1]:in_polygon_bounds[3]] = in_polygon
+
+ return mask
+
+ @staticmethod
+ def gt_snake_intersection(snake, gt):
+ yx = snake.in_polygon_yx
+ in_polygon = snake.in_polygon
+ in_polygon_bounds = np.array([yx, np.array(yx) + in_polygon.shape]).flatten()
+ intersection_local = gt.binary_mask[in_polygon_bounds[0]:in_polygon_bounds[2],
+ in_polygon_bounds[1]:in_polygon_bounds[3]] * in_polygon
+ return np.count_nonzero(intersection_local)
+
+ @staticmethod
+ def out_of_gt_penalty(snake_area, gt_snake_area, intersection):
+ snake_less_gt = snake_area - intersection
+ snake_less_gt_percent = snake_less_gt / gt_snake_area * 100
+ if snake_less_gt_percent < 20:
+ return 1
+ elif snake_less_gt_percent < 80:
+ return 1.3
+ else:
+ return 2
+
+ @staticmethod
+ def fitness_with_gt(snake, gt_snake):
+ intersection = PFSnake.gt_snake_intersection(snake, gt_snake)
+ return intersection / (
+ gt_snake.area + (snake.area - intersection) * PFSnake.out_of_gt_penalty(snake.area, gt_snake.area,
+ intersection))
+
+ def multi_fitness(self, gt_snake):
+ return max([PFSnake.fitness_with_gt(pf_snake, gt_snake) for pf_snake in self.snakes])
+
+
+class GTSnake(object):
+ def __init__(self, binary_mask, seed=None):
+ self.binary_mask = binary_mask
+ self.eroded_mask = morph.binary_erosion(binary_mask, np.ones((3, 3)))
+ self.area = np.count_nonzero(self.binary_mask)
+ if seed is not None:
+ self.seed = seed
+ self.centroid_x, self.centroid_y = seed.x, seed.y
+ else:
+ self.calculate_centroids(binary_mask)
+ self.seed = Seed(self.centroid_x, self.centroid_y, "gt_snake")
+
+ def calculate_centroids(self, binary_mask):
+ self.centroid_y, self.centroid_x = measure.center_of_mass(binary_mask, binary_mask, [1])[0]
+
+ def set_erosion(self, size):
+ self.eroded_mask = morph.binary_erosion(self.binary_mask, np.ones((size, size)))
+
+ def is_inside(self, x, y):
+ """
+ Check if seed is inside of eroded mask.
+ @type seed: Seed
+ """
+ x = to_int(x)
+ y = to_int(y)
+ if x < 0 or x >= self.eroded_mask.shape[1] or y < 0 or y >= self.eroded_mask.shape[0]:
+ return False
+ return self.eroded_mask[y, x]
diff --git a/cellstar/segmentation.py b/cellstar/segmentation.py
new file mode 100644
index 00000000..303b6ab9
--- /dev/null
+++ b/cellstar/segmentation.py
@@ -0,0 +1,246 @@
+# -*- coding: utf-8 -*-
+"""
+Segmentation is a main entry point for CellStar segmentation.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import ast
+import logging
+
+logger = logging.getLogger(__name__)
+from copy import copy
+
+from cellstar.utils.params_util import *
+from cellstar.core.image_repo import ImageRepo
+from cellstar.utils.params_util import default_parameters
+from cellstar.utils import image_util, debug_util
+from cellstar.core.seeder import Seeder
+from cellstar.core.snake import Snake
+from cellstar.core.snake_filter import SnakeFilter
+from cellstar.core.polar_transform import PolarTransform
+from cellstar.parameter_fitting.pf_auto_params import rank_parameters_range as rank_auto_params
+from cellstar.parameter_fitting.pf_auto_params import parameters_range as snake_auto_params
+
+
+class Segmentation(object):
+ def __init__(self, segmentation_precision=9, avg_cell_diameter=35):
+ self.parameters = default_parameters(segmentation_precision, avg_cell_diameter)
+ self.images = None
+ self.all_seeds = []
+ self.seeds = []
+ self.grown_seeds = set() # seeds from which we already have snakes
+ self.snakes = []
+ self.new_snakes = []
+ self._seeder = None
+ self._filter = None
+ self.polar_transform = PolarTransform.instance(self.parameters["segmentation"]["avgCellDiameter"],
+ self.parameters["segmentation"]["stars"]["points"],
+ self.parameters["segmentation"]["stars"]["step"],
+ self.parameters["segmentation"]["stars"]["maxSize"])
+ self.debug_output_image_path = None
+
+ @property
+ def seeder(self):
+ if self._seeder is None:
+ self.init_seeder()
+
+ return self._seeder
+
+ @property
+ def filter(self):
+ if self._filter is None:
+ self.init_filter()
+
+ return self._filter
+
+ def clear_lists(self):
+ self.all_seeds = []
+ self.seeds = []
+ self.grown_seeds = set()
+ self.snakes = []
+ self.new_snakes = []
+
+ def set_frame(self, frame):
+ # Extract previous background
+ prev_background = None
+ if self.images is not None:
+ prev_background = self.images.background
+ # Initialize new image repository for new frame
+ self.images = ImageRepo(frame, self.parameters)
+ # One background per whole segmentation
+ if prev_background is not None:
+ self.images.background = prev_background
+
+ # Clear all temporary results
+ self.clear_lists()
+
+ def set_background(self, background):
+ self.images._background = background
+
+ def set_mask(self, ignore_mask):
+ if ignore_mask is not None:
+ self.images.apply_mask(ignore_mask)
+
+ def init_seeder(self):
+ self._seeder = Seeder(self.images, self.parameters)
+
+ def init_filter(self):
+ self._filter = SnakeFilter(self.images, self.parameters)
+
+ def decode_auto_params(self, text):
+ """
+ Decode automatic parameters from text and apply to self.
+
+ @param text: parameters denoted as python list
+ @return true if parsing was successful
+ """
+ return Segmentation.decode_auto_params_into(self.parameters, text)
+
+ @staticmethod
+ def decode_auto_params_into(complete_params, text):
+ """
+ Decode automatic parameters from text and apply to self.
+
+ @param text: parameters denoted as python list
+ @return true if parsing was successful
+ """
+ new_stars = copy(complete_params["segmentation"]["stars"])
+ new_ranking = copy(complete_params["segmentation"]["ranking"])
+ try:
+ all_params = ast.literal_eval(text)
+ snake_params = all_params[0]
+ rank_params = all_params[1]
+ if len(snake_params) != len(snake_auto_params) or len(rank_params) != len(rank_auto_params):
+ raise Exception("text invalid: list size not compatible")
+
+ for name in sorted(snake_auto_params.keys()):
+ val = snake_params[0]
+ if name == "sizeWeight": # value to list
+ original = complete_params["segmentation"]["stars"]["sizeWeight"]
+ val = list(np.array(original) * (val/np.mean(original)))
+
+ new_stars[name] = val
+ snake_params = snake_params[1:]
+
+ for name in sorted(rank_auto_params.keys()):
+ new_ranking[name] = rank_params[0]
+ rank_params = rank_params[1:]
+ except:
+ return False
+
+ complete_params["segmentation"]["stars"] = new_stars
+ complete_params["segmentation"]["ranking"] = new_ranking
+ return True
+
+ @staticmethod
+ def encode_auto_params_from_all_params(parameters):
+ snake_auto_params_values = []
+ for name in sorted(snake_auto_params.keys()):
+ val = parameters["segmentation"]["stars"][name]
+ if name == "sizeWeight": # list to mean value
+ original = parameters["segmentation"]["stars"]["sizeWeight"]
+ val = np.mean(original)
+ snake_auto_params_values.append(val)
+
+ rank_auto_params_values = [parameters["segmentation"]["ranking"][name]
+ for name in sorted(rank_auto_params.keys())]
+ auto_values_list = [snake_auto_params_values, rank_auto_params_values]
+
+ return str(auto_values_list)
+
+ def encode_auto_params(self):
+ return Segmentation.encode_auto_params_from_all_params(self.parameters)
+
+ def pre_process(self):
+ # background getter is never None but it creation background only if not existant
+ if self.images.background is None:
+ self.images.calculate_background()
+
+ self.images.calculate_brighter_original()
+ self.images.calculate_darker_original()
+ self.images.calculate_clean_original()
+ self.images.calculate_forebackground_masks()
+
+ self.images.calculate_clean()
+ self.images.calculate_brighter()
+ self.images.calculate_darker()
+
+ self.images.calculate_cell_border_content_mask()
+
+ def find_seeds(self, exclude):
+ self.seeds = self.seeder.find_seeds(self.snakes, self.all_seeds, exclude_current_segments=exclude)
+ self.all_seeds += self.seeds
+
+ def snakes_from_seeds(self):
+ self.new_snakes = [
+ Snake.create_from_seed(
+ self.parameters, seed, self.parameters["segmentation"]["stars"]["points"], self.images
+ )
+ for seed in self.seeds if seed not in self.grown_seeds
+ ]
+ for seed in self.seeds:
+ if seed not in self.grown_seeds:
+ self.grown_seeds.add(seed)
+
+ def grow_snakes(self):
+ grown_snakes = []
+ size_weights = self.parameters["segmentation"]["stars"]["sizeWeight"]
+ logger.debug("%d snakes seeds to grow with %d weights options -> %d snakes to calculate"%(len(self.new_snakes), len(size_weights), len(self.new_snakes) * len(size_weights)))
+ for snake in self.new_snakes:
+ best_snake = None
+ for weight in size_weights:
+ curr_snake = copy(snake)
+
+ curr_snake.grow(weight, self.polar_transform)
+ curr_snake.evaluate(self.polar_transform)
+
+ if best_snake is None:
+ best_snake = curr_snake
+ else:
+ if curr_snake.rank < best_snake.rank:
+ best_snake = curr_snake
+
+ grown_snakes.append(best_snake)
+
+ self.new_snakes = grown_snakes
+
+ def evaluate_snakes(self):
+ for snake in self.new_snakes:
+ snake.evaluate(self.polar_transform)
+
+ def filter_snakes(self):
+ self.snakes = self.filter.filter(self.snakes + self.new_snakes)
+ self.new_snakes = []
+
+ def debug_images(self):
+ if self.debug_output_image_path is not None:
+ debug_util.debug_image_path = self.debug_output_image_path
+ debug_util.images_repo_save(self.images)
+
+ def debug_seeds(self, step):
+ if self.debug_output_image_path is not None:
+ image_util.debug_image_path = self.debug_output_image_path
+ debug_util.draw_seeds(self.all_seeds, self.images.image, title=str(step))
+
+ def run_one_step(self, step):
+ logger.debug("find_seeds")
+ self.find_seeds(step > 0)
+ self.debug_seeds(step)
+ logger.debug("snake_from_seeds")
+ self.snakes_from_seeds()
+ logger.debug("grow_snakes")
+ self.grow_snakes()
+ logger.debug("filter_snakes")
+ debug_util.draw_snakes(self.images.image, self.snakes + self.new_snakes, it=step)
+ self.filter_snakes()
+ logger.debug("done")
+
+ def run_segmentation(self):
+ logger.debug("preproces...")
+ self.pre_process()
+ self.debug_images()
+ debug_util.explore_cellstar(self)
+ for step in range(self.parameters["segmentation"]["steps"]):
+ self.run_one_step(step)
+ return self.images.segmentation, self.snakes
\ No newline at end of file
diff --git a/cellstar/utils/__init__.py b/cellstar/utils/__init__.py
new file mode 100644
index 00000000..db21e07f
--- /dev/null
+++ b/cellstar/utils/__init__.py
@@ -0,0 +1,7 @@
+# -*- coding: utf-8 -*-
+"""
+Utilities package contain a number of general functions used in Cell Star segmentation process.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+__all__ = ["calc_util", "cluster_util", "image_util", "params_util", "python_util"]
diff --git a/cellstar/utils/calc_util.py b/cellstar/utils/calc_util.py
new file mode 100644
index 00000000..1d6f33be
--- /dev/null
+++ b/cellstar/utils/calc_util.py
@@ -0,0 +1,243 @@
+# -*- coding: utf-8 -*-
+"""
+Calculation package contains a number of functions used in contour grow and evaluation.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import math
+
+import numpy as np
+import scipy.ndimage as sp_image
+
+from index import Index
+
+
+def euclidean_norm((x1, y1), (x2, y2)):
+ return math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2)
+
+
+def interpolate_radiuses(values_mask, length, values):
+ """
+ Fill values with linear interpolation using values_mask values.
+ @type values_mask: np.ndarray
+ @param values_mask: mask of existing values
+ @type values: np.ndarray
+ @type length: int
+ """
+ cumlengths = np.where(values_mask)[0]
+ if len(cumlengths) > 0:
+ cumlengths_loop = np.append(cumlengths, cumlengths[0] + int(length))
+ for i in range(len(cumlengths)):
+ # Find left and right boundary in existing values.
+ left_interval_boundary = cumlengths_loop[i]
+ right_interval_boundary = cumlengths_loop[i + 1] % length
+
+ # Length of the interpolated interval.
+ interval_length = cumlengths_loop[i + 1] - left_interval_boundary - 1
+
+ # Interpolate for every point in the interval.
+ for k in range(1, interval_length + 1):
+ interpolated = (left_interval_boundary + k) % length
+
+ new_val = round(values[left_interval_boundary] +
+ (values[right_interval_boundary] - values[left_interval_boundary]) *
+ k / (interval_length + 1.0)) # (interval_length + 1.0)
+
+ values[interpolated] = new_val # min(values[interpolated], new_val)
+
+
+def loop_connected_components(mask):
+ """
+ @type mask: np.ndarray
+ @rtype (np.ndarray, np.ndarray, np.ndarray)
+ """
+
+ c = np.array([])
+ init = np.array([])
+ fin = np.array([])
+
+ if mask.sum() > 0:
+ labeled = sp_image.label(mask)[0]
+ components = sp_image.measurements.find_objects(labeled)
+ c_fin = [(s[0].stop - s[0].start, s[0].stop - 1) for s in components]
+ if len(c_fin) > 1 and mask[0] and mask[-1]:
+ c_fin[0] = c_fin[0][0] + c_fin[-1][0], c_fin[0][1]
+ c_fin = c_fin[0:-1]
+
+ c, fin = zip(*c_fin)
+ c = np.array(c, dtype=int)
+ fin = np.array(fin, dtype=int)
+ init = (fin - c + 1) % mask.shape[0]
+ return c, init, fin
+
+
+def unstick_contour(edgepoints, unstick_coeff):
+ """
+ Removes edgepoints near previously discarded points.
+ @type edgepoints: list[bool]
+ @param edgepoints: current edgepoint list
+ @type unstick_coeff: float
+ @param unstick_coeff
+ @return: filtered edgepoints
+ """
+ (n, init, end) = loop_connected_components(np.logical_not(edgepoints))
+ filtered = np.copy(edgepoints)
+ n_edgepoint = len(edgepoints)
+ for size, s, e in zip(n, init, end):
+ for j in range(1, int(size * unstick_coeff + 0.5) + 1):
+ filtered[(e + j) % n_edgepoint] = 0
+ filtered[(s - j) % n_edgepoint] = 0
+ return filtered
+
+
+def sub2ind(dim, (x, y)):
+ return x + y * dim
+
+
+def get_gradient(im, index, border_thickness_steps):
+ """
+ Fun. calc. radial gradient including thickness of cell edges
+ @param im: image (for which grad. will be calc.)
+ @param index: indices of pixes sorted by polar coordinates (alpha, radius)
+ @param border_thickness_steps: number of steps to cop. grad. - depends on cell border thickness
+ @return: gradient matrix for cell
+ """
+ # index of axis used to find max grad.
+ max_gradient_along_axis = 2
+
+ # preparing the image limits (called subimage) for which grad. will be computed
+ radius_lengths, angles = index.shape[:2]
+
+ # matrix init
+ # for each single step for each border thick. separated grad. is being computed
+ # at the end the max. grad values are returned (for all steps of thickness)
+ border_thickness_steps = int(border_thickness_steps)
+ gradients_for_steps = np.zeros((radius_lengths, angles, border_thickness_steps), dtype=np.float64)
+
+ # for every step of thickness:
+ for border_thickness_step in range(1, int(border_thickness_steps) + 1):
+ # find beg. and end indices of input matrix for which the gradient will be computed
+ matrix_end = radius_lengths - border_thickness_step
+ matrix_start = border_thickness_step
+
+ # find beg. and end indices of pix. for which the gradient will be computed
+ starting_index = index[:matrix_end, :]
+ ending_index = index[matrix_start:, :]
+
+ # find internal in matrix where computed gradient will go
+ intersect_start = int(math.ceil(border_thickness_step / 2.0))
+ intersect_end = int(intersect_start + matrix_end)
+
+ # comp. current gradient for selected (sub)image
+ current_step_gradient = im[Index.to_numpy(ending_index)] - im[Index.to_numpy(starting_index)]
+ current_step_gradient /= np.sqrt(border_thickness_step)
+
+ # save gradient to previously determined place in results matrix
+ gradients_for_steps[intersect_start:intersect_end, :, border_thickness_step - 1] = current_step_gradient
+
+ return gradients_for_steps.max(axis=max_gradient_along_axis)
+
+
+def extend_slices(my_slices, extension):
+ def extend_slice(my_slice, extend):
+ max_len = 100000
+ ind = (max(0, my_slice.indices(max_len)[0] - extend), my_slice.indices(max_len)[1] + extend)
+ return slice(*ind)
+
+ return extend_slice(my_slices[0], extension), extend_slice(my_slices[1], extension)
+
+
+def inslice_point(point_yx_in_slice, slices):
+ y = point_yx_in_slice[0]
+ x = point_yx_in_slice[1]
+ max_len = 1000000
+ return y - slices[0].indices(max_len)[0], x - slices[1].indices(max_len)[0]
+
+
+def unslice_point(point_yx_in_slice, slices):
+ y = point_yx_in_slice[0]
+ x = point_yx_in_slice[1]
+ max_len = 1000000
+ return y + slices[0].indices(max_len)[0], x + slices[1].indices(max_len)[0]
+
+
+def get_cartesian_bounds(polar_coordinate_boundary, origin_x, origin_y, polar_transform):
+ polygon_x, polygon_y = polar_to_cartesian(polar_coordinate_boundary, origin_x, origin_y, polar_transform)
+ x1 = int(max(0, math.floor(min(polygon_x))))
+ x2 = int(math.ceil(max(polygon_x)) + 1)
+ y1 = int(max(0, math.floor(min(polygon_y))))
+ y2 = int(math.ceil(max(polygon_y)) + 1)
+ return slice(y1, y2), slice(x1, x2)
+
+
+def polar_to_cartesian(polar_coordinate_boundary, origin_x, origin_y, polar_transform):
+ t = polar_transform.t
+ step = polar_transform.step
+ px = origin_x + step * polar_coordinate_boundary * np.cos(t.T)
+ py = origin_y + step * polar_coordinate_boundary * np.sin(t.T)
+
+ return px, py
+
+
+def mask_with_pil(ys, xs, yslice, xslice):
+ from PIL import Image
+ rxs = np.round(xs) - xslice[0]
+ rys = np.round(ys) - yslice[0]
+
+ lx = xslice[1] - xslice[0]
+ ly = yslice[1] - yslice[0]
+ rxys = zip(rxs, rys)
+
+ img = Image.new('L', (lx, ly), 0)
+ draw = Image.core.draw(img.im, 0)
+ ink = draw.draw_ink(1, "white")
+ draw.draw_polygon(rxys, ink, 1)
+ draw.draw_polygon(rxys, ink, 0)
+ return np.array(img) != 0
+
+
+def star_in_polygon((max_y, max_x), polar_coordinate_boundary, seed_x, seed_y, polar_transform):
+ polygon_x, polygon_y = polar_to_cartesian(polar_coordinate_boundary, seed_x, seed_y, polar_transform)
+
+ polygon_x_bounded = np.maximum(0, np.minimum(max_x - 1, polygon_x))
+ polygon_y_bounded = np.maximum(0, np.minimum(max_y - 1, polygon_y))
+
+ x1 = int(math.floor(np.min(polygon_x_bounded)))
+ x2 = int(math.ceil(np.max(polygon_x_bounded)) + 1)
+ y1 = int(math.floor(np.min(polygon_y_bounded)))
+ y2 = int(math.ceil(np.max(polygon_y_bounded)) + 1)
+
+ small_boolean_mask = mask_with_pil(polygon_y_bounded, polygon_x_bounded, (y1, y2), (x1, x2))
+
+ boolean_mask = np.zeros((max_y, max_x), dtype=bool)
+ boolean_mask[y1:y2, x1:x2] = small_boolean_mask
+
+ yx = [y1, x1]
+
+ return boolean_mask, small_boolean_mask, yx
+
+
+def multiply_list(ls, times):
+ list_length = len(ls)
+ integer_times = int(times)
+ fraction_elements = int((times - int(times)) * list_length)
+ res = ls * integer_times
+ res += ls[:fraction_elements]
+ return res
+
+
+def to_int(num):
+ return int(num)
+
+
+def fast_power(a, n):
+ mn = a
+ res = 1
+ n = int(n)
+ while n > 0:
+ if (n % 2 == 1):
+ res *= mn
+ mn = mn * mn
+ n /= 2
+ return res
\ No newline at end of file
diff --git a/cellstar/utils/debug_util.py b/cellstar/utils/debug_util.py
new file mode 100644
index 00000000..e2a81b09
--- /dev/null
+++ b/cellstar/utils/debug_util.py
@@ -0,0 +1,260 @@
+# -*- coding: utf-8 -*-
+"""
+Utilities with tools that can help with debuging / profiling CellStar
+Date: 2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import os
+from os import makedirs
+from os.path import exists
+
+import numpy as np
+import scipy as sp
+
+from cellstar.core import image_repo
+
+debug_image_path = "debug"
+
+# profiling switches
+PROFILE_SPEED = False
+PROFILE_MEMORY = False
+
+# main switch which turn off all debugging utils (always deploy with False)
+DEBUGING = False
+SHOW = False
+
+SNAKE_PROPERTIES = False
+# allow the user to inspect cell star results before segmentation (only for debugging)
+EXPLORE = False
+
+# pyplot import can fail if cellstar used as a plugin
+try:
+ import matplotlib
+ import matplotlib.pyplot as plt
+
+ if not SHOW:
+ try:
+ matplotlib.use('Agg')
+ except:
+ pass
+except:
+ DEBUGING = False # turn off debugging images if unavailable (e.g. in CP 2.2 BETA)
+
+# try to load user32.dll for key lock state
+user32_dll = None
+try:
+ import ctypes
+ user32_dll = ctypes.WinDLL ("User32.dll")
+ user32_dll.GetKeyState.restype = ctypes.c_short
+except:
+ pass
+
+
+def check_caps_scroll_state():
+ if user32_dll is None:
+ return None, None
+
+ VK_CAPITAL = 0x14
+ VK_SCROLL = 0X91
+ return user32_dll.GetKeyState(VK_CAPITAL), user32_dll.GetKeyState(VK_SCROLL)
+
+
+def prepare_debug_folder():
+ if not exists(debug_image_path):
+ makedirs(debug_image_path)
+
+
+def draw_seeds_on_axes(seeds, axes):
+ if DEBUGING:
+ return axes.plot([s.x for s in seeds], [s.y for s in seeds], 'bo', markersize=3)
+
+
+def draw_seeds(seeds, background, title="some_source"):
+ if DEBUGING:
+ prepare_debug_folder()
+ fig = plt.figure("draw_seeds")
+ fig.frameon = False
+ plt.imshow(background, cmap=plt.cm.gray)
+ plt.plot([s.x for s in seeds], [s.y for s in seeds], 'bo', markersize=3)
+ plt.savefig(os.path.join(debug_image_path, "seeds_" + title + ".png"), pad_inches=0.0)
+ fig.clf()
+ plt.close(fig)
+
+
+def images_repo_save(images):
+ """
+ @type images: image_repo.ImageRepo
+ """
+ image_save(images.background, "background")
+ image_save(images.brighter, "brighter")
+ image_save(images.brighter_original, "brighter_original")
+ image_save(images.darker, "darker")
+ image_save(images.darker_original, "darker_original")
+ image_save(images.cell_content_mask, "cell_content_mask")
+ image_save(images.cell_border_mask, "cell_border_mask")
+ image_save(images.foreground_mask, "foreground_mask")
+ image_save(images.mask, "image_mask")
+ image_save(images.image_back_difference, "image_back_difference")
+
+
+def image_save(image, title):
+ """
+ Displays image with title using matplotlib.pyplot
+ @param image:
+ @param title:
+ """
+ if DEBUGING:
+ prepare_debug_folder()
+ file_path = os.path.join(debug_image_path, title + '.tif')
+ sp.misc.imsave(file_path, image)
+ return file_path
+ return None
+
+
+def image_show(image, title, override=False):
+ """
+ Displays image with title using matplotlib.pyplot
+ @param image:
+ @param title:
+ """
+ if DEBUGING and (SHOW or override):
+ prepare_debug_folder()
+ fig = plt.figure(title)
+ plt.imshow(image, cmap=plt.cm.gray, interpolation='none')
+ plt.show()
+ fig.clf()
+ plt.close(fig)
+
+
+def draw_overlay(image, x, y):
+ if DEBUGING and SHOW:
+ prepare_debug_folder()
+ fig = plt.figure()
+ plt.imshow(image, cmap=plt.cm.gray, interpolation='none')
+ plt.plot(x, y)
+ plt.show()
+ fig.clf()
+ plt.close(fig)
+
+
+def explorer_expected():
+ if DEBUGING:
+ check_state = check_caps_scroll_state()
+ if EXPLORE or check_state[0] == True and check_state[1] == True:
+ return True
+ return False
+
+
+def explore_cellstar(cellstar=None, seeds=[], snakes=[], images=None, image=None, params=None):
+ if explorer_expected():
+ value = 0
+ try:
+ app = None
+ try:
+ import wx
+ app = wx.App(0)
+ except:
+ pass
+
+ import utils.explorer as exp
+ if image is None:
+ image = cellstar.images.image
+ if images is None:
+ images = cellstar.images
+
+ explorer_ui = exp.ExplorerFrame(images)
+ explorer = exp.Explorer(image, images, explorer_ui, cellstar, params)
+ explorer.stick_seeds = seeds
+ explorer.stick_snakes = snakes
+ value = explorer_ui.ShowModal()
+
+ #if app is not None:
+ # app.MainLoop()
+ except Exception as ex:
+ print ex
+ pass
+
+ if value == exp.ExplorerFrame.ABORTED:
+ raise Exception("Execution aborted")
+
+
+def draw_snakes_on_axes(snakes, axes, outliers=.1):
+ if DEBUGING and len(snakes) >= 1:
+ snakes = sorted(snakes, key=lambda ss: ss.rank)
+ snakes_tc = snakes[:max(1, int(len(snakes) * (1 - outliers)))]
+
+ max_rank = snakes_tc[-1].rank
+ min_rank = snakes_tc[0].rank
+
+ rank_range = max_rank - min_rank
+ if rank_range == 0: # for example there is one snake
+ rank_range = max_rank
+
+ rank_ci = lambda rank: 999 * ((rank - min_rank) / rank_range) if rank <= max_rank else 999
+ colors = plt.cm.jet(np.linspace(0, 1, 1000))
+ s_colors = [colors[int(rank_ci(s.rank))] for s in snakes]
+
+ # we want the best on top
+ for snake, color in reversed(zip(snakes, s_colors)):
+ axes.plot(snake.xs, snake.ys, c=color, linewidth=1.0)
+
+
+def draw_snakes(image, snakes, outliers=.1, it=0):
+ if DEBUGING and len(snakes) > 1:
+ prepare_debug_folder()
+
+ fig = plt.figure("draw_snakes")
+ plt.imshow(image, cmap=plt.cm.gray, interpolation='none')
+ draw_snakes_on_axes(snakes, plt)
+
+ plt.savefig(os.path.join(debug_image_path, "snakes_rainbow_" + str(it) + ".png"), pad_inches=0.0)
+ if SHOW:
+ plt.show()
+
+ fig.clf()
+ plt.close(fig)
+
+try:
+ import line_profiler
+ import memory_profiler
+except:
+ pass
+
+
+def speed_profile(func):
+ def profiled_func(*args, **kwargs):
+ try:
+ profiler = line_profiler.LineProfiler()
+ profiler.add_function(func)
+ profiler.enable_by_count()
+ return func(*args, **kwargs)
+ finally:
+ profiler.print_stats()
+
+ if PROFILE_SPEED:
+ return profiled_func
+ else:
+ return func
+
+
+def memory_profile(func):
+ if not PROFILE_MEMORY:
+ return func
+ else:
+ if func is not None:
+ def wrapper(*args, **kwargs):
+ if PROFILE_MEMORY:
+ prof = memory_profiler.LineProfiler()
+ val = prof(func)(*args, **kwargs)
+ memory_profiler.show_results(prof)
+ else:
+ val = func(*args, **kwargs)
+ return val
+
+ return wrapper
+ else:
+ def inner_wrapper(f):
+ return memory_profiler.profile(f)
+
+ return inner_wrapper
diff --git a/cellstar/utils/image_util.py b/cellstar/utils/image_util.py
new file mode 100644
index 00000000..73e31eaf
--- /dev/null
+++ b/cellstar/utils/image_util.py
@@ -0,0 +1,323 @@
+# -*- coding: utf-8 -*-
+"""
+Image util module contains additional methods for easy image manipulations.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import numpy as np
+import scipy as sp
+import scipy.misc
+import scipy.ndimage
+from numpy import argwhere
+from numpy.fft import rfft2, irfft2
+from scipy.ndimage.filters import *
+
+try:
+ # for older version of scipy
+ from scipy.signal.signaltools import _next_regular as next_fast_len
+except:
+ # for 0.19 version of scipy
+ from scipy.fftpack.helper import next_fast_len
+
+
+from calc_util import extend_slices, fast_power, to_int
+
+
+def fft_convolve(in1, in2, times):
+ def _centered(arr, newsize):
+ # Return the center newsize portion of the array.
+ currsize = np.array(arr.shape)
+ startind = (currsize - newsize) // 2
+ endind = startind + newsize
+ myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
+ return arr[tuple(myslice)]
+
+ if times == 0:
+ return in1.copy()
+
+
+ s1 = np.array(in1.shape)
+ s2 = np.array(in2.shape)
+ shape = s1 + (s2 - 1) * times
+
+ # Speed up FFT by padding to optimal size for FFTPACK
+ fshape = [next_fast_len(int(d)) for d in shape]
+ fslice = tuple([slice(0, int(sz)) for sz in shape])
+
+ resfft = fast_power(rfft2(in2, fshape), times)
+ resfft = resfft * rfft2(in1, fshape)
+ ret = irfft2(resfft, fshape)[fslice].copy()
+ ret = ret.real
+
+ return _centered(ret, s1)
+
+
+def get_bounding_box(image_mask):
+ """
+ Calculates the minimal bounding box for non zero elements.
+ @returns [ystart,ystop), [xstart,xstop) or None, None
+ """
+ non_zero_points = argwhere(image_mask)
+ if len(non_zero_points) == 0:
+ return None
+ (ystart, xstart), (ystop, xstop) = non_zero_points.min(0), non_zero_points.max(0) + 1
+ return (ystart, ystop), (xstart, xstop)
+
+
+def get_circle_kernel(radius):
+ """
+ Creates radius x radius bool image of the circle.
+ @param radius: radius of the circle
+ """
+ y, x = np.ogrid[np.floor(-radius):np.ceil(radius) + 1, np.floor(-radius):np.ceil(radius) + 1]
+ return x ** 2 + y ** 2 <= radius ** 2
+
+
+def image_dilate(image, radius):
+ image = np.copy(image)
+ if radius <= 1:
+ return image
+
+ box = get_bounding_box(image)
+ if box is None:
+ return image
+ ys, xs = box
+ lp, hp = contain_pixel(image.shape, (ys[0] - radius, xs[0] - radius)), \
+ contain_pixel(image.shape, (ys[1] + radius, xs[1] + radius))
+ ys, xs = (lp[0], hp[0]), (lp[1], hp[1])
+ morphology_element = get_circle_kernel(radius)
+ dilated_part = sp.ndimage.morphology.binary_dilation(image[ys[0]:ys[1], xs[0]:xs[1]], morphology_element)
+ image[ys[0]:ys[1], xs[0]:xs[1]] = dilated_part
+ return image
+
+
+def image_dilate_with_element(image, n):
+ return sp.ndimage.morphology.grey_dilation(image, size=(n, n))
+
+
+def image_erode(image, radius):
+ morphology_element = get_circle_kernel(radius)
+ return sp.ndimage.morphology.binary_erosion(image, morphology_element)
+
+
+def fill_foreground_holes(mask, kernel_size, minimal_hole_size, min_cluster_area_scaled, mask_min_radius_scaled):
+ filled_black_holes = fill_holes(mask, kernel_size, minimal_hole_size)
+
+ holes_remaining = np.logical_not(filled_black_holes)
+ filled_small_holes = mark_small_areas(holes_remaining, min_cluster_area_scaled, filled_black_holes)
+
+ morphology_enhanced = image_erode(filled_small_holes, mask_min_radius_scaled)
+ morphology_enhanced = image_dilate(morphology_enhanced, mask_min_radius_scaled)
+
+ dilated_mask = dilate_big_areas(morphology_enhanced, min_cluster_area_scaled, kernel_size)
+
+ return dilated_mask
+
+
+def mark_small_areas(mask, max_hole_size, result_mask):
+ components, num_components = sp.ndimage.label(mask, np.ones((3, 3)))
+ slices = sp.ndimage.find_objects(components)
+ for label, slice in zip(range(1, num_components + 1), slices):
+ components_slice = components[slice] == label
+ if np.count_nonzero(components_slice) < max_hole_size:
+ result_mask[slice][components_slice] = True
+ return result_mask
+
+
+def dilate_big_areas(mask, min_area_size, dilate_radius):
+ components, num_components = sp.ndimage.label(mask, np.ones((3, 3)))
+ component = np.zeros(mask.shape, dtype=bool)
+ for label in range(1, num_components + 1):
+ np.equal(components, label, component)
+ if np.count_nonzero(component) > min_area_size:
+ tmp_mask = image_dilate(component, dilate_radius)
+ mask = mask | tmp_mask
+
+ return mask
+
+
+def fill_holes(mask, kernel_size, minimal_hole_size):
+ """
+ Fills holes in a given mask using iterative close + dilate morphological operations and filtering small patches.
+ @param mask: mask which holes are to be filled
+ @param kernel_size: size of the morphological element used to dilate/erode mask
+ @param minimal_hole_size: holes with area smaller than param are to be removed
+ """
+ nr = 1
+ morphology_element = get_circle_kernel(kernel_size)
+ while True:
+ new_mask = mask.copy()
+ # find connected components
+ components, num_components = sp.ndimage.label(np.logical_not(new_mask), np.ones((3, 3)))
+ slices = sp.ndimage.find_objects(components)
+ for label, slice in zip(range(1, num_components + 1), slices):
+ slice = extend_slices(slice, to_int(kernel_size * 2))
+ components_slice = components[slice] == label
+ # filter small components
+ if np.count_nonzero(components_slice) < minimal_hole_size:
+ new_mask[slice] |= components_slice
+ else:
+ # shrink components and check if they fell apart
+ # close holes
+ components_slice = sp.ndimage.morphology.binary_closing(components_slice, morphology_element)
+
+ # erode holes
+ components_slice = sp.ndimage.morphology.binary_erosion(components_slice, morphology_element)
+
+ # don't invade masked pixels
+ components_slice &= np.logical_not(new_mask[slice])
+
+ # recount connected components and check sizes
+ mark_small_areas(components_slice, minimal_hole_size, new_mask[slice])
+
+ # check if it is the fix point
+ if (mask == new_mask).all():
+ break
+ else:
+ mask = new_mask
+
+ nr += 1
+
+ return mask
+
+
+def contain_pixel(shape, pixelYX):
+ """
+ Trims pixel to given dimentions, converts pixel position to int
+ @param shape: size (height, width) exclusive
+ @param pixel: pixel to push inside shape
+ """
+ (py, px) = pixelYX
+ (py, px) = ((np.minimum(np.maximum(py + 0.5, 0), shape[0] - 1)).astype(int),
+ (np.minimum(np.maximum(px + 0.5, 0), shape[1] - 1)).astype(int))
+ return py, px
+
+
+def find_maxima(image):
+ """
+ Finds local maxima in given image
+ @param image: image for which maxima will be found
+ """
+ height = image.shape[0]
+ width = image.shape[1]
+
+ right = np.zeros(image.shape, dtype=bool)
+ left = np.zeros(image.shape, dtype=bool)
+ up = np.zeros(image.shape, dtype=bool)
+ down = np.zeros(image.shape, dtype=bool)
+
+ right[:, 0:width - 1] = image[:, 0:width - 1] > image[:, 1:width]
+ left[:, 1:width] = image[:, 1:width] > image[:, 0:width - 1]
+ up[0:height - 1, :] = image[0:height - 1, :] > image[1:height, :]
+ down[1:height, :] = image[1:height, :] > image[0:height - 1, :]
+
+ return right & left & up & down
+
+
+def exclude_segments(image, segments, val):
+ """
+ Sets exclusion value for given segments in given image
+ @param image: image from which segments will be excluded
+ @param segments: segments to be excluded from image
+ @param val: value to be set in segments as exclusion value
+ """
+ segment_mask = segments > 0
+ inverted_segment_mask = np.logical_not(segment_mask)
+ image_segments_zeroed = image * inverted_segment_mask
+ image_segments_valued = image_segments_zeroed + (segment_mask * val)
+
+ return image_segments_valued
+
+
+def image_median_filter(image, size):
+ if size < 1:
+ return image
+
+ size = to_int(size)
+ return median_filter(image, (size, size))
+
+
+def image_blur(image, times):
+ """
+ Performs image blur with kernel: [[2, 3, 2], [3, 12, 3], [2, 3, 2]] / 32
+ @param image: image to be blurred (assumed as numpy.array of values from 0 to 1)
+ @param times: specifies how many times blurring will be performed
+ """
+ kernel = np.array([[2, 3, 2], [3, 12, 3], [2, 3, 2]]) / 32.0
+
+ if times >= 8:
+ return fft_convolve(image, kernel, times)
+ else:
+ blurred = convolve(image, kernel)
+ for _ in xrange(int(times) - 1):
+ blurred = convolve(blurred, kernel)
+ return blurred
+
+
+def image_smooth(image, radius, fft_use=True):
+ """
+ Performs image blur with circular kernel.
+ @param image: image to be blurred (assumed as numpy.array of values from 0 to 1)
+ @param radius: radius of the kernel
+ """
+ if radius < 1:
+ return image
+
+ kernel = get_circle_kernel(radius).astype(float)
+ kernel /= np.sum(kernel)
+ image = np.array(image, dtype=float)
+
+ if radius >= 8 and fft_use:
+ image_2 = np.pad(image, int(radius), mode='reflect')
+ res = fft_convolve(image_2, kernel, 1)
+ radius_round = to_int(radius)
+ return res[radius_round:-radius_round, radius_round:-radius_round]
+ else:
+ return convolve(image, kernel, mode='reflect', cval=0.0)
+
+
+def image_normalize(image):
+ """
+ Performs image normalization (vide: matlab mat2gray)
+ @param image: image to be normalized (assumed as numpy.array of values from 0 to 1)
+ """
+ minimum = np.amin(image)
+ maximum = np.amax(image)
+
+ delta = 1
+ if maximum != minimum:
+ delta = 1 / (maximum - minimum)
+ shift = - minimum * delta
+
+ image_normalized = delta * image + shift
+
+ return np.minimum(np.maximum(image_normalized, 0), 1)
+
+
+def set_image_border(image, val):
+ """
+ Sets pixel values at image borders to given value
+ @param image: image that borders will be set to given value
+ @param val: value to be s et
+ """
+ image[0, :] = val
+ image[:, 0] = val
+ image[image.shape[0] - 1, :] = val
+ image[:, image.shape[1] - 1] = val
+
+ return image
+
+
+def load_image(filename, scaling=True):
+ if filename == '':
+ return None
+ image = scipy.misc.imread(filename)
+ if image.max() > 1 and scaling:
+ image2d = np.array(image, dtype=float) / np.iinfo(image.dtype).max
+ else:
+ image2d = image.astype(float)
+
+ if image2d.ndim == 3:
+ image2d = np.sum(image, 2) / image.shape[2]
+ return image2d
diff --git a/cellstar/utils/index.py b/cellstar/utils/index.py
new file mode 100644
index 00000000..7d50aea9
--- /dev/null
+++ b/cellstar/utils/index.py
@@ -0,0 +1,23 @@
+# -*- coding: utf-8 -*-
+"""
+Index is a convenient structure for processing calculate all angle x radius point in one numpy operation.
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import numpy as np
+
+
+class Index(object):
+ @classmethod
+ def create(cls, px, py):
+ return np.column_stack((py.flat, px.flat)).astype(np.int64)
+
+ @staticmethod
+ def to_numpy(index):
+ if len(index.shape) == 2:
+ return index[:, 0], index[:, 1]
+ elif len(index.shape) == 3:
+ return index[:, :, 0], index[:, :, 1]
+ else:
+ return index
diff --git a/cellstar/utils/params_util.py b/cellstar/utils/params_util.py
new file mode 100644
index 00000000..09df6953
--- /dev/null
+++ b/cellstar/utils/params_util.py
@@ -0,0 +1,92 @@
+# -*- coding: utf-8 -*-
+"""
+Params util module contains methods for manipulating parameters and precision to parameters mapping
+Date: 2013-2016
+Website: http://cellstar-algorithm.org/
+"""
+
+import numpy as np
+
+from cellstar.core.config import default_config
+
+
+def default_parameters(segmentation_precision=-1, avg_cell_diameter=-1):
+ parameters = default_config()
+ if avg_cell_diameter != -1:
+ parameters["segmentation"]["avgCellDiameter"] = avg_cell_diameter
+
+ if segmentation_precision is None:
+ return parameters
+ else:
+ return parameters_from_segmentation_precision(parameters, segmentation_precision)
+
+
+def create_size_weights(size_weight_average, length):
+ if length == 1:
+ size_weight_multiplier = np.array([1])
+ elif length == 2:
+ size_weight_multiplier = np.array([0.8, 1.25])
+ elif length == 3:
+ size_weight_multiplier = np.array([0.6, 1, 1.6])
+ elif length == 4:
+ size_weight_multiplier = np.array([0.5, 0.8, 1.3, 2])
+ elif length == 5:
+ size_weight_multiplier = np.array([0.5, 0.8, 1, 1.3, 2])
+ elif length == 6:
+ size_weight_multiplier = np.array([0.35, 0.5, 0.8, 1.3, 2, 3])
+ else:
+ size_weight_multiplier = np.array([0.25, 0.35, 0.5, 0.8, 1.3, 2, 3, 5, 8])
+
+ return size_weight_average * size_weight_multiplier / np.average(size_weight_multiplier)
+
+
+def parameters_from_segmentation_precision(parameters, segmentation_precision):
+ sfrom = lambda x: max(0, segmentation_precision - x)
+ segmentation_precision = min(20, segmentation_precision)
+ if segmentation_precision <= 0:
+ parameters["segmentation"]["steps"] = 0
+ elif segmentation_precision <= 6:
+ parameters["segmentation"]["steps"] = 1
+ else:
+ parameters["segmentation"]["steps"] = min(10, segmentation_precision - 5)
+
+ parameters["segmentation"]["stars"]["points"] = 8 + max(segmentation_precision - 2, 0) * 4
+
+ parameters["segmentation"]["maxFreeBorder"] = \
+ max(0.4, 0.7 * 16 / max(16, parameters["segmentation"]["stars"]["points"]))
+
+ parameters["segmentation"]["seeding"]["from"]["cellBorder"] = int(segmentation_precision >= 2)
+ parameters["segmentation"]["seeding"]["from"]["cellBorderRandom"] = sfrom(14)
+ parameters["segmentation"]["seeding"]["from"]["cellContent"] = int(segmentation_precision >= 11)
+ parameters["segmentation"]["seeding"]["from"]["cellContentRandom"] = min(4, sfrom(12))
+ parameters["segmentation"]["seeding"]["from"]["cellBorderRemovingCurrSegments"] = \
+ int(segmentation_precision >= 11)
+ parameters["segmentation"]["seeding"]["from"]["cellBorderRemovingCurrSegmentsRandom"] = max(0, min(4, sfrom(16)))
+ parameters["segmentation"]["seeding"]["from"]["cellContentRemovingCurrSegments"] = \
+ int(segmentation_precision >= 7)
+ parameters["segmentation"]["seeding"]["from"]["cellContentRemovingCurrSegmentsRandom"] = max(0, min(4, sfrom(12)))
+ parameters["segmentation"]["seeding"]["from"]["snakesCentroids"] = int(segmentation_precision >= 9)
+ parameters["segmentation"]["seeding"]["from"]["snakesCentroidsRandom"] = max(0, min(4, sfrom(14)))
+
+ parameters["segmentation"]["stars"]["step"] = 0.0067 * max(1, (1 + (15 - segmentation_precision) / 2.0))
+
+ if segmentation_precision <= 9:
+ weight_length = 1
+ elif segmentation_precision <= 11:
+ weight_length = 2
+ elif segmentation_precision <= 13:
+ weight_length = 3
+ elif segmentation_precision <= 15:
+ weight_length = 4
+ elif segmentation_precision <= 17:
+ weight_length = 6
+ else:
+ weight_length = 9
+
+ parameters["segmentation"]["stars"]["sizeWeight"] = list(
+ create_size_weights(np.average(parameters["segmentation"]["stars"]["sizeWeight"]), weight_length)
+ )
+
+ parameters["segmentation"]["foreground"]["pickyDetection"] = segmentation_precision > 8
+
+ return parameters
diff --git a/identifyyeastcells.py b/identifyyeastcells.py
new file mode 100644
index 00000000..72c24397
--- /dev/null
+++ b/identifyyeastcells.py
@@ -0,0 +1,1148 @@
+#!/usr/bin/env python
+import threading
+
+import cellprofiler.icons
+try:
+ # version 3.0.0
+ from cellprofiler.modules._help import PROTIP_RECOMEND_ICON, PROTIP_AVOID_ICON, TECH_NOTE_ICON
+except:
+ from cellprofiler.modules._help import PROTIP_RECOMMEND_ICON as PROTIP_RECOMEND_ICON
+ from cellprofiler.modules._help import PROTIP_AVOID_ICON, TECH_NOTE_ICON
+
+__doc__ = """IdentifyYeastCells identifies yeast (or other round) objects in the image. This module was designed
+to work on brightfield images. However in some cases it can also be efficient with fluorescent imagery.
+
+The important part of the module is parameter fitting mechanism which allows to adapt it to many different types of imagery.
+
+There is a number of various usage example available on website
+CellStar algorithm which includes ready to use pipelines.
+
+
+
+What do I need as input?
+To use this module, you will need to make sure that your input image has the following qualities:
+
+- The image should be grayscale.
+
+
+What do the settings mean?
+See below for help on the individual settings. The following icons are used to call attention to
+key items:
+
+
Our recommendation or example use case
+for which a particular setting is best used.
Indicates a condition under which
+a particular setting may not work well.
Technical note. Provides more
+detailed information on the setting.
+
+What do I get as output?
+A set of primary objects are produced by this module, which can be used in downstream modules
+for measurement purposes or other operations.
+See the section "Available measurements" below for
+the measurements that are produced by this module.
+
+Once the module has finished processing, the module display window
+will show the following panels:
+
+- Left: The raw, original image.
+- Right: The identified objects shown as a color
+image where connected pixels that belong to the same object are assigned the same
+color (label image). It is important to note that assigned colors are
+arbitrary; they are used simply to help you distingush the various objects.
+
+
+
+Available measurements
+Image measurements:
+
+- Count: The number of primary objects identified.
+
+
+Object measurements:
+
+- Location_X, Location_Y: The pixel (X,Y) coordinates of the primary
+object centroids. The centroid is calculated as the center of mass of the binary
+representation of the object.
+- Features_Quality: The module internal quality measure describing how well objects fits the model.
+It can be used to filter out the weakest objects.
+
+
+Publications
+
+This module was prepared by Filip Mroz, Adam Kaczmarek and Szymon Stoma. The method uses CellStar approach
+developed by Versari et al. It is also available with a dedicated GUI for semi-automatic segmentation.
+
+"Long-term Tracking of Budding Yeast Cells in Brightfield Microscopy: CellStar and the Evaluation Platform."
+Journal of The Royal Society Interface 14.127 (2017)
+
+Please reach us at CellStar algorithm website for inquires and manuals.
+
+Credits
+Filip Mroz, Adam Kaczmarek, Szymon Stoma,
+Artemis Llamosi,
+Kirill Batmanov,
+Cristian Versari,
+Cedric Lhoussaine,
+Gabor Csucs,
+Simon F. Noerrelykke,
+Gregory Batt,
+Pawel Rychlikowski,
+Pascal Hersen.
+
+""".format(PROTIP_RECOMEND_ICON=PROTIP_RECOMEND_ICON, PROTIP_AVOID_ICON=PROTIP_AVOID_ICON,
+ TECH_NOTE_ICON=TECH_NOTE_ICON)
+
+# Module documentation variables:
+__authors__ = """Filip Mroz,
+Adam Kaczmarek,
+Szymon Stoma
+"""
+__contact__ = "fafafft@gmail.com"
+__date__ = "2017.05.24"
+__version__ = "1.3.0"
+__docformat__ = "restructuredtext en"
+__revision__ = "$Id$"
+
+#################################
+#
+# Imports from useful Python libraries
+#
+#################################
+from os.path import join as pj
+from os.path import isfile
+import logging
+
+logger = logging.getLogger(__name__)
+
+import math
+import numpy as np
+import scipy.ndimage
+
+#################################
+#
+# Imports from CellProfiler
+#
+##################################
+
+import bioformats
+
+import cellprofiler.modules
+import cellprofiler.measurement as cpmeas
+import cellprofiler.setting as cps
+import cellprofiler.preferences as pref
+
+#################################
+#
+# Specific imports
+#
+##################################
+
+from cellstar.utils.params_util import default_parameters, create_size_weights
+from cellstar.segmentation import Segmentation
+from cellstar.parameter_fitting.pf_runner import run_pf, run_rank_pf
+from cellstar.utils.debug_util import memory_profile, explorer_expected
+
+###################################
+#
+# Constants
+#
+###################################
+
+DEBUG = 1
+F_BACKGROUND = "Background"
+BKG_CURRENT = 'Actual image'
+BKG_FIRST = 'First image'
+BKG_FILE = 'File'
+
+C_OBJECT_FEATURES = "Features"
+
+FTR_OBJECT_QUALITY = "Quality"
+'''The object quality - floating number the higher the better'''
+M_OBJECT_FEATURES_OBJECT_QUALITY = '%s_%s' % (C_OBJECT_FEATURES, FTR_OBJECT_QUALITY)
+
+
+def hack_add_from_file_into_EditObjects(dialog_box):
+ import wx
+ def on_load(event):
+ with wx.FileDialog(None,
+ message="Select image with labels",
+ wildcard="Image file (*.tif,*.tiff,*.jpg,*.jpeg,*.png,*.gif,*.bmp)|*.tif;*.tiff;*.jpg;*.jpeg;*.png;*.gif;*.bmp|*.* (all files)|*.*",
+ style=wx.FD_OPEN) as dlg:
+ if dlg.ShowModal() != wx.ID_OK:
+ return
+ path = dlg.Path
+
+ labels_loaded = (bioformats.load_image(path) * 255).astype(int)
+ if labels_loaded.ndim == 3:
+ labels_loaded = np.sum(labels_loaded, 2) / labels_loaded.shape[2]
+
+ for l in np.unique(labels_loaded):
+ if l != 0:
+ dialog_box.add_label(l == labels_loaded)
+ # no shape check, no float check
+ dialog_box.init_labels()
+ dialog_box.display()
+ dialog_box.record_undo()
+
+ ID_ACTION_LOAD_FROM_FILE = wx.NewId()
+ sizer = dialog_box.toolbar.ContainingSizer
+ load_file_button = wx.Button(dialog_box, ID_ACTION_LOAD_FROM_FILE, "Add from file")
+ list(sizer.Children)[-1].Sizer.Add(load_file_button, 0, wx.ALIGN_CENTER)
+ dialog_box.Bind(wx.EVT_BUTTON, on_load, load_file_button)
+
+
+###################################
+#
+# The module class
+#
+###################################
+
+
+class IdentifyYeastCells(cellprofiler.module.ImageSegmentation):
+ module_name = "IdentifyYeastCells"
+ category = "Yeast Toolbox"
+ variable_revision_number = 9
+ current_workspace = ''
+ fitting_image_set = None
+ param_fit_progress = 0
+ param_fit_progress_partial = 0
+
+ def create_settings(self):
+ super(IdentifyYeastCells, self).create_settings()
+
+ self.input_image_name = self.x_name
+ self.input_image_name.doc = "How do you call the images you want to use to identify objects?"
+
+ self.object_name = self.y_name
+ self.object_name.doc = "How do you want to call the objects identified by this module?"
+
+ self.background_image_name = cps.ImageNameSubscriber(
+ "Select the empty field image", doc="""
+ (Used only when you select "loaded from file" background calculation strategy)
+ How do you call the image you want to use as background
+ image (same image will be used for every image in the workflow)?
+ """)
+
+ self.ignore_mask_image_name = cps.ImageNameSubscriber(
+ "Select ignore mask image", doc="""
+ You can provide a ignore mark with regions in the image which are to be ignored by the algorithm in segmentation.
+ """, can_be_blank=True)
+
+ self.background_elimination_strategy = cps.Choice(
+ 'Select the background calculation mode',
+ [BKG_CURRENT, BKG_FILE], doc="""
+ You can choose from the following options:
+
+ - loaded from file: Use this option if your background does not change at all for all images in the series.
+ - computed from actual image: This is default option. The algorithm will try to automatically compute the background for
+ each individual image. In some specific cases it is better to use manually precomputed background loaded from file.
+
""")
+
+ self.average_cell_diameter = cps.Float(
+ "Average cell diameter in pixels",
+ 30.0, minval=10, doc='''\
+ The average cell diameter is used to scale many algorithm parameters.
+ Please use e.g. ImageJ to measure the average cell size in pixels.
+ '''
+ )
+
+ self.advanced_cell_filtering = cps.Binary(
+ 'Do you want to filter cells by area?', False, doc="""
+ This parameter is used for cell filtering. The algorithm creates many more cell candidates than finally accepted cells (these
+ cells overlap with each other). At the final phase of algorithm the cells are selected from the ensemble of cell candidates based on the
+ "quality" measure. Use this option if you want to absolutely prohibit too small (or too big) cells to be chosen (regardless of "quality").
+ """)
+
+ self.min_cell_area = cps.Integer(
+ "Minimal area of accepted cell in pixels",
+ 900, minval=10, doc='''\
+ (Used only when you want to filter cells based on area)
+ The minimum cell area is used while final filtering of cells.
+ Please use e.g. ImageJ to measure the average cell size in pixels.
+ '''
+ )
+
+ self.max_cell_area = cps.Integer(
+ "Maximum area of accepted cell in pixels",
+ 5 * 900, minval=10, doc='''\
+ (Used only when you want to filter cells based on area)
+ The maximum cell area is used while final filtering of cells.
+ Please use e.g. ImageJ to measure the average cell size in pixels.
+ '''
+ )
+
+ self.background_brighter_then_cell_inside = cps.Binary(
+ 'Is the area without cells (background) brighter then cell interiors?', True, doc="""
+ Please check if the area inside of the cells is darker than the area without the cells (background). Use e.g. ImageJ to measure
+ average intensity.
+ """
+ )
+
+ self.bright_field_image = cps.Binary(
+ 'Do you want to segment brightfield images?', True, doc="""
+ Choose this option if you want to segment a brightfield image. For segmentation of fluorescent images please answer "No".
+ """
+ )
+
+ self.advanced_parameters = cps.Binary(
+ 'Use advanced configuration parameters?', False, doc="""
+ Do you want to use advanced parameters to configure plugin? They allow for more flexibility, and require
+ understanding of the algorithm to configure well.
+ """
+ )
+
+ self.segmentation_precision = cps.Integer(
+ "Segmentation precision",
+ 2, minval=1, maxval=5, doc='''\
+ (Used only when you want to specify advanced parameters)
+ Describes how thoroughly the algorithm searches for cells. Higher values should
+ make it easier to find smaller cells because the more parameters sets are searched.
+ The cost is longer runtime.
+ '''
+ )
+
+ self.specify_precision_details = cps.Binary(
+ "Do you want to edit details of segmentation precision?", False,
+ doc='''(Used only when you want to specify advanced parameters)
+ Do you want to edit details of segmentation precision?
+ These are low level parameters that can be modified
+ to further tweak result.
+ '''.format(TECH_NOTE_ICON=TECH_NOTE_ICON))
+
+ self.iterations = cps.Integer(
+ "Iterations",
+ 6, minval=1, maxval=15, doc='''\
+ (Used only when you want to specify precision details)
+ Number of segmentation iterations done by CellStar.''')
+
+ self.seeds_border = cps.Float(
+ "Seeds from border",
+ 1.2, minval=0, maxval=5, doc='''\
+ (Used only when you want to specify precision details)
+ How many seeds are to be extracted from this source:
+ val = 0 - no seeds from this source
+ val in (0.0, 0.5] - seeds only from second iteration
+ val in (0.5, 1) - seeds only in first iteration
+ val = 1 - seeds in every iteration
+ val in (1, 5.0) - additional random seeding: val*seeds_num total
+ ''')
+
+ self.seeds_content = cps.Float(
+ "Seeds from content",
+ 1.2, minval=0, maxval=5, doc=self.seeds_border.doc)
+
+ self.seeds_centroid = cps.Float(
+ "Seeds from centroid",
+ 1.2, minval=0, maxval=5, doc=self.seeds_border.doc)
+
+ self.contour_points = cps.Integer(
+ "Contour points",
+ 44, minval=36, maxval=70, doc='''\
+ (Used only when you want to specify precision details)
+ Number of points in every contour.''')
+
+ self.contour_precision = cps.Float(
+ "Contour precision",
+ 49.75, minval=25, maxval=200, doc='''\
+ (Used only when you want to specify precision details)
+ Number of space points covering average cell diameter.''')
+
+ self.weights_number = cps.Integer(
+ "Cell size variants",
+ 2, minval=1, maxval=4, doc='''\
+ (Used only when you want to specify precision details)
+ Number of tested difference size weights in every cell grow.''')
+
+ self.maximal_cell_overlap = cps.Float(
+ "Maximal overlap allowed while final filtering of cells",
+ 0.2, minval=0, maxval=1, doc='''\
+ (Used only when you want to specify advanced parameters)
+ This parameter is used for cell filtering. The algorithm creates many more cell candidates than finally accepted cells (these
+ cells overlap with each other). At the final phase of algorithm the cells are selected from the ensamble of cell candidates based on the
+ "quality" measure. One of the conditions checked while filtering is the overlap of the current candidate cell with already chosen cells.
+ Use this parameter if you want to allow to choose cells even if they overlap. Important: at the end cells do not overlap - they are
+ trimmed in such a way that the cell of higher "quality" will "borrow" the area of lower "quality" cells.
+ '''
+ )
+
+ self.autoadaptation_steps = cps.Integer(
+ "Number of steps in the autoadaptation procedure",
+ 1, minval=1, maxval=1000, doc='''
+ Describes how thoroughly we want to adapt the algorithm to current image sets. Higher values should
+ make it easier to correctly discover cells, however you will have to wait longer for the autoadaptation
+ procedure to finish. Remember that you do it once for all images, and you can copy the values from other
+ pipelines, if you have already found the parameters before.
+ '''
+ )
+
+ self.use_ground_truth_to_set_params = cps.DoSomething("", "Autoadapt parameters",
+ self.ground_truth_editor, doc="""
+ Use this option to autoadapt parameters required for correct contour identification. This procedure should be run once
+ for one image in the series. Using your input the algorithm "learns" to recognize cells. When you click this button the window will open.
+ Please select one of the images which you would like to segment. If you use background or ignore mask images you will have to provide them as well.
+ On the input image you should draw few cells (3-6) and click "Done". Then the
+ "learning" procedure will start. Please be patient. Usually single iteration lasts 1-3 min. The more iterations you will choose, the more likely
+ it is that the algorithm will work better on your images.
+
+ There is an experimental alternative to selecting existing image files which may be useful when images are "produced" by the pipeline:
+
+ - Start Test Run and step down to IdentifyYeastCells module
+ - Step through IdentifyYeastCells - it will run segmentation and remember input images.
+ - Now you can adapt without specifing any images.
+
+
+
+ - Recommendations:
+
+ - Start with single iteration, then check how well the algorithm is able to recognize your cells
+ - If you have more time, use more iterations. You can always cancel (it will last 1-5 min). The algorithm will
+ then remember the best parameters learned during the session.
+ - Parameters are stored as a text in "Autoadapted parameters". If you found good parameters for your datasets
+ you can copy them to another pipeline.<\li>
+
+
+ """.format(PROTIP_RECOMEND_ICON=PROTIP_RECOMEND_ICON))
+
+ self.show_autoadapted_params = cps.Binary(
+ 'Do you want to see autoadapted parameters?', False, doc="""
+ (Used only when you want to specify advanced parameters)
+ Use this option to display autoadapted parameters.""")
+
+ self.autoadapted_params = cps.Text(text="Autoadapted parameters: ",
+ value="[[0.0442, 304.45, 15.482, 189.40820000000002], [300, 10, 0, 18, 10]]",
+ doc="""
+ (Used only when you want to specify advanced and autoadapted parameters)
+ Autoadapted parameters are pasted here from the "learning" preocedure. These parameters are used to characterize cell borders.
+ Edit them only if you know what you are doing. If you found good parameters for your datasets
+ you can copy them to another pipeline.
+ """)
+
+ PRECISION_PARAMS_START = 18
+ PRECISION_PARAMS_END = 24
+
+ def settings(self):
+ return [self.input_image_name,
+ self.object_name,
+ self.average_cell_diameter,
+ self.segmentation_precision,
+ self.maximal_cell_overlap,
+ self.background_image_name,
+ self.advanced_parameters,
+ self.background_brighter_then_cell_inside,
+ self.bright_field_image,
+ self.min_cell_area,
+ self.max_cell_area,
+ self.advanced_cell_filtering,
+ self.background_elimination_strategy,
+ self.show_autoadapted_params,
+ self.autoadapted_params,
+ self.autoadaptation_steps,
+ self.ignore_mask_image_name,
+ self.specify_precision_details,
+ self.iterations,
+ self.seeds_border,
+ self.seeds_content,
+ self.seeds_centroid,
+ self.contour_points,
+ self.contour_precision,
+ self.weights_number,
+ ]
+
+ def visible_settings(self):
+ visible_settings = [self.input_image_name,
+ self.object_name,
+ self.average_cell_diameter,
+ ]
+
+ visible_settings += [self.bright_field_image,
+ self.background_brighter_then_cell_inside,
+ self.autoadaptation_steps,
+ self.use_ground_truth_to_set_params,
+ self.show_autoadapted_params]
+
+ if self.show_autoadapted_params:
+ visible_settings.append(self.autoadapted_params)
+
+ visible_settings.append(self.advanced_parameters)
+
+ #
+ # Show the user the background only if self.provide_background is checked
+ #
+ if self.advanced_parameters:
+ visible_settings.append(self.segmentation_precision)
+
+ visible_settings.append(self.specify_precision_details)
+ if self.specify_precision_details:
+ visible_settings.append(self.iterations)
+ visible_settings.append(self.seeds_border)
+ visible_settings.append(self.seeds_content)
+ visible_settings.append(self.seeds_centroid)
+ visible_settings.append(self.contour_points)
+ visible_settings.append(self.contour_precision)
+ visible_settings.append(self.weights_number)
+
+ visible_settings.append(self.maximal_cell_overlap)
+ visible_settings.append(self.advanced_cell_filtering)
+ if self.advanced_cell_filtering:
+ visible_settings.append(self.min_cell_area)
+ visible_settings.append(self.max_cell_area)
+ # Show the user the background only if self.provide_background is checked
+ visible_settings.append(self.background_elimination_strategy) #
+
+ if self.background_elimination_strategy == BKG_FILE:
+ visible_settings.append(self.background_image_name)
+
+ visible_settings.append(self.ignore_mask_image_name)
+
+ return visible_settings
+
+ def on_setting_changed(self, setting, pipeline):
+ '''If precision is changed then update all the related settings'''
+ if setting == self.segmentation_precision or \
+ setting == self.specify_precision_details and not self.specify_precision_details.value:
+ self.set_ui_from_precision(self.segmentation_precision.value)
+
+ def is_interactive(self):
+ return False
+
+ def prepare_group(self, workspace, grouping, image_numbers):
+ '''Erase module information at the start of a run'''
+ d = self.get_dictionary(workspace.image_set_list)
+ d.clear()
+ return True
+
+ def get_ws_dictionary(self, workspace):
+ return self.get_dictionary(workspace.image_set_list)
+
+ def __get(self, field, workspace, default):
+ if self.get_ws_dictionary(workspace).has_key(field):
+ return self.get_ws_dictionary(workspace)[field]
+ return default
+
+ def __set(self, field, workspace, value):
+ self.get_ws_dictionary(workspace)[field] = value
+
+ def upgrade_settings(self, setting_values, variable_revision_number, module_name, from_matlab):
+ '''Adjust setting values if they came from a previous revision
+
+ setting_values - a sequence of strings representing the settings
+ for the module as stored in the pipeline
+ variable_revision_number - the variable revision number of the
+ module at the time the pipeline was saved. Use this
+ to determine how the incoming setting values map
+ to those of the current module version.
+ module_name - the name of the module that did the saving. This can be
+ used to import the settings from another module if
+ that module was merged into the current module
+ from_matlab - True if the settings came from a Matlab pipeline, False
+ if the settings are from a CellProfiler 2.0 pipeline.
+
+ Overriding modules should return a tuple of setting_values,
+ variable_revision_number and True if upgraded to CP 2.0, otherwise
+ they should leave things as-is so that the caller can report
+ an error.
+ '''
+ if variable_revision_number == 2:
+ setting_values = setting_values + [False]
+ variable_revision_number = 3
+ if variable_revision_number == 3:
+ setting_values = setting_values + [True]
+ variable_revision_number = 4
+
+ if variable_revision_number < 7:
+ setting_values = setting_values + ['Leave blank'] # ignore mask
+ # decode precision index 3
+ setting_values[3] = str(self.precision_to_ui_map[int(setting_values[3])])
+ variable_revision_number = 7
+ if variable_revision_number == 7:
+ # fill new ones based on precision
+ setting_values = setting_values + [False]
+ params_from_precision = self.get_ui_params_from_precision(int(setting_values[3]))
+ setting_values[20:26 + 1] = params_from_precision
+ variable_revision_number = 8
+ if variable_revision_number == 8:
+ # two settings were removed
+ setting_values = setting_values[:6] + setting_values[8:]
+ variable_revision_number = 9
+ return setting_values, variable_revision_number, from_matlab
+
+ def display(self, workspace, figure=None):
+ if self.show_window:
+ figure.set_subplots((2, 1))
+
+ title = "Input image, cycle #%d" % (workspace.measurements.image_number,)
+ image = workspace.display_data.input_pixels
+ labeled_image = workspace.display_data.segmentation_pixels
+
+ ax = figure.subplot_imshow_grayscale(0, 0, image, title)
+ figure.subplot_imshow_labels(1, 0, labeled_image,
+ self.object_name.value,
+ sharexy=ax)
+
+ #
+ # Measuremets:
+ # - objects count
+ # - objects location
+
+ def get_measurement_columns(self, pipeline):
+ '''Column definitions for measurements made by IdentifyPrimAutomatic'''
+ columns = super(IdentifyYeastCells, self).get_measurement_columns(pipeline)
+ columns += [(self.object_name.value, M_OBJECT_FEATURES_OBJECT_QUALITY, cpmeas.COLTYPE_FLOAT)]
+
+ return columns
+
+ def get_categories(self, pipeline, object_name):
+ """Return the categories of measurements that this module produces
+
+ object_name - return measurements made on this object (or 'Image' for image measurements)
+ """
+ result = super(IdentifyYeastCells, self).get_categories(pipeline, object_name)
+ result += [C_OBJECT_FEATURES]
+ return result
+
+ def get_measurements(self, pipeline, object_name, category):
+ """Return the measurements that this module produces
+
+ object_name - return measurements made on this object (or 'Image' for image measurements)
+ category - return measurements made in this category
+ """
+ result = super(IdentifyYeastCells, self).get_measurements(pipeline, object_name, category)
+ if category == C_OBJECT_FEATURES:
+ result += [FTR_OBJECT_QUALITY]
+ return result
+
+ ui_to_precision_map = dict([(1, 9), (2, 11), (3, 12), (4, 13), (5, 15)])
+ precision_to_ui_map = dict([(v, k) for k, v in ui_to_precision_map.items()] + [(10, 2)])
+
+ @property
+ def decoded_segmentation_precision_value(self):
+ return self.ui_to_precision_map[self.segmentation_precision.value]
+
+ @memory_profile
+ def run(self, workspace):
+ input_image_name = self.input_image_name.value
+ self.current_workspace = workspace
+ image_set = workspace.image_set
+
+ # same image_set for fitting
+ self.fitting_image_set = image_set
+
+ #
+ # Load images from workspace
+ #
+ input_pixels, background_pixels, ignore_mask_pixels = self.get_input_images_from_image_set(image_set)
+
+ input_image = image_set.get_image(input_image_name, must_be_grayscale=True)
+ self.input_image_file_name = input_image.file_name
+ self.current_workspace.display_data.input_pixels = input_pixels
+
+ #
+ # Preprocessing
+ #
+ input_pixels, background_pixels, ignore_mask_pixels = self.preprocess_images(input_pixels, background_pixels,
+ ignore_mask_pixels)
+
+ #
+ # Segmentation
+ #
+ objects, objects_qualities, background_pixels = self.segmentation(input_pixels, background_pixels,
+ ignore_mask_pixels)
+ objects.parent_image = input_image
+
+ if self.__get(F_BACKGROUND, workspace, None) is None and self.background_elimination_strategy == BKG_FIRST:
+ self.__set(F_BACKGROUND, workspace, background_pixels)
+
+ workspace.object_set.add_objects(objects, self.object_name.value)
+
+ # Save measurements
+ workspace.measurements.add_measurement(self.object_name.value, M_OBJECT_FEATURES_OBJECT_QUALITY,
+ objects_qualities)
+ super(IdentifyYeastCells, self).add_measurements(workspace)
+
+ def set_params_from_ui(self, params):
+ def update_params(next_name, first_name, random_name, ui_value):
+ params["segmentation"]["seeding"]["from"][next_name] = ui_value >= 1 or ui_value <= 0.5
+ params["segmentation"]["seeding"]["from"][first_name] = ui_value > 0.5
+ params["segmentation"]["seeding"]["from"][random_name] = max(0, ui_value - 1)
+
+ params_seeding = params["segmentation"]["seeding"]["from"]
+ params["segmentation"]["steps"] = self.iterations.value
+ update_params("cellBorderRemovingCurrSegments", "cellBorder", "cellBorderRandom", self.seeds_border.value)
+ update_params("cellContentRemovingCurrSegments", "cellContent", "cellContentRandom", self.seeds_content.value)
+ params_seeding["cellBorderRemovingCurrSegmentsRandom"] = params_seeding["cellBorderRandom"]
+ params_seeding["cellContentRemovingCurrSegmentsRandom"] = params_seeding["cellContentRandom"]
+
+ params["segmentation"]["seeding"]["from"]["snakesCentroids"] = self.seeds_centroid.value > 0.0
+ params["segmentation"]["seeding"]["from"]["snakesCentroidsRandom"] = max(0, self.seeds_centroid.value - 1)
+
+ params["segmentation"]["stars"]["points"] = self.contour_points.value
+ params["segmentation"]["stars"]["step"] = 1.0 / self.contour_precision.value
+
+ default_size_weight_average = np.average(params["segmentation"]["stars"]["sizeWeight"])
+ params["segmentation"]["stars"]["sizeWeight"] = list(
+ create_size_weights(default_size_weight_average, self.weights_number.value)
+ )
+
+ def get_ui_params_from_precision(self, ui_precision):
+ def params_to_ui(params_seeding, next_name, first_name, random_name):
+ random = params_seeding[random_name]
+ next = params_seeding[next_name]
+ first = params_seeding[first_name]
+
+ if random != 0:
+ return random + 1
+ elif next and first:
+ return 1
+ elif first:
+ return 0.7
+ elif next:
+ return 0.5
+ return 0.0
+
+ def override_random_seeds_number(precision, value):
+ if precision == 2:
+ return max(value, 1.2)
+ elif precision >= 3:
+ return max(value, 1.5)
+ return value
+
+ params = default_parameters(self.ui_to_precision_map[ui_precision], 30)
+
+ ui_params = [params["segmentation"]["steps"]]
+ params_seeding = params["segmentation"]["seeding"]["from"]
+ ui_params.append(override_random_seeds_number(ui_precision,
+ params_to_ui(params_seeding, "cellBorderRemovingCurrSegments", "cellBorder", "cellBorderRandom")))
+ ui_params.append(override_random_seeds_number(ui_precision,
+ params_to_ui(params_seeding, "cellContentRemovingCurrSegments", "cellContent", "cellContentRandom")))
+ ui_params.append(override_random_seeds_number(ui_precision,
+ params_to_ui(params_seeding, "snakesCentroids", "snakesCentroids", "snakesCentroidsRandom")))
+
+ ui_params.append(params["segmentation"]["stars"]["points"])
+ ui_params.append(1.0 / params["segmentation"]["stars"]["step"])
+ ui_params.append(len(params["segmentation"]["stars"]["sizeWeight"]))
+
+ return ui_params
+
+ def set_ui_from_precision(self, ui_precision):
+ ui_precision_values = self.get_ui_params_from_precision(ui_precision)
+ ui_precision_settings = self.settings()[self.PRECISION_PARAMS_START:self.PRECISION_PARAMS_END + 1]
+ for setting, value in zip(ui_precision_settings, ui_precision_values):
+ setting.value = value
+
+ def prepare_cell_star_object(self, segmentation_precision):
+ cellstar = Segmentation(segmentation_precision, self.average_cell_diameter.value)
+ cellstar.parameters["segmentation"]["maxOverlap"] = self.maximal_cell_overlap.value
+ self.set_params_from_ui(cellstar.parameters)
+
+ if self.advanced_cell_filtering.value:
+ def calculate_area_multiplier(area):
+ return 4.0 * area / self.average_cell_diameter.value ** 2 / math.pi
+
+ # def calculate_size_multiplier(area):
+ # return calculate_area_multiplier(area) ** 0.5
+
+ areas_range = self.min_cell_area.value, self.max_cell_area.value
+ cellstar.parameters["segmentation"]["minArea"] = max(cellstar.parameters["segmentation"]["minArea"],
+ calculate_area_multiplier(areas_range[0]))
+ cellstar.parameters["segmentation"]["maxArea"] = calculate_area_multiplier(areas_range[1])
+ # to some extent change length of rays
+ # cellstar.parameters["segmentation"]["stars"]["maxSize"] = max(cellstar.parameters["segmentation"]["stars"]["maxSize"], min(2.5, calculate_size_multiplier(areas_range[1])))
+
+ success = cellstar.decode_auto_params(self.autoadapted_params.value)
+ if not success: # if current value is invalid overwrite it with current settings
+ self.autoadapted_params.value = cellstar.encode_auto_params()
+ return cellstar
+
+ @staticmethod
+ def normalized_log(image, sigma):
+ '''Perform the Laplacian of Gaussian transform on the image
+
+ image - 2-d image array
+ mask - binary mask of significant pixels
+ size - length of side of square kernel to use
+ sigma - standard deviation of the Gaussian
+ '''
+ size = int(sigma * 4) + 1
+
+ half_size = size / 2
+ i, j = np.mgrid[-half_size:half_size + 1,
+ -half_size:half_size + 1].astype(float) / float(sigma)
+ distance = (i ** 2 + j ** 2) / 2
+ gaussian = np.exp(-distance)
+ #
+ # Normalize the Gaussian
+ #
+ gaussian = gaussian / np.sum(gaussian)
+
+ log = (distance - 1) * gaussian
+ #
+ # Normalize the kernel to have a sum of zero
+ #
+ log_norm = log - np.mean(log)
+ output = scipy.ndimage.convolve(image, log_norm, mode='nearest')
+ return output
+
+ def preprocess_images(self, input_image, background_image, ignore_mask):
+ # Invert images if required.
+ if not self.background_brighter_then_cell_inside:
+ input_image = 1 - input_image
+ if background_image is not None:
+ background_image = 1 - background_image
+
+ # support for fluorescent images
+ # here it is design question: we assume that the user *should* say
+ # truth about background and inside of cells: insides are brighter
+ # than background (so the bkg and image for fluorescent will be
+ # inverted at this stage)
+ if not self.bright_field_image:
+ sigma = 4
+ factor = 10
+ edge_pixels = IdentifyYeastCells.normalized_log(input_image, sigma)
+ input_image = np.subtract(input_image, factor * edge_pixels)
+
+ return input_image, background_image, ignore_mask
+
+ #
+ # Segmentation of the image into yeast cells.
+ # Returns: yeast cells, yeast cells qualities, background
+ #
+ def segmentation(self, normalized_image, background_pixels, ignore_mask_pixels=None):
+ cellstar = self.prepare_cell_star_object(self.decoded_segmentation_precision_value)
+
+ if self.input_image_file_name is not None:
+ dedicated_image_folder = pj(pref.get_default_output_directory(), self.input_image_file_name)
+ if dedicated_image_folder is not None:
+ cellstar.debug_output_image_path = dedicated_image_folder
+
+ cellstar.set_frame(normalized_image)
+ cellstar.set_background(background_pixels)
+ cellstar.set_mask(ignore_mask_pixels)
+ segmented_image, snakes = cellstar.run_segmentation()
+
+ objects = cellprofiler.object.Objects()
+ objects.segmented = segmented_image
+ objects.unedited_segmented = segmented_image
+ objects.small_removed_segmented = np.zeros(normalized_image.shape)
+
+ self.current_workspace.display_data.segmentation_pixels = objects.segmented
+
+ raw_qualities = [-s.rank for s in snakes]
+ if not raw_qualities == []:
+ raw_interval = min(raw_qualities), max(raw_qualities)
+ logger.info("Qualities are in interval [%.3f,%.3f]." % raw_interval)
+
+ return objects, np.array(raw_qualities), cellstar.images.background
+
+ def fit_parameters(self, input_image, background_image, ignore_mask_image, ground_truth_labels, number_of_steps,
+ update_callback, wait_callback):
+ """
+
+ :param wait_callback: function that wait and potentially updates UI
+ :param update_callback: function that take number of steps completed and return if fitting should be continued
+ """
+ keep_going = True
+
+ self.param_fit_progress = 0
+ self.best_snake_score = 10
+ self.best_rank_score = 1000000000
+ aft_active = []
+ adaptations_stopped = False
+ cellstar = self.prepare_cell_star_object(min(11, self.decoded_segmentation_precision_value))
+ self.best_parameters = cellstar.parameters
+ self.autoadapted_params.value = cellstar.encode_auto_params()
+
+ try:
+ while (keep_going or not adaptations_stopped) and self.param_fit_progress < number_of_steps:
+ # here put one it. of fitting instead
+ wait_callback(0.5)
+
+ # Thread ended with exception so optimisation have to be stopped.
+ if any(aft_active) and aft_active[0].exception is not None:
+ break
+
+ # Clean aft_active from dead threads.
+ while any(aft_active) and (not aft_active[0].is_alive() and aft_active[0].started):
+ aft_active = aft_active[1:]
+
+ # If thread in line start first of them.
+ if any(aft_active) and not aft_active[0].started:
+ # update parameters which may already be changed
+ aft_active[0].update_params(self.best_parameters)
+ aft_active[0].start()
+
+ adaptations_stopped = aft_active == []
+
+ if adaptations_stopped and keep_going and self.param_fit_progress < number_of_steps:
+ aft_active.append(
+ AutoFitterThread(run_pf, self.update_snake_params,
+ input_image, background_image, ignore_mask_image, ground_truth_labels,
+ self.best_parameters,
+ self.decoded_segmentation_precision_value, self.average_cell_diameter.value,
+ self.update_partial_iteration_progress))
+
+ aft_active.append(
+ AutoFitterThread(run_rank_pf, self.update_rank_params,
+ input_image, background_image, ignore_mask_image, ground_truth_labels,
+ self.best_parameters,
+ self.update_partial_iteration_progress))
+
+ # here update params. in the GUI
+ keep_going_update = update_callback(self.param_fit_progress + self.param_fit_progress_partial)
+ keep_going = keep_going and keep_going_update
+
+ finally:
+ update_callback(number_of_steps)
+
+ def fit_parameters_with_ui(self, input_image, background_image, ignore_mask_image, ground_truth_labels):
+ import wx
+ # reading GT from dialog_box.labels[0] and image from self.pixel
+ progress_max = self.autoadaptation_steps.value * 2 # every step consists of: snake params and ranking params fitting
+
+ with wx.ProgressDialog("Fitting parameters..", "Iterations remaining", progress_max * 100,
+ # show percents of change
+ style=wx.PD_CAN_ABORT | wx.PD_ELAPSED_TIME | wx.PD_REMAINING_TIME) as dialog:
+ def update(steps):
+ return dialog.Update(steps * 100)[0]
+
+ def wait(time):
+ return wx.Sleep(int(time + 0.5))
+
+ self.fit_parameters(input_image, background_image, ignore_mask_image, ground_truth_labels, progress_max,
+ update, wait)
+
+ def get_param_fitting_input_images_from_image_set(self):
+ """
+ Try to load images from current workspace. Can be used when fitting is called in test run after segmentation has been
+ run at least once.
+ """
+ if self.fitting_image_set is None:
+ return None
+
+ images = self.get_input_images_from_image_set(self.fitting_image_set)
+ if images is None:
+ return None
+ return images + (None,)
+
+ def get_input_images_from_image_set(self, image_set):
+ try:
+ background_needed = self.background_elimination_strategy == BKG_FILE
+ ignore_mask_needed = self.ignore_mask_image_name.value != cps.LEAVE_BLANK
+
+ background_image = None
+ ignore_mask = None
+
+ # load images from workspace
+ input_image = image_set.get_image(self.input_image_name.value, must_be_grayscale=True).pixel_data
+ if background_needed:
+ background_image = image_set.get_image(self.background_image_name.value,
+ must_be_grayscale=True).pixel_data
+ if ignore_mask_needed:
+ ignore_mask_image = image_set.get_image(self.ignore_mask_image_name)
+ ignore_mask = ignore_mask_image.pixel_data > 0
+
+ return input_image, background_image, ignore_mask
+ except Exception as ex:
+ logger.info("Could not use image from workspace.image_set because: " + str(ex))
+ return None
+
+ @staticmethod
+ def load_image_grayscale(path):
+ image = bioformats.load_image(path)
+ if image.ndim == 3:
+ image = np.sum(image, 2) / image.shape[2]
+ return image
+
+ def get_param_fitting_input_images_from_user(self):
+ import wx
+ def get_file_path(message):
+ with wx.FileDialog(None,
+ message=message,
+ wildcard="Image file (*.tif,*.tiff,*.jpg,*.jpeg,*.png,*.gif,*.bmp)|*.tif;*.tiff;*.jpg;*.jpeg;*.png;*.gif;*.bmp|*.* (all files)|*.*",
+ style=wx.FD_OPEN) as dlg:
+ if dlg.ShowModal() == wx.ID_OK:
+ return dlg.Path
+ else:
+ return None
+
+ input_image = None
+ background_image = None
+ ignore_mask = None
+ labels = None
+
+ background_needed = self.background_elimination_strategy == BKG_FILE
+ ignore_mask_needed = self.ignore_mask_image_name.value != cps.LEAVE_BLANK
+
+ message = "In order to autoadapt parameters you need to provide:\n- input image"
+ if background_needed:
+ message += "\n- background image"
+
+ if ignore_mask_needed:
+ message += "\n- ignore mask image"
+
+ with wx.MessageDialog(None, message, "Select images for autoadapt", wx.OK | wx.ICON_INFORMATION) as dlg:
+ dlg.ShowModal()
+
+ input_image_path = get_file_path("Select sample input image")
+ if input_image_path is None:
+ return None
+ else:
+ input_image = self.load_image_grayscale(input_image_path)
+ label_path = input_image_path + ".lab.png" # if file attached load labels from file
+ if isfile(label_path):
+ labels = (self.load_image_grayscale(label_path) * 255).astype(int)
+
+ if background_needed:
+ background_image_path = get_file_path("Select background image")
+ if background_image_path is None:
+ return None
+ else:
+ background_image = self.load_image_grayscale(background_image_path)
+
+ if ignore_mask_needed:
+ ignore_mask_path = get_file_path("Select ignore mask image")
+ if ignore_mask_path is None:
+ return None
+ else:
+ ignore_mask = self.load_image_grayscale(ignore_mask_path) > 0
+
+ return input_image, background_image, ignore_mask, labels
+
+ def ground_truth_editor(self):
+ import wx
+ from cellprofiler.gui.editobjectsdlg import EditObjectsDialog
+
+ '''Display a UI for GT editing'''
+ ### check if user want to use last pipeline run images
+ pipeline_imagery = self.get_param_fitting_input_images_from_image_set()
+ if pipeline_imagery is not None:
+ # ask if user want to use it
+ with wx.MessageDialog(None,
+ "Images used in previous pipeline run are available. Do you want to use them for autoadapting?",
+ "Using pipeline images", wx.YES_NO | wx.ICON_QUESTION) as dlg:
+ use_pipeline = dlg.ShowModal()
+ if use_pipeline == wx.ID_NO:
+ pipeline_imagery = None
+
+ ### opening file dialogs
+ input_data = pipeline_imagery or self.get_param_fitting_input_images_from_user()
+ if input_data is None:
+ return
+ input_image, background_image, ignore_mask, labels = input_data
+
+ ### opening GT editor
+ title = "Please mark few representative cells to allow for autoadapting parameters. \n"
+ title += " \n"
+ title += 'Press "F" to being freehand drawing.\n'
+ title += "Click Help for full instructions."
+
+ if labels is None or not labels.any():
+ edit_labels = [np.zeros(input_image.shape[:2], int)]
+
+ if getattr(self, "last_labeling", None) is not None:
+ if self.last_labeling[0] == hash(abs(np.sum(input_image))):
+ edit_labels = [self.last_labeling[1]]
+
+ ## two next lines are hack from Lee
+ edit_labels[0][0, 0] = 1
+ edit_labels[0][-2, -2] = 1
+ with EditObjectsDialog(
+ input_image, edit_labels, False, title) as dialog_box:
+ hack_add_from_file_into_EditObjects(dialog_box)
+ result = dialog_box.ShowModal()
+ if result != wx.OK:
+ return None
+ labels = dialog_box.labels[0]
+ ## two next lines are hack from Lee
+ labels[0, 0] = 0
+ labels[-2, -2] = 0
+
+ self.last_labeling = (hash(abs(np.sum(input_image))), labels)
+
+ # check if the user provided GT
+ if not labels.any():
+ with wx.MessageDialog(None,
+ "Please correctly select at least one cell. Otherwise, parameters can not be autoadapted!",
+ "Warning!", wx.OK | wx.ICON_WARNING) as dlg:
+ dlg.ShowModal()
+ return
+
+ input_processed, background_processed, ignore_mask_processed = self.preprocess_images(input_image,
+ background_image,
+ ignore_mask)
+
+ self.fit_parameters_with_ui(input_processed, background_processed, ignore_mask_processed, labels)
+
+ def update_partial_iteration_progress(self, fraction):
+ self.param_fit_progress_partial = min(0.99, max(self.param_fit_progress_partial, fraction))
+
+ def update_snake_params(self, new_parameters, new_snake_score):
+ if new_snake_score < self.best_snake_score:
+ self.best_snake_score = new_snake_score
+ self.best_parameters = new_parameters
+ self.best_rank_score = 1000000000 # clear best ranking params as it no longer valid for different snakes
+ if self.autoadapted_params.value != Segmentation.encode_auto_params_from_all_params(new_parameters):
+ self.autoadapted_params.value = Segmentation.encode_auto_params_from_all_params(new_parameters)
+ logger.info("New auto parameters applied.")
+ else:
+ logger.info(
+ "New auto parameters (%f) are not better than current (%f)." % (new_snake_score, self.best_snake_score))
+ self.param_fit_progress += 1
+ self.param_fit_progress_partial = 0
+
+ def update_rank_params(self, new_parameters, new_rank_score):
+ if new_rank_score < self.best_rank_score:
+ self.best_rank_score = new_rank_score
+ self.best_parameters = new_parameters
+ if self.autoadapted_params.value != Segmentation.encode_auto_params_from_all_params(new_parameters):
+ self.autoadapted_params.value = Segmentation.encode_auto_params_from_all_params(new_parameters)
+ logger.info("New auto ranking parameters applied.")
+ else:
+ logger.info("New auto ranking parameters (%f) are not better than current (%f)." % (
+ new_rank_score, self.best_rank_score))
+ self.param_fit_progress += 1
+ self.param_fit_progress_partial = 0
+
+
+class AutoFitterThread(threading.Thread):
+ def __init__(self, target, callback, *args):
+ self._target = target
+ self._args = list(args)
+ self._callback = callback
+ self.started = False
+ self.mock_alive = False
+ self.exception = None
+ super(AutoFitterThread, self).__init__()
+ pass
+
+ def update_params(self, new_params):
+ self._args[4] = new_params
+
+ def run(self):
+ try:
+ self.started = True
+ params, score = self._target(*self._args)
+ self._callback(params, score)
+ except:
+ self.exception = True
+ raise
+
+ def start(self):
+ # check if call on the same thread because of explorer
+ if not explorer_expected():
+ self.started = True
+ super(AutoFitterThread, self).start()
+ else:
+ self.mock_alive = True
+ self.run()
+ self.mock_alive = False
+
+ def is_alive(self):
+ if not explorer_expected():
+ return super(AutoFitterThread, self).is_alive()
+ else:
+ return self.mock_alive
+
+ def kill(self):
+ pass
diff --git a/tests/test_identifyyeastcells.py b/tests/test_identifyyeastcells.py
new file mode 100644
index 00000000..b4bce9a6
--- /dev/null
+++ b/tests/test_identifyyeastcells.py
@@ -0,0 +1,724 @@
+"""test_identifyyeastcells.py: test the IdentifyYeastCells module
+
+CellProfiler is distributed under the GNU General Public License.
+See the accompanying file LICENSE for details.
+
+Copyright (c) 2003-2009 Massachusetts Institute of Technology
+Copyright (c) 2009-2015 Broad Institute
+All rights reserved.
+
+Credits (coding)
+Filip Mroz, Adam Kaczmarek, Szymon Stoma.
+
+Website: http://www.cellprofiler.org
+"""
+
+import StringIO
+import ast
+import time
+import unittest
+
+import numpy as np
+import scipy.ndimage
+
+# run in headless mode so wx is not required
+import cellprofiler.preferences
+cellprofiler.preferences.set_headless()
+
+import cellprofiler.image as cpi
+import cellprofiler.measurement as cpmeas
+import identifyyeastcells as YS
+import cellprofiler.object as cpo
+import cellprofiler.pipeline
+import cellprofiler.setting
+from cellprofiler.workspace import Workspace
+
+IMAGE_NAME = "my_image"
+BACKGROUND_IMAGE_NAME = "background_image"
+MASK_IMAGE_NAME = "mask_image"
+OBJECTS_NAME = "my_objects"
+BINARY_IMAGE_NAME = "binary_image"
+MASKING_OBJECTS_NAME = "masking_objects"
+MEASUREMENT_NAME = "my_measurement"
+
+
+class test_IdentifyYeastCells(unittest.TestCase):
+ def load_error_handler(self, caller, event):
+ if isinstance(event, cellprofiler.pipeline.LoadExceptionEvent):
+ self.fail(event.error.message)
+
+ def assertRange(self, range_start, range_end, value):
+ self.assertGreaterEqual(value, range_start)
+ self.assertLessEqual(value, range_end)
+
+ def make_workspace(self, image,
+ mask=None,
+ labels=None,
+ binary_image=None):
+ '''Make a workspace and IdentifyPrimaryObjects module
+
+ image - the intensity image for thresholding
+
+ mask - if present, the "don't analyze" mask of the intensity image
+
+ labels - if thresholding per-object, the labels matrix needed
+
+ binary_image - if thresholding using a binary image, the image
+ '''
+ module = YS.IdentifyYeastCells()
+ module.module_num = 1
+ module.input_image_name.value = IMAGE_NAME
+ module.object_name.value = OBJECTS_NAME
+ module.binary_image.value = BINARY_IMAGE_NAME
+ module.masking_objects.value = MASKING_OBJECTS_NAME
+
+ pipeline = cellprofiler.pipeline.Pipeline()
+ pipeline.add_module(module)
+ m = cpmeas.Measurements()
+ cpimage = cpi.Image(image, mask=mask)
+ m.add(IMAGE_NAME, cpimage)
+ if binary_image is not None:
+ m.add(BINARY_IMAGE_NAME, cpi.Image(binary_image))
+ object_set = cpo.ObjectSet()
+ if labels is not None:
+ o = cpo.Objects()
+ o.segmented = labels
+ object_set.add_objects(o, MASKING_OBJECTS_NAME)
+ workspace = cellprofiler.workspace.Workspace(
+ pipeline, module, m, object_set, m, None)
+ return workspace, module
+
+ def test_00_00_init(self):
+ x = YS.IdentifyYeastCells()
+
+ def test_01_00_test_zero_objects(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 5
+ img = np.zeros((25, 25))
+ image = cpi.Image(img, file_name="test_01_00_test_zero_objects")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ self.assertEqual(len(object_set.object_names), 1)
+ self.assertTrue(OBJECTS_NAME in object_set.object_names)
+ objects = object_set.get_objects(OBJECTS_NAME)
+ segmented = objects.segmented
+ self.assertTrue(np.all(segmented == 0))
+ self.assertTrue("Image" in measurements.get_object_names())
+ self.assertTrue(OBJECTS_NAME in measurements.get_object_names())
+ self.assertTrue("Count_" + OBJECTS_NAME in measurements.get_feature_names("Image"))
+ count = measurements.get_current_measurement("Image", "Count_" + OBJECTS_NAME)
+ self.assertEqual(count, 0)
+ self.assertTrue("Location_Center_X" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ self.assertTrue(isinstance(location_center_x, np.ndarray))
+ self.assertEqual(np.product(location_center_x.shape), 0)
+ self.assertTrue("Location_Center_Y" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_y = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_Y")
+ self.assertTrue(isinstance(location_center_y, np.ndarray))
+ self.assertEqual(np.product(location_center_y.shape), 0)
+
+ def test_01_01_test_one_object(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 10
+ img = convert_to_brightfield(get_one_cell_mask(), False)
+ image = cpi.Image(img, file_name="test_01_01_test_one_object")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ self.assertEqual(len(object_set.object_names), 1)
+ self.assertTrue(OBJECTS_NAME in object_set.object_names)
+ objects = object_set.get_objects(OBJECTS_NAME)
+ segmented = objects.segmented
+ self.assertTrue(is_segmentation_correct(get_one_cell_mask(), segmented))
+ self.assertTrue("Image" in measurements.get_object_names())
+ self.assertTrue(OBJECTS_NAME in measurements.get_object_names())
+ self.assertTrue("Features_Quality" in measurements.get_feature_names(OBJECTS_NAME))
+ quality = measurements.get_current_measurement(OBJECTS_NAME, "Features_Quality")
+ self.assertTrue(quality[0] > 0)
+ self.assertTrue("Location_Center_Y" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_y = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_Y")
+ self.assertTrue(isinstance(location_center_y, np.ndarray))
+ self.assertEqual(np.product(location_center_y.shape), 1)
+ self.assertTrue(location_center_y[0] > 8)
+ self.assertTrue(location_center_y[0] < 12)
+ self.assertTrue("Location_Center_X" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ self.assertTrue(isinstance(location_center_x, np.ndarray))
+ self.assertEqual(np.product(location_center_x.shape), 1)
+ self.assertTrue(location_center_x[0] > 13)
+ self.assertTrue(location_center_x[0] < 16)
+ columns = x.get_measurement_columns(pipeline)
+ for object_name in (cpmeas.IMAGE, OBJECTS_NAME):
+ ocolumns = [x for x in columns if x[0] == object_name]
+ features = measurements.get_feature_names(object_name)
+ self.assertEqual(len(ocolumns), len(features))
+ self.assertTrue(all([column[1] in features for column in ocolumns]))
+
+ def test_01_02_test_two_bright_objects(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 10
+ img = convert_to_brightfield(get_two_cell_mask(), False)
+ image = cpi.Image(img, file_name="test_01_02_test_two_bright_objects")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ self.assertEqual(len(object_set.object_names), 1)
+ self.assertTrue(OBJECTS_NAME in object_set.object_names)
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertTrue(is_segmentation_correct(get_two_cell_mask(), objects.segmented))
+ self.assertTrue("Image" in measurements.get_object_names())
+ self.assertTrue(OBJECTS_NAME in measurements.get_object_names())
+ self.assertTrue("Features_Quality" in measurements.get_feature_names(OBJECTS_NAME))
+ quality = measurements.get_current_measurement(OBJECTS_NAME, "Features_Quality")
+ self.assertTrue(len(quality) == 2)
+ location_center_y = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_Y")
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ positions = sorted(zip(location_center_x, location_center_y))
+ self.assertEqual(2, len(positions))
+
+ self.assertRange(3, 18, positions[0][0])
+ self.assertRange(25, 45, positions[0][1])
+
+ self.assertRange(20, 40, positions[1][0])
+ self.assertRange(5, 25, positions[1][1])
+
+ def test_01_03_test_two_dark_objects(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = True
+ x.average_cell_diameter.value = 10
+ img = convert_to_brightfield(get_two_cell_mask(), True)
+ image = cpi.Image(img, file_name="test_01_03_test_two_dark_objects")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ self.assertEqual(len(object_set.object_names), 1)
+ self.assertTrue(OBJECTS_NAME in object_set.object_names)
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertTrue(is_segmentation_correct(get_two_cell_mask(), objects.segmented))
+ self.assertTrue("Image" in measurements.get_object_names())
+ self.assertTrue(OBJECTS_NAME in measurements.get_object_names())
+ self.assertTrue("Features_Quality" in measurements.get_feature_names(OBJECTS_NAME))
+ quality = measurements.get_current_measurement(OBJECTS_NAME, "Features_Quality")
+ self.assertTrue(len(quality) == 2)
+ self.assertTrue("Location_Center_Y" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_y = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_Y")
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ positions = sorted(zip(location_center_x, location_center_y))
+ self.assertEqual(2, len(positions))
+
+ self.assertRange(3, 18, positions[0][0])
+ self.assertRange(25, 45, positions[0][1])
+
+ self.assertRange(20, 40, positions[1][0])
+ self.assertRange(5, 25, positions[1][1])
+
+ def test_01_04_test_two_flu_bright_objects(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.bright_field_image.value = False
+ x.average_cell_diameter.value = 10
+ img = convert_to_fluorescent(get_two_cell_mask(), False)
+ image = cpi.Image(img, file_name="test_01_04_test_two_flu_bright_objects")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ self.assertEqual(len(object_set.object_names), 1)
+ self.assertTrue(OBJECTS_NAME in object_set.object_names)
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertTrue(is_segmentation_correct(get_two_cell_mask(), objects.segmented))
+ self.assertTrue("Image" in measurements.get_object_names())
+ self.assertTrue(OBJECTS_NAME in measurements.get_object_names())
+ self.assertTrue("Features_Quality" in measurements.get_feature_names(OBJECTS_NAME))
+ quality = measurements.get_current_measurement(OBJECTS_NAME, "Features_Quality")
+ self.assertTrue(len(quality) == 2)
+ self.assertTrue("Location_Center_X" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_y = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_Y")
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ positions = sorted(zip(location_center_x, location_center_y))
+ self.assertEqual(2, len(positions))
+
+ self.assertRange(3, 18, positions[0][0])
+ self.assertRange(25, 45, positions[0][1])
+
+ self.assertRange(20, 40, positions[1][0])
+ self.assertRange(5, 25, positions[1][1])
+
+ def test_01_05_test_two_flu_dark_objects(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = True
+ x.bright_field_image.value = False
+ x.average_cell_diameter.value = 10
+ img = convert_to_fluorescent(get_two_cell_mask(), True)
+ image = cpi.Image(img, file_name="test_01_05_test_two_flu_dark_objects")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ self.assertEqual(len(object_set.object_names), 1)
+ self.assertTrue(OBJECTS_NAME in object_set.object_names)
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertTrue(is_segmentation_correct(get_two_cell_mask(), objects.segmented))
+ self.assertTrue("Image" in measurements.get_object_names())
+ self.assertTrue(OBJECTS_NAME in measurements.get_object_names())
+ self.assertTrue("Features_Quality" in measurements.get_feature_names(OBJECTS_NAME))
+ quality = measurements.get_current_measurement(OBJECTS_NAME, "Features_Quality")
+ self.assertTrue(len(quality) == 2)
+ self.assertTrue("Location_Center_X" in measurements.get_feature_names(OBJECTS_NAME))
+ location_center_y = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_Y")
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ positions = sorted(zip(location_center_x, location_center_y))
+ self.assertEqual(2, len(positions))
+
+ self.assertRange(3, 18, positions[0][0])
+ self.assertRange(25, 45, positions[0][1])
+
+ self.assertRange(20, 40, positions[1][0])
+ self.assertRange(5, 25, positions[1][1])
+
+ def test_01_06_fill_holes(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 5
+
+ img = np.zeros((40, 40))
+ draw_disc(img, (10, 10), 7, .5)
+ draw_disc(img, (30, 30), 7, .5)
+ img[10, 10] = 0
+ img[30, 30] = 0
+ image = cpi.Image(img, file_name="test_01_06_fill_holes")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertTrue(objects.segmented[10, 10] > 0)
+ self.assertTrue(objects.segmented[30, 30] > 0)
+
+ def test_01_07_extreme_params(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 14
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 77
+ img = get_two_cell_mask()
+ draw_disc(img, (5, 5), 2, 0.7)
+ draw_disc(img, (35, 11), 3, 0.2)
+ img = convert_to_brightfield(img, False)
+ image = cpi.Image(img, file_name="test_01_07_extreme_params")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ segmented = objects.segmented
+ self.assertTrue(np.all(segmented == 0)) # no found because of parameters (no foreground)
+ self.assertEqual(0, objects.count)
+
+ def test_02_01_discard_large(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 30
+ x.min_cell_area.value = 100
+ x.max_cell_area.value = 1000
+ x.advanced_cell_filtering.value = True
+
+ img = np.ones((200, 200)) * 0.5
+ draw_brightfield_cell(img, 100, 100, 20, False)
+ draw_brightfield_cell(img, 25, 25, 10, False)
+ draw_brightfield_cell(img, 150, 150, 15, False)
+ image = cpi.Image(img, file_name="test_02_01_discard_large")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertEqual(objects.segmented[25, 25] > 0, 1, "The small object was not there")
+ self.assertEqual(objects.segmented[150, 150] > 0, 1, "The medium object was not there")
+ self.assertEqual(objects.segmented[100, 100] > 0, 0, "The large object was not filtered out")
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ self.assertTrue(isinstance(location_center_x, np.ndarray))
+ self.assertEqual(np.product(location_center_x.shape), 2)
+
+ def test_02_02_discard_small(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.segmentation_precision.value = 11
+ x.input_image_name.value = IMAGE_NAME
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 30
+ x.min_cell_area.value = 500
+ x.max_cell_area.value = 5000
+ x.advanced_cell_filtering.value = True
+
+ img = np.ones((200, 200)) * 0.5
+ draw_brightfield_cell(img, 100, 100, 20, False)
+ draw_brightfield_cell(img, 25, 25, 10, False)
+ draw_brightfield_cell(img, 150, 150, 15, False)
+ image = cpi.Image(img, file_name="test_02_02_discard_small")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertEqual(objects.segmented[25, 25] > 0, 0, "The small object was not filtered out")
+ self.assertEqual(objects.segmented[150, 150] > 0, 1, "The medium object was not there")
+ self.assertEqual(objects.segmented[100, 100] > 0, 1, "The large object was not there")
+ location_center_x = measurements.get_current_measurement(OBJECTS_NAME, "Location_Center_X")
+ self.assertTrue(isinstance(location_center_x, np.ndarray))
+ self.assertEqual(np.product(location_center_x.shape), 2)
+
+ def test_02_03_use_background_image(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.input_image_name.value = IMAGE_NAME
+ x.segmentation_precision.value = 11
+ x.background_image_name.value = BACKGROUND_IMAGE_NAME
+ x.background_elimination_strategy.value = YS.BKG_FILE
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 30
+
+ img = np.ones((200, 200)) * 0.5
+ draw_brightfield_cell(img, 100, 100, 20, False)
+ draw_brightfield_cell(img, 25, 25, 10, False)
+ draw_brightfield_cell(img, 150, 150, 15, False) # background blob
+
+ bkg = np.ones((200, 200)) * 0.5
+ draw_brightfield_cell(bkg, 150, 150, 15, False) # background blob
+
+ image = cpi.Image(img, file_name="test_02_03_use_background_image")
+ background = cpi.Image(bkg)
+
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ image_set.providers.append(cpi.VanillaImageProvider(BACKGROUND_IMAGE_NAME, background))
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertEqual(objects.segmented[25, 25] > 0, 1, "The small object was not there")
+ self.assertEqual(objects.segmented[100, 100] > 0, 1, "The large object was not there")
+ self.assertEqual(objects.segmented[150, 150] > 0, 0, "The background blob was not filtered out")
+
+ def test_02_04_mask_input_image(self):
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.input_image_name.value = IMAGE_NAME
+ x.segmentation_precision.value = 11
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 30
+
+ img = np.ones((200, 200)) * 0.5
+ draw_brightfield_cell(img, 100, 100, 20, False)
+ draw_brightfield_cell(img, 25, 25, 10, False)
+ img[0:10, 0:10] = 1
+ img[180:200, 180:200] = 0
+
+ msk = np.zeros((200, 200))
+ msk[0:10, 0:10] = 1
+ msk[180:200, 180:200] = 1
+
+ image = cpi.Image(img, file_name="test_02_04_mask_input_image")
+ mask = cpi.Image(msk)
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ image_set.providers.append(cpi.VanillaImageProvider(MASK_IMAGE_NAME, mask))
+
+ # first try without masking
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertEqual(0, objects.segmented.max(), "Cells should not be found due to the distractors")
+
+ # now if we use masking option we should find these cells
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.ignore_mask_image_name.value = MASK_IMAGE_NAME
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertEqual(objects.segmented[25, 25] > 0, 1, "The small object was not there")
+ self.assertEqual(objects.segmented[100, 100] > 0, 1, "The large object was not there")
+
+ def test_03_01_simple_fitting(self):
+ # reduce depth of fitting to speed up testing
+ import cellstar.parameter_fitting.pf_process as process
+ import cellstar.parameter_fitting.pf_runner as runner
+ process.SEARCH_LENGTH_NORMAL = 20
+ # test multicore but only two cores
+ process.get_max_workers = lambda: 2
+ runner.get_max_workers = lambda: 2
+ # make it deterministic
+ np.random.seed(1)
+
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+
+ x.input_image_name.value = IMAGE_NAME
+ x.segmentation_precision.value = 9 # so that it is faster for fitting
+ x.maximal_cell_overlap.value = 0.4
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 30
+ x.autoadaptation_steps.value = 1
+
+ img = np.ones((200, 200)) * 0.5
+ draw_brightfield_cell(img, 100, 100, 15, False)
+ draw_brightfield_cell(img, 120, 120, 15, False)
+ draw_brightfield_cell(img, 110, 70, 15, False)
+ draw_brightfield_cell(img, 160, 160, 10, False)
+ draw_disc(img, (100, 100), 15, .65)
+ draw_disc(img, (120, 120), 15, .65)
+ draw_disc(img, (110, 70), 15, .65)
+ draw_disc(img, (160, 160), 10, .65)
+ img = img + np.random.normal(3., 0.01, img.shape)
+ img = scipy.ndimage.gaussian_filter(img, 3)
+
+ label = np.zeros((200, 200), dtype=int)
+ draw_disc(label, (100, 100), 15, 1)
+ draw_disc(label, (110, 70), 15, 2)
+
+ image = cpi.Image(img, file_name="test_03_01_simple_fitting")
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+
+ old_params = ast.literal_eval(x.autoadapted_params.value)
+ input_processed, background_processed, ignore_mask_processed = x.preprocess_images(img, None, None)
+ x.fit_parameters(input_processed, background_processed, ignore_mask_processed, label,
+ x.autoadaptation_steps.value * 2, lambda x: True, lambda secs: time.sleep(secs))
+ new_params = ast.literal_eval(x.autoadapted_params.value)
+ self.assertNotEqual(old_params[0], new_params[0])
+ self.assertNotEqual(old_params[1], new_params[1])
+
+ # now if we use new parameters option we should find these cells
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.segmentation_precision.value = 11
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ self.assertEqual(4, objects.segmented.max())
+ colours = sorted([objects.segmented[100, 100], objects.segmented[120, 120], objects.segmented[110, 70],
+ objects.segmented[160, 160]])
+ self.assertEqual(colours[0], 1)
+ self.assertEqual(colours[1], 2)
+ self.assertEqual(colours[2], 3)
+ self.assertEqual(colours[3], 4)
+
+ def test_03_02_fitting_background_masked(self):
+ # reduce depth of fitting to speed up testing
+ import cellstar.parameter_fitting.pf_process as process
+ import cellstar.parameter_fitting.pf_runner as runner
+ process.SEARCH_LENGTH_NORMAL = 20
+ # test one core
+ process.get_max_workers = lambda: 1
+ runner.get_max_workers = lambda : 1
+ # make it deterministic
+ np.random.seed(1)
+
+ x = YS.IdentifyYeastCells()
+ x.object_name.value = OBJECTS_NAME
+ x.input_image_name.value = IMAGE_NAME
+ x.segmentation_precision.value = 9 # so that it is faster for fitting
+ x.maximal_cell_overlap.value = 0.4
+ x.background_brighter_then_cell_inside.value = False
+ x.average_cell_diameter.value = 12
+ x.autoadaptation_steps.value = 1
+
+ x.background_image_name.value = BACKGROUND_IMAGE_NAME
+ x.background_elimination_strategy.value = YS.BKG_FILE
+ x.ignore_mask_image_name.value = MASK_IMAGE_NAME
+
+ img = np.ones((50, 50)) * 0.5
+ draw_brightfield_cell(img, 7, 7, 5, False)
+ draw_brightfield_cell(img, 25, 28, 5, False)
+ draw_brightfield_cell(img, 15, 16, 5, False)
+ draw_brightfield_cell(img, 40, 40, 4, False)
+ draw_disc(img, (7, 7), 5, .65)
+ draw_disc(img, (25, 28), 5, .65)
+ draw_disc(img, (15, 16), 5, .65)
+ draw_disc(img, (40, 40), 4, .65)
+ img = img + np.random.normal(0.5, 0.01, img.shape)
+ img = scipy.ndimage.gaussian_filter(img, 2)
+
+ # bright flares
+ draw_disc(img, (40, 10), 10, 1.5)
+ ignore_mask = np.zeros((50, 50), dtype=bool)
+ draw_disc(ignore_mask, (40, 10), 10, 1)
+
+ # dark areas
+ draw_disc(img, (50, 25), 10, 0.0)
+ background_mask = np.ones((50, 50)) * 0.5
+ draw_disc(background_mask, (50, 25), 10, 0.0)
+
+ label = np.zeros((50, 50), dtype=int)
+ draw_disc(label, (10, 10), 7, 1)
+ draw_disc(label, (25, 20), 7, 2)
+
+ image = cpi.Image(img, file_name="test_03_02_fitting_background_masked")
+ mask = cpi.Image(ignore_mask)
+ background = cpi.Image(background_mask)
+
+ image_set_list = cpi.ImageSetList()
+ image_set = image_set_list.get_image_set(0)
+ image_set.providers.append(cpi.VanillaImageProvider(IMAGE_NAME, image))
+ image_set.providers.append(cpi.VanillaImageProvider(MASK_IMAGE_NAME, mask))
+ image_set.providers.append(cpi.VanillaImageProvider(BACKGROUND_IMAGE_NAME, background))
+
+ old_params = ast.literal_eval(x.autoadapted_params.value)
+ input_processed, background_processed, ignore_mask_processed = x.preprocess_images(img, background_mask, ignore_mask)
+ x.fit_parameters(input_processed, background_processed, ignore_mask_processed, label,
+ x.autoadaptation_steps.value * 2, lambda x: True, lambda secs: time.sleep(secs))
+ new_params = ast.literal_eval(x.autoadapted_params.value)
+ self.assertNotEqual(old_params[0], new_params[0])
+ self.assertNotEqual(old_params[1], new_params[1])
+
+ # now if we use new parameters option we should find these cells
+ object_set = cpo.ObjectSet()
+ measurements = cpmeas.Measurements()
+ pipeline = cellprofiler.pipeline.Pipeline()
+ x.segmentation_precision.value = 11
+ x.run(Workspace(pipeline, x, image_set, object_set, measurements, None))
+ objects = object_set.get_objects(OBJECTS_NAME)
+ # 3 or four objects are acceptable
+ self.assertLessEqual(3, objects.segmented.max())
+ self.assertLessEqual(objects.segmented.max(), 4)
+
+ colours = sorted([objects.segmented[10, 10], objects.segmented[25, 25], objects.segmented[40, 40]])
+ self.assertEqual(colours[0], 1)
+ self.assertEqual(colours[1], 2)
+ self.assertEqual(colours[2], 3)
+
+
+def add_noise(img, fraction):
+ '''Add a fractional amount of noise to an image to make it look real'''
+ np.random.seed(0)
+ noise = np.random.uniform(low=1 - fraction / 2, high=1 + fraction / 2,
+ size=img.shape)
+ return img * noise
+
+
+def get_one_cell_mask():
+ img = np.zeros((30, 30))
+ draw_disc(img, (10, 15), 5, 1)
+ return img
+
+
+def get_two_cell_mask():
+ img = np.zeros((50, 50))
+ draw_disc(img, (10, 25), 5, 1)
+ draw_disc(img, (25, 15), 5, 1)
+ return img
+
+
+def convert_to_brightfield(img, content_dark):
+ if (content_dark):
+ img *= 0.3
+ else:
+ img *= 0.6
+ # get ring with dilation (5x5 radius)
+ ring = (scipy.ndimage.morphology.binary_dilation(img, np.ones((3, 3))) - (img > 0))
+ img[ring] = .8
+ # fill rest with background
+ img[img == 0] = 0.5
+ return add_noise(img, 0.000)
+
+
+def draw_brightfield_cell(img, x, y, radius, content_dark=True):
+ draw_disc(img, (x, y), radius + 2, .8)
+ if (content_dark):
+ draw_disc(img, (x, y), radius, .3)
+ else:
+ draw_disc(img, (x, y), radius, .6)
+ return img
+
+
+def convert_to_fluorescent(img, content_dark):
+ if content_dark:
+ img = 1 - img + (img * 0.1)
+ else:
+ img *= 0.9
+ img = scipy.ndimage.gaussian_filter(img, sigma=2)
+ return add_noise(img, .000)
+
+
+def is_segmentation_correct(ground_truth, segmentation):
+ return are_masks_similar(segmentation > 0, ground_truth > 0)
+
+
+def are_masks_similar(a, b):
+ return 1.0 - (a & b).sum() / float((a | b).sum()) < 0.5
+
+
+def draw_disc(img, center, radius, value):
+ x, y = np.mgrid[0:img.shape[0], 0:img.shape[1]]
+ distance = np.sqrt((x - center[0]) * (x - center[0]) + (y - center[1]) * (y - center[1]))
+ img[distance <= radius] = value