# Load R packages

library(raster)

library(rgdal)

library(Rcpp)


# Translate compiled code from Rcpp file into the R working environment  

sourceCpp("C:/Users/35799/Desktop/Manuscript Tracy/For publication/RootHairFinder.cpp")




# Import rhizotron image tiff file as raster stack

## This is the utility of the "raster" package ##

Image <- stack("C:/Users/35799/Desktop/Rhizo/mine_run2/Raw/R15_new2.tif")

# Greyscale conversion (if not already)

if ((dim(Image)[3]) == 3) {

# Express red spectrum in the form of a matrix

red_matrix <- as.matrix(Image[[1]])

# Express green spectrum in the form of a matrix

green_matrix <- as.matrix(Image[[2]])

# Express blue spectrum in the form of a matrix

blue_matrix <- as.matrix(Image[[3]])

# Convert to greyscale

grey_matrix <- greyscale_par(red_matrix, green_matrix, blue_matrix)

rm("Image", "red_matrix", "green_matrix", "blue_matrix")

} else {

# Express greyscale in the form of a matrix

grey_matrix <- as.matrix(Image[[1]])

rm("Image") }

# Reduce matrix dimensions to even numbers

if ((dim(grey_matrix)[1] %% 2) != 0) { grey_matrix <- grey_matrix[-1, ] }
if ((dim(grey_matrix)[2] %% 2) != 0) { grey_matrix <- grey_matrix[, -1] }

# Identify foreground pixel id

foreground_id <- filter_indices(grey_matrix)

# Apply Frangi filter for root detection

vessels_root <- frangi_filter_par(grey_matrix, foreground_id, sigma = 10)

rm("foreground_id")

# Root binarization

root_binary <- root_binarization(vessels_root)

rm("vessels_root")

# Apply Frangi filter for root hair detection

vessels_root_hair <- root_hair_enhance(grey_matrix, root_binary, roothair_pixel_radius = 100)

# Compute output root image and parameters 

results <- root_outline(x = root_binary, y = grey_matrix, vessels_matrix = vessels_root_hair, dpi = 1200, roothair_pixel_radius = 100)

rm("root_binary", "grey_matrix", "vessels_root_hair")

# Print parameters table

data.frame("Parameter" = c("x coord min", "x coord max", "y coord min", "y coord max", "root length x", "root length y", "root length total", "root area", "root hair area", "root hair to root area ratio"), "Value" = c(unlist(results[2:11], use.names = FALSE)))

# View output root image

r <- raster(ncol = nrow(results[[1]]), nrow = ncol(results[[1]]), xmx = ncol(results[[1]]), xmn = 1, ymn = 1, ymx = nrow(results[[1]]))
values(r) = as.vector(results[[1]])
plot(r)

# Export output root image (tiff)

x <- writeRaster(r, 'output.tif', overwrite = TRUE)




Image <- stack("C:/Users/35799/Desktop/Rhizo/mine_run2/Binary/R15_Binary_new2.tif")

my_binary <- as.matrix(Image[[1]])

if ((dim(my_binary)[1] %% 2) != 0) { my_binary <- my_binary[-1, ] }
if ((dim(my_binary)[2] %% 2) != 0) { my_binary <- my_binary[, -1] }

eval_binary(x = results[[1]], binary_matrix = my_binary)