1. Wstęp

Projekt stanowi badawczy, nienadzorowany system analizy obrazów oka zewnętrznego, ukierunkowany na wykrywanie odchyleń morfologicznych i teksturalnych widocznych w obrębie spojówki, twardówki oraz powierzchni rogówki. Celem nie jest diagnostyka nozologiczna, lecz identyfikacja nieprawidłowości fenotypu obrazowego (np. przekrwienie, zaburzenia mikrokrążenia, nieregularność powierzchni), które mogą stanowić przesłankę do dalszej oceny klinicznej.

Metodyka łączy klasyczne przetwarzanie obrazu z podejściami inspirowanymi CNN: lokalne pola recepcyjne, filtracja przestrzenna oraz analiza gradientowa pozwalają wydzielić mapy cech o potencjalnie wysokiej wartości diagnostycznej. Wykorzystanie tekstur GLCM i miar statystycznych umożliwia opis zarówno składowej barwnej (zaczerwienienie), jak i mikrostruktury tkanek.

Zastosowanie standaryzacji Z‑score oraz skoringu anomalii tworzy spójne ramy wczesnego wykrywania zmian w danych obrazowych, co jest zgodne z paradygmatem medycznych systemów wspomagania decyzji i monitorowania stanu powierzchni oka.

Zbiór danych został przygotowany na podstawie obrazów pozyskanych z wyszukiwarki Google Grafika (obrazy ogólnodostępne), wyselekcjonowanych i ręcznie zweryfikowanych pod kątem jakości oraz zgodności z tematyką badania.

2. Podstawy teoretyczne

2.1 Konwolucja i lokalne pola recepcyjne

Operacja splotu stanowi podstawę lokalnej ekstrakcji informacji w obrazach biomedycznych, umożliwiając uwydatnienie krawędzi, gradientów oraz drobnych zmian strukturalnych:

\[ S(i,j) = \sum_{m=-k}^{k} \sum_{n=-k}^{k} I(i+m, j+n) \cdot K(m,n) \]

W kontekście obrazów oka splot jest użyteczny w detekcji przejść tonalnych związanych z naczyniami, mikrouszkodzeniami nabłonka oraz lokalnymi nieregularnościami powierzchni.

2.2 Analiza częstotliwości przestrzennej

Zmiany patologiczne w tkankach oka często skutkują wzrostem komponentów wysokoczęstotliwościowych w obrazie, co odpowiada:

  • mikrostrukturze naczyń
  • nieregularności powierzchni
  • zwiększonej chropowatości teksturalnej
  • lokalnym zmianom kontrastu

Miary częstotliwości przestrzennej pełnią rolę wskaźników chropowatości i niejednorodności obrazu, co ma znaczenie w detekcji wczesnych zmian zapalnych lub dystroficznych.

2.3 Deskryptory statystyczne

System wykorzystuje zestaw deskryptorów statystycznych opisujących rozkład intensywności oraz złożoność informacji:

  • Energy – intensywność gradientów
  • Variance – rozproszenie intensywności
  • Entropy – złożoność informacyjna
  • Laplace Energy – uwydatnianie struktur drugiej pochodnej

W połączeniu z teksturami GLCM deskryptory te zapewniają stabilny opis zarówno globalnej, jak i lokalnej organizacji obrazu.

2.4 Nienadzorowana detekcja anomalii

Anomalia definiowana jest jako wektor cech istotnie odbiegający od rozkładu populacyjnego:

\[ Z = \frac{x - \mu}{\sigma} \]

oraz poprzez agregację wagową:

\[ AnomalyScore = \sum w_i Z_i \]

W praktyce klinicznej takie podejście pozwala na wykrywanie obrazów odstających bez potrzeby definiowania klas chorobowych, co jest szczególnie użyteczne w warunkach ograniczonej dostępności etykiet eksperckich.

3. Implementacja systemu

setwd("C:/Users/luck0/Documents/Kamil/UMN final project blok 2") # Ścieżka do folderu z danymi
# ================================
# Biblioteki
# ================================
required_pkgs <- c(
  "jpeg", "png", "webp", "imager", "entropy", "DT", "ggplot2",
  "logger", "purrr", "patchwork", "GLCMTextures", "dbscan", "Rtsne", "umap",
  "htmltools", "cluster"
)

#' Instaluje brakujące pakiety CRAN.
#' @param pkgs Wektor nazw pakietów.
#' @return Brak wartości zwrotnej.
ensure_packages <- function(pkgs) {
  missing_pkgs <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
  if (length(missing_pkgs) > 0) {
    install.packages(missing_pkgs)
  }
}

ensure_packages(required_pkgs)

library(jpeg)
library(png)
library(webp)
library(imager)
library(entropy)
library(DT)
library(ggplot2)
library(logger)
library(purrr)
library(patchwork)
library(GLCMTextures)
library(dbscan)
library(Rtsne)
library(umap)
library(htmltools)
library(cluster)

