彩色图像空间转换代码

1、图像篡改检测Matlab源代码说明（1） 彩色图像空间转换代码（RGB空间转换为YUV空间）clc;clear all;close all;img = imread(3.jpg);%figure;imshow(img);title(原图);img_5=img;img_ycbcr = rgb2ycbcr(img); % rgb-yuvfigure;subplot(121);imshow(img);title(原始图像);subplot(122);imshow(img_ycbcr);title(YUV空间图);图1 彩色图像空间转换（2）图像分块（分为16*16子块）row,col,i=size(

2、img_ycbcr);%对图像进行扩展row_expand=ceil(row/16)*16; %行数上取整再乘16，及扩展成16的倍数if mod(row,16)=0 %行数不是16的倍数，用最后一行进行扩展 for i=row:row_expand img_ycbcr(i,:,:)=img_ycbcr(row,:,:); endendcol_expand=ceil(col/16)*16; %列数上取整if mod(col,16)=0 %列数不是16的倍数，用最后一列进行扩展 for j=col:col_expand img_ycbcr(:,j,:)=img_ycbcr(:,col,:); e

3、ndend（3）对各个分量（Y,Cr，Cb）进行离散余弦变换量化%对各分量进行4:2:0采样Y=img_ycbcr(:,:,1); %Y分量figure;subplot(231);imshow(Y);title(原Y分量);Cb=img_ycbcr(:,:,2); %Cb分量Cr=img_ycbcr(:,:,3); %Cr分量 Cb = zeros(row_expand/2,col_expand/2); Cr = zeros(row_expand/2,col_expand/2);for i=1:row_expand/2 for j=1:2:col_expand/2-1 %奇数 Cb(i,j)=

