Nearest Mean Interpolation: a method to scale up images

Here is a nice and clean method of scaling up an image to double its size:

screenshot of a comment in pepperpickle's source code

I found it in the research paper that describes the HCRHide image steganography algorithm: https://www.researchgate.net/publication/272307744_A_high-capacity_reversible_data_hiding_method_HCRHide

The main quality of this algorithm is that the original image can be reconstructed without any loss of color information. Pixels a, b, c and d are taken from the original image. the three lambda-pixels we get for each original pixel can be used to store information steganographically. Here is my golang implementation of the NMI scaling-algorithm:

// scaleNMI scales up an image using the NMI method (Nearest Mean
// Interpolation)
/*
Here is how NMI works:
four neighbouring pixels of the original image are obtained:
 _________
|    |    |
| a  | b  |
|____|____|
|    |    |
| c  | d  |
|____|____|

between those pixels, new ones are inserted:
 ______________
|    |    |    |
| a  | λ₁ | b  |
|____|____|____|
|    |    |    |
| λ₂ | λ₃ | λ₄ |
|____|____|____|
|    |    |    |
| c  | λ₅ | d  |
|____|____|____|

The color values for the new pixels are calculated using the following
formulae:

λ₁=(a+b)/2
λ₂=(a+c)/2
λ₃=(a+b+c+d)/4

Those are calculated in the next pixel set as λ₁ and λ₂:
λ₄=(b+d)/2
λ₅=(c+d)/2

The original formula for λ₃ is:
λ₃=(a+λ₁+λ₂)/3

Per pass, the image is scaled up to double its size.
*/
func scaleNMI(src image.Image) image.Image {
        bounds := src.Bounds()
        w, h := bounds.Max.X, bounds.Max.Y
        dest := image.NewRGBA(image.Rect(0, 0, w*2, h*2))
        dX, dY := 0, 0

        for x := 0; x < w; x++ {
        for y := 0; y < h; y++ {
        a := src.At(x, y)
        b := src.At(x+1, y)
        c := src.At(x, y+1)
        d := src.At(x+1, y+1)

        if x > w-2 {
        b = a
        d = a
        }
        if y > h-2 {
        c = a
        d = a
        }

        aR, aG, aB, aA := a.RGBA()
        bR, bG, bB, bA := b.RGBA()
        cR, cG, cB, cA := c.RGBA()
        dR, dG, dB, dA := d.RGBA()

        l1R, l1G, l1B, l1A :=
        ((aR + bR) / 2), // >> 1: divide by 2
        ((aG + bG) / 2),
        ((aB + bB) / 2),
        ((aA + bA) / 2)
        l2R, l2G, l2B, l2A :=
        ((aR + cR) / 2),
        ((aG + cG) / 2),
        ((aB + cB) / 2),
        ((aA + cA) / 2)
        l3R, l3G, l3B, l3A :=
        ((aR + bR + cR + dR) / 4), // >> 2: divide by 4
        ((aG + bG + cG + dG) / 4),
        ((aB + bB + cB + dB) / 4),
        ((aA + bA + cA + dA) / 4)

        dest.Set(dX, dY, color.RGBA{
        uint8(aR / 256), // >> 8 divide by 256
        uint8(aG / 256),
        uint8(aB / 256),
        uint8(aA / 256)})
        dest.Set(dX+1, dY, color.RGBA{
        uint8(l1R / 256),
        uint8(l1G / 256),
        uint8(l1B / 256),
        uint8(l1A / 256)})
        dest.Set(dX, dY+1, color.RGBA{
        uint8(l2R / 256),
        uint8(l2G / 256),
        uint8(l2B / 256),
        uint8(l2A / 256)})
        dest.Set(dX+1, dY+1, color.RGBA{
        uint8(l3R / 256),
        uint8(l3G / 256),
        uint8(l3B / 256),
        uint8(l3A / 256)})

        dY += 2
        }
        dY = 0
        dX += 2
        }

        return dest
}