R/vis_volume.R

In fsbrain: Managing and Visualizing Brain Surface Data

# Functions for volume manipulation and rendering.
# The 3D imaging functions are designed to work on gray-scale (single channel) images.


#' @title Extract a slice of a 3D image stack.
#'
#' @description Extracts one or more 2D slices from a 3D image (or a frame of a 4D image). To display the result, you can use \code{\link[fsbrain]{volvis.lightbox}}.
#'
#' @param volume a 3D or 4D image volume. Note that empty dimensions will be dropped before any processing, and the remaining volume must have 3 or 4 dimensions.
#'
#' @param slice_index positive integer or vector of positive integers, the index into the slices (for the axis). A *slice* in the sense of this function is any 2D image plane extracted from the 3D volume (no matter the axis). If NULL, the slice in the middle of the volume is used. One can pass the magic character string 'all' to use all slice indices along the axis.
#'
#' @param frame positive integer, optional. The frame (time point) to use, only relevant for 4D volumes. The last (i.e. 4th) dimension is assumed to be the time dimension in that case.
#'
#' @param axis positive integer, the axis to use when indexing the slices. Defaults to 1.
#'
#' @param rotation integer, rotation in degrees. Defaults to 0 (no ratation). Must be a multiple of 90L if given.
#'
#' @param flip NULL or one of the character strings 'vertically' or 'horizontally'. Note that flipping *horizontally* means that the image will be mirrored along the central *vertical* axis. If `NULL` is passed, nothing is flipped. Flipping occurs after rotation.
#'
#' @return slice data. If `slice_index` is a scalar, a numerical 2D matrix (a 2D image from the stack). Otherwise, a numerical 3D array that contains the selected 2D images.
#'
#' @family volume utility
#'
#' @export
vol.slice <- function(volume, slice_index=NULL, frame=1L, axis=1L, rotation=0L, flip=NULL) {
    if(axis < 1 | axis > 3) {
        stop(sprintf("Axis must be integer with value 1, 2 or 3 but is %d.\n", axis));
    }

    if(length(dim(volume))==4) {
        vol3d = volume[,,,frame];
    } else if(length(dim(volume))==3) {
        vol3d = volume;
    } else {
        stop("Data passed as parameter 'volume' must have 3 or 4 dimensions.");
    }

    if(is.null(slice_index)) {
        # Select a middle slice, the first one is often (almost) empty.
        slice_index = as.integer(round(dim(vol3d)[axis] / 2));
    }

    if(is.character(slice_index)) {
        if(slice_index == 'all') {
            slice_index = seq_len(dim(vol3d)[axis]);
        }
    }

    if(!is.numeric(slice_index)) {
        stop("Could not determine a valid slice index.");
    }

    # Limit the index to the range 1..axis dim
    slice_index = ifelse(slice_index < 1L, 1L, slice_index);
    slice_index = ifelse(slice_index > dim(vol3d)[axis], dim(vol3d)[axis], slice_index);

    # Select requested axis
    if(axis == 1L) {
        slice = vol3d[slice_index,,];
    } else if(axis == 2L) {
        slice = vol3d[,slice_index,];
    } else {
        slice = vol3d[,,slice_index];
    }

    if(rotation != 0L) {
        if(length(slice_index) == 1) {
            slice = freesurferformats::rotate2D(slice, rotation);
        } else {
            slice = freesurferformats::rotate3D(slice, axis=axis, degrees=rotation);
        }
    }

    if(!is.null(flip)) {
        if(length(slice_index) == 1) {
            slice = freesurferformats::flip2D(slice, how = flip);
        } else {
            slice = freesurferformats::flip3D(slice, axis = axis, how = flip);
        }
    }

    return(slice);
}




#' @title Compute foreground pixels over the whole 3D imagestack.
#'
#' @description Compute, over all images in a stack along an axis, the foreground and background pixels as a binary mask. A pixel is a `foreground` pixel iff its value is greater than the `threshold` parameter in at least one of the slices. A pixel is a `background` pixel iff its value is below or euqal to the `threshold` in all slices.
#'
#' @param volume a 3D image volume
#'
#' @param plane integer vector of length 2 or something that will be turned into one by \code{\link[fsbrain]{vol.plane.axes}}.
#'
#' @param threshold numerical, the threshold intensity used to separate background and foreground. All voxels with intensity values greater than this value will be considered `foreground` voxels.
#'
#' @return integer 2D matrix with dimensions of a slice of the volume. Positions set to 1 are `foreground` pixels and positions set to 0 are `background` pixels (see `Details` section).
#'
#' @keywords internal
vol.boundary.mask <- function(volume, plane=1L, threshold=0L) {
    axes = vol.plane.axes(plane);
    if(length(dim(volume)) != 3) {
        stop(sprintf("Volume must have exactly 3 dimensions but has %d.\n", length(dim(volume))));
    }

    height = dim(volume)[axes[1]];
    width = dim(volume)[axes[2]];

    if(is.null(threshold)) {
        return(matrix(rep(1L, width*height), nrow=height));
    }

    maximg = apply(volume, axes, max);
    foreground_indices = which(maximg > threshold, arr.ind=TRUE);

    mask = matrix(rep(0L, width*height), nrow=height);
    mask[foreground_indices] = 1L;
    #cat(sprintf(" Using axes '%s', dimension of mask is '%s'.\n", paste(axes, collapse=", "), paste(dim(mask), collapse=", ")));
    return(mask);
}


#' @title Compute 3D bounding box of a volume.
#'
#' @description Compute the axis-aligned foreground bounding box of a 3D volume, i.e. the inner foreground area that must be retained if you want to remove all background from the corners of the volume. The foreground is determined by thresholding, such that all values greater than 0 are considered foreground. See \code{\link[fsbrain]{vol.boundary.mask}} for details.
#'
#' @param volume a 3D image volume
#'
#' @param threshold numerical, the threshold intensity used to separate background and foreground. All voxels with intensity values greater than this value will be considered `foreground` voxels.
#'
#' @param apply logical, whether to directly apply the bounding box and return the resulting volume instead.
#'
#' @return named list with 2 entries: `from` is an integer vector of length 3, defining the minimal (x,y,z) foreground indices. `to` is an integer vector of length 3, defining the maximal (x,y,z) foreground indices.
#'
#' @family volume utility
#'
#' @export
vol.boundary.box <- function(volume, threshold=0L, apply=FALSE) {
    if(length(dim(volume)) != 3) {
        stop(sprintf("Volume must have exactly 3 dimensions but has %d.\n", length(dim(volume))));
    }
    min_index_per_axis = rep(-1L, 3);
    max_index_per_axis = rep(-1L, 3);
    for(axis in c(1L, 2L, 3L)) {
        mmask = vol.boundary.mask(volume, axis, threshold=threshold);
        colmax = apply(mmask, 1, max);
        min_index_per_axis[axis] = Position(function(x) x >= 1L, colmax, right=FALSE);
        max_index_per_axis[axis] = Position(function(x) x >= 1L, colmax, right=TRUE);
    }
    coords = boxcoords.from.bbox(min_index_per_axis, max_index_per_axis);
    bbox = list("from"=min_index_per_axis, "to"=max_index_per_axis, "edge_coords"=coords);
    if(apply) {
        return(vol.boundary.box.apply(volume, bbox));
    } else {
        return(bbox);
    }
}


#' @title Apply a boundary box to a volume, returning the inner volume part
#'
#' @param volume a 3D image volume
#'
#' @param bbox the boundary box
#'
#' @return a 3D image volume, the inner volume part, resulting from the application of the boundary box
#'
#' @export
vol.boundary.box.apply <- function(volume, bbox) {
    return(volume[bbox$from[1]:bbox$to[1], bbox$from[2]:bbox$to[2], bbox$from[3]:bbox$to[3]]);
}


#' @title Compute the coordinates of the 8 corners of a 3D box.
#'
#' @description Given the extreme values (min and max) along the 3 axes, compute the coordinates of the 8 corners of a 3D box.
#'
#' @param axes_min numerical vector of length 3, the min values of the 3 axes
#'
#' @param axes_max numerical vector of length 3, the max values of the 3 axes
#'
#' @return numerical matrix with 3 columns and 8 rows, the edge coordinates
#'
#' @keywords internal
boxcoords.from.bbox <- function(axes_min, axes_max) {
    xlen = axes_max[1] - axes_min[1];
    ylen = axes_max[2] - axes_min[2];
    zlen = axes_max[3] - axes_min[3];

    coords = matrix(rep(0.0, 8*3), ncol=3);
    coords[1,] = axes_min;
    coords[2,] = c(axes_min[1] + xlen, axes_min[2], axes_min[3]);
    coords[3,] = c(axes_min[1], axes_min[2] + ylen, axes_min[3]);
    coords[4,] = c(axes_min[1] + xlen, axes_min[2] + ylen, axes_min[3]);

    coords[5,] = c(axes_min[1], axes_min[2], axes_max[3]);
    coords[6,] = c(axes_min[1] + xlen, axes_min[2], axes_max[3]);
    coords[7,] = c(axes_min[1], axes_min[2] + ylen, axes_max[3]);
    coords[8,] = axes_max;
    return(coords);
}


