function dFusionMaxEntropy = FusionDetectorsMaxEntropy(HyperCube);
% 基于最大熵準則的多檢測器融合檢測算法
% INPUT -- 多個檢測器檢測結果組成的立方體
% OUTPUT -- 根據最大熵準則得到的檢測器融合結果，二維圖像矩陣
[NumRows,NumCols,NumBands] = size(HyperCube);
indexLength = fix(NumBands/3);
SubHyperCube1 = HyperCube(:,:,1:indexLength);
SubHyperCube2 = HyperCube(:,:,indexLength+1:indexLength*2);
SubHyperCube3 = HyperCube(:,:,indexLength*2+1:end);


dCube(:,:,1) = RXcov_detector(SubHyperCube1); %%% dRXcov
dCube(:,:,2) = RXcov_detector(SubHyperCube2);
dCube(:,:,3) = RXcov_detector(SubHyperCube3);

[num_row,num_col,NumDetectors]= size(dCube);

for iDetector = 1 : NumDetectors
    temp = reshape(dCube(:,:,iDetector),1,num_row*num_col);
%     estPdCube(iDetector,:) = funcAdaptiveKernel_LEN(temp,10,temp,0.5,1,'KDE',10);%估計每個檢測器輸出統計量的概率密度--邊緣分布Marginal Distribution，每個輸出值都得到一個對應的概率值
    for j = 1:num_row*num_col
        estPdCube(iDetector,j) = KernelDenEst_LEN(temp,temp(j),20,1);
    end
    estCDFCube(iDetector,:) = funcEstCDF(temp,temp); %計算累計概率(Cumulative distribution function,CDF),對每一個檢測器的每一個輸出量都得到對應的累計概率值
    zGaussSpaceTAdjust(iDetector,:)  = norminv(estCDFCube(iDetector,:)*(num_row*num_col/(num_row*num_col+1)),0,1); % 計算將檢測器輸出向高斯分布投影空間投影后對應的z值，投影到高斯分布空間是為了滿足最大熵原則
end


zCovMatrix = zGaussSpaceTAdjust*zGaussSpaceTAdjust';
detC = (det(zCovMatrix))^(-1);
IMatrix = diag(ones(1,NumDetectors));


for iPixel = 1:num_row*num_col
    iPixel
    JointMarginal(iPixel) = prod(estPdCube(:,iPixel));% 將第iPixel個象素點對應的多個檢測器的輸出統計量的邊緣概率密度連乘
    pz(iPixel) = exp(-0.5*zGaussSpaceTAdjust(:,iPixel)'*pinv(zCovMatrix)*zGaussSpaceTAdjust(:,iPixel));
    eItem(iPixel) = exp(-0.5*zGaussSpaceTAdjust(:,iPixel)'*pinv(zCovMatrix-IMatrix)*zGaussSpaceTAdjust(:,iPixel));
    dFusionMaxEntropy(iPixel) = detC*eItem(iPixel)*JointMarginal(iPixel);
end

dFusionMaxEntropy = reshape(dFusionMaxEntropy,num_row,num_col);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function EstimatedPDF = KernelDenEst_LEN(x,t,h,lag)
[dim NumSample] = size(x);

EstimatedPDF=0;
for i=1:NumSample%scalar
   temp = KerNPDEScalar(x(i),t,h,lag)/(NumSample*h);
   EstimatedPDF=EstimatedPDF+temp;
   temp=0;
end
%*************************************************************    
    function yKer = KerNPDEScalar(Observation,t,BandWidth,lag)
    %calculate the likelihood of random varible to a observation
    %Observation:scalar
    %BandWidth:scalar
     switch lag
    
        case 1 %Gaussian kernel
             yKer = (1/((2*pi)^0.5))*exp(-0.5*((Observation-t)/BandWidth)^2);
         
         case 2
        
             
     end

     
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     function estCDF = funcEstCDF(Observations, NewObservations,varargin)
     
if length(varargin)==1
    FullObservationSpace = varargin{1};
    FullObservationSpace = sort(FullObservationSpace); 
else
    ObservationMin = min([min(Observations) min(NewObservations)]);
    ObservationMax = max([max(Observations) max(NewObservations)]);
    FullObservationSpace = linspace(ObservationMin,ObservationMax,10000);
end

for i = 1:length(NewObservations)
     estPFullObseSpace = KernelDenEst_LEN(Observations,NewObservations(i), 20,1); 
     FrequencyFullObseSpace = sum(estPFullObseSpace);
end

for i = 1:length(NewObservations)
    index = find(FullObservationSpace<=NewObservations(i)); 
    estPsubObservations = estPFullObseSpace(index);  
    FrequencySubObservations = sum(estPsubObservations);
    estCDF(i) = FrequencySubObservations/FrequencyFullObseSpace; 
    clear index;
    clear FrequencySubObservations;
    clear estPsubObservations;
end

function dRXcov = RXcov_detector(X)
% The Abnormaly Detector RX using Covariance Matrix
% Input: X -- The HyperCube 
% Output: D -- The quatratic statistics of the RX detector
[num_row,num_col,num_bands] = size(X);
for i = 1:num_bands
    Y(i,:) = reshape(X(:,:,i),1,num_row*num_col);
end
Y = Y';
CovMat = cov(Y); %spactral covariance matrix
InvCovMat = pinv(CovMat); %對協方差陣求逆
MeanY = mean(Y);
for i=1:num_row*num_col
 t=Y(i,:);   
 tempD(i)=(t-MeanY)*InvCovMat*(t-MeanY)';
end
dRXcov = reshape(tempD,num_row,num_col);