message("Just wait for a bit...")
Sys.sleep(3)
message("Done.")

# Load the libraries
library(romero.gateway)
library(EBImage, warn.conflicts = FALSE)

user_name = readline('Username: ')
user_password <- getPass::getPass('OMERO password: ')

server <- OMEROServer(host = 'wss://workshop.openmicroscopy.org/omero-ws', port = 443L, username=user_name, password=user_password)
server <- connect(server)
paste('Successfully logged in as', server@user$getUserName())

imageId <- REPLACE_WITH_IMAGE_ID 
image <- loadObject(server, "ImageData", imageId)
paste("Image", imageId, "loaded.")

# There is just one plane, so z = 1 and t = 1
z <- 1
t <- 1

# Load the second channel
channelIndex <- 2

pixels <- getPixelValues(image, z, t, channelIndex)

ebimage <- EBImage::Image(data = pixels, colormode = 'Grayscale')
ebimage <- normalize(ebimage)
EBImage::display(ebimage)

# Threshold
threshImg <- thresh(ebimage, w=10, h=10, offset=0.01)

# Remove noise
threshImg <- medianFilter(threshImg, size=3)

# Fill holes
threshImg <- fillHull(threshImg)          

# Distinguish objects
threshImg <- bwlabel(threshImg)

# Show the segmented image
EBImage::display(colorLabels(threshImg))

shapes = computeFeatures.shape(threshImg)
moments = computeFeatures.moment(threshImg)

# merge the two dataframes together into one 'features' dataframe
features <- merge(shapes, moments, by=0, all=TRUE)
features

#' Create ROIs from a features table
#'
#' @param features The shape and moments generated by EBImage computeFeatures
#' @return A dataframe specifying the x, y, rx, ry, w, h, theta and text parameters of the ROIs
createROIs <- function(features) {
    rois <- data.frame(x = c(0), y = c(0), rx = c(0), ry = c(0), w = c(0), h = c(0), theta = c(0), text = c('remove'), stringsAsFactors = FALSE)
    for (index in 1:length(features[,1])) {
        x <- features[index,8]
        y <- features[index,9]
        r1 <- features[index,10]
        ecc <- features[index,11]
        r2 <- sqrt(r1^2 * (1 - ecc^2))
        theta <- features[index,12]
        rois <- rbind(rois, c(x, y, r1, r2, NA, NA, -theta, as.character(index)))
    }
    rois <- rois[-1,]
    rownames(rois) <- c()
    return(rois)
}
print("Function 'createROIs' defined.")

rois <- createROIs(features)
rois

# addROIs creates an ROI for each entry in the dataframe specified by the 'coords' parameter
addROIs(image, coords = rois)
print("ROIs successfully created.")

# 'attachDataframe' directly attaches an R dataframe to an OMERO image (project, dataset, etc.)
# ('invisible' not strictly necessary, just suppresses some unnecessary output)
invisible(attachDataframe(image, features, "ROI features"))
print("Data frame successfully attached.")

csvFile <- "/tmp/ROI_features.csv"
write.csv(features, file = csvFile)
fileannotation <- attachFile(image, csvFile)
print("CSV file successfully attached.")

#' Performs an image segmentation to find the largest ROI of the image
#' and returns some features of interest (area, perimeter and diameter).
#' Optionally: Creates an ROI for each object detected by the segmentation.
#'
#' @param image The image to segment
#' @param channelName The channel to be taken into account
#' @param df The dataframe to add the features to (channelName, imageName, ImageID, ROIIndex, area, perimeter, diameter)
#' @param saveROIs Flag if ROIs should be created and attached to the image (default: FALSE)
#' @return The dataframe with the features
analyzeImage <- function(image, channelName, df, saveROIs = FALSE) {
    # Find the channel index
    chnames <- getChannelNames(image)
    chindex <- match(channelName, chnames, nomatch = 0)
    if (chindex == 0) {
      message (paste("Could not resolve channel name, skipping ", image@dataobject$getId()))
      return(df)
    }
    
    # Load the pixels
    pixels <- getPixelValues(image, 1, 1, chindex)
    ebi <- EBImage::Image(data = pixels, colormode = 'Grayscale')
    ebi <- normalize(ebi)
     
    # this is our segmentation workflow from above
    threshImg <- thresh(ebi, w=15, h=15, offset=0.1)
    threshImg <- medianFilter(threshImg, size=3)
    threshImg <- fillHull(threshImg)          
    threshImg <- bwlabel(threshImg)
    
    # Calculate the features
    shapes = suppressMessages(computeFeatures.shape(threshImg))
    moments = suppressMessages(computeFeatures.moment(threshImg))
    features <- merge(shapes, moments, by=0, all=TRUE)
    
    if (length(features[,1])>1) {
        # Add the ROIs to the image
        if (saveROIs) {
            rois <- createROIs(features)
            addROIs(image, coords = rois)
        }
        
        # Add the interesting properties (area, perimeter and diameter)
        # of the largest object together with channel name, image name, image id 
        # and roi index to the dataframe
        features <- features[order(-features[,2]),]
        diameter <- features[1,4]*2
        df <- rbind(df, c(channelName, image@dataobject$getName(), image@dataobject$getId(), features[1,1], features[1,2], features[1,3], diameter))
    }
    return(df)
}
print("Function 'analyzeImage' defined.")