#' @title Get indices of the axes defining the given plane.
#'
#' @description When using plane names, this function assumes that the volume is in the standard FreeSurfer orientation, as returned by reading a conformed volume with functions like \code{\link[fsbrain]{subject.volume}}.
#'
#' @param plane integer or string. If a string, one of "axial", "coronal", or "sagittal". If this is an integer vector of length 2 already, it is returned as given. If it is a single integer, it is interpreted as an axis index, and the plane orthogonal to the axis is returned. A warning on using the plane names: these only make sense if the volume is in the expected orientation, no checking whatsoever on this is performed.
#'
#' @return integer vector of length 2, the axes indices.
#'
#' @keywords internal
vol.plane.axes <- function(plane) {
    if(is.double(plane)) {
        plane = as.integer(plane);
    }
    if(is.integer(plane)) {
        if(length(plane) == 2L) {   # Already done. We do not check whether the contents makes sense in this case.
            return(plane);
        } else if(length(plane) == 1L) {
            # Treat it as an axis, and return the plane that is orthogonal to the axis
            if(plane == 1L) {
                return(c(1L, 2L));
            } else if(plane == 2L) {
                return(c(2L, 3L));
            } else if(plane == 3L) {
                return(c(3L, 1L));
            } else {
                stop("If plane is an integer (vector), the values must be in range 1..3");
            }
        } else {
            stop("If plane is an integer vector, it must have length 1 or 2.");
        }
    }

    # The plane must be a character string when we hit this.
    if(!(plane %in% c("axial", "coronal", "sagittal"))) {
        stop(sprintf("Parameter 'plane' must be on of c('axial', 'coronal', 'sagittal') but is '%s'.\n", plane));
    }
    if(plane == "sagittal") {
        return(c(2L, 3L));
    } else if(plane == "coronal") {
        return(c(1L, 2L))
    } else { # axial
        return(c(3L, 1L));
    }
}


#' @title Translate names and indices of planes.
#'
#' @description Translate names and indices of 3D image planes. The names only make sense if the data in the volume is in the default FreeSurfer conformed orientation.
#'
#' @param plane NULL, a plane index, or a plane name.
#'
#' @return if `plane` is NULL, all available planes and their indices as a named list. If `plane` is an integer (a plane index), its name. If `plane` is an characters string (a plane name), its index.
#'
#' @family volume utility
#'
#' @export
vol.planes <- function(plane=NULL) {
    planes = list("coronal"=1, "sagittal"=2, "axial"=3);
    if(is.null(plane)) {
        return(planes);
    } else {
        if(is.character(plane)) {
            return(planes[[plane]]);
        } else {
            return(names(planes)[plane]);
        }
    }
}


#' @title Turn volume into an ImageMagick image stack.
#'
#' @description Create an image from each slice along the axis, then stack those into an ImageMagick image stack.
#'
#' @param volume a 3D image volume. Can be numeric, or something that can be read directly by \code{magick::image_read} in 2D matrices (slices along the axis), e.g., a 3D array of color strings. If a 2D matrix is passed, the resulting stack will contain a single image.
#'
#' @param axis positive integer in range 1L..3L or an axis name, the axis to use.
#'
#' @param intensity_scale integer, value by which to scale the intensities in the volume to the range `[0, 1]`. Only used for numeric volumes. Set to NULL for data that can be read directly by \code{magick::image_read}, and to 1 for intensity data that requires no scaling. Defaults to 255, which is suitable for 8 bit image data.
#'
#' @return a vectorized ImageMagick image, containing one subimage per slice. This can be interpreted as an animation or whatever.
#'
#' @family volume utility
#'
#' @export
vol.imagestack <- function(volume, axis=1L, intensity_scale=255) {
    if(is.character(axis)) {
        axis = vol.planes(axis);
    }
    axis = as.integer(axis);
    if(axis < 1L | axis > 3L) {
        stop(sprintf("Axis must be integer with value 1, 2 or 3 but is %d.\n", axis));
    }

    if(length(dim(volume)) == 2) {
        if(is.null(intensity_scale) | !is.numeric(volume)) {
            return(magick::image_read(volume));
        } else {
            return(magick::image_read(grDevices::as.raster(volume / intensity_scale)));
        }
    } else if(length(dim(volume)) == 3) {
        if(is.null(intensity_scale) | !is.numeric(volume)) {
            image_list = apply(volume, axis, function(x){magick::image_read(x)});
        } else {
            image_list = apply(volume, axis, function(x){magick::image_read(grDevices::as.raster(x / intensity_scale))});
        }
        image_stack = Reduce(c, image_list);
        return(image_stack);
    } else{
        stop("The image must have exactly 2 or 3 dimensions.");
    }
}


#' @title Generate colors for a 3D volume, based on the activation data and a colormap.
#'
#' @description Applies the colormap function to the data, then sets the alpha value (transparency) to full in all areas without any activation. Feel free to clip data or whatever before passing it, so that all your no-activation data has the same value.
#'
#' @param volume a 3D array, the activation data (or p-values, effect sizes, or whatever)
#'
#' @param colormap_fn function, a colormap function
#'
#' @param no_act_source_value numerical scalar, the value from the data in 'volume' that means no activation. The output colors for this value will be set to `NA`. Set to NULL to not change anything.
#'
#' @return a 3D matrix of color strings, with the same dimensions as the input volume
#'
#' @importFrom squash makecmap blueorange cmap
#' @importFrom grDevices adjustcolor
#'
#' @family volume utility
#'
#' @export
vol.overlay.colors.from.activation <- function(volume, colormap_fn=squash::blueorange, no_act_source_value=0) {
    col = squash::cmap(volume, map = squash::makecmap(volume, colFn = colormap_fn));
    if( ! is.null(no_act_source_value)) {
        no_act_indices = which(volume == no_act_source_value, arr.ind = TRUE);
        col[no_act_indices] = NA;
    }
    return(col);
}


#' @title Draw a lightbox view from volume slices.
#'
#' @description A lightbox is a single image that holds a set of subimages, arranged in a grid. The images can have a small border or spacing between them. Consecutive subimages will be appear the same row of the grid.
#'
#' @param volume 3D array, can be numeric (gray-scale intensity values) or color strings. If numeric, the intensity values must be in range `[0, 1]`.
#'
#' @param slices slice index definition. If a vector of integers, interpreted as slice indices. If a single negative interger `-n`, interpreted as every `nth` slice, starting at slice 1. The character string 'all' or the value `NULL` will be interpreted as *all slices*.
#'
#' @param axis positive integer in range 1L..3L, the axis to use.
#'
#' @param scale integer >= 1, the factor by which to upscale the slice images in both dimensions using nearest-neighbor interpolation. This keeps the MRI pixels sharp (no smoothing) while making the images larger. Defaults to 1L (no upscaling).
#'
#' @param per_row positive integer, the number of subimages per row in the output image. If `NULL`, automatically computed from the number of slices and the `per_col` parameter.
#'
#' @param per_col positive integer, the number of subimages per column in the output image. If `NULL`, automatically computed from the number of slices and the `per_row` parameter.
#'
#' @param border_geometry string, a geometry string passed to \code{magick::image_border} to define the borders to add to each image tile. The default value adds 5 pixels, both horizontally and vertically.
#'
#' @param background_color string, a valid ImageMagick color string such as "white" or "#000080". The color to use when extending images (e.g., when creating the border). Defaults to black.
#'
#' @param arrange_single_image logical, whether to apply the given arrangement (from parameters `per_row` and `per_column`) even if a single slice (a 2D image) is passed as `volume`. Defaults to FALSE, which prevents that background tiles are added to fill the row up to `per_row` images. This also prevents the border from getting added to a single image, so all you see is the raw image. Set to `TRUE` if you want to arrange even a single image in a row with a border.
#'
#' @description If overlay_colors are given, the volume will be used as the background, and it will only be visible where overlay_colors has transparency.
#'
#' @return a magick image instance
#'
#' @family volume visualization
#'
#' @note You should, in most cases, not call this function directly. Use \code{\link[fsbrain]{volvis.lb}} instead, which has a more intuitive interface.
#'
#' @seealso \code{\link[fsbrain]{volvis.lb}}
#'
#' @export
volvis.lightbox <- function(volume, slices=-5, axis=1L, per_row=5L, per_col=NULL, border_geometry="5x5", background_color = "#000000", arrange_single_image=FALSE, scale=1L) {

    skip_border = FALSE;
    if(length(dim(volume)) == 2) {
        just_a_slice = volume;
        message("Inflating single 2D slice to volume, new axis added at position 1.");
        volume = array(just_a_slice, dim = c(1, dim(just_a_slice)[1], dim(just_a_slice)[2]));
        if(!arrange_single_image) {
            per_row=1L;
            skip_border = TRUE;
        }
    }

    if(is.character(axis)) {
        axis = vol.planes(axis);
    }
    axis = as.integer(axis);
    if(axis < 1L | axis > 3L) {
        stop(sprintf("Axis must be integer with value 1, 2 or 3 but is %d.\n", axis));
    }

    scale = as.numeric(scale);
    if(length(scale) != 1L || is.na(scale) || scale < 1 || scale != round(scale)) {
        stop("Parameter 'scale' must be a single integer >= 1.");
    }
    scale = as.integer(scale);

    if(is.numeric(volume)) {
        volume = vol.intensity.to.color(volume);
    }

    if(length(dim(volume)) != 3) {
        stop("Volume must have exactly 3 dimensions.");
    }

    # Compute the slice indices from the slice definition
    slice_indices = get.slice.indices(dim(volume), axis, slices);

    # Get the subset of requested slices as a 3D image
    img_slices = vol.slice(volume, slice_index=slice_indices, axis=axis);

    # Transform the slices into an ImageMagick stack of 2D images
    images = vol.imagestack(img_slices, axis=axis);

    # Optionally upscale the slices with nearest-neighbor (keeps pixels sharp).
    if(scale > 1L) {
        img_info_all = magick::image_info(images);
        images = magick::image_resize(images, magick::geometry_size_pixels(img_info_all$width[1] * scale, img_info_all$height[1] * scale), filter = "point");
    }

    # Add tiny border
    if(!skip_border) {
        if(!((is.null(border_geometry) | is.null(background_color)))) {
            #message("Applying border to individual images.");
            images = magick::image_border(images, background_color, border_geometry);
        }
    }

    # Arrange the stack of images
    #merged_img = magick::image_append(images);
    merged_img = magick.grid(images, per_row=per_row, per_col=per_col, background_color = background_color);

    return(merged_img);
}


