2017-05-06 00:42:13 baimafujinji 阅读数 9297

数字图像处理对数学的要求颇高,这不禁令很多学习者望而却步。在阅读图像处理方面的论文时,面对梯度、散度、黑塞矩阵、傅里叶变换等这些本该在微积分中早已耳熟能详的概念时,很多人仍然感觉一筹莫展。为了弭平图像处理道路上的数学险阻,帮助更多人学好数字图像处理,并更快地具备深入研究的能力。笔者特别撰写了这本《图像处理中的数学修炼》(该书现已由清华大学出版社正式出版)。欲了解《图像处理中的数学修炼》的更多详细内容,你可以参考本书的目录

通常,我不喜欢翻开一本技术书,里面满篇满篇的都是代码。我也不希望用代码去凑页数或者挤占宝贵的篇幅。就《图像处理中的数学修炼》一书而言,更是如此。本书的重点在于全面系统地对图像处理中可能用到的各种数学基础加以总结归纳,所以数学原理才是我们所关注的重点!

但是如果把数学原理同图像处理实践割裂开来,又可能令某些读者觉得“英雄无用武之地”!为了真正帮读者建立起从数学原理到图像处理实践的桥梁,本书特意在后半部分安排了一些比较实用且非常经典的图像处理算法介绍,而书籍的前半部分则是纯数学部分,后半部分中的经典算法大量使用了前半部分里的数学知识,我们希望以这种方式来帮助读者夯实基础、巩固所学。

在后半部分的图像处理算法介绍部分,偶尔为了便于理解,本书配备了必要的MATLAB代码。但所用到的代码也非常有限的。因为这并不是一本教你怎么使用MATLAB的书!全书的代码不过六百多行,主要是为了对一些比较经典但又比较复杂的算法加以解释。这些“比较经典但又比较复杂的算法”主要包括(但不限于):

  • 基于暗通道的图像去雾实现
  • 对比度有限的自适应直方图均衡(CLAHE)
  • 自适应直方图均衡(AHE)
  • 基于小波变换的图像融合(失焦图像修复)
  • 基于泊松方程的泊松融合算法(基于迭代)
  • 基于泊松方程的泊松融合算法(有限差分方法解偏微分方程)
  • 基于主成分提取(PCA)的图像压缩

特别感谢前期试读本书的读者给予笔者的意见反馈,以及CSND博客上【图像处理中的数学原理详解】专栏的读者提出的宝贵见解。下面所列之代码已经结合之前收到的意见反馈进行了调整和优化,并修正了此前存在的一些小问题。非常感谢这部分读者为提高本书质量所做之努力。此书文字部分的勘误可以参见图像处理中的数学修炼》一书之勘误表

重要的事情说三遍:
代码或者编程不是这本书的重点!
代码或者编程不是这本书的重点!
代码或者编程不是这本书的重点!

所以全本书也就只有这么点代码。但如果你——本书的读者——对下面这些代码实现有疑问,或者发现那里有缺失,可以直接发邮件(邮件地址见博客左侧联系方式)给作者进行咨询或者获取补充代码。你也可以在“图像处理书籍读者群”中直接联系作者,提出你的疑问,或者得到你想得到的内容。面向本书读者的沟通渠道永远都是畅通的。


P270

i=double(imread('vase.tif'));
[C,S]=wavedec2(i,2,'db1');
a2=appcoef2(C,S,'db1',2);
dh1=detcoef2('h',C,S,1);
dv1=detcoef2('v',C,S,1);
dd1=detcoef2('d',C,S,1);
dh2=detcoef2('h',C,S,2);
dv2=detcoef2('v',C,S,2);
dd2=detcoef2('d',C,S,2);
[x,y]=size(i);
img = zeros(x,y);
img(1:x/4,1:y/4) =im2uint8(mat2gray(a2));
img(((x/4)+1):y/2,1:y/4) = im2uint8(mat2gray(dv2));
img(((x/4)+1):x/2,1:y/4) = im2uint8(mat2gray(dv2));
img(1:x/4,((y/4)+1):y/2) = im2uint8(mat2gray(dh2));
img(((x/4)+1):x/2,((y/4)+1):y/2) = im2uint8(mat2gray(dd2));
img(((x/2)+1):x,1:y/2) = im2uint8(mat2gray(dv1));
img(1:x/2,((y/2)+1):y) = im2uint8(mat2gray(dh1));
img(((x/2)+1):x,((y/2)+1):y) = im2uint8(mat2gray(dd1));
imshow(img,[]);

P272

X1 = imread('cathe1.bmp');
X2 = imread('cathe2.bmp');
XFUS = wfusimg(X1,X2,'sym4',5,'mean','max');
imshow(XFUS,[]);

P273

X1 = imread('cathe1.bmp');
X2 = imread('cathe2.bmp');
M1 = double(X1) / 256;
M2 = double(X2) / 256;
N = 4;
wtype = 'sym4';
[c0,s0] = wavedec2(M1, N, wtype);
[c1,s1] = wavedec2(M2, N, wtype);
length = size(c1);
Coef_Fusion = zeros(1,length(2));
%低频系数的处理,取平均值
Coef_Fusion(1:s1(1,1)) = (c0(1:s1(1,1))+c1(1:s1(1,1)))/2;
%处理高频系数,取绝对值大者,这里用到了矩阵乘法
MM1 = c0(s1(1,1)+1:length(2));
MM2 = c1(s1(1,1)+1:length(2));
mm = (abs(MM1)) > (abs(MM2));
Y  = (mm.*MM1) + ((~mm).*MM2);
Coef_Fusion(s1(1,1)+1:length(2)) = Y;
%重构
Y = waverec2(Coef_Fusion,s0,wtype);
imshow(Y,[]);

P274

I = imread('noise_lena.bmp');
[thr,sorh,keepapp] = ddencmp('den','wv',I);
de_I = wdencmp('gbl',I,'sym4',2,thr,sorh,keepapp);
imwrite(im2uint8(mat2gray(de_I)), 'denoise_lena.bmp');