log_layout(layout_glue_generator("[{level}] {msg}"))

# ================================
# Funkcje pomocnicze
# ================================

# Zmienne sterujące obliczeniami GLCM
glcm_enabled <- TRUE
glcm_warned <- FALSE

#' Instaluje pakiet Bioconductor (jeśli jest wymagany).
#' @param pkg Nazwa pakietu Bioconductor.
#' @return Brak wartości zwrotnej.
ensure_bioc_package <- function(pkg) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    if (!requireNamespace("BiocManager", quietly = TRUE)) {
      install.packages("BiocManager")
    }
    BiocManager::install(pkg, ask = FALSE, update = FALSE)
  }
}

# Dostępność EBImage warunkuje użycie CLAHE
ensure_bioc_package("EBImage")

#' Bezpiecznie konwertuje obiekt cimg na macierz RGB.
#' @param img_cimg Obiekt klasy cimg.
#' @return Tablica 3D w układzie [x, y, kanał].
cimg_to_rgb_array <- function(img_cimg) {
  arr <- as.array(img_cimg)
  if (length(dim(arr)) == 4) {
    arr <- arr[,,1,]
  }
  arr
}

#' Wczytuje obraz oka z obsługą wielu formatów.
#' @param path Ścieżka do pliku obrazu.
#' @return Macierz RGB lub NULL przy błędzie.
load_eye_image <- function(path) {
  if (!file.exists(path)) {
    log_warn(sprintf("Plik nie istnieje: %s", path))
    return(NULL)
  }

  tryCatch({
    header <- readBin(path, what = "raw", n = 4)
    is_webp <- any(header == as.raw(c(0x52, 0x49, 0x46, 0x46)))

    img <- if (is_webp) {
      webp::read_webp(path)
    } else {
      ext <- tolower(tools::file_ext(path))
      switch(ext,
             "jpg"  = jpeg::readJPEG(path),
             "jpeg" = jpeg::readJPEG(path),
             "png"  = png::readPNG(path),
             webp::read_webp(path))
    }

    # Standaryzacja wymiarów do przestrzeni RGB
    if (length(dim(img)) == 2) {
      img <- array(img, dim = c(dim(img)[1], dim(img)[2], 1))
    }
    if (length(dim(img)) == 3 && dim(img)[3] == 4) {
      img <- img[,,1:3]
    }
    if (length(dim(img)) == 3 && dim(img)[3] == 1) {
      img <- array(img, dim = c(dim(img)[1], dim(img)[2], 3))
      img[,,2] <- img[,,1]
      img[,,3] <- img[,,1]
    }

    log_info(sprintf("Wczytano obraz: %s", basename(path)))
    img
  }, error = function(e) {
    log_error(sprintf("Błąd wczytywania obrazu %s: %s", path, e$message))
    NULL
  })
}

#' Normalizuje obraz w sposób odporny na wartości odstające.
#' @param img Obiekt cimg w skali 0-1.
#' @param q_low Dolny kwantyl.
#' @param q_high Górny kwantyl.
#' @return Obraz po normalizacji.
#' @details W razie braku rozrzutu stosuje normalizację min-max.
robust_normalize <- function(img, q_low = 0.01, q_high = 0.99) {
  v <- as.numeric(img)
  v <- v[is.finite(v)]
  qs <- stats::quantile(v, probs = c(q_low, q_high), na.rm = TRUE)
  if (diff(qs) < .Machine$double.eps) {
    return(normalize_minmax(img))
  }
  img_norm <- (img - qs[1]) / (qs[2] - qs[1])
  pmin(pmax(img_norm, 0), 1)
}

#' Normalizuje obraz metodą min-max.
#' @param img Obiekt cimg.
#' @return Obraz w skali 0-1.
normalize_minmax <- function(img) {
  v <- as.numeric(img)
  v <- v[is.finite(v)]
  if (length(v) == 0) return(img)
  rng <- range(v, na.rm = TRUE)
  if (diff(rng) < .Machine$double.eps) return(img)
  img_norm <- (img - rng[1]) / (rng[2] - rng[1])
  pmin(pmax(img_norm, 0), 1)
}