#' @title Arrange a multi-frame ImageMagick image into a grid.
#'
#' @description Arrange all subimages of the given ImageMagick image into a single 2D image, that contains the subimages arranged in a grid-like structure. Consecutive subimages will be appear the same row.
#'
#' @param magickimage an ImageMagick image
#'
#' @param per_row positive integer, the number of subimages per row in the output image. If `NULL`, automatically computed from the number of slices and the `per_col` parameter.
#'
#' @param per_col positive integer, the number of subimages per column in the output image. If `NULL`, automatically computed from the number of slices and the `per_row` parameter.
#'
#' @param background_color string, a valid ImageMagick color string such as "white" or "#000080". The color to use when extending images (e.g., when creating the border). Defaults to black.
#'
#' @keywords internal
magick.grid <- function(magickimage, per_row=5L, per_col=NULL, background_color = "#000000") {
    images = magickimage;
    num_subimages = length(images);
    if(is.null(per_row) & is.null(per_col)) {
        # If both are none, just return one horizontal strip of images.
        return(magick::image_append(images));
    }

    if(!(is.null(per_row) | is.null(per_col))) {
        # Both are given, let's see whether they make sense.
        expected_per_row = ceiling(num_subimages / as.double(per_col));
        if(per_row != expected_per_row) {
            warning(sprintf("Changing 'per_row' value from %d to %d based on %d subimages and %d per column.\n", per_row, expected_per_row, num_subimages, per_col));
            per_row = expected_per_row;
        }
    } else {
        if(is.null(per_row)) {
            per_row = ceiling(num_subimages / as.double(per_col));
        } else {
            per_col = ceiling(num_subimages / as.double(per_row));
        }
    }

    num_rows = per_col;
    num_columns = per_row;

    #message(sprintf("Distributing %d images over %d rows and %d colums.\n", num_subimages, num_rows, num_columns));

    img_start_indices_each_row = seq.int(from = 1L, to = num_subimages, by=per_row);
    image_data = NULL;
    for(row_idx in seq_len(length(img_start_indices_each_row))) {
        row_start_index = img_start_indices_each_row[[row_idx]];
        row_end_index = row_start_index + (per_row -1);
        num_this_row = per_row;
        if(row_end_index > num_subimages) {  # Can only happen in last row.
            out_of_bounds_row_end_index = row_end_index;
            row_end_index = row_start_index + (num_subimages %% per_row) -1;
            num_this_row = (row_end_index - row_start_index) + 1;
            #message(sprintf("At row %d of %d, handling the %d subimages %d to %d of %d.\n", row_idx, num_rows, num_this_row, row_start_index, row_end_index, num_subimages));

            if(!is.null(background_color)) {
                # Fill the rest of the background with the background color. Otherwise, we will get a white background for the remainder of the last row.
                num_missing = out_of_bounds_row_end_index - row_end_index;
                background_tile = magick::image_blank(magick::image_info(images[1])$width, magick::image_info(images[1])$height, background_color);
                for(tile_idx in 1:num_missing) {
                    images = c(images, background_tile);
                }
                #message(sprintf("-At row %d of %d: row had %d images, missed %d background tiles to a full row of %d images, added them.\n", row_idx, num_rows, num_this_row, num_missing, per_row));
                #message(sprintf("-After adding the  %d missing background tiles, changes row_end_index from %d to %d.\n", num_missing, row_end_index, out_of_bounds_row_end_index));
                row_end_index = out_of_bounds_row_end_index; # It's not out of bounds anymore since we added the background tiles.

            }
        }

        if(is.null(image_data)) {
            image_data = magick::image_append(images[row_start_index:row_end_index], stack=FALSE);
        } else {
            image_row = magick::image_append(images[row_start_index:row_end_index], stack=FALSE);
            image_data = magick::image_append(c(image_data, image_row), stack = TRUE);
        }
    }
    return(image_data);

}


#' @title Compute slice indices from slice definition.
#'
#' @param voldim integer vector, the dimension of the volume
#'
#' @param axis integer, the axis
#'
#' @param slices slice index definition. If a vector of integers, interpreted as slice indices. If a single negative interger `-n`, interpreted as every `nth` slice, starting at slice 1. The character string 'all' or the value `NULL` will be interpreted as *all slices*.
#'
#' @return integer vector, the computed slice indices. They are guaranteed to be valid indices into the volume.
#'
#' @keywords internal
get.slice.indices <- function(voldim, axis, slices) {
    if(is.character(axis)) {
        axis = vol.planes(axis);
    }
    axis = as.integer(axis);
    num_slices_in_volume = voldim[axis];    # along the requested axis
    if(is.numeric(slices)) {
        if(length(slices) == 1L && slices[1L] < 0L) {
            # every nth slice
            return(seq.int(from=1L, to=num_slices_in_volume, by=abs(slices)))
        } else {
            if(any(slices > num_slices_in_volume)) {
                stop(sprintf("The %d requested slice indices include %d which are out of bounds, volume has %d slices along axis %d.\n", length(slices), length(which(slices > num_slices_in_volume)), num_slices_in_volume, axis));
            }
            if(any(slices < 0)) {
                stop(sprintf("The %d requested slice indices include %d negative ones which are invalid.\n", length(slices), length(which(slices < 0))));
            }
            return(slices);
        }
    } else if(is.character(slices)) {
        if(slices == "all") {
            return(seq_len(num_slices_in_volume));
        } else {
            stop(sprintf("Slice definition character string '%s' not recognized.", slices));
        }
    } else if (is.null(slices)) {
        return(seq_len(num_slices_in_volume));
    } else {
        stop("Invalid slice definition.");
    }
}


#' @title Extract subset from a volume by value.
#'
#' @description Extract subset from a volume by value, set all other voxel values to `NA`. Typically used to extract a brain structure (like corpus callosum) from a volume segmentation (like the `mri/aseg.mgz` file of a subject). You should consider passing the volume and the include values as integers.
#'
#' @param volume numeric 3D array
#'
#' @param include_values numerical vector, the intensity values which qualify a voxel to be part of the result (without being set to NA)
#'
#' @return numerical array with same dimensions as the input volume. All values which are not part of `include_values` replaced with `NA`.
#'
#' @export
vol.mask.from.segmentation <- function(volume, include_values) {
    if(length(dim(volume)) < 3) {
        stop("Volume must have 3 or more dimensions.");
    }

    mask = array(rep(NA, prod(dim(volume))), dim(volume));

    for(inc_val in include_values) {
        act_indices = which(volume == inc_val, arr.ind = TRUE);
        mask[act_indices] = volume[act_indices];
    }
    return(mask);
}


#' @title Merge background volume and overlay to new colors.
#'
#' @param volume 3D array, can be numeric (gray-scale intensity values) or color strings. If numeric, the intensity values must be in range `[0, 1]`.
#'
#' @param overlay_colors 3D array of color strings, values which are not part of the overlay (and should display background in the result) must have `NA` instead of a color string. Must have same dimensions as the `volume`.
#'
#' @param bbox_threshold numerical, the threshold intensity used to separate background and foreground. All voxels with intensity values greater than this value in the background `volume` will be considered `foreground` voxels. Background-only slices at the borders of the volume will be discarded (in the merged, final image). Pass `NULL` to use the full image without applying any bounding box.
#'
#' @param forced_overlay_color NULL or an rgb color string, like '#FF0000' for red. If NULL, the activation colors will be used as foreground colors. Otherwise, the given color will be for all foreground vertices.
#'
#' @return 3D array of color strings, the merged colors
#'
#' @family volume utility
#'
#' @importFrom grDevices rgb
#' @export
vol.merge <- function(volume, overlay_colors, bbox_threshold=0L, forced_overlay_color=NULL) {
    if(length(dim(volume)) != 3) {
        stop("Volume must have exactly 3 dimensions.");
    }

    if(!(all.equal(dim(volume), dim(overlay_colors)))) {
        stop("If 'overlay_colors' are given, they must have the same dimension as the 'volume'. Hint: use RGB color strings.");
    }

    if(!is.null(bbox_threshold)) {
        bbox = vol.boundary.box(volume, threshold=bbox_threshold);
        volume = volume[bbox$from[1]:bbox$to[1], bbox$from[2]:bbox$to[2], bbox$from[3]:bbox$to[3]];
        overlay_colors = overlay_colors[bbox$from[1]:bbox$to[1], bbox$from[2]:bbox$to[2], bbox$from[3]:bbox$to[3]];
    }

    # Compute background volume color strings from intensity values if needed.
    if(is.numeric(volume)) {
        merged = vol.intensity.to.color(volume);
    } else {
        merged = volume;
    }

    # Same for the overlay
    if(is.numeric(overlay_colors)) {
        message("Overlay is numerical and will be gray-scale. It may be hard to discern from the background volume.");
        overlay_colors = vol.intensity.to.color(overlay_colors);
    }

    if(!(all.equal(dim(volume), dim(overlay_colors)))) {
        stop("Bug: dimensions of 'volume' do not match 'overlay_colors' anymore");
    }

    # Copy background colors into NA voxels of the activation.
    overlay_idc = which(!is.na(overlay_colors), arr.ind = TRUE);
    #message(sprintf("Setting %d activated vertices to colors.\n", nrow(overlay_idc)));

    if(is.null(forced_overlay_color)) {
        merged[overlay_idc] = overlay_colors[overlay_idc];
    } else {
        if(!is.character(forced_overlay_color)) {
            stop("Parameter 'forced_overlay_color' must be a color string (like '#FF0000') or NULL.");
        } else {
            merged[overlay_idc] = forced_overlay_color;
        }
    }

    return(merged);
}


