Python skimage.feature.blob_doh() Examples

def blob_image(image):
    #img_path    =   '../1.3.6.1.4.1.14519.5.2.1.6279.6001.100332161840553388986847034053.mhd'
    # slice = numpy_image[240,:,:]
    # normalized = norm.normalize(return_surrounding([240,240,240],numpy_image, 240))
    # thresholded = threshold_by_histogram(normalized)
    # blobs = label_image(thresholded)
    # show_images([blobs, normalized, thresholded])

    # normalized3d = norm.normalize(numpy_image)
    # thresholded3d = threshold_by_histogram(normalized3d)
    list = []
    image = norm.normalize(image)
    #print "normalized and thresholded"

    for z, slice in tqdm(enumerate(image)):
        blobs = feature.blob_doh(slice)
        #print blobs.shape
        blob_coords = np.zeros((len(blobs),3))
        for i, blob in enumerate(blobs):
            blob_coords[i] = [z, blob[0], blob[1]]
        list.append(blob_coords)

    #print list
    return list 

def __init__(
            self,
            min_sigma: Union[Number, Tuple[Number, ...]],
            max_sigma: Union[Number, Tuple[Number, ...]],
            num_sigma: int,
            threshold: Number,
            overlap: float = 0.5,
            measurement_type='max',
            is_volume: bool = True,
            detector_method: str = 'blob_log',
            exclude_border: Optional[int] = None,
    ) -> None:

        self.min_sigma = min_sigma
        self.max_sigma = max_sigma
        self.num_sigma = num_sigma
        self.threshold = threshold
        self.overlap = overlap
        self.is_volume = is_volume
        self.measurement_function = self._get_measurement_function(measurement_type)
        self.exclude_border = exclude_border
        try:
            self.detector_method = blob_detectors[detector_method]
        except ValueError:
            raise ValueError("Detector method must be one of {blob_log, blob_dog, blob_doh}") 

def blob_detection(self, img):
        blobs = blob_doh(img, max_sigma=80, threshold=.001)
        return blobs 

def main():
    numberOfImages = 11;

    # TODO: AUTOMATICALLY GET NUMBER OF IMAGES
    # Get number of images. Remeber to divide by 2 as for every relevant image,
    # theres also the comparison image.
    # if ".DS_Store" in os.listdir("Wheat_ROIs"):
    #     numberOfImages = (len(os.listdir("Wheat_ROIs")) - 1)/2;
    # else:
    #     numberOfImages = len(os.listdir("Wheat_ROIs"))/2;

    # For each ROI image in folder
    for i in tqdm.tqdm(range(1, numberOfImages+1)):
        # Load image
        filename = "../Wheat_ROIs/{:03d}_ROI.png".format(i);
        img = misc.imread(filename);
        img_gray = rgb2gray(img);

        # Detect blobs. See http://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.blob_doh
        # for function documentation
        blobs = blob_doh(img_gray, min_sigma=1, max_sigma=100, threshold=.01)

        # Display blobs on image and save image
        fig, ax = plt.subplots()
        plt.title("Number of Blobs Detected: {}".format(blobs.shape[0]))
        plt.grid(False)
        ax.imshow(img, interpolation='nearest')
        for blob in blobs:
            y, x, r = blob
            c = plt.Circle((x, y), r, color='red', linewidth=2, fill=False)
            ax.add_patch(c)
        fig.savefig("../Wheat_ROIs/{:03d}_Blob.png".format(i)) 

def blob_image_multiscale2(image, type=0,scale=2):
    # function that return a list of blob_coordinates, 0 = dog, 1 = doh, 2 = log
    list = []
    image = norm.normalize(image)
    for z, slice in tqdm(enumerate(image)):
        # init list of different sigma/zoom blobs
        featureblobs = []
        # x = 0,1,2,3,4
        if scale == 2:
            # for x in xrange(0,6):
            #     if type == 0:
            #         featureblobs.append(feature.blob_dog(slice, 2**x, 2**x))
            #     if type == 1:
            #         featureblobs.append(feature.blob_doh(slice, 2**x, 2**x))
            #     if type == 2:
            #         featureblobs.append(feature.blob_log(slice, 2**x, 2**x))
            for x in xrange(0,5):
                if type == 0:
                    featureblobs.append(feature.blob_dog(slice, 2**x, 2**(x+1)))
                if type == 1:
                    featureblobs.append(feature.blob_doh(slice, 2**x, 2**(x+1)))
                if type == 2:
                    featureblobs.append(feature.blob_log(slice, 2**x, 2**(x+1),16,.1))
        else:
            for x in xrange(0,4):
                if type == 0:
                    featureblobs.append(feature.blob_dog(slice, 3**x, 3**x))
                if type == 1:
                    featureblobs.append(feature.blob_doh(slice, 3**x, 3**x))
                if type == 2:
                    featureblobs.append(feature.blob_log(slice, 3**x, 3**x))
        # init list of blob coords
        blob_coords = []
        #print featureblobs
        # start at biggest blob size
        for featureblob in reversed(featureblobs):
            # for every blob found of a blobsize

            for blob in enumerate(featureblob):
                # if that blob is not within range of another blob, add it
                blob = blob[1]
                if not within_range(blob, blob_coords):
                    blob_coords.append([z, blob[0], blob[1], blob[2]])
        list.append(blob_coords)
    return list 

