raw_debayer.py

#!/usr/bin/python3

import numpy as np
from numpy.lib.stride_tricks import as_strided
from PIL import Image

infilename = "/tmp/foo/xac"
inwidth = 3280
inheight = 2464
bitspp = 10
bayer_pattern = "RGGB"

def read_raw_10_mipi():
	inbitspp = 10

	data = np.fromfile(infilename, np.uint8).reshape((inheight, inwidth * inbitspp / 8))

	# Horizontally, each row consists of 10-bit values. Every four bytes are
	# the high 8-bits of four values, and the 5th byte contains the packed low
	# 2-bits of the preceding four values. In other words, the bits of the
	# values A, B, C, D and arranged like so:
	#
	#  byte 1   byte 2   byte 3   byte 4   byte 5
	# AAAAAAAA BBBBBBBB CCCCCCCC DDDDDDDD AABBCCDD
	#
	# Here, we convert our data into a 16-bit array, shift all values left by
	# 2-bits and unpack the low-order bits from every 5th byte in each row,
	# then remove the columns containing the packed bits

	data = data.astype(np.uint16) << 2
	for byte in range(4):
		data[:, byte::5] |= ((data[:, 4::5] >> ((4 - byte) * 2)) & 0b11)
	data = np.delete(data, np.s_[4::5], 1)

	return data


