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

## ----internet-check, include=FALSE--------------------------------------------
knitr::opts_chunk$set(error = FALSE)
if (!requireNamespace("curl", quietly = TRUE) || !curl::has_internet()) {
  message("No internet connection detected — skipping sections requiring downloads.")
  knitr::knit_exit()
}

## ----setup--------------------------------------------------------------------
library(BoneDensityMapping)
library(rgl)
library(testthat)

## -----------------------------------------------------------------------------
# CT scan file
url_scan <- "https://github.com/Telfer/BoneDensityMapping/releases/download/v1.0.1/test_CT_hip.nii.gz"
scan_file <- tempfile(fileext = ".nii.gz")
scan_success <- tryCatch({
  download.file(url_scan, scan_file, mode = "wb", quiet = TRUE)
  TRUE
}, error = function(e) {
  message("Failed to download scan: ", e$message)
  FALSE
})

skip_if_not(scan_success, "Scan download failed")

CTscan <- tryCatch({
  import_scan(scan_file)
}, error = function(e) {
  skip(paste("import_scan failed:", e$message))
})

# STL mesh file
url_mesh <- "https://github.com/Telfer/BoneDensityMapping/releases/download/v1.0.2/test_CT_femur.stl"
mesh_file <- tempfile(fileext = ".stl")
mesh_success <- tryCatch({
  download.file(url_mesh, mesh_file, mode = "wb", quiet = TRUE)
  TRUE
}, error = function(e) {
  message("Failed to download mesh: ", e$message)
  FALSE
})

skip_if_not(mesh_success, "Mesh download failed")

surface_mesh <- tryCatch({
  import_mesh(mesh_file)
}, error = function(e) {
  skip(paste("import_mesh failed:", e$message))
})

landmark_path <- system.file("extdata", "test_femur.mrk.json",
                             package = "BoneDensityMapping")
landmarks <- import_lmks(landmark_path)

## -----------------------------------------------------------------------------
landmark_check(surface_mesh, landmarks, threshold = 1.0)

## -----------------------------------------------------------------------------
bone_scan_check(surface_mesh, CTscan)

## -----------------------------------------------------------------------------
surface_density <- surface_normal_intersect(surface_mesh, normal_dist = 3.0, 
                                            nifti = CTscan, ct_params = c(1, 0.74))

## -----------------------------------------------------------------------------
surface_colors <- color_mapping(surface_density)

## -----------------------------------------------------------------------------
plot_mesh(surface_mesh, surface_colors, title = 'Bone surface', 
          legend_maxi = max(surface_density), legend_mini = min(surface_density))

## -----------------------------------------------------------------------------
internal_fill <- fill_bone_points(surface_mesh, 1)

## -----------------------------------------------------------------------------
internal_density <- voxel_point_intersect(internal_fill, CTscan)

## -----------------------------------------------------------------------------
internal_colors <- color_mapping(internal_density)

## -----------------------------------------------------------------------------
plot_cross_section_bone(surface_mesh, surface_colors, IncludeSurface = FALSE, 
                        internal_fill, internal_colors, slice_axis = 'x', 
                        slice_val = 0.5, legend_maxi = max(internal_density), 
                        legend_mini = min(internal_density))

## -----------------------------------------------------------------------------
patients <- list()

base_url <- "https://github.com/Telfer/BoneDensityMapping/releases/download/v1.0.2/"
patients <- list()

for (i in 1:6) {
  id <- paste0("SCAP00", i)
  
  # STL mesh file
  url_mesh <- paste0(base_url, id, ".stl")
  mesh_file <- tempfile(fileext = ".stl")
  mesh_success <- tryCatch({
    download.file(url_mesh, mesh_file, mode = "wb", quiet = TRUE)
    TRUE
  }, error = function(e) {
    message("Failed to download mesh: ", e$message)
    FALSE
  })
  
  skip_if_not(mesh_success, "Mesh download failed")
  
  surface_mesh <- tryCatch({
    import_mesh(mesh_file)
  }, error = function(e) {
    skip(paste("import_mesh failed:", e$message))
  })
  
  # CT scan file
  url_scan <- paste0(base_url, id, "_resampled.nii.gz")
  scan_file <- tempfile(fileext = ".nii.gz")
  scan_success <- tryCatch({
    download.file(url_scan, scan_file, mode = "wb", quiet = TRUE)
    TRUE
  }, error = function(e) {
    message("Failed to download scan: ", e$message)
    FALSE
  })
  
  skip_if_not(scan_success, "Scan download failed")
  
  nifti <- tryCatch({
    import_scan(scan_file)
  }, error = function(e) {
    skip(paste("import_scan failed:", e$message))
  })
  
  landmark_path <- system.file("extdata", paste0(id, "_landmarks.fcsv"), package = "BoneDensityMapping")
  
  patients[[id]] <- list(
    mesh = surface_mesh,
    landmarks = import_lmks(landmark_path),
    ct = nifti
  )
}

