What is the correct method to demosaic raw images from a Lytro camera?

Question

I have a Lytro camera and images retrieved from it, and am trying to read the image. The end goal is that I want to view the images on a Looking Glass display, so I will eventually have to deal with the light field data too.

So far, I have read the greyscale image from the file, and am trying to convert it to colour.

This web page says of the file format,

Do not forget that the resulting image is grayscale with Bayer color filter over it, so it needs to be demosaiced to obtain the color information.

All well and good, but the page doesn't go on to explain how this is supposed to be done, which Bayer filter it's using, and so on, and so forth.

The code I have for this so far:

import java.awt.image.BufferedImage
import kotlin.math.roundToInt

enum class MosaicCell { R, GR, GB, B }

data class MosaicInfo(
    val cell: MosaicCell,
    val xOffset: Int,
    val yOffset: Int,
    val pixelFormatBlack: Int,
    val pixelFormatWhite: Int,
)

data class IntOffset(
    val x: Int,
    val y: Int,
)

enum class Direction(val offset: IntOffset) {
    LEFT(IntOffset(-1, 0)),
    RIGHT(IntOffset(1, 0)),
    UP(IntOffset(0, -1)),
    DOWN(IntOffset(0, 1)),
    UP_LEFT(IntOffset(-1, -1)),
    UP_RIGHT(IntOffset(1, -1)),
    DOWN_LEFT(IntOffset(-1, 1)),
    DOWN_RIGHT(IntOffset(1, 1)),
    ;

    companion object {
        val Horizontal get() = listOf(LEFT, RIGHT)
        val Vertical get() = listOf(UP, DOWN)
        val Orthogonal get() = Horizontal + Vertical
        val Diagonal get() = listOf(UP_LEFT, UP_RIGHT, DOWN_LEFT, DOWN_RIGHT)
    }
}

/**
 * Removes Bayer filtering from an image, retrieving the full colour image from the greyscale.
 *
 * This method allows poking the mosaic values in directly, which is more convenient in
 * unit test code.
 * 
 * @param input the input greyscale image.
 * @param mosaicValues a list of programmatically provided mosaic values.
 * @return the demosaiced image.
 */
fun demosaicImage(input: BufferedImage, mosaicValues: List<MosaicInfo>): BufferedImage {
    val mosaicCellsByModOffset = mosaicValues
        .associateBy { IntOffset(it.xOffset, it.yOffset) }
        .mapValues { (_, v) -> v.cell }

    val width = input.width
    val xRange = 0..<width
    val height = input.height
    val yRange = 0..<height

    val result = BufferedImage(width, height, BufferedImage.TYPE_INT_RGB)

    for (y in yRange) {
        for (x in xRange) {
            // Which mosaic cell are we pointing at?
            val cell = mosaicCellsByModOffset[IntOffset(x % 2, y % 2)]
                ?: throw IllegalStateException("Missing map entry!")

            fun sample(x: Int, y: Int) = input.getRGB(x, y).and(0xFF)

            fun currentPixel() = sample(x, y)

            fun averageOfNeighbours(directions: List<Direction>): Int {
                val v = directions.mapNotNull { direction ->
                    val dx = direction.offset.x
                    val dy = direction.offset.y
                    val x1 = x + dx
                    val y1 = y + dy
                    if (x1 in xRange && y1 in yRange) {
                        sample(x1, y1)
                    } else {
                        null
                    }
                }
                check(v.isNotEmpty()) { "No neighbours were in range somehow. x = $x, y = $y, directions = $directions" }
                return v.average().roundToInt()
            }

            val r: Int
            val g: Int
            val b: Int

            when (cell) {
                MosaicCell.R -> {
                    r = currentPixel()
                    g = averageOfNeighbours(Direction.Orthogonal)
                    b = averageOfNeighbours(Direction.Diagonal)
                }

                MosaicCell.GR -> {
                    r = averageOfNeighbours(Direction.Horizontal)
                    g = currentPixel()
                    b = averageOfNeighbours(Direction.Vertical)
                }

                MosaicCell.GB -> {
                    r = averageOfNeighbours(Direction.Vertical)
                    g = currentPixel()
                    b = averageOfNeighbours(Direction.Horizontal)
                }

                MosaicCell.B -> {
                    r = averageOfNeighbours(Direction.Diagonal)
                    g = averageOfNeighbours(Direction.Orthogonal)
                    b = currentPixel()
                }
            }

            assert(r in 0..255 && g in 0..255 && b in 0..255)

            result.setRGB(x, y, 0xFF000000.toInt().or(r.shl(16)).or(g.shl(8)).or(b))
        }
    }

    return result
}

My input file (tiny fragment of the full image because SO's image upload is so restrictive):

My result:

As we can see, the current result is overall too yellow, and there still appears to be a lot of noise.

I also found another question related to demosaic which included a simpler sample greyscale image and the expected result - I compared my result with theirs, and mine is the same except for producing sharper detail at the edges of the image.

Additionally, you can see from the greyscale image, there is some kind of pattern every 4 pixels, but the Bayer filter can only explain a repeating pattern every 2 pixels. So the next path I can investigate is double and triple checking the greyscale reading code, to make sure that's correct, but I don't have any good way to verify it because I don't have any known good data to compare my results with.

Additional data:

• The original raw LFR file