P298

I = imread('baboon.bmp');  
I1 = double(I);  
T = hadamard(8);  
myFun1 = @(block_struct)T*block_struct.data*T/64;  
H = blockproc(I1, [8 8], myFun1);  
H(abs(H)<3.5)=0;  
myFun2 = @(block_struct)T*block_struct.data*T;  
I2 = blockproc(H, [8 8], myFun2);  
subplot(121), imshow(I1,[]), title('original image');  
subplot(122), imshow(I2,[]), title('zipped image');  

P299

I = imread('baboon.bmp');  
I1 = double(I);  
[m n] =size(I);  
sizi = 8;  
num = 16;  
%分块进行离散沃尔什变换  
T = hadamard(sizi);  
myFun1 = @(block_struct)T*block_struct.data*T/(sizi.^2);  
hdcoe = blockproc(I1, [sizi, sizi], myFun1);  
%重新排列系数  
coe = im2col(hdcoe,  [sizi, sizi], 'distinct');  
coe_t = abs(coe);  
[Y, ind] = sort(coe_t);  
%舍去绝对值较小的系数  
[m_c, n_c] = size(coe);  
for i = 1:n_c  
coe(ind(1:num, i), i)=0;  
end  
%重建图像  
re_hdcoe = col2im(coe, [sizi, sizi], [m, n], 'distinct');  
myFun2 = @(block_struct)T*block_struct.data*T;  
re_s = blockproc(re_hdcoe, [sizi, sizi], myFun2);  
subplot(121), imshow(I1,[]), title('original image');  
subplot(122), imshow(re_s,[]), title('compressed image');  

P307

I = imread('baboon.bmp');  
x = double(I)/255;  
[m,n]=size(x);  
y =[];  
%拆解图像  
for i = 1:m/8;  
    for j = 1:n/8;  
        ii = (i-1)*8+1;  
        jj = (j-1)*8+1;  
        y_app = reshape(x(ii:ii+7,jj:jj+7),1,64);  
        y=[y;y_app];  
    end  
end  
  
%KL变换  
[COEFF,SCORE,latent] = princomp(y);  
kl = y * COEFF;  
  
kl1 = kl;  
kl2 = kl;  
kl3 = kl;  
  
%置零压缩过程  
kl1(:, 33:64)=0;  
kl2(:, 17:64)=0;  
kl3(:, 9:64)=0;  
  
%KL逆变换  
kl_i = kl*COEFF';  
kl1_i = kl1*COEFF';  
kl2_i = kl2*COEFF';  
kl3_i = kl3*COEFF';  
  
image = ones(256,256);  
image1 = ones(256,256);  
image2 = ones(256,256);  
image3 = ones(256,256);  
  
k=1;  
%重组图像  
for i = 1:m/8;  
    for j = 1:n/8;  
  
        y = reshape(kl_i(k, 1:64),8,8);  
        y1 = reshape(kl1_i(k, 1:64),8,8);  
        y2 = reshape(kl2_i(k, 1:64),8,8);  
        y3 = reshape(kl3_i(k, 1:64),8,8);  
  
        ii = (i-1)*8+1;  
        jj = (j-1)*8+1;  
  
        image(ii:ii+7,jj:jj+7) = y;  
        image1(ii:ii+7,jj:jj+7) = y1;  
        image2(ii:ii+7,jj:jj+7) = y2;  
        image3(ii:ii+7,jj:jj+7) = y3;  
  
        k=k+1;  
    end  
end  

##P356-1

mountains = double(imread('./img/mountain.jpg'));
moon = double(imread('./img/moon.png'));
sizeSrc = size(mountains);
sizeDst = size(moon);

gradient_inner = moon(1:sizeDst(1)-2,2:sizeDst(2)-1,:)...
    + moon(3:sizeDst(1),2:sizeDst(2)-1,:)...
    + moon(2:sizeDst(1)-1,1:sizeDst(2)-2,:)...
    + moon(2:sizeDst(1)-1,3:sizeDst(2),:)...
    - 4*moon(2:sizeDst(1)-1,2:sizeDst(2)-1,:);

P356-2

Lap = [0, 1, 0;1, -4, 1;0, 1, 0];

I1 = conv2(double(moon(:,:,1)), double(Lap));
I2 = conv2(double(moon(:,:,2)), double(Lap));
I3 = conv2(double(moon(:,:,3)), double(Lap));
gradient_inner(:, :, 1) = I1(3:sizeDst(1),3:sizeDst(2));
gradient_inner(:, :, 2) = I2(3:sizeDst(1),3:sizeDst(2));
gradient_inner(:, :, 3) = I3(3:sizeDst(1),3:sizeDst(2));

P357

dstX = 350;dstY = 100;
rebuilt = mountains(dstY:dstY+sizeDst(1)-1,dstX:dstX+sizeDst(2)-1,:);

for n = [1:1000]
    rebuilt(2:2:sizeDst(1)-1,2:2:sizeDst(2)-1,:)= ...
        (rebuilt(1:2:sizeDst(1)-2 , 2:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1) , 2:2:sizeDst(2)-1,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 1:2:sizeDst(2)-2,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 3:2:sizeDst(2),:)...
        -gradient_inner(1:2:sizeDst(1)-2 , 1:2:sizeDst(2)-2,:))/4;
     rebuilt(3:2:sizeDst(1)-1,3:2:sizeDst(2)-1,:)= ...
        (rebuilt(2:2:sizeDst(1)-2 , 3:2:sizeDst(2)-1,:)...
        +rebuilt(4:2:sizeDst(1) , 3:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 2:2:sizeDst(2)-2,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 4:2:sizeDst(2),:)...
        -gradient_inner(2:2:sizeDst(1)-2 , 2:2:sizeDst(2)-2,:))/4;
     rebuilt(3:2:sizeDst(1)-1,2:2:sizeDst(2)-1,:)= ...
        (rebuilt(2:2:sizeDst(1)-2 , 2:2:sizeDst(2)-1,:)...
        +rebuilt(4:2:sizeDst(1) , 2:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 1:2:sizeDst(2)-2,:)...
        +rebuilt(3:2:sizeDst(1)-1 , 3:2:sizeDst(2),:)...
        -gradient_inner(2:2:sizeDst(1)-2 , 1:2:sizeDst(2)-2,:))/4;
    rebuilt(2:2:sizeDst(1)-1 , 3:2:sizeDst(2)-1,:)= ...
        (rebuilt(1:2:sizeDst(1)-2 , 3:2:sizeDst(2)-1,:)...
        +rebuilt(3:2:sizeDst(1) , 3:2:sizeDst(2)-1,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 2:2:sizeDst(2)-2,:)...
        +rebuilt(2:2:sizeDst(1)-1 , 4:2:sizeDst(2),:)...
        -gradient_inner(1:2:sizeDst(1)-2 , 2:2:sizeDst(2)-2,:))/4;
