{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"\n\n\n\n\n","metadata":{}},{"cell_type":"markdown","source":"# **How you might be able to clean the data**","metadata":{}},{"cell_type":"markdown","source":"* The first part until the data cleaning part was taken from @https://www.kaggle.com/its7171/mmdetection-for-segmentation-training\n* The basic idea is to use cell segmentor to get a data set where each image shows a single cell\n* Then training a classifier on single cell images\n* Such a dataset will contain faulty data because:\n    1. the labels correspond to a bunch of cells without\n    2. there is no guarantee that each cell actually expresses the protein of interest\n    3. even if the cell expresses it there is no guarantee that it is visible properly through the green marker for example due to photobleaching\n    4. overexposure might lead to bright signal in the green channel but the characteristic feature for this protein might not be visible due to loss of information \n    5. there are erros in the segmentation masks\n \n \n* In order to get rid of some of the faulty data the following criterion is applied:\n    1. from the pixels corresponding to no cell, the background level of the green channel is calculated\n    2. for each cell the area of pixels above this background level is divided by the cell's entire area (this quantity is refered to as fraction in the code)    \n    3. each cell with a fraction below a threshold and above another threshold is discarded as faulty\n    4. the thresholds are calulated by calculating two percentile of the distribution of the fractions of cells with the same protein of interest\n    \n   \n    ","metadata":{}},{"cell_type":"code","source":"from itertools import groupby\nimport numpy as np\nfrom tqdm import tqdm\nimport pandas as pd\nimport os\nimport pickle\nimport cv2\nfrom multiprocessing import Pool\nimport matplotlib.pyplot as plt","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exp_name = \"v3\"\nconf_name = \"mask_rcnn_s101_fpn_syncbn-backbone+head_mstrain_1x_coco\"\ncell_mask_dir = '../input/hpa-mask/hpa_cell_mask'    \nROOT = '../input/hpa-single-cell-image-classification/'\ntrain_or_test = 'train'\nimg_dir = f'../work/mmdet_{exp_name}_{train_or_test}'\n!mkdir -p {img_dir}\ndf = pd.read_csv(os.path.join(ROOT, 'train.csv'))\n\n# this script takes more than 9hours for full data.\ndebug = False\nif debug:\n    df = df[:4]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# helper funcs","metadata":{}},{"cell_type":"code","source":"# convert segmentation mask image to run length encoding\nMAX_GREEN = 64 # filter out dark green cells\ndef get_rles_from_mask(image_id, class_id):\n    mask = np.load(f'{cell_mask_dir}/{image_id}.npz')['arr_0']\n    if class_id != '18':\n        green_img = read_img(image_id, 'green')\n    rle_list = []\n    mask_ids = np.unique(mask)\n    for val in mask_ids:\n        if val == 0:\n            continue\n        binary_mask = np.where(mask == val, 1, 0).astype(bool)\n        if class_id != '18':\n            masked_img = green_img * binary_mask\n            #print(val, green_img.max(),masked_img.max())\n            if masked_img.max() < MAX_GREEN:\n                continue\n        rle = coco_rle_encode(binary_mask)\n        rle_list.append(rle)\n    return rle_list, mask.shape[0], mask.shape[1]\n\ndef coco_rle_encode(mask):\n    rle = {'counts': [], 'size': list(mask.shape)}\n    counts = rle.get('counts')\n    for i, (value, elements) in enumerate(groupby(mask.ravel(order='F'))):\n        if i == 0 and value == 1:\n            counts.append(0)\n        counts.append(len(list(elements)))\n    return rle\n\n# mmdet custom dataset generator\ndef mk_mmdet_custom_data(image_id, class_id):\n    rles, height, width = get_rles_from_mask(image_id, class_id)\n    if len(rles) == 0:\n        return {\n            'filename': image_id+'.jpg',\n            'width': width,\n            'height': height,\n            'ann': {}\n        }\n    rles = mutils.frPyObjects(rles, height, width)\n    bboxes = mutils.toBbox(rles)\n    bboxes[:, 2] += bboxes[:, 0]\n    bboxes[:, 3] += bboxes[:, 1]\n    return {\n        'filename': image_id+'.jpg',\n        'width': width,\n        'height': height,\n        'ann':\n            {\n                'bboxes': np.array(bboxes, dtype=np.float32),\n                'labels': np.zeros(len(bboxes)), # dummy data.(will be replaced later)\n                'masks': rles\n            }\n    }\n\n# print utility from public notebook\ndef print_masked_img(image_id, mask):\n    img = load_RGBY_image(image_id, train_or_test)\n    \n    plt.figure(figsize=(15, 15))\n    plt.subplot(1, 3, 1)\n    plt.imshow(img)\n    plt.title('Image')\n    plt.axis('off')\n    \n    plt.subplot(1, 3, 2)\n    plt.imshow(mask)\n    plt.title('Mask')\n    plt.axis('off')\n    \n    plt.subplot(1, 3, 3)\n    plt.imshow(img)\n    plt.imshow(mask, alpha=0.6)\n    plt.title('Image + Mask')\n    plt.axis('off')\n    plt.show()\n    \n# image loader, using rgb only here\ndef load_RGBY_image(image_id, train_or_test='train', image_size=None):\n    red = read_img(image_id, \"red\", train_or_test, image_size)\n    green = read_img(image_id, \"green\", train_or_test, image_size)\n    blue = read_img(image_id, \"blue\", train_or_test, image_size)\n    #yellow = read_img(image_id, \"yellow\", train_or_test, image_size)\n    stacked_images = np.transpose(np.array([red, green, blue]), (1,2,0))\n    return stacked_images\n\n# \ndef read_img(image_id, color, train_or_test='train', image_size=None):\n    filename = f'{ROOT}/{train_or_test}/{image_id}_{color}.png'\n    assert os.path.exists(filename), f'not found {filename}'\n    img = cv2.imread(filename, cv2.IMREAD_UNCHANGED)\n    if image_size is not None:\n        img = cv2.resize(img, (image_size, image_size))\n    if img.max() > 255:\n        img_max = img.max()\n        img = (img/255).astype('uint8')\n    return img\n\n# make annotation helper called multi processes\ndef mk_ann(idx):\n    image_id = df.iloc[idx].ID\n    class_id = df.iloc[idx].Label\n    anno = mk_mmdet_custom_data(image_id, class_id)\n    img = load_RGBY_image(image_id, train_or_test)\n    cv2.imwrite(f'{img_dir}/{image_id}.jpg', img)\n    return anno, idx, image_id","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# checking segment mask\nTo extract the each cells, [CellSegmentator](https://github.com/CellProfiling/HPA-Cell-Segmentation) can be used.\nAnd The extracted segment masks are stored in [this dataset](https://www.kaggle.com/its7171/hpa-mask).\n\nThis mask files are made as follows:\n<pre>\nnucl_mask, cell_mask = segmentCell(im, segmentator)\nnp.savez_compressed(f'{cell_dir}/{image_id}', cell_mask)\nnp.savez_compressed(f'{nucl_dir}/{image_id}', nucl_mask)\n</pre>\nSo you can load the mask as follows:\n<pre>\ncell_mask = np.load(f'{cell_dir}/{image_id}.npz')['arr_0']\nnucl_mask = np.load(f'{nucl_dir}/{image_id}.npz')['arr_0']\n</pre>\n","metadata":{}},{"cell_type":"code","source":"cell_mask_dir = '../input/hpa-mask/hpa_cell_mask'    \nfor idx in range(2):\n    image_id = df.iloc[idx].ID\n    cell_mask = np.load(f'{cell_mask_dir}/{image_id}.npz')['arr_0']\n    print_masked_img(image_id, cell_mask)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# **Data Cleaning**","metadata":{}},{"cell_type":"markdown","source":"# **Calculating fraction as metric for expression level of target protein**","metadata":{}},{"cell_type":"code","source":"def background(mask, ch):\n    bm = binary_mask(mask)\n    background = np.logical_not(bm) * ch\n    return background\n\ndef threshold(bg, p):\n    return np.quantile(bg, p)\n\n\ndef bbox(img):\n    rows = np.any(img, axis=1)\n    cols = np.any(img, axis=0)\n    rmin, rmax = np.where(rows)[0][[0, -1]]\n    cmin, cmax = np.where(cols)[0][[0, -1]]\n\n    return rmin, rmax, cmin, cmax\n\ndef crop(img, cell_mask, index):\n    bb = bbox(cell_mask == index)\n    return crop(img, bb)\n\ndef crop(img, bb):   \n    return img[bb[0] : bb[1], bb[2] : bb[3]]\n    \ndef binary_mask(mask):\n    bm = mask > 0\n    return bm\n\ndef fraction(croppedimg, croppedmask, th):\n    filtered = croppedimg * (croppedmask > 0)\n    s1 = np.sum(filtered > th)   \n    s2 = (np.sum(croppedmask > 0))\n    \n    frac = s1 / s2\n    return frac\n\ndef cell_mask_to_fraction(cell_mask, img, thp, chindex):\n    unique = np.unique(cell_mask)\n    ch = img[:,:,chindex]\n    bg = background(cell_mask, ch)\n    th = threshold(bg, thp)\n    \n    result =[]\n    for index in unique:\n        cri = crop(ch, cell_mask, index)\n        crm = crop(cell_mask, cell_mask, index)\n        \n        frac = fraction(cri, crm == index, th)\n        result.append(frac)\n    \n    return result     \n    \n    \ndef get_fraction(cell_mask_dir, image_id, chindex=1, thp=0.99):\n    idext = []\n    bbs = []\n    fracs = []\n    \n    cell_mask =load_cell_mask(image_id)\n    unique = np.unique(cell_mask)\n    img = load_RGBY_image(image_id)\n    \n    ch = img[:,:,chindex]\n    bg = background(cell_mask, ch)\n    th = threshold(bg, thp)  \n    \n    for i in unique[1:]:\n        cmtemp = cell_mask == i\n        bb = bbox(cmtemp)                \n        \n        idext.append(str(image_id) + '_' + str(i))\n        bbs.append(bb)\n        \n        cri = crop(ch, bb)\n        crm = crop(cell_mask, bb)\n        \n        frac = fraction(cri, crm == i, th)\n        fracs.append(frac)\n        \n    return idext, bbs, fracs     \n\ndef load_cell_mask(image_id):\n    cell_mask = np.load(f'{cell_mask_dir}/{image_id}.npz')['arr_0']\n    return cell_mask\n\n\ndef calculate_fractions(chindex=1, thp=0.99):\n    cell_mask_dir = '../input/hpa-mask/hpa_cell_mask' \n    idext, bbs, fracs = [], [], []\n    n = len(df)\n    for idx in range(n):\n        if idx % 100 == 0: print(str(idx) + '|' + str(n))    \n        image_id = df.iloc[idx].ID\n        if(df.iloc[idx].Label.find('|') != -1): continue    \n        idexttemp, bbstemp, fracstemp = get_fraction(cell_mask_dir, image_id, chindex, thp)\n        idext.extend(idexttemp)\n        bbs.extend(bbstemp)\n        fracs.extend(fracstemp)\n\n\n    df_fractions = {'cell_ID' : idext, 'bbox' : bbs, 'fraction' : fracs}    \n    df_fractions = pd.DataFrame(data=df_fractions)\n    \n    labels = []\n    i = 0\n    \n    print(\"Number of cells: \")\n    print(len(df_fractions))\n    print(\"Number of calculated fractions: \")\n    \n    for cid in df_fractions['cell_ID']:\n        if i % 100 == 0:   \n            print(i)\n        l = df.loc[df['ID'] == cid.split('_')[0]].Label\n        labels.append(l)\n        i += 1\n\n    df_fractions['Label'] = [k.to_list()[0] for k in labels]   \n    \n    df_fractions.to_csv('df_fractions.csv',index=False)\n    \n    return df_fractions","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"calculate_Fractions = False","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if calculate_Fractions:\n    df_fractions = calculate_fractions()\nelse:\n    df_fractions = pd.read_csv('../input/hpafractions/hpa_fractions.csv')\n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def rgby_from_frac_id(frac_id, bbox):\n    split = frac_id.split(\"_\")\n    image_id = split[0]\n    cell_mask_id = int(split[1])   \n    cell_mask = load_cell_mask(image_id)\n    cell_mask_crop = crop(cell_mask, bbox)\n    cell_mask_crop = cell_mask_crop == cell_mask_id\n    cell_mask_crop = cell_mask_crop.reshape(cell_mask_crop.shape[0], cell_mask_crop.shape[1], 1)\n    \n    img = crop(load_RGBY_image(image_id), bbox) * cell_mask_crop\n    \n    return img \n    \n        \ndef bbox_from_string(bboxstring):\n    bbox = bboxstring.split(\",\")\n    bbox[0] = bbox[0][1:]\n    bbox[-1] = bbox[-1][:-1]\n    \n    return [int(c) for c in bbox]\n\n\ndef get_entries_in_percentiles(label, pl, pu):\n    df_label = df_fractions[df_fractions['Label'] == label]\n    qu = np.quantile(df_label.fraction.to_list(), pu)    \n    ql = np.quantile(df_label.fraction.to_list(), pl)\n    \n    extremel = df_label[df_label['fraction'] < ql]    \n    extremeu = df_label[df_label['fraction'] > qu]\n    middle = df_label[df_label['fraction'] < qu]\n    middle = middle[middle['fraction'] > ql]    \n    \n    return extremel, middle, extremeu \n\ndef get_percentiles(label, pl, pu):\n    df_label = df_fractions[df_fractions['Label'] == label]\n    qu = np.quantile(df_label.fraction.to_list(), pu)    \n    ql = np.quantile(df_label.fraction.to_list(), pl)\n    \n    return ql, qu \n\ndef get_extremes(label, pl, pu):\n    df_label = df_fractions[df_fractions['Label'] == label]\n    qu = np.quantile(df_label.fraction.to_list(), pu)    \n    ql = np.quantile(df_label.fraction.to_list(), pl)\n    \n    extremel = df_label[df_label['fraction'] < ql]    \n    extremeu = df_label[df_label['fraction'] > qu]\n    \n    return extremel, extremeu \n\n\n\ndef visualize_random(df_sub, n, ci = None):\n    n_ev = n - (n % 2)\n    randindsl = ((len(df_sub) - 1) * np.random.uniform(0,1,n_ev)).astype(int)   \n    fig, axs = plt.subplots(int((len(randindsl) / 2)), 2, figsize=(15,15))\n    plotix = 0\n    plotiy = 0\n    for i in randindsl:\n        frac_id, bbox = df_sub.cell_ID.values[i], bbox_from_string(df_sub.bbox.values[i])\n        im = rgby_from_frac_id(frac_id, bbox)\n        \n        if ci is None:\n            axs[plotiy, plotix % 2].imshow(im)      \n        else: \n            axs[plotiy, plotix % 2].imshow(im[:,:,ci])\n\n        if plotix % 2 == 1: plotiy += 1\n        plotix += 1    \n        \n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Look at the distribution and visualize some examples","metadata":{}},{"cell_type":"code","source":"plt.hist(df_fractions[df_fractions['Label'] == 0].fraction.to_list())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"lowerpercent = 0.001\nupperpercent = 1 - 0.001\n\nextremel, middle, extremeu = get_entries_in_percentiles(0, lowerpercent, upperpercent)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Show some random example from below lower percentile:","metadata":{}},{"cell_type":"code","source":"visualize_random(extremel, 20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Show some random example between lower and upper percentile:","metadata":{}},{"cell_type":"code","source":"visualize_random(middle, 20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Show some random example above upper percentile:","metadata":{}},{"cell_type":"code","source":"visualize_random(extremeu, 20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ql, qu = get_percentiles(0, lowerpercent, upperpercent)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# **Output**","metadata":{}},{"cell_type":"code","source":"df_fractions[\"clean\"] = True\nlabels = np.unique(df_fractions.Label)\n\nfor l in labels:\n    ql, qu = get_percentiles(l, lowerpercent, upperpercent)\n    df_fractions[\"clean\"] = True\n    df_fractions.loc[np.bitwise_and(df_fractions['Label'] == l, df_fractions['fraction'] < ql), ['clean']] = False\n    df_fractions.loc[np.bitwise_and(df_fractions['Label'] == l, df_fractions['fraction'] > ql), ['clean']] = False","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_fractions.to_csv(\"df_fractions_clean.csv\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# **Outlook**","metadata":{}},{"cell_type":"markdown","source":"You could try different metrics for detecting faulty data like image entropy or brenner and do it for the other channels for example in order to detect if a nucleus is present","metadata":{}}]}