Skip to content
Snippets Groups Projects
count_foci.R 3.26 KiB
Newer Older
count_foci <- function(file_list, img_path)
{
  # input :
  
  # output : a bunch of output jpegs? Or save them all?
  
  
  
  #BiocManager::install("EBImage")
  library(EBImage)
  cell_count <- 0
  image_count <-0
  pair <- 0
  foci_counts <- 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(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){
      
      #print("I have a pair of images")
      #display(img_orig)
      #display(img_orig_foci)
      
      new_img<-img_orig
      #### now see which have the right amount of strands 
      disc = makeBrush(21, "disc")
      disc = disc / sum(disc)
      localBackground = filter2(new_img, disc)
      offset = 0.1
      nucBadThresh_crop = (new_img - localBackground > offset)
      strands <- bwlabel(nucBadThresh_crop)
      color_img_strands<- colorLabels(strands, normalize = TRUE)
      num_strands <- computeFeatures.shape(strands)
      num_strands <- data.frame(num_strands)
      ### define a foci mask
      w = makeBrush(size = 51, shape = 'gaussian', sigma = 10)
      img_flo_foci = filter2(img_orig_foci, w)
      ## do adaptive thresholding on the image.. this should solve the issue...
      disc = makeBrush(21, "disc")
      disc = disc / sum(disc)
      localBackground = filter2(img_orig_foci, disc)
      offset = 0.1
      nucBadThresh_foci = (img_orig_foci - localBackground > offset)
      foci_mask <- bwlabel(nucBadThresh_foci)
      #display(foci_mask)
      
      foci_mask_crop <- foci_mask
      
      
      bg <- mean(img_orig_foci)
      if(bg < 0.05 && bg > 0.02){
        display(foci_mask)
        offset = 4*bg
        foci_th = foci_mask_crop > bg + offset
        ### smooth it
        ### maybe up the contrast first??
        img_tmp_contrast = foci_mask_crop
        print("cell counter is")
        print(cell_count)
        display(foci_mask_crop)
        w = makeBrush(size = 3, shape = 'gaussian', sigma = 1)
        img_flo = filter2(img_tmp_contrast, w)
        ## only choose objects above bright pixel value
        #foci_th = foci_mask_crop > 0.2
        #foci_th = img_flo > 0.2
        
        #display(foci_th)
        foci_label = bwlabel(foci_th)
        foci_label <- channel(foci_label, "grey")
        display(strands)
        display(foci_label)
        
        ### multiply strands by foci_label
        display(rgbImage(strands,foci_label,foci_label))
        coincident_foci <- bwlabel(foci_label*strands)
        display(colorLabels(coincident_foci))
        overlap_no = table(coincident_foci)
        foci_per_cell <-  length(overlap_no)
        #print(file)
        #print("foci count was:")
        #print(foci_per_cell)
        if (foci_per_cell > 0){
          foci_counts <- append(foci_counts,foci_per_cell)
        }
        
        
        
      }
      
      
      ###
    }
  
  }
  hist(foci_counts)
  print(mean(foci_counts))
  print(median(foci_counts))
  print(sd(foci_counts))
  return(foci_counts)

}