% This programme is to find the frequency structure and the multiplexing 
% matrix (T) for a CFA pattern given in a matrix.
% The demosaicking matrix is the inverse (or pseudo-inverse) of T, which
% can be used for CFA demosaicking with Dubois's method:
%         http://www.site.uottawa.ca/~edubois/CFA/
%
%   For any commercial use of this program, a written permission can be 
%   obtained from the authors. please email: phao@eecs.qmul.ac.uk
%
%   If this program is used for research, please cite in your publications:
%
%   * Pengwei Hao, Yan Li, Zhouchen Lin, and Eric Dubois, "A Geometric
%     Method for Optimal Design of Color Filter Arrays", IEEE Transactions
%     on Image Processing, Vol.20, No.3, pp.709-722, March 2011.
%
%   * Yan Li, Pengwei Hao, and Zhouchen Lin, "The Frequency Structure Matrix: 
%     A Representation of Color Filter Arrays", International Journal 
%     of Imaging Systems and Technology, Vol.21, No.1, pp.101-106, 2011.
%
%   By Pengwei Hao and Yan Li, 2007-2012

syms R G B C M Y W real

% A collection of CFA patterns
sCFA1 = [G,R; B,G]; % Bayer CFA, 1976
sCFA2 = [W,R,W,B; B,W,R,W; W,B,W,R; R,W,B,W]; % Dillon, 1977
sCFA3 = [G,R,G,B; G,B,G,R]; % Yamanaka, 1977
sCFA4 = [G,R; G,B; B,G; R,G]; % Takashi, et al, 1985
sCFA5 = [M,C,G,Y; M,Y,G,C]; % Morimura, et al, 1986
sCFA6 = [W,G,W,G; R,W,B,W; W,G,W,G; B,W,R,W]; % Yamagami, et al, 1994
sCFA7 = [M,W; W,Y]; % Bawolek, et al, 1999
sCFA8 = [M,Y; C,G]; % Hamilton, et al, 2001
sCFA9 = [W,R; B,G]; % Gindele & Gallagher, 2002
sCFA10 = [C,M; B,Y]; % Bean, 2003
sCFA11 = [G,R; G,B; R,G; B,G]; % Lukac & Plataniotis, 2005
sCFA12 = [R,G,B,C; G,R,C,B; B,C,R,G; C,B,G,R]; % Roddy, et al, 2006
sCFA13 = [M,Y,C,G; G,C,Y,M]; % Hirakawa, et al, 2007
sCFA14 = [W,B,W,G; B,W,G,W; W,G,W,R; G,W,R,W]; % Kodak CFA2.0, 2007
sCFA15 = [R,B,G; G,R,B; B,G,R]; % Fillfactory diagonal
sCFA16 = [R+G/2, B+G/2; B+R/2, G+R/2; G+B/2, R+B/2]*2/3; % Condat, 2009
sCFA17 = [B,Y,C,Y; Y,C,M,C; C,M,G,M; Y,C,M,C]; % Hao, Li, Lin and Dubois (ours), TIP, 2011
sCFA18 = [W,R,B,W,G; W,G,W,R,B; R,B,W,G,W; G,W,R,B,W; B,W,G,W,R]; % Wang, Zhang and Hao (ours), ICIP11, 40% white 5x5

sCFA = sCFA18; % change this to show the frequency structures of a CFA pattern

disp('The CFA pattern :');
disp(sCFA);
sCFA = subs( sCFA, {C,M,Y,W}, {(G+B)/2,(R+B)/2,(R+G)/2,(R+G+B)/3} );

[fCFA, bCFA, T] = fs(sCFA);

disp('The frequency structure :');
disp( fCFA );
disp('The frequency structure (1 for nonzeros) :');
disp( bCFA );
disp('The multiplexing matrix (T) :');
disp( sym(T) );
disp('The 2-norm of T :');
disp( norm(T) );
disp('The demosaicking matrix (D) :');
disp( sym(threshold(pinv(T),2^-16)) );
disp('The 2-norm of D :');
disp( norm(pinv(T)) );