#' @title Convert integer intensity image to RGB color string form.
#'
#' @description Convert a gray-scale image defined by intensity values in range `[0, 1]` to an image with identical dimensions that contains an R color string (like `#222222`) at each position. The color strings are computed from the intensities, by taking the intensity value as the value for all three RGB channels. I.e., the output is still gray-scale, but defined in RGB space. To make it clear, this function does **not** apply a colormap. It only changes the representation of the data, not the resulting colors.
#'
#' @param volume numeric array, typically a 3D image with intensities in range `[0, 1]`. This function now also supports numeric matrices (2D images, slices) and numeric vectors (1D).
#'
#' @param scale numeric or character string, a scaling to apply to the values. Defaults to NULL, which means *no scaling* and requires the values in `volume` to be in rage `[0, 1]`. You can pass a number like 255 or the string 'normalize' to scale based on the data. You can pass the string 'normalize_if_needed' to scale only if the data is *outside* the range `[0, 1]`, so that data in range `[0.3, 0.5]` would **not** be rescaled to `[0, 1]`.
#'
#' @return array (or matrix, or vector) of RGB color strings. All of them will represent gray values.
#'
#' @examples
#'    vol.intensity.to.color(c(0.0, 0.5, 1.0));
#'    # output: "#000000" "#808080" "#FFFFFF"
#'    vol.intensity.to.color(c(20, 186, 240), scale="normalize");
#'    vol.intensity.to.color(c(20, 186, 240), scale=255);
#'    vol.intensity.to.color(c(0.0, 0.5, 0.8), scale="normalize");
#'    vol.intensity.to.color(c(0.0, 0.5, 0.8), scale="normalize_if_needed");
#'
#' @importFrom grDevices rgb
#' @export
vol.intensity.to.color <- function(volume, scale=NULL) {
    if(is.numeric(volume)) {

        if(!is.null(scale)) {
            rng = range(volume);
            if(scale == 'normalize') {
                volume = normalize(volume);
            } else if(scale == 'normalize_if_needed') {
                if(rng[1] < 0.0 | rng[2] > 1.0) {
                    volume = normalize(volume);
                }
            } else {
                if(is.numeric(scale) & length(scale) == 1) {
                    volume = volume / scale;
                } else {
                    stop("Parameter 'scale' must be exactly the character string 'normalize', 'normalize_if_needed', or a numeric scalar.");
                }
            }
        }

        rng = range(volume);
        if(rng[1] < 0.0 | rng[2] > 1.0) {
            warning(sprintf("Intensity values of volume are in range range [%.2f, %.2f], please scale the intensity values to range [0, 1] before passing them to this function or use 'scale' parameter.\n", rng[1], rng[2]));
        }

        num_dims = length(dim(volume));
        if(num_dims == 3L) {
            return(array(grDevices::rgb(volume, volume, volume), dim(volume))); # try magick::image_read(vol.slice(return_value)) or volvis.lightbox(return_value)
        } else if (num_dims == 2L) {
            return(matrix(grDevices::rgb(as.vector(volume), as.vector(volume), as.vector(volume)), nrow=nrow(volume)));
        } else if (is.vector(volume)) {
            return(grDevices::rgb(volume, volume, volume));
        } else {
            stop("Parameter 'volume' must have 1, 2, or 3 dimensions.");
        }
    } else {
        stop("Parameter 'volume' must be numeric.");
    }
}


#' @title Normalize data.
#'
#' @description Scales data to the range `[0, 1]` based on min and max values.
#'
#' @param x the data
#'
#' @return the scaled data
#'
#' @keywords internal
normalize <- function(x) {
    return((x- min(x)) /(max(x)-min(x)));
}


#' @title Return triangles for a 3D cube or cuboid.
#'
#' @description Each row of the returned matrix encodes a point (the x, y, and z coordinates), and 3 consecutive rows encode a triangle. Obvisouly, a point will occur several times (as part of several triangles). The result can be passed to \code{\link[rgl]{triangles3d}} to render a 3D box. The defaults for the parameters will create a cube with edge length 1 centered at (0, 0, 0).
#'
#' @param xmin numeric, minimal x coordinate
#'
#' @param xmax numeric, maximal x coordinate
#'
#' @param ymin numeric, minimal y coordinate
#'
#' @param ymax numeric, maximal y coordinate
#'
#' @param zmin numeric, minimal z coordinate
#'
#' @param zmax numeric, maximal z coordinate
#'
#' @param center numeric vector of length 3 or NULL, coordinates where to center a cube with the edge length defined in parameter `edge_length`. If this is not `NULL`, the parameters `xmin`, `xmax`, ... will be ignored, and their values will be computed for a cube based on the `center` and `edge_length`. Note that you can only create cubes using `center` and `edge_length`, while the min/max methods allows the construction of cuboids.
#'
#' @param edge_length numeric, the edge length of the cube. Only used if parameter `center` is used, ignored otherwise.
#'
#' @return numerical matrix with 36 rows and 3 columns, the 3D coordinates. Each row encodes a point (the x, y, and z coordinates), and 3 consecutive rows encode a triangle.
#'
#' @examples
#'    # Create a cube with edge length 2, centered at (3,4,5):
#'    cube_coords = cube3D.tris(center=c(3,4,5), edge_length=2.0);
#'    # Create the same cube using the min/max method:
#'    cube_coords = cube3D.tris(xmin=2, xmax=4, ymin=3, ymax=5, zmin=4, zmax=6);
#'    # Create a cuboid:
#'    cuboid_coords = cube3D.tris(xmin=2, xmax=4, ymin=3, ymax=9, zmin=4, zmax=5);
#'    # To render the cuboid:
#'    #rgl::triangles3d(cuboid_coords, col="red");
#'
#' @export
cube3D.tris <- function(xmin=-0.5, xmax=0.5, ymin=-0.5, ymax=0.5, zmin=-0.5, zmax=0.5, center=NULL, edge_length=1.0) {

    if(!is.null(center)) {
        if(is.numeric(center) & length(center) == 3) {
            hel = edge_length / 2.0;   # half edge length
            xmin = center[1] - hel;
            xmax = center[1] + hel;
            ymin = center[2] - hel;
            ymax = center[2] + hel;
            zmin = center[3] - hel;
            zmax = center[3] + hel;

        } else {
            stop("If given, 'center' must be a numeric vector of length 3: the x, y, and z coordinates of the cube center.");
        }
    }

    tris_cube = matrix(c(xmin, ymin, zmin, # the 2 front tris start
                         xmin, ymin, zmax,
                         xmax, ymin, zmin,
                         # tris 2
                         xmax, ymin, zmin,
                         xmin, ymin, zmax,
                         xmax, ymin, zmax,
                         # tris 3: the back starts (shifted by +1.0 on y axis compared to front)
                         xmin, ymax, zmin,
                         xmin, ymax, zmax,
                         xmax, ymax, zmin,
                         # tris 4
                         xmax, ymax, zmin,
                         xmin, ymax, zmax,
                         xmax, ymax, zmax,
                         # tris 5: the left side starts
                         xmin, ymax, zmin,
                         xmin, ymax, zmax,
                         xmin, ymin, zmin,
                         # tris 6
                         xmin, ymin, zmin,
                         xmin, ymax, zmax,
                         xmin, ymin, zmax,
                         # tris 7: the right side starts (shifted by +1.0 on x axis compared to left side)
                         xmax, ymax, zmin,
                         xmax, ymax, zmax,
                         xmax, ymin, zmin,
                         # tris 8
                         xmax, ymin, zmin,
                         xmax, ymax, zmax,
                         xmax, ymin, zmax,
                         # tris 9: the bottom starts
                         xmin, ymin, zmin,
                         xmin, ymax, zmin,
                         xmax, ymin, zmin,
                         # tris 10
                         xmax, ymin, zmin,
                         xmin, ymax, zmin,
                         xmax, ymax, zmin,
                         # tris 11: the top starts (shifted by +1.0 on z axis compared to the bottom)
                         xmin, ymin, zmax,
                         xmin, ymax, zmax,
                         xmax, ymin, zmax,
                         # tris 12
                         xmax, ymin, zmax,
                         xmin, ymax, zmax,
                         xmax, ymax, zmax
    ), ncol=3, byrow = TRUE);
    return(tris_cube);
}



#' @title Vectorized version of cube3D.tris
#'
#' @param centers numerical matrix with 3 columns. Each column represents the x, y, z coordinates of a center at which to create a cube.
#'
#' @param edge_length numerical vector or scalar, the edge length. Must have length 1 (same edge length for all cubes), or the length must be identical to the number of rows in parameter `centers`.
#'
#' @return matrix of triangle coordinates, see \code{\link[fsbrain]{cube3D.tris}}.
#'
#' @examples
#'    # Plot a 3D cloud of 20000 voxels:
#'    centers = matrix(rnorm(20000*3)*100, ncol=3);
#'    rgl::triangles3d(cubes3D.tris(centers));
#'
#' @export
cubes3D.tris <- function(centers, edge_length=1) {
    if(is.vector(centers)) {
        centers = matrix(centers, ncol = 3);
    }

    num_tris = nrow(centers);

    if(length(edge_length) == 1) {
        edge_length = rep(edge_length, num_tris);
    }

    points_per_rect = 36L;    # Each 3D rectangle consists of 12 triangles, each of which has 3 points.

    tris = matrix(rep(NA, num_tris*points_per_rect*3), ncol=3);
    for(row_idx in seq_len(nrow(centers))) {
        start_idx = (row_idx - 1) * points_per_rect + 1;
        end_idx = start_idx + points_per_rect - 1;
        tris[start_idx:end_idx,] = cube3D.tris(center=centers[row_idx,], edge_length = edge_length[row_idx]);
    }
    return(tris);
}


