R/colourramp.R

In pvermees/IsoplotR: Statistical Toolbox for Radiometric Geochronology

delta <- 6 / 29
delta2 <- delta * delta
delta3 <- delta2 * delta

curve <- function(x) {
    if (delta3 < x) out <- x ^ (1 / 3)
    else out <- x / (3 * delta2) + 4 / 29
    out
}

anticurve <- function(y) {
    if (delta < y) out <- y ^ 3
    else out <- (y - 4 / 29) * (3 * delta2)
    out
}

rgb2xyz <- cbind(c(0.49 / 95, 0.0018, 0),
                 c(0.31 / 95, 0.0081, 0.01 / 109),
                 c(0.2 / 95, 0.0001, 0.99 / 109)
                 )

xyz2rgb <- solve(rgb2xyz)

# R, g and b from 0 to 1
rgb2ceilab <- function(rgb) {
    v <- rgb2xyz %*% rgb
    curvex <- curve(v[1])
    curvey <- curve(v[2])
    curvez <- curve(v[3])
    c(
        116 * curvey - 16,
        500 * (curvex - curvey),
        200 * (curvey - curvez)
    )
}

ceilab2rgb <- function(lab) {
    l <- lab[1]
    a <- lab[2]
    b <- lab[3]
    curvey <- (l + 16) / 116
    curved <- c(a / 500 + curvey, curvey, curvey - b / 200)
    xyz <- vapply(curved, anticurve, 0)
    xyz2rgb %*% xyz
}

clamp01 <- function(x) pmin(1, pmax(0, x))

rgbv <- function(v) {
    c <- clamp01(v)
    rgb(c[1], c[2], c[3])
}

colour_ramp <- function(
                        colour_list,
                        collapse_red_green = FALSE,
                        collapse_blue_yellow = FALSE
                        ) {
    collapse <- c(1,
                  if (collapse_red_green) 0.001 else 1,
                  if (collapse_blue_yellow) 0.001 else 1
                  )
    last_colour <- colour_list[[length(colour_list)]]
    colour_count <- length(colour_list)
    rgb_matrix <- grDevices::col2rgb(colour_list) / 255
    lab_list <- lapply(1:colour_count, function(i) {
        rgb2ceilab(rgb_matrix[, i]) * collapse
    })
    lab_diffs <- mapply(
        function(u, v) (u - v),
        lab_list[2:colour_count],
        lab_list[1:colour_count - 1]
    )
    dists <- unlist(lapply(
        1:(colour_count - 1),
        function(i) sqrt(sum(lab_diffs[, i]^2))
    ))
    dist_order <- order(dists, decreasing = TRUE)
    total_dist <- sum(dists)
    function(n) {
        max_dist <- total_dist / (n - 2)
        min_counts <- floor(dists / max_dist)
        total <- sum(min_counts)
        extras <- n - total - 1
        which_sides <- dist_order <= extras
        counts <- min_counts + which_sides
        # which colours are we interpolating for each result colour?
        i2coli <- rep(1:(colour_count - 1), times = counts)
        # how many colours are on each interpolation?
        i2count <- rep(counts, times = counts)
        isstart <- c(TRUE, i2coli[1:(n - 2)] != i2coli[2:(n - 1)])
        # how far are we interpolating for each result colour?
        i2sidei <- (1:(n - 1)) - cummax((1:(n - 1)) * isstart)
        results <- mapply(
            function(coli, sidei, count) {
                lab_list[[coli]] + sidei * lab_diffs[, coli] / count
            },
            i2coli,
            i2sidei,
            i2count,
            SIMPLIFY = FALSE
        )
        rgb_results <- lapply(results, function(v) {
            ceilab2rgb(v / collapse)
        })
        c(vapply(rgb_results, rgbv, ""), last_colour)
    }
}

colour_iridesce <- colour_ramp(
    list("black", "dark magenta", "dark blue", "green", "yellow"),
    collapse_red_green = TRUE
)

colour_salmon <- colour_ramp(
    list("black", "red", "pink"),
    collapse_red_green = TRUE
)

colour_rainbow_plus <- colour_ramp(
    list("black", "red", "yellow", "green", "blue", "magenta")
)

colour_red_blue_black <- colour_ramp(
    list("red", "magenta", "blue", "black"),
    collapse_red_green = TRUE
)