{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Image Preprocessing\n\nThis notebook demonstrates how to:\n- read image data from a parquet file via a batch generator function\n- normalize, reshape, resize, denoise, center and crop images using Numpy  \n- use tensorflow ImageDataGenerator to create generated data\n\nWe will be using the (Bengali.AI Handwritten Grapheme)[https://www.kaggle.com/c/bengaliai-cv19] dataset. \n\nFor an example of how this code can be used as part of a data pipeline for a CNN Neural Network, see the my writeup of my competition entry:\n- https://www.kaggle.com/jamesmcguigan/bengali-ai-cnn-data-pipeline-problem-solving"},{"metadata":{"trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport tensorflow as tf\nimport matplotlib.pyplot as plt\nimport glob2\nimport gc\nimport pyarrow\nimport time\nimport math\nimport skimage\nimport skimage.measure\nimport skimage.filters\nimport scipy\nfrom timeit import timeit\nfrom time import sleep\nfrom pyarrow.parquet import ParquetFile\nfrom typing import Callable\nfrom pandas import DataFrame\n\n\nfrom tensorflow.keras.preprocessing.image import ImageDataGenerator\nfrom collections import Counter\nfrom itertools import chain\nfrom functools import reduce\n\npd.set_option('display.max_columns',   500)\npd.set_option('display.max_colwidth',   -1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plt_imshow_batch( batch_X, cols=8, cmap='gray' ):\n    imgs = (batch_X.squeeze() * 255).astype('uint8')\n    rows = math.ceil(len(imgs)/cols)\n        \n    fig, ax = plt.subplots(figsize=(20, rows*2))    \n    plt.axis('off')        \n    for index, img in enumerate(imgs):\n        plt.subplot(rows, cols, index+1)\n        plt.imshow(img, cmap=cmap)    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Source: https://www.kaggle.com/jamesmcguigan/reading-parquet-files-ram-cpu-optimization/\ndef parquet_batch_generator(batch_size=128, reads_per_file=2):\n    filenames = sorted(glob2.glob('../input/bengaliai-cv19/train_image_data_*.parquet')); filenames    \n    for filename in filenames:\n        num_rows    = ParquetFile(filename).metadata.num_rows\n        cache_size  = math.ceil( num_rows / batch_size / reads_per_file ) * batch_size\n        batch_count = math.ceil( cache_size / batch_size )\n        for n_read in range(reads_per_file):\n            cache = pd.read_parquet(filename).iloc[ cache_size * n_read : cache_size * (n_read+1) ].copy()\n            gc.collect(); sleep(1);  # sleep(1) is required to allow measurement of the garbage collector\n            for n_batch in range(batch_count):            \n                yield cache[ batch_size * n_batch : batch_size * (n_batch+1) ].copy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Source: https://stackoverflow.com/a/31402351/748503\ndef bounding_box(img):\n    r = np.any(img, axis=(1, 2))\n    c = np.any(img, axis=(0, 2))\n    z = np.any(img, axis=(0, 1))\n\n    rmin, rmax = np.where(r)[0][[0, -1]]\n    cmin, cmax = np.where(c)[0][[0, -1]]\n    zmin, zmax = np.where(z)[0][[0, -1]]\n    \n    return rmin, rmax, cmin, cmax, zmin, zmax    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Source: https://codereview.stackexchange.com/questions/132914/crop-black-border-of-image-using-numpy\n\n# This is the fast method that simply remove all empty rows/columns\n# NOTE: assumes inverted\ndef crop_image(img, tol=0):\n    mask = img > tol\n    return img[np.ix_(mask.any(1),mask.any(0))]\n\n\n# DOCS: https://docs.scipy.org/doc/numpy/reference/generated/numpy.pad.html\n# NOTE: assumes inverted\ndef crop_center_image(img, tol=0):\n    org_shape   = img.shape\n    img_cropped = crop_image(img)\n    pad_x       = (org_shape[0] - img_cropped.shape[0])/2\n    pad_y       = (org_shape[1] - img_cropped.shape[1])/2\n    padding     = (\n        (math.floor(pad_x), math.ceil(pad_x)),\n        (math.floor(pad_y), math.ceil(pad_y))\n    )\n    img_center = np.pad(img_cropped, padding, 'constant', constant_values=0)\n    return img_center\n\n\n# Source: https://www.kaggle.com/jamesmcguigan/bengali-ai-image-processing/\n# noinspection PyArgumentList\ndef transform_X(train: DataFrame, denoise=True, normalize=True, center=True, invert=True, resize=2, resize_fn=None) -> np.ndarray:\n    train = (train.drop(columns='image_id', errors='ignore')\n             .values.astype('uint8')                   # unit8 for initial data processing\n             .reshape(-1, 137, 236)                    # 2D arrays for inline image processing\n    )\n    gc.collect(); sleep(1)\n\n    # Invert for processing\n    # Colors   |   0 = black      | 255 = white\n    # invert   |   0 = background | 255 = line\n    # original | 255 = background |   0 = line\n    train = (255-train)\n    \n    if denoise:                     \n        # Rescale lines to maximum brightness, and set background values (less than 2x mean()) to 0\n        train = np.array([ train[i] + (255-train[i].max())              for i in range(train.shape[0]) ])        \n        train = np.array([ train[i] * (train[i] >= np.mean(train[i])*2) for i in range(train.shape[0]) ])                                  \n            \n    if isinstance(resize, bool) and resize == True:\n        resize = 2    # Reduce image size by 2x\n    if resize and resize != 1:                  \n        # NOTEBOOK: https://www.kaggle.com/jamesmcguigan/bengali-ai-image-processing/\n        # Out of the different resize functions:\n        # - np.mean(dtype=uint8) produces produces fragmented images (needs float16 to work properly - but RAM intensive)\n        # - np.median() produces the most accurate downsampling\n        # - np.max() produces an enhanced image with thicker lines (maybe slightly easier to read)\n        # - np.min() produces a  dehanced image with thiner lines (harder to read)\n        resize_fn = resize_fn or (np.max if invert else np.min)\n        cval      = 0 if invert else 255\n        train = skimage.measure.block_reduce(train, (1, resize,resize), cval=cval, func=resize_fn)\n            \n    if center:\n        # NOTE: crop_center_image() assumes inverted\n        train = np.array([\n            crop_center_image(train[i,:,:])\n            for i in range(train.shape[0])\n        ])\n        \n    # Un-invert if invert==False\n    if not invert: train = (255-train)\n\n    if normalize:\n        train = train.astype('float16') / 255.0   # prevent division cast: int -> float64\n\n    train = train.reshape(*train.shape, 1)        # 4D ndarray for tensorflow CNN\n\n    gc.collect(); sleep(1)\n    return train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Without preprocessing\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=False, center=False, invert=False, resize=False)\n    plt_imshow_batch(batch_X)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=False, center=False, resize=False)\n    plt_imshow_batch(batch_X)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=False,  resize=False)\n    plt_imshow_batch(batch_X)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert + Center\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=True,  resize=False)\n    plt_imshow_batch(batch_X, cols=4)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"Out of the different resize functions:\n- np.mean(dtype=uint8) produces distored results (needs float16 to work properly - but RAM intensive)\n- np.median() produces the most accurate downsampling\n- np.max() produces an enhanced image with thicker lines (maybe slightly easier to read)\n- np.min() produces a  dehanced image with thiner lines (harder to read)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Center + Resize Default (without Invert)\nmean = lambda x, axis: np.mean(x, axis=axis, dtype=np.uint8)  # func_kwargs={} arg not in pip version    \nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=False, center=True, resize=2)\n    plt_imshow_batch(batch_X, cols=8)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert + Center + Resize Mean uint8\nmean = lambda x, axis: np.mean(x, axis=axis, dtype=np.uint8)  # func_kwargs={} arg not in pip version    \nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=True, resize=2, resize_fn=mean)\n    plt_imshow_batch(batch_X, cols=8)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert + Center + Resize Median\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=True, resize=2, resize_fn=np.median)\n    plt_imshow_batch(batch_X, cols=8)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert + Center + Resize Median\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=True, resize=2)\n    plt_imshow_batch(batch_X, cols=8)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert + Center + Resize Max\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=True, resize=2, resize_fn=np.max)\n    plt_imshow_batch(batch_X, cols=8)\n    break","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Denoise + Invert + Center + Resize Min\nfor batch_raw in parquet_batch_generator(16):\n    batch_X = transform_X(batch_raw, denoise=True, invert=True, center=True, resize=2, resize_fn=np.min)\n    plt_imshow_batch(batch_X, cols=8)\n    break","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# ImageDataGenerator"},{"metadata":{"trusted":true},"cell_type":"code","source":"for train_batch in parquet_batch_generator(8):\n    train_batch = transform_X( train_batch, normalize=False )  # ImageDataGenerator(rescale=1./255) saves float16 conversion\n    break\n( train_batch.dtype, train_batch.shape )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# DOCS: https://keras.io/preprocessing/image/#imagedatagenerator-class\ntime_start = time.time()\ntrain_datagen = ImageDataGenerator(\n    rescale=1./255,\n#     featurewise_center=True,             # No visible effect in plt.imgshow() \n#     samplewise_center=True,              # No visible effect in plt.imgshow()\n#     featurewise_std_normalization=True,  # No visible effect in plt.imgshow() | requires .fit()\n#     samplewise_std_normalization=True,   # No visible effect in plt.imgshow() | requires .fit()\n#     zca_whitening=True,                   # Kaggle, insufficent memory\n\n    zoom_range=0.2,\n    rotation_range=45/2,\n    shear_range=45/2,\n    width_shift_range=0.1,    # we already have centering\n    height_shift_range=0.1,   # we already have centering\n    brightness_range=(0.5,2),\n    fill_mode='constant',\n    cval=0,        \n)\ntrain_datagen.fit( train_batch )\nprint('time: datagen.fit()', time.time() - time_start)\n\ntrain_generator = train_datagen.flow(\n    train_batch,\n    subset=\"training\",\n    shuffle=False,  \n    batch_size=32,  # cant be larger than train_batch\n)\ngc.collect()\n\ncount = 0\nfor datagen_batch in train_generator:\n    plt_imshow_batch(datagen_batch, cols=8)\n    count += datagen_batch.shape[0]\n    if count >= 8*8: break\n        \nprint('time: datagen.fit() + flow()', time.time() - time_start)        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note that ZCA is very slow, and only ever seems to produce a single image"},{"metadata":{"trusted":true},"cell_type":"code","source":"time_start = time.time()\nzca_datagen = ImageDataGenerator(\n    rescale=1./255,\n    featurewise_center=True,             # No visible effect in plt.imgshow() \n    samplewise_center=True,              # No visible effect in plt.imgshow()    \n    zca_whitening=True, # Kaggle, insufficent memory\n)\nzca_datagen.fit( train_batch )\nprint('time: zca.fit()', time.time() - time_start)\n\nzca_generator = zca_datagen.flow(\n    train_batch,\n    subset=\"training\",\n    shuffle=False,  \n    batch_size=8,  # cant be larger than train_batch\n)\ncount = 0\nfor datagen_batch in zca_generator:\n    plt_imshow_batch(datagen_batch, cols=8)\n    count += datagen_batch.shape[0]\n    if count >= 8*8: break\n        \nprint('time: zca.fix() + flow()', time.time() - time_start)        ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}