#' Skaluje obraz z zachowaniem proporcji i dopełnia do docelowego rozmiaru.
#' @param img Obiekt cimg.
#' @param target_size Wektor c(liczba_kolumn, liczba_wierszy).
#' @return Obraz o rozmiarze docelowym.
resize_with_aspect <- function(img, target_size = c(256, 256)) {
  target_w <- target_size[1]
  target_h <- target_size[2]
  w <- imager::width(img)
  h <- imager::height(img)
  if (!is.finite(w) || !is.finite(h) || w <= 0 || h <= 0) {
    log_warn("Nieprawidłowy rozmiar obrazu: stosuję bezpośrednie skalowanie")
    return(imager::resize(img, size_x = target_w, size_y = target_h))
  }
  scale <- min(target_w / w, target_h / h)
  new_w <- max(1, round(w * scale))
  new_h <- max(1, round(h * scale))

  resized <- tryCatch(
    imager::resize(img, size_x = new_w, size_y = new_h),
    error = function(e) {
      log_warn(sprintf("Błąd skalowania z zachowaniem proporcji: %s", e$message))
      imager::resize(img, size_x = target_w, size_y = target_h)
    }
  )
  pad_x <- target_w - new_w
  pad_y <- target_h - new_h

  if (pad_x > 0 || pad_y > 0) {
    pad_left <- floor(pad_x / 2)
    pad_right <- ceiling(pad_x / 2)
    pad_top <- floor(pad_y / 2)
    pad_bottom <- ceiling(pad_y / 2)

    arr <- as.array(resized)
    if (length(dim(arr)) == 4) {
      arr <- arr[,,1,]
    }
    if (length(dim(arr)) == 2) {
      arr <- array(arr, dim = c(dim(arr)[1], dim(arr)[2], 1))
    }

    w_arr <- dim(arr)[1]
    h_arr <- dim(arr)[2]
    ch <- dim(arr)[3]

    padded <- array(0, dim = c(
      w_arr + pad_left + pad_right,
      h_arr + pad_top + pad_bottom,
      ch
    ))
    padded[(pad_left + 1):(pad_left + w_arr),
           (pad_top + 1):(pad_top + h_arr), ] <- arr

    padded <- padded[1:target_w, 1:target_h, ]
    resized <- imager::as.cimg(padded)
  }

  resized
}

#' Wykonuje CLAHE na obrazie w skali szarości.
#' @param img_gray Obiekt cimg w skali szarości.
#' @return Obraz po CLAHE lub oryginał przy braku wsparcia.
apply_clahe <- function(img_gray) {
  if ("clahe" %in% getNamespaceExports("imager")) {
    return(imager::clahe(img_gray))
  }
  if (requireNamespace("EBImage", quietly = TRUE)) {
    img_mat <- as.matrix(img_gray)
    clahe_args <- list(x = img_mat, ntiles = 8, clip.limit = 0.01)
    valid_args <- intersect(names(clahe_args), names(formals(EBImage::clahe)))
    clahe_mat <- do.call(EBImage::clahe, clahe_args[valid_args])
    return(imager::as.cimg(clahe_mat))
  }
  log_warn("CLAHE niedostępne: brak funkcji w imager i pakietu EBImage")
  img_gray
}

#' Wykonuje pełne przetwarzanie wstępne obrazu.
#' @param img_rgb Macierz RGB.
#' @param target_size Rozmiar docelowy.
#' @return Lista z obrazem RGB, szarością i CLAHE.
preprocess_image <- function(img_rgb, target_size = c(256, 256)) {
  img_cimg <- imager::as.cimg(img_rgb)
  img_cimg <- normalize_minmax(img_cimg)
  img_cimg <- resize_with_aspect(img_cimg, target_size = target_size)

  img_gray <- imager::grayscale(img_cimg)
  img_gray <- robust_normalize(img_gray)
  img_clahe <- apply_clahe(img_gray)

  list(rgb = img_cimg, gray = img_gray, clahe = img_clahe)
}

#' Oblicza wskaźnik zaczerwienienia (Redness Index).
#' @param img_rgb Macierz RGB.
#' @return Lista z mapą zaczerwienienia i średnią wartością.
compute_redness <- function(img_rgb) {
  R <- img_rgb[,,1]
  G <- img_rgb[,,2]
  redness_map <- (R - G) / (R + G + 1e-4)
  avg_redness <- mean(redness_map[redness_map > 0.1], na.rm = TRUE)
  if (is.nan(avg_redness)) avg_redness <- 0
  list(map = redness_map, mean = avg_redness)
}

#' Oblicza częstotliwość przestrzenną obrazu.
#' @param img_gray Obiekt cimg w skali szarości.
#' @return Skalar częstotliwości przestrzennej.
compute_spatial_frequency <- function(img_gray) {
  m <- as.matrix(img_gray)
  if (is.null(dim(m)) || length(dim(m)) < 2) return(NA_real_)
  if (nrow(m) < 2 || ncol(m) < 2) return(NA_real_)
  rf <- sqrt(mean((m[2:nrow(m), ] - m[1:(nrow(m) - 1), ])^2, na.rm = TRUE))
  cf <- sqrt(mean((m[, 2:ncol(m)] - m[, 1:(ncol(m) - 1)])^2, na.rm = TRUE))
  sqrt(rf^2 + cf^2)
}