end

mountains(dstY:sizeDst(1)+dstY-1,dstX:sizeDst(2)+dstX-1,:) = rebuilt;
figure,imshow(uint8(mountains));

P360-P361

TargetImg   = imread('pool-target.jpg');  
SourceImg   = imread('bear.jpg');  
SourceMask  = im2bw(imread('bear-mask.jpg'));  

[SrcBoundry, ~] = bwboundaries(SourceMask, 8);
figure, imshow(SourceImg), axis image  
hold on  
for k = 1:length(SrcBoundry)  
    boundary = SrcBoundry{k};  
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2)  
end  
title('Source image intended area for cutting from');

position_in_target = [10, 225];%xy
[TargetRows, TargetCols, ~] = size(TargetImg);
[row, col] = find(SourceMask);
start_pos = [min(col), min(row)];
end_pos  = [max(col), max(row)];
frame_size = end_pos - start_pos;

if (frame_size(1) + position_in_target(1) > TargetCols)
    position_in_target(1) = TargetCols - frame_size(1);
end

if (frame_size(2) + position_in_target(2) > TargetRows)
    position_in_target(2) = TargetRows - frame_size(2);
end

MaskTarget = zeros(TargetRows, TargetCols);
MaskTarget(sub2ind([TargetRows, TargetCols], row-start_pos(2)+ ...
position_in_target(2), col-start_pos(1)+position_in_target(1))) = 1;

TargBoundry = bwboundaries(MaskTarget, 8);
figure, imshow(TargetImg), axis image
hold on
for k = 1:length(TargBoundry)
    boundary = TargBoundry{k};
    plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1)
end

P362

templt = [0 -1 0; -1 4 -1; 0 -1 0];  
LaplacianSource = imfilter(double(SourceImg), templt, 'replicate');  
VR = LaplacianSource(:, :, 1);  
VG = LaplacianSource(:, :, 2);  
VB = LaplacianSource(:, :, 3);

TargetImgR = double(TargetImg(:, :, 1));  
TargetImgG = double(TargetImg(:, :, 2));  
TargetImgB = double(TargetImg(:, :, 3));  
  
TargetImgR(logical(MaskTarget(:))) = VR(SourceMask(:));  
TargetImgG(logical(MaskTarget(:))) = VG(SourceMask(:));  
TargetImgB(logical(MaskTarget(:))) = VB(SourceMask(:));  
  
TargetImgNew = cat(3, TargetImgR, TargetImgG, TargetImgB);  
figure, imagesc(uint8(TargetImgNew)), axis image;

AdjacencyMat = calcAdjancency(MaskTarget);

ResultImgR=PoissonSolver(TargetImgR,MaskTarget,AdjacencyMat,TargBoundry);
ResultImgG=PoissonSolver(TargetImgG,MaskTarget,AdjacencyMat,TargBoundry);
ResultImgB=PoissonSolver(TargetImgB,MaskTarget,AdjacencyMat,TargBoundry);

ResultImg = cat(3, ResultImgR, ResultImgG, ResultImgB);

%% Show the final results
figure;
imshow(uint8(ResultImg));

特别说明:函数calcAdjancency和PoissonSolver的实现请在图像处理学习群中直接联系群主获取,群号见文末。


P393

image = imread('Unequalized_Hawkes_Bay_NZ.jpg');  
Img = rgb2gray(image);  
[height,width]=size(Img);

NumPixel = zeros(1,256);%统计各灰度数目,共256个灰度级
for i = 1:height  
    for j = 1: width  
    %对应灰度值像素点数量增加一
    %因为NumPixel的下标是从1开始,但是图像像素的取值范围是0~255,
    %所以用NumPixel(Img(i,j) + 1)  
    NumPixel(Img(i,j) + 1) = NumPixel(Img(i,j) + 1) + 1;  
    end  
end  

ProbPixel = zeros(1,256);
for i = 1:256  
    ProbPixel(i) = NumPixel(i) / (height * width * 1.0);
end 

CumuPixel = cumsum(ProbPixel);  
CumuPixel = uint8(255 .* CumuPixel + 0.5);  

for i = 1:height  
    for j = 1: width  
        Img(i,j) = CumuPixel(Img(i,j));  
    end  
end

imshow(Img)

P395-1

a = imread('couple.tiff');
R = a(:,:,1);
G = a(:,:,2);
B = a(:,:,3);
      
R = histeq(R, 256);
G = histeq(G, 256);
B = histeq(B, 256);
      
a(:,:,1) = R;
a(:,:,2) = G;
a(:,:,3) = B;
imshow(a)

P395-2

Img = imread('couple.tiff');
hsvImg = rgb2hsv(Img);
V = hsvImg(:,:,3);
[height,width] = size(V);

V = uint8(V*255);
NumPixel = zeros(1,256);
for i = 1:height
    for j = 1: width
        NumPixel(V(i,j) + 1) = NumPixel(V(i,j) + 1) + 1;
    end  
end

ProbPixel = zeros(1,256);
for i = 1:256
    ProbPixel(i) = NumPixel(i) / (height * width * 1.0);
