{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Reading these Ginormous images in R\n\nIn this notebook, we cover how to open and render the very large images. We also take a crack at subsetting the images down into smaller areas, better suited to running computer vision algorithms on.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"### First, lets take a look at our training set","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#initial package load:\nlibrary(magrittr)\n\ntrain=utils::read.csv(\"../input/prostate-cancer-grade-assessment/train.csv\", stringsAsFactors=F)\nhead(train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Entires where the `isup_grade` and `gleason_score` are equal to 0 or 0+0 respectively are graded as non-cancerous. Values above 0 in either scoring system indicate increasing 'severity'. The ISUP grading system is the easiest to understand, with simple scores ranging from 0-5 (healthy-high severity). Lets read in the first healthy and first high-severity tiff in, and take a look:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Find the image name of the first helathy image (and construct the image path)\nhealthyIndex=which(train$isup_grade==0)\nhealthy=train$image_id[healthyIndex[1]]\nhealthyPath=paste0(\"../input//prostate-cancer-grade-assessment//train_images/\", healthy, \".tiff\")\n\n# now, the first entry in the training set with a very severe ISUP grade:\nsevereIndex=which(train$isup_grade==5)\nsevere=train$image_id[severeIndex[1]]\nseverePath=paste0(\"../input//prostate-cancer-grade-assessment//train_images/\", severe, \".tiff\")\n\n# \nprint(c(healthy, severe))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Using `stack` in the raster package, we can read all 3 bands in at once. Lets do that:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"healthyRaster=raster::stack(x = healthyPath)\nsevereRaster=raster::stack(x=severePath)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"---\n**Viewing the healthy image first:**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"raster::plotRGB(healthyRaster)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"---\n**Now the severely-graded image:**","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"raster::plotRGB(severeRaster)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"These images look pretty different, but the most obvious difference is that the second image looks less saturated than the first. This is because the images are from different data sources:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"c(healthy=train$data_provider[healthyIndex[1]], severe=train$data_provider[severeIndex[1]])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's look at a severe image from the Karolinska set:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"severeK=train$image_id[train$data_provider=='karolinska' & train$isup_grade==5][1]\nsevereKPath=paste0(\"../input//prostate-cancer-grade-assessment//train_images/\", severeK, \".tiff\")\nsevereKRaster=raster::stack(severeKPath)\nraster::plotRGB(severeKRaster)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"That looks much more comparable to the first image, at least as far as the color saturation. \n> Note: We can see kinda see how the unhealthy image above differs from the healthy image in how densely packed the tissue looks.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"---\n### Sub-sampling our rasters\nThe severe image from the Karolinska set is _pretty big_:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"nCols=ncol(severeKRaster)\nnRows=nrow(severeKRaster)\nc(totalColumns=nCols, totalRows=nRows)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Note that the number of columns/rows shown above is for one channel in the raster. There are 3x as many columns and rows to examine accross all the channels. Plus, the raster is made larger by the orientation of the tissue sample. The diagonal layout means there is quite a bit more whitespace than we'd like.\n\nClearly, to work effiecintly in R, we'll need to extract areas of this raster to feed into our CV algorithm.\n\n## 100 CELLS!\n10 x 10 is good. Lets go with that. Why not?\n> Note: A reason why not is smaller or bigger cell sizes may be approproate based on the image (consider lower resolution images, or a rectangular aspect ratio)\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"cellSizeX=floor(nCols/10)\ncellSizeY=floor(nRows/10)\nprint(c(cellSizeX, cellSizeY))\n\n# This technically makes 1-pixel overlaps\nxMins=seq(from = 0, by = cellSizeX, length.out = 10)\nxMaxs=seq(from = cellSizeX, by = cellSizeX, length.out = 10)\nyMins=seq(from = 0, by = cellSizeY, length.out = 10)\nyMaxs=seq(from = cellSizeY, by = cellSizeX, length.out = 10)\n\ngridBounds=cbind(expand.grid(xMins, yMins) %>% `names<-`(value=c(\"xMin\", \"yMin\")), \n      expand.grid(xMaxs, yMaxs) %>% `names<-`(value=c(\"xMax\", \"yMax\")), \n     cellNum=seq(1:100))\n\nraster::plotRGB(severeKRaster)\nabline(h = gridBounds$yMin, v=gridBounds$xMin)\n\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The grid pattern above shows roughly how our 100 cells are laid out on the image. Obviously, most of the cells are white, and we can ignore them. Let's have a funciton to check that the cell contains a color other than white:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# given a raster and the x/y extent of the raster you want to check,\n# this function returns a logical TRUE/FALSE for if every pixel is white.\nallWhite=function(raster, xMin, xMax, yMin, yMax){\n    # Kinda fast return of values with extract\n    cell=raster::extract(x=raster, y=raster::extent(c(xMin, xMax, yMin, yMax)))\n    #cell=raster::crop(x=raster, y=raster::extent(c(xMin, xMax, yMin, yMax)))\n    \n    # Check- is everything in cell white\n    # This check is relatively fast\n    allWhite=all(cell==255)\n    return(allWhite)\n}\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we have our function, we can check each cell we've defined for it's conctents, and only keep cells with some color in them. \n\n## _NOTE this is kinda slow. A few minutes..._\n**Call your mother while it runs.** Or spend some time pondering how to make it faster.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#Initialize an empty vector to add to\nkeepCell=c()\n\n# Loop thru the grid cells, returning T/F for if the cell is all white\nfor(r in 1:nrow(gridBounds)){\n    #print(r)\n    #s=Sys.time()\n    keepCell=c(keepCell,\n              !allWhite(raster=severeKRaster, \n                       xMin = gridBounds$xMin[r],\n                       xMax = gridBounds$xMax[r],\n                       yMin = gridBounds$yMin[r],\n                       yMax = gridBounds$yMax[r])\n              )\n    #e=Sys.time()\n    #message(difftime(e, s, \"secs\"))\n}\n# Note:\n# This block takes ~250 secs (4 min) to run on a large image. \n# The extract function in raster is faster than crop, but still not that fast","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**NOW WE CAN CONFIDENTLY ZOOM TO TISSUE ONLY REGIONS FOR ANALYSIS.**  \n\nLets store the 'good' cells in `tissueCells` (confusing name?), and take a look what is in grid cell 35","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"tissueCells=gridBounds[keepCell,]\n# Row number 8 has the info for grid cell #35\nraster::plotRGB(raster::crop(x=severeKRaster, y=raster::extent(c(tissueCells$xMin[8], \n                                                  tissueCells$xMax[8], \n                                                  tissueCells$yMin[8], \n                                                  tissueCells$yMax[8]))))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# WE DID IT\n![ma](https://www.zerohedge.com/sites/default/files/images/user5/imageroot/mission%20accomplished%202.jpg)\n\nOnto a new notebook where we devise a CV grading algorithm...\n","execution_count":null}],"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.6.0","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":4}