{"cells":[{"metadata":{"_cell_guid":"ebba0809-6591-4386-8e4a-ba8d1935dd57","_uuid":"af1935f19e6ffeca47d2933c39bfcaf5a816535a"},"cell_type":"markdown","source":"Original kernel is [Teaching notebook for total imaging newbies]((https://www.kaggle.com/stkbailey/teaching-notebook-for-total-imaging-newbies).\n\nThis kernel will implement classical image techniques and will hopefully serve as a useful primer to people who have never worked with image data before. Ultimately, we will develop a simple pipeline using `scipy` and `numpy` (and a little bit of `scikit-image`) that we can apply to the test images -- in fact, we won't even use the training images except to optimize parameters.\n\nI'll keep updating this notebook to try and improve it - user Gabro Vecsei takes a similar approach and scores 0.22 in [this kernel]((https://www.kaggle.com/gaborvecsei/basic-pure-computer-vision-segmentation-lb-0-229). My main intention here, though, is to help out people who are new to analyses, not to score highly. "},{"metadata":{"_cell_guid":"2b168289-7746-4255-b933-356bd04914b7","_uuid":"25fa830192ea11807891ec2f71487dbeb18364e3","trusted":true},"cell_type":"code","source":"import pathlib\nimport imageio\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\n# Glob the training data and load a single image path\ntraining_paths = pathlib.Path('../input/stage1_train').glob('*/images/*.png')\n\ntraining_sorted = sorted([x for x in training_paths])\nim_path = training_sorted[45]\nim = imageio.imread(str(im_path))\nprint(\"Total training images:%i\"%len(training_sorted))\nprint('im :'+str(format(im.shape)))","execution_count":1,"outputs":[]},{"metadata":{"_cell_guid":"aca64183-1c16-474c-88d6-6574a9dcf56e","_uuid":"1e41cf0c93629a559d0e62423b60b1f67691d9fc"},"cell_type":"markdown","source":"## Lectura de las máscaras\nFunciones varias"},{"metadata":{"_cell_guid":"bb23d8c6-8474-40aa-9113-10708b2c746c","_kg_hide-input":true,"_uuid":"de680bf8740e070c7daa01d4f74acd88bd75983e","trusted":true},"cell_type":"code","source":"from matplotlib.colors import ListedColormap\nnp.random.seed(1)\nlc=(np.random.rand(256,3)+0.4)/1.4\nlc[0]=[0,0,0]\nrand_cmap = ListedColormap(lc)\n\n# Función para obtener la región donde hay células\n# en una  máscara a partir del directorio de la imagen\ndef getMask(im_path):\n    yy=(im_path.parents[1]/'masks').glob('*.png')\n    yy=[x for x in yy]\n    im_mask=imageio.imread(str(yy[0]))\n    for dirIm in yy:\n        im_t=imageio.imread(str(dirIm))\n        im_mask=im_t | im_mask\n    return im_mask\n# Función para obtener las regiones donde hay células\n# en una sola imagen a partir del directorio de la imagen\ndef getRegion(im_path):\n    yy=(im_path.parents[1]/'masks').glob('*.png')\n    yy=[x for x in yy]\n    im_labels=imageio.imread(str(yy[0]))\n    im_labels= im_labels//255\n    for i,dirIm in enumerate(yy[1:]):\n        im_t=imageio.imread(str(dirIm))\n        im_labels=(im_t//255)*(i+1) + im_labels\n    return (im_labels, len(yy))\n\n#Función para graficar una célula dada la dirección de la imagen\n# la dirección es un path (pathlib)\ndef plotImg(im_path, figsize=(12,9),cmaps=('gnuplot', 'tab20', 'brg')):\n    im = imageio.imread(str(im_path))\n    im_mask= getMask(im_path)\n    im_regiones, nlabels = getRegion(im_path)\n    fondo= im[:,:,1] & (~ im_mask )\n    plt.figure(figsize=figsize)\n    plt.subplot(2,2,1)\n    plt.imshow(im[:,:,1],cmap=cmaps[0],norm=None,vmin=0,vmax=255)\n    plt.axis('off')\n    plt.title('Original Image')\n    plt.subplot(2,2,2)\n    plt.imshow(fondo, cmap=cmaps[0],norm=None, vmin=0,vmax=255)\n    plt.axis('off')\n    plt.title('Background')\n    plt.subplot(2,2,3)\n    plt.imshow(im_regiones, cmap=rand_cmap)\n    plt.axis('off')\n    plt.title('Regions (%s Nuclei)'%nlabels)\n    plt.subplot(2,2,4)\n    plt.imshow(im_mask, cmap=rand_cmap)\n    plt.axis('off')\n    plt.title('Mask')\nprint('Done')","execution_count":2,"outputs":[]},{"metadata":{"_cell_guid":"505d8f16-1f54-4132-8a60-d73a145fee89","_uuid":"092bf1fb8a05ca9b975180d305928ab250306765"},"cell_type":"markdown","source":"## Visualización de una imagen"},{"metadata":{"_cell_guid":"0a33baae-7708-40ff-add2-2d76063e8981","_uuid":"ac61bf401692b8d4a1933b4bcbe7ec08b8f6f07b","trusted":true},"cell_type":"code","source":"im_path = training_sorted[0]\nplotImg(im_path)","execution_count":3,"outputs":[]},{"metadata":{"_cell_guid":"4f6cff77-6108-427e-a555-ac12020b2c90","_uuid":"2c17d6e81a47698ad91cccef3800241539c8b185"},"cell_type":"markdown","source":"# Normalizando las imágenes\n\nEn el conjunto de datos se hallan diferentes tipos de imágenes, RGB,  y en escala de grises. Se ha estudiado que el canal que tiene mejores resultadose es el verde en imágenes a color, por lo tanto, se tomará el segundo canal."},{"metadata":{"_cell_guid":"1a227d66-1063-41d4-ae6a-d8a5ae8bac2a","_kg_hide-input":true,"_uuid":"2e1da694c4c9ee337215d2b88e723963ed87410b","trusted":true},"cell_type":"code","source":"from skimage.color import rgb2gray\nfrom scipy import signal\nfrom skimage.filters import threshold_otsu\n\n# Obtiene la imagen normalizada en escala de grises\ndef getGrayImg(im_path):\n    im = imageio.imread(str(im_path))\n    im_gray = rgb2gray(im)\n    if im_gray.mean() > 0.5:\n        return 1-im_gray\n    return im_gray\n#Descripción de un array\ndef describe(array, name='array', show=True, otsu=False):\n    datos=(1,1,1)\n    if otsu:\n        datos=(threshold_otsu(array), array.min(),array.max())\n    else:\n        datos=(array.mean(), array.min(),array.max())\n    if show:\n        print(name, '\\n  mean:', datos[0],'\\n  min:', datos[1],'\\n  max:', datos[2])\n    return datos\n#Obtiene un kernel ponderando por distancias\ndef getDistKernel(radio=1):\n    row=[i for i in range(-radio,radio+1)]\n    a=np.array([row for i in range(radio*2+1)])\n    b=a.transpose()\n    d= np.multiply(a,a)+np.multiply(b,b)\n    vmax=d.max()\n    d=((d/vmax)-1)*(-1)    \n    return d\n# binarización de imágenes con una imágen de referencia\ndef binarizaMedias(im_original, im_umbral,factor=0.5,otsu=False):\n    meanA,maxA,minA=describe(im_original,show=False,otsu=otsu)\n    meanB,maxB,minB=describe(im_umbral,show=False,otsu=otsu)\n    imA= ((im_original-minA)/(maxA-minA))-(meanA-minA)/(maxA-minA)\n    imB=((im_umbral-minB)/(maxB-minB))-(meanB-minB)/(maxB-minB)\n    #describe(imA,'Original')\n    #describe(imB,'Umbral')\n    return np.where(imA>imB*factor,0,255)\n# Mide la exactitud de coincidencia entre dos imagenes binarizadas\ndef mideExactitud(im_ideal, im_calculada, imagen=False, mensaje=False):\n    dif=np.zeros(im_ideal.shape,dtype=np.uint)\n    dif=im_ideal^im_calculada\n    AreaMascara=np.sum(im_calculada==255)\n    v= np.sum(dif==255)/AreaMascara\n    if mensaje:\n        print('Exactitud', 1-v)\n    if imagen:\n        return(v, dif)\n    return v\n\nprint(\"Done\")","execution_count":134,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"cad5e637281b18ca76abbe136af5d3dce806218d"},"cell_type":"code","source":"# Realiza analisis en multiples parametros para binarizar una imagen\n# se da una tupla con los pares de factor, radio del kernel\ndef analizaParametros(im_path,factores=(0.4,0.5), radios=(1,2,3,4)):\n    im_gray= getGrayImg(im_path)\n    im_mask= getMask(im_path)\n    resultados=[]\n    for radio in radios:\n        kernel=getDistKernel(radio)\n        normal_result = signal.convolve2d(im_gray,kernel, boundary='symm', mode='same')\n        for factor in factores:\n            img_bin=binarizaMedias(im_gray,normal_result,factor=factor,otsu=False)\n            img_binOtsu=binarizaMedias(im_gray,normal_result,factor=factor,otsu=True)\n            #print('Shape img_bin:',img_bin.shape, 'Shape im_mask:', im_mask.shape)\n            exact=mideExactitud(im_mask,img_bin)\n            exactOtsu=mideExactitud(im_mask,img_binOtsu)\n            resultados.append([factor,radio,exact,exactOtsu]) \n    return resultados\nprint(\"Done\")","execution_count":59,"outputs":[]},{"metadata":{"_cell_guid":"be1c69c3-66ea-43ed-a4ab-d1c7a0711320","_uuid":"9e058a0e7ecf866fd7b34f2cae9341f9f51b27a7"},"cell_type":"markdown","source":"Ahora se obtiene una imagen del perfil de umbral"},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"af23991fb4c0b7296bc90745434735218b70027f"},"cell_type":"code","source":"def getMejor(resultados, mensaje=True):\n    rMean=[ex[2]for ex in resultados]\n    rOtsu=[ex[3]for ex in resultados]\n    iMean=rMean.index(min(rMean))\n    iOtsu=rOtsu.index(min(rOtsu))\n    \n    mejorMean=(resultados[iMean][1],resultados[iMean][0] , 1-resultados[iMean][2])\n    mejorOtsu=(resultados[iOtsu][1],resultados[iOtsu][0] , 1-resultados[iOtsu][3])\n    if mensaje:\n        print('Mejor exactitud (radio,factor):\\n','\\tMean: (%i,%f) = %f \\n'%mejorMean,'\\tOtsu: (%i,%f) = %f '%mejorOtsu, )\n    return (mejorMean,mejorOtsu)\nprint('Done')","execution_count":121,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true,"_uuid":"7d95ce946e9e37ec034263068bff083361fdad02"},"cell_type":"code","source":"im_path = training_sorted[5]\nfactores=np.linspace(0.0,0.9,20)\nradios=[1,2,3,4,5,6,7,8,9,10]\n#factores=[0,0.1,0.2,0.3,0.4]\nexact= analizaParametros(im_path,factores=factores ,radios=radios)\nexactitud=[(ex[-2],ex[-1])for ex in exact]\ngetMejor(exact)\nif len(radios)>1:\n    im_exact=np.zeros((len(factores),len(radios)))\n    im_exactOtsu=np.zeros((len(factores),len(radios)))\n    for f,r,em,eo in exact:\n        ix=np.where(factores==f)\n        iy=radios.index(r)\n        im_exact[ix,iy]=1-em\n        im_exactOtsu[ix,iy]=1-eo\n    valoresEjes=[min(radios),max(radios),max(factores),min(factores)]\n    fig = plt.figure(1,figsize=(10,5))\n    ax = fig.add_subplot(1,2,1)#Impares\n    fig.colorbar(ax.imshow(im_exact,cmap='inferno',extent=valoresEjes,aspect=15))\n    ax.set_adjustable('box-forced')\n    ax.set_title('Exactitud por media')\n    plt.ylabel('factor')\n    plt.xlabel('radio')\n    ax = fig.add_subplot(1,2,2)#Impares\n    fig.colorbar(ax.imshow(im_exactOtsu,cmap='inferno',extent=valoresEjes,aspect=15))\n    #ax[1].imshow(, cmap=plt.cm.gray)\n    ax.set_title('Exactitud por Otsu')\n    ax.set_adjustable('box-forced')\n    plt.ylabel('Factor')\n    plt.xlabel('Fadio')\n    \nelse:\n    #print(pares)\n    plt.figure(figsize=(10,5))\n    plt.plot(factores,exactitud, label=['a','b'])\n    plt.grid()\n    plt.show()","execution_count":138,"outputs":[]},{"metadata":{"_uuid":"a203deecd0143fba98f465d3adbe47b4c3b39b0a"},"cell_type":"markdown","source":"Hallados los mejores parámetros, se analiza como es el resultado obtenido."},{"metadata":{"_cell_guid":"afb74250-212f-446c-89e2-f13604f8b8e8","_kg_hide-input":true,"_uuid":"645d99d00a57e5f4e9909f977fec2caed5de4b75","trusted":true},"cell_type":"code","source":"\nim_path = training_sorted[5]\ncmaps=('gnuplot','gray')\nfigsize=(10,10)\nimgray= getGrayImg(im_path)\nim_mask= getMask(im_path)\n\nradio=1\nfactor=0.7\notsu=True\nkernel=getDistKernel(radio)\nnormal_result = signal.convolve2d(imgray,kernel, boundary='symm', mode='same')\n\nimg_bin=binarizaMedias(imgray,normal_result,factor=factor,otsu=otsu)\nex,img_dif=mideExactitud(im_mask,img_bin,imagen=True)\n\nfig, axes = plt.subplots(nrows=2, ncols=2, figsize=figsize,\n                         sharex=True, sharey=True)\nax = axes.ravel()\ntitles = ['Original', 'Ideal', 'Diferencia','Binarizado']\nimgs = [imgray, im_mask, img_dif, img_bin  ]\nfor n in range(0, len(imgs)):\n    ax[n].imshow(imgs[n], cmap=plt.cm.gray)\n    ax[n].set_title(titles[n])\n    ax[n].set_adjustable('box-forced')\n    #ax[n].axis('off')\n\nplt.tight_layout()\nprint(imgray.shape, normal_result.shape)\nprint('  Factor:', factor, '  radio:', radio, '  Otsu:',otsu , ' Exactitud:',1-ex)\n\nplt.show()\n","execution_count":139,"outputs":[]},{"metadata":{"_uuid":"b758e31fb8783034c0806e925e08e08fbad18ed2"},"cell_type":"markdown","source":"## Creando una base de datos con los mejores parámetros\n"},{"metadata":{"trusted":true,"_uuid":"a53011734e0d9476de3c14dc6da2b345f978d0ef"},"cell_type":"code","source":"factores=np.linspace(0.0,0.9,20)\nradios=[1,2,3,4,5,6,7,8,9,10]\ncolumns=['Otsu','Radio','Factor','Exactitud']\ndf = pd.DataFrame(columns=columns)\nfor i,im_path in enumerate(training_sorted[:20]):\n    #factores=[0,0.1,0.2,0.3,0.4]\n    exact= analizaParametros(im_path,factores=factores ,radios=radios)\n    mejorMean,mejorOtsu=getMejor(exact, mensaje=False)\n    if mejorMean[-1]>mejorOtsu[-1]:\n        tipoRef=(0,mejorMean[0],mejorMean[1],mejorMean[2])\n    else:\n        tipoRef=(1,mejorOtsu[0],mejorOtsu[1],mejorOtsu[2])\n    df.loc[i]= tipoRef\n\npath=pathlib.Path('featuresA.csv')\ndf.to_csv(path)\ndf.head(5)","execution_count":145,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}