{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"daecbdb6-470a-3efb-c2ec-28a2c399e860"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"6e16c8e3-de87-918e-b59d-b8643256ae5a"},"outputs":[],"source":"import pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib.mlab as mlab\nimport matplotlib.image as mpimg\nimport numpy as np\nfrom PIL import Image\nfrom scipy import sparse\nfrom scipy import ndimage\nfrom scipy.ndimage import gaussian_filter\nfrom skimage import data\nfrom skimage import img_as_float\nfrom skimage import morphology, measure\nfrom skimage.color import label2rgb\n\n%matplotlib inline"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"7d9aae62-d2a1-50e4-6541-3107a18009fd"},"outputs":[],"source":"raw_image = mpimg.imread('../input/Train/0.jpg')\ndot_image = mpimg.imread('../input/TrainDotted/0.jpg')\n\n# Convert to floats. Will save us headache later.\nraw_image = raw_image.astype(float)\ndot_image = dot_image.astype(float)\nraw_image = raw_image / raw_image.max()\ndot_image = dot_image / dot_image.max()\n\n# Spotcheck the data.\nplt.figure(figsize=(20,20))\nplt.subplot(1, 3, 1)\nimgplot = plt.imshow(raw_image[:, :, 0], cmap=plt.cm.gray)\nplt.subplot(1, 3, 2)\nimgplot = plt.imshow(raw_image[:, :, 1], cmap=plt.cm.gray)\nplt.subplot(1, 3, 3)\nimgplot = plt.imshow(raw_image[:, :, 2], cmap=plt.cm.gray)\nplt.show()\n\ndef img_to_float(m):\n    m = m.astype(float)\n    return m / m.max()\n\ndef float_to_img(m):\n    return (m * 255.).astype(int)\n\ndef spot_check(img, show_pil=False):\n    imgplot = plt.imshow(img)\n    plt.show()\n\n    if show_pil:\n        # We can use this to see a higher resolution image, useful for checking things.\n        pil_image = Image.fromarray((img*255.0).astype('uint8'))\n        pil_image.show()"},{"cell_type":"markdown","metadata":{"_cell_guid":"2c3db247-c9c1-d85c-3cb7-1e0b84b0b883"},"source":""},{"cell_type":"markdown","metadata":{"_cell_guid":"08d19185-ef84-843e-a811-1285b4e14670"},"source":""},{"cell_type":"markdown","metadata":{"_cell_guid":"a02bf0b5-4cc8-958a-7894-fdc4e5d2b373"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"73e59ef0-9fd8-b8a3-e4a5-3e91829e9488"},"outputs":[],"source":"# Dot images have some black artifacts. Let's mask those out.\ndot_norm = np.linalg.norm(dot_image, axis=2)\nthreshold = (dot_norm.max() - dot_norm.min()) * 0.02\ninitial_mask_1d = dot_norm <= threshold\n\n# Broadcast to 3d for true image mask.\n_, initial_mask = np.broadcast_arrays(dot_image, initial_mask_1d[..., None])\n\nplt.imshow(initial_mask.astype(float))\nplt.show()"},{"cell_type":"markdown","metadata":{"_cell_guid":"5435aeb7-28e7-337e-90c6-3dc20c3daa0c"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"60dc3d66-3411-9670-dcb6-8a6cbd628f7e"},"outputs":[],"source":"##### Remove the background. Mask with the original image.\ndot_diff = np.abs(dot_image - raw_image)\ndot_diff[initial_mask] = 0\n#dot_diff = np.linalg.norm(dot_diff, axis=2)\ndot_diff = np.max(dot_diff, axis=2)\n\n# Remove jpeg artifact noise.\ndot_diff[dot_diff < 0.3] = 0\n\n# Other ways to consider removing noise.\n#dot_diff = ndimage.binary_opening(dot_diff)\n#dot_diff = ndimage.binary_closing(dot_diff)\n\nplt.figure(figsize=(20,20))\nplt.subplot(1, 2, 1)\nimgplot = plt.imshow(dot_diff)\n\ndef mask_image(m, mask, c):\n    mask_1d = mask < 0.01\n    _, mask_3d = np.broadcast_arrays(m, mask_1d[..., None])\n    m[mask_3d] = c\n\ndots = dot_image.copy()\nmask_image(dots, dot_diff, 0)\nplt.subplot(1, 2, 2)\nimgplot = plt.imshow(dots)\nplt.show()\n\ndef channel_histogram(m, idx, color):\n    c = (m[:, :, idx] * 255).astype('uint8')\n    c = c[c > 0]\n    n, bins, patches = plt.hist(c.ravel(), bins=256, facecolor=color)\n    plt.plot(bins)\n\ndef print_histograms(m):\n    plt.figure(figsize=(20,5))\n    plt.subplot(1, 3, 1)\n    channel_histogram(m, 0, 'red')\n    plt.subplot(1, 3, 2)\n    channel_histogram(m, 1, 'blue')\n    plt.subplot(1, 3, 3)\n    channel_histogram(m, 2, 'green')\n    plt.show()\n\nprint_histograms(dots)\n"},{"cell_type":"markdown","metadata":{"_cell_guid":"2dc186b4-81bb-1920-5c03-d3beb7214b20"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"3f620fb1-09b7-e073-5253-7b96152f7683"},"outputs":[],"source":"# Other ways to consider removing noise.\ndot_diff_eroded = ndimage.binary_erosion(dot_diff)\n#dot_diff_eroded = ndimage.binary_erosion(dot_diff_eroded)\ndots = dot_image.copy()\nmask_image(dots, dot_diff_eroded, 0)\n\nspot_check(dots, False)\n\nprint_histograms(dots)"},{"cell_type":"markdown","metadata":{"_cell_guid":"a552e19e-2efc-76ce-dc27-6f57758f9ea6"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"7552cc0d-494d-7d82-5bb6-df36a84a4be7"},"outputs":[],"source":"labeled_dots, label_count = morphology.label(dot_diff_eroded, return_num=True, connectivity=2)\nrgb_labels = label2rgb(labeled_dots, raw_image)\nimgplot = plt.imshow(rgb_labels)\nprint(\"Num of labels = {0}\".format(label_count))"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"54809ecb-5e6b-7f74-e208-3f93cddeb855"},"outputs":[],"source":"from skimage import measure\nplt.figure(figsize=(20,20))\nplt.imshow(dot_image)\naxis = plt.gcf().gca()\nprops = measure.regionprops(labeled_dots)\nprint(len(props))\nfor region in props:\n    center = region.centroid\n    c = plt.Circle((center[1], center[0]), 30, color='y', fill=True, clip_on=False, alpha=0.5)\n    axis.add_artist(c)\nplt.show()"},{"cell_type":"markdown","metadata":{"_cell_guid":"0e5181f9-e77c-4929-3d1f-1983f3917bb5"},"source":""},{"cell_type":"markdown","metadata":{"_cell_guid":"fbcb7624-bf5f-1aa7-631c-3d10587ae963"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"fd135510-1a4b-7804-369b-4d2fcab41167"},"outputs":[],"source":"def anisotropic_diffusion(img, niter=1, kappa=50, gamma=0.1, voxelspacing=None, option=1):\n    import numpy\n    r\"\"\"\n    Edge-preserving, XD Anisotropic diffusion.\n\n\n    Parameters\n    ----------\n    img : array_like\n        Input image (will be cast to numpy.float).\n    niter : integer\n        Number of iterations.\n    kappa : integer\n        Conduction coefficient, e.g. 20-100. ``kappa`` controls conduction\n        as a function of the gradient. If ``kappa`` is low small intensity\n        gradients are able to block conduction and hence diffusion across\n        steep edges. A large value reduces the influence of intensity gradients\n        on conduction.\n    gamma : float\n        Controls the speed of diffusion. Pick a value :math:`<= .25` for stability.\n    voxelspacing : tuple of floats\n        The distance between adjacent pixels in all img.ndim directions\n    option : {1, 2}\n        Whether to use the Perona Malik diffusion equation No. 1 or No. 2.\n        Equation 1 favours high contrast edges over low contrast ones, while\n        equation 2 favours wide regions over smaller ones. See [1]_ for details.\n\n    Returns\n    -------\n    anisotropic_diffusion : ndarray\n        Diffused image.\n        \n    Notes\n    -----\n    Original MATLAB code by Peter Kovesi,\n    School of Computer Science & Software Engineering,\n    The University of Western Australia,\n    pk @ csse uwa edu au,\n    <http://www.csse.uwa.edu.au>\n\n    Translated to Python and optimised by Alistair Muldal,\n    Department of Pharmacology,\n    University of Oxford,\n    <alistair.muldal@pharm.ox.ac.uk>\n    \n    Adapted to arbitrary dimensionality and added to the MedPy library Oskar Maier,\n    Institute for Medical Informatics,\n    Universitaet Luebeck,\n    <oskar.maier@googlemail.com>\n    \n    June 2000  original version. -\n    March 2002 corrected diffusion eqn No 2. - \n    July 2012 translated to Python -\n    August 2013 incorporated into MedPy, arbitrary dimensionality -    \n\n    References\n    ----------\n    .. [1] P. Perona and J. Malik. \n       Scale-space and edge detection using ansotropic diffusion.\n       IEEE Transactions on Pattern Analysis and Machine Intelligence, \n       12(7):629-639, July 1990.\n    \"\"\"\n    # define conduction gradients functions\n    if option == 1:\n        def condgradient(delta, spacing):\n            return numpy.exp(-(delta/kappa)**2.)/float(spacing)\n    elif option == 2:\n        def condgradient(delta, spacing):\n            return 1./(1.+(delta/kappa)**2.)/float(spacing)\n\n    # initialize output array\n    out = numpy.array(img, dtype=numpy.float32, copy=True)\n    \n    # set default voxel spacong if not suppliec\n    if None == voxelspacing:\n        voxelspacing = tuple([1.] * img.ndim)\n\n    # initialize some internal variables\n    deltas = [numpy.zeros_like(out) for _ in range(out.ndim)]\n\n    for _ in range(niter):\n\n        # calculate the diffs\n        for i in range(out.ndim):\n            slicer = [slice(None, -1) if j == i else slice(None) for j in range(out.ndim)]\n            deltas[i][slicer] = numpy.diff(out, axis=i)\n        \n        # update matrices\n        matrices = [condgradient(delta, spacing) * delta for delta, spacing in zip(deltas, voxelspacing)]\n\n        # subtract a copy that has been shifted ('Up/North/West' in 3D case) by one\n        # pixel. Don't as questions. just do it. trust me.\n        for i in range(out.ndim):\n            slicer = [slice(1, None) if j == i else slice(None) for j in range(out.ndim)]\n            matrices[i][slicer] = numpy.diff(matrices[i], axis=i)\n\n        # update the image\n        out += gamma * (numpy.sum(matrices, axis=0))\n\n    return out"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"70457820-e795-8b4f-be78-de9711a384d4"},"outputs":[],"source":"img = raw_image.astype(float)\ndc = np.abs(img[:, :, 0] - img[:, :, 2])\naniso = anisotropic_diffusion(dc, 10, 200, .01)\n#dc[dc < 0] = 0\nmod = aniso\nplt.figure(figsize=(10,10))\nimgplot = plt.imshow(dc, cmap=plt.cm.gray)\nplt.show()\nplt.figure(figsize=(10,10))\nimgplot = plt.imshow(aniso, cmap=plt.cm.gray)\nplt.show()\nplt.figure(figsize=(10,10))\nimgplot = plt.imshow(aniso > .3, cmap=plt.cm.gray)\nplt.show()"},{"cell_type":"markdown","metadata":{"_cell_guid":"603425ea-4026-1a41-7950-85730e5e066e"},"source":"## Thoughts\nOn the first image this seems like a reasonable approach. On the second image we aren't getting nearly the same value out of such a simple calculation. I do think I'm going to go back to the first image and see if I can get some type of prediction. But for now I think I'll move on and start brainstorming more complex solutions."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"f5f055ee-637c-524a-8b74-f9cda8b1b0ee"},"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}