#' @title Compute voxel colors based on colortable.
#'
#' @description Use the intensity values of the voxels in volume and lookup the respective colors in a colortable.
#'
#' @param volume numeric 3D array, the values should be integers present in the `struct_index` column of the colortable. All other values will be assigned `NA` as a color.
#'
#' @param colortable a colortable, as returned by \code{\link[freesurferformats]{read.fs.colortable}}, or a character string representing a path to a colortable file.
#'
#' @param ignored_struct_indices integer vector, `struct_index` entries in the colortable that should be ignored
#'
#' @param ignored_struct_names vector of character strings, `struct_name` entries in the colortable that should be ignored. Can be combined with `ignored_struct_indices`.
#'
#' @return character string 3D array, the colors. Voxels in the volume which were not matched by the colortable are set to `NA` in it.
#'
#' @export
vol.overlay.colors.from.colortable <- function(volume, colortable, ignored_struct_indices=c(), ignored_struct_names=c('unknown', 'Unknown')) {
    if(length(dim(volume)) != 3) {
        stop("Volume must have exactly 3 dimensions.");
    }

    overlay_colors = array(rep(NA, length(volume)), dim(volume));

    if(is.character(colortable) && length(colortable) == 1L) {
        if(file.exists(colortable)) {
            colortable = freesurferformats::read.fs.colortable(colortable);
        } else {
            stop(sprintf("Received character string  '%s' for parameter 'colortable', but no such file exists.", colortable));
        }
    }

    occuring_color_rows = which(colortable$struct_index %in% volume);

    for(row_idx in occuring_color_rows) {
        ct_structure = colortable[row_idx, ];
        if(ct_structure$struct_index %in% ignored_struct_indices | ct_structure$struct_name %in% ignored_struct_names) {
            #cat(sprintf("Skipping colortable entry with struct_index %d, named '%s'.\n", ct_structure$struct_index, ct_structure$struct_name));
            next;
        }
        struct_voxels = which(volume==ct_structure$struct_index, arr.ind = TRUE);
        voxel_color = grDevices::rgb(ct_structure$r/255., ct_structure$g/255., ct_structure$b/255.);
        if(nrow(struct_voxels) > 0) {
            #cat(sprintf("Assigning %d voxels to color %s from colortable entry with struct_index %d, named '%s'.\n", nrow(struct_voxels), voxel_color, ct_structure$struct_index, ct_structure$struct_name));
            overlay_colors[struct_voxels] = voxel_color;
        }
    }
    return(overlay_colors);
}


#' @title Try to extract a 3D volume from the input argument.
#'
#' @description Check whether it already is such an array, whether it is a filename that can be loaded with \code{freesurferformats::read.fs.volume} into such an array, etc.
#'
#' @param stats a 3D array, 4D array or a string that can be treated as a filename to a volume image containing such an array.
#'
#' @param silent logical, whether to suppress file reading messages.
#'
#' @note This function stops with an error if the input cannot be returned as a 3D array.
#'
#' @return the obtained 3D array
#'
#' @keywords internal
extract.volume.3D <- function(stats, silent = getOption("fsbrain.silent", default = FALSE)) {
    if(! (is.array(stats) && length(dim(stats)) == 3L)) {
        if(is.character(stats) && is.vector(stats) && length(stats) == 1L) {
            if(! file.exists(stats)) {
                stop(sprintf("Received single character vector '%s', but no such file exists.", stats));
            }
            if(!silent) {
                message(sprintf("Reading volume file '%s'...\n", stats));
            }

            stats = freesurferformats::read.fs.volume(stats, drop_empty_dims = TRUE);
            if((is.array(stats) && length(dim(stats)) == 4L)) {
                message("Using first frame of loaded 4D image (assuming time is the 4th dimension in the volume).");
                stats = stats[,,,1]; # select first frame
            }
            if(! ((is.array(stats) && length(dim(stats)) == 3L))) {
                stop(sprintf("Volume data loaded from file '%s' does not have 3 non-empty dimensions, not a 3D volume.", stats));
            }
        } else {
            stop("Parameter 'stats' must be a 3D array.");
        }
    }
    if((is.array(stats) && length(dim(stats)) == 4L)) {
        message("Using first frame of input 4D image (assuming time is the 4th dimension in the volume).");
        stats = stats[,,,1]; # select first frame
    }
    if(! (is.array(stats) && length(dim(stats)) == 3L)) {
        stop("Could not extract 3D volume from input.");
    } else {
        return(stats);
    }
}


#' @title Show continuous 3D voxel/volume data as a lightbox, optionally with a background brain volume and colormap.
#'
#' @description This function is the main way to visualize 3D volume images that contain raw MRI scans or statistical results.
#'
#' @param volume numerical 3D array of per-voxel data, typically activation data, a raw MRI image, or a segmentation to show. Can also be a filename if the file can be loaded as such a volume with \code{\link[freesurferformats]{read.fs.volume}}.
#'
#' @param background numerical 3D array or 3D array of color strings, the background volume. Typically a raw brain volume. Dimensions and space must match those of the 'volume' for an array. Can also be a single file name as a character string. Can also be a single color name, like '#FEFEFE' but the string then must start with '#' (color names like 'red' are not allowed, they would be treated as file names). If a color string, be sure to use the \code{...} parameter to set the same color as \code{background_color} for the tiles.
#'
#' @param colFn a colormap function, passed to \code{vol.overlay.colors.from.activation} and used as colormap for the 'volume' data. Pass NULL to derive gray-scale values from the raw data (only recommended with single-color backgrounds). Note that the colormap is not used for the the background data (if any), which will be shown in grayscale (unless it is a 3D array of color strings).
#'
#' @param colortable optional, only makes sense for categorical 'volume' data like segmentations. If not NULL, a colortable as returned by \code{\link[freesurferformats]{read.fs.colortable}}, or a character string representing a path to a colortable file (like \code{"FREESURFER_HOME/FreeSurferColorLUT.txt"]}).
#'
#' @param no_act_source_value numerical value, passed to \code{vol.overlay.colors.from.activation}. Specifies the value which is treated as transparent in the 'volume' parameter data (where you will see the background). If you need more control, e.g., you want to treat one or morge ranges of values as NA, you should load the 'volume' data first, modify it as needed, as pass it to this function afterwards. Set this parameter to \code{NULL} to disable it. Only for 'colFn', ignored if a 'colortable' is used.
#'
#' @param bbox_threshold numerical scalar, passed on to \code{vol.merge}. If set, voxels with intensities smaller than this threshold will be dropped at the outside of the image. If \code{bbox_of_volume} parameter is \code{TRUE} (the default), this threshold applies to the 'volume', otherwise to the 'background'. Set to \code{NULL} to disable bounding box and show the full image.
#'
#' @param bbox_of_volume logical, whether the bounding box is computed on the volume (foreground), which typically is what you want. Leave alone if in doubt.
#'
#' @param ... extra parameters to be passed to \code{\link[fsbrain]{volvis.lightbox}}, can be used to select specific slices, set the \code{backgroud_color} for the border between and around the image tiles, etc.
#'
#' @note This function should be preferred over manually calling \code{\link[fsbrain]{volvis.lightbox}}.
#'
#' @examples
#' \dontrun{
#' volume = subject.volume(subjects_dir, subject_id, 'brain');
#' volvis.lb(volume);
#' volvis.lb("~/study1/subject1/mri/brain.mgz");
#' volvis.lb("~/study1/subject1/mri/brain.mgz", bbox_threshold = 1L);
#' volvis.lb("~/study1/subject1/mri/brain.mgz", background = "~/data/study1/subject1/mri/T1.mgz");
#' volvis.lb("~/study1/subject1/mri/brain.mgz", background = "#FEFEFE", background_color="#FEFEFE");
#' ct = file.path(find.freesurferhome(mustWork = T), "FreeSurferColorLUT.txt"); # ct = "color table"
#' volvis.lb("~/study1/subject1/mri/aseg.mgz", background="~/study1/subject1/mri/T1.mgz",
#'  colortable = ct, colFn=NULL, axis=2L);
#' volvis.lb("~/study1/subject1/mri/aseg.mgz", background = "~/study1/subject1/mri/T1.mgz",
#'  colortable = ct, colFn=NULL, bbox_threshold = 0);
#' }
#'
#' @family volume visualization
#'
#' @export
volvis.lb <- function(volume, background=NULL, colFn=viridis::viridis, colortable=NULL, no_act_source_value = 0, bbox_threshold = NULL, bbox_of_volume = TRUE, ...) {


    volume = extract.volume.3D(volume);

    #volume = clip.data(volume);
    #cat(sprintf("volume dim: %s\n", paste(dim(volume), collapse = ", ")));

    bbox = NULL;
    if(!is.null(bbox_threshold)) {

        if(bbox_of_volume) {
            bbox_of = volume;
        } else {
            if(!(is.array(background) && is.numeric(background))) {
                if(is.vector(background) && is.character(background) && length(background) == 1L && (!startsWith(background, "#"))) {
                    background = extract.volume.3D(background);
                } else {
                    stop("Setting parameter 'bbox_of_volume' to FALSE is only possible if parameter 'background' is a numerical 3D array.");
                }
            }
            bbox_of = background;
        }
        bbox = vol.boundary.box(bbox_of, threshold=bbox_threshold);
        volume = volume[bbox$from[1]:bbox$to[1], bbox$from[2]:bbox$to[2], bbox$from[3]:bbox$to[3]];
    }

    if(! is.null(colFn)) {
        if(!is.null(colortable)) {
            stop("Both parameters 'colFn' and 'colortable' are given, only one supported. Set 'colFn' to NULL to use a colortable.");
        }
        overlay_colors = vol.overlay.colors.from.activation(volume, colormap_fn = colFn, no_act_source_value = no_act_source_value);
    } else {
        if(! is.null(colortable)) {
            overlay_colors = vol.overlay.colors.from.colortable(volume, colortable);
        } else {
            overlay_colors = scale01(volume);
        }
    }

    if(! is.null(background)) {
        if(is.vector(background) && is.character(background) && length(background) == 1L) {
            if(startsWith(background, "#")) {
                bg_col = background;
                background = array(rep(bg_col, sum(dim(volume))), dim = dim(volume));
            }
        }

        background = extract.volume.3D(background);

        # Apply bounding box to loaded background file as well.
        if(!is.null(bbox_threshold)) {
            background = background[bbox$from[1]:bbox$to[1], bbox$from[2]:bbox$to[2], bbox$from[3]:bbox$to[3]];
        }


        if(length(dim(volume)) != length(dim(background))) {
            stop(sprintf("The dimensions of the volume and background parameters must be identical, but they are %s versus %s.", paste(dim(volume), collapse = ", "), paste(dim(background), collapse = ", ")));
        }
        for(dim_idx in seq_along(dim(volume))) {
            if(dim(volume)[dim_idx] != dim(background)[dim_idx]) {
                stop(sprintf("The dimensions of the volume and background parameters must be identical, but they are %s versus %s.", paste(dim(volume), collapse = ", "), paste(dim(background), collapse = ", ")));
            }
        }
        if(is.numeric(background)) {
            background = scale01(background);
        }
    } else {
        # use black background
        bg_col = "#000000"; # should get this from function parameters to allow customization.
        background = array(rep(bg_col, sum(dim(volume))), dim = dim(volume));
        # No need to apply bounding box, background dimensions are based on already bounded volume value.
    }


    color_vol = vol.merge(background, overlay_colors, bbox_threshold = NULL);
    return(volvis.lightbox(color_vol, ...));
}


