{"cells":[{"metadata":{},"cell_type":"markdown","source":"# OSIC Pulmonary Fibrosis Progression: EDA - An R perspective.\nIn this competition I am starting to see the young sprouts of R candidates building models and trying to compete. Sadly one quite quickly hits a bump in the road when it comes to importing images from R. As most would be aware a few of these images are spectecled with some problems, which causes [`pydicom`](https://pydicom.github.io/) in python and [`oro.dicom`](https://cran.r-project.org/package=oro.dicom) in R to crash, and several a few images have [missing features](https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/165723#930391) such as slice position. The latter can be fixed or incorporated into the model procedure itself, but for a complete model one will have to analyse the scans themselves where the first step is being able to import every single one into R.\n\nThis is the purpose of this notebook, to give a method to import **all** the ct-scans, and thus be able to further analyse these in a familiar (cough superior) environment. \n\nThe first problem becomes apparent quite quickly when reading the (new dataset discussion)[https://www.kaggle.com/c/osic-pulmonary-fibrosis-progression/discussion/165723] posted about a month ago. `pydicom` contains a binding to the `gdcm` python package (which wraps around the [`gdcm`](http://gdcm.sourceforge.net/) c++ library) while `oro.dicom` does not. I have not got the time (nor the interest) in implement a wrapper through [`Rcpp`](http://www.rcpp.org/) to `gdcm`, and as such a natural thought is to use `pydicom` and `gdcm` through the [`reticulate`](https://rstudio.github.io/reticulate/) package. However if one tries to execute\n\n```R\nlibrary(reticulate)\npy_install('gdcm') \npy_install('pydicom')\n```\nreturns the fantasticly unhelpful error message:\n\n>Error: Error installing package(s): 'gdcm'\n>Traceback:\n>\n>1. py_install(\"gdcm\")\n>2. virtualenv_install(envname = envname, packages = packages, ...)\n>3. pip_install(python, packages, ignore_installed = ignore_installed, \n> .     pip_options = pip_options)\n>4. stop(msg, call. = FALSE)\n\nSo our very first problem becomes getting acces to `gdcm` . Luckily this doesn't become too complex, and a very helpful post on [`pydicom`s github page](https://github.com/pydicom/pydicom/wiki/Installing-the-Python-GDCM-bindings-without-Conda) gives us an ubuntu specific way of installing the package.\n\n```bash\nsudo apt install python3-gdcm\n```\n\nSuch commands can be executes from R using `system` (or simpler) `system2` functions. And in this way we can gain access to all of the nice functionality from `pydicom` library.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"library(reticulate)\nres <- system2('sudo', c('apt', 'install', 'python3-gdcm'))\npy_install('pydicom')\npydicom = import('pydicom')\nlibrary(oro.dicom)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that we've got the packages, the next step search our scan files for images that cannot be imported by `oro.dicom` but may be imported by `pydicom` with `gdcm`. For this I'll use a simple `for-loop`.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"path <- '../input/osic-pulmonary-fibrosis-progression/'\nfiles <- list.files(path, pattern = '*.dcm', full.names = TRUE, recursive = TRUE)\nlength(files)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Faster alternative: Read only 1 file per directory.\n#allFiles <- files\n#dirs <- c(list.dirs('../input/osic-pulmonary-fibrosis-progression/test/', full.names = TRUE),\n#          list.dirs('../input/osic-pulmonary-fibrosis-progression/train/', full.names = TRUE))\n#files <- sapply(dirs, function(x){list.files(x, recursive = TRUE, full.names = TRUE)[1]})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"readAll <- compiler::cmpfun(\n    function(f){\n        out <- vector('list', length(f))\n        names(out) <- f\n        for(z in seq_along(f)){\n            i <- f[z]\n            j <- try(readDICOM(i), silent = TRUE)\n            if(inherits(j, 'try-error')){\n                j <- try(pydicom$dcmread(i))\n                if(inherits(j, 'try-error')){\n                    out[[i]] <- 'error'\n                }else\n                    out[[i]] <- 'pydicom'\n            }else\n                out[[i]] <- 'oro.dicom'\n        }\n        out\n    }\n)\nreaders <- readAll(files)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This takes quite a while so in order to avoid doing this twice, we'll save the result.\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"library(data.table)\ndf <- data.table(file = names(readers), reader = unlist(readers))\nfwrite(df, 'readers.csv')\ndf[, .N, by = reader]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we see only a few files require us to read through `pydicom` which makes our life simpler. \n\nNow that we know which images require which reader, next we have to \n\n1. Import the standout images into R\n1. Extract the image matrix itself\n1. Extract all attributes into the same format as returned by `oro.dicom`\n\n\nImporting the images and extracting the image matrix is a simple job. \n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pydicom_file <- df[reader == 'pydicom'][1, file]\npydicom_class <- pydicom$dcmread(pydicom_file)\npydicom_image_matrix <- pydicom_class$pixel_array\npydicom_image_matrix[1:6, 1:6]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"oro.dicom_file <- df[reader == 'oro.dicom'][1, file]\noro.dicom_image <- readDICOM(oro.dicom_file)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Our real problem comes when we want to import the attributes contained within each image.\n\nFor this we need to extract\n1. The attribute value \n1. The attribute binary identifier eg. (0002, 0000)\n1. The attribute descriptor eg. File Meta Information Group Length\n1. The attribute descriptor (short) eg. UL\n\nNow all of these values can be extracted seperately, but are hidden underneath the hood of the `pydicom` dicom class.\n\nSo for illustration our goal is to achieve something similar to (but with varying fields)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"head(oro.dicom_image$hdr[[1]], 2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There is (to my knowledge) no truely simple way to extract all this information from the python object itself. \n\nHere I simply default to a loop over the `pydicom` object fields, which I came up with after some regretting even starting the process. The result of my personal self-imposed torture is composed in the extract_attributes function below. Note the `make_exact` argument decides whether to return the attributes in a list (values of different classes) or as a `data.frame` (similar to `oro.dicom`).  \n\n\nI'd imagine that using `pydicom` for all images and just keeping the attributes as a list (or tibble) would simplify some pre-processing. \n\nAs a side note I've left out private members, those not contained in `pydicom_class$dir()`, but these only seem to contain only irrelevant fields for this analysis.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Required packages and R stuff.\nbuiltins <- import_builtins()\nlibrary(stringr)\nextract_attributes <- compiler::cmpfun(function(pydicom_im, make_exact = TRUE, builtins = reticulate::import_builtins()){\n    # Change to iterate only over values in pydicom_im$dir()\n    # Seems private members are contained in values.\n    fields <- pydicom_im$dir()\n    n <- length(fields)\n    value <- vector('list', n)\n    sequence <- code <- name <- element <- group <- character(n)\n    length <- numeric(n)\n    convVal <- function(val){\n        clv <- class(val)\n        if(is.raw(val))\n            value <- rawToChar(val)\n        else if('python.builtin.float' %in% clv)\n            value <- builtins$float(val)\n        else if('python.builtin.bytes' %in% clv) # Bytes fuck us over quite a bit if they contain a NULL string \\\\x00, and these are terrible to remove.. Here's a quite manual way of fixing it.\n            value <- gsub('^(b\\')|(\\')$', '', py_to_r(r_to_py(builtins$str(val))$replace('\\\\x00', '')))\n        else if('python.builtin.int' %in% clv)\n            value <- builtins$int(val)\n        else if('python.builtin.str' %in% clv || grepl('pydicom.valuerep', clv))\n            value <- builtins$str(val)\n        else if('collections.abc.Sequence' %in% clv)\n            value <- lapply(builtins$list(val), convVal)\n        else \n            value <- val\n        value\n    }\n    for(i in seq_len(n)){\n        field_i <- pydicom_im$get_item(fields[i])\n        tag_i <- field_i$tag\n        group[i] <- tag_i$group\n        element[i] <- tag_i$element\n        name[i] <- fields[i]\n        code[i] <- field_i$VR\n        sequence[i] <- ''\n        length[i] <- tag_i$bit_length()\n        value[[i]] <- convVal(field_i$value)\n    }\n    if('PixelData' %in% name)\n        value[[which('PixelData' == name)]] <- 'PixelData'\n    out <- list(group = str_pad(group, 4, 'left', '0')\n                      , element = str_pad(element, 4, 'left', '0')\n                      , name = name\n                      , code = code\n                      , length = length\n                      , value = value\n                      , sequence = rep('', n)\n    )\n    if(make_exact){\n        # Convert values into same format as oro.dicom\n        lists <- which(lengths(out[['value']]) > 1)\n        out[['value']][lists] <- lapply(out[['value']][lists], paste0, collapse = ' ')\n        out[['value']] <- unlist(out[['value']])\n        out <- as.data.frame(out)[order(out[['group']], out[['element']]), ]\n        rownames(out) <- seq_len(n)\n    }\n    out\n})","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pydicom_attributes <- extract_attributes(pydicom_class, builtins = builtins)\nhead(pydicom_attributes, 2)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now to obtain the exact same result as oro.dicom we have to combine the attributes and image into a nested list, and we can obtain this using","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pydicom_image <- list(hdr = pydicom_attributes, img = list(pydicom_image_matrix))\nnames(pydicom_image[2]) <- pydicom_file","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Finally we finish this up with conveniently wrapping this all into a single function, that can be called for any of the images within the competition.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"#' Wrapper for oro.dicom and pydicom to read \n#'\n#' @param path path for the file that should be imported\n#' @param ... extra arguments pass to oro.dicom::readDICOM\n#' @param method method for importing the dcm file (auto/oro.dicom/pydicom)\n#' @param make_exact for method == 'pydicom', should the output be returned as oro.dicom format?\n#' @param pyDicom binding to the imported pyDicom library. gdcm is assumed to be installed.\n#'\n#' @return a list with fields hdr and img.\nreadDICOM_V2 <- compiler::cmpfun(function(path, ..., method = 'auto', make_exact = TRUE, builtins = NULL, pyDicom = NULL){\n    if(!is.character(method) || !(method <- tolower(method)) %in% c('auto', 'pydicom', 'oro.dicom') )\n        rlang::abort('Method not in c(\"auto\", \"oro.dicom\", \"pydicom\")')\n    if(method == 'auto'){\n        # Hacky way to implement an auto method. Prolly should do this properly, but I can't be bothered.\n        out <- try(readDICOM_V2(path, ..., method = 'oro.dicom'), silent = TRUE)\n        if(inherits(out, 'try-error'))\n            method <- 'pydicom'\n    }\n    if(method %in% 'oro.dicom'){\n        out <- readDICOM(path, ...)\n    }else if(method == 'pydicom'){\n        pyDicomMissing <- !is.null(pyDicom) && ('python.builtin.module' %in% class(pyDicom)) && (pyDicom$'__name__' != 'pydicom')\n        if(pyDicomMissing)\n            rlang::abort('pyDicom missing while trying to import image. Use reticulate::import(pydicom) as argument to pyDicom, and make sure the module is installed.')\n        \n        out <- pyDicom$dcmread(path)\n        arr <- out$pixel_array\n        hdr <- extract_attributes(out, make_exact = make_exact, builtins = builtins)\n        out <- list(hdr = hdr, img = list(arr))\n        names(out$img) <- path\n    }\n    out\n})\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And last we can test that the function works as intended.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"pydicom$dcmread(pydicom_file)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"py_read <- readDICOM_V2(pydicom_file, method = 'pydicom', pyDicom = pydicom, builtins = builtins)\noro_read <- readDICOM_V2(oro.dicom_file)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As an end note I'd like to mention, that this is definitely not the most efficient way to store the data for analysis purposes. As can be seen this notebook took about 1 hour to complete, and this is what one would say 'not viable' for batched learning purposes. It is likely a better method is to import the data, and save it in a more convenient and faster format such as `parquet` or `fst` files, the latter likely being the fastest alternative. ","execution_count":null}],"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"3.6.3"}},"nbformat":4,"nbformat_minor":4}