datasetId <- REPLACE_WITH_DATASET_ID

channelName <- 'CDK5RAP2-C'

dataset <- loadObject(server, "DatasetData", datasetId)

# Keep the channel name, image name, image id, area, perimeter, and diameter of the largest ROIs
result <- data.frame(Channel = c('remove'), ImageName = c('remove'), Image = c(0), ROIIndex = c(0), Area = c(0), Perimeter = c(0), Diameter = c(0), stringsAsFactors = FALSE)

images <- getImages(dataset)
for (image in images) {
    result <- tryCatch({
        analyzeImage(image, channelName, result, saveROIs = TRUE)
    }, warning = function(war) {
        message(paste("WARNING:  ", image@dataobject$getId(),war))
        return(result)
    }, error = function(err) {
        message(paste("ERROR:  ", image@dataobject$getId() ,err))
        return(result)
    }, finally = {
    })
}

result <- result[-1,]
rownames(result) <- c()

# set the correct datatypes
result$Channel <- as.factor(result$Channel)
result$Area <- as.numeric(result$Area)
result$Perimeter <- as.numeric(result$Perimeter)
result$Diameter <- as.numeric(result$Diameter)

result

pxSizeInNM <- NA # REPLACE WITH PIXEL SIZE
if (!is.na(pxSizeInNM)) {
    # Replace result$Area, Perimeter, Diameter with the
    # size in nanometers, using the variable 'pxSizeInNM'
    # declared above.
    
    # Enter code here...
    
    # Just print out the data frame again to check result:
    result
}

projectId <- REPLACE_WITH_PROJECT_ID 
project <- loadObject(server, "ProjectData", projectId)
dataframes <- availableDataframes(project)
print(dataframes)

# Load the dataframe
dfID <- dataframes$ID[1]
centrioles <- loadDataframe(project, dfID)

# Make sure the data types are correct
centrioles$Dataset <- as.factor(centrioles$Dataset)
centrioles$Diameter <- as.numeric(centrioles$Diameter)
centrioles

# Order the datasets ascending by mean diameter 
ag <-aggregate(centrioles$Diameter ~ centrioles$Dataset, centrioles, median)
ordered <- factor(centrioles$Dataset, levels=ag[order(ag$`centrioles$Diameter`), 'centrioles$Dataset'])

# Draw the plot 
plot(centrioles$Diameter ~ ordered, ylab='Diameter (nm)', xlab="Protein", cex.axis=0.5)

fit <- aov(centrioles$Diameter ~ centrioles$Dataset)
summary(fit)

# Two-sample Wilcoxon test ('Mann-Whitney') of all pairwise combinations:
combins <- combn(levels(centrioles$Dataset), 2)
params_list <- split(as.vector(combins), rep(1:ncol(combins), each = nrow(combins)))
testResults <- data.frame()
for (param in params_list) {
  testdf <- subset(centrioles, centrioles$Dataset %in% param)
  pval <- wilcox.test(formula = Diameter ~ Dataset, data = testdf)$p.value
  testResults <- rbind(testResults, data.frame(Protein_1=param[1], Protein_2=param[2], p_value=pval))
}
testResults

disconnect(server)

# Solution for the optional exercise ('set pixel size in nanometer'):

pxSizeInNM <- 40
result$Area <- result$Area * pxSizeInNM ^ 2
result$Perimeter <- result$Perimeter * pxSizeInNM
result$Diameter <- result$Diameter * pxSizeInNM