Use ImageArray package for images (and labels)

.Rbuildignore

@@ -6,3 +6,5 @@ docs/
 vignettes/*html
 vignettes/*cache
 vignettes/*files
+^\.positai$
+^\.claude$

.gitignore

@@ -7,3 +7,4 @@ R/_*
 vignettes/*html
 vignettes/*cache
 vignettes/*files
+.positai

DESCRIPTION

@@ -45,7 +45,8 @@ Imports:
     SingleCellExperiment,
     SummarizedExperiment,
     SparseArray,
-    ZarrArray
+    ZarrArray,
+    ImageArray
 Suggests:
     BiocStyle,
     EBImage,

NAMESPACE

@@ -124,6 +124,7 @@ importFrom(BiocGenerics,table)
 importFrom(BiocGenerics,transform)
 importFrom(DBI,dbIsValid)
 importFrom(DelayedArray,DelayedArray)
+importFrom(ImageArray,ImageArray)
 importFrom(Matrix,sparseMatrix)
 importFrom(Matrix,sparseVector)
 importFrom(Matrix,summary)
@@ -209,5 +210,3 @@ importFrom(sf,st_polygon)
 importFrom(sf,st_sf)
 importFrom(sf,st_sfc)
 importFrom(utils,.DollarNames)
-importFrom(utils,head)
-importFrom(utils,tail)

R/AllClasses.R

@@ -3,26 +3,26 @@
 #' @export
 #' @rdname SpatialData
 .SpatialData <- setClass(
-    Class="SpatialData",
-    contains=c("list", "Annotated"),
-    representation(
-        images="list",  # 'SpatialDataImage's
-        labels="list",  # 'SpatialDataLabel's
-        points="list",  # 'SpatialDataPoint's
-        shapes="list",  # 'SpatialDataShape's
-        tables="list")) # 'SingleCellExperiment's
+  Class="SpatialData",
+  contains=c("list", "Annotated"),
+  representation(
+    images="list",  # 'SpatialDataImage's
+    labels="list",  # 'SpatialDataLabel's
+    points="list",  # 'SpatialDataPoint's
+    shapes="list",  # 'SpatialDataShape's
+    tables="list")) # 'SingleCellExperiment's
 
 .LAYERS <- `names<-`(. <- c("images","labels","points","shapes","tables"), .)
 .SpatialDataAttrs <- setClass("SpatialDataAttrs", contains="list")
 setOldClass("duckspatial_df")
 
 setClass("SpatialDataArray", 
-    contains=c("Annotated", "VIRTUAL"),
-    slots=list(data="list", meta="SpatialDataAttrs"))
+         contains=c("Annotated", "VIRTUAL"),
+         slots=list(data="ImageArray", meta="SpatialDataAttrs"))
 
 setClass("SpatialDataFrame",
-    contains=c("Annotated", "VIRTUAL"),
-    slots=list(data="duckspatial_df", meta="SpatialDataAttrs"))
+         contains=c("Annotated", "VIRTUAL"),
+         slots=list(data="duckspatial_df", meta="SpatialDataAttrs"))
 
 .SpatialDataImage <- setClass("SpatialDataImage", contains="SpatialDataArray")
 .SpatialDataLabel <- setClass("SpatialDataLabel", contains="SpatialDataArray")
@@ -31,5 +31,5 @@ setClass("SpatialDataFrame",
 .SpatialDataShape <- setClass("SpatialDataShape", contains="SpatialDataFrame")
 
 setClassUnion("SpatialDataElement", c(
-    "SpatialDataImage", "SpatialDataLabel", 
-    "SpatialDataPoint", "SpatialDataShape"))
+  "SpatialDataImage", "SpatialDataLabel", 
+  "SpatialDataPoint", "SpatialDataShape"))

R/crop.R

@@ -168,7 +168,8 @@ setMethod("crop", "SpatialDataArray", \(x, y, j=1, ...) {
     .check_box(y)
     z <- .box2rev(x, y, j)
     # offset current origin
-    wh <- metadata(x)$wh
+    # wh <- metadata(x)$wh
+    wh <- extent(x)
     if (!is.null(wh)) {
         z$xmin <- z$xmin - wh[[1]][1]
         z$xmax <- z$xmax - wh[[1]][1]
@@ -181,22 +182,14 @@ setMethod("crop", "SpatialDataArray", \(x, y, j=1, ...) {
     z$ymin <- floor(max(z$ymin, 0))
     z$xmax <- ceiling(min(z$xmax, d[n]))
     z$ymax <- ceiling(min(z$ymax, d[n-1]))
-    # update origin
-    if (is.null(wh)) {
-        # set from bounding box
-        wh <- list(
-            c(z$xmin, z$xmax), 
-            c(z$ymin, z$ymax))
-    } else {
-        # offset current origin
-        wh[[1]] <- wh[[1]][1] + c(z$xmin, z$xmax)
-        wh[[2]] <- wh[[2]][1] + c(z$ymin, z$ymax)
-    }
-    metadata(x)$wh <- wh
     # subset array
     i <- seq(z$ymin+1, z$ymax)
     j <- seq(z$xmin+1, z$xmax) 
-    if (n == 3) x[, i, j] else x[i, j]
+    x <- if (n == 3) x[, i, j] else x[i, j]
+    # update origin
+    # TODO: what is the order of translation
+    x <- translation(x, t = c(z$ymin, z$xmin))
+    x
 })
 
 #' @export

R/extent.R

@@ -37,7 +37,7 @@ setMethod("extent", "SpatialData", \(x, i=1) {
 #' @rdname extent
 setMethod("extent", "SpatialDataArray", \(x, i=1) {
     x <- transform(x, i)
-    wh <- metadata(x)$wh %||% {
+    wh <- ImageArray::extent(data(x, NULL)) %||% {
         n <- length(d <- dim(x))
         if (n == 3) d <- d[-1]
         d <- rev(d)

R/misc.R

@@ -88,9 +88,9 @@ setMethod("show", "SpatialData", .showSpatialData)
 
 #' @importFrom S4Vectors coolcat
 .showArray <- function(object) {
-    n.object <- length(object@data)
+    n.object <- length(levels(object@data))
     cat("class: ", class(object), ifelse(n.object > 1, "(MultiScale)", ""),"\n")
-    scales <- vapply(object@data, \(x) paste0(dim(x), collapse=","), character(1))
+    scales <- vapply(levels(object@data), \(x) paste0(dim(x), collapse=","), character(1))
     coolcat("Scales (%d): (%s)", scales)
 }
 

R/read.R

@@ -47,17 +47,15 @@ NULL
 .readArray <- function(x, ...) {
     md <- read_zarr_attributes(x)
     mdattr <- SpatialDataAttrs(md)
-    # TODO: paths to datasets have to be validated properly in the future
-    # https://ngff.openmicroscopy.org/specifications/0.5/index.html#images
-    # The name of the array is arbitrary with the ordering defined by
-    # by the "multiscales" metadata, but is often a sequence starting at 0.
     ds <- .validate_multiscales_paths(x, datasets(mdattr))
     ds <- file.path(x, as.character(ds))
-    as <- lapply(ds, ZarrArray)
+    as <- ImageArray(levels = lapply(ds, ZarrArray),
+                     meta = list(axes = vapply(axes(mdattr), \(.) .$name, character(1))))
     list(array=as, mdattr=mdattr)
 }
 
 #' @rdname readSpatialData
+#' @importFrom ImageArray ImageArray
 #' @export
 readImage <- function(x, ...) {
     l <- .readArray(x, ...)

R/sdArray.R

@@ -67,8 +67,8 @@ NULL
 #' @rdname SpatialDataArray
 #' @importFrom methods new
 #' @importFrom S4Vectors metadata<-
-SpatialDataImage <- function(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
-    x <- .SpatialDataImage(data=data, meta=meta, ...)
+SpatialDataImage <- function(data = list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
+    x <- .SpatialDataImage(data=as_imgarray(data, 3), meta=meta, ...)
     metadata(x) <- metadata
     return(x)
 }
@@ -77,12 +77,24 @@ SpatialDataImage <- function(data=list(), meta=SpatialDataAttrs(), metadata=list
 #' @rdname SpatialDataArray
 #' @importFrom methods new
 #' @importFrom S4Vectors metadata<-
-SpatialDataLabel <- function(data=list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
-    x <- .SpatialDataLabel(data=data, meta=meta, ...)
+SpatialDataLabel <- function(data = list(), meta=SpatialDataAttrs(), metadata=list(), ...) {
+    x <- .SpatialDataLabel(data=as_imgarray(data, 2), meta=meta, ...)
     metadata(x) <- metadata
     return(x)
 }
 
+#' @noRd
+as_imgarray <- function(data, dim){
+  if(is(data, "ImageArray")){
+    data
+  } else {
+    if(!is.list(data))
+      data <- list(data)
+    ImageArray(levels = data, 
+               meta = list(axes = c(if(dim == 3) "c" else NULL, "y", "x")))
+  }
+}
+
 # utils ----
 
 #' @rdname SpatialDataArray
@@ -146,58 +158,11 @@ setMethod("channels", "SpatialDataElement", \(x, ...) stop("only 'images' have c
 
 # sub ----
 
-.check_jk <- \(x, .) {
-    if (isTRUE(x)) return()
-    tryCatch(
-        stopifnot(
-            is.numeric(x), x == round(x),
-            diff(range(x)) == length(x)-1,
-            (y <- abs(x)) == seq(min(y), max(y))
-        ),
-        error=\(e) stop(sprintf("invalid '%s'", .))
-    )
-}
-
-#' @exportMethod [
-#' @rdname SpatialDataArray
-#' @importFrom utils head tail
-setMethod("[", "SpatialDataImage", \(x, i, j, k, ..., drop=FALSE) {
-    if (missing(i)) i <- TRUE
-    if (missing(j)) j <- TRUE else if (isFALSE(j)) j <- 0 else .check_jk(j, "j")
-    if (missing(k)) k <- TRUE else if (isFALSE(k)) k <- 0 else .check_jk(k, "k")
-    ijk <- list(i, j, k)
-    n <- length(data(x, NULL))
-    d <- dim(data(x))
-    data(x) <- lapply(seq_len(n), \(.) {
-        j <- if (isTRUE(j)) seq_len(d[2]) else j
-        k <- if (isTRUE(k)) seq_len(d[3]) else k
-        jk <- lapply(list(j, k), \(jk) {
-            fac <- 2^(.-1)
-            seq(floor(head(jk, 1)/fac), 
-                ceiling(tail(jk, 1)/fac))
-        })
-        data(x, .)[i, jk[[1]], jk[[2]], drop=FALSE]
-    })
-    x
-})
-
 #' @exportMethod [
 #' @rdname SpatialDataArray
-#' @importFrom utils head tail
-setMethod("[", "SpatialDataLabel", \(x, i, j, ..., drop=FALSE) {
-    if (missing(i)) i <- TRUE else if (isFALSE(i)) i <- 0 else .check_jk(i, "i")
-    if (missing(j)) j <- TRUE else if (isFALSE(j)) j <- 0 else .check_jk(j, "j")
-    n <- length(data(x, NULL))
-    d <- dim(data(x, 1))
-    data(x) <- lapply(seq_len(n), \(.) {
-        i <- if (isTRUE(i)) seq_len(d[1]) else i
-        j <- if (isTRUE(j)) seq_len(d[2]) else j
-        ij <- lapply(list(i, j), \(ij) {
-            fac <- 2^(.-1)
-            seq(floor(head(ij, 1)/fac), 
-                ceiling(tail(ij, 1)/fac))
-        })
-        data(x, .)[ij[[1]], ij[[2]], drop=FALSE]
-    })
+setMethod("[", "SpatialDataArray", \(x, i, j,..., drop=FALSE) {
+    cl <- sys.call()
+    cl[[2]] <- substitute(data(x, NULL))
+    data(x) <- eval(cl, parent.frame())
     x
-})
+})

R/trans.R

@@ -86,7 +86,8 @@ setMethod("sequence", "SpatialDataElement", \(x, t, ..., rev=FALSE) {
 .mirror <- \(x, t, k=1) {
     d <- length(dim(x)) == 3
     i <- if (d) c(1, 3, 2) else c(2, 1)
-    data(x) <- list(aperm(data(x, k), i))
+    # data(x) <- list(aperm(data(x, k), i))
+    data(x) <- ImageArray::aperm(data(x, NULL), perm = i)
     rotate(x, t, k=1)
 }
 
@@ -97,11 +98,13 @@ setMethod("mirror", "SpatialDataArray", \(x, t=c("v", "h"), k=1, ...)
 
 #' @export
 #' @rdname trans
-setMethod("flip", "SpatialDataArray", \(x, k=1, ...) .mirror(x, -90, k))
+setMethod("flip", "SpatialDataArray", \(x, ...) .mirror(x, 270))
+# TODO: allow -90 as angle in ImageArray
+# setMethod("flip", "SpatialDataArray", \(x, ...) .mirror(x, -90))
 
 #' @export
 #' @rdname trans
-setMethod("flop", "SpatialDataArray", \(x, k=1, ...) .mirror(x, 90, k))
+setMethod("flop", "SpatialDataArray", \(x, ...) .mirror(x, 90))
 
 # rotation matrix to rotate points counter-clockwise through an angle 't'
 .R <- \(t) matrix(c(cos(t), -sin(t), sin(t), cos(t)), 2, 2)
@@ -111,58 +114,58 @@ setMethod("flop", "SpatialDataArray", \(x, k=1, ...) .mirror(x, 90, k))
 #' @importFrom methods as
 #' @importFrom BiocGenerics rotate
 #' @importFrom S4Vectors metadata<-
-setMethod("rotate", "SpatialDataArray", \(x, t, k=1, ..., rev=FALSE) {
-    if (!requireNamespace("EBImage", quietly=TRUE))
-        stop("install 'EBImage' to use this function")
-    # negate angle since 'EBImage' rotates clockwise
+setMethod("rotate", "SpatialDataArray", \(x, t, ..., rev=FALSE) {
+  # complement angle with 360 to turn counterclockwise
     stopifnot(length(t) == 1, is.finite(t))
     if (t %% 360 == 0) return(x)
-    if (rev) t <- -t
-    if (length(d <- dim(data(x, k))) == 3) d <- d[-1]
-    metadata(x)$wh <- lapply(rev(d), \(.) c(c(0, .) %*% .R(t*pi/180)))
-    f <- \(.) EBImage::rotate(., -t) 
-    a <- f(aperm(as.array(data(x, k))))
-    metadata(x)$data_type <- data_type(x)
-    data(x) <- list(as(aperm(a), "SparseArray"))
+    if (rev) t <- 360-t
+    data(x) <- ImageArray::rotate(data(x, NULL), t)
     return(x)
 })
 
-.trans_a <- \(x, t, f=c("scale", "translation"), k=1, rev=FALSE) {
-    f <- match.arg(f)
-    n <- length(d <- dim(data(x, k)))
-
-    # setup: identity, operator
-    map <- list(
-        ids=c(scale=1, translation=0),
-        ops=c(scale="*", translation="+"))
-
-    # validation & identity check
-    stopifnot(is.numeric(t), is.finite(t), length(t) == n)
-    if (all(t == map$ids[f])) return(x)
-    if (rev) t <- if (f == "scale") 1/t else -t
-
-    # project to spatial (XY) dims
-    if (n == 3) { t <- t[-1]; d <- d[-1] }
-    t <- rev(t); d <- rev(d)
-
-    # update 'wh' metadata
-    wh <- metadata(x)$wh %||% list(c(0, d[1]), c(0, d[2]))
-    op <- get(map$ops[f])
-    metadata(x)$wh <- mapply(op, t, wh, SIMPLIFY=FALSE)
-    return(x)
+.trans_a_scale <- \(x, t, rev=FALSE) {
+  n <- length(d <- dim(x))
+
+  # validation & identity check
+  stopifnot(is.numeric(t), is.finite(t), length(t) == 2)
+  if (all(t == 0)) return(x)
+
+  # scale
+  if (rev) t <- 1/t
+
+  # project to spatial (XY) dims
+  if (n == 3) { d <- d[-1] }
+  t <- rev(t); d <- rev(d)
+
+  data(x) <- ImageArray::scale(data(x, NULL), 
+                               output.dim = as.integer(round(d*t)))
+  x
 }
 
 #' @export
 #' @rdname trans
 #' @importFrom BiocGenerics scale
 setMethod("scale", "SpatialDataArray",
-    \(x, t, ...) .trans_a(x, t, "scale", ...))
+    \(x, t, ...) .trans_a_scale(x, t, ...))
+
+.trans_a_trans <- \(x, t, rev=FALSE) {
 
+  # validation & identity check
+  stopifnot(is.numeric(t), is.finite(t), length(t) == 2)
+  if (all(t == 0)) return(x)
+
+  # project to spatial (XY) dims
+  t <- rev(t)
+
+  data(x) <- ImageArray::translation(data(x, NULL), shift = t)
+  x
+}
+
 #' @export
 #' @rdname trans
 setMethod("translation", 
     c("SpatialDataArray", "numeric"), 
-    \(x, t, ...) .trans_a(x, t, "translation", ...))
+    \(x, t, ...) .trans_a_trans(x, t, ...))
 
 # point/shape ----