{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"0370e2f9-aabc-ca3f-3a4c-9c01194a1f4c"},"source":"This kernel introduced preprocessing procedures for cervical images, including ROI (region of interest) croping and Specular Reflections removing. The main ideas came from two papers, please read them to fine more detailed infomations): \n\n- Ref. 1 [Automatic Detection of Anatomical Landmarks in Uterine Cervix Images](https://www.researchgate.net/profile/Sameer_Antani/publication/24041301_Automatic_Detection_of_Anatomical_Landmarks_in_Uterine_Cervix_Images/links/0deec51cdce9261312000000/Automatic-Detection-of-Anatomical-Landmarks-in-Uterine-Cervix-Images.pdf)\n\n- Ref.2 [Automatic Detection of Specular Reflections in Uterine Cervix Images](http://paloma.eng.tau.ac.il/research/laboratories/mip_lab/publications/cervix_sr_spie06.pdf)"},{"cell_type":"markdown","metadata":{"_cell_guid":"90631bfb-b756-5e13-5ae4-139d00e02caa"},"source":"## ROI Cropping ##"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"2fda0b4e-383c-8458-181a-c9f0f06fe537"},"outputs":[],"source":"# packeges importing\nimport pandas as pd\nimport numpy as np\nimport skimage\nfrom skimage import io, transform, morphology, segmentation, measure\nfrom sklearn.mixture import GaussianMixture\nfrom sklearn.preprocessing import MinMaxScaler\nimport glob\nimport os\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nDIR_BASE = '../input'\nimgs = glob.glob(os.path.join(DIR_BASE, 'train/Type_1', '*.jpg'))\nnum_imgs = len(imgs)\nprint('Number of images of type_1 in train size: {}'.format(num_imgs))"},{"cell_type":"markdown","metadata":{"_cell_guid":"f6b79a48-9392-c7dc-880c-18c691d6d7cc"},"source":"Let's first define a function which reads in an image as numpy array and resize it to smaller size, then show a sample image."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"2b5d468f-d4dd-69e1-4021-41f8bad47082"},"outputs":[],"source":"SHORTER_EDGE = 400  # image's shorter edge after resizing\ndef read_img(img_name):\n    img = skimage.io.imread(img_name)\n    h0, w0, _ = img.shape\n    if h0 >= w0:\n        w = SHORTER_EDGE\n        h = int(h0/w0*w)\n    else:\n        h = SHORTER_EDGE\n        w = int(w0/h0*h)\n    img = transform.resize(img, (h,w))\n    img_lab = skimage.color.rgb2lab(img) # Later the LAB color space will be used.\n\n    return img, img_lab\n\nimg_name = imgs[60]\nimg, img_lab = read_img(img_name)\nplt.imshow(img)"},{"cell_type":"markdown","metadata":{"_cell_guid":"0e14c642-cbeb-4aa5-e242-27a87442a1a9"},"source":"As paper *Ref. 1* suggested: \n\n> The cervix region is a relatively pink region located near the image center.  We use two features: the *A* color channel of the *LAB* color space (the higher the value of *A*, the “redder” the pixel color) and the distance *R* of a pixel from the image center. The *R* feature provides spatial information and supports the extraction of continuous regions within the image plane. The image is separated next into two clusters in the 2-D (*A*-*R*) feature space; we use Gaussian mixture modeling, initialized by a K-means procedure, as a statistical clustering methodology. The cluster that has the highest *A*-mean and the lowest *R*-mean is selected. The ROI is chosen as the largest connected component within the pixels associated with this cluster. \n\nLet's implement this idea in PYTHON. The following function *find_roi_by_gsmix* clusters the images' pixels (all or in the mask region if mask is provided) into two classes, and return a mask represents the ROI."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"75a5a2a0-8660-5919-85cc-5bc4d48dd47d"},"outputs":[],"source":"# find ROI by using Gaussian mixture modeling\ndef find_roi_by_gsmix(img, mask=None):\n    h,w,_ = img.shape\n    x_coor = np.repeat(range(h), w)  # for calculating the R\n    y_coor = np.tile(range(w), h)\n    if mask is None:\n        center = [h/2, w/2]\n    else:\n        mask = mask.reshape(-1)\n        center = [np.mean(x_coor[mask==1]), np.mean(y_coor[mask==1])]\n    R = np.sqrt((x_coor-center[0])**2 + (y_coor-center[1])**2)  # R\n    A = img_lab[:,:,1].reshape(-1)  # A\n    Ra = np.vstack([R, A]).T  # concat R and A\n\n    scaler = MinMaxScaler()\n    Ra = scaler.fit_transform(Ra)\n    gs_mix = GaussianMixture(n_components=2, random_state=42, init_params='kmeans') # Gaussian mixture modele\n    gs_mix.fit(Ra)\n    labels = gs_mix.predict(Ra)\n    \n    # Cluster with lowest R-mean will be chosen as ROI\n    means = gs_mix.means_\n    if means[0,0] < means[1,0]:\n        labels = 1 - labels\n    mask = labels.reshape(h, w)\n\n    return mask\n\nroi_mask = find_roi_by_gsmix(img)\nplt.imshow(roi_mask, 'gray')"},{"cell_type":"markdown","metadata":{"_cell_guid":"dee9823a-c63e-0076-90b4-559d95c91bef"},"source":"See, it seems fairly reasonable. But, it's not enough. We then need perform some postprocessing. The first thing we should do is removal of the small separated regions in the mask shown above, hence only the largest continuous region is keeping alive.  The following function *postprocess_mask* let it come true. Well, if you want to perform some morphology operations, just set the *morp* parameter to *True* when you call this function."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"03c14790-90e6-da4d-2ae3-2bd6fd6a76a9"},"outputs":[],"source":"def postprocess_mask(mask, morp=False):    \n    \n    selem = skimage.morphology.disk(5)\n    mask = skimage.morphology.binary_erosion(mask,selem)\n    \n    labels_mask = measure.label(mask)\n    regions = measure.regionprops(labels_mask)\n    regions.sort(key=lambda x: x.area, reverse=True)\n    if len(regions) > 1:\n        for rg in regions[1:]:\n            labels_mask[rg.coords[:,0], rg.coords[:,1]] = 0\n    labels_mask[labels_mask!=0] = 1\n    mask = labels_mask\n    \n    selem = skimage.morphology.disk(5)\n    mask = skimage.morphology.binary_dilation(mask,selem)\n\n    if morp:\n        selem = skimage.morphology.disk(10)\n        mask = skimage.morphology.binary_erosion(mask,selem)\n        mask = skimage.morphology.binary_dilation(mask,selem)\n\n    return mask\n\nmask = postprocess_mask(roi_mask, True)\nplt.imshow(mask, 'gray')"},{"cell_type":"markdown","metadata":{"_cell_guid":"c5530f6b-4116-36e2-2662-c3f63dbe6975"},"source":"Now we can crop the ROI off the original image. Firstly, we need a function to define the rectangle which will be used to crop the image. The smallest circumscribed rectangle is our interest."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"53bac652-f07a-eb6e-eb0e-da7726a6244d"},"outputs":[],"source":"def find_rectangle_border(mask):\n    h, w = mask.shape\n    left = np.argmax(mask, axis=1)\n    left = left[left!=0].min()\n    right = np.argmax(mask[:,::-1], axis=1)\n    right = w - right[right!=0].min()\n    up = np.argmax(mask, axis=0)\n    up = up[up!=0].min()\n    down = np.argmax(mask[::-1,:], axis=0)\n    down = h - down[down!=0].min()\n\n    return up, down, left, right\n\ndef segment_img(img, rectangle):\n    up, down, left, right = rectangle\n    img_out = img[up:down, left:right]\n\n    return img_out\n\nrect = find_rectangle_border(mask)\nimg_out = segment_img(img, rect)\nplt.imshow(img_out)"},{"cell_type":"markdown","metadata":{"_cell_guid":"d660bdb9-ef59-16cc-080c-9540cd17b6c4"},"source":"Yes, we did it! You may notice that some bling bling points in the image. These are the specular reflections which interfere with numerous computer vision tasks. It’s better to identificate them and remove them. "},{"cell_type":"markdown","metadata":{"_cell_guid":"2f8fdbd8-f0c4-6abc-337a-1d73713a5e14"},"source":"## Specular Reflections Removing ##"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c26cef2e-72fc-1f77-961d-b2fcdb71ee79"},"outputs":[],"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0ea46da4-c43e-526a-4e3d-dcc78e1ec9ea"},"outputs":[],"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"232b68e6-1b13-84d9-ac59-28a679f28da2"},"outputs":[],"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"ad5cf13a-43f8-b3fc-89b8-fb82b795c284"},"outputs":[],"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"fa62ea68-8b4a-5ce5-a1ae-8d886e5743db"},"outputs":[],"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"e54ce07e-1a14-8843-f02c-7ef5184ab234"},"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}