# ══════════════════════════════════════════════════════════════════════════════════════
# Surface contour overlay on 2D volume slices (lightbox)
# ══════════════════════════════════════════════════════════════════════════════════════


#' @title Transform surface vertices from surface RAS to 0-based volume CRS space.
#'
#' @description Applies the inverse of the FreeSurfer vox2ras_tkr matrix to all surface vertices, converting them from surface RAS coordinates to 0-based CRS (column, row, slice) indices. The resulting CRS coordinates can be used to directly index into a brain volume array (after adding 1 for R's 1-based indexing).
#'
#' @param surface an `fs.surface` instance, as returned by \code{\link[freesurferformats]{read.fs.surface}}.
#'
#' @return a modified copy of the input surface with `vertices` in 0-based CRS space.
#'
#' @keywords internal
mesh.ras2crs <- function(surface) {
    ras2vox <- solve(vox2ras_tkr());
    verts_ras <- surface$vertices;
    verts_homog <- cbind(verts_ras, 1);
    verts_crs <- t(ras2vox %*% t(verts_homog))[, 1:3, drop = FALSE];
    surface$vertices <- verts_crs;
    return(surface);
}


#' @title Compute intersection of a triangular surface mesh with an axis-aligned plane.
#'
#' @description For a given axis and coordinate, finds all triangles that straddle the plane and computes the line segment where each crossing triangle intersects the plane. The input surface must have vertices in the same coordinate space as the slice plane definition (typically 0-based CRS).
#'
#' @param surface an `fs.surface` instance with vertices in CRS space (0-based), as produced by \code{\link[fsbrain]{mesh.ras2crs}}.
#'
#' @param axis integer, 1, 2, or 3. The axis perpendicular to the slice plane (1=sagittal, 2=coronal, 3=axial in CRS convention).
#'
#' @param slice_crs_coord numeric scalar, the coordinate along `axis` where the slice plane is located. In CRS space (0-based).
#'
#' @return a list of 2×2 numeric matrices. Each matrix represents one line segment with rows `[start, end]` and columns `[axis_val_1, axis_val_2]` giving the coordinates in the two axes orthogonal to the slice axis, in 0-based CRS units. Returns an empty list if no triangles intersect the plane.
#'
#' @keywords internal
mesh.slice.intersection <- function(surface, axis, slice_crs_coord) {
    verts <- surface$vertices;
    faces <- surface$faces;
    ax_val <- verts[, axis];

    # Signed distance from plane for all 3 triangle vertices (fully vectorized)
    d1 <- ax_val[faces[, 1]] - slice_crs_coord;
    d2 <- ax_val[faces[, 2]] - slice_crs_coord;
    d3 <- ax_val[faces[, 3]] - slice_crs_coord;

    # Which triangles cross the plane?
    crosses <- (d1 * d2 < 0) | (d2 * d3 < 0) | (d3 * d1 < 0);
    # Vertex exactly on plane, other two split:
    crosses <- crosses | (d1 == 0 & d2 * d3 < 0);
    crosses <- crosses | (d2 == 0 & d1 * d3 < 0);
    crosses <- crosses | (d3 == 0 & d1 * d2 < 0);

    crossing_idx <- which(crosses);
    other_axes <- setdiff(1:3, axis);
    n_cross <- length(crossing_idx);
    segs <- vector("list", n_cross);
    seg_count <- 0L;

    for (k in seq_len(n_cross)) {
        fi <- crossing_idx[k];
        idx <- faces[fi, ];
        vv  <- verts[idx, , drop = FALSE];
        dd  <- vv[, axis] - slice_crs_coord;

        pts <- matrix(NA_real_, 0, 2);
        for (ei in 1:3) {
            i <- ei;
            j <- if (ei == 3) 1 else ei + 1;
            if (dd[i] * dd[j] < 0 || (dd[i] == 0 && dd[j] != 0)) {
                t_val <- -dd[i] / (dd[j] - dd[i]);
                pt <- vv[i, other_axes] + t_val * (vv[j, other_axes] - vv[i, other_axes]);
                pts <- rbind(pts, pt);
            } else if (dd[i] == 0 && dd[j] == 0) {
                pts <- rbind(pts, vv[i, other_axes], vv[j, other_axes]);
            }
        }

        if (nrow(pts) >= 2) {
            seg_count <- seg_count + 1L;
            segs[[seg_count]] <- pts[1:2, , drop = FALSE];
        }
    }

    if (seg_count == 0) return(list());
    return(segs[1:seg_count]);
}


#' @title Draw contour segments onto a magick image.
#'
#' @description Draws a set of line segments (as returned by \code{\link[fsbrain]{mesh.slice.intersection}}) onto a magick image. Individual mesh triangles produce small segments, which collectively trace the contour boundary when drawn together.
#'
#' @param img a magick image instance (a single 2D slice image).
#'
#' @param segments list of 2×2 matrices, each a line segment in 0-based CRS coordinates on the two orthogonal axes.
#'
#' @param slice_axis integer, the axis (1,2,3) perpendicular to the slice plane.
#'
#' @param row_axis integer, which CRS axis (1, 2, or 3) maps to the image row dimension.
#'
#' @param col_axis integer, which CRS axis (1, 2, or 3) maps to the image column dimension.
#'
#' @param color character string, the color for the contour lines.
#'
#' @param lwd numeric, line width passed to \code{\link[graphics]{segments}}. Defaults to 1.
#'
#' @param coord_scale positive integer, the factor by which the CRS coordinates must be multiplied to map them to image pixel coordinates. Defaults to 1L, which maps coordinates 1:1 to pixels (as is the case when the image is at native resolution). Use larger values when drawing onto a scaled-up image or overlay.
#'
#' @return the modified magick image (invisibly).
#'
#' @keywords internal
draw.segments.on.image <- function(img, segments, slice_axis, row_axis, col_axis, color = "#FF0000", lwd = 1, coord_scale = 1L) {
    if (length(segments) == 0) return(invisible(img));

    other_axes <- setdiff(1:3, slice_axis);
    row_col <- match(row_axis, other_axes);
    col_col <- match(col_axis, other_axes);

    img_info <- magick::image_info(img);
    img_width  <- img_info$width;
    img_height <- img_info$height;

    img <- magick::image_draw(img);
    for (seg in segments) {
        x0 <- max(1, min(img_width,  seg[1, col_col] * coord_scale + 1));
        y0 <- max(1, min(img_height, seg[1, row_col] * coord_scale + 1));
        x1 <- max(1, min(img_width,  seg[2, col_col] * coord_scale + 1));
        y1 <- max(1, min(img_height, seg[2, row_col] * coord_scale + 1));
        graphics::segments(x0, y0, x1, y1, col = color, lwd = lwd);
    }
    dev.off();

    return(invisible(img));
}