end

CumuPixel = cumsum(ProbPixel);
CumuPixel = uint8(255 .* CumuPixel + 0.5);
  
for i = 1:height  
    for j = 1: width  
        if V(i,j)==0
			V(i,j) = CumuPixel(V(i,j)+1);
		else 
			V(i,j) = CumuPixel(V(i,j));
		end 
    end  
end  
   
V = im2double(V);
hsvImg(:,:,3) = V;
outputImg = hsv2rgb(hsvImg);
imshow(outputImg);

P397

img = imread('space.jpg');
rimg = img(:,:,1);
gimg = img(:,:,2);
bimg = img(:,:,3);
resultr = adapthisteq(rimg);
resultg = adapthisteq(gimg);
resultb = adapthisteq(bimg);
result = cat(3, resultr, resultg, resultb);
imshow(result);

P398-1

clear;  
img = imread('space.jpg');  
cform2lab = makecform('srgb2lab');  
LAB = applycform(img, cform2lab);  
L = LAB(:,:,1);  
LAB(:,:,1) = adapthisteq(L);  
cform2srgb = makecform('lab2srgb');  
J = applycform(LAB, cform2srgb);  
imshow(J);  

P398-2

clc;
clear all;
Img = rgb2gray(imread('space.jpg'));
[h,w] = size(Img);  
minV = double(min(min(Img)));
maxV = double(max(max(Img)));
imshow(Img);

P399

NrX = 8;
NrY = 4;
HSize = ceil(h/NrY);
WSize = ceil(w/NrX);
      
deltay = NrY*HSize - h;  
deltax = NrX*WSize - w;  
      
tmpImg = zeros(h+deltay,w+deltax);  
tmpImg(1:h,1:w) = Img;  

new_w = w + deltax;
new_h = h + deltay;
NrPixels = WSize * WSize;

% NrBins - Number of greybins for histogram ("dynamic range")
NrBins = 256;

LUT = zeros(maxV+1,1);

for i=minV:maxV  
    LUT(i+1) = fix(i - minV);%i+1
end  

Bin = zeros(new_h, new_w);  
for m = 1 : new_h  
    for n = 1 : new_w  
        Bin(m,n) = 1 + LUT(tmpImg(m,n) + 1);
    end  
end  

Hist = zeros(NrY, NrX, 256);
for i=1:NrY  
    for j=1:NrX  
        tmp = uint8(Bin(1+(i-1)*HSize:i*HSize, 1+(j-1)*WSize:j*WSize));
        %tmp = tmpImg(1+(i-1)*HSize:i*HSize,1+(j-1)*WSize:j*WSize);
        [Hist(i, j, :), x] = imhist(tmp, 256);  
    end  
end

Hist = circshift(Hist,[0, 0, -1]);

ClipLimit = 2.5;  
ClipLimit = max(1,ClipLimit * HSize * WSize/NrBins);
Hist = clipHistogram(Hist,NrBins,ClipLimit,NrY,NrX);
Map=mapHistogram(Hist, minV, maxV, NrBins, NrPixels, NrY, NrX);

yI = 1;  
for i = 1:NrY+1  
    if i == 1  
        subY = floor(HSize/2);
        yU = 1;
        yB = 1;
    elseif i == NrY+1  
        subY = floor(HSize/2);
        yU = NrY;
        yB = NrY;
    else
        subY = HSize;
        yU = i - 1;
        yB = i;
    end
    xI = 1;
    for j = 1:NrX+1
        if j == 1  
            subX = floor(WSize/2);
            xL = 1;
            xR = 1;
        elseif j == NrX+1
            subX = floor(WSize/2);
            xL = NrX;
            xR = NrX;
        else
            subX = WSize;
            xL = j - 1;
            xR = j;
        end
        UL = Map(yU,xL,:);
        UR = Map(yU,xR,:);
        BL = Map(yB,xL,:);
        BR = Map(yB,xR,:);
        subImage = Bin(yI:yI+subY-1,xI:xI+subX-1);  

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        sImage = zeros(size(subImage));
        num = subY * subX;
        for i = 0:subY - 1
            inverseI = subY - i;
            for j = 0:subX - 1
                inverseJ = subX - j;
                val = subImage(i+1,j+1);
                sImage(i+1, j+1)=(inverseI*(inverseJ*UL(val)+j*UR(val))...
                                   + i*(inverseJ*BL(val)+j*BR(val)))/num;
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        output(yI:yI+subY-1, xI:xI+subX-1) = sImage;
        xI = xI + subX;
    end
    yI = yI + subY;
end  

output = output(1:h, 1:w);
figure, imshow(output, []);

P404

img = rgb2gray(imread('theatre.jpg'));
img_ref = rgb2gray(imread('rpic.jpg'));
[hgram, x] = imhist(img_ref);
J = histeq(img, hgram);
subplot(2,3,1), imshow(img), title('original image');
subplot(2,3,4), imhist(img), title('original image');
subplot(2,3,2), imshow(img_ref), title('reference image');
subplot(2,3,5), imhist(img_ref), title('reference image');
subplot(2,3,3), imshow(J), title('output image');
subplot(2,3,6), imhist(J), title('output image');

P409

%求一幅图像的暗通道图,窗口大小为15*15
imageRGB = imread('picture.bmp');
imageRGB = double(imageRGB);
imageRGB = imageRGB./255;
dark = darkChannel(imageRGB);

% 选取暗通道图中最亮的0.1%像素,从而求得大气光
[m, n, ~] = size(imageRGB);
imsize = m * n;
numpx = floor(imsize/1000);
JDarkVec = reshape(dark,imsize,1);
ImVec = reshape(imageRGB,imsize,3);

[JDarkVec, indices] = sort(JDarkVec);
indices = indices(imsize-numpx+1:end);

atmSum = zeros(1,3);
for ind = 1:numpx
    atmSum = atmSum + ImVec(indices(ind),:);
end

atmospheric = atmSum / numpx;