## -----------------------------------------------------------------------------
patients$SCAP001$surface_points <- surface_points_template(patients$SCAP001$mesh, patients$SCAP001$landmarks, 5000)  

## -----------------------------------------------------------------------------
patient_sides <- data.frame(
  ID = sprintf("SCAP%03d", 1:6),
  Side = c("LEFT", "RIGHT", "LEFT", "LEFT", "RIGHT", "RIGHT"),
  stringsAsFactors = FALSE
)

for (i in 2:6) {
  id <- sprintf("SCAP%03d", i)

  mirror_side <- if (patient_sides$Side[i] == "RIGHT") "x" else FALSE

  patients[[id]]$surface_points <- surface_points_new(
    patients[[id]]$mesh,
    patients[[id]]$landmarks,
    patients$SCAP001$surface_points,
    mirror = mirror_side,
    plot_check = FALSE
  ) 
}

## -----------------------------------------------------------------------------
for (i in 1:6) {
  id <- sprintf("SCAP%03d", i)

  patients[[id]]$surface_density <- surface_normal_intersect(
    patients[[id]]$mesh,
    patients[[id]]$surface_points,
    normal_dist = 1.0,
    patients[[id]]$ct, ct_params = c(1, 0.74)
  )
}

## -----------------------------------------------------------------------------
for (i in 1:6) {
  id <- sprintf("SCAP%03d", i)
  
  patients[[id]]$colored_mesh <- color_mesh(
    patients$SCAP001$mesh,
    patients$SCAP001$surface_points,
    patients[[id]]$surface_density,
    maxi = 2000, 
    mini = 0
  )
}

plot_mesh(patients$SCAP001$colored_mesh, title = 'SCAP001')
plot_mesh(patients$SCAP002$colored_mesh, title = 'SCAP002')
plot_mesh(patients$SCAP003$colored_mesh, title = 'SCAP003')
plot_mesh(patients$SCAP004$colored_mesh, title = 'SCAP004')
plot_mesh(patients$SCAP005$colored_mesh, title = 'SCAP005')
plot_mesh(patients$SCAP006$colored_mesh, title = 'SCAP006')


## -----------------------------------------------------------------------------
# Collect all surface density vectors into a matrix
density_matrix <- do.call(cbind, lapply(1:6, function(i) {
  id <- sprintf("SCAP%03d", i)
  patients[[id]]$surface_density
}))

# Average across columns (i.e. across patients), row by row
average_density <- rowMeans(density_matrix)

average_color_mesh <- color_mesh(
    patients$SCAP001$mesh,
    patients$SCAP001$surface_points,
    average_density
  )

plot_mesh(average_color_mesh, title = 'Average density')

## -----------------------------------------------------------------------------
# Define patient IDs
young_ids <- c(1, 4, 6)
old_ids <- c(2, 3, 5)

# Helper function to get surface densities
get_densities <- function(indices) {
  do.call(cbind, lapply(indices, function(i) {
    id <- sprintf("SCAP%03d", i)
    patients[[id]]$surface_density
  }))
}

average_young_density <- rowMeans(get_densities(young_ids))
average_old_density <- rowMeans(get_densities(old_ids))

# Create colored meshes for each group
young_colored_mesh <- color_mesh(patients$SCAP001$mesh, patients$SCAP001$surface_points, average_young_density)
old_colored_mesh <- color_mesh(patients$SCAP001$mesh, patients$SCAP001$surface_points, average_old_density)

close3d()
plot_mesh(young_colored_mesh, title = "Younger Group")
plot_mesh(old_colored_mesh, title = "Older Group")


density_diff <- average_young_density - average_old_density
max_abs <- max(abs(density_diff))

# Color scheme: red = young > old, blue = old > young
my_colors <- c("blue", "white", "red")

# Create colored mesh of differences
diff_colored_mesh <- color_mesh(
  patients$SCAP001$mesh,
  patients$SCAP001$surface_points,
  density_diff,
  maxi = max_abs,
  mini = -max_abs,
  color_sel = my_colors
)

# Plot it
plot_mesh(diff_colored_mesh, title = "Younger - Older Density Difference", legend_maxi = max_abs, legend_mini = -max_abs, legend_color_sel = my_colors)