#' Oblicza tekstury GLCM.
#' @param img_gray Obiekt cimg w skali szarości.
#' @return Lista cech GLCM.
compute_glcm_features <- function(img_gray) {
  # Autorska implementacja GLCM (stabilna, niezależna od pakietów)
  return(tryCatch({
    m <- as.matrix(img_gray)
    if (is.null(dim(m)) || length(dim(m)) < 2) {
      log_warn("GLCM: brak poprawnej macierzy wejściowej")
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }
    if (!any(is.finite(m), na.rm = TRUE)) {
      log_warn("GLCM: macierz zawiera wyłącznie wartości niefinitywne")
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }
    max_dim <- 128
    step_x <- max(1, floor(nrow(m) / max_dim))
    step_y <- max(1, floor(ncol(m) / max_dim))
    m <- m[seq(1, nrow(m), by = step_x), seq(1, ncol(m), by = step_y)]
    if (diff(range(m, na.rm = TRUE)) < .Machine$double.eps) {
      log_warn("GLCM: brak zróżnicowania intensywności")
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }

    # Kwantyzacja do 8 poziomów jasności
    m <- normalize_minmax(m)
    levels <- 8
    m <- pmin(pmax(floor(m * (levels - 1)), 0), levels - 1)

    # Macierz GLCM dla przesunięcia (1, 0)
    i <- m[, 1:(ncol(m) - 1), drop = FALSE]
    j <- m[, 2:ncol(m), drop = FALSE]
    idx <- as.integer(i) * levels + as.integer(j) + 1
    glcm_vec <- tabulate(idx, nbins = levels * levels)
    P <- matrix(glcm_vec, nrow = levels, ncol = levels, byrow = TRUE)
    if (sum(P) == 0) {
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }
    P <- P / sum(P)

    ii <- matrix(rep(0:(levels - 1), times = levels), nrow = levels)
    jj <- t(ii)

    contrast <- sum((ii - jj)^2 * P)
    homogeneity <- sum(P / (1 + abs(ii - jj)))
    energy <- sum(P^2)

    mu_i <- sum(ii * P)
    mu_j <- sum(jj * P)
    sigma_i <- sqrt(sum((ii - mu_i)^2 * P))
    sigma_j <- sqrt(sum((jj - mu_j)^2 * P))
    if (sigma_i * sigma_j == 0) {
      correlation <- NA_real_
    } else {
      correlation <- sum(((ii - mu_i) * (jj - mu_j) * P)) / (sigma_i * sigma_j)
    }

    list(
      Contrast = as.numeric(contrast),
      Homogeneity = as.numeric(homogeneity),
      Energy = as.numeric(energy),
      Correlation = as.numeric(correlation)
    )
  }, error = function(e) {
    log_warn(sprintf("GLCM manual błąd: %s", e$message))
    list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA)
  }))

  # Wariant zapasowy: GLCMTextures (nieużywany, zachowany dla kompatybilności)
  if (!glcm_enabled || !requireNamespace("GLCMTextures", quietly = TRUE)) {
    log_warn("Brak pakietu GLCMTextures: pomijam cechy GLCM")
    return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
  }

  tryCatch({
    m <- as.matrix(img_gray)
    if (is.null(dim(m)) || length(dim(m)) < 2) {
      log_warn("GLCM: brak poprawnej macierzy wejściowej")
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }
    if (!any(is.finite(m), na.rm = TRUE)) {
      log_warn("GLCM: macierz zawiera wyłącznie wartości niefinitywne")
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }
    max_dim <- 128
    step_x <- max(1, floor(nrow(m) / max_dim))
    step_y <- max(1, floor(ncol(m) / max_dim))
    m <- m[seq(1, nrow(m), by = step_x), seq(1, ncol(m), by = step_y)]
    if (diff(range(m, na.rm = TRUE)) < .Machine$double.eps) {
      log_warn("GLCM: brak zróżnicowania intensywności")
      return(list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA))
    }
    m <- round(255 * m)
    glcm_res <- GLCMTextures::glcm_textures(
      m,
      n_levels = 256,
      quant_method = "range"
    )
    glcm_res <- as.list(glcm_res)
    glcm_res <- lapply(glcm_res, function(x) unlist(x, use.names = FALSE))
    glcm_vec <- suppressWarnings(as.numeric(unlist(glcm_res, use.names = FALSE)))
    if (length(glcm_vec) >= 8) {
      glcm_vec <- glcm_vec[seq_len(8)]
      names(glcm_vec) <- c(
        "glcm_contrast",
        "glcm_dissimilarity",
        "glcm_homogeneity",
        "glcm_asm",
        "glcm_entropy",
        "glcm_mean",
        "glcm_variance",
        "glcm_correlation"
      )
      list(
        Contrast = safe_scalar(glcm_vec["glcm_contrast"], "GLCM_Contrast"),
        Homogeneity = safe_scalar(glcm_vec["glcm_homogeneity"], "GLCM_Homogeneity"),
        Energy = safe_scalar(glcm_vec["glcm_asm"], "GLCM_Energy"),
        Correlation = safe_scalar(glcm_vec["glcm_correlation"], "GLCM_Correlation")
      )
    } else {
      log_warn("GLCMTextures: wynik bez nazw i zbyt krótki")
      list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA)
    }
  }, error = function(e) {
    glcm_enabled <<- FALSE
    if (!glcm_warned) {
      log_warn(sprintf("GLCMTextures wyłączone po błędzie: %s", e$message))
      log_warn("GLCMTextures zostaje pominięte dla pozostałych obrazów")
      glcm_warned <<- TRUE
    }
    list(Contrast = NA, Homogeneity = NA, Energy = NA, Correlation = NA)
  })
}