%求解透射率,并通过omega参数来选择保留一定程度的雾霾,以免损坏真实感
omega = 0.95; 
im = zeros(size(imageRGB));

for ind = 1:3 
    im(:,:,ind) = imageRGB(:,:,ind)./atmospheric(ind);
end

dark_2 = darkChannel(im);
t = 1-omega*dark_2;

%通过导向滤波来获得更为精细的透射图
r = 60;
eps = 10^-6;
refined_t = guidedfilter_color(imageRGB, t, r, eps);
refinedRadiance = getRadiance(atmospheric, imageRGB, refined_t);

P410

function dark = darkChannel(imRGB)

r=imRGB(:,:,1);
g=imRGB(:,:,2);
b=imRGB(:,:,3);

[m n] = size(r);
a = zeros(m,n);
for i = 1: m     
    for j = 1: n
         a(i,j) = min(r(i,j), g(i,j));
         a(i,j)= min(a(i,j), b(i,j));
    end
end

d = ones(15,15);
fun = @(block_struct)min(min(block_struct.data))*d;
dark = blockproc(a, [15 15], fun); 

dark = dark(1:m, 1:n);

本书源起

数字图像处理对数学的要求颇高,这不禁令很多学习者望而却步。在阅读图像处理方面的论文时,面对梯度、散度、黑塞矩阵、傅里叶变换等这些本该在微积分中早已耳熟能详的概念时,很多人仍然感觉一筹莫展。

大约三年前,我开始在博客上介绍一些数字图像处理中的数学基础和数学原理,并专门开设了一个“图像处理中的数学原理”专栏,其中大概收录了三十几篇文章。很多人开始留言询问我,有没有这样一本书?彼时我并没有将该系列文章付梓的想法。但是耐不住众多热情的读者不断敦促我将这个系列相对完整和系统地整理成一本可读性更强的书籍。于是便有了这本书——《图像处理中的数学修炼》。

大约一年前清华出版社的编辑同我联系,我想或许时机已到,遂开始着手将千头万绪已经发布的和未曾发布的资料整理成书。期间为了能够更好地呈现这本书,出版日期也是一拖再拖。从最初的十一,到双十一,再到春节前,直到现在的烟花三月。这本《图像处理中的数学修炼》可真算是“千呼万唤始出来”。期间很多热情的读者三番五次地问询该书的出版进度,着实令我感觉身上压力不小。无奈好事毕竟多磨,所幸终于瓜熟蒂落。现在这本书终于同大家见面了!


2015-12-14 14:18:12 baimafujinji 阅读数 7193

欢迎关注我的博客专栏“图像处理中的数学原理详解

全文目录请见 图像处理中的数学原理详解(总纲)

http://blog.csdn.net/baimafujinji/article/details/48467225



1.1.3 函数的极限


本小节介绍两个重要的函数极限,并讨论它们的应用。

重要极限1:


此外,该重要极限的另一种形式也常常被用到,即


综上,结论得证。
由此,也很容易推出如下结论,证明从略,有兴趣的读者可以自行尝试推导



我的“图像处理中的数学原理”专栏中之系列文章已经以《图像处理中的数学修炼》为名结集出版(清华大学出版社)。该书详细介绍图像处理中的数学原理,为你打开一道通往图像世界的数学之门,详细内容及目录请见 http://blog.csdn.net/baimafujinji/article/details/48467225



2015-09-15 15:10:15 baimafujinji 阅读数 37711

数字图像处理技术的研究与开发对数学基础的要求很高,一些不断涌现的新方法中,眼花缭乱的数学推导令很多期待深入研究的人望而却步。一个正规理工科学生大致已经具备了包括微积分、线性代数、概率论在内的数学基础。但在分析一些图像处理算法的原理时,好像还是感觉无从入手。实际所牵涉出来的问题主要可归结为如下几个原因:1)微积分、线性代数、概率论这些是非常重要的数学基础,但显示不是这些课程中所有的内容都在图像处理算法中有直接应用;2)当你将图像处理和数学分开来学的时候,其实并没有设法建立它们二者的联系;3)一些新方法或者所谓的高大上的算法之基础已经超过了上面三个数学课程所探讨的基本领域,这又涉及到偏微分方程、变分法、复变函数、实变函数、泛函分析等等;4)如果你不是数学科班出身,要想自学上面所谈到所有内容,工作量实在太过浩繁,恐怕精力也难以顾及。

业余时间,笔者结合自己对图像处理的学习和实践,大致总结了一部分图像处理研究中所需的数学原理基础。这些内容主要涉及微积分、向量分析、场论、泛函分析、偏微分方程、复变函数、变分法等。线性代数和概率论笔者认为比较基础,于是并没有将其收入。我总结、归纳、提取了上面这些数学课程中,在研究图像处理最容易碰到也最需要知道的一些知识点,然后采取一种循序渐进的方式将它们重新组织到了一起。并结合具体的图像处理算法讨论来讲解这些数学知识的运用。从而建立数学知识与图像处理之间的一座桥梁。

这部分内容主要是笔者日常研究和学习的一个总结,我原本并没有将其出版的计划(毕竟这个Topic还是有点小众而且可能还有点艰深)。但之前我撷取了其中的一小部分发到了网上,已经有读者表现出了浓厚的兴趣。再后来亦有多位出版社的编辑联系到我,希望可以将该书付梓。而且不知不觉中,这个系列专栏的内容日积月累,个人感觉确实已经形成了一个相对比较完整的体系,于是便应下了出版社的合作意向。


我的“图像处理中的数学原理”专栏中之系列文章已经以《图像处理中的数学修炼》为名结集出版(清华大学出版社)。该书详细介绍图像处理中的数学原理,为你打开一道通往图像世界的数学之门。以下是最新版本的该书的完整目录,方便各位网友查阅以及确定本书是否符合你的选购目标:


第1章 必不可少的数学基础

1.1  极限及其应用
    1.1.1  数列的极限
    1.1.2  级数的敛散
    1.1.3  函数的极限
    1.1.4  极限的应用

