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