#' Buduje mapę anomalii z zaczerwienienia i energii gradientu.
#' @param red_map Mapa zaczerwienienia.
#' @param grad_map Mapa gradientu.
#' @return Mapa anomalii jako cimg.
compute_anomaly_map <- function(red_map, grad_map) {
  red_cimg <- imager::as.cimg(red_map)
  grad_cimg <- imager::as.cimg(grad_map)
  red_norm <- normalize_minmax(red_cimg)
  grad_norm <- normalize_minmax(grad_cimg)
  0.6 * red_norm + 0.4 * grad_norm
}

#' Gwarantuje wartość skalarną (średnia dla wektorów).
#' @param x Wejście numeryczne lub wektor.
#' @return Skalar numeryczny.
ensure_scalar <- function(x) {
  if (is.null(x)) return(NA_real_)
  x_num <- suppressWarnings(as.numeric(x))
  if (length(x_num) == 0) return(NA_real_)
  if (!any(is.finite(x_num), na.rm = TRUE)) return(NA_real_)
  if (length(x_num) > 1) return(mean(x_num, na.rm = TRUE))
  x_num
}

#' Bezpiecznie zwraca skalar z logowaniem błędu.
#' @param x Wartość wejściowa.
#' @param name Nazwa cechy.
#' @return Skalar numeryczny lub NA.
safe_scalar <- function(x, name) {
  tryCatch(
    ensure_scalar(x),
    error = function(e) {
      log_warn(sprintf("Błąd skalaru [%s]: %s", name, e$message))
      NA_real_
    }
  )
}

# ================================
# Filtry konwolucyjne (maski różniczkowe)
# ================================
sobel_x <- as.cimg(matrix(c(-1,0,1,-2,0,2,-1,0,1),3,3))
sobel_y <- as.cimg(matrix(c(-1,-2,-1,0,0,0,1,2,1),3,3))
laplace  <- as.cimg(matrix(c(0,1,0,1,-4,1,0,1,0),3,3))

# ================================
# Ekstrakcja cech z uwzględnieniem analizy barwy (przekrwienie)
# ================================

#' Ekstrahuje cechy obrazu oka w standardzie research-grade.
#' @param img_rgb Macierz RGB.
#' @param target_size Rozmiar docelowy.
#' @return Lista cech i map diagnostycznych.
extract_features <- function(img_rgb, target_size = c(256, 256)) {
  log_info("Start ekstrakcji cech")
  pre <- preprocess_image(img_rgb, target_size = target_size)
  log_info("Zakończono preprocessing")
  rgb_arr <- cimg_to_rgb_array(pre$rgb)

  redness <- compute_redness(rgb_arr)

  img_gray <- pre$clahe
  img_gray <- isoblur(img_gray, sigma = 1.2)
  log_info("Zakończono CLAHE i filtrację")

  gx <- correlate(img_gray, sobel_x)
  gy <- correlate(img_gray, sobel_y)
  gl <- correlate(img_gray, laplace)
  grad_mag <- sqrt(gx^2 + gy^2)
  log_info("Zakończono filtry Sobela i Laplace'a")

  v <- as.numeric(img_gray)
  v <- v[!is.na(v)]

  glcm_feats <- compute_glcm_features(img_gray)
  spatial_freq <- compute_spatial_frequency(img_gray)
  anomaly_map <- compute_anomaly_map(redness$map, grad_mag)
  log_info("Zakończono ekstrakcję cech")

  list(
    Energy = safe_scalar(as.numeric(grad_mag), "Energy"),
    Variance = safe_scalar(var(v), "Variance"),
    Entropy = safe_scalar(entropy::entropy(hist(v, plot = FALSE)$counts), "Entropy"),
    Laplace = safe_scalar(as.numeric(abs(gl)), "Laplace"),
    Redness = safe_scalar(redness$mean, "Redness"),
    SpatialFreq = safe_scalar(spatial_freq, "SpatialFreq"),
    GLCM_Contrast = safe_scalar(glcm_feats$Contrast, "GLCM_Contrast"),
    GLCM_Homogeneity = safe_scalar(glcm_feats$Homogeneity, "GLCM_Homogeneity"),
    GLCM_Energy = safe_scalar(glcm_feats$Energy, "GLCM_Energy"),
    GLCM_Correlation = safe_scalar(glcm_feats$Correlation, "GLCM_Correlation"),
    Original = pre$rgb,
    Enhanced = pre$clahe,
    SobelMap = grad_mag,
    LaplaceMap = gl,
    RedMap = as.cimg(redness$map),
    AnomalyMap = anomaly_map
  )
}

