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

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

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

在后半部分的图像处理算法介绍部分，偶尔为了便于理解，本书配备了必要的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, []);

特别说明： 如果本书的读者需要函数clipHistogram和mapHistogram的具体实现，可以从“图像处理书籍读者群”中直接获取。你也可以参考MATLAB中CLAHE的实现代码中相关部分来自行编写这部分内容。

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);

特别说明：函数getRadiance和guidedfilter_color的实现可以从“图像处理书籍读者群”中直接获取（群号见文末）。你也可以直接使用MATLAB中提供的内置导向滤波函数来替代guidedfilter_color，这样实现代码会更加简洁！

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);

本书中用于实验的图片下载【百度网盘】

本书源起

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

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

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