### FODs para Histología

En este tutorial veremos como calcular los Fiber Orientation Distributions (FODs) en imágenes histologicas. Este tipo de análisis es común para las imágenes de resonancia magnética sensibles a difusión, sin embargo, extrapolarlo de alguna manera a las imágenes histológicas puede ser de gran ayuda para la validación directa.

Este acercamiento ha sido anteriormente descrito en el articulo:
> Budde MD, Annese J. Quantification of anisotropy and fiber orientation in human brain histological sections. Front Integr Neurosci. 2013 Feb 1;7:3. doi: 10.3389/fnint.2013.00003. PMID: 23378830; PMCID: PMC3561729.

Sin embargo, el código y la manera exacta en como calculan los FODs no es completamente clara. Aqui nos dimos a la tarea de crear unos códigos usando una combinacion entre `javascript` y `R`, lo cual veremos a continuación.

Para este ejercicio utilizamos la siguiente tincion de plata de la corteza, la cual podemos ver que tiene las fibras que ascienden y descienden, asi como aquellas que van tangencialmente en el plano:


Ahora bien, el primer paso es abrir `FIJI` y la imágen a analizar. Despues corremos `Pluing > Macros > New` con el código que se proporciona aqui abajo:
> Puedes modificar la dimension del parche al cambiar el número de pixeles de altura y ancho.

// Get image dimensions
width = getWidth();
height = getHeight();

// Suggested number of patches
suggestedPatchesX = Math.ceil(width / 400);
suggestedPatchesY = Math.ceil(height / 400);

// Prompt the user to input the number of patches or use the suggested value
Dialog.create("Especifica la cantidad de parches por imagen");
Dialog.addNumber("Numero de parches en eje X (Sugerido: " + suggestedPatchesX+ ")", suggestedPatchesX);
Dialog.addNumber("Numero de parches en eje Y (Sugerido: " + suggestedPatchesY+ ")", suggestedPatchesY);
Dialog.show();

numPatchesX = Dialog.getNumber();
numPatchesY = Dialog.getNumber();

// Calculate patch size based on the number of patches
patchSizeX = Math.floor(width / numPatchesX);
patchSizeY = Math.floor(height / numPatchesY);

// Create folders for storing patches and results
patchFolder = getDirectory("current") + "patches/";
resultFolder = getDirectory("current") + "results/";
File.makeDirectory(patchFolder);
File.makeDirectory(resultFolder);

// Loop through the image in patches
for (x = 0; x < width; x += patchSizeX) {
    for (y = 0; y < height; y += patchSizeY) {
        // Define the size of the current patch
        currentPatchWidth = patchSizeX;
        currentPatchHeight = patchSizeY;

        // Ensure the patch fits within the image boundaries
        if (x + currentPatchWidth > width) {
            currentPatchWidth = width - x;
        }
        if (y + currentPatchHeight > height) {
            currentPatchHeight = height - y;
        }

        // Define ROIs
        makeRectangle(x, y, currentPatchWidth, currentPatchHeight);
        
        // Duplicate the ROIs and save them
        title = "ROI_" + x + "_" + y;
        run("Duplicate...", "title=" + title);
        
        // Save the patch
        patchFilePath = patchFolder + title + ".tif";
        saveAs("Tiff", patchFilePath);
        
        // Close the duplicated ROI window
        close();
    }
}

// Apply analysis to each patch and save the results
fileList = getFileList(patchFolder);
for (i = 0; i < fileList.length; i++) {
    // Open the patch image
    open(patchFolder + fileList[i]);

    // Perform analysis (OrientationJ Distribution on the duplicated patch)
    run("OrientationJ Distribution", "tensor=5");

    // Name the new table file (output of 0J) according to the patch file name
    fileName = fileList[i].replace(".tif", ".csv");

    // Save the results
    filePath = resultFolder + fileName;
    saveAs("Results", filePath);

    // Close the patch image
    close();
}

El codigo va a crear dos carpetas, una llamada **patches** donde va a almacenar cada recuadro al que se le hizo `OrientationJ Distribution` , y la segunda carpeta es la de **results** que va a contener los archivos `.csv` con la información de la orientación y anisotropía de cada uno de los parches. **MUY importarte es que tanto los parches como los archivos CSV se guardan de acuerdo a las coordenadas que les corresponde en la imágen. Esto nos será realmente util para poder realizar la recontrucción y saber quién es quién**.