1.2  微分中值定理
    1.2.1  罗尔中值定理
    1.2.2  拉格朗日中值定理
    1.2.3  柯西中值定理
    1.2.4  泰勒公式

    1.2.5  海塞矩阵与多元函数极值

1.3  向量代数与场论
    1.3.1  牛顿-莱布尼茨公式
    1.3.2  内积与外积
    1.3.3  方向导数与梯度
    1.3.4  曲线积分
    1.3.5  格林公式
    1.3.6  积分与路径无关条件
    1.3.7  曲面积分
    1.3.8  高斯公式与散度
    1.3.9  斯托克斯公式与旋度
本章参考文献

第2章 更进一步的数学内容

2.1  傅立叶级数展开
    2.1.1  函数项级数的概念
    2.1.2  函数项级数的性质
    2.1.3  傅立叶级数的概念
    2.1.4  傅立叶变换的由来
    2.1.5  卷积定理及其证明

2.2  复变函数论初步
    2.2.1  解析函数
    2.2.2  复变积分
    2.2.3  基本定理
    2.2.4  级数展开

2.3  凸函数与詹森不等式

    2.3.1  凸函数的概念

    2.3.2  詹森不等式及其证明

    2.3.3  詹森不等式的应用

2.4  常用经典数值解法

    2.4.1  牛顿迭代法

    2.4.2  雅各比迭代

    2.4.3  高斯迭代法

    2.4.4  托马斯算法

本章参考文献


第3章  泛函分析以及变分法

3.1  勒贝格积分理论
    3.1.1  点集的勒贝格测度
    3.1.2  可测函数及其性质
    3.1.3  勒贝格积分的定义
    3.1.4  积分序列极限定理
3.2  泛函与抽象空间
    3.2.1  线性空间
    3.2.2  距离空间
    3.2.3  赋范空间
    3.2.4  巴拿赫空间
    3.2.5  内积空间
    3.2.6  希尔伯特空间
    3.2.7  索伯列夫空间
3.3  从泛函到变分法
    3.3.1  理解泛函的概念
    3.3.2  关于的变分概念
    3.3.3  变分法的基本方程
    3.3.4  理解哈密尔顿原理
    3.3.5  等式约束下的变分
    3.3.6  巴拿赫不动点定理
    3.3.7  有界变差函数空间

本章参考文献


第4章 概率论与统计学基础

4.1  概率论的基本概念

4.2  随机变量数字特征

    4.2.1  期望

    4.2.2  方差

    4.2.3  矩与矩母函数

    4.2.4  协方差与协方差矩阵

4.3  基本概率分布模型

    4.3.1  离散概率分布

    4.3.2  连续概率分布

4.4  概率论中的重要定理

    4.4.1  大数定理

    4.4.2  中央极限定理

4.5  随机采样

    4.5.1  随机采样分布

    4.5.2  蒙特卡洛采样

4.6  参数估计

4.7  假设检验

    4.7.1  基本概念

    4.7.2  两类错误

    4.7.3  均值检验

4.8  极大似然估计

    4.8.1  极大似然法的基本原理

    4.8.2  求极大似然估计的方法

4.9  贝叶斯推断

    4.9.1  先验概率与后验概率

    4.9.2  共轭分布

参考文献

第5章 子带编码与小波变换

5.1  图像编码的理论基础
    5.1.1  率失真函数
    5.1.2  香农下边界
    5.1.3  无记忆高斯信源
    5.1.4  有记忆高斯信源
5.2  子带编码基本原理
    5.2.1  数字信号处理基础
    5.2.2  多抽样率信号处理
    5.2.3  图像信息子带分解
5.3  哈尔函数及其变换
    5.3.1  哈尔函数的定义
    5.3.2  哈尔函数的性质
    5.3.3 酉矩阵与酉变换
    5.3.4 二维离散线性变换
    5.3.5 哈尔基函数
    5.3.6 哈尔变换
5.4  小波及其数学原理
    5.4.1  小波的历史
    5.4.2  理解小波的概念
    5.4.3  多分辨率分析
    5.4.4  小波函数的构建
    5.4.5  小波序列展开
    5.4.6  离散小波变换
    5.4.7  连续小波变换
    5.4.8  小波的容许条件与基本特征
5.5  快速小波变换算法
    5.5.1  快速小波正变换
    5.5.2  快速小波逆变换
    5.5.3  图像的小波变换

5.6  小波在图像处理中的应用

本章参考文献


第6章 正交变换与图像压缩
6.1  傅立叶变换
    6.1.1  信号处理中的傅立叶变换

        1. 连续时间,连续频率——傅立叶变换

        2. 连续时间,离散频率——傅立叶级数

        3. 离散时间,连续频率——序列的傅立叶变换

        4. 离散时间,离散频率——离散的傅立叶变换

    6.1.2  数字图像的傅立叶变换
    6.1.3  快速傅立叶变换的算法
6.2  离散余弦变换
    6.2.1  基本概念及数学描述
    6.2.2  离散余弦变换的快速算法
    6.2.3  离散余弦变换的意义与应用
6.3  沃尔什-阿达马变换
    6.3.1  沃尔什函数
    6.3.2  离散沃尔什变换及其快速算法
    6.3.3  沃尔什变换的应用
6.4  卡洛南-洛伊变换
    6.4.1  一些必备的基础概念
    6.4.2  主成分变换的推导
    6.4.3  主成分变换的实现
    6.4.4  基于K-L变换的图像压缩
本章参考文献

第7章 无所不在的高斯分布
7.1  卷积积分与邻域处理
    7.1.1  卷积积分的概念
    7.1.2  模板与邻域处理
    7.1.3  图像的高斯平滑
7.2  边缘检测与微分算子
    7.2.1  哈密尔顿算子
    7.2.2  拉普拉斯算子
    7.2.3  高斯-拉普拉斯算子
    7.2.4  高斯差分算子
7.3  保持边缘的平滑处理
    7.3.1  双边滤波算法应用
    7.3.2  各向异性扩散滤波
    7.3.3  基于全变差的方法
