{"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":"# Preamble\n\nTraining on manually labelled ink has one downside:\n\nWhen label is not accurate model is training to recognize patterns that are not ink as containing ink.\nThis can significantly increase amount of false positives or affect overall ability of model to train.\n\nLuckily we are given a more reliable source of truth - Infra Red images, that were used for manual labelling.\nLet's see if we can use some automated data processing on same infrared images to improve quality of both ink labels and training mask.","metadata":{"papermill":{"duration":0.004278,"end_time":"2023-05-17T05:47:32.879468","exception":false,"start_time":"2023-05-17T05:47:32.875190","status":"completed"},"tags":[]}},{"cell_type":"code","source":"import numpy as np\nimport glob\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom io import StringIO\nfrom sklearn.metrics import fbeta_score\nfrom skimage.util import view_as_windows\nfrom scipy.ndimage import distance_transform_edt, gaussian_filter\nfrom numba import jit","metadata":{"papermill":{"duration":3.062464,"end_time":"2023-05-17T05:47:35.945878","exception":false,"start_time":"2023-05-17T05:47:32.883414","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-05-27T07:44:43.177927Z","iopub.execute_input":"2023-05-27T07:44:43.178337Z","iopub.status.idle":"2023-05-27T07:44:43.186151Z","shell.execute_reply.started":"2023-05-27T07:44:43.178304Z","shell.execute_reply":"2023-05-27T07:44:43.185027Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Constants\nPREFIX = '/kaggle/input/vesuvius-challenge-ink-detection/train/3/'\n\n# Load mask image\nmask = np.array(Image.open(PREFIX+\"mask.png\").convert('1'))\n\n# Load label image\nlabel = (np.array(Image.open(PREFIX+\"inklabels.png\")) > 0).astype(np.float32)\n\n# Load infrared image\nir = np.array(Image.open(PREFIX+\"ir.png\"))","metadata":{"papermill":{"duration":1.169647,"end_time":"2023-05-17T05:47:37.119688","exception":false,"start_time":"2023-05-17T05:47:35.950041","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-05-27T07:44:43.188188Z","iopub.execute_input":"2023-05-27T07:44:43.189240Z","iopub.status.idle":"2023-05-27T07:44:44.410180Z","shell.execute_reply.started":"2023-05-27T07:44:43.189196Z","shell.execute_reply":"2023-05-27T07:44:44.409019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Image value\n\nFirst, let's see how values are distributed in the infrared image.","metadata":{}},{"cell_type":"code","source":"def plot_horizontal(images):\n    fig, ax = plt.subplots(1, len(images), figsize=(14,5))\n    for i, image in enumerate(images):\n        ax[i].imshow(image, cmap='gray')\n    plt.show()\n\nplot_horizontal([mask, ir, label])\nprint({ 'min': ir.min(), 'mean': ir.mean(), 'max': ir.max() })\nplot_horizontal([ir < 255, ir < 224, ir < 196, ir < 164])","metadata":{"execution":{"iopub.status.busy":"2023-05-27T07:44:44.411421Z","iopub.execute_input":"2023-05-27T07:44:44.411821Z","iopub.status.idle":"2023-05-27T07:44:54.223802Z","shell.execute_reply.started":"2023-05-27T07:44:44.411785Z","shell.execute_reply":"2023-05-27T07:44:54.222986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Interesting facts:\n- there are no values near zero. this can be useful if later we will want to use zero values as no-data\n- values above mean value can for the most part be found only outside of the papirus surface. Only exceptions are occasional bright spots on the papyrus with relatively small areas\n\nBy setting area around scroll to zero we will distinguish between data and non-data on IR image, that will help in further analysis.\n\nWe can use medium size gausian filter to distinguish between small bright spots on the surface from wide bright areas on edges and around papyrus.\nAs a bonus some part of border area of papyrus, where data is not reliable, will also be filtered out.","metadata":{}},{"cell_type":"code","source":"ir_mean = 180\nnew_mask = gaussian_filter(ir,sigma=40) < ir_mean\n\nplot_horizontal([ir < ir_mean, new_mask, mask])\n\nir[new_mask == 0] = 0\nlabel[new_mask == 0] = 0\nprint(ir[new_mask].min(), ir[new_mask].mean(), ir[new_mask].max())\nplot_horizontal([label, ir, new_mask])","metadata":{"execution":{"iopub.status.busy":"2023-05-27T07:44:54.226202Z","iopub.execute_input":"2023-05-27T07:44:54.226529Z","iopub.status.idle":"2023-05-27T07:45:22.314613Z","shell.execute_reply.started":"2023-05-27T07:44:54.226501Z","shell.execute_reply":"2023-05-27T07:45:22.313739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Nice! By this simple manipulation we improved data mask to only include surface of papyrus and updated our Infra red image so that values from outside of papyrus don't interfere with values inside papyrus.\n\nTime to try value to label correlation analysis!","metadata":{}},{"cell_type":"markdown","source":"# Image value to label correlation\n\nLet's see if we can rely solely on brightness value of infrared image to distinguish ink from non-ink.\n\nI'll take some code from my previous <> Kernel.","metadata":{}},{"cell_type":"code","source":"@jit(nopython=True)\ndef get_value_ink_ratio(value_count_ink, value_count_all, a, label):\n    for v, l in zip(a.ravel(), label.ravel()):\n        value_count_all[v] += 1\n        if l:\n            value_count_ink[v] += 1","metadata":{"papermill":{"duration":0.290446,"end_time":"2023-05-17T05:47:37.416342","exception":false,"start_time":"2023-05-17T05:47:37.125896","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-05-27T07:45:22.316087Z","iopub.execute_input":"2023-05-27T07:45:22.316693Z","iopub.status.idle":"2023-05-27T07:45:22.397732Z","shell.execute_reply.started":"2023-05-27T07:45:22.316648Z","shell.execute_reply":"2023-05-27T07:45:22.396518Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_ink_ratio(value_count_ink, value_count_all, img, label, plot_th = None):\n    value_ink_ratio = np.where(value_count_all == 0, 0, value_count_ink/value_count_all)\n    x = np.arange(len(value_ink_ratio))\n    \n    # plot ink ratio distribution\n    fig, ax = plt.subplots(1, 1, figsize=(14,1))\n    ax.plot(x, value_ink_ratio, linestyle='', marker='.')\n    ax.set_title('Ink ratio by value')\n    plt.show()\n    \n    # select \n    sorted_by_ink = np.argsort(value_ink_ratio)\n    sorted_ink_ratio = value_ink_ratio[sorted_by_ink]\n    truth = value_count_all.sum()\n    f05 = np.zeros(101)\n    best_f05 = 0\n    best_th = 0\n    for th in range(101):\n        high_ink_ratio = sorted_by_ink[sorted_ink_ratio > th / 100]\n        tp = value_count_ink[high_ink_ratio].sum()\n        fp = value_count_all[high_ink_ratio].sum() - tp\n        fn = truth - tp\n        f05[th] = 1.25 * tp / (1.25 * tp + fp + 0.25 * fn)\n        if best_f05 < f05[th] + th / 1000:\n            best_th = th / 100\n            best_f05 = f05[th] + th / 1000\n            \n    fig, ax = plt.subplots(1, 1, figsize=(14,1))\n    ax.plot(np.arange(101) / 100, f05, linestyle='', marker='.')\n    ax.set_title('F05 score by threshold')\n    plt.show()\n            \n    high_ink_ratio = sorted_by_ink[sorted_ink_ratio > (plot_th if plot_th is not None else best_th)]\n    high_ink = np.isin(img, high_ink_ratio)\n    ratio_ink = value_ink_ratio[img]\n        \n    print('Best scoring threshold:', best_th)\n    print('Best score:', best_f05 - best_th / 10)\n    print('% of values:', len(high_ink_ratio)/len(value_ink_ratio))\n\n    fig, ax = plt.subplots(1, 4, figsize=(14,5))\n    ax[0].imshow(label, cmap='binary')\n    ax[1].imshow(img, cmap='binary')\n    ax[2].imshow(ratio_ink, cmap='binary')\n    ax[3].imshow(high_ink, cmap='binary')\n    ax[0].set_title('Ink labels')\n    ax[1].set_title('Source image')\n    ax[2].set_title('Ink ratio by pixel')\n    ax[3].set_title('Best prediction')\n    plt.show()\n    \n    return high_ink, ratio_ink, best_f05\n\nNUM_VALUES = 256    \nvalue_count_ink = np.zeros(NUM_VALUES, dtype=int)\nvalue_count_all = np.zeros(NUM_VALUES, dtype=int)\nget_value_ink_ratio(value_count_ink, value_count_all, ir, label)\n_ = plot_ink_ratio(value_count_ink, value_count_all, ir, label)","metadata":{"papermill":{"duration":9.002093,"end_time":"2023-05-17T05:47:46.422368","exception":false,"start_time":"2023-05-17T05:47:37.420275","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2023-05-27T07:45:22.399453Z","iopub.execute_input":"2023-05-27T07:45:22.400101Z","iopub.status.idle":"2023-05-27T07:45:31.467804Z","shell.execute_reply.started":"2023-05-27T07:45:22.400071Z","shell.execute_reply":"2023-05-27T07:45:31.466861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There is clear correlation between pixel values of Infrared image and manually created labels.\n\nBut this correlation not very strong.\nCurrently, due to different amount of light falling on different parts of the image, same values might correspond to ink on darker parts of the papirus and non-ink on its brighter parts.\n\nTo fix it we need to make brightness of papirus uniform across entire image.\n\nUltimately it should allow us to improve labels to correlate better with an infrared image.","metadata":{}},{"cell_type":"markdown","source":"# Uniform brightness\n\nOur infrared image already has areas not representing surface of papyrus removed.\n\nNext step is to also rule out values that are most likely to be ink, because ink has noticably different value on infrared image, which can affect brightness calculation.\n\nEasiest way to do so is by setting to zero pixels manually labelled as ink.\n\nI'm using small trick to apply wide gausian filter to the image while excluding zero valued data from calculations.\nIt is important to select gausian kernel to be wider than inked areas, so that letters are completely covered on resulting light map","metadata":{}},{"cell_type":"code","source":"# Weight of pixels with data that particitated in computations\nW = gaussian_filter(1.*new_mask*(label==0),200)\n\n# papyrus without ink\nlight = np.where(label, 0, ir)\n# wide gausian filter to leave only average brightness\nlight = np.where(W == 0, 0, gaussian_filter(light,200) / W)\nmin_light = light[new_mask].min()\nprint(light[new_mask].min(), light[new_mask].mean(), light[new_mask].max())\n\nplot_horizontal([label, W, ir, light, new_mask])","metadata":{"execution":{"iopub.status.busy":"2023-05-27T07:45:31.469261Z","iopub.execute_input":"2023-05-27T07:45:31.469561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Pretty good. Next step is to remove spare light from the image in order to get uniformly lighted version.\n\nNote that brightness calcualtion by gauss interpolation is not perfect and only partly represents light map\nI use small amplifying multiplier to compensate for not 100% of light being \"detected\".","metadata":{}},{"cell_type":"code","source":"min_light = light[new_mask].min()\n\nnew_ir = (ir - (light - min_light) * 1.05)\nnew_ir[new_mask == 0] = 0\nprint(new_ir[new_mask].min(), new_ir[new_mask].mean(), new_ir[new_mask].max())\nnew_ir = np.around(new_ir).astype(np.uint8)\n\nplot_horizontal([label, new_ir, ir, new_mask])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Nice, image looks even better then original, is not it?\n\nLet see if correlation with a label has increased!","metadata":{}},{"cell_type":"code","source":"value_count_ink = np.zeros(NUM_VALUES, dtype=int)\nvalue_count_all = np.zeros(NUM_VALUES, dtype=int)\nget_value_ink_ratio(value_count_ink, value_count_all, new_ir, label)\nhigh_ink, ratio_ink, _ = plot_ink_ratio(value_count_ink, value_count_all, new_ir, label)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Much better. Other source of error can be from noisy image, let's try to blur image a little bit and see if that helps.","metadata":{}},{"cell_type":"code","source":"blurred_ir = gaussian_filter(new_ir, 6)\nblurred_ir[new_mask == 0] = 0\nprint(blurred_ir[new_mask].min(), blurred_ir[new_mask].mean(), blurred_ir[new_mask].max())\n\nvalue_count_ink = np.zeros(NUM_VALUES, dtype=int)\nvalue_count_all = np.zeros(NUM_VALUES, dtype=int)\nget_value_ink_ratio(value_count_ink, value_count_all, blurred_ir, label)\nhigh_ink, ratio_ink, _ = plot_ink_ratio(value_count_ink, value_count_all, blurred_ir, label)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's zoom in a little bit...","metadata":{}},{"cell_type":"code","source":"slice = np.s_[2500:3500,1500:3500]\n_ = plot_ink_ratio(value_count_ink, value_count_all, blurred_ir[slice], label[slice])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Results\n\nArguably our prediction already better represents ink than manual labeling. While ink was likely continuous originally some of it might be already lost and prediction based on IR image can better represent where it happened.","metadata":{}},{"cell_type":"code","source":"\nblurred_ink = gaussian_filter(1.*ratio_ink, 5)\n\n_ = plot_horizontal([label, ir, blurred_ink, ratio_ink])\nprint(blurred_ink[blurred_ink > 0.01].min(), blurred_ink[blurred_ink > 0.01].mean(), blurred_ink[blurred_ink > 0.01].max())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"new_label = label * (blurred_ink > 0.2)\n_ = plot_horizontal([ir, blurred_ink > 0.2, new_label, label])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_horizontal([ir[slice], (blurred_ink > 0.2)[slice], new_label[slice], label[slice]])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Unfortunately we still have to rely on manual labels to filter out noise, however using only part of manual labels that are also supported by data extracted directly from IR image can slightly improve quality of predictions.\n\n# Save\n\nNow it's time to save updated mask and model to be used for training","metadata":{}},{"cell_type":"code","source":"Image.fromarray(new_mask).save('/kaggle/working/mask.png')\nImage.fromarray(new_label.astype(bool)).save('/kaggle/working/inklabel.png')","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}