基于奇异值分解的图像压缩matlab

话不多说上才艺

大家每天都在使用的jpg，png，gif其实都是压缩过的图片格式，他们在底层的数据并不是挨个记录每个像素的rgb值的，而真正这么干的格式，大家比较熟悉的大概就只有bmp格式了，所以一般bmp格式的图片都比较大。

嘛是图像压缩‍♂️

目前图像压缩的方法主要分为两大类：无损压缩和有损压缩。无损压缩利用数据的统计冗余进行压缩，可完全恢复原始数据而不引入任何失真，但无损压缩受图像本身信息熵的制约，压缩比受限，但由于压缩比的限制，仅使用无损压缩的方法已经不能解决海量图像的存储和传输问题。有损压缩方法利用了人类视觉对图像中的某些频率成分不敏感的特性，允许压缩过程中损失一定的信息：虽然不能完全恢复原始数据，但是采用适当的方法可以使所损失的部分对原始图像的影响较小，进而得到较大的压缩比。

• 图像数据的冗余主要表现为：

1. 图像中相邻像素间的相关性引起的空间冗余。
2. 图像序列中不同帧之间存在相关性引起的时间冗余。
3. 不同彩色平面或频谱带的相关性引起的频谱冗余。

聊聊图像格式‍

以常用的PNG 和 JPG 格式的图像为例,读取它们的图像数据可得到一个 3 维的数据矩阵.第1维表示的横向的像素,第2维表示纵向的像素,第3维表示图像的通道,用 0~255 范围的数据表示数据值.PNG 和 JPG 两种格式不同的是, PNG 格式图像的第3维有4个通道,分别表示R(Red,红色)、G(Green,绿色)、B(Blue,蓝色)、A(Alpha,透明度);而JPG格式图像只有 R、G、B 这 3 个通道,没有A这个通道。这也是JPG格式图像比 PNG 格式图像占用空间更小的根本原因

• 一张rgb彩色通道的图像矩阵：

奇异值分解

特征值分解evd

在这里，选择一种特殊的矩阵——对称阵（酉空间中叫hermite矩阵）。对称阵有一个性质：它总能相似对角化，对称阵不同特征值对应的特征向量两两正交。一个矩阵能相似对角化即说明其特征子空间即为其列空间，若不能对角化则其特征子空间为列空间的子空间。现在假设存在mxm的满秩对称矩阵A

a 1 x 1 = λ 1 x 1 a 2 x 2 = λ 2 x 2 ⋮ a n x n = λ n x n A U = U Σ a_1x_1=\\lambda_1x_1\\\\ a_2x_2=\\lambda_2x_2\\\\ \\vdots\\\\ a_nx_n=\\lambda_nx_n\\\\ AU=U\\Sigma\\\\ a1​x1​=λ1​x1​a2​x2​=λ2​x2​⋮an​xn​=λn​xn​AU=UΣ
其中
U = [ x 1 , x 2 , ⋯   , x n ] Σ = [ λ 1 0 ⋯ 0 0 λ 2 ⋯ 0 ⋮ λ n 0 0 ⋯ 0 ] t h e n A = U Σ U − 1 U=[x_1,x_2,\\cdots,x_n]\\\\ \\Sigma=\\left[ \\begin{matrix} \\lambda_1&0&\\cdots&0\\\\ 0&\\lambda_2&\\cdots&0\\\\ &\\vdots&\\lambda_n\\\\ 0&0&\\cdots&0 \\end{matrix} \\right]\\\\ then\\\\ A=U\\Sigma U^{-1} U=[x1​,x2​,⋯,xn​]Σ=⎣⎢⎢⎢⎡​λ1​00​0λ2​⋮0​⋯⋯λn​⋯​000​⎦⎥⎥⎥⎤​thenA=UΣU−1
因为正交矩阵转置就是他的逆： A = U Σ U T A=U\\Sigma U^T A=UΣUT

奇异值分解svd

特征值分解是一个提取矩阵特征的好方法，但是它只是对方阵而言的，对于非方阵描述其特征就要用到奇异值分解了。想详细了解请观看矩阵论中奇异值分解那章。