7.4  数学物理方程的应用
    7.4.1  泊松方程的推导
    7.4.2  图像的泊松编辑
    7.4.3  离散化数值求解

    7.4.4  基于稀疏矩阵的解法

7.5  多尺度空间及其构建
    7.5.1  高斯滤波与多尺度空间的构建
    7.5.2  基于各向异性扩散的尺度空间
本章参考文献


 

第8章 处理彩色图像

8.1  从认识色彩开始

    8.1.1  什么是颜色

    8.1.2  颜色的属性

        1. 色相

        2. 亮度

        3. 纯度

    8.1.3  光源能量分布图

8.2  CIE色度图

    8.2.1  CIE色彩模型的建立

    8.2.2  CIE色度图的理解

        1. 确定互补颜色

        2. 确定色光主波

        3. 定义颜色区域

    8.2.3  CIE色度图的后续发展

8.3  常用的色彩空间

    8.3.1  RGB颜色空间

    8.3.2 CMY/CMYK颜色空间

    8.3.3  HSV/HSB颜色空间

    8.3.4  HSI/HSL颜色空间

    8.3.5  Lab颜色空间

    8.3.6 YUV/YCbCr颜色空间

8.4  色彩空间的转换方法

    8.4.1  RGB转换到HSV的方法

    8.4.2  RGB转换到HSI的方法

    8.4.3  RGB转换到YUV的方法

    8.4.4  RGB转换到YCbCr的方法

8.5  基于直方图的色彩增强

    8.5.1     普通直方图均衡

    8.5.2    CLAHE算法

    8.5.3     直方图规定化

8.6  暗通道先验的去雾算法

    8.6.1  暗通道的概念与意义

    8.6.2  暗通道去雾霾的原理

    8.6.3  算法实现与应用

本章参考文献




2018-04-27 14:41:53 jayandchuxu 阅读数 342

矩阵的特征值、特征向量的概念

这里,我们讨论的是n阶的方阵A

定义

从向量的定义可知,它是方向和长度的结合体。当一个线性变换A作用在n维线性空间V中的某一非零向量x上时,便是对该向量的长度和方向进行变化。然而,存在一些向量,线性变换A并没有改变其方向,而只是改变了长度,这种向量,叫做线性变换A特征向量,它在变换中被改变的倍数,叫做它的特征值。用数学公式表示这一概念,即:

Ax=λx(1)

其中,λ的个线性变换A的某一个特征值。从公式上可以轻易发现,如果某一向量x是线性变换A的特征向量,那么与其方向相同的任意长度(不为零)向量,都是A的特征向量,并且属于同一个特征值λ。由于相同方向的特征向量具有相同特征值,我们可以同特征值来描述这一族向量(同一个特征值,可以有多个方向的特征向量)。
从公式(1)中,可以看出特征向量和特征值的计算方法:
|λEA|=0(2)
(λEA)x=0(3)

对应于同一个线性变换A,可以有多个特征向量(方向不同),但是有多个特征向量可以对应同一个特征值。 一个向量是一个方向,两个不同方向的向量就可以张成一个空间。在相同特征值的特征向量张成的空间内,任何一个向量在变换A下,都有相同的放大倍数λ

公式(2)的左侧,总可以展成如下形式的多项式:

|λEA|=λn(a11+a22++ann)λn1++(1)n|A|

所以求特征值就是求下面方程的解:

λn(a11+a22++ann)λn1++(1)n|A|=0(4)

关于从方程(4)得到的特征值,有几个比较重要的结论(参考资料1):

  1. n阶矩阵在复数范围内,一定有n个特征值(重特征值按重数计算个数)。
  2. n阶矩阵在实数范围内有多少个特征值是不一定的
  3. n阶实对称矩阵可以看成是一个特例,因为它一定有n个实特征值(重特征值按重数计算个数)。如果其中一个特征值λ=0,矩阵的秩r(A)=k,(0<k<nk是正整数),则λ=0恰为Ank重特征值。
  4. 如果 n阶矩阵A不是对称矩阵,那么,λ=0至少为Ank重特征值。

关于特征值和特征向量的理解,参考资料3写的也很好,知乎上有很多大神的回答直击要害,对问题的理解很有帮助。

作用

参考资料4中,认为矩阵的变换有三个作用:旋转拉伸投影
A是一个n×n方阵时,只涉及到旋转,拉伸。如果矩阵在实数域内可以得到n个特征值(重特征值按重数计算个数),那么利用其对应的特征向量(单位化后)组成矩阵正交Q可以使得:

A=QΛQ1

其中正交矩阵Q起到旋转作用(旋转矩阵都是正交矩阵,且行列式都为1),对角矩阵Λ起拉伸作用。
当矩阵不是方阵而是m×n时,可以对其进行SVD分解
在之前,想研究一下正交矩阵。

正交矩阵

按照定义,正交矩阵是QQT=E,它的行列式为1或者-1。

正交矩阵的性质

了解正交矩阵的性质,在很多计算方面,能够更深入了解所进行的运算的意义。
我们这里说的都是有限维欧式空间内的正交矩阵

  1. 正交矩阵的转置伴随矩阵、之间的积矩阵都是正交矩阵;
  2. .每一行(列)都是单位向量
  3. 任意两行或两列相互垂直
  4. 其行列式等于±1

假设n维欧式空间Rn中,一个正交变换Q存在一个一一对应的正交矩阵Q。所以研究正交变换的性质,可以转为研究正交矩阵。
根据正交矩阵行列式的值,将其分为两类:|Q|=1,为第一类;|Q|=1,为第二类。
第一类正交矩阵,当其左乘一个向量时,几何意义是使该量在Oxyz坐标系下旋转;
第一类正交矩阵,当其左乘一个向量时,几何意义是使该向量沿Oxyz某一轴(点)进行反射;[5]
无论是哪一类正交矩阵,其左乘向量,均不会改变向量的长度,即|Qv|=|Q||v|=|v|
所以上面将矩阵A拆成A=QΛQ1的形式后,由于Q是正交矩阵,它作用在向量v上,只是对其进行旋转或者反射,而对向量长度有影响的是矩阵Λ。其上的元素由矩阵A的特征值组成。
需要注意的是,只有当矩阵A能够在实数域内有n个特征值(重数也计入)的情况下,可以如此分解。但是,对更多的一般性矩阵A,在实数域内没有n个实特征根,或者是更一般的m×n维矩阵,此时,我们用SVD分解的方法进行研究。