def read_raw_12_mipi():
	inbitspp = 12

	data = np.fromfile(infilename, np.uint8)

	# XXX we have packed 12 bit raw mipi, but in a buffer with stride for 16 bit pixel
	# XXX delete the empty part on the right
	data = data.reshape((inheight, inwidth * 16 // 8))
	data = np.delete(data, np.s_[inwidth * inbitspp // 8:], 1)

	#data = data.reshape((inheight, inwidth * inbitspp // 8))

	# Horizontally, each row consists of 10-bit values. Every four bytes are
	# the high 8-bits of four values, and the 5th byte contains the packed low
	# 2-bits of the preceding four values. In other words, the bits of the
	# values A, B, C, D and arranged like so:
	#
	#  byte 1   byte 2   byte 3   byte 4   byte 5
	# AAAAAAAA BBBBBBBB CCCCCCCC DDDDDDDD AABBCCDD
	#
	# Here, we convert our data into a 16-bit array, shift all values left by
	# 2-bits and unpack the low-order bits from every 5th byte in each row,
	# then remove the columns containing the packed bits

	# AAAAAAAA BBBBBBBB AAAABBBB

	data = data.astype(np.uint16) << (inbitspp - 8)

	for byte in range(2):
		data[:, byte::3] |= (data[:, 2::3] >> (4 * byte + 4)) & 0b1111

	data = np.delete(data, np.s_[2::3], 1)

	return data

def read_raw_12_unpacked():
	data = np.fromfile(infilename, np.uint16).reshape((inheight, inwidth))

	return data

idx = bayer_pattern.find("R")
assert(idx != -1)
r0 = (idx % 2, idx // 2)

idx = bayer_pattern.find("B")
assert(idx != -1)
b0 = (idx % 2, idx // 2)

idx = bayer_pattern.find("G")
assert(idx != -1)
g0 = (idx % 2, idx // 2)

idx = bayer_pattern.find("G", idx + 1)
assert(idx != -1)
g1 = (idx % 2, idx // 2)

def separate_components(data):
	# Now to split the data up into its red, green, and blue components. The
	# Bayer pattern of the OV5647 sensor is BGGR. In other words the first
	# row contains alternating green/blue elements, the second row contains
	# alternating red/green elements, and so on as illustrated below:
	#
	# GBGBGBGBGBGBGB
	# RGRGRGRGRGRGRG
	# GBGBGBGBGBGBGB
	# RGRGRGRGRGRGRG
	#
	# Please note that if you use vflip or hflip to change the orientation
	# of the capture, you must flip the Bayer pattern accordingly

	rgb = np.zeros(data.shape + (3,), dtype=data.dtype)
	rgb[r0[1]::2, r0[0]::2, 0] = data[r0[1]::2, r0[0]::2] # Red
	rgb[g0[1]::2, g0[0]::2, 1] = data[g0[1]::2, g0[0]::2] # Green
	rgb[g1[1]::2, g1[0]::2, 1] = data[g1[1]::2, g1[0]::2] # Green
	rgb[b0[1]::2, b0[0]::2, 2] = data[b0[1]::2, b0[0]::2] # Blue

	return rgb

def nodemosaic(rgb):
	#output = np.empty(rgb.shape, dtype=rgb.dtype)
	return rgb

def demosaic(rgb):
	# At this point we now have the raw Bayer data with the correct values
	# and colors but the data still requires de-mosaicing and
	# post-processing. If you wish to do this yourself, end the script here!
	#
	# Below we present a fairly naive de-mosaic method that simply
	# calculates the weighted average of a pixel based on the pixels
	# surrounding it. The weighting is provided b0[1] a b0[1]te representation of
	# the Bayer filter which we construct first:

	bayer = np.zeros(rgb.shape, dtype=np.uint8)
	bayer[r0[1]::2, r0[0]::2, 0] = 1 # Red
	bayer[g0[1]::2, g0[0]::2, 1] = 1 # Green
	bayer[g1[1]::2, g1[0]::2, 1] = 1 # Green
	bayer[b0[1]::2, b0[0]::2, 2] = 1 # Blue

	# Allocate an array to hold our output with the same shape as the input
	# data. After this we define the size of window that will be used to
	# calculate each weighted average (3x3). Then we pad out the rgb and
	# bayer arrays, adding blank pixels at their edges to compensate for the
	# size of the window when calculating averages for edge pixels.

	output = np.empty(rgb.shape, dtype=rgb.dtype)
	window = (3, 3)
	borders = (window[0] - 1, window[1] - 1)
	border = (borders[0] // 2, borders[1] // 2)

	#rgb_pad = np.zeros((
	#	rgb.shape[0] + borders[0],
	#	rgb.shape[1] + borders[1],
	#	rgb.shape[2]), dtype=rgb.dtype)
	#rgb_pad[
	#	border[0]:rgb_pad.shape[0] - border[0],
	#	border[1]:rgb_pad.shape[1] - border[1],
	#	:] = rgb
	#rgb = rgb_pad
	#
	#bayer_pad = np.zeros((
	#	bayer.shape[0] + borders[0],
	#	bayer.shape[1] + borders[1],
	#	bayer.shape[2]), dtype=bayer.dtype)
	#bayer_pad[
	#	border[0]:bayer_pad.shape[0] - border[0],
	#	border[1]:bayer_pad.shape[1] - border[1],
	#	:] = bayer
	#bayer = bayer_pad

	# In numpy >=1.7.0 just use np.pad (version in Raspbian is 1.6.2 at the
	# time of writing...)
	#
	rgb = np.pad(rgb, [
		(border[0], border[0]),
		(border[1], border[1]),
		(0, 0),
		], 'constant')
	bayer = np.pad(bayer, [
		(border[0], border[0]),
		(border[1], border[1]),
		(0, 0),
		], 'constant')

	# For each plane in the RGB data, we use a nifty numpy trick
	# (as_strided) to construct a view over the plane of 3x3 matrices. We do
	# the same for the bayer array, then use Einstein summation on each
	# (np.sum is simpler, but copies the data so it's slower), and divide
	# the results to get our weighted average:

	for plane in range(3):
		p = rgb[..., plane]
		b = bayer[..., plane]
		pview = as_strided(p, shape=(
		p.shape[0] - borders[0],
		p.shape[1] - borders[1]) + window, strides=p.strides * 2)
		bview = as_strided(b, shape=(
		b.shape[0] - borders[0],
		b.shape[1] - borders[1]) + window, strides=b.strides * 2)
		psum = np.einsum('ijkl->ij', pview)
		bsum = np.einsum('ijkl->ij', bview)
		output[..., plane] = psum // bsum

	return output

data = read_raw_12_unpacked()
#data = read_raw_10_mipi()
#data = read_raw_12_mipi()

rgb = separate_components(data)
output = demosaic(rgb)
#output = nodemosaic(rgb)

output = (output >> (bitspp - 8)).astype(np.uint8)

Image.fromarray(output).save("out.png")