{"cells":[{"metadata":{"_uuid":"1ad8ea4d5a8e2618b42457bd05b2710c2a6ba16e"},"cell_type":"markdown","source":"**README**\n\nWhat follows is a background extraction script for Data Science Bowl 2018 (identification of cell nuclei in images) based on simple image analysis (thresholding).\n\nThe first part, defining the algorithm used, is followed by a few examples. It may be helpful to simply execute the definition cells first and jump to the examples -- I tried to comment and visualize the procedure where relevant, so it should give one a good feel of what is happening.\n\nNote: the background extraction is still far from perfect (see the last example), but I will be very glad if you find something useful for you in the code. Enjoy! :)"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# useful when you move it to your machine:\npath_stem = '../input/'\n\nverbose = True\nuse_agg = False # turn on when run without a graphics head\n\nMINIMAL_BLURRING_KERNEL_SIZE = 4.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"bb7436f0d00da9c8fc33938f63601a1e4a08f679"},"cell_type":"code","source":"# ----------------------------------------------------------------------------\n# IMPORTING\n\n# import train data\ndef import_data (folder='stage1_train'):\n    print('Importing train data...')\n    from tqdm import tqdm\n    import os\n\n    DATA = {}; DATA['images'] = {};# DATA['masks'] = {}\n    for dataset in tqdm(os.listdir(path_stem+folder)):\n        for typ in ['images']: #, 'masks']:\n            DATA[typ][dataset] = os.listdir(path_stem+folder+'/'+dataset+'/'+typ)\n        # report if anything unusual found\n        if len(DATA['images'][dataset]) != 1:\n            print('   %i images for %s' % (len(DATA['images']),dataset))\n        #if len(DATA['masks'][dataset]) < 1:\n        #    print('   no masks for %s' % dataset)\n    print('done. Imported %i datasets' % len(DATA['images']))\n\n    return DATA\n\n# image normalization functions\ndef normalize_color (img):\n    import numpy as np\n    maxs = np.max(img, axis=(0,1))\n    img = np.array(img) / max(maxs)\n    return img\ndef normalize_gray (img):\n    import numpy as np\n    maxx = np.max(img)\n    img = np.array(img) / maxx\n    return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c8c0a8f816fd90de66cf02eef5be4b43e51d5a4f"},"cell_type":"code","source":"# ----------------------------------------------------------------------------\n# BACKGROUND EXTRACTION\n\ndef detect_background (dataset, imgname, visualize=False, folder='stage1_train', n_channels=128):\n\n    import matplotlib\n    if use_agg:\n        matplotlib.use('Agg')\n    import matplotlib.pyplot as plt\n    import numpy as np\n    import cv2\n\n    from scipy.optimize import curve_fit, bisect\n\n    # definitions need to be copied here for multithreading to work\n    def normalize_gray (img):\n        import numpy as np\n        maxx = np.max(img)\n        img = np.array(img) / maxx\n        return img\n    def gaussian (x, A, x0, log_sigma):\n        sigma = np.exp(log_sigma)\n        return A * np.exp(-(x-x0)**2/(2.0*sigma**2))\n    def lorentzian (x, A, x0, log_sigma):\n        sigma = np.exp(log_sigma)\n        return A * sigma / ( (x-x0)**2 + sigma**2)\n\n    try:\n        if verbose:\n            print('Reading the image... ', end='', flush=True)\n        imgpath = path_stem+folder+'/'+dataset+'/images/'+imgname\n        img = cv2.imread(imgpath, 0)\n        img = normalize_gray(img)\n        if len(set(np.array(img).flatten())) == 1:\n            print(' Image empty, moving on.')\n            return -1, -1\n        if verbose:\n            print('done.', flush=True)\n        img0 = 1. * img\n        img_range = np.min(img), np.max(img)\n        n_bins = max(int(len(set(np.array(img).flatten()))/3.), 10)\n        force_gmm = False\n        force_gmm_brightest = False\n        blurring_kernel_size = MINIMAL_BLURRING_KERNEL_SIZE\n\n        for trial_idx in range(10): # redo if noisy image detected, see below\n\n            try:\n\n                if visualize:\n                    plt.clf()\n                    plt.imshow(img)\n                    plt.title('Current image state')\n                    plt.show()\n                    plt.clf()\n                vals_sparse, _, _ = plt.hist(np.array(img).flatten(), 32)\n                plt.clf()\n                vals, bins, _ = plt.hist(np.array(img).flatten(), n_bins)\n                bins = np.array(bins); vals = np.array(vals)\n                bins = 0.5*(np.roll(bins, -1)[:-1] + bins[:-1])\n                if visualize:\n                    plt.title('Pixel brightness histogram. Black vertical solid line: threshold.')\n                    plt.plot(bins, vals, linewidth=2, color='red')\n                    plt.xlim(0.,1.)\n\n                # manual fitting of a lorentzian (bg) and a gaussian (cells) performs better if the histogram is background-dominated:\n                if max(vals_sparse) > 0.35 * np.sum(vals_sparse) and not force_gmm:\n                    # fit the main maximum (usually: background)\n                    if verbose:\n                        print('Finding the maximum... ', end='', flush=True)\n                    posmax = bins[np.where(vals == max(vals))[0][0]]\n                    first_gauss, cov = curve_fit(lambda x, A, log_sigma : lorentzian(x, A, posmax, log_sigma), bins, vals, p0=[max(vals), np.log(0.1)], bounds=([0.,-np.inf], [np.inf, 0.]))\n                    A, log_sigma = first_gauss\n                    if verbose:\n                        print('done.', flush=True)\n                    #correct for position\n                    if verbose:\n                        print('Correcting for position... ', end='', flush=True)\n                    first_gauss, cov = curve_fit(lambda x, x0 : lorentzian(x, A, x0, log_sigma), bins, vals, p0=[posmax], bounds=img_range)\n                    posmax = first_gauss[0]\n                    first_gauss, cov = curve_fit(lambda x, A, log_sigma : lorentzian(x, A, posmax, log_sigma), bins, vals, p0=[max(vals), np.log(0.1)], bounds=([0.,-np.inf], [np.inf, np.inf]))\n                    A, log_sigma = first_gauss\n                    if verbose:\n                        print('done.', flush=True)\n                    #plot\n                    if visualize:\n                        plt.axvline(posmax, color='g', linewidth=2)\n                        plt.plot(bins, lorentzian(bins, A, posmax, log_sigma), color='g', linewidth=2)\n                    # save the fit\n                    first_gauss = [A, posmax, log_sigma]\n\n                    # fit the secondary maximum (usually: cells)\n                    if verbose:\n                        print('Fitting the secondary maximum... ', end='', flush=True)\n                    primary_vals = lorentzian(bins, A, posmax, log_sigma)\n                    primary_max_idx = np.argmax(primary_vals)\n                    sigma_idx = max(int(len(primary_vals) * np.exp(log_sigma)), 2)\n                    mask = max(0, primary_max_idx-2*sigma_idx), min(len(bins), primary_max_idx+2*sigma_idx)\n                    vals = vals / (1. + primary_vals)\n                    vals[mask[0]:mask[1]] = 0. #0.5*(vals[mask[0]]+vals[mask[1]])\n                    if visualize:\n                        plt.plot(bins, vals, color='k', linewidth=2)\n                    posmax_secondary = bins[np.argmax(vals)]\n                    print(posmax_secondary)\n                    if abs(posmax - posmax_secondary) < 0.1:\n                        if posmax > 0.5:\n                            posmax_secondary = 0.5 * posmax\n                        else:\n                            posmax_secondary = 1. - 0.5*(1.-posmax)\n                    print(\"max vals\", max(vals))\n                    second_gauss, cov = curve_fit(lambda x, A, log_sigma : gaussian(x, A, posmax_secondary, log_sigma), bins, vals, p0=[max(vals), np.log(0.4)], bounds=([5.,-np.inf], [max(vals), np.log(1.)]))\n                    A, log_sigma = second_gauss\n                    if verbose:\n                        print('done.', flush=True)\n                    #correct for position\n                    if verbose:\n                        print('Correcting for position... ', end='', flush=True)\n                    second_gauss, cov = curve_fit(lambda x, x0 : gaussian(x, A, x0, log_sigma), bins, vals, p0=[posmax_secondary], bounds=img_range)\n                    posmax_secondary = second_gauss[0]\n                    second_gauss, cov = curve_fit(lambda x, A, log_sigma : gaussian(x, A, posmax_secondary, log_sigma), bins, vals, p0=[max(vals), np.log(0.1)], bounds=([5.,-np.inf], [max(vals), np.log(1.)]))\n                    A, log_sigma = second_gauss\n                    if verbose:\n                        print('done.', flush=True)\n                    #plot\n                    if visualize:\n                        plt.axvline(posmax_secondary, color='orange', linewidth=2)\n                        plt.plot(bins, gaussian(bins, A, posmax_secondary, log_sigma), color='orange', linewidth=2)\n                    # save the fit\n                    second_gauss = [A, posmax_secondary, log_sigma]\n\n                    # find where the two fits are equal\n                    if verbose:\n                        print('Splitting... ', end='', flush=True)\n                    try:\n                        splitter = bisect(lambda x : lorentzian(x, *first_gauss) - gaussian(x, *second_gauss), first_gauss[1], second_gauss[1])\n                    except Exception as e:\n                        splitter = 0.5 * (first_gauss[1] + second_gauss[1])\n                        print('Bisection failed, defaulting to average of maxima.')\n                    if verbose:\n                        print('done.', flush=True)\n                    if visualize:\n                        plt.axvline(splitter, color='k', linewidth=2)\n                        plt.show()\n                else: # use a Gaussian mixture model if cells contribute significantly to the histogram\n                    from sklearn.mixture import GaussianMixture\n                    X = np.array(img).flatten().reshape((-1,1)) # samples\n                    models = [GaussianMixture(n).fit(X) for n in range(2,5)]\n                    AIC = [m.aic(X) for m in models]\n                    model = models[np.argmin(AIC)] # best model\n                    x_axis = np.linspace(0.,1.,32)\n                    p = model.predict_proba(x_axis.reshape((-1,1))).transpose()\n                    if visualize:\n                        for y in p:\n                            plt.plot(x_axis, y)\n                    splitters = []\n                    for y in p:\n                        pos = np.where(y > 0.5)[0]\n                        if len(pos) > 0:\n                            splitters.append(x_axis[max(np.where(y > 0.5)[0])])\n                        else:\n                            splitters.append(x_axis[np.argmax(y)])\n                    splitters = np.array(sorted(splitters))[:-1]\n                    if visualize:\n                        for splitter in splitters:\n                            plt.axvline(splitter)\n                    if not force_gmm_brightest:\n                        splitter = splitters[0]\n                    else:\n                        tail = 1\n                        splitter = splitters[-tail]\n                        while len(np.where(img0 > splitter)[0]) < 0.05 * np.product(img0.shape):\n                            tail += 1\n                            if tail == len(splitters):\n                                break\n                            splitter = splitters[-tail]\n\n\n                if verbose:\n                    print(\"Threshold: \", splitter)\n\n                if visualize:\n                    plt.show()\n                    plt.clf()\n                    plt.subplot(121)\n                    img_to_show = 1. * img\n                    mask = np.where(img > splitter)\n                    img_to_show[mask] = 0.\n                    plt.imshow(img_to_show)\n                    plt.subplot(122)\n                    img_to_show = 1. * img\n                    mask = np.where(img < splitter)\n                    img_to_show[mask] = 0.\n                    plt.imshow(img_to_show)\n                    plt.suptitle('Background and foreground (yet unclassified)')\n                    plt.show()\n\n                # Perform statistics of the contiguous regions in the flattened image thresholded by splitter. If the sizes are only few pixels, the image is likely noise-dominated -- convolve it with a gaussian kernel and redo the analysis. Note: I added this feature to the code after viewing the cases that failed during 1st pass of stage2_test_final.\n                flattened_cells_image = np.array(1. * img)\n                mask = np.where(img < splitter)\n                flattened_cells_image[mask] = 0\n                mask = np.where(img > splitter)\n                flattened_cells_image[mask] = 1\n                flattened_cells_image = flattened_cells_image.astype(np.int)\n                flattened_cells_image = np.array(flattened_cells_image).flatten()\n                if np.sum(flattened_cells_image) > 0.85 * len(flattened_cells_image) and not force_gmm:\n                    if verbose:\n                        print(\"Threshold too low, forcing GMM.\")\n                    force_gmm = True\n                    continue\n                flattened_cells_image = np.array([flattened_cells_image[:-1], flattened_cells_image[1:]]).transpose()\n                contiguous_starts = list(np.where(list(map(lambda x : x[0] == 0 and x[1] == 1, flattened_cells_image)))[0])\n                contiguous_stops = list(np.where(list(map(lambda x : x[0] == 1 and x[1] == 0, flattened_cells_image)))[0])\n                if flattened_cells_image[0][0] > 0: contiguous_starts.insert(0,-1)\n                if flattened_cells_image[-1][1] > 0: contiguous_stops.append(len(flattened_cells_image))\n                contiguous_sizes = np.array(contiguous_stops) - np.array(contiguous_starts)\n                if (np.median(contiguous_sizes)) < 3.:\n                    if len(contiguous_sizes) < 10 and force_gmm == False:\n                        if verbose:\n                            print(\"Simple background detection failed. Forcing Gaussian Mixture Model.\")\n                        force_gmm = True\n                    else:\n                        if verbose:\n                            print(\"Noisy image detected. Restarting background detection with a blurred image.\")\n                        # convolve with a flat kernel and redo\n                        blurring_kernel_size = int(1.5 * blurring_kernel_size)\n                        img = cv2.filter2D(img0, -1, np.ones((blurring_kernel_size,blurring_kernel_size), np.float32)/(blurring_kernel_size**2))\n                        #img = cv2.GaussianBlur(img, (5,5), 0)\n                        #n_bins = 32\n                        if not force_gmm:\n                            force_gmm_brightest = True # in noise dominated images, the nuclei seem to be the brightest objects (imaging method?)\n                            force_gmm = True # force the Gaussian mixture model background (two-element model is too simple for blurred images)\n                else:\n                    break\n\n            except Exception as e:\n                print(\"Exception caught: %s. \" % e, end='')\n                if n_bins < 129:\n                    print('Attempting to use more bins.')\n                    n_bins = n_bins * 2\n                else:\n                    print('Forcing GMM.')\n                    n_bins = 32\n                    force_gmm = True\n                continue\n\n        return splitter, blurring_kernel_size\n        \n    except Exception as e:\n        print(\"Could not process image: %s. Error: %s\" % (imgname, e))\n        return 0.5\ndef detect_all_backgrounds (DATA, dataset, folder='stage1_train', visualize=False):\n    print('Detecting backgrounds for %s' % dataset)\n    thresholds = {}\n    for imgname in DATA['images'][dataset]:\n        thresholds[imgname] = detect_background(dataset, imgname, folder=folder, visualize=visualize)\n    return thresholds","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a94b05e12b75899e27b2008a71e43359cadbbe33"},"cell_type":"markdown","source":"**EXAMPLES START HERE**"},{"metadata":{"trusted":true,"_uuid":"a6093c7e4b4ed4e918bf2ad35f6277a9ddac8c48"},"cell_type":"code","source":"folder = 'stage2_test_final' # directory holding the dataset to be analysed\n\n# import data\nDATA = import_data(folder=folder)\ndatasets = list(DATA['images'].keys())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"5002385295e000967fab819bcfab8bee9e443c2e"},"cell_type":"code","source":"# case 1: pixel brightness histogram fit as a Lorentzian and a Gaussian\ndataset = '0154f82b5f214a91e3475d0b14d9bbe93960b7b1f925da9bb8e2b4aa65ec006a'\nthresholds = detect_all_backgrounds(DATA, dataset, folder=folder, visualize=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"db05c2a500ecf9d405e728c9f6c53f4e3e5f4d78"},"cell_type":"code","source":"# case 2: using blurring to fix noisy images\ndataset = '2ac1e92c749e54b5bad111be1220d0ea07dd176a48c51fac4c212f2b7aeec8af'\nsplitters = detect_all_backgrounds(DATA, dataset, folder=folder, visualize=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"19b5c2d72bd40fb309b7c69e87657c3385516c96"},"cell_type":"code","source":"# case 4: using a Gaussian Mixture Model decomposition\ndataset = '63bc2a0c564882b414349b3485208e50b47aa56f8f0e9cc1fe1e0fe8cf60d1f6'\nthresholds = detect_all_backgrounds(DATA, dataset, folder=folder, visualize=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dfd463e4e5aee31ead51964f3b45179116e9fd83"},"cell_type":"code","source":"# case 5: it's still far from perfect\ndataset = 'ec54bb85006c1a4ebfcc87aa0d369c6ce1bf85ef898c613b490c9dcce379e096'\nthresholds = detect_all_backgrounds(DATA, dataset, folder=folder, visualize=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"e9d026ff96503575240530e287799c0b8bf09690"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}