%function coeffs = decdemo( im, option )
% DECDEMO demonstrates nonsubsampled Contourlet decomposition and reconstruction.
disp('Welcome to the nonsubsampled Contourlet decomposition demo! :)');
disp('Type help decdemo for help' ) ;
disp('You can also view decdemo.m for details.') ;
disp(' ');
% Input image
if ~exist('im', 'var')
% Zoneplate image: good for illustrating multiscale and directional
% decomposition
im = imread ('sar.png') ;
elseif isstr(im)
im = imread ( im ) ;
else
error('You shall input valid image name!');
end
% % % Show the input image
% % disp( 'Displaying the input image...');
% % clf;
% % imagesc(im, [0, 255]);
% % title( 'Input image' ) ;
% % axis image off;
% % colormap(gray);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Image decomposition by nonsubsampled contourlet transform (NSCT).
% This is the iterated filter bank that computes the nonsubsampled
% contourlet transform.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 分解参数设置
nlevels = [1, 1, 1] ; % Decomposition level
pfilter = 'maxflat' ; % Pyramidal filter
dfilter = 'dmaxflat7' ; % Directional filter
% nsct分解
%coeffs = nsctdec( double(im), nlevels, dfilter, pfilter );
y = nsctdec( double(im), nlevels, dfilter, pfilter );
% Display the coefficients
disp('Displaying the contourlet coefficients...') ;
%显示各层系数
%shownsct( coeffs ) ;
% Level of decomposition.
clevels = length( y ) ;
% Show the subband images.
% % for i=1:clevels
% % figure;
% % if iscell( y{i} ) %如果是单元数组
% % % The number of directional subbands.
% % csubband = length( y{i} ) ; %数组个数
% % if csubband > 7
% % col = 4 ;
% % else
% % col = 2 ;
% % end
% % row = csubband / col ;
% % for j = 1:csubband
% % subplot( row, col, j ) ;
% % imshow( uint8(y{i}{j}) ); %显示各层各方向系数图片
% % title( sprintf('NSSC coefficients: level %d, %d', i,length(y{i}{j})) );
% % end
% % else
% % imshow ( uint8(y{i}) ) ; %显示不是单元数组的低通系数图片
% % title( sprintf('Nonsubsampled Contourlet coefficients level %d, %d', i,length(y{i})) );
% % end
% % end
%显示完毕
%系数处理
%计算所有系数层标准差%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for yi=2:4;
for yj=1:2;
in{yi}{yj}=double(zeros(256,256)) ;
in1{yi}{yj}=double(zeros(256,256)) ;
in{yi}{yj}=y{yi}{yj};
zs=median(in{yi}{yj}(:));
for zi=1:256;
for zj=1:256;
in1{yi}{yj}(zi,zj)=abs(in{yi}{yj}(zi,zj)-zs) ;
end
end
zs2{yi}{yj}=median(in1{yi}{yj}(:))/0.6745; %标准差
end
end
%标准差计算完毕
%假设第2层为无噪系数层，依次估计第3层掩膜%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ym31=double(zeros(256,256)) ;
ym32=double(zeros(256,256)) ;
yc31=double(zeros(256,256)) ;
yc32=double(zeros(256,256)) ;
for ymi=1:256
for ymj=1:256
yc31=abs(y{2}{1}(ymi,ymj)*y{3}{1}(ymi,ymj)); %相邻层乘积绝对值
yc32=abs(y{2}{2}(ymi,ymj)*y{3}{2}(ymi,ymj));
if yc31 < zs2{3}{1}^2
ym31(ymi,ymj)=0;
else
ym31(ymi,ymj)=1;
end
if yc32 < zs2{3}{2}^2
ym32(ymi,ymj)=0;
else
ym32(ymi,ymj)=1;
end
end
end
%掩膜估计完毕
%计算第3层缩减因子%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算似然比
%xs310=double(zeros(256,256)) ; % 掩码为0的3层1区样本集合
%xs311=double(zeros(256,256)) ;
%xs320=double(zeros(256,256)) ;
%xs321=double(zeros(256,256)) ;
pr311=double(zeros(256,256)) ; %1似然概率
pr321=double(zeros(256,256)) ;
pr310=double(zeros(256,256)) ; %0似然概率
pr320=double(zeros(256,256)) ;
sirb31=double(zeros(256,256)) ;
sirb32=double(zeros(256,256)) ;
%xs310=y{3}{1};
xx1=1;
xx2=1;
xx3=1;
xx4=1;
for xsi=1:256
for xsj=1:256
if ym31(xsi,xsj)==1 %掩码为1的3层1区样本集合
xs311(xx1)=abs(y{3}{1}(xsi,xsj));
xx1=xx1+1;
% else
% xs311(xsi,xsj)=0;
end
if ym31(xsi,xsj)==0 %掩码为0的3层1区样本集合
xs310(xx2)=abs(y{3}{1}(xsi,xsj));
xx2=xx2+1;
% else
% xs310(xsi,xsj)=0;
end
if ym32(xsi,xsj)==1 %掩码为1的3层2区样本集合
xs321(xx3)=abs(y{3}{2}(xsi,xsj));
xx3=xx3+1;
% else
% xs321(xsi,xsj)=0;
end
if ym32(xsi,xsj)==0 %掩码为0的3层2区样本集合
xs320(xx4)=abs(y{3}{2}(xsi,xsj));
xx4=xx4+1;
% else
% xs320(xsi,xsj)=0;
end
end
end
junz310=mean2(xs310); %掩膜为0或者为1的系数均值
junz311=mean2(xs311);
junz320=mean2(xs320);
junz321=mean2(xs321);
biaofc310=std2(xs310); %掩膜为0或者为1的系数方差
biaofc311=std2(xs311);
biaofc320=std2(xs320);
biaofc321=std2(xs321);
xingz311=(junz311/biaofc311)^2 ; %形状参数
chidu311=junz311/(biaofc311)^2 ; %尺度参数
xingz321=(junz321/biaofc321)^2 ; %形状参数
chidu321=junz321/(biaofc321)^2 ; %尺度参数
zhis320=1/junz320; %指数分布参数
zhis310=1/junz310;
ga311=gamma(xingz311) ;
ga321=gamma(xingz321) ;
for pi=1:256
for pj=1:256
%1似然分布
pr311(pi,pj)=( (chidu311^xingz311)/ga311 ) * ( ( abs(y{3}{1}(pi,pj)) )^(xingz311-1) ) * exp( (-1)*chidu311 * abs(y{3}{1}(pi,pj)) );
pr321(pi,pj)=( (chidu321^xingz321)/ga321 ) * ( ( abs(y{3}{2}(pi,pj)) )^(xingz321-1) ) * exp( (-1)*chidu321 * abs(y{3}{2}(pi,pj)) );
%0似然分布
pr310(pi,pj)=zhis310*exp(-zhis310 * abs(y{3}{1}(pi,pj)) ) ;
pr320(pi,pj)=zhis320*exp(-zhis320 * abs(y{3}{2}(pi,pj)) ) ;
sirb31(pi,pj)=pr311(pi,pj)/pr310(pi,pj); %似然比
sirb32(pi,pj)=pr321(pi,pj)/pr320(pi,pj);
end
end
%第3层似然比计算完毕
%计算第3层先验比%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xy31=double(zeros(256,256)) ; %定义先验系数集合
xy32=double(zeros(256,256)) ;
xy31(1,:)=exp(0.4*8);
xy31(256,:)=exp(0.4*8);
xy31(:,1)=exp(0.4*8);
xy31(:,256)=exp(0.4*8);
xy32(1,:)=exp(0.4*8);
xy32(256,:)=exp(0.4*8);
xy32(:,1)=exp(0.4*8);
xy32(:,256)=exp(0.4*8);
contr=0.4 %缩减控制因子
for xyi=2:255
for xyj=2:255
sum2=(2*ym31(xyi,xyj+1)-1) + (2*ym31(xyi-1,xyj+1)-1) +(2*ym31(xyi-1,xyj)-1)+(2*ym31(xyi-1,xyj-1)-1)+(2*ym31(xyi,xyj-1)-1) + (2*ym31(xyi+1,xyj-1)-1) + (2*ym31(xyi+1,xyj)-1) + (2*ym31(xyi+1,xyj+1)-1) ; %系数周围掩码和
if ym31(xyi,xyj)==1&&ym31(xyi,xyj-1)==1&&ym31(xyi,xyj+1)==1 %&&ym31(xyi-1,xyj-1)==0&&ym31(xyi-1,xyj)==0&&ym31(xyi-1,xyj+1)==0 ) || ( ym31(xyi,xyj)==1&&ym31(xyi,xyj-1)==1&&ym31(xyi,xyj+1)==1&&ym31(xyi+1,xyj-1)==0&&ym31(xyi+1,xyj)==0&&ym31(xyi+1,xyj+1)==0 )
xy31(xyi,xyj)=exp(32) ;
elseif ym31(xyi,xyj)==1&&ym31(xyi-1,xyj)==1&&ym31(xyi+1,xyj)==1 %&&ym31(xyi,xyj-1)==0&&ym31(xyi-1,xyj-1)==0&&ym31(xyi+1,xyj-1)==0 ) || (ym31(xyi,xyj)==1&&ym31(xyi-1,xyj)==1&&ym31(xyi+1,xyj)==1&&ym31(xyi,xyj+1)==0&&ym31(xyi-1,xyj+1)==0&&ym31(xyi+1,xyj+1)==0 )
xy31(xyi,xyj)=exp(32) ;
elseif ym31(xyi,xyj)==1&&ym31(xyi+1,xyj+1)==1&&ym31(xyi-1,xyj-1)==1 %&&ym31(xyi,xyj+1)==0&&ym31(xyi-1,xyj)==0&&ym31(xyi-1,xyj+1)==0 ) || (ym31(xyi,xyj)==1&&ym31(xyi+1,xyj+1)==1&&ym31(xyi-1,xyj-1)==1&&ym31(xyi,xyj-1)==0&&ym31(xyi+1,xyj)==0&&ym31(xyi+1,xyj-1)==0 )
xy31(xyi,xyj)=exp(32) ;
elseif ym31(xyi,xyj)==1
