matlab开发-高光谱解混和去噪

matlab开发-高光谱解混和去噪。高光谱混合噪声解混演示。

function demo（）
%% This program is demo for unmixing. Please see the readme file for details.
%This program assume following noise model: 
%Y=MA+S+G  and solve the following problem
% min_{AS} ||Y-MA-S||_F^2 + lam1*||DhA‘||_21 + lam1*||DvA‘||_1
% +lam2||S||_1 + lam3||A||_{21}

% M: mixing matrix; 
% Y: Noisy Image ; 
% A: Abundance matrix
% S : sparse noise; Dh and Dv are Total variation operators.

% m: number of rows in one band.
% n: number of columns in one band.
% b: number of bands.
% e: total number of endmemebrs available.
% k: number of endmembers used to make the image.


% Refrence paper :
% Hyperspectral Unmixing in the Presence of Mixed Noise using Joint-Sparsity and Total-Variation 
% Hemant Kumar Aggarwal Angshul Majumdar
% IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing （JSTARS） 2016
% link : http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=7414394
% paper PDF is available from author‘s homepage as well.

%% create HS image using abundance maps & endmember matrix
%clc; clear all; close all;

M=rand（5070）; %take some random mixing matrix
M=normc（M）;     % normalize the each endmember spectrum
[be]=size（M）;        
load abundance  %Ground truth abundance maps
[mnk]=size（abundance）;  
Atemp=reshape（abundancem*nk）‘;
%idx=randperm（ek）;
idx=[6    41    11     5    39];
A=zeros（em*n）;
A（idx:）=Atemp;
X=M*A;
img=reshape（X‘mnb）; %
img2=myhisteq（img）; %just for display purpose
%% Add Gaussian noise + impulse noise + line strips
snr=50;  %add Gaussian noise of this much SNR
noisy=addGaussianNoise（imgsnr）; 
noisy2=addGaussianNoise（img2snr）; %noisy2 is for display purpose only

%add vertical lines
    noisy（ :[10 24  46][1  5]）=min（img（:））;
    noisy2（ :[10 24  46][1  5]）=min（img2（:））;
 
%add  some impulse noise
noisy=imnoise（noisy‘salt & pepper‘.01）;
noisy2=imnoise（noisy2‘salt & pepper‘.01）;

% Check PSNR of noisy image
noisyPSNR=myPSNR（imgnoisy1）;
fprintf（‘PSNR of noisy image= %f \n ‘noisyPSNR）;

X=reshape（noisym*nb）‘;

%% Set parameters and Run the algorithm

  opts.m=m;
  opts.n=n;
 % These  parameters can be controlled as required.
 opts.lambda1=.5; % for total variation
 opts.lambda2=.5; % for sparse term
 opts.lambda3=1; % for joint sparsity i.e. L21 term
     opts.mu1=.35; % for TV regularization term
     opts.mu2=10;  % for L21 regularization term
  

 abun=abundance;     
 opts.iter=50;
 tic;
 Ahat=funJSTV（MXopts）;
 toc;
%% Calculate PSNR of denoised image and abundance maps 
%get image by multiplying with recovered abundance maps in variable Ahat 
Xhat=M*Ahat; 
 rec=reshape（Xhat‘mnb）;
Ahat1=reshape（Ahat（idx:）‘mnnumel（idx））;

recAbunPSNR=myPSNR（abunAhat11）;
recPSNR=myPSNR（imgrec1）;

fprintf（‘\nReconstructed Abundance PSNR=%f \n‘recAbunPSNR）;
fprintf（‘Reconstructed Image PSNR= %f \n ‘recPSNR）;

%% showing abundance map
%figure;
  % These are five original abundance maps. 
figure（‘position‘ [100 100 700 500]） 
 h1=subplot（251）;    imagesc（abun

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        2512  2016-06-26 07:14  UnmixingCode\ReadMe.txt
     文件        1335  2016-06-26 07:14  UnmixingCode\abundance.mat
     文件        8267  2016-06-26 07:14  UnmixingCode\demo.m
     文件        1521  2016-06-26 07:14  license.txt