svd(Singular Value Decomposition)：
M = U Σ V T [ x 1 y 1 x 2 y 2 x 3 y 3 ] = [ − − − − − − − − − ] [ a 0 0 b 0 0 ] [ − − − − ] M=U\\Sigma V^T\\\\ \\left[ \\begin{matrix} x_1&y_1\\\\ x_2&y_2\\\\ x_3&y_3\\\\ \\end{matrix} \\right]= \\left[ \\begin{matrix} -&-&-\\\\ -&-&-\\\\ -&-&-\\\\ \\end{matrix} \\right] \\left[ \\begin{matrix} a&0\\\\ 0&b\\\\ 0&0\\\\ \\end{matrix} \\right] \\left[ \\begin{matrix} -&-\\\\ -&-\\\\ \\end{matrix} \\right] M=UΣVT⎣⎡​x1​x2​x3​​y1​y2​y3​​⎦⎤​=⎣⎡​−−−​−−−​−−−​⎦⎤​⎣⎡​a00​0b0​⎦⎤​[−−​−−​]
上面式子M做线性变换分解为：旋转(U) 拉伸( Σ \\Sigma Σ) 旋转( V T V^T VT)

图像压缩

由上面的式子可以看出来中间那个绿色方块里面对角线放的就是特征值，而且这特征值从左上角开始是从大到小排列的，那么下面几行（天蓝色this这一部分）都是0元素，即可以把他们舍掉，以减小数据量，这样就是对图像矩阵A做压缩了。

那么用数学式子表示：
A = λ 1 u 1 v 1 ′ + λ 2 u 2 v 2 ′ + ⋯ + λ k u k v k ′ + ⋯ + λ n u n v n ′ + 0 + ⋯ A=\\lambda_1u_1v_1\'+\\lambda_2u_2v_2\'+\\cdots+\\lambda_ku_kv_k\'+\\cdots+\\lambda_nu_nv_n\'+0+\\cdots A=λ1​u1​v1′​+λ2​u2​v2′​+⋯+λk​uk​vk′​+⋯+λn​un​vn′​+0+⋯
只取前k项来表示矩阵A。

• k的取值

可以根据压缩率计算出来。我们假设压缩率是rate，矩阵A总共有m*n个像素点，压缩后只有m*n*rate个像素点，上面式子一项 λ i u i v i ′ \\lambda_iu_iv_i\' λi​ui​vi′​包含像素点有 m + n + 1 m+n+1 m+n+1个，那么
k = m ∗ n ∗ r a t e m + n + 1 k=\\frac{m*n*rate}{m+n+1} k=m+n+1m∗n∗rate​
这里提一嘴，压缩率跟压缩比不一样。

图像重构

要知道将 RGB 图像拆开，每个通道单独拿出来都是一个与原图大小相同的灰度图，灰度图中像素的灰度值越大，该像素的颜色就越浅，灰度值越小颜色就越深。当要显示彩色图像时，将RGB三个分量分别放入红色、绿色和蓝色通道里，经过处理后就得到了彩色图像。

其实RGB三个通道的灰度图和颜色没有关系，如果愿意它们还可以随便定义为红橙黄绿青蓝紫，当然三原色还是最好的，其他颜色不也是可以用RGB表示出来咩。交换一下 R 和 B 通道，也就是把 R 分量送给蓝色通道，把 B 分量送给红色通道，出的图也不错喔~

matlab中有可以直接用的函数

cat(k,a1,a2,...)%k表示维数，ai表示要串联的矩阵，前提是维数相同

matlab代码‍