SVD分解

协方差

假设有两个变量XY,他们的取值为X={x1,x2,,xn}Y={y1,y2,,yn},他们的平均值(期望)记为是E(X)=X¯=1nxiE(Y)=Y¯=1nyiXY的协方差就是研究两个变量之间的相关关系。比如当一个的取值不断增大时,另一个变量的取值如何变化,知乎上的一篇文章( 如何通俗易懂地解释「协方差」与「相关系数」的概念?)讲的很好。

根据奇异值分解

n×n阶矩阵按特征值分解相似,任一m×n阶矩阵A也可以写成类似的形式:

A=UΣVT(5)

那么得到的U是一个m×m的方阵(里面的向量是正交的,U里面的向量称为左奇异向量),Σ是一个m×n的矩阵(除了对角线的元素都是0,对角线上的元素称为奇异值),VT(V的转置)是一个n×n的矩阵,里面的向量也是正交的,V里面的向量称为右奇异向量)。

求法

利用如下公式可以计算出各矩阵:

(ATA)vi=λivi(6)

δi=λi(7)

ui=1δiAvi(8)

对方程计算,得到的v,就是的右奇异向量σ就是奇异值u就是左奇异向量

作用

当矩阵A=UΣVT左乘一个向量v时,只有Σ对向量的长度进行了拉伸(收缩),而矩阵UV都只是对其进行旋转或反射。当作用在图像上时,也只有Σ对图像进行了各个方向上的伸缩改变。

用奇异值分解图像

这里写图片描述
用奇异值的方法,将这幅图像进行分解,得形如A=UΣVT格式的矩阵。其中Σ是由矩阵奇异值由大到小排列组成的对角矩阵。
我们分别保留前10,30,100,300个奇异值,其余奇异值设为0,比较图像的变化:
奇异值保留前10
奇异值保留前10
奇异值保留前30
奇异值保留前30
奇异值保留前100
奇异值保留前100
奇异值保留前300
奇异值保留前300。

代码

import cv2
import numpy as np
from numpy import linalg as la  # 用到别名
from scipy.misc import imsave

import scipy
im = cv2.imread('lena512.bmp')
print(im.shape)
gray = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
U, Sigma, VT = la.svd(gray)
print('矩阵U的形状:', U.shape, '    矩阵Sigma的形状:',
      Sigma.shape, '    矩阵VT的形状:', VT.shape)
se = np.eye(512, dtype=np.float64)
n = 512
i = 0
k = 30  # 保留特征值数目
# 改变特征值
while i < n:
    if i > k - 1:
        Sigma[i] = 0
    se[i, i] = Sigma[i]
    i += 1

svt = np.dot(se, VT)
usvt = np.dot(U, svt)
imsave('USVT_Sigma=30m.bmp', usvt)

参考资料

  1. 秦川, 李小飞. 方阵的秩与特征值的关系[J]. 课程教育研究:学法教法研究, 2015(27):120-120.
  2. 如何通俗易懂地解释「协方差」与「相关系数」的概念?
  3. 如何理解矩阵特征值?马同学的回答
  4. 矩阵的特征值分解与奇异值分解的几何意义
  5. 杜美华, 孙建英. 正交变换的几何意义及其应用[J]. 哈尔滨师范大学自然科学学报, 2014, 30(3):36-39.
2016-01-16 16:38:16 baimafujinji 阅读数 12033

 

欢迎关注我的博客专栏“图像处理中的数学原理详解

全文目录请见 图像处理中的数学原理详解(总纲)

http://blog.csdn.net/baimafujinji/article/details/48467225

 

 

有段时间没继续更新我的“图像处理中的数学原理详解”专栏了。因为前面基础的部分已经发布的差不多了,现在已经进入 “深水区”。一方面现在文章的长度都有所增加,所以我写起来就更加麻烦了。另一方面,现在的话题进入了微分方程和泛函分析领域,这部分内容对于非数学专业的人来说实在太难了:(  。而且如果是初学图像处理的人,也就是还在琢磨高斯平滑和中值滤波的人,基本上也用不到这些东西。所以如果你属于这种情况,我不建议你看本文。还有,如果你的目的是写一个类似Magic House的程序(http://blog.csdn.net/baimafujinji/article/details/50500757),那么你也不用学这么深。而且研究下面这些话题的人,更适合用MATLAB。

我的意见是那些每天要看Paper的图像处理研究者来学微分方程和泛函分析。如果你每天都要看一些满篇公式的图像处理paper,却不懂L1范数和L2范数、不懂希尔伯特空间、不懂泊松方程或者格林函数法,我真的无法想象你得有多痛苦。当然这部分内容同样是以我前面发布的文章为基础的,比如你会再次看到我在傅立叶变换里讲过的帕塞瓦尔等式,希望你还记得这是什么东西:)

 

2.3.6 希尔伯特空间


定义:在由内积所定义的范数意义下完备的内积空间称为希尔伯特(Hilbert)空间
希尔伯特空间是一类性质非常好的线性赋范空间,在工程上有着非常广泛的应用,而且在希尔伯特空间中最佳逼近问题可以得到比较完满的解决。

 

 

我的“图像处理中的数学原理”专栏中之系列文章已经以《图像处理中的数学修炼》为名结集出版(清华大学出版社)。该书详细介绍图像处理中的数学原理,为你打开一道通往图像世界的数学之门,详细内容及目录请见 http://blog.csdn.net/baimafujinji/article/details/48467225

 

 

 

没有更多推荐了,返回首页