LUT Analysis
Yuyang Zhao
Motivation
Given a Look Up Table (LUT) that represents some mapping \(f:\mathbb{R}^3\to\mathbb{R}^3\), we want to analyze some of its properties. One idea is to transform the LUT into some 2 dimensional plot. Some ways to do this include projecting all the (r,g,b) vectors onto the \(\mathbb{R}^2\) plane, or examining a cross-section plane of the LUT. The former method will be cluttered with all data points, while the latter is an infinitely thin slice that intersects with very few data points.
To get the best of both worlds, we will take a slice of the LUT bounded by two parallel planes and project this slice onto a 2D plot. We can setup the orientation of the slice by defining a normal vector \(\vec{n}\) orthogonal to a plane \(P\) crossing the origin. We can define the slice itself by translating \(P\) along \(P^\perp=\{t\vec{n}:t\in\mathbb{R}\}\) to get two new planes that bound the slice. Finally We can setup the location of the slice by translating these two planes by some vector in space.
Implementation
rm(list=ls())
vector_dot <- function(x, y) {
return(sum(x * y))
}
vector_norm <- function(x) {
return(sqrt(vector_dot(x, x)))
}
vector_proj_fn <- function(n) {
return(function(x) {
return(vector_dot(x, n)/vector_dot(n, n) * n)
})
}
vector_perp_fn <- function(n) {
return(function(x) {
return(x - vector_proj_fn(n)(x))
})
}
cross_product <- function(x, y) {
v1 <- x[2] * y[3] - x[3] * y[2]
v2 <- x[3] * y[1] - x[1] * y[3]
v3 <- x[1] * y[2] - x[2] * y[1]
return(c(v1, v2, v3))
}Let’s implement a special case where we’re not interested in relocating the slice. Notice that a slice is bounded by 2 planes orthogonal to \(\vec{n}\). Thus a point \(\vec{v}\) lies in the slice if \(proj_n(\vec{v})\) lies in the slice. We can define a helper functional to return functions that subset our data in this way relative to some normal vector \(\vec{n}\).
slice_perp_to_n_fn <- function(n) {
proj_n <- vector_proj_fn(n)
return(function(data, bounds) {
mapply(FUN=function(x, y, z) {
v <- proj_n(c(x, y, z))
norm_v <- vector_norm(v)
return(bounds[1] <= norm_v && norm_v <= bounds[2])
}, data[,1], data[,2], data[,3])
})
}Later we will set \(\vec{n}:=(1,1,1)\) which is an interesting normal vector because \[ proj_{(1,1,1)}(r,g,b) =\frac{r+g+b}{3}\begin{bmatrix}1\\1\\1\end{bmatrix} =\begin{bmatrix}\bar{v}\\\bar{v}\\\bar{v}\end{bmatrix},\ \bar{v}=mean(r,g,b) \] Thus each component of the projection becomes the arithmetic mean of (r,g,b) i.e the brightness. The goal would be to find all the points that lie in a certain brightness range bounded by bounds. What we will essentially do is tilt the cube on its point and slice up through the brightness range; best represented in the following picture:
knitr::include_graphics("grapher.png")
The blue line is the span of brightness perpendicular to the white plane. The slice will be defined within two of these planes and projected into \(\mathbb{R}^2\). The Wikipedia page on HSL also has a nice graphic for what we’re trying to achieve: 
The first step is the tilt transformation we’ll be doing to the cube. If our RGB data was transformed into HSL data, then we can further define a slice perpendicular to the lightness axis, which is similar to but not the same as brightness.
slice_perp_to_lightness <- function(hsl_data, bounds) {
bounds[1] <= hsl_data[,3] & hsl_data[,3] <= bounds[2]
}Next we will define a new projection transformation that “flattens” our slice into the plane \(P\) orthogonal to \(\vec{n}\). It’s easier to define this as a diagonal matrix \(D\) in terms of some basis for \(P\). We’ll use a change of basis matrix \(R\) to apply the transformation. To preserve angles, we want \(R\) to be a rotation matrix. A natural choice for \(R\) would be the matrix that rotates \(\hat{n}\) to \(\vec{e_3}=(0,0,1)\). This source helps us find such a matrix.
rotate_matrix <- function(D, a, b) {
I <- cbind(c(1,0,0), c(0,1,0), c(0,0,1))
a_hat <- a / vector_norm(a)
b_hat <- b / vector_norm(b)
v <- cross_product(a_hat, b_hat)
c <- vector_dot(a_hat, b_hat)
v_cross <- rbind(
c(0, -v[3], v[2]),
c(v[3], 0, -v[1]),
c(-v[2], v[1], 0)
)
R <- I + v_cross + v_cross %*% v_cross / (1 + c)
return(D %*% R)
}Now we can define a plotting function that plots some desired brightness range of RGB data. For comparison, we will also define a plotting function that plots some desired lightness range of HSL data. We can define a helper function to transform RGB data to HSL data using this source.
rgb_to_hsl <- function(lut) {
t(apply(lut, MARGIN=1, FUN=function(rgb_vec) {
M <- max(rgb_vec)
m <- min(rgb_vec)
C <- M - m
if (C == 0) {
return(c(0, 0, 0))
}
r <- rgb_vec[1]
g <- rgb_vec[2]
b <- rgb_vec[3]
if (M == r) {
H <- (g - b) / C %% 6
} else if (M == g) {
H <- (b - r) / C + 2
} else if (M == b) {
H <- (r - g) / C + 4
}
H <- 60 * H * pi / 180
if (H < 0) {
H <- 2 * pi + H
}
L <- (M + m) / 2
if (L == 1 || L == 0) {
S <- 0
} else {
S <- C / (1 - abs(2 * L - 1))
}
return(c(H, S, L))
}))
}
plot_rgb_slice_perp_to_brightness <- function(lut, pos, thickness, pointSize=1) {
n <- c(1, 1, 1)
norm_n <- vector_norm(n)
slice_perp_to_brightness <- slice_perp_to_n_fn(n)
bounds <- c(pos-thickness, pos+thickness) * norm_n
slice <- lut[slice_perp_to_brightness(lut, bounds),]
D <- cbind(c(1,0,0), c(0,1,0), c(0,0,0))
A <- rotate_matrix(D, n, c(0,0,1))
slice_flattened <- mapply(FUN=function(r, g, b) {A %*% c(r, g, b)}, slice$r, slice$b, slice$g)
slice_flattened <- matrix(slice_flattened, ncol=3, byrow=TRUE)
colouring <- mapply(FUN=function(r, g, b) {
colour_range <- 0.8
colour <- colour_range * c(r, g, b)
return(rgb(colour[1], colour[2], colour[3]))
}, slice$r, slice$g, slice$b)
plot(slice_flattened[,1], slice_flattened[,2],
main="RGB Slice", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, col=colouring, pch=20, cex=pointSize)
text(-1, -0.9, paste("brightness: ", round(pos, digits=2)))
}
plot_hsl_slice_perp_to_lightness <- function(lut, pos, thickness, pointSize=1) {
hsl_lut <- rgb_to_hsl(lut)
bounds <- c(pos-thickness, pos+thickness)
slice_hsl <- hsl_lut[slice_perp_to_lightness(hsl_lut, bounds),]
colour_range <- 0.8
slice_rgb <- colour_range * lut[slice_perp_to_lightness(hsl_lut, bounds),]
colouring <- mapply(FUN=rgb, slice_rgb[,1], slice_rgb[,2], slice_rgb[,3])
if (length(slice_hsl) == 0) {
slice_xyz <- cbind(c(), c(), c())
} else {
slice_xyz <- t(apply(slice_hsl, MAR=1, FUN=function(hsl_vec) {
x <- hsl_vec[2] * cos(hsl_vec[1])
y <- hsl_vec[2] * sin(hsl_vec[1])
z <- hsl_vec[3]
return(c(x, y, z))
}))
}
plot(slice_xyz[,1], slice_xyz[,2],
main="HSL Slice", xlim=c(-1, 1), ylim=c(-1, 1), asp=1, col=colouring, pch=20, cex=pointSize)
text(-1, -0.9, paste("lightness: ", round(pos, digits=2)))
}We now load some film LUTs and compare with the identity LUT (identity transform) to see how the data points differ. The two film LUTs we will be looking at were generated from a software called FilmConvert. Each represents a mapping from the colour space S-Gamut3.Cine/S-Log3 to the desired film colour space. k5207_neg_lut is an emulation of a Cineon scanned colour negative motion picture film stock Kodak Vision3 5207. k5207_print_lut is an emulation of 5207 printed onto motion picture print film stock. The visual difference between the two is that the former is logarithmically scaled and low contrast meant for post production work, while the latter emulates the final look of the film projected on screen. Since S-Gamut3.Cine/S-Log3 is already logarithmically scaled, we should see little difference between the identity LUT and k5207_neg_lut, and big difference with the k5207_print_lut.
d <- seq(0, 1, length.out=64)
identity_lut <- Reduce(f=function(m, r) {
gb_plane <- Reduce(f=function(m, g) {
b_line <- Reduce(f=function(m, b) {
return(rbind(m, c(r, g, b)))
}, d, init=matrix(, nrow=0, ncol=3))
return(rbind(m, b_line))
}, d, init=matrix(, nrow=0, ncol=3))
return(rbind(m, gb_plane))
}, d, init=matrix(, nrow=0, ncol=3))
identity_lut <- data.frame(r=identity_lut[,1], g=identity_lut[,2], b=identity_lut[,3])
k5207_neg_lut <- read.csv("slog3cine-k5207.cube",
sep=" ", header=FALSE, skip=5, col.names=c("r", "g", "b"))
k5207_print_lut <- read.csv("slog3cine-k5207-print.cube",
sep=" ", header=FALSE, skip=5, col.names=c("r", "g", "b"))
#lut <- k5207_print_lut[sample(nrow(k5207_print_lut), 20000),]Cross-sections of the identity RGB LUT
par(mfrow=c(3, 3), mar=rep(0.1, 4))
# Build images -> save them at .png format
#png(file="slice%02d.png", width=640, height=640)
num_images <- 9
offset <- ceiling(num_images / 2) - num_images / 2
for (i in c(1:num_images)) {
t <- (offset + i - 1) / num_images
plot_rgb_slice_perp_to_brightness(identity_lut, t, 0.01)
}
# Use image magick
#system("convert -delay 10 *.png identity_rgb_slices.gif")
# Remove png files
#file.remove(list.files(pattern=".png"))Cross-sections of the negative film RGB LUT
par(mfrow=c(3, 3), mar=rep(0.1, 4))
# Build images -> save them at .png format
#png(file="slice%02d.png", width=640, height=640)
num_images <- 9
offset <- ceiling(num_images / 2) - num_images / 2
for (i in c(1:num_images)) {
t <- (offset + i - 1) / num_images
plot_rgb_slice_perp_to_brightness(k5207_neg_lut, t, 0.01)
}
# Use image magick
#system("convert -delay 10 *.png identity_rgb_slices.gif")
# Remove png files
#file.remove(list.files(pattern=".png"))Cross-sections of the print film RGB LUT
par(mfrow=c(3, 3), mar=rep(0.1, 4))
# Build images -> save them at .png format
#png(file="slice%02d.png", width=640, height=640)
num_images <- 9
offset <- ceiling(num_images / 2) - num_images / 2
for (i in c(1:num_images)) {
t <- (offset + i - 1) / num_images
plot_rgb_slice_perp_to_brightness(k5207_print_lut, t, 0.01)
}
# Use image magick
#system("convert -delay 10 *.png identity_rgb_slices.gif")
# Remove png files
#file.remove(list.files(pattern=".png"))Cross-sections of the identity HSL LUT
par(mfrow=c(3, 3), mar=rep(0.1, 4))
# Build images -> save them at .png format
#png(file="slice%02d.png", width=640, height=640)
num_images <- 9
offset <- ceiling(num_images / 2) - num_images / 2
for (i in c(1:num_images)) {
t <- (offset + i - 1) / num_images
plot_hsl_slice_perp_to_lightness(identity_lut, t, 0.01)
}
# Use image magick
#system("convert -delay 10 *.png identity_hsl_slices.gif")
# Remove png files
#file.remove(list.files(pattern=".png"))Interestingly, you can see that mapping RGB to HSL has spread the points near origin (black) out across the plane of the cylinder. The same effect is seen on lighter colour points.
Cross-sections of the negative film HSL LUT
par(mfrow=c(3, 3), mar=rep(0.1, 4))
# Build images -> save them at .png format
#png(file="slice%02d.png", width=640, height=640)
num_images <- 9
offset <- ceiling(num_images / 2) - num_images / 2
for (i in c(1:num_images)) {
t <- (offset + i - 1) / num_images
plot_hsl_slice_perp_to_lightness(k5207_neg_lut, t, 0.01)
}
# Use image magick
#system("convert -delay 10 *.png identity_hsl_slices.gif")
# Remove png files
#file.remove(list.files(pattern=".png"))Cross-sections of the print film HSL LUT
par(mfrow=c(3, 3), mar=rep(0.1, 4))
# Build images -> save them at .png format
#png(file="slice%02d.png", width=640, height=640)
num_images <- 9
offset <- ceiling(num_images / 2) - num_images / 2
for (i in c(1:num_images)) {
t <- (offset + i - 1) / num_images
plot_hsl_slice_perp_to_lightness(k5207_print_lut, t, 0.01)
}
# Use image magick
#system("convert -delay 10 *.png identity_hsl_slices.gif")
# Remove png files
#file.remove(list.files(pattern=".png"))The 1D colour curves of negative film and print film LUTs
Another useful metric is how the LUT changes brightness. We can compute the brightness of all colour points and compare with the brightness of the identity. Furthermore, we can focus our attention on a specific 1 dimensional span of the LUT, the achromatic colours. These are the grayscale colours on the diagonal line from the origin of the cube to \((1,1,1)\).
identity_brightness <- apply(identity_lut, MAR=1, FUN=mean)
neg_brightness <- apply(k5207_neg_lut, MAR=1, FUN=mean)
print_brightness <- apply(k5207_print_lut, MAR=1, FUN=mean)
logical_achromatic_span <- apply(identity_lut, MAR=1, FUN=function(vec) {
brightness <- mean(vec)
return(brightness == vec[1] && brightness == vec[2] && brightness == vec[3])
})
identity_achromatic_brightness <- identity_lut[logical_achromatic_span, 1]
neg_achromatic_brightness <- k5207_neg_lut[logical_achromatic_span, 1]
print_achromatic_brightness <- k5207_print_lut[logical_achromatic_span, 1]
plot(identity_brightness, print_brightness,
xlab="S-Gamut3.Cine/S-Log3 brightness", ylab="output brightness", pch=".", col="darkgray")
points(identity_brightness, neg_brightness, pch=".", col="darkslateblue")
points(identity_achromatic_brightness, print_achromatic_brightness, pch=20, col="black")
points(identity_achromatic_brightness, neg_achromatic_brightness, pch=20, col="darkslateblue")
lines(identity_achromatic_brightness, identity_achromatic_brightness, col="red", lwd=3)
lines(identity_achromatic_brightness, print_achromatic_brightness, col="black", lwd=3)
lines(identity_achromatic_brightness, neg_achromatic_brightness, col="darkslateblue", lwd=3)
legend(0.55, 0.3,
legend=c("identity brightness",
"negative film brightness",
"negative achromatic film brightness",
"print film brightness",
"print achromatic film brightness"),
col=c("red", "darkslateblue", "darkslateblue", "darkgray", "black"),
pch=c(NA, 20, 20, 20, 20), lty=c(1, NA, 1, NA, 1), lwd=3, cex=0.8)
If we want to directly compare the output vectors vs input vectors (identity) on the same graph, we’d need a 6 dimensional plot. The best we can do is 3. Another way to analyze a LUT is to treat it as a vector valued function like so: \[ S=[0,1]\subseteq\mathbb{R},\ f:S^3\to S^3,\ f_r,f_g,f_b:S^3\to S\\ f(r,g,b)=\begin{bmatrix} f_r(r,g,b)\\f_g(r,g,b)\\f_b(r,g,b) \end{bmatrix} \] Then we can analyze the individual component functions separately.