Ahora bien, ya que tenemos los resultados del tensor de estructura por parche, podemos proceder a la segunda parte que es calcular los FODs. Para esto utilizaremos el siguiente código en `R`.

Lo unico que necesitamos para este paso es la carpeta en donde se encuentran los archivos `.csv`, que contienen las distribuciones que queremos graficar.

## Load the libraries
library(tidyverse)
library(dplyr)
library(readr)
library(stringr)
library(ggplot2)
library(xml2)

## Check directory
getwd()
## Change directory where my csv are located
setwd("/Users/paulinav/Desktop/results/")
## list all files from the current directory
list.files()
## Crate a list from these files 
list_files <- list.files()
list_files
## Create an empty list that will serve as a container to the receiving files
files_all <- list()
files_all
## Loop to read all files
for (i in 1:length(list_files)) {
  files_all[[i]] <- read_csv(list_files[i])
}
files_all
## Add the name of the files 
names(files_all) <- list_files
files_all
## Merged the files into a dataframe by rows
df <- do.call(rbind, files_all)
## Convert rownames to a first column (name of the file)
dataframe <- cbind(rownames(df), data.frame(df, row.names = NULL)) 
colnames(dataframe)[1] <- c("filename")
## Normalize the negatives to positive 
dataframe$adjOrientation <- ifelse(
  dataframe$Orientation < 0, 180 + dataframe$Orientation, 
  dataframe$Orientation) 

## Do the mirror of FODs

dataframe$reversed2 <- ifelse(
  dataframe$adjOrientation > 0, 180 + dataframe$adjOrientation, 
  dataframe$adjOrientation)

## Create new variables for coordinates
# Use str_split_fixed() to split the filenames and extract X and Y
dataframe$split_names <- str_split_fixed(dataframe$filename, "_", 3)
dataframe$X <- as.numeric(dataframe$split_names[, 2])
dataframe$Y <- as.numeric(str_replace(dataframe$split_names[, 3], ".csv.*", ""))

dataframe <- dataframe[, -5]
## Create a variable "group" so we can differentiate between patches prior to plotting
paste(sample(LETTERS[cur_group_id()], 10, replace=TRUE), collapse="")

dataframe <- dataframe %>% group_by(X, Y) %>%
 mutate(group = paste(sample(LETTERS, 3, replace=TRUE), collapse="")) 
 # mutate(group = rep.len(LETTERS[cur_group_id()]))


"Lets now plot in bulk for each combination of X and Y"

## We will use ggplot and geom_line + coord_polar
dataframe

group_unique <- unique(dataframe$group)
X_unique <- unique(dataframe$X)
Y_unique <- unique(dataframe$Y)

plot_list <- list()


for (i in seq_along(group_unique)) {
  grp <- group_unique[i]
  
  # Subset the dataframe for the current group
  df_sub <- subset(dataframe, group == grp)
  
  # Create the plot
  p <- ggplot(data = df_sub) +
    geom_line(aes(x=adjOrientation, y=Slice1), colour = "white", size = 2) +
    geom_line(aes(x=reversed2, y=Slice1), colour = "white", size = 2) +
    coord_polar()  +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          legend.text = element_blank(),
          legend.title = element_blank(),
          legend.position = "none",
          axis.line = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank(),
          rect = element_rect(fill = "transparent"),
          plot.background = element_rect(fill = "transparent",
                                         color = NA)
    )  
    #ggtitle(paste0("X: ", df_sub$X[1], ", Y: ", df_sub$Y[1])) 
  
  # Add the plot to the plot list
  plot_list[[paste0("Group_", grp)]] <- p
  
  # Create a filename based on X and Y
  filename <- paste0(df_sub$X[1], "_", df_sub$Y[1], ".png")
  
  # Save the plot as an SVG file
  ggsave(filename = filename, 
         plot = p, 
         device = "png", 
        # width = 8, 
        # height = 6, 
         bg = "transparent")
  
}

Como resultado final, tendremos una carpeta en donde se almacenarán todos los gráficos que representan los FODs de cada uno de los parches. Y como podrás notar, los FODs se guardaron también de acuerdo a las coordenadas de la imágen de histologia original. Lo cual nos ayudará a la siguiente fase que es la recontrucción de los parches y el overlay de los FODs.