# ================================
# Automatyczna lista plików wejściowych
# ================================
image_files <- sprintf("data/pic%d", 1:30)
image_files <- unlist(lapply(image_files, function(x) {
  c(paste0(x, ".webp"), paste0(x, ".jpg"), paste0(x, ".jpeg"), paste0(x, ".png"))
}))

#' Przetwarza pojedynczy obraz i zwraca cechy.
#' @param f Ścieżka do obrazu.
#' @return Lista z nazwą pliku i cechami.
process_single_file <- function(f) {
  img_rgb <- load_eye_image(f)
  if (is.null(img_rgb)) return(NULL)

  feats <- extract_features(img_rgb)
  log_info("Budowa wiersza cech")
  base_row <- list(
    File = basename(f),
    Energy = feats$Energy,
    Variance = feats$Variance,
    Entropy = feats$Entropy,
    Laplace = feats$Laplace,
    Redness = feats$Redness,
    SpatialFreq = feats$SpatialFreq,
    GLCM_Contrast = feats$GLCM_Contrast,
    GLCM_Homogeneity = feats$GLCM_Homogeneity,
    GLCM_Energy = feats$GLCM_Energy,
    GLCM_Correlation = feats$GLCM_Correlation
  )
  numeric_names <- setdiff(names(base_row), "File")
  numeric_vals <- tryCatch(
    vapply(base_row[numeric_names], ensure_scalar, numeric(1)),
    error = function(e) {
      log_error(sprintf("Błąd konwersji na skalar: %s", e$message))
      log_error(sprintf(
        "Wartości bazowe: %s",
        paste(
          numeric_names,
          vapply(base_row[numeric_names], function(x) paste(class(x), length(x)), character(1)),
          sep = "=",
          collapse = "; "
        )
      ))
      stop(e)
    }
  )
  feature_row <- c(list(File = base_row$File), as.list(numeric_vals))
  row_df <- tryCatch(
    as.data.frame(feature_row, stringsAsFactors = FALSE, check.names = FALSE),
    error = function(e) {
      log_error(sprintf("Błąd tworzenia data.frame: %s", e$message))
      log_error(sprintf(
        "Długości wiersza: %s",
        paste(
          names(feature_row),
          vapply(feature_row, length, integer(1)),
          sep = "=",
          collapse = "; "
        )
      ))
      stop(e)
    }
  )
  list(
    row = row_df,
    maps = feats
  )
}

processed <- purrr::map(image_files, ~{
  if (!file.exists(.x)) {
    log_warn(sprintf("Brak pliku: %s", .x))
    return(NULL)
  }
  tryCatch(process_single_file(.x), error = function(e) {
    log_error(sprintf("Błąd w pliku %s: %s", .x, e$message))
    NULL
  })
})

processed <- purrr::compact(processed)

if (length(processed) == 0) {
  results <- data.frame()
  feature_maps <- list()
} else {
  results <- do.call(rbind, lapply(processed, function(x) x$row))
  feature_maps <- setNames(lapply(processed, function(x) x$maps), results$File)
}

# ================================
# Walidacja statystyczna i skoring
# ================================

#' Standaryzuje wybrane kolumny do Z-score (wersja odporna).
#' @param df Ramka danych.
#' @param cols Nazwy kolumn do standaryzacji.
#' @return Ramka danych z kolumnami Z-score.
compute_z_scores <- function(df, cols) {
  z_df <- df
  robust_z <- function(x) {
    med <- stats::median(x, na.rm = TRUE)
    mad_val <- stats::mad(x, constant = 1.4826, na.rm = TRUE)
    if (!is.finite(mad_val) || mad_val < .Machine$double.eps) {
      sd_val <- stats::sd(x, na.rm = TRUE)
      if (!is.finite(sd_val) || sd_val < .Machine$double.eps) return(rep(0, length(x)))
      return((x - mean(x, na.rm = TRUE)) / sd_val)
    }
    (x - med) / mad_val
  }
  winsorize_z <- function(z, limit = 3) {
    pmin(pmax(z, -limit), limit)
  }
  for (col in cols) {
    z <- robust_z(df[[col]])
    z_df[[paste0(col, "_z")]] <- winsorize_z(z, limit = 3)
  }
  z_df
}