#' @title Compute surface contour slice images (internal).
#'
#' @description Internal helper that loads a volume and surface meshes, computes the intersection contours of each surface with each slice plane, and draws them onto individual magick slice images. Used by both \code{\link{volvis.lb.with.surface}} (which grids the results) and \code{\link{volvis.slices.with.surface}} (which exports them individually).
#'
#' @param subjects_dir character string, the FreeSurfer SUBJECTS_DIR.
#'
#' @param subject_id character string, the subject identifier.
#'
#' @param volume numeric 3D array or character string. Either a 3D brain volume or the name of a volume file to load from the subject's \code{mri/} directory. Defaults to \code{"brain"}.
#'
#' @param surface character string or vector of strings, the surface(s) to use for contour extraction.
#'
#' @param hemi character string, one of \code{'lh'}, \code{'rh'}, or \code{'both'}.
#'
#' @param surface_color character string or vector, color(s) for the contour lines.
#'
#' @param surface_lwd numeric, line width for contours.
#'
#' @param slices passed to \code{\link[fsbrain]{volvis.lightbox}}. A negative integer N means "use every Nth slice". A numeric vector gives explicit slice indices (1-based). Defaults to \code{-5}.
#'
#' @param axis integer, the slice axis. 1 = sagittal, 2 = coronal, 3 = axial. Defaults to \code{1L}.
#'
#' @param scale integer >= 1, the factor by which to upscale the slice images in both dimensions using nearest-neighbor interpolation (keeps the MRI pixels sharp, which is good for QA). The surface contour lines are drawn on a high-resolution transparent overlay which is downsampled to the final size, so the lines appear thin and smooth. Defaults to \code{1L} (no upscaling, previous behavior).
#'
#' @param silent logical, whether to suppress messages.
#'
#' @return named list with entries: \code{images} (list of magick images, one per slice), \code{slice_indices} (integer vector of 1-based slice indices used), \code{has_contour} (logical vector, TRUE if at least one contour segment was drawn on the slice). Also includes metadata entries \code{axis}, \code{surfaces}, \code{hemis}, \code{subject_id} for downstream use.
#'
#' @keywords internal
compute.surface.contour.slices <- function(subjects_dir, subject_id,
    volume = "brain",
    surface = "white",
    hemi = "both",
    surface_color = "#FF0000",
    surface_lwd = 1,
    slices = -5,
    axis = 1L,
    scale = 1L,
    silent = TRUE
) {
    # ── Validate ──────────────────────────────────────────────────────────
    if(!(hemi %in% c("lh", "rh", "both"))) {
        stop("Parameter 'hemi' must be one of 'lh', 'rh', or 'both'.");
    }
    if(!(axis %in% 1:3)) {
        stop("Parameter 'axis' must be 1, 2, or 3.");
    }
    scale <- as.numeric(scale);
    if(length(scale) != 1L || is.na(scale) || scale < 1 || scale != round(scale)) {
        stop("Parameter 'scale' must be a single integer >= 1.");
    }
    scale <- as.integer(scale);

    # ── Load volume ───────────────────────────────────────────────────────
    if(is.character(volume) && length(volume) == 1L && !startsWith(volume, "#")) {
        vol_data <- subject.volume(subjects_dir, subject_id, volume);
    } else {
        vol_data <- volume;
    }
    if(!(is.array(vol_data) && length(dim(vol_data)) == 3L)) {
        stop("Parameter 'volume' must be a 3D numeric array or a volume file name.");
    }
    voldim <- dim(vol_data);

    # ── Normalize surface and surface_color to parallel vectors ───────────
    surfaces <- as.character(surface);
    nsurf <- length(surfaces);

    hemis <- if(hemi == "both") c("lh", "rh") else hemi;
    nhemi <- length(hemis);

    # Build per-surface, per-hemisphere color map
    nc <- length(surface_color);
    if(nc == 1L) {
        surf_col_map <- matrix(surface_color, nrow = nsurf, ncol = nhemi);
    } else if(nc == nsurf) {
        surf_col_map <- matrix(rep(surface_color, each = nhemi),
            nrow = nsurf, ncol = nhemi, byrow = TRUE);
    } else if(nc == nsurf * nhemi) {
        surf_col_map <- matrix(surface_color, nrow = nsurf, ncol = nhemi, byrow = TRUE);
    } else {
        stop(sprintf(paste0("Length of 'surface_color' (%d) must be 1, the ",
            "number of surfaces (%d), or the number of surfaces times the ",
            "number of hemispheres (%d)."), nc, nsurf, nsurf * nhemi));
    }
    rownames(surf_col_map) <- surfaces;
    colnames(surf_col_map) <- hemis;

    # ── Load all surfaces and transform to CRS ────────────────────────────
    surf_crs <- list();
    for (surf_name in surfaces) {
        surf_crs[[surf_name]] <- list();
        for (h in hemis) {
            surf_path <- file.path(subjects_dir, subject_id, "surf",
                paste0(h, ".", surf_name));
            if(!file.exists(surf_path)) {
                warning(sprintf("Surface file not found, skipping '%s': %s", h, surf_path));
                next;
            }
            s <- freesurferformats::read.fs.surface(surf_path);
            surf_crs[[surf_name]][[h]] <- mesh.ras2crs(s);
        }
    }

    if(length(surf_crs) == 0) {
        stop("No surface files could be loaded. Check subjects_dir, subject_id, and surface parameters.");
    }

    # ── Compute slice indices ─────────────────────────────────────────────
    slice_indices <- get.slice.indices(voldim, axis, slices);
    if(!silent) {
        message(sprintf("Computing contours for %d slices along axis %d...",
            length(slice_indices), axis));
    }

    # ── Mapping from contour CRS axes to image pixel axes ─────────────────
    if(axis == 1L) {
        row_axis <- 2L; col_axis <- 3L;
    } else if(axis == 2L) {
        row_axis <- 1L; col_axis <- 3L;
    } else {
        row_axis <- 1L; col_axis <- 2L;
    }

    # ── Convert volume to grayscale RGB and extract slices ────────────────
    vol_rgb <- vol.intensity.to.color(vol_data, scale = "normalize");

    # Supersampling factor for smooth contour lines: when the image is upscaled
    # (scale > 1), the lines are drawn on a high-resolution transparent overlay
    # which is then downsampled, yielding thin and smooth lines, while the
    # background MRI stays blocky and sharp (nearest-neighbor upscaling).
    line_ss <- 4L;

    # ── For each slice: compute contours, draw on image ──────────────────
    magick_images <- list();
    has_contour <- logical(length(slice_indices));
    for (si in seq_along(slice_indices)) {
        slice_r_idx <- slice_indices[si];       # 1-based R index
        slice_crs   <- slice_r_idx - 1;          # 0-based CRS

        slice_2d <- vol.slice(vol_rgb, slice_index = slice_r_idx, axis = axis);
        img <- magick::image_read(slice_2d);

        # Upscale the background slice with nearest-neighbor interpolation and
        # prepare a high-resolution transparent overlay for the contour lines.
        overlay <- NULL;
        if(scale > 1L) {
            nat_w <- magick::image_info(img)$width;
            nat_h <- magick::image_info(img)$height;
            img <- magick::image_resize(img, magick::geometry_size_pixels(nat_w * scale, nat_h * scale), filter = "point");
            overlay <- magick::image_blank(nat_w * scale * line_ss, nat_h * scale * line_ss, color = "transparent");
        }

        any_contour <- FALSE;
        for (surf_name in names(surf_crs)) {
            for (h in names(surf_crs[[surf_name]])) {
                segs <- mesh.slice.intersection(surf_crs[[surf_name]][[h]],
                    axis, slice_crs);
                if(length(segs) > 0) {
                    any_contour <- TRUE;
                    if(scale > 1L) {
                        overlay <- draw.segments.on.image(overlay, segs,
                            slice_axis = axis,
                            row_axis = row_axis, col_axis = col_axis,
                            color = surf_col_map[surf_name, h],
                            lwd = surface_lwd * line_ss,
                            coord_scale = scale * line_ss);
                    } else {
                        img <- draw.segments.on.image(img, segs,
                            slice_axis = axis,
                            row_axis = row_axis, col_axis = col_axis,
                            color = surf_col_map[surf_name, h], lwd = surface_lwd);
                    }
                }
            }
        }

        # Downsample the overlay to the final image size (smoothing the contour
        # lines) and composite it onto the blocky background.
        if(scale > 1L && any_contour) {
            overlay_small <- magick::image_resize(overlay, magick::geometry_size_pixels(nat_w * scale, nat_h * scale), filter = "lanczos");
            img <- magick::image_composite(img, overlay_small);
        }

        magick_images[[si]] <- img;
        has_contour[si] <- any_contour;
    }

    return(list(
        images = magick_images,
        slice_indices = slice_indices,
        has_contour = has_contour,
        axis = axis,
        surfaces = surfaces,
        hemis = hemis,
        subject_id = subject_id
    ));
}


#' @title Visualize volume slices with surface mesh contours overlaid in lightbox view.
#'
#' @description Creates a lightbox view of 2D MRI slices with the cortical surface boundary contours drawn on top. This is the fsbrain equivalent of what FreeSurfer's `freeview` or `tkmedit` show — useful for quality assessment of surface reconstruction against the underlying MRI volume. The surface mesh is loaded from the subject's FreeSurfer directory and its intersection with each displayed slice plane is computed and drawn as colored contour lines.
#'
#' @param subjects_dir character string, the FreeSurfer SUBJECTS_DIR.
#'
#' @param subject_id character string, the subject identifier.
#'
#' @param volume numeric 3D array or character string. Either a 3D brain volume (e.g., from T1.mgz or brain.mgz), or the name of a volume file to load from the subject's `mri/` directory. Defaults to \code{"brain"}.
#'
#' @param surface character string or vector of strings, the surface(s) to use for contour extraction. One or more of \code{"white"}, \code{"pial"}, or \code{"inflated"}. If a vector, contours for all surfaces are drawn. Defaults to \code{"white"}.
#'
#' @param hemi character string, one of \code{'lh'}, \code{'rh'}, or \code{'both'}. Which hemisphere surface(s) to overlay. Defaults to \code{"both"}.
#'
#' @param surface_color character string or character vector, the color(s) for the surface contour lines. If a single color, all surfaces and hemispheres use it. If a vector of length \code{length(surface)}, each surface gets its own color (both hemispheres use it). If a vector of length \code{2 * length(surface)}, colors are assigned as \code{(lh_surf1, rh_surf1, lh_surf2, rh_surf2, ...)}. Defaults to \code{"#FF0000"} (red).
#'
#' @param surface_lwd numeric, line width for the contour lines (passed as `lwd` to \code{\link[graphics]{segments}}). Defaults to 1.
#'
#' @param slices passed to \code{\link[fsbrain]{volvis.lightbox}}. A negative integer N means "use every Nth slice". A numeric vector gives explicit slice indices (1-based). Defaults to \code{-5} (every 5th slice).
#'
#' @param axis integer, the slice axis. 1 = sagittal, 2 = coronal, 3 = axial (in volume CRS convention). Defaults to \code{1L}.
#'
#' @param scale integer >= 1, the factor by which to upscale the slice images in both dimensions using nearest-neighbor interpolation (keeps the MRI pixels sharp, which is good for QA). The surface contour lines are drawn on a high-resolution overlay and downsampled to the final size, so they appear thin and smooth. Defaults to \code{1L} (no upscaling, previous behavior).
#'
#' @param per_row integer, number of slice images per row in the lightbox grid. Defaults to \code{5L}.
#'
#' @param per_col integer, number of rows. If \code{NULL}, computed from \code{per_row}. Defaults to \code{NULL}.
#'
#' @param border_geometry character string, geometry for borders between tiles, passed to \code{\link[magick]{image_border}}. Defaults to \code{"5x5"}.
#'
#' @param background_color character string, background color for borders and empty tiles. Defaults to \code{"#000000"}.
#'
#' @param silent logical, whether to suppress messages. Defaults to \code{TRUE}.
#'
#' @return a magick image instance containing the lightbox with surface contours overlaid. Can be saved with \code{\link[magick]{image_write}} or displayed interactively.
#'
#' @examples
#' \dontrun{
#'    fsbrain::download_optional_data();
#'    subjects_dir <- fsbrain::get_optional_data_filepath("subjects_dir");
#'
#'    # Axial slices with white surface contours in red:
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=3L,
#'       surface_color="#FF0000");
#'    magick::image_write(img, "~/fsbrain_qa_axial.png");
#'
#'    # Coronal slices in green:
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=2L,
#'       surface_color="#00FF00");
#'
#'    # Sagittal slices: left hemisphere red, right hemisphere blue:
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=1L,
#'       surface_color=c("#FF0000", "#0000FF"));
#'
#'    # Single hemisphere, pial surface, every 5th slice:
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="pial", hemi="lh",
#'       axis=3L, slices=-5, surface_color="#FF8800");
#'
#'    # Both white and pial surfaces overlaid: white in red, pial in yellow:
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface=c("white", "pial"), axis=3L,
#'       surface_color=c("#FF0000", "#FFFF00"));
#'
#'    # White (lh red, rh blue) + pial (lh green, rh orange):
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface=c("white", "pial"), axis=3L,
#'       surface_color=c("#FF0000", "#0000FF", "#00FF00", "#FF8800"));
#'
#'    # 2x upscaled slices for larger images with thin, smooth contours:
#'    img <- volvis.lb.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=3L, scale=2L);
#' }
#'
#' @family volume visualization
#'
#' @export
volvis.lb.with.surface <- function(subjects_dir, subject_id,
    volume = "brain",
    surface = "white",
    hemi = "both",
    surface_color = "#FF0000",
    surface_lwd = 1,
    slices = -5,
    axis = 1L,
    scale = 1L,
    per_row = 5L,
    per_col = NULL,
    border_geometry = "5x5",
    background_color = "#000000",
    silent = TRUE
) {
    # ── Compute per-slice images with surface contours ────────────────────────
    result <- compute.surface.contour.slices(
        subjects_dir = subjects_dir,
        subject_id = subject_id,
        volume = volume,
        surface = surface,
        hemi = hemi,
        surface_color = surface_color,
        surface_lwd = surface_lwd,
        slices = slices,
        axis = axis,
        scale = scale,
        silent = silent
    );

    # ── Combine into magick image stack and arrange in grid ──────────────────
    img_stack <- Reduce(c, result$images);

    if(!is.null(border_geometry) && !is.null(background_color)) {
        img_stack <- magick::image_border(img_stack, background_color, border_geometry);
    }

    merged_img <- magick.grid(img_stack, per_row = per_row, per_col = per_col,
        background_color = background_color);

    return(merged_img);
}


