遥感影像变化检测算法(IR-MAD、MAD、CVA、PCA)

clc clear %% load the image data % imt1 = imread('04.bmp'); % imt2 = imread('05.bmp'); % [X,p]=freadenvi('data\\Whu\\2002roi',1); % [Y,p2]=freadenvi('data\\Whu\\2005roi',1); [X,p]=freadenvi('data\\Taizhou\\2000TM',1); [Y,p2]=freadenvi('data\\Taizhou\\2003TM',1); % [X,p]=freadenvi('data\\Nanjing\\LT51200382000124BJC00_Nanjing',1); % [Y,p2]=freadenvi('data\\Nanjing\\LT51200382002193BJC00_Nanjing',1); rows=p(1); cols=p(2); bands=p(3); N =size(X,1); %行数 rows*cols X=double(X); Y=double(Y); %滤波 h2 = fspecial('gaussian',5,1); %高斯低通滤波 hsize表示模板尺寸，默认值为【3 3】，sigma为滤波器的标准值，单位为像素，默认值为0.5.这里是1 X = imfilter(X, h2, 'symmetric'); %滤波器 Y = imfilter(Y, h2, 'symmetric'); lastrou=0; xpixs=[];%记录迭代次数信息 rouvalues=[];%记录CCA最大特征值信息 allrouvalues=[]; iteraNums=50;% 10 30 50 ||70 90 110 %初始权重 weight = ones(N,1); % originalWeight=reshape(weight,rows,cols,1); % imwrite(originalWeight','Test result\\NanjingChange_110_iterations_ChangeInageWeight_original.tif'); % enviwrite2(originalWeight,'Change Result\\NanjingChange_110_iterations_ChangeInage_weight_original'); for i=1:iteraNums disp(['--------------',num2str(i),'次迭代-------------']); % 中心化 % X=zscore(X); % Y=zscore(Y); % X=X-repmat(mean(X),N,1); % Y=Y-repmat(mean(Y),N,1); % disp(size(X)); %% Calculate Cov 11 12 21 22 % covxy = cov([X Y]); %8*8 [covxy,meanw] = covw([X Y],weight); meanwX = meanw(1:bands); meanwY = meanw(bands+1:end); %disp(covxy); cov11 = covxy(1:bands,1:bands); cov22 = covxy(bands+1:end,bands+1:end); cov12 = covxy(1:bands,bands+1:end); cov21 = cov12'; %% Calculate a and b invcov22=inv(cov22); [v,d]=eig(cov12*invcov22*cov21,cov11); %式32 v 特征向量表示a d 特征值表示rou方 d=diag(d); [~,idx] = sort(d); %idx 索引位置 v = v(:,idx); %参数a d = diag(d(idx)); %rou2 rou=diag(sqrt(d))';%求出bands个rou allrouvalues=[allrouvalues;rou]; disp(num2str(rou,' %0.6g')); nowmaxrou=max(rou); %终止条件 两次最大特征值相差10-6 % disp('----------最大特征值差值--------'); % disp(num2str(nowmaxrou-lastrou,'%0.6g')); xpixs=[xpixs,i]; if(abs(nowmaxrou-lastrou)<10^(-6)) break; end; lastrou=nowmaxrou; %满足条件 a'*cov11*a=1 aux1 = v'*cov11*v; aux2 = 1./sqrt(diag(aux1)); aux3 = repmat(aux2',bands,1); v1 = v.*aux3; % invstderr = diag(1./std(X)); sgn = diag(sign(sum(invstderr*cov11*v1))); v1 = v1*sgn; %test b v2 = invcov22*cov21*v1; %./repmat(rho(end,:),nvar1,1); % scaling doesn't matter % scale v2 to give CVs with unit variance aux1 = v2'*cov22*v2; % dispersion of CVs aux2 = 1./sqrt(diag(aux1)); aux3 = repmat(aux2',bands,1); v2 = v2.*aux3; % now dispersion is unit matrix %MAD=X*v1-Y*v2; MAD = (X-repmat(meanwX,N,1))*v1 - (Y-repmat(meanwY,N,1))*v2; % save('MAD','MAD'); %% MAD % rou=fliplr(rou); %公式11 倒置波段 sigma2=2*(1-rou); tj = sum((MAD./repmat(sqrt(sigma2),N,1)).^2,2); %公式12 weight=1-gammainc(0.5*tj,0.5*bands); % w'=1-a gammainc(x,a)=x^a 论文2 公式6 % % % % if(i==1) % % firstWeight=reshape(weight,rows,cols,1); % % imwrite(firstWeight','TaizhouChange_50_iterations_ChangeInageWeight_first.tif'); % % enviwrite2(firstWeight,'TaizhouChange_50_iterations_ChangeInage_weight_first'); % % end; % save(['MAD\',num2str(i)],'MAD'); % save(['Weight\',num2str(i)],'weight'); % a=i; % xpixs=[xpixs,i]; % rouvalues=[rouvalues,nowmaxrou]; end; % createfigure(xpixs, allrouvalues); 显示迭代过程中各个波段特征值的变化趋势 save('LastWeight','weight'); % %权重影像 % weight=1-weight; % weight=reshape(weight,rows,cols,1); % imwrite(weight,'TaizhouChange_50_iterations_ChangeInageWeight_last.tif'); % enviwrite2(weight,'TaizhouChange_50_iterations_ChangeInage_weight_last'); %IR-MAD特征差值影像 DifferMAD=reshape(MAD,rows,cols,bands); % enviwrite(DifferMAD,'Change Result\\TaizhouChange_test_iterations'); % enviwrite(MAD,'Update_TaizhouChange_testBandsRou_iterations'); %转换为变化强度图 方法有欧式距离、卡方距离 % 1. 欧式距离 2.卡方距离 EuclideanMAD=sqrt(sum(MAD.*MAD,2)); %像素在各个波段统一变化强度 %k-means 二值聚类 CM=KmeansMap(EuclideanMAD,rows,cols); ChangeMap=CM; ChangeMap(find(ChangeMap==1))=255; % disp(size(ChangeMap)); % imwrite(ChangeMap,'Test result\\Taizhou_Kmeans_test_TwoValue DifferImage.bmp'); EuclideanMAD=reshape(EuclideanMAD,rows,cols,1); %c存储变化强度图 % enviwrite2(EuclideanMAD,'Change Result\\TaihzoutestChange_Intensity_Image'); %% Accuracy %load Change Image Ref % ChangeRI_Nanjing=GetChangeOfNanjing(rows,cols); %400*400 ChangeRI=imread('data\\Taizhou\\TaizhouChange_blackWhite.bmp'); ChangeRI=ChangeRI(:,:); %load unChange Image Ref UnchangeRI=imread('data\\Taizhou\\TaizhouChange_blackWhite_unchange.bmp'); UnchangeRI=UnchangeRI(:,:); AllRef=zeros(rows,cols); for i=1:N if(UnchangeRI(i)==255) AllRef(i)=1; end; if(ChangeRI(i)==255) AllRef(i)=2; end; end %ROC AUC disp(sprintf('检测方法指标评价 IR-MAD Kappa And Accuracy:')); [IRMAD_test,IRMAD_lastRef,IRMAD_Accuracy,IRMAD_Kappa,IRMAD_canshu,IRMAD_AUC,IRMAD_Result]=DAcom(AllRef,CM,0); % % ------------CVA------------ disp(sprintf('检测方法指标评价 CVA Kappa And Accuracy:')); CVA_CM=CVADemo(); [CVA_test,CVA_lastRef,CVA_Accuracy,CVA_Kappa,CVA_canshu,CVA_AUC,CVA_Result]=DAcom(AllRef,CVA_CM,1) ; %------------PCA------------- disp(sprintf('检测方法指标评价 PCA Kappa And Accuracy:')); PCA_CM=PCADemo(); [PCA_test,PCA_lastRef,PCA_Accuracy,PCA_Kappa,PCA_canshu,PCA_AUC,PCA_Result]=DAcom(AllRef,PCA_CM,1) ; %------------MAD------------- disp(sprintf('检测方法指标评价 MAD Kappa And Accuracy:')); MADCM=MADGet(); [MAD_test,MAD_lastRef,MAD_Accuracy,MAD_Kappa,MAD_canshu,MAD_AUC,MAD_Result]=DAcom(AllRef,MADCM,0) ; %% ROC AUC [val,ind] = sort(CVA_test,'descend'); roc_y = CVA_lastRef(ind); stack_x = cumsum(roc_y == 0)/sum(roc_y == 0); stack_y = cumsum(roc_y == 1)/sum(roc_y == 1); CVA_auc = sum((stack_x(2:length(roc_y),1)-stack_x(1:length(roc_y)-1,1)).*stack_y(2:length(roc_y),1)); [val,ind] = sort(PCA_test,'descend'); roc_y = PCA_lastRef(ind); stack_x1 = cumsum(roc_y == 0)/sum(roc_y == 0); stack_y1 = cumsum(roc_y == 1)/sum(roc_y == 1); PCA_auc = sum((stack_x1(2:length(roc_y),1)-stack_x1(1:length(roc_y)-1,1)).*stack_y1(2:length(roc_y),1)); [val,ind] = sort(MAD_test,'descend'); roc_y = MAD_lastRef(ind); stack_x2 = cumsum(roc_y == 0)/sum(roc_y == 0); stack_y2 = cumsum(roc_y == 1)/sum(roc_y == 1); MAD_auc = sum((stack_x2(2:length(roc_y),1)-stack_x2(1:length(roc_y)-1,1)).*stack_y2(2:length(roc_y),1)); [val,ind] = sort(IRMAD_test,'descend'); roc_y = IRMAD_lastRef(ind); stack_x3 = cumsum(roc_y == 0)/sum(roc_y == 0); stack_y3 = cumsum(roc_y == 1)/sum(roc_y == 1); IRMAD_auc = sum((stack_x3(2:length(roc_y),1)-stack_x3(1:length(roc_y)-1,1)).*stack_y3(2:length(roc_y),1)); plot(stack_x,stack_y,stack_x1,stack_y1,stack_x2,stack_y2,stack_x3,stack_y3); xlabel('False Positive Rate'); ylabel('True Positive Rate'); title(['ROC curve of (IRMAD-AUC = ' num2str(IRMAD_auc) '、MAD-AUC =' num2str(MAD_auc) '、CVA-AUC = ' num2str(CVA_auc) '、PCA-AUC = ' num2str(PCA_auc) ')' ]);