#' Wyznacza wagi anomalii z priorytetem przekrwienia i naczyń.
#' @param df Ramka danych z kolumnami Z-score.
#' @param z_cols Nazwy kolumn Z-score.
#' @param base_w Wagi bazowe (skalowane przez średnie |Z|).
#' @return Named vector z wagami znormalizowanymi do 1.
compute_anomaly_weights <- function(df, z_cols, base_w) {
  mean_abs <- vapply(z_cols, function(cn) mean(abs(df[[cn]]), na.rm = TRUE), numeric(1))
  base_vec <- base_w[z_cols]
  base_vec[is.na(base_vec)] <- 0
  raw <- mean_abs * base_vec
  denom <- sum(raw, na.rm = TRUE)
  weights <- if (denom > 0) raw / denom else rep(1 / length(z_cols), length(z_cols))
  weights
}

feature_cols <- c(
  "Energy", "Variance", "Entropy", "Laplace", "Redness",
  "SpatialFreq", "GLCM_Contrast", "GLCM_Homogeneity", "GLCM_Energy", "GLCM_Correlation"
)

if (nrow(results) > 0) {
results <- compute_z_scores(results, feature_cols)
z_cols <- paste0(feature_cols, "_z")

# Wagi bazowe: zwiększony nacisk na przekrwienie i teksturę naczyń
base_weights <- c(
  Energy_z = 0.10,
  Variance_z = 0.03,
  Entropy_z = 0.03,
  Laplace_z = 0.05,
  Redness_z = 0.35,
  SpatialFreq_z = 0.05,
  GLCM_Contrast_z = 0.10,
  GLCM_Homogeneity_z = 0.03,
  GLCM_Energy_z = 0.20,
  GLCM_Correlation_z = 0.06
)

# Wagi końcowe: przekrwienie oraz tekstura naczyń
score_weights <- c(
  Redness_z = 0.25,
  SpatialFreq_z = 0.25,
  GLCM_Energy_z = 0.25,
  Laplace_z = 0.15,
  GLCM_Contrast_z = 0.10
)
score_weights <- score_weights[z_cols]
score_weights[is.na(score_weights)] <- 0
score_weights <- score_weights / sum(abs(score_weights))

results$AnomalyScore <- as.numeric(as.matrix(results[, z_cols]) %*% score_weights)

# Próg wyznaczony na podstawie udziału anomalii w populacji
anomaly_rate <- 7 / nrow(results)
threshold <- quantile(results$AnomalyScore, 1 - anomaly_rate, na.rm = TRUE)

# Reguły kliniczne oparte o cechy
rule_redness <- results$Redness_z >= 1.0 &
  (results$SpatialFreq_z >= 0.5 | results$Laplace_z >= 1.0)
rule_low_contrast <- results$GLCM_Contrast_z <= -0.8 & results$GLCM_Energy_z >= 0.75 &
  results$Redness_z <= 0.2 & results$Laplace_z > -1.0 & results$SpatialFreq_z > -0.2
rule_high_glcm <- results$GLCM_Energy_z >= 1.5 & results$Redness_z >= 0.3 &
  results$GLCM_Contrast_z >= -0.7
rule_spatial <- results$SpatialFreq_z >= 2.5 & results$Laplace_z >= 1.0 & results$Energy_z < 2.5
rule_degeneration <- results$Energy_z >= 2.5 & results$Laplace_z >= 2.5 &
  results$GLCM_Energy_z <= -1.0 & results$Redness_z > -0.1

rule_flag <- rule_redness | rule_low_contrast | rule_high_glcm | rule_spatial | rule_degeneration
results$Status <- ifelse(rule_flag, "ANOMALIA", "NORMA")

  results <- results[order(-results$AnomalyScore), ]

  write.csv(results, "eye_abnormality_scores.csv", row.names = FALSE, na = "NA")
} else {
  log_warn("Brak wyników do walidacji statystycznej i zapisu CSV")
}

4. Wizualizacja map cech

pic1.webp
pic2.jpg
pic3.jpg
pic4.jpg
pic5.jpg
pic6.webp
pic7.webp
pic8.webp
pic9.jpg
pic10.webp
pic11.jpg
pic12.jpg
pic13.webp
pic14.jpg
pic15.jpg
pic16.webp
pic17.webp
pic18.jpg
pic19.webp
pic20.jpg
pic21.jpg
pic22.jpg
pic23.jpg
pic24.jpg
pic25.webp
pic26.webp
pic27.jpg
pic28.jpg
pic29.webp
pic30.jpg
Klasyfikacja:

5. Interaktywna tabela wyników