4、double(img_ycbcr(i*2-1,j*2-1,2); Cr(i,j)=double(img_ycbcr(i*2-1,j*2+1,3); endendfor i=1:row_expand/2 %偶数 for j=2:2:col_expand/2 Cb(i,j)=double(img_ycbcr(i*2-1,j*2-2,2); Cr(i,j)=double(img_ycbcr(i*2-1,j*2,3); endendsubplot(232);imshow(uint8(Cb);title(原Cb分量); %Cb分量subplot(233);imshow(uint8(Cr);title(原

5、Cr分量); %Cr分量Y_Table=16 11 10 16 24 40 51 61 12 12 14 19 26 58 60 55 14 13 16 24 40 57 69 56 14 17 22 29 51 87 80 62 18 22 37 56 68 109 103 77 24 35 55 64 81 104 113 92 49 64 78 87 103 121 120 101 72 92 95 98 112 100 103 99;%亮度量化表 CbCr_Table=17, 18, 24, 47, 99, 99, 99, 99; 18, 21, 26, 66, 99, 99, 99,

6、 99; 24, 26, 56, 99, 99, 99, 99, 99; 47, 66, 99 ,99, 99, 99, 99, 99; 99, 99, 99, 99, 99, 99, 99, 99; 99, 99, 99, 99, 99, 99, 99, 99; 99, 99, 99, 99, 99, 99, 99, 99; 99, 99, 99, 99, 99, 99, 99, 99;%色差量化表Qua_Factor=0.2;%量化因子%对Y分量进行处理Qua_Matrix=Qua_Factor*Y_Table; %Y量化矩阵Y = blkproc(Y,8 8,dct2(x); %DCT变

7、换Y = blkproc(Y,8 8,round(x./P1),Qua_Matrix); %量化Y = blkproc(Y,8 8,x.*P1,Qua_Matrix);%反量化Y= blkproc(Y,8 8,idct2(x);%反DCT变换subplot(234);imshow(uint8(Y);title(Y分量量化后);%对Cb分量进行处理Qua_Matrixcc=0.2*CbCr_Table; % Cb质量因子选取为5，并不采用Qua_Factor的数值Cb = blkproc(Cb,8 8,dct2(x); %DCT变换Cb = blkproc(Cb,8 8,round(x./P1)

8、,Qua_Matrixcc); %量化Cb = blkproc(Cb,8 8,x.*P1,Qua_Matrixcc);%反量化Cb = blkproc(Cb,8 8,idct2(x);%反DCT变换subplot(235);imshow(uint8(Cb);title(Cb分量量化后); %对Cr分量进行处理Qua_Matrixcc=0.2*CbCr_Table; %Cr质量因子选取为5，并不采用Qua_Factor的数值Cr = blkproc(Cr,8 8,dct2(x); %DCT变换Cr = blkproc(Cr,8 8,round(x./P1),Qua_Matrixcc); %量化C

9、r = blkproc(Cr,8 8,x.*P1,Qua_Matrixcc);%反量化Cr = blkproc(Cr,8 8,idct2(x);%反DCT变换subplot(236);imshow(uint8(Cr);title(Cr分量量化后);图2 各个分量量化图像（4）对各分量（Y,Cr，Cb）进行二值化处理 %对Y分量处理%差分处理dh_Y = diff(Y,2);%水平差分dvt_Y = diff(Y); %纵向差分dv_Y= dvt_Y; %显示功率谱密度fs=1000;n=0:1/fs:1;nfft=1024;window=blackman(100);noverlap=20;ra

10、nge=half;pxx1,f=pwelch(dh_Y,window,noverlap,nfft,fs,range);pxx2,f=pwelch(dvt_Y,window,noverlap,nfft,fs,range);plot_pxx1=10*log10(pxx1);plot_pxx2=10*log10(pxx2); figure,subplot(121);plot(f,plot_pxx1);title(水平功率谱密度);subplot(122);plot(f,plot_pxx2);title(垂直功率谱密度); %计算功率谱密度dh1_Y= blkproc(dh_Y,8 8,welch);

11、 %水平方向差分图像的区域功率谱密度dv1_Y= blkproc(dv_Y,8 8,welch); %垂直方向差分图像的区域功率谱密度b1_Y= blkproc(dh1_Y,513 1,sum(x(:);%水平方向差分图像的区域功率谱密度总和b2_Y= blkproc(dv1_Y,513 1,sum(x(:);%垂直方向差分图像的区域功率谱密度总和b_Y=b1_Y+b2_Y;%整幅图像的区域功率谱密度总和 for i=1:size(b_Y,1) for j=1:size(b_Y,2) B_Y(i,j)=exp(log(sum(b_Y(i,j)+1)/(sum(b_Y(i,j)+1);%整幅图像

12、块效应评价 endend %二值化处理T=exp(1.5)*sum(sum(B_Y)/(row_expand/8)*(col_expand/8);%阈值c_Y=zeros(size(Y);for i=1:size(Y,1) for j=1:size(Y,2) if b_Y(ceil(i/8),ceil(j/8)T c_Y(i,j)=1; end endendfigure;subplot(121);imshow(int8(b_Y); title(Y分量离散余弦系数图);subplot(122);imshow(c_Y); title(Y分量二值图像);图3 Y分量离散余弦变换二值化图像图4 Cb分

13、量离散余弦变换二值化图像图5 Cr分量离散余弦变换二值化图像(5) 提取各个子图特征function featureVec = GetBlockFeature(Block)featureVec=zeros(1,7); %特征数组，共设置7个特征R = Block(:,:, 1);G = Block(:,:, 2);B = Block(:,:, 3); %cl,c2, c3三个特征值featureVec(1)=mean2(R(:);featureVec(2)=mean2(G(:);featureVec(3)=mean2(B(:); %c4.5.6.7四个特征值,Y通道公式Y = im2doubl

14、e(0.299 * R + 0.587 * G + 0.114 * B); %块内所有像素灰度之和totalGray = sum(Y(:); %c4sum_partI= sum(sum(Y(1 : 8,:);featureVec(4) = sum_partI / totalGray; %c5sum_partI = sum(sum(Y(:,1 : 8);featureVec(5) = sum_partI/ totalGray; %c6sumPart11 = 0.0;for r = 1 : 16 for c = r : 16 sumPart1 = sumPart11 + Y(r, c); ende

15、ndfeatureVec(6)=sumPart1/totalGray; %c7sumPart11 = 0.0;for r = 1 : 16 for c = 1 : (16 - r + 1) sumPart1 = sumPart11 + Y(r, c); endendfeatureVec(7)= sumPart1 / totalGray;return;(6) 特征排序FeatureVecThresh =2.5,1.5,3.0,0.006,0.005,0.005,0.005;BlockDistThresh = 10;BlockPairsNum = NumBlocks - 1 ;ShiftVecto

16、rs = zeros(1,2);candidateSimilarPairs.indices = zeros(1,2);candidateSimilarPairs.frequence(1) = 0;features, indices = sortrows(Blocks.features);vi = 0;for i = 1 : NumBlocks - 1 ind1 = indices(i); ind2 = indices(i + 1); deltaCoord = Blocks.topLeft(ind1,:)-Blocks.topLeft(ind2,:); deltaFeature = Blocks

17、.features(ind1,:)-Blocks.features(ind2,:); pixelDist = norm(deltaCoord); %向景取模 if (pixelDistBlockDistThresh) & (sum(abs(deltaFeature)FeatureVecThresh) = 7) vi = vi + 1; candidateSimilarPairs.indices(vi, :)=ind1 ind2; ShiftVectors(vi,:) = Blocks.topLeft(ind1,:)-Blocks.topLeft(ind2,:); endend(7) 查找并删除

18、重复图像块SimilarPairs.indices=zeros(1,2);pairs.vector1=ShiftVectors(1,:); %第一行pairs.freq(1)=1;ind=1;for i=2:size(ShiftVectors,1); vec=ShiftVectors(i,:); isExist=false; for j=1:ind if norm(vec-pairs.vectorj)2 pairs.freq(j)=pairs.freq(j)+1; isExist=true; break; end end if isExist ind=ind+1; pairs.vectorind=vec; pairs.freq(ind)=1; end