#' @title Export individual volume slices with surface contours to image files.
#'
#' @description Creates individual 2D slice images with surface contour overlays and exports them to separate files — ideal for browsing through slices one by one. This is the single-slice export counterpart of \code{\link{volvis.lb.with.surface}}, which arranges slice images in a lightbox grid. Supports optional slice index labelling and automatic skipping of empty slices (i.e., slices where no surface contour intersects the plane).
#'
#' @param subjects_dir character string, the FreeSurfer SUBJECTS_DIR.
#'
#' @param subject_id character string, the subject identifier.
#'
#' @param volume numeric 3D array or character string. Either a 3D brain volume, or the name of a volume file to load from the subject's \code{mri/} directory. Defaults to \code{"brain"}.
#'
#' @param surface character string or vector of strings, the surface(s) to use for contour extraction. One or more of \code{"white"}, \code{"pial"}, or \code{"inflated"}. Defaults to \code{"white"}.
#'
#' @param hemi character string, one of \code{'lh'}, \code{'rh'}, or \code{'both'}. Which hemisphere surface(s) to overlay. Defaults to \code{"both"}.
#'
#' @param surface_color character string or character vector, the color(s) for the surface contour lines. See \code{\link{volvis.lb.with.surface}} for the color assignment rules. Defaults to \code{"#FF0000"} (red).
#'
#' @param surface_lwd numeric, line width for the contour lines. Defaults to 1.
#'
#' @param slices passed to \code{\link[fsbrain]{volvis.lightbox}}. A negative integer N means "use every Nth slice". A numeric vector gives explicit slice indices (1-based). Defaults to \code{-5}.
#'
#' @param axis integer, the slice axis. 1 = sagittal, 2 = coronal, 3 = axial (in volume CRS convention). Defaults to \code{1L}.
#'
#' @param scale integer >= 1, the factor by which to upscale the slice images in both dimensions using nearest-neighbor interpolation (keeps the MRI pixels sharp, which is good for QA). The surface contour lines are drawn on a high-resolution overlay and downsampled to the final size, so they appear thin and smooth. Defaults to \code{1L} (no upscaling, previous behavior).
#'
#' @param silent logical, whether to suppress messages. Defaults to \code{TRUE}.
#'
#' @param output_dir character string, directory in which to write the slice image files. Created if it does not exist. Defaults to \code{"."} (current working directory).
#'
#' @param output_prefix character string, prefix for output filenames. If \code{NULL}, auto-derived as \code{"<subject_id>"}. Defaults to \code{NULL}.
#'
#' @param image_format character string, image format for output files (e.g., \code{"png"}, \code{"jpg"}). Passed to \code{\link[magick]{image_write}}. Defaults to \code{"png"}.
#'
#' @param label_slices logical, whether to annotate each slice image with its slice index (e.g., "slice 42") in the top-left corner. Uses \code{\link[magick]{image_annotate}}. Defaults to \code{FALSE}.
#'
#' @param label_color character string, color for the slice label text. Defaults to \code{"white"}.
#'
#' @param label_size integer, font size for the slice label. Defaults to 20.
#'
#' @param label_gravity character string, gravity for label placement, passed to \code{\link[magick]{image_annotate}}. Defaults to \code{"northwest"}.
#'
#' @param skip_empty logical, whether to skip slices where no surface contour intersects the slice plane. Defaults to \code{FALSE}.
#'
#' @return invisible character vector of file paths that were written.
#'
#' @examples
#' \dontrun{
#'    fsbrain::download_optional_data();
#'    subjects_dir <- fsbrain::get_optional_data_filepath("subjects_dir");
#'
#'    # Export all sagittal slices with white surface contours:
#'    volvis.slices.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=1L,
#'       output_dir="~/qa_slices");
#'
#'    # Export with slice labels, skipping empty slices:
#'    volvis.slices.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=3L,
#'       output_dir="~/qa_slices", label_slices=TRUE,
#'       skip_empty=TRUE);
#'
#'    # Both white and pial surfaces, different colors:
#'    volvis.slices.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface=c("white","pial"), axis=3L,
#'       surface_color=c("#FF0000","#FFFF00"),
#'       output_dir="~/qa_slices");
#'
#'    # 2x upscaled slices for larger images with thin, smooth contours:
#'    volvis.slices.with.surface(subjects_dir, "subject1",
#'       volume="brain", surface="white", axis=3L, scale=2L,
#'       output_dir="~/qa_slices");
#' }
#'
#' @family volume visualization
#'
#' @export
volvis.slices.with.surface <- function(subjects_dir, subject_id,
    volume = "brain",
    surface = "white",
    hemi = "both",
    surface_color = "#FF0000",
    surface_lwd = 1,
    slices = -5,
    axis = 1L,
    scale = 1L,
    silent = TRUE,
    output_dir = ".",
    output_prefix = NULL,
    image_format = "png",
    label_slices = FALSE,
    label_color = "white",
    label_size = 20,
    label_gravity = "northwest",
    skip_empty = FALSE
) {
    # ── Compute per-slice images with surface contours ────────────────────
    result <- compute.surface.contour.slices(
        subjects_dir = subjects_dir,
        subject_id = subject_id,
        volume = volume,
        surface = surface,
        hemi = hemi,
        surface_color = surface_color,
        surface_lwd = surface_lwd,
        slices = slices,
        axis = axis,
        scale = scale,
        silent = silent
    );

    # ── Filter empty slices if requested ──────────────────────────────────
    keep_idx <- seq_along(result$images);
    if(skip_empty) {
        keep_idx <- which(result$has_contour);
        if(length(keep_idx) == 0) {
            warning("All slices are empty (no surface contour intersects any slice plane), nothing exported.");
            return(invisible(character(0)));
        }
        if(!silent) {
            n_skipped <- length(result$images) - length(keep_idx);
            message(sprintf("Skipping %d empty slice(s), exporting %d slice(s).",
                n_skipped, length(keep_idx)));
        }
    }

    # ── Build output path components ──────────────────────────────────────
    if(is.null(output_prefix)) {
        output_prefix <- subject_id;
    }

    surf_str <- paste(result$surfaces, collapse = "_");
    hemi_str <- paste(result$hemis, collapse = "_");
    scale_str <- sprintf("scl%d", as.integer(scale));

    if(!dir.exists(output_dir)) {
        dir.create(output_dir, recursive = TRUE, showWarnings = !silent);
    }

    # ── Write each slice ──────────────────────────────────────────────────
    out_files <- character(length(keep_idx));
    for (i in seq_along(keep_idx)) {
        ki <- keep_idx[i];
        img <- result$images[[ki]];
        sl_idx <- result$slice_indices[ki];

        if(label_slices) {
            label_text <- sprintf("slice %d", sl_idx);
            img <- magick::image_annotate(img, label_text,
                size = label_size, color = label_color,
                gravity = label_gravity);
        }

        fname <- sprintf("%s_axis%d_slice%04d_%s_%s_%s.%s",
            output_prefix, result$axis, sl_idx, surf_str, hemi_str,
            scale_str, image_format);
        fpath <- file.path(output_dir, fname);
        magick::image_write(img, path = fpath, format = image_format);
        out_files[i] <- fpath;
    }

    if(!silent) {
        message(sprintf("Exported %d slice image(s) to '%s'.",
            length(out_files), output_dir));
    }

    return(invisible(out_files));
}