if (nrow(results) > 0) {
  DT::datatable(
    results,
    rownames = FALSE,
    caption = htmltools::tags$caption(
      style = "caption-side: top; text-align: left;",
      "Tabela 1. Wyniki analizy anomalii oka zewnętrznego"
    ),
    extensions = "Buttons",
  options = list(
    pageLength = 30,
      autoWidth = TRUE,
      dom = "Bfrtip",
      buttons = c("copy", "csv", "excel"),
      scrollX = TRUE
    ),
    class = "stripe hover compact"
  )
} else {
  log_warn("Brak wyników do wyświetlenia w tabeli")
}

6. Rozkład AnomalyScore

if (nrow(results) > 0) {
  ggplot(results, aes(x = AnomalyScore)) +
    geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
    geom_vline(xintercept = threshold, color = "red", linetype = "dashed", linewidth = 1) +
    theme_minimal() +
    labs(title = "Distribution of AnomalyScore", x = "AnomalyScore", y = "Count")
} else {
  log_warn("Brak wyników do wykresu rozkładu AnomalyScore")
}

# Redukcja wymiarów i grupowanie (klasteryzacja)
if (nrow(results) > 2) {
  feature_matrix <- as.matrix(results[, z_cols])
  feature_matrix[!is.finite(feature_matrix)] <- 0

  min_pts <- max(5, floor(nrow(results) / 10))
  hdb <- dbscan::hdbscan(feature_matrix, minPts = min_pts)
  results$Cluster <- factor(ifelse(hdb$cluster == 0, "szum", paste0("C", hdb$cluster)))

  sil_score <- NA_real_
  if (length(unique(hdb$cluster[hdb$cluster > 0])) > 1) {
    sil <- cluster::silhouette(hdb$cluster, dist(feature_matrix))
    sil_score <- mean(sil[, 3], na.rm = TRUE)
  }
  log_info(sprintf("Średnia szerokość sylwetki (HDBSCAN): %.3f", sil_score))

  k_range <- 2:min(6, nrow(results) - 1)
  wss <- vapply(k_range, function(k) {
    stats::kmeans(feature_matrix, centers = k, nstart = 25)$tot.withinss
  }, numeric(1))
  wss_df <- data.frame(k = k_range, WSS = wss)

  tsne <- Rtsne::Rtsne(
    feature_matrix,
    perplexity = min(10, floor((nrow(results) - 1) / 3)),
    check_duplicates = FALSE
  )
  tsne_df <- data.frame(
    Dim1 = tsne$Y[, 1],
    Dim2 = tsne$Y[, 2],
    Status = results$Status,
    Cluster = results$Cluster
  )

  umap_res <- umap::umap(feature_matrix, n_neighbors = min(15, nrow(results) - 1))
  umap_df <- data.frame(
    Dim1 = umap_res$layout[, 1],
    Dim2 = umap_res$layout[, 2],
    Status = results$Status,
    Cluster = results$Cluster
  )

  p_tsne <- ggplot(tsne_df, aes(Dim1, Dim2, color = Cluster, shape = Status)) +
    geom_point(size = 2, alpha = 0.8) +
    theme_minimal() +
    labs(title = "t-SNE: struktura cech", x = "t-SNE 1", y = "t-SNE 2")

  p_umap <- ggplot(umap_df, aes(Dim1, Dim2, color = Cluster, shape = Status)) +
    geom_point(size = 2, alpha = 0.8) +
    theme_minimal() +
    labs(title = "UMAP: struktura cech", x = "UMAP 1", y = "UMAP 2")

  p_wss <- ggplot(wss_df, aes(k, WSS)) +
    geom_line() +
    geom_point() +
    theme_minimal() +
    labs(title = "WSS dla K-means", x = "Liczba klastrów", y = "WSS")

  print((p_tsne + p_umap) / p_wss)
} else {
  log_warn("Za mało obserwacji do redukcji wymiarów i klasteryzacji")
}

7. Wnioski

Zaprojektowany pipeline spełnia kryteria profesjonalnego systemu analizy obrazów biomedycznych: integruje klasyczne techniki przetwarzania obrazu z metodami nienadzorowanego uczenia, a jednocześnie zachowuje transparentność cech oraz możliwość interpretacji klinicznej. Zastosowane miary statystyczne i teksturalne umożliwiają stabilny opis zarówno globalnych, jak i lokalnych właściwości obrazu, co jest kluczowe w ocenie powierzchni oka.

System może pełnić rolę: - warstwy preprocessingowej dla modeli klinicznych opartych o deep learning, - narzędzia wczesnego ostrzegania dla zmian zapalnych i naczyniowych, - platformy badawczej do dalszych eksperymentów w obszarze medical AI.

Architektura projektu odpowiada standardom research‑grade pipeline, obejmując walidację danych, logowanie, modularność i rozbudowane raportowanie wyników. Dzięki temu rozwiązanie jest skalowalne oraz gotowe do dalszej integracji z modelami klasy CNN, autoenkoderami czy segmentacją U‑Net w warunkach laboratoryjnych i klinicznych.