{"cells":[{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"802e2ae5-a614-639d-384c-1cf2bba7be72"},"outputs":[],"source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nfrom subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n\n# Any results you write to the current directory are saved as output."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b4ce0b51-f359-7cc3-2185-b4f851f975c3"},"outputs":[],"source":"\nimport numpy as np\nfrom scipy.ndimage import gaussian_filter, gaussian_laplace\nimport itertools as itt\nimport math\nfrom math import sqrt, hypot, log\nfrom numpy import arccos\nfrom ..util import img_as_float\nfrom .peak import peak_local_max\nfrom ._hessian_det_appx import _hessian_matrix_det\nfrom ..transform import integral_image\nfrom .._shared.utils import assert_nD\n\n\n# This basic blob detection algorithm is based on:\n# http://www.cs.utah.edu/~jfishbau/advimproc/project1/ (04.04.2013)\n# Theory behind: http://en.wikipedia.org/wiki/Blob_detection (04.04.2013)\n\n\ndef _blob_overlap(blob1, blob2):\n    \"\"\"Finds the overlapping area fraction between two blobs.\n\n    Returns a float representing fraction of overlapped area.\n\n    Parameters\n    ----------\n    blob1 : sequence\n        A sequence of ``(y,x,sigma)``, where ``x,y`` are coordinates of blob\n        and sigma is the standard deviation of the Gaussian kernel which\n        detected the blob.\n    blob2 : sequence\n        A sequence of ``(y,x,sigma)``, where ``x,y`` are coordinates of blob\n        and sigma is the standard deviation of the Gaussian kernel which\n        detected the blob.\n\n    Returns\n    -------\n    f : float\n        Fraction of overlapped area.\n    \"\"\"\n    root2 = sqrt(2)\n\n    # extent of the blob is given by sqrt(2)*scale\n    r1 = blob1[2] * root2\n    r2 = blob2[2] * root2\n\n    d = hypot(blob1[0] - blob2[0], blob1[1] - blob2[1])\n\n    if d > r1 + r2:\n        return 0\n\n    # one blob is inside the other, the smaller blob must die\n    if d <= abs(r1 - r2):\n        return 1\n\n    ratio1 = (d ** 2 + r1 ** 2 - r2 ** 2) / (2 * d * r1)\n    ratio1 = np.clip(ratio1, -1, 1)\n    acos1 = arccos(ratio1)\n\n    ratio2 = (d ** 2 + r2 ** 2 - r1 ** 2) / (2 * d * r2)\n    ratio2 = np.clip(ratio2, -1, 1)\n    acos2 = arccos(ratio2)\n\n    a = -d + r2 + r1\n    b = d - r2 + r1\n    c = d + r2 - r1\n    d = d + r2 + r1\n    area = r1 ** 2 * acos1 + r2 ** 2 * acos2 - 0.5 * sqrt(abs(a * b * c * d))\n\n    return area / (math.pi * (min(r1, r2) ** 2))\n\n\ndef _prune_blobs(blobs_array, overlap):\n    \"\"\"Eliminated blobs with area overlap.\n\n    Parameters\n    ----------\n    blobs_array : ndarray\n        A 2d array with each row representing 3 values, ``(y,x,sigma)``\n        where ``(y,x)`` are coordinates of the blob and ``sigma`` is the\n        standard deviation of the Gaussian kernel which detected the blob.\n    overlap : float\n        A value between 0 and 1. If the fraction of area overlapping for 2\n        blobs is greater than `overlap` the smaller blob is eliminated.\n\n    Returns\n    -------\n    A : ndarray\n        `array` with overlapping blobs removed.\n    \"\"\"\n\n    # iterating again might eliminate more blobs, but one iteration suffices\n    # for most cases\n    for blob1, blob2 in itt.combinations(blobs_array, 2):\n        if _blob_overlap(blob1, blob2) > overlap:\n            if blob1[2] > blob2[2]:\n                blob2[2] = -1\n            else:\n                blob1[2] = -1\n\n    # return blobs_array[blobs_array[:, 2] > 0]\n    return np.array([b for b in blobs_array if b[2] > 0])\n\n\ndef blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=2.0,\n             overlap=.5,):\n    \"\"\"Finds blobs in the given grayscale image.\n\n    Blobs are found using the Difference of Gaussian (DoG) method [1]_.\n    For each blob found, the method returns its coordinates and the standard\n    deviation of the Gaussian kernel that detected the blob.\n\n    Parameters\n    ----------\n    image : ndarray\n        Input grayscale image, blobs are assumed to be light on dark\n        background (white on black).\n    min_sigma : float, optional\n        The minimum standard deviation for Gaussian Kernel. Keep this low to\n        detect smaller blobs.\n    max_sigma : float, optional\n        The maximum standard deviation for Gaussian Kernel. Keep this high to\n        detect larger blobs.\n    sigma_ratio : float, optional\n        The ratio between the standard deviation of Gaussian Kernels used for\n        computing the Difference of Gaussians\n    threshold : float, optional.\n        The absolute lower bound for scale space maxima. Local maxima smaller\n        than thresh are ignored. Reduce this to detect blobs with less\n        intensities.\n    overlap : float, optional\n        A value between 0 and 1. If the area of two blobs overlaps by a\n        fraction greater than `threshold`, the smaller blob is eliminated.\n\n    Returns\n    -------\n    A : (n, 3) ndarray\n        A 2d array with each row representing 3 values, ``(y,x,sigma)``\n        where ``(y,x)`` are coordinates of the blob and ``sigma`` is the\n        standard deviation of the Gaussian kernel which detected the blob.\n\n    References\n    ----------\n    .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach\n\n    Examples\n    --------\n    >>> from skimage import data, feature\n    >>> feature.blob_dog(data.coins(), threshold=.5, max_sigma=40)\n    array([[  45.      ,  336.      ,   16.777216],\n           [  52.      ,  155.      ,   16.777216],\n           [  52.      ,  216.      ,   16.777216],\n           [  54.      ,   42.      ,   16.777216],\n           [  54.      ,  276.      ,   10.48576 ],\n           [  58.      ,  100.      ,   10.48576 ],\n           [ 120.      ,  272.      ,   16.777216],\n           [ 124.      ,  337.      ,   10.48576 ],\n           [ 125.      ,   45.      ,   16.777216],\n           [ 125.      ,  208.      ,   10.48576 ],\n           [ 127.      ,  102.      ,   10.48576 ],\n           [ 128.      ,  154.      ,   10.48576 ],\n           [ 185.      ,  347.      ,   16.777216],\n           [ 193.      ,  213.      ,   16.777216],\n           [ 194.      ,  277.      ,   16.777216],\n           [ 195.      ,  102.      ,   16.777216],\n           [ 196.      ,   43.      ,   10.48576 ],\n           [ 198.      ,  155.      ,   10.48576 ],\n           [ 260.      ,   46.      ,   16.777216],\n           [ 261.      ,  173.      ,   16.777216],\n           [ 263.      ,  245.      ,   16.777216],\n           [ 263.      ,  302.      ,   16.777216],\n           [ 267.      ,  115.      ,   10.48576 ],\n           [ 267.      ,  359.      ,   16.777216]])\n\n    Notes\n    -----\n    The radius of each blob is approximately :math:`\\sqrt{2}sigma`.\n    \"\"\"\n    assert_nD(image, 2)\n\n    image = img_as_float(image)\n\n    # k such that min_sigma*(sigma_ratio**k) > max_sigma\n    k = int(log(float(max_sigma) / min_sigma, sigma_ratio)) + 1\n\n    # a geometric progression of standard deviations for gaussian kernels\n    sigma_list = np.array([min_sigma * (sigma_ratio ** i)\n                           for i in range(k + 1)])\n\n    gaussian_images = [gaussian_filter(image, s) for s in sigma_list]\n\n    # computing difference between two successive Gaussian blurred images\n    # multiplying with standard deviation provides scale invariance\n    dog_images = [(gaussian_images[i] - gaussian_images[i + 1])\n                  * sigma_list[i] for i in range(k)]\n    image_cube = np.dstack(dog_images)\n\n    # local_maxima = get_local_maxima(image_cube, threshold)\n    local_maxima = peak_local_max(image_cube, threshold_abs=threshold,\n                                  footprint=np.ones((3, 3, 3)),\n                                  threshold_rel=0.0,\n                                  exclude_border=False)\n    # Convert local_maxima to float64\n    lm = local_maxima.astype(np.float64)\n    # Convert the last index to its corresponding scale value\n    lm[:, 2] = sigma_list[local_maxima[:, 2]]\n    local_maxima = lm\n    return _prune_blobs(local_maxima, overlap)\n\n\ndef blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,\n             overlap=.5, log_scale=False):\n    \"\"\"Finds blobs in the given grayscale image.\n\n    Blobs are found using the Laplacian of Gaussian (LoG) method [1]_.\n    For each blob found, the method returns its coordinates and the standard\n    deviation of the Gaussian kernel that detected the blob.\n\n    Parameters\n    ----------\n    image : ndarray\n        Input grayscale image, blobs are assumed to be light on dark\n        background (white on black).\n    min_sigma : float, optional\n        The minimum standard deviation for Gaussian Kernel. Keep this low to\n        detect smaller blobs.\n    max_sigma : float, optional\n        The maximum standard deviation for Gaussian Kernel. Keep this high to\n        detect larger blobs.\n    num_sigma : int, optional\n        The number of intermediate values of standard deviations to consider\n        between `min_sigma` and `max_sigma`.\n    threshold : float, optional.\n        The absolute lower bound for scale space maxima. Local maxima smaller\n        than thresh are ignored. Reduce this to detect blobs with less\n        intensities.\n    overlap : float, optional\n        A value between 0 and 1. If the area of two blobs overlaps by a\n        fraction greater than `threshold`, the smaller blob is eliminated.\n    log_scale : bool, optional\n        If set intermediate values of standard deviations are interpolated\n        using a logarithmic scale to the base `10`. If not, linear\n        interpolation is used.\n\n    Returns\n    -------\n    A : (n, 3) ndarray\n        A 2d array with each row representing 3 values, ``(y,x,sigma)``\n        where ``(y,x)`` are coordinates of the blob and ``sigma`` is the\n        standard deviation of the Gaussian kernel which detected the blob.\n\n    References\n    ----------\n    .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_Laplacian_of_Gaussian\n\n    Examples\n    --------\n    >>> from skimage import data, feature, exposure\n    >>> img = data.coins()\n    >>> img = exposure.equalize_hist(img)  # improves detection\n    >>> feature.blob_log(img, threshold = .3)\n    array([[ 113.        ,  323.        ,    1.        ],\n           [ 121.        ,  272.        ,   17.33333333],\n           [ 124.        ,  336.        ,   11.88888889],\n           [ 126.        ,   46.        ,   11.88888889],\n           [ 126.        ,  208.        ,   11.88888889],\n           [ 127.        ,  102.        ,   11.88888889],\n           [ 128.        ,  154.        ,   11.88888889],\n           [ 185.        ,  344.        ,   17.33333333],\n           [ 194.        ,  213.        ,   17.33333333],\n           [ 194.        ,  276.        ,   17.33333333],\n           [ 197.        ,   44.        ,   11.88888889],\n           [ 198.        ,  103.        ,   11.88888889],\n           [ 198.        ,  155.        ,   11.88888889],\n           [ 260.        ,  174.        ,   17.33333333],\n           [ 263.        ,  244.        ,   17.33333333],\n           [ 263.        ,  302.        ,   17.33333333],\n           [ 266.        ,  115.        ,   11.88888889]])\n\n    Notes\n    -----\n    The radius of each blob is approximately :math:`\\sqrt{2}sigma`.\n    \"\"\"\n\n    assert_nD(image, 2)\n\n    image = img_as_float(image)\n\n    if log_scale:\n        start, stop = log(min_sigma, 10), log(max_sigma, 10)\n        sigma_list = np.logspace(start, stop, num_sigma)\n    else:\n        sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)\n\n    # computing gaussian laplace\n    # s**2 provides scale invariance\n    gl_images = [-gaussian_laplace(image, s) * s ** 2 for s in sigma_list]\n    image_cube = np.dstack(gl_images)\n\n    local_maxima = peak_local_max(image_cube, threshold_abs=threshold,\n                                  footprint=np.ones((3, 3, 3)),\n                                  threshold_rel=0.0,\n                                  exclude_border=False)\n\n    # Convert local_maxima to float64\n    lm = local_maxima.astype(np.float64)\n    # Convert the last index to its corresponding scale value\n    lm[:, 2] = sigma_list[local_maxima[:, 2]]\n    local_maxima = lm\n    return _prune_blobs(local_maxima, overlap)\n\n\ndef blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,\n             overlap=.5, log_scale=False):\n    \"\"\"Finds blobs in the given grayscale image.\n\n    Blobs are found using the Determinant of Hessian method [1]_. For each blob\n    found, the method returns its coordinates and the standard deviation\n    of the Gaussian Kernel used for the Hessian matrix whose determinant\n    detected the blob. Determinant of Hessians is approximated using [2]_.\n\n    Parameters\n    ----------\n    image : ndarray\n        Input grayscale image.Blobs can either be light on dark or vice versa.\n    min_sigma : float, optional\n        The minimum standard deviation for Gaussian Kernel used to compute\n        Hessian matrix. Keep this low to detect smaller blobs.\n    max_sigma : float, optional\n        The maximum standard deviation for Gaussian Kernel used to compute\n        Hessian matrix. Keep this high to detect larger blobs.\n    num_sigma : int, optional\n        The number of intermediate values of standard deviations to consider\n        between `min_sigma` and `max_sigma`.\n    threshold : float, optional.\n        The absolute lower bound for scale space maxima. Local maxima smaller\n        than thresh are ignored. Reduce this to detect less prominent blobs.\n    overlap : float, optional\n        A value between 0 and 1. If the area of two blobs overlaps by a\n        fraction greater than `threshold`, the smaller blob is eliminated.\n    log_scale : bool, optional\n        If set intermediate values of standard deviations are interpolated\n        using a logarithmic scale to the base `10`. If not, linear\n        interpolation is used.\n\n    Returns\n    -------\n    A : (n, 3) ndarray\n        A 2d array with each row representing 3 values, ``(y,x,sigma)``\n        where ``(y,x)`` are coordinates of the blob and ``sigma`` is the\n        standard deviation of the Gaussian kernel of the Hessian Matrix whose\n        determinant detected the blob.\n\n    References\n    ----------\n    .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian\n\n    .. [2] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool,\n           \"SURF: Speeded Up Robust Features\"\n           ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf\n\n    Examples\n    --------\n    >>> from skimage import data, feature\n    >>> img = data.coins()\n    >>> feature.blob_doh(img)\n    array([[ 121.        ,  271.        ,   30.        ],\n           [ 123.        ,   44.        ,   23.55555556],\n           [ 123.        ,  205.        ,   20.33333333],\n           [ 124.        ,  336.        ,   20.33333333],\n           [ 126.        ,  101.        ,   20.33333333],\n           [ 126.        ,  153.        ,   20.33333333],\n           [ 156.        ,  302.        ,   30.        ],\n           [ 185.        ,  348.        ,   30.        ],\n           [ 192.        ,  212.        ,   23.55555556],\n           [ 193.        ,  275.        ,   23.55555556],\n           [ 195.        ,  100.        ,   23.55555556],\n           [ 197.        ,   44.        ,   20.33333333],\n           [ 197.        ,  153.        ,   20.33333333],\n           [ 260.        ,  173.        ,   30.        ],\n           [ 262.        ,  243.        ,   23.55555556],\n           [ 265.        ,  113.        ,   23.55555556],\n           [ 270.        ,  363.        ,   30.        ]])\n\n    Notes\n    -----\n    The radius of each blob is approximately `sigma`.\n    Computation of Determinant of Hessians is independent of the standard\n    deviation. Therefore detecting larger blobs won't take more time. In\n    methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation\n    of Gaussians for larger `sigma` takes more time. The downside is that\n    this method can't be used for detecting blobs of radius less than `3px`\n    due to the box filters used in the approximation of Hessian Determinant.\n    \"\"\"\n\n    assert_nD(image, 2)\n\n    image = img_as_float(image)\n    image = integral_image(image)\n\n    if log_scale:\n        start, stop = log(min_sigma, 10), log(max_sigma, 10)\n        sigma_list = np.logspace(start, stop, num_sigma)\n    else:\n        sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)\n\n    hessian_images = [_hessian_matrix_det(image, s) for s in sigma_list]\n    image_cube = np.dstack(hessian_images)\n\n    local_maxima = peak_local_max(image_cube, threshold_abs=threshold,\n                                  footprint=np.ones((3, 3, 3)),\n                                  threshold_rel=0.0,\n                                  exclude_border=False)\n\n    # Convert local_maxima to float64\n    lm = local_maxima.astype(np.float64)\n    # Convert the last index to its corresponding scale value\n    lm[:, 2] = sigma_list[local_maxima[:, 2]]\n    local_maxima = lm\n    return _prune_blobs(local_maxima, overlap)\n"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d126101f-d93d-af73-a3b8-69c7f07ec213"},"outputs":[],"source":""}],"metadata":{"_change_revision":0,"_is_fork":false,"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.0"}},"nbformat":4,"nbformat_minor":0}