• 主函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%测试用例：cramp01(\'mjjy.bmp\',0.15)
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function cramp01(photo,rate)
    img=imread(photo);
    img_g=double(img);
    %颜色分离
    img_gr = img_g(:,:,1);
    img_gg = img_g(:,:,2);
    img_gb = img_g(:,:,3);
    figure(1)
    subplot(2,2,1)
    image(img_gr),title(\'颜色分离r压缩前\');
    subplot(2,2,2)
    image(img_gg),title(\'颜色分离g压缩前\');
    subplot(2,2,3)
    image(img_gb),title(\'颜色分离b压缩前\');
    %cramp_all奇异值分解
    [cramp_lst01,row_left01,col_left01,row_right01]=cramp_all(img_gr,rate);
    [cramp_lst02,row_left02,col_left02,row_right02]=cramp_all(img_gg,rate);
    [cramp_lst03,row_left03,col_left03,row_right03]=cramp_all(img_gb,rate);
    
    %cramp_lst_all =[cramp_lst01,cramp_lst02,cramp_lst03];
    %保存数据
    fileid01 = cramp_save(cramp_lst01,\'a1.svd\');
    fileid02 = cramp_save(cramp_lst02,\'a2.svd\');
    fileid03 = cramp_save(cramp_lst03,\'a3.svd\');
    %读取数据
    cramp_newlst01 = cramp_open(fileid01,row_left01,col_left01,row_right01);
    cramp_newlst02 = cramp_open(fileid02,row_left02,col_left02,row_right02);
    cramp_newlst03 = cramp_open(fileid03,row_left03,col_left03,row_right03);
    
    %图像重构
    last_lst01 = cat(3,cramp_newlst01,cramp_newlst02,cramp_newlst03);
    figure
    image(uint8(last_lst01));
    %重新rgb通道位置
    last_lst02 = cat(3,cramp_newlst02,cramp_newlst01,cramp_newlst03);
    figure
    image(uint8(last_lst02));
    
    %误差
    delta = intnum(img_g,last_lst01);
    fprintf(\'平均每个像素点误差=%.6f\',delta);   
end

• 奇异值分解

%pho:待压缩的图片
%rate:压缩率(0~1)
function [cramp_lst,row_left,col_left,row_right]=cramp_all(pho,rate) 
    %颜色分离
    %传入的pho默认已经是bouble型
    [u,s,v]=svd(pho);
    [row,col] = size(s);
    k=floor((row*col)*rate/(row+col+1));
    left=u(:,1:k);
    val=s(1:k,1:k);
    right=v(:,1:k);
    %需要把处理后的数据约束在0-255范围内
    low=min(min(pho));
    high=max(max(pho));
    value=left*val*right\';
    figure(2);
    for i=1:3
        subplot(2,2,i);
        %显示灰度图像 I，以二元素向量 [low high] 形式指定显示范围。
        imshow(value,[low,high]),title(\'压缩后图像\');
    end

    val_new = diag(val);
    val_new = val_new\';
    [row_left,col_left] = size(left);
    lala = size(right);
    row_right = lala(1);

    cramp_lst = [left;val_new;right];

end

• 数据存储

%filename：保存的文件的名字
function fileid = cramp_save(cramp_lst,filename)
    %m，n分别是cramp_lst的行列数：1551*56
    [m,n] = size(cramp_lst);
    fid = fopen(filename,\'w\');
    for i=1:m
        for j=1:n
            if j==n
                fprintf(fid,\'%g\\n\',cramp_lst(i,j));
            else
                fprintf(fid,\'%g\\t\',cramp_lst(i,j));
            end
        end
    end
    fclose(fid);
    fileid = filename;
end

• 数据读取

%数据读取
%filename是保存数据的文件
function cramp_newlst = cramp_open(filename,row_left,col_left,row_right)
    %读取数据,前面这种写入数据的方式读出来的数据是按行读取显示成一列
    fidye = fopen(filename,\'r\');
    xiaoye = fscanf(fidye,\'%g\');
    fclose(fidye);
    rowdim_all = row_left+row_right+1;
    yeye = reshape(xiaoye,[col_left,rowdim_all]);%56*1151
    piao = yeye\';
    new_left = piao(1:row_left,:);%1:633
    new_val = piao(row_left+1,:);%634
    new_right = piao(row_left+2:rowdim_all,:);%635:1551
    cramp_newlst = new_left*diag(new_val)*new_right\';
end

• 误差

function res=intnum(x,y)
    [row01,col01]=size(x);
    for i=1:row01
        for j=1:col01
            res=((abs(x(i,j)-y(i,j)))^2)/(row01*col01);
        end
    end
end

结果展示

原图：

我们用mjjy.bmp来做图像压缩，压缩率0.15：

cramp01(\'mjjy.bmp\',0.15)

分离颜色通道：

奇异值分解后：

重构：

我们这个截断奇异值分解属于有损压缩，原图有2268kb大小，原图‘mjjy.bmp’像素个数：1741383。保存的.svd格式文件像素点个数总共：260568，这么来看确实小不少哈~

jio的还阔以记得关注点起来~