auto_crop.R

crop <- function(file_list, img_path, crop_method = "regular")
{
  # input :

  # output : a bunch of output jpegs? Or save them all?



  #BiocManager::install("EBImage")
  library(EBImage)
  cell_count <- 0
  image_count <-0
  pair <- 0

  ## for each image that is *-dna.jpeg,
  for (file in file_list){
    setwd(img_path)
    if(grepl("*dna.jpeg$", file)){
      file_dna = file
      image_count <- image_count +1
      image <- readImage(file_dna)
      img_orig <- channel(2*image, "grey")
      pair <- 0
    }
    if(grepl("*foci.jpeg$", file)){
      file_foci = file
      image <- readImage(file_foci)
      img_orig_foci <- channel(image, "gray")
      # call functions: get
      pair <- 1
    }
    if(pair ==1){

      #### function: blur the image
      ## call it on img_orig, optional offset
      blob_th <- get_blobs(img_orig)

      blob_label = bwlabel(blob_th)
      blob_label <- channel(blob_label, "gray")
      candidate <- bwlabel(blob_label)

      ## function: remove things that aren't cells

      removed <- keep_cells(candidate)

      ### crop over each cell
      ## function: crop around every single viable cell
      ### crop foci channel here
      ## Loop over all objects in this final "removed" image


      x_final <- computeFeatures.shape(removed)
      x_final <- data.frame(x_final)
      OOI_final <- width(x_final)

      counter_final <- 0
      # looping through each object to crop
      while(counter_final<OOI_final){
        counter_final <- counter_final+1


        cell_count <- cell_count +1
        crop_single_object(removed,OOI_final,counter_final,img_orig,img_orig_foci,file_dna,cell_count)

        print("cell count for above crop is")
        print(cell_count)
      }
    }
    pair <- 0
  }

print("out of")
print(image_count)
print("we got")
print(cell_count)
print("viable cells")
}



### add all the functions used here

#################################### new function ####################################
get_blobs <- function(img_orig, crop_method = "regular"){

  # input:

  # output:

  ## subfunction: get signals and make BW mask
  ### default offset
  nucBadThresh <- 10*img_orig
  # subfunction: big blur to blobs
  img_tmp_dna <- img_orig
  img_tmp <- nucBadThresh
  w = makeBrush(size = 51, shape = 'gaussian', sigma = 15)
  img_flo = filter2(img_tmp, w)
  ## default amplification
  bg <- mean(10*img_tmp)
  offset = 0.2
  if(crop_method == "regular"){
    blob_th = 10*img_flo > bg + offset
  }

  if(crop_method == "watershed"){
    blob_th = 10*img_flo > bg + offset
  }
  return(blob_th)

}

#################################### new function ####################################
keep_cells <- function(candidate){

  # input:

  # output:



  # delete everything that's too small
  colorimg<- colorLabels(candidate, normalize = TRUE)
  x <- computeFeatures.shape(candidate)
  x <- data.frame(x)
  OOI <- width(x)
  counter <- 0
  removed <- candidate

  while(counter<OOI){
    counter <- counter+1
    pixel_area = x$s.area[counter]
    semi_maj <- x$s.radius.max[counter]
    semi_min <- x$s.radius.min[counter]
    # if statement checking if it's the wrong area
    if(pixel_area> 20000 | pixel_area < 5000){
      removed <- as.numeric(removed)*rmObjects(candidate, counter, reenumerate = TRUE)
    }
    ## if statement checking that it's not too long i.e. not at edge.
    if(semi_maj/semi_min > 2 & is.na(semi_maj/semi_min)==FALSE){
      removed <- as.numeric(removed)*rmObjects(candidate, counter, reenumerate = TRUE)
    }

  }
  removed <- bwlabel(removed)
  return(removed)
}


crop_single_object <- function(removed, OOI_final,counter_final,img_orig,img_orig_foci,file,cell_count){
  tmp_img <- removed
  ## have a single object
  ### delete all other objects
  counter_single <- 0
  # looping over all other objects to crop
  while(counter_single < OOI_final){
    counter_single <- counter_single + 1
    # iteratively remove all other objects
    if(counter_single != counter_final){
      tmp_img <- as.numeric(tmp_img)*rmObjects(bwlabel(removed), counter_single, reenumerate = TRUE)

    }

  }

  ## function: remove noise
  noise_gone <- bwlabel(tmp_img)*as.matrix(img_orig)
  noise_gone_foci <- bwlabel(tmp_img)*as.matrix(img_orig_foci)
  ## first get the row and column list that has a one in it.
  row_list <- c()
  col_list <- c()
  # I think this is quick enough for now.. takes less than 10s...
  xx = data.frame(as.numeric(tmp_img))
  xx <- data.frame(bwlabel(tmp_img))
  xm2 = t(as.matrix(xx))
  i <- 0
  my_matrix <- xm2
  ### now loop over matrix
  for(row in 1:nrow(my_matrix)) {
    for(col in 1:ncol(my_matrix)) {
      if(my_matrix[row, col]==1){
        row_list[i] <- row
        col_list[i] <- col
        i <- i+1
      }
    }
  }

  if (length(col_list>1)){
    cy <- mean(row_list)
    cx <- mean(col_list)

    x_maj <- max(col_list)
    x_min <- min(col_list)

    y_maj <- max(row_list)
    y_min <- min(row_list)
    ## radius
    ## total number of pixels given by list length. Find radius of area.
    ## crops to a square for the moment. can change this.
    max_r <- max(x_maj-cx,y_maj-cy)
    crop_r <- floor(max_r)
    # might want to do this as a matrix
    top_left_x <- floor(cx-crop_r)
    top_left_y <- floor(cy-crop_r)

    bottom_left_x <- floor(cx-crop_r)
    bottom_left_y <- floor(cy+crop_r)

    bottom_right_x <- floor(cx+crop_r)
    bottom_right_y <- floor(cy+crop_r)

    top_right_x <- floor(cx-crop_r)
    top_right_y <- floor(cy+crop_r)

    ## crop image
    ix <- bottom_left_x:bottom_right_x
    iy <- top_left_y:bottom_left_y

    # determine the dimensions, 2 in this case

    ## cropping part
    tryCatch({
      new_img <- noise_gone[ix, iy]
      ## want all images to have the same mean to 0.1
      orig_mean <- mean(new_img)
      mean_factor <- 0.08/orig_mean
      new_img <- new_img*mean_factor
      filename_crop = paste0("./crops/", file,"-crop",cell_count,"-dna.jpeg")
      writeImage(new_img, filename_crop)

      new_img_foci <- noise_gone_foci[ix, iy]
      filename_crop_foci = paste0("./crops/", file,"-crop",cell_count,"-foci.jpeg")
      writeImage(new_img_foci, filename_crop_foci)


      #### strand related stuff here
    },
    error = function(e) {
      #what should be done in case of exception?
      str(e) # #prints structure of exception
      print("couldn't crop it")

    }
    )
  }

}