0001 %BAYER_DEMOSAIC - Converts a Bayer raw image to an MxNx3 image matrix with
0002 %                 separate red, green, and blue channels and interpolates the
0003 %                 unknown components of a channel from neighbors of the same
0004 %                 channel.
0005 %
0006 % Syntax:  R = bayer_demosaic(img, pattern, interpolate)
0007 %
0008 % Inputs:
0009 %    img - Bayer raw image's filename or an image matrix (MxNx1)
0010 %    pattern - the Bayer filter pattern of the camera that took img.
0011 %              default: uses value returned from access_bayer_pattern().
0012 %    interpolate - a boolean that if true, will result in unknown R, G, or B
0013 %                  component values being interpolated from it's most proximate
0014 %                  neighbors of the same color, and if false, will only separate
0015 %                  the Bayer image into channels, setting the unknown components
0016 %                  to zero. default: true.
0017 %
0018 % Outputs:
0019 %    R - MxNx3 image matrix with separate red, green, and blue channels
0020 %
0021 % Example:
0022 %    R = bayer_demosaic('path/myfilename.PGM', 'gbrg', true)
0023 %
0024 %    pattern = determine_bayer_pattern('ref/red.PGM', 101, 5);
0025 %    access_bayer_pattern(pattern);
0026 %    R = bayer_raw_read(I, [], false)
0027 %
0028 % Other m-files required: access_bayer_pattern, myimread
0029 % Subfunctions: bayer_interp
0030 % MAT-files required: none
0031 %
0032 % See also: BAYER_MOSAIC, ACCESS_BAYER_PATTERN
0033 %
0034 % Author: Jonathan Thomson
0035 % Work:
0036 % email:
0037 % Website: http://jethomson.wordpress.com
0038 %
0039 
0040 function R = bayer_demosaic(img, pattern, interpolate)
0041 
0042     if (nargin < 1 || nargin > 3 || isempty(img))
0043         usage('bayer_demosaic(img, pattern, interpolate)');
0044     end
0045 
0046     if (~exist('pattern','var') || isempty(pattern))
0047         pattern = access_bayer_pattern();
0048         if isempty(pattern)
0049             error(['bayer_demosaic: neither parameter pattern nor ' ...
0050                    'access_bayer_pattern() are set.']);
0051         end
0052     end
0053 
0054     if (~exist('interpolate','var') || isempty(interpolate))
0055         interpolate = true;
0056     end
0057 
0058     if ischar(img)
0059         I = myimread(img);
0060     else
0061         I = double(img);
0062     end
0063 
0064     [h, w, nd] = size(I);
0065 
0066     if (nd ~= 1)
0067         error('bayer_demosaic: image data must be MxNx1.');
0068     end
0069 
0070     % Transform pattern string into 2x2 cell array, one letter per cell.
0071     % This gives us easy to use start indices through strcmp.
0072     bfp = mat2cell(reshape(pattern, [2,2]).', [1,1], [1,1]);
0073 
0074     %red
0075     Red = zeros(h, w);
0076     [bi,bj]=find(strcmpi(bfp, 'r'));
0077     ri = bi:2:h;
0078     cj = bj:2:w;
0079     Red(ri, cj) = I(ri, cj);
0080 
0081     %green
0082     Green1 = zeros(h, w);
0083     Green2 = Green1;
0084     [bi,bj]=find(strcmpi(bfp, 'g')); % two green filters
0085     ri = bi(1):2:h;
0086     cj = bj(1):2:w;
0087     Green1(ri, cj) = I(ri, cj);
0088 
0089     ri = bi(2):2:h;
0090     cj = bj(2):2:w;
0091     Green2(ri, cj) = I(ri, cj);
0092 
0093     %blue
0094     Blue = zeros(h, w);
0095     [bi,bj]=find(strcmpi(bfp, 'b'));
0096     ri = bi:2:h;
0097     cj = bj:2:w;
0098     Blue(ri, cj) = I(ri, cj);
0099 
0100     if (interpolate == true)
0101         Red = bayer_interp(Red);     % interpolate columnwise
0102         Red = bayer_interp(Red.').'; % interpolate rowwise
0103         Green1 = bayer_interp(Green1);
0104         Green1 = bayer_interp(Green1.').';
0105         Green2 = bayer_interp(Green2);
0106         Green2 = bayer_interp(Green2.').';
0107         Green = (Green1+Green2)/2;
0108         Blue = bayer_interp(Blue);
0109         Blue = bayer_interp(Blue.').';
0110     else
0111         Green = Green1 + Green2;
0112     end
0113 
0114     R = reshape([Red Green Blue], [h w 3]);
0115 
0116 end
0117 
0118 % Interpolates out zeros by replacing them with the average of the pixel
0119 % directly above and the pixel directly below. Zeros along the top/bottom
0120 % edge are replaced with the value of the pixel that's directly below/above.
0121 % This function only interpolates columnwise. To interpolate rowwise
0122 % rotate the input 90 degrees.
0123 function Ri = bayer_interp(img_chan)
0124 
0125     h = size(img_chan, 1);
0126     B = [img_chan(1:2,:); img_chan; img_chan((h-1):h,:)];
0127     Ri = filter(triang(3), 1, B);
0128     Ri = Ri(4:(4+h-1),:);
0129 
0130 end
0131