FCLS matlab源码 丰度估计

   function [alpha rmse] = fcls（totalpixelendmember）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  function: use FCLS algorithm to unmix the image
%%  totalpixel - pixel need to be unmixed
%%  endmember - the endmember matrix
%%  alpha - the abundance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear all
% [imagept]=freadenvi（‘D:/data/Fractal1_noise_30‘）;
% numofend = 9;
% [U11P11] = OSP（imagenumofend）;
% totalpixel=image;endmember=U11;

[rowlineb] = size（totalpixel）;
[band class] = size（endmember）;       % obtain the band number and class number
 
pp = zeros（band1）;
aa = zeros（class1）;
endm = endmember;

dita = 100;
pixel = dita*ones（band+11）;
alpha = zeros（rowlineclass）;
M = dita*ones（band+1class）;
M（1:b:） = endmember（1:b:）;
endmember = M;
for i=1:row
    for j=1:line
        m = endmember;
        pixel（1:band） = totalpixel（ij1:band）;
        flag = 1;
        [l p]=size（m）;      %l-波段数；p-端元光谱数
        order = 1:p;     %纪录求解出的系数的顺序
        counter = 0;            %纪录当前求出的系数
        while flag ~= 0
            [l p]=size（m）;      %l-波段数；p-端元光谱数
            one = ones（p1）;
            a = zeros（band1）;
            invmm=pinv（m‘*m）;    %inv（m*m‘）
            a_ls = invmm*m‘*pixel;      %使用无约束最小二乘求得系数
%             s = invmm*one;
%             lamda = （1-one‘*a_ls）/（one‘*s）;
%             a_fcls = a_ls-lamda*s;
            a_fcls = a_ls;
            re = （a_fcls>=0）;
            if sum（re）==p
                flag=0;
            else
%                 a_temp = a_fcls./s;
                a_temp = a_fcls;
                temp_max = 0;
                flag_zero = 0;
                %%寻找最大负系数
                for k=1:p
                    if temp_max                        temp_max = abs（a_temp（k））;
                        flag_zero = k;
                    end
                end
                counter = counter+1;
                order（flag_zero:p-1） = order（flag_zero+1:p）;
                %%删除最大负系数对应的项
                matrx_new = zeros（lp-1）;
                if flag_zero==1
                    matrx_new（:1:p-1） = m（:2:p）;
                else
                    if flag_zero==p
                        matrx_new（:1:p-1） = m（:1:p-1）;
                    else
                        matrx_new（:1:flag_zero-1） = m（:1:flag_zero-1）;
                        matrx_new（:flag_zero:p-1） = m（:flag_zero+1:p）;
                    end
                end
                m = matrx_new;
            end
        end
        if p==1
            invmm=inv（m‘*m）;    
            a_ls = invmm*m‘*pixel;
            alpha（ijorder（1）） = a_ls;
        else
            for k=1:p
                alpha（ijorder（k）） = a_fcls（k）;       %最终系数
            end
        end    
        
        pp（1:band） = totalpixel（ij1:band）;
        aa（1:class） = alpha（ij1:class）;