def image_to_spots(
            self, data_image: np.ndarray,
    ) -> PerImageSliceSpotResults:
        """
        Find spots using a gaussian blob finding algorithm

        Parameters
        ----------
        data_image : np.ndarray
            image containing spots to be detected

        Returns
        -------
        PerImageSpotResults :
            includes a SpotAttributes DataFrame of metadata containing the coordinates, intensity
            and radius of each spot, as well as any extra information collected during spot finding.

        """

        spot_finding_args = {
            "min_sigma": self.min_sigma,
            "max_sigma": self.max_sigma,
            "threshold": self.threshold,
            "exclude_border": self.exclude_border,
            "overlap": self.overlap,
            "num_sigma": self.num_sigma
        }
        if self.detector_method is not blob_doh:
            spot_finding_args["exclude_border"] = self.exclude_border

        fitted_blobs_array: np.ndarray = self.detector_method(
            data_image,
            **spot_finding_args
        )

        if fitted_blobs_array.shape[0] == 0:
            empty_spot_attrs = SpotAttributes.empty(
                extra_fields=[Features.INTENSITY, Features.SPOT_ID])
            return PerImageSliceSpotResults(spot_attrs=empty_spot_attrs, extras=None)

        # measure intensities
        z_inds = fitted_blobs_array[:, 0].astype(int)
        y_inds = fitted_blobs_array[:, 1].astype(int)
        x_inds = fitted_blobs_array[:, 2].astype(int)
        radius = np.round(fitted_blobs_array[:, 3] * np.sqrt(3))
        data_image = np.asarray(data_image)
        intensities = data_image[tuple([z_inds, y_inds, x_inds])]

        # construct dataframe
        spot_data = pd.DataFrame(
            data={
                Features.INTENSITY: intensities,
                Axes.ZPLANE.value: z_inds,
                Axes.Y.value: y_inds,
                Axes.X.value: x_inds,
                Features.SPOT_RADIUS: radius,
            }
        )
        spots = SpotAttributes(spot_data)
        spots.data[Features.SPOT_ID] = np.arange(spots.data.shape[0])
        return PerImageSliceSpotResults(spot_attrs=spots, extras=None) 

def run(
            self,
            image_stack: ImageStack,
            reference_image: Optional[ImageStack] = None,
            n_processes: Optional[int] = None,
            *args,
    ) -> SpotFindingResults:
        """
        Find spots in the given ImageStack using a gaussian blob finding algorithm.
        If a reference image is provided the spots will be detected there then measured
        across all rounds and channels in the corresponding ImageStack. If a reference_image
        is not provided spots will be detected _independently_ in each channel. This assumes
        a non-multiplex imaging experiment, as only one (ch, round) will be measured for each spot.

        Parameters
        ----------
        image_stack : ImageStack
            ImageStack where we find the spots in.
        reference_image : Optional[ImageStack]
            (Optional) a reference image. If provided, spots will be found in this image, and then
            the locations that correspond to these spots will be measured across each channel.
        n_processes : Optional[int] = None,
            Number of processes to devote to spot finding.
        """
        spot_finding_method = partial(self.image_to_spots, *args)
        if reference_image:
            data_image = reference_image._squeezed_numpy(*{Axes.ROUND, Axes.CH})
            if self.detector_method is blob_doh and data_image.ndim > 2:
                raise ValueError("blob_doh only support 2d images")
            reference_spots = spot_finding_method(data_image)
            results = spot_finding_utils.measure_intensities_at_spot_locations_across_imagestack(
                data_image=image_stack,
                reference_spots=reference_spots,
                measurement_function=self.measurement_function)
        else:
            if self.detector_method is blob_doh and self.is_volume:
                raise ValueError("blob_doh only support 2d images")
            spot_attributes_list = image_stack.transform(
                func=spot_finding_method,
                group_by=determine_axes_to_group_by(self.is_volume),
                n_processes=n_processes
            )
            results = SpotFindingResults(imagestack_coords=image_stack.xarray.coords,
                                         log=image_stack.log,
                                         spot_attributes_list=spot_attributes_list)
        return results 

def blob_detection(
    input_image,
    method="log",
    threshold=0.5,
    min_sigma=3,
    max_sigma=20,
    overlap=0.5,
    return_sigmas=False,
    **kwargs,
):
    """
    Interface to the `feature.blob_*`_ methods implemented in scikit-image. A
    blob is defined as local a maximum of a Gaussian-filtered image.

    .. _`feature.blob_*`:\
    https://scikit-image.org/docs/dev/auto_examples/features_detection/plot_blob.html

    Parameters
    ----------
    input_image : array_like
        Array of shape (m, n) containing the input image. Nan values are ignored.
    method : {'log', 'dog', 'doh'}, optional
        The method to use: 'log' = Laplacian of Gaussian, 'dog' = Difference of
        Gaussian, ''
    threshold : float, optional
        Detection threshold.
    min_sigma : float, optional
        The minimum standard deviation for the Gaussian kernel.
    max_sigma : float, optional
        The maximum standard deviation for the Gaussian kernel.
    overlap : float, optional
        A value between 0 and 1. If the area of two blobs overlaps by a fraction
        greater than threshold, the smaller blob is eliminated.
    return_sigmas : bool, optional
        If True, the return array has a third column indicating the standard
        deviations of the Gaussian kernels that detected the blobs.
    """
    if method not in ["log", "dog", "doh"]:
        raise ValueError("unknown method %s, must be 'log', 'dog' or 'doh'" % method)

    if method == "log":
        detector = feature.blob_log
    elif method == "dog":
        detector = feature.blob_dog
    else:
        detector = feature.blob_doh

    blobs = detector(
        input_image,
        min_sigma=min_sigma,
        max_sigma=max_sigma,
        threshold=threshold,
        overlap=overlap,
        **kwargs,
    )

    if not return_sigmas:
        blobs = blobs[:, :2]

    return np.column_stack([blobs[:, 1], blobs[:, 0]])