## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----installation, include = TRUE, eval=FALSE---------------------------------
# if (!requireNamespace("BiocManager")) {
#     install.packages("BiocManager")
# }
# BiocManager::install("sosta")

## ----setup, message=FALSE, warning=FALSE--------------------------------------
library("sosta")
library("dplyr")
library("tidyr")
library("ggplot2")
library("sf")
library("SpatialExperiment")

theme_set(theme_bw())

## ----loading, echo=TRUE, message=FALSE----------------------------------------
# load the data
data("sostaSPE")
sostaSPE

## -----------------------------------------------------------------------------
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = cellType)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal()

## -----------------------------------------------------------------------------
shapeIntensityImage(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    imageId = "image1",
    markSelect = "A"
)

## -----------------------------------------------------------------------------
n <- estimateReconstructionParametersSPE(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    markSelect = "A",
    plotHist = TRUE
)

## -----------------------------------------------------------------------------
(thresSPE <- mean(n$thres))
(bndwSPE <- mean(n$bndw))

## -----------------------------------------------------------------------------
(struct <- reconstructShapeDensityImage(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    imageId = "image1",
    markSelect = "A",
    bndw = bndwSPE,
    dim = 500,
    thres = thresSPE
))

## -----------------------------------------------------------------------------
cbind(
    colData(sostaSPE[, sostaSPE$imageName == "image1"]),
    spatialCoords(sostaSPE[, sostaSPE$imageName == "image1"])
) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = cellType)) +
    geom_point(size = 0.5) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = struct,
        fill = NA,
        color = "darkblue",
        inherit.aes = FALSE, # this is important
        linewidth = 0.75
    )

## -----------------------------------------------------------------------------
struct2 <- reconstructShapeDensityImage(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    imageId = "image1",
    markSelect = "A",
    dim = 500
)

## -----------------------------------------------------------------------------
cbind(
    colData(sostaSPE[, sostaSPE$imageName == "image1"]),
    spatialCoords(sostaSPE[, sostaSPE$imageName == "image1"])
) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = cellType)) +
    geom_point(size = 0.5) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = struct2,
        fill = NA,
        color = "darkblue",
        inherit.aes = FALSE, # this is important
        linewidth = 0.75
    )

## ----eval=TRUE----------------------------------------------------------------
allStructs <- reconstructShapeDensitySPE(
    sostaSPE,
    marks = "cellType",
    imageCol = "imageName",
    markSelect = "A",
    bndw = bndwSPE,
    thres = thresSPE,
    nCores = 1
)

## -----------------------------------------------------------------------------
# Define colnames by numbering the cells
colnames(sostaSPE) <- paste0("cell_", c(1:dim(sostaSPE)[2]))

## -----------------------------------------------------------------------------
assign <- assingCellsToStructures(sostaSPE, allStructs,
    imageCol = "imageName", nCores = 1
)
# Assign using the correct order of the columns in the spe object
sostaSPE$structAssign <- assign[colnames(sostaSPE)]

## -----------------------------------------------------------------------------
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = structAssign)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal()

## -----------------------------------------------------------------------------
cellTypeProportions(sostaSPE, "structAssign", "cellType")

## -----------------------------------------------------------------------------
shapeMetrics <- totalShapeMetrics(allStructs)
head(shapeMetrics)

## -----------------------------------------------------------------------------
cbind(allStructs, t(shapeMetrics)) |>
    ggplot() +
    geom_sf(aes(fill = Area)) +
    facet_wrap(~imageName)

## -----------------------------------------------------------------------------
sostaSPE$minDist <- minBoundaryDistances(
    spe = sostaSPE, imageCol = "imageName",
    structColumn = "structAssign", allStructs = allStructs
)

cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = minDist)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    scale_colour_gradient2() +
    geom_sf(
        data = allStructs,
        fill = NA,
        inherit.aes = FALSE
    ) +
    facet_wrap(~imageName)

## -----------------------------------------------------------------------------
sostaSPE$border <- ifelse(abs(sostaSPE$minDist) < 3, TRUE, FALSE)


cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = border)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = allStructs,
        fill = NA,
        inherit.aes = FALSE
    ) +
    facet_wrap(~imageName)

## -----------------------------------------------------------------------------
borders <- lapply(
    st_geometry(allStructs),
    \(x) st_difference(st_buffer(x, 3), st_buffer(x, -3))
) |>
    st_as_sfc() |>
    st_as_sf() # both functions needed for proper conversion

borders$imageName <- allStructs$imageName
borders$borderID <- paste0("border_", allStructs$structID)

borderAssign <- assingCellsToStructures(sostaSPE,
    borders,
    imageCol = "imageName",
    uniqueId = "borderID",
    nCores = 1
)
# Assign using the correct order of the columns in the spe object
sostaSPE$borderSf <- borderAssign[colnames(sostaSPE)]

## -----------------------------------------------------------------------------
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = borderSf)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    geom_sf(
        data = borders,
        fill = NA,
        inherit.aes = FALSE
    ) +
    facet_wrap(~imageName)

## -----------------------------------------------------------------------------
# create the fov bounding box using st_bbox
fov_bbox <- st_bbox(c(xmin = 0, ymin = 0, xmax = 128, ymax = 128))

# convert to polygon
fov <- st_as_sfc(fov_bbox)

# one fov per image
fovBorder <- rep(fov, length(unique(sostaSPE$imageName))) |>
  st_as_sf()

# name the polygons accordingly
fovBorder$imageName <- unique(sostaSPE$imageName)

fovBorder

## -----------------------------------------------------------------------------
sostaSPE$fovBorderDist <- minBoundaryDistances(
  spe = sostaSPE,
  imageCol = "imageName",
  structColumn = "imageName", # here any variable that is true for all cells
  allStructs = fovBorder
) |> 
  abs()

## -----------------------------------------------------------------------------
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = fovBorderDist)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    facet_wrap(~imageName)

## -----------------------------------------------------------------------------
sostaSPE$borderCorrected <- 
  ifelse(sostaSPE$fovBorderDist > 5, sostaSPE$borderSf, NA)

## -----------------------------------------------------------------------------
cbind(colData(sostaSPE), spatialCoords(sostaSPE)) |>
    as.data.frame() |>
    ggplot(aes(x = x, y = y, color = borderCorrected)) +
    geom_point(size = 0.25) +
    facet_wrap(~imageName) +
    coord_equal() +
    facet_wrap(~imageName)

## ----sessionInfo--------------------------------------------------------------
sessionInfo()