function [c,s,varagout] = maxchosenimf(c1,s1,c2,s2,consistency,app)
%  MAXCHOSENIMF利用最大值選擇法對圖像進行融合
%
%  [c,s] = MAXCHOSENIMF(C1,S1,C2,S2,CONSISTENCY)對利用最大值選擇法對
%  輸入進行融合,并返回融合小波系數和簿記矩陣
%  [C,S,MAP] = MAXCHOSENIMF(C1,S1,C2,S2,CONSISTENCY)返回融合小波系數,
%  簿記矩陣和融合map
%  輸出map記錄融合小波系數來自那幅圖像,若來自圖像1,則記0;
%  否則記1
%
%  app的默認值為0
%    app = 0 將兩幅輸入圖像的近似系數進行平均
%    app = 1 取兩幅圖像中近似系數的絕對值較小的作為融合系數
%    app = 2 取兩幅圖像中近似系數的絕對值較大的作為融合系數
%  consistency的默認值為1
%    consistency = 0 不進行一致性檢驗
%    consistency = 1 進行一致性檢驗
%
% 參見minchosenimf

%  檢查輸入和輸出的合理性
error(nargchk(2,3,nargout));
error(nargchk(4,6,nargin));

sz1 = size(s1);
nmax = sz1(1) - 2;
tmap = cell(nmax,3);

elements1 = prod(s1,2);
elements2 = prod(s2,2);
if (ndims(c1) ~= 2)|(size(c1,1) ~= 1)|...
        (ndims(c2) ~= 2)|(size(c2,1) ~= 1)
    error('C1和C2必須是行向量');
end 

if (ndims(s1) ~= 2)|~isreal(s1)|~isnumeric(s1)|(size(s1,2) ~= 2)|...
        (ndims(s2) ~= 2)|~isreal(s2)|~isnumeric(s2)|(size(s2,2) ~= 2)
    error('S1和S2必須是二維實數陣列'); 
end     
    
if (length(c1) < elements1(end))| ...
        ~(elements1(1)+3*sum(elements1(2:end - 1)) >= elements1(end))| ...
        (length(c2) < elements2(end))| ...
        ~(elements2(1)+3*sum(elements2(2:end - 1)) >= elements2(end))
    error('[C1,S1]和[C2,S2]必須是標準小波分解結構');
end

if (length(c1) ~= length(c2))|(s1 ~= s2)
    error('輸入的源圖像不符合要求:大小不同!')
end

if(nargin > 4) &(consistency ~= 0) & (consistency ~= 1) 
    error('consistency(是否進行一致性檢驗)輸入錯誤!');
end

if (nargin > 5) &  (~((app == 0)|(app == 1)|(app == 2)))
    error('輸入錯誤,請重新輸入對近似圖像的處理方法');
end

if nargin == 5
    app = 0;
end
 
if nargin == 4
    consistency = 1;
    app = 0;
end

%選取絕對值最大的作為融合系數
decision1 = abs(c1) > abs(c2);
decision2 = abs(c1) <= abs(c2);
decision1 = im2double(decision1);
decision2 = im2double(decision2);


s = s1;
elements = prod(s,2);
%創建map若融合系數來自圖1,則寫入0,若來自圖2則寫入1
map = ones(1,length(c1));
map = map - decision1;
mask = ones(3,3)/9;



if app == 1
    decision1(1:elements1(1)) = 1 - decision1(1:elements1(1));
    decision2(1:elements2(1)) = 1 - decision2(1:elements2(1));
end

if app == 2
    decision1(1:elements1(1)) = decision1(1:elements1(1));
    decision2(1:elements2(1)) = decision2(1:elements2(1));
end

if consistency == 1
    for i =1:nmax
        for j = 1:3
            tmap{i,j} = repmat(0,s(i + 1,:));
            tmap{i,j}(:) = map(elements1(1) + (j - 1)*elements1(2) + 1:elements1(1) + j*elements1(2));
            tmap{i,j} = impad(double(tmap{i,j}),1);
            tmap{i,j} = round(conv2(tmap{i,j},mask,'valid'));     
            map(elements1(1) + (j - 1)*elements1(2) + 1:elements1(1) + j*elements1(2)) = tmap{i,j}(:);  
        end
        elements1(1) = elements1(1) + 3*elements1(i + 1);
        elements1(2) = elements1(i + 2);   
    end
    if app == 1
        mapa = repmat(0,s(1,:));
        a = map(1:elements(1));
        mapa(:) = map(1:elements(1));
        mapa = impad(double(mapa),1);
        mapa = round(conv2(mapa,mask,'valid')); 
        mapa = ~mapa;
        mapa = impad(double(mapa),1);
        mapa = round(conv2(mapa,mask,'valid'));
        mapa = ~mapa;
        map(1:elements(1)) = mapa(:);          
    end
    decision2 = map;
    decision1 = 1 - map;
end

if app == 0
    decision1(1:elements1(1)) = 0.5;
    decision2(1:elements2(1)) = 0.5;
end

c = c1.*decision1 + c2.*decision2;