{"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":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Executive Summary\n\nIt was found ................ from results and conclusion.\n\n# Objectives\n\nThis model serves two objectives. \n\n1. Build and test an image classification model - able to predict blood clot origins from ischeamic stroke pathology images.\n2. Provide evidence of competency for an academic capstone project.","metadata":{}},{"cell_type":"markdown","source":"# Setup\n\nWorking with images and machine learning in the R Programming language - the following packages will be of use:\n\n- Tidyverse, goes without saying.\n- Magick, image processing library.\n- Keras, Nueral Net library.\n- Caret, Comprehensive Machine Learning library.\n- Tiff, manipulation of tiff images. \n- Raster, image analysis with layering.\n- RcolorBrewer, for pretty plots.\n- jpeg, useful for image compression.","metadata":{}},{"cell_type":"code","source":"#Online Package Setup.\n\n#if(!require(tidyverse)) install.packages(\"tidyverse\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(magick)) install.packages(\"magick\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(keras)) install.packages(\"keras\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(caret)) install.packages(\"caret\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(tiff)) install.packages(\"tiff\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(raster)) install.packages(\"raster\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(RColorBrewer)) install.packages(\"RColorBrewer\", repos = \"http://cran.us.r-project.org\")\n\n#if(!require(BiocManager)) install.packages(\"BiocManager\")\n\n#if(!require(jpeg)) install.packages(\"jpeg\", repos = \"http://cran.us.r-project.org\")","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:30.071861Z","iopub.execute_input":"2022-10-08T23:50:30.074937Z","iopub.status.idle":"2022-10-08T23:50:30.117283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Meta: Analysis | Method\n\nPrior to committing to any particular course of action. Let's first observe and analyse the meta data attached to our images.","metadata":{}},{"cell_type":"code","source":"#First lets look at the meta.data\n\nlibrary(tidyverse)\n\ntrain_meta <- read_csv(\"../input/mayo-clinic-strip-ai/train.csv\")\n\ntrain_meta <- train_meta %>% data.frame()\n\nstr(train_meta)\n\ntest_meta <- read_csv(\"../input/mayo-clinic-strip-ai/test.csv\")\n\ntest_meta <- test_meta %>% data.frame()\n\nstr(test_meta)\n\n#Top and tail observation.\n\nhead(train_meta)\n\ntail(train_meta)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:30.121062Z","iopub.execute_input":"2022-10-08T23:50:30.161312Z","iopub.status.idle":"2022-10-08T23:50:31.47515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Variables:\n\n- image_id - unique TIFF identifier.\n- center_id - Hospital or Lab where clot was reported.\n- patient_id - unique patient identifier.\n- image_num - if multiple images taken of the same clot. This is the sequence number [starting at 0]\n- label - clot aetiology (if known) provided from medical analysis: LAA = Large Artery Atherosclerosis and CE = Cardioembolic.","metadata":{}},{"cell_type":"markdown","source":"## Visualization\n\nLet's visualize the distributions of Aetiology (label) and Location (center_id).","metadata":{}},{"cell_type":"code","source":"#Distribution of Aetiology and Location.\n\nbar1 <- as.data.frame(table(train_meta$label))\n\ncolnames(bar1) <- c(\"Aetiology\",\"Frequency\")\n\nbar2 <- as.data.frame(table(train_meta$center_id))\n\ncolnames(bar2) <- c(\"Center_ID\", \"Frequency\")\n\nplot <- function(x){\n    ggplot(x, aes(x = x[,1], y = x[,2])) + \n    geom_bar(stat = \"identity\") +\n    xlab(colnames(x[1])) +\n    ylab(colnames(x[2]))\n    }\n\nplot(bar1)\n\nplot(bar2)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:31.478668Z","iopub.execute_input":"2022-10-08T23:50:31.480611Z","iopub.status.idle":"2022-10-08T23:50:32.094214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Relative Proportions.\n\nProp <- function(tbl){\n    \n    Prop <- data.frame(Class = tbl[,1], Proportion = tbl[,2]/sum(tbl[,2]))\n    \n    Prop %>% arrange(desc(tbl[,2]))\n}\n\nProp(bar1)\n\nProp(bar2)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.097995Z","iopub.execute_input":"2022-10-08T23:50:32.099605Z","iopub.status.idle":"2022-10-08T23:50:32.156481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Let's evaluate the proportion of Aetiology/labels for each hospital.\n\nbar3 <- train_meta %>% group_by(center_id) %>% \nsummarize(Total = n(), LAA = sum(label == \"LAA\"), CE = sum(label == \"CE\")) %>% \npivot_longer(cols = LAA:CE, names_to = \"Aetiology\", values_to = \"Count\")\n\nbar3 %>% head()\n\nbar3 %>% ggplot(aes(x = center_id, y = Count, fill = Aetiology)) +\n       geom_bar(stat = \"identity\")\n","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.160037Z","iopub.execute_input":"2022-10-08T23:50:32.16165Z","iopub.status.idle":"2022-10-08T23:50:32.508114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#facet\n\nbar3 %>% group_by(Aetiology) %>% \nggplot(aes(x = center_id, y = Count)) +\n       geom_bar(stat = \"identity\") +\nfacet_grid(~Aetiology)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.511967Z","iopub.execute_input":"2022-10-08T23:50:32.513556Z","iopub.status.idle":"2022-10-08T23:50:32.804934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Meta Data Summary\n\nThrough the analysis and visualization of our training meta data the following can be inferred:\n\n- Cardioembolic clots make up 73% of the total observations within the training data set. (Large Artery Atherosclerosis - 27%) \n- Cardioembolic clots have the majority across all centers, with the exception being center_id - 3.\n- Center_id = 11. Had the greatest magnitude of observations submitted (33% of the total).","metadata":{}},{"cell_type":"markdown","source":"## Method: Meta Model","metadata":{}},{"cell_type":"markdown","source":"If one was to derive a simplistic classification model from meta data alone. Predictions could be made as to the Aetiology of clots based off the following algorithm.\n\n---\n\n$$ With \\: CE = 1 \\: and \\: LAA = 0 $$ \n\n$$ Prediction = IF\\: (Pr(CE \\: | \\: Center Id)\\: \\geq .5),\\: THEN \\:1 ,\\: ELSE \\: 0 $$\n\n---\n\n*NOTE - this is an overly simplistic model, however parameters could be added (naively) in an iterative fashion to increase accuracy and decrease variability.*","metadata":{}},{"cell_type":"code","source":"#Method: Step 1.\n\n#Evaluate the proportion of observed CE clots for given center_id x.\n\ncenter_bias <- train_meta %>% group_by(center_id) %>% \nsummarize(CE_prop = sum(label == \"CE\")/n())\n\ncenter_bias","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.809015Z","iopub.execute_input":"2022-10-08T23:50:32.810649Z","iopub.status.idle":"2022-10-08T23:50:32.846244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Method: Step 2.\n\n#Join the above proportions to our test data.\n\n#Utilize an IFELSE logic statement to evaluate proportion against threshold (.5).\n\n#return CE if greater than or equal to threshold, else return LAA.\n\nPred_Meta <- test_meta %>% \nleft_join(center_bias, by = \"center_id\") %>%\nmutate(Pred_Meta = ifelse(CE_prop >= .5, \"CE\", \"LAA\"))\n\nPred_Meta","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.849847Z","iopub.execute_input":"2022-10-08T23:50:32.851693Z","iopub.status.idle":"2022-10-08T23:50:32.888276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Meta Model Results\n\nThe meta model has returned a prediction that each test observation is 69% likely to be of cardioembolic (CE) origin. This is due to the fact that our test observations are derived from the same center_id. For brevity and to complete objective 2, i will submit the meta model and continue to work on the image classification model in slow time. ","metadata":{}},{"cell_type":"markdown","source":"## Meta Model Submission","metadata":{}},{"cell_type":"code","source":"#Read in test csv.. and assign to DF.\n\ntest_sub <- read.csv('../input/mayo-clinic-strip-ai/test.csv')\n\ntest_sub","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.892443Z","iopub.execute_input":"2022-10-08T23:50:32.895287Z","iopub.status.idle":"2022-10-08T23:50:32.92995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Join training CE proportions based for each center_id in test data.\n\ntest_sub <- test_sub %>% left_join(center_bias, by = \"center_id\") %>%\nmutate(CE = CE_prop, LAA = 1 - CE_prop)\n\ntest_sub\n\n#If a new center_id is provided within the test data set, NAs may be produced and will need to be replaced with default values.\n\ntest_sub$CE <- test_sub$CE %>% replace_na(.7)\n\ntest_sub$LAA <- test_sub$LAA %>% replace_na(.3)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.932314Z","iopub.execute_input":"2022-10-08T23:50:32.933689Z","iopub.status.idle":"2022-10-08T23:50:32.968269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Remove redundant variables\n\ntest_sub <- test_sub[,c(\"patient_id\",\"CE\",\"LAA\")]\n\n#Remove duplicate patient IDs\n\ntest_sub <- test_sub[!duplicated(test_sub),]\n\ntest_sub","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:32.970618Z","iopub.execute_input":"2022-10-08T23:50:32.972038Z","iopub.status.idle":"2022-10-08T23:50:32.997678Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Write to output file \"submission.csv\"\n\nwrite_delim(test_sub, file='submission.csv', delim=',')\n\nsessionInfo()","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:33.000076Z","iopub.execute_input":"2022-10-08T23:50:33.00146Z","iopub.status.idle":"2022-10-08T23:50:33.054019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Image Classification: Analysis | Model","metadata":{}},{"cell_type":"markdown","source":"## Image Analysis","metadata":{}},{"cell_type":"code","source":"library(EBImage)\n\nlibrary(raster)\n\nlibrary(tiff)\n\nlibrary(jpeg)\n\n#To reduce memory consumption of large rasters, i will change a setting by which, writing data to disk as opposed to storing in memory.\n\nrasterOptions(todisk = TRUE)\n\n#Sample tiff image.\n\ntrain_file_list <- list.files(\"../input/mayo-clinic-strip-ai/train\")\n\nsample <- as.character(sample(train_file_list,size = 1))\n\nsample_path <- paste(\"../input/mayo-clinic-strip-ai/train/\",sample, sep = \"\")\n\n#Visualise a sample Tiff in train data.set\n\nimage_sample_raster <- raster(sample_path, band = 1)\n\nimage_sample_raster","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:33.056603Z","iopub.execute_input":"2022-10-08T23:50:33.058125Z","iopub.status.idle":"2022-10-08T23:50:38.891402Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Visualize image\n\nimage(image_sample_raster, zlim = c(0,250))\n\nhist(values(image_sample_raster))","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:50:38.893962Z","iopub.execute_input":"2022-10-08T23:50:38.895397Z","iopub.status.idle":"2022-10-08T23:51:03.744779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Image Pre-Processing","metadata":{}},{"cell_type":"markdown","source":"Since we are dealing with very large and dimensionally varied tiff files. We will need to standardize the images to a single format, which can then be ingested by our image classification model. The model will most probably be a nueral network accepting image / tile vectors with dimensions 256 x 256.\n\nThis process is called image 'pre-processing' and will look something like this:\n\n- Read in TIFF image as raster.\n- Aggregate, if required.\n- Initial resample of raster to standard template, avoiding unacceptable quality loss.\n- Crop to grab single clot from image, if required.\n- Confirm image quality.\n- Vectorize image and add to training data frame.\n- Write resampled image to output directory with the same unique identifier + \"processed\", reducing memory cost.\n- Loop the above for all training TIFFs.","metadata":{}},{"cell_type":"code","source":"#Aggregate to decrease prominence of the mode value; prior to resampling.\n\nimage_sample_aggregated <- aggregate(image_sample_raster, fact = 5, fun = mean, na.rm = TRUE)\n\nimage(image_sample_aggregated, zlim = c(0,250))\n\nimage_sample_aggregated\n\nhist(image_sample_aggregated)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:51:03.748639Z","iopub.execute_input":"2022-10-08T23:51:03.750144Z","iopub.status.idle":"2022-10-08T23:52:27.547209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Develop function to binarize raster values.\n\n# 0 > x <= 245, 1 else 0\n\nset.seed(2)\n\nbinarize <- function (x) {\n    vec <- as.vector(x)\n    print(vec[1:5])\n    vec <- replace_na(vec, 0)\n    vec[vec >= 245] <- 0\n    print(vec[1:5])\n    vec <- ifelse(vec == 0, 0, 1)\n    print(vec[1:5])    \n}\n\n#Testing of the function\n\ntest_vec <- sample(240:250, 50, replace = TRUE)\n\ntest_values <- binarize(test_vec)\n\nstr(test_values)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:52:27.549543Z","iopub.execute_input":"2022-10-08T23:52:27.55093Z","iopub.status.idle":"2022-10-08T23:52:27.579564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Refine and incorporate binarize function\n\nbinarize_inc <- function (x) {\n    vec <- as.vector(x)\n    vec <- replace_na(vec, 0)\n    vec[vec > 245] <- 0\n    vec <- ifelse(vec == 0, 0, 1)  \n}","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:52:27.581925Z","iopub.execute_input":"2022-10-08T23:52:27.583281Z","iopub.status.idle":"2022-10-08T23:52:27.595692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Validate binarize function on image sample..\n\nvals <- binarize_inc(values(image_sample_aggregated))\n\nmat <- matrix(vals, nrow(image_sample_aggregated), ncol(image_sample_aggregated), byrow = TRUE)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:52:27.598043Z","iopub.execute_input":"2022-10-08T23:52:27.599374Z","iopub.status.idle":"2022-10-08T23:52:28.362965Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Plot binarized values.\n\nimage(mat)","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:52:28.366374Z","iopub.execute_input":"2022-10-08T23:52:28.367851Z","iopub.status.idle":"2022-10-08T23:52:43.851783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Initial resample...\n\nr1 <- raster(nrows = 1000, ncols = 1000, xmn = 0, xmx = 1000, ymn = 0, ymx = 1000)\n\nr1_1 <- resample(image_sample_aggregated, r1, method = \"ngb\")\n\nr1_1","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:52:43.85403Z","iopub.execute_input":"2022-10-08T23:52:43.855414Z","iopub.status.idle":"2022-10-08T23:52:47.48018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Second Resample to fixed dimensions 256 x 256.\n\nr2 <- raster(nrows = 256, ncols = 256, xmn = 0, xmx = 256, ymn = 0, ymx = 256)\n\nr3 <- resample(image_sample_raster, r2, method = \"bilinear\")\n\nr3\n\nextent(r3)\n\nextent(image_sample_raster)\n\n#Vectorize\n\nvect <- values(r3)\n\nvect[1:20]\n\n#Generate matrix from vector\n\nmatrx <- matrix(vect, 256, 256, byrow = TRUE)\n\n#plot matrx\n\nimage(matrx)\n        #Step XX - Confirm file size.\n\n    #Step X - Identify and isolate clots in images.","metadata":{"execution":{"iopub.status.busy":"2022-10-08T23:54:30.318198Z","iopub.execute_input":"2022-10-08T23:54:30.319975Z","iopub.status.idle":"2022-10-08T23:54:32.202534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"    #Step X - Confirm image can be represented as a Matrix (R m x n). M = total pixel width (226) and N = total pixel height (226) = Total of 50625 pixels.","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Example of process - utilising a randomly generated image with 10w x 15h pixels.\n\nRandom_vector <- rnorm(150)\n\nRandom_image <- matrix(Random_vector, nrow = 10, byrow = FALSE)\n\nRandom_image[2,] <- 5\n\nRandom_image[8,] <- 5\n\nRandom_image[,5] <- 5\n\nRandom_image[,12] <- 5\n\ndim(Random_image)\n\nRandom_image\n\nimage(Random_image)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"    #Step X - Vectorize Image for Neural Net. (Later building the train_image_matrix) \n\n#Image/Matrix requires \"unrolling/flattening\". Converting the Matrix into a N dimensional vector. Each image vector can then be fed to a Neural Nets Input Layer. \n\nImage_vector <- as.numeric(Random_image)\n\n#Confirm the original random image vector matches the vector produced from the unrolling process.\n\nConfirmation <- data.frame(Image = Random_vector, Unrolled = Image_vector)\n\nhead(Confirmation)\n\nhead(Confirmation[,1] == Confirmation[,2])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Image_vector","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}