Demosaicing

Bilinear interpolation

% demosaicing photos in Matlab

% using bilinear interpolation

%

% input variables: image_in (mosaiced image (1 layer ), any data type)

% output variables: image_out ( RGB image (3 layers), double )

% settings : sens_coord

%

% coordinates of color sensors in elementary array: [v h] for R, G1, G2 and B:

% (usual pattern is R G2; in this case, sens_coord is [1 1]

% G1 B [2 1]

% [1 2]

% [2 2])

sens_coord = [1 1; 2 1; 1 2; 2 2];

%

% filter: 3*3 square without the center:

filter_ring8 = [1 1 1; 1 0 1; 1 1 1];

% output image:

image_out = zeros([size(image_in) 3], 'double');

% converting part-by-part to prevent memory overflow:

size_part_v = 200;

for v_beg_out = 1 : size_part_v : size(image_in, 1)

v_end_out = min(v_beg_out+size_part_v-1, size(image_in, 1));

v_beg_in = max(v_beg_out-4, 1);

v_end_in = min(v_end_out+4, size(image_in, 1));

img_part = double(image_in(v_beg_in : v_end_in, :));

for z_out = 1 : 3

% current color layer:

out_part{z_out} = zeros(size(img_part), 'double');

% map of original pixels:

map_n_part = zeros(size(img_part), 'double');

% picking of original pixels into the current color layer and the map:

z_in_list = {[1] [2 3] [4]};

for z_in = z_in_list{z_out}

out_part{z_out}(sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end) =...

img_part( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end);

map_n_part( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end) = 1;

end

% interpolated data (values for original pixels are wrong, but will not be used):

out_part_interp =...

imfilter(out_part{z_out}, filter_ring8, 'replicate') ./...

imfilter(map_n_part, filter_ring8, 'replicate');

% where original values are absent, filling by interpolated values:

map_n_part = (map_n_part == 0);

out_part{z_out}(map_n_part) = out_part_interp(map_n_part);

clear out_part_interp map_n_part

end

img_part = cat(3, out_part{1}, out_part{2}, out_part{3});

% removing of added borders:

img_part = img_part(v_beg_out - v_beg_in + 1 : end - (v_end_in - v_end_out), :, :);

image_out(v_beg_out : v_end_out, :, :) = img_part;

clear img_part out_part

end

Interpolation with Matlab's function griddata

(triangulation-based cubic interpolation and other methods)

% demosaicing photos in Matlab

% using function griddata

% possible methods of interpolation:

% "triangulation-based linear interpolation" ('linear')

% "triangulation-based nearest neighbor interpolation" ('nearest')

% "triangulation-based natural neighbor interpolation" ('natural')

% "triangulation-based cubic interpolation" ('cubic')

% "biharmonic spline interpolation (MATLAB 4 griddata method)" ('v4')

% ('natural' is almost identical to bilinear interpolation, and 'linear' differs from both;

% 'cubic' is the best one).

%

% input variables: image_in (mosaiced image (1 layer ), any data type)

% output variables: image_out ( RGB image (3 layers), double )

% settings : sens_coord, interpolation_method

%

% coordinates of color sensors in elementary array: [v h] for R, G1, G2 and B:

% (usual pattern is R G2; in this case, sens_coord is [1 1]

% G1 B [2 1]

% [1 2]

% [2 2])

sens_coord = [1 1; 2 1; 1 2; 2 2];

% interpolation method (for Matlab's function griddata):

% (possible values: 'linear', 'nearest', 'natural', 'cubic', 'v4')

interpolation_method = 'cubic';

% output image:

image_out = zeros([size(image_in) 3], 'double');

% converting part-by-part to prevent memory overflow:

size_part_v = 200;

for v_beg_out = 1 : size_part_v : size(image_in, 1)

v_end_out = min(v_beg_out+size_part_v-1, size(image_in, 1));

v_beg_in = max(v_beg_out-4, 1);

v_end_in = min(v_end_out+4, size(image_in, 1));

img_part = double(image_in(v_beg_in : v_end_in, :));

for z_out = 1 : 3

% current color layer:

out_part{z_out} = zeros(size(img_part), 'double');

% map of original pixels:

map_n_part = false(size(img_part));

% picking of original pixels into the current color layer and the map:

z_in_list = {[1] [2 3] [4]};

for z_in = z_in_list{z_out}

out_part{z_out}(sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end) =...

img_part( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end);

map_n_part( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end) = true;

end

% replacing of missing values with interpolated ones:

[v_in_list h_in_list ] = find( map_n_part);

[v_out_list h_out_list] = find(~map_n_part);

out_part{z_out}(~map_n_part(:)) = griddata(...

h_in_list, v_in_list, out_part{z_out}(map_n_part(:)),...

h_out_list, v_out_list, interpolation_method);

clear map_n_part v_in_list h_in_list v_out_list h_out_list

end

img_part = cat(3, out_part{1}, out_part{2}, out_part{3});

% replacing of NaNs (generated by some interpolation algorithms at the borders)

% with neighbouring values:

part_border = img_part(:, [1 end ], :);

part_subborder = img_part(:, [2 end-1], :);

part_border(isnan(part_border)) = part_subborder(isnan(part_border));

img_part(:, [1 end], :) = part_border; % on the left and right border

part_border = img_part([1 end ], :, :);

part_subborder = img_part([2 end-1], :, :);

part_border(isnan(part_border)) = part_subborder(isnan(part_border));

img_part([1 end], :, :) = part_border; % on the upper and lower border

clear part_border part_subborder

% removing of added borders:

img_part = img_part(v_beg_out - v_beg_in + 1 : end - (v_end_in - v_end_out), :, :);

image_out(v_beg_out : v_end_out, :, :) = img_part;

clear img_part out_part

end

Gradient-corrected linear interpolation

after Malvar et al., 2004

% demosaicing photos in Matlab

% using gradient-corrected linear interpolation after Malvar et al., 2004:

% Malvar, Henrique S., Li-wei He, and Ross Cutler. “High-Quality Linear Interpolation

% for Demosaicing of Bayer-Patterned Color Images.” Microsoft Research, May 2004.

% goo.gl/NnL1XU

%

% input variables: image_in (mosaiced image (1 layer ), any data type)

% output variables: image_out ( RGB image (3 layers), double )

% settings : sens_coord

%

% coordinates of color sensors in elementary array: [v h] for R, G1, G2 and B:

% (usual pattern is R G2; in this case, sens_coord is [1 1]

% G1 B [2 1]

% [1 2]

% [2 2])

sens_coord = [1 1; 2 1; 1 2; 2 2];

%

% filters:

% first index: place of interpolation

% (1 = R, 2 = green in B row, R column, 3 = green in R row, B column, 4 = B)

% second index: result of interpolation

% (1 = R, 2 = G, 3 = B)

%

% filters for G output:

filter{1, 2} = [...

0 0 -1 0 0;...

0 0 2 0 0;...

-1 2 4 2 -1;...

0 0 2 0 0;...

0 0 -1 0 0] / 8;

filter{4, 2} = filter{1, 2};

% filters for R output:

filter{3, 1} = [...

0 0 0.5 0 0;...

0 -1 0 -1 0;...

-1 4 5 4 -1;...

0 -1 0 -1 0;...

0 0 0.5 0 0] / 8;

filter{2, 1} = filter{3, 1}';

filter{4, 1} = [...

0 0 -1.5 0 0;...

0 2 0 2 0;...

-1.5 0 6 0 -1.5;...

0 2 0 2 0;...

0 0 -1.5 0 0] / 8;

% filters for B output:

filter{2, 3} = filter{3, 1};

filter{3, 3} = filter{2, 1};

filter{1, 3} = filter{4, 1};

% output image:

image_out = zeros([size(image_in) 3], 'double');

% converting part-by-part to prevent memory overflow:

size_part_v = 200;

for v_beg_out = 1 : size_part_v : size(image_in, 1)

v_end_out = min(v_beg_out+size_part_v-1, size(image_in, 1));

v_beg_in = max(v_beg_out-4, 1);

v_end_in = min(v_end_out+4, size(image_in, 1));

img_part = double(image_in(v_beg_in : v_end_in, :));

for z_out = 1 : 3

% current color layer:

out_part{z_out} = zeros(size(img_part), 'double');

% picking of original values:

z_in_list = {[1] [2 3] [4]};

for z_in = z_in_list{z_out}

out_part{z_out}( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end) =...

img_part( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end);

end

% where original values are absent, filling by interpolated values:

z_in_list = {[2 3 4] [1 4] [1 2 3]};

for z_in = z_in_list{z_out}

color_part_curr = imfilter(img_part, filter{z_in, z_out}, 'replicate');

out_part{z_out}( sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end) =...

color_part_curr(sens_coord(z_in, 1) : 2 : end, sens_coord(z_in, 2) : 2 : end);

clear color_part_curr

end

end

img_part = cat(3, out_part{1}, out_part{2}, out_part{3});

% removing of added borders:

img_part = img_part(v_beg_out - v_beg_in + 1 : end - (v_end_in - v_end_out), :, :);

image_out(v_beg_out : v_end_out, :, :) = img_part;

clear img_part out_part

end