Matlab 蒙太奇马赛克拼图（小图拼成大图）

0. 介绍

马赛克照片拼图的原则是用小的图片拼成一张大的图片，特点是保留了大图的信息，放大后仍能看到每一张小图的内容。本文把马赛克照片拼图分为3类，分别为图1的框架型、图2的叠加型和图3的组合型。

图1. 框架型拼图

图2. 叠加型

图3. 组合型

框架型只关注图像主体的形状，不关心形状的填充颜色，如图1的心形的颜色不是重点。叠加型保留大图作为背景，小图依次叠加在大图上。按我的理解，“叠加型”并不能算作马赛克拼图，因为，首先没有对图像进行马赛克处理；其次没有对图像进行拼接；最后小图仅仅是以一定比例覆盖在大图上，小图的颜色不影响大图的内容。最后一种是真正意义上的马赛克拼图，每一张小图的主体颜色和大图的一个像素块的颜色相近。

人脸数据集来自于香港中文大学，地址http://mmlab.ie.cuhk.edu.hk/projects/TCDCN.html

1. 框架型

框架型的实现方法是：

1. 对图像进行马赛克处理
2. 在相应的网格内填充图片

首先定义马赛克处理函数

% mosaic.m
function img = mosaic(img, gridRow, gridColumn) 
    [row, column, channel] = size(img);
    for i = floor(gridRow/2)+1 : gridRow : row 
        for j = floor(gridColumn/2)+1 : gridColumn : column
            boundRow = i + floor(gridRow/2);
            boundColumn = j + floor(gridColumn/2);
            if boundRow > row
                boundRow = row;
            end
            if boundColumn > column
                boundColumn = column;
            end
            if channel == 1
                img(i-floor(gridRow/2):boundRow, j-floor(gridColumn/2):boundColumn) = img(i,j);
            else
                startRow = i - floor(gridRow/2);
                startColumn = j-floor(gridColumn/2);
                img(startRow:boundRow, startColumn:boundColumn, 1) = img(i, j, 1);
                img(startRow:boundRow, startColumn:boundColumn, 2) = img(i, j, 2);
                img(startRow:boundRow, startColumn:boundColumn, 3) = img(i, j, 3);
            end
        end
    end
end

然后依次读取图片文件和二值化。二值化时要选择合适的阈值，使得图像主体与背景尽量分开。最后进行马赛克处理，效果图如下：

小图填充的完整代码如下

% matlab
clear;

img = imread('heart.jpg');
img = rgb2gray(img);
img = img < 240;
fillImgRow = 25;
fillImgColumn = 25; % 填充小图的行和列的大小
img = mosaic(img, fillImgRow, fillImgColumn); % 根据设定的行和列大小进行马赛克处理
[row, column] = size(img);

% 首先生成和原图大小一致，具有3通道的矩阵，背景默认是255（即白色）
% 注意： outputImg的类型是double，imshow时把像素值当做[0,1]处理
outputImg = 255*ones(row, column, 3);

d = dir('MTFL/AFLW'); % 人脸数据集的地址
numImgs = length(d);

for i = 1 : fillImgRow : row
    for j = 1 : fillImgColumn : column
        if img(i, j) % 如果是需要填充的位置
        	% 随机取出一张图片，如果不是3通道，重新选取
            subimg = imread(['MTFL/AFLW/' d(int32((numImgs-3)*rand()+3)).name]);
            while length(size(subimg))<3
                subimg = imread(['MTFL/AFLW/' d(int32((numImgs-2)*rand()+3)).name]);
            end
            % 修改取出图片的大小，设置为设定值
            subimg = imresize(subimg, [fillImgRow, fillImgColumn]);
            r = fillImgRow;
            c = fillImgColumn;
            % 到达图像边界时，切割小图
            if i + fillImgRow - 1 > row
                r = row-i+1;
                subimg = subimg(1:r, :, :);
            end
            if j + fillImgColumn - 1 > column
                c = column-j+1;
                subimg = subimg(:, 1:c, :);
            end
            % 分别赋值
            outputImg(i:i+r-1, j:j+c-1, 1) = subimg(:, :, 1);
            outputImg(i:i+r-1, j:j+c-1, 2) = subimg(:, :, 2);
            outputImg(i:i+r-1, j:j+c-1, 3) = subimg(:, :, 3);
        end
    end
end
% 将outputImg的类型转换为uint8，使图片正常显示
outputImg = uint8(outputImg);
imwrite(outputImg, 'mosaicJigsawHeart.jpg');

效果图：

2. 叠加型

叠加型不需要对图像进行马赛克处理，仅仅要求两张图片以一定权重叠加。假设有两张大小相同的图片$I_1(x,y)$和$I_2(x,y)$，并决定$I_1(x,y)$以$w$的权重添加到输出图像$I(x,y)$中，则有
$$I(x,y)=w\times I_1(x,y)+(1-w)\times I_2(x,y)$$输出图像将能观察到$I_1(x,y)$和$I_2(x,y)$的图像内容。
在这个例子中，本文采用大小不一致的图片叠加，这意味着小图不再是间隔一定像素规则地叠加到背景图上。这里采用掩模版的方法：创建一个大小和大图任意一个通道一致的矩阵，值取0或1。遍历掩模版上所有的点，当$mask(i,j)=1$时，从$(i,j)$开始取出大小和小图一致的掩模版，并与小图相乘，最后以一定权重和背景图叠加并更新掩模版。下图是$(500,500)$时的掩模版，

可以看到蓝色部分是不平整的。
掩模版使用图解如下：

1. mask的实时情况如左图所示，黑色表示已经被使用过(等于0），现在发现绿色点的位置$mask(i,j)=1$，可以看到从绿色点向行和列递增方向的可使用空间都不规则，无法满足小图的填充
   
2. 根据小图的尺寸，取出对应大小的掩模版，并与之相乘，最后与背景图叠加，相当于队列操作，最先进入的显示在最上层，越后面显示越底层。
   
   完整代码如下：

% matlab

clear;

minRow = 30; % 小图最小行数
maxRow = 50; % 小图最大行数，如果小图行数大于最大行数，则调整大小
weight = 0.3; % 小图的权重

img = imread('photo.jpg');
% 把img转化为带权重的double型。
% img后面与小图可以直接相加，避免多个小图重叠处，
% img多次乘以权重后接近于0 （即可能出现0.7*0.7*0.7*img)
img = (1-weight)*double(img);
[row, column, ~] = size(img);
mask = uint8(ones(row, column));
d = dir('MTFL/AFLW');
numImgs = length(d);
for i = 1 : row
    for j = 1 : column
        if mask(i, j)
            subimg = imread(['MTFL/AFLW/' d(int32((numImgs-3)*rand()+3)).name]);
            sz = size(subimg);
            while length(sz) < 3
                subimg = imread(['MTFL/AFLW/' d(int32((numImgs-3)*rand()+3)).name]);
                sz = size(subimg);
            end
            r = sz(1);
            % 如果小图的行数大于50，则随机调整行数在[30,50]范围内
            if r > maxRow 
                subimg = imresize(subimg, ((maxRow-minRow)*rand()+minRow)/r);
            end
            
            [r, c, ~] = size(subimg);
            if i+r-1 > row
                r = row-i+1;
                subimg = subimg(1:r,:,:);
            end
            if j+c-1 > column
                c = column-j+1;
                subimg = subimg(:, 1:c,:);
            end
                
            submask = mask(i:i+r-1, j:j+c-1);
            subimg(:,:,1) = subimg(:,:,1) .* submask;
            subimg(:,:,2) = subimg(:,:,2) .* submask;
            subimg(:,:,3) = subimg(:,:,3) .* submask;
            img(i:i+r-1, j:j+c-1, :) = weight*double(subimg) + img(i:i+r-1, j:j+c-1, :);
            mask(i:i+r-1, j:j+c-1) = 0;
        end
    end
end
img = uint8(img);
imwrite(img, 'output.jpg');

效果图：

3. 组合图

程序所用素材

组合图完全由图片拼合而成，大图的颜色由小图的主体颜色决定。因此在一个位置，像素块不再是由随机选取的小图填充，而是由大量图片中颜色最接近像素块的小图填充。
第一步，对所有图片按照颜色进行分类。对颜色进行分类至少有两种方法——平均分类和区间分类。平均分类是分别对图像的RGB进行平均处理得到$(R_{avg},G_{avg},B_{avg})$，按照$(R_{avg},G_{avg},B_{avg})$分类。区间分类是分别统计图像RGB灰度值的分布，选择包含像素最多的灰度$(R_m,G_m,B_m)$作为该图像的代表，以$(R_m,G_m,B_m)$对图像分类。本文采用区间分类的方法。
**备注：**图片素材数量很少，可以通过调整素材的RGB直方图来达到不同的颜色，这种方法实现更简单而且得到的结果颜色上的效果更好，但是素材单调。

马赛克处理效果如下图所示

3.1 resize

首先要重新调整数据集中的所有图像的大小，编写函数datasetReisze.m：

% datasetResize.m
function numImgs = datasetResize(fromFolder, toFolder, row, column)
% fromFolder: 图像数据集目录地址
% toFolder：调整大小后存放图像数据集的目录
% row, column：图像要调整的大小
% numImgs：函数返回图像总数
    d = dir(fromFolder);
    n = length(d);
    numImgs = 0;
    for i = 3 : n
        img = imread([fromFolder, '/', d(i).name]);
        sz = size(img);
        if length(sz) < 3 || sz(1) < row || sz(2) < column
            continue;
        end
        img = imresize(img, [row, column]);
        imwrite(img, [toFolder, '/', num2str(numImgs+1), '.jpg']);
        numImgs = numImgs + 1;
    end
end

3.2 区间分类

3.2.1 灰度级映射

首先编写生成单个图像三个分量的灰度分布的函数countRGB.m。根据传递参数interval计算出落入在每个区间的像素数，计算得到的区间具有如下形式：

1~interval, interval+1~2*interval, ···, ~256
记为区间1, 2, ···, ceil(256/interval)

因为matlab的索引从1开始，灰度级由0~255调整为1~256。具体实现

% countRGB.m
function [R, G, B] = countRGB(img, interval)
	% img：输入图像
	% interval：灰度区间大小
	% 返回值R,G,B是向量，表示落入到每一个区间的像素数目
    r = zeros(256, 1);
    g = zeros(256, 1);
    b = zeros(256, 1);
    [row, column, ~] = size(img);
    for i = 1 : row
        for j = 1 : column
            r(img(i, j, 1) + 1) = r(img(i, j, 1) + 1) + 1;
            g(img(i, j, 2) + 1) = g(img(i, j, 2) + 1) + 1;
            b(img(i, j, 3) + 1) = b(img(i, j, 3) + 1) + 1;
        end
    end
    len = ceil(256/interval);
    R = zeros(len, 1);
    G = zeros(len, 1);
    B = zeros(len, 1);
    for i = 1 : len
        startIndex = (i-1)*interval + 1;
        endIndex   = i*interval;
        if endIndex > 256
            endIndex = 256;
        end
        R(i) = sum(r(startIndex: endIndex));
        G(i) = sum(g(startIndex: endIndex));
        B(i) = sum(b(startIndex: endIndex));
    end
end

接下来统计图像每一个分量中哪一个区间的灰度占主导（即像素数目最大），编写函数dominantRGB.m统计数据集中所有图像的“主导”灰度级。

% dominantRGB.m
function dominantRGB = countDatasetRGB(addr, interval)
% addr：图像数据集的目录地址
% interval：灰度区间大小
% 返回值dominantRGB是一个包含RGB分量主导区间的矩阵
% dominantRGB的行是图像的编号，图像的命名是以1~numImgs的方式进行
% 每行的第一个元素表示R分量的主导区间，以此类推
% 
    %  #img    dominantValueR dominantValueG dominantValueB
    %    1    [       8             4             3  ]
    %    2    [       11            1             1  ]
    %         [        .........
    % numImgs [                                      ]
    
    d = dir(addr); 
    numImgs = length(d)-2;
    dominantRGB = zeros(numImgs, 3);
    %dominantRGB(:, 1) = 1: numImgs;
    for i = 1 : numImgs
        img = imread([addr, '/', num2str(i), '.jpg']);
        [R, G, B] = countRGB(img, interval);
        [val, index] = max([R, G, B]);
        dominantRGB(i, :) = index;
        %if mod(i, 100)==0
        %    disp(i)
        %end
    end
end

最后定义一个函数，功能是生成每一个分量中主导灰度落在每一个灰度区间的图像编号的矩阵。例如红色分量的灰度区间2，表示灰度范围在interval+1~3*interval，如果一个图像的红色分量的主导灰度属于区间2，则把该图像的编号添加到红色分量的区间2中。编写数据生成函数mapDatasetRGB.m

% mapDatasetRGB.m
function [R, G, B] = mapDatasetRGB(dominantRGB, interval)
    level = ceil(256/interval);
    R = cell(level, 1);
    G = cell(level, 1);
    B = cell(level, 1);
    [n, ~] = size(dominantRGB);
    for i = 1 : n
        rgb = dominantRGB(i, :);
        r = rgb(1);
        g = rgb(2);
        b = rgb(3);
        R{r} = [R{r}, i];
        G{g} = [G{g}, i];
        B{b} = [B{b}, i];
    end
end

函数返回值R、G或者B的数据类型是cell。以R为例，R的大小是ceil(256/interval)，R{i}的数据是向量，表示红色分量主导灰度区间是i的图像的集合，如 R { i } = [ 1 , 200 , 3001 , 5000 ] R\{i\} = [1, 200, 3001, 5000] R{i}=[1,200,3001,5000]代表1.jpg、200.jpg、3001.jpg和5000.jpg的红色分量落在第$i$个区间的像素最多。

3.2.2 根据RGB选取合适的图像

根据原图像的RGB，计算该RGB落入的灰度区间，从mapDatasetRGB的返回向量中，找到对应颜色近似的图像。假设该像素的RGB属于的灰度区间分别是$l,m,n$，那么，颜色与该像素接近的图像编号的集合可以通过取交集的方法得到 I = R { l } ∩ G { m } ∩ B { n } I=R\{l\}\cap G\{m\}\cap B\{n\} I=R{l}∩G{m}∩B{n}，最后从$I$集合里随机选出一张图像。如果$I$是空集，扩大区间的搜索范围， I = R { l − 1 : l + 1 } ∩ G { m − 1 : m + 1 } ∩ B { n − 1 : n + 1 } I=R\{l-1:l+1\}\cap G\{m-1:m+1\}\cap B\{n-1:n+1\} I=R{l−1:l+1}∩G{m−1:m+1}∩B{n−1:n+1}如果$I$还是空集，重复上述操作，直到灰度范围覆盖所有灰度取值。编写chooseImage.m实现上述功能。

% chooseImage.m
function [id, n, ii] = chooseImage(R, G, B, Rset, Gset, Bset)
    r = Rset{R};
    g = Gset{G};
    b = Bset{B};
    id = intersect(intersect(r, g), b);
    ii = 1;
    len = length(Rset);
    while isempty(id)
        sR = R - ii;
        eR = R + ii;
        sG = G - ii;
        eG = G + ii;
        sB = B - ii;
        eB = B + ii;
        if sR < 1
            sR = 1;
        end
        if eR > len
            eR = len;
        end
        if sG < 1
            sG = 1;
        end
        if eG > len
            eG = len;
        end
        if sB < 1
            sB = 1;
        end
        if eB > len
            eB = len;
        end
        for i = sR : eR
            r = union(r, Rset{i});
        end
        for i = sG : eG
            g = union(g, Gset{i});
        end
        for i = sB : eB
            b = union(b, Bset{i});
        end
        id = intersect(intersect(r, g), b);
        ii = ii + 1;
        if sR == 1 && sG == 1 && sB == 1 && eR == len && eG == len && eB == len
            break;
        end
    end
    n = length(id);
    %disp(n);
    %disp(ii);
    id = id(floor(n*rand()+1));
end

显然当图像很大时，chooseImage函数时间效率很低，因为对每一个像素块都要进行多次判断、取交集的操作。如果图像的分辨率是40004000，马赛克取4040像素，那么将有10000个像素块需要调用chooseImage函数进行判断。对于这种情况，牺牲空间来节省时间十分有效。每次出现一种RGB组合$(r,g,b)$时，把它添加到字典结构中，调用chooseImage，选取$(r,g,b)$对应的图像。下次再出现$(r,g,b)$时，直接选取上次的图像即可。

// 伪代码
if (r, g, b) in RGBDict
	subimg = RGBDict[(r, g, b)]
else
	subimg = chooseImage((r, g, b))
	RGBDict[(r, g, b)] = subimg

即使所有的RGB组合都出现，以灰度区间大小20为例，总共的组合方式有131313=2197。与10000次调用相比，2197次调用小得多。

3.2.3 灰度分类主程序

% matlab
clear;
rowInter = 40;
columnInter = 40; % 马赛克大小
levelInterval = 20; % 灰度区间大小
img = imread('pic.jpg');
img = imresize(img, 3000/1686); % 如果图像较小，适当用imsize插值
img = mosaic(img, rowInter, columnInter);
[row, column, ~] = size(img);

% 以下代码的作用是判断目录中是否有数据，如果有跳过或者load，如果没有，则计算
addr = ['image-dataset_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist(addr, 'dir')
    mkdir(addr);
    datasetResize('MTFL/AFLW', addr, rowInter, columnInter);
end
dominantRGBFilename = ['dominantRGB_', num2str(rowInter), 'x', num2str(columnInter), '_', num2str(levelInterval)];
if ~exist([dominantRGBFilename, '.mat'], 'file')
    dominantRGB = countDatasetRGB(addr, levelInterval);
    save(dominantRGBFilename, 'dominantRGB', 'levelInterval');
else
    load(dominantRGBFilename);
end
RGBMapFilename = ['RGBMap_', num2str(rowInter), 'x', num2str(columnInter), '_', num2str(levelInterval)];
if ~exist([RGBMapFilename, '.mat'], 'file')
    [Rmap, Gmap, Bmap] = mapDatasetRGB(dominantRGB, levelInterval);
    save(RGBMapFilename, 'Rmap', 'Gmap', 'Bmap', 'levelInterval');
else
    load(RGBMapFilename);
end

for i = 1 : rowInter : row
    for j = 1 : columnInter : column
        Rlevel = ceil((double(img(i, j, 1)) + 1)/levelInterval);
        Glevel = ceil((double(img(i, j, 2)) + 1)/levelInterval);
        Blevel = ceil((double(img(i, j, 3)) + 1)/levelInterval);
        subimg = imread([addr, '/', num2str(chooseImage(Rlevel, Glevel, Blevel, Rmap, Gmap, Bmap)), '.jpg']);
        rowBound = i + rowInter - 1;
        columnBound = j + columnInter - 1;
        if rowBound > row
            rowBound = row;
            subimg = subimg(1:rowBound-i+1, :, :);
        end 
        if columnBound > column
            columnBound = column;
            subimg = subimg(:, 1:columnBound-j+1, :);
        end
        img(i:rowBound,j:columnBound,:) = subimg;
    end
end

3.2.4 结果

实现效果图1：
图像重新调整大小行数等于2500，马赛克大小是3030，图像数据集共有2967张

实现效果图2：
图像重新调整大小行数等于3000，马赛克大小是4040，图像数据集共有16414张

感官上，第二张的效果比第一张好，原因是数据集更大了，存在的RGB组合的可能性也更大。
实现效果图3：
图像重新调整大小行数等于3000，马赛克大小是2020，灰度区间大小10，图像数据集共有16414张
实现效果图4：极端情形
图像重新调整大小行数等于3000，马赛克大小是2020，灰度区间大小1，图像数据集共有16414张

3.3 平均分类

3.3.1 定义函数

平均分类与上面的区间分类基本一致，主要区别是从统计每一张图片出现最多的像素区间改为统计该图像的平均灰度$$av{g}_K=\sum _i\sum _j{I}_K(i,j)$$其中，$K=R,G,B$。
对应区间分类中的countRGB、countDatasetRGB、mapDatasetRGB和chooseImage函数，平均分类中定义calcAvgRGB、calcDatasetAvgRGB、mapDatasetAvgRGB和chooseImageAvg函数。

% calcAvgRGB.m
function [avgR, avgG, avgB] = calcAvgRGB(img)
    [row, column, ~] = size(img);
    img = int32(img);
    numPixels = row * column;
    avgR = uint8(sum(sum(img(:, :, 1)))/numPixels);
    avgG = uint8(sum(sum(img(:, :, 2)))/numPixels);
    avgB = uint8(sum(sum(img(:, :, 3)))/numPixels);
end

% calcDatasetAvgRGB.m
function avgRGB = calcDatasetAvgRGB(addr)
    d = dir(addr); 
    numImgs = length(d)-2;
    avgRGB = zeros(numImgs, 3);
    for i = 1 : numImgs
        img = imread([addr, '/', num2str(i), '.jpg']);
        avgRGB(i, :) = calcAvgRGB(img);
        if mod(i, 100) == 0
            disp(i/numImgs);
        end
    end
end

% mapDatasetAvg.m
function [Rmap, Gmap, Bmap] = mapDatasetAvgRGB(avgRGB)
    Rmap = cell(256, 1);
    Gmap = cell(256, 1);
    Bmap = cell(256, 1);
    [n, ~] = size(avgRGB);
    for i = 1 : n
        Rmap{avgRGB(i, 1)} = [Rmap{avgRGB(i, 1)}, i];
        Gmap{avgRGB(i, 2)} = [Gmap{avgRGB(i, 2)}, i];
        Bmap{avgRGB(i, 3)} = [Bmap{avgRGB(i, 3)}, i];
    end
end

% chooseImageAvg.m
function [id, n, ii] = chooseImageAvg(R, G, B, Rmap, Gmap, Bmap)
    r = Rmap{R+1};
    g = Gmap{G+1};
    b = Bmap{B+1};
    id = intersect(intersect(r, g), b);
    ii = 1;
    len = length(Rmap);
    while isempty(id)
        sR = R - ii;
        eR = R + ii;
        sG = G - ii;
        eG = G + ii;
        sB = B - ii;
        eB = B + ii;
        if sR < 1
            sR = 1;
        end
        if eR > len
            eR = len;
        end
        if sG < 1
            sG = 1;
        end
        if eG > len
            eG = len;
        end
        if sB < 1
            sB = 1;
        end
        if eB > len
            eB = len;
        end
        for i = sR : eR
            r = union(r, Rmap{i});
        end
        for i = sG : eG
            g = union(g, Gmap{i});
        end
        for i = sB : eB
            b = union(b, Bmap{i});
        end
        id = intersect(intersect(r, g), b);
        ii = ii + 1;
        if sR == 1 && sG == 1 && sB == 1 && eR == len && eG == len && eB == len
            break;
        end
    end
    n = length(id);
    %disp(n);
    %disp(ii);
    id = id(floor(n*rand()+1));

同样地，chooseImageAvg函数的效率很低，需要优化。

3.3.2 平均分类主程序

平均主程序和区间分类的主程序大体一致，不过是把相应的函数替换。

% matlab
% mosaic jigsaw combination average

clear;

rowInter = 20;
columnInter = 20;
img = imread('pic.jpg');
img = imresize(img, 3000/1686);
img = mosaic(img, rowInter, columnInter);
[row, column, ~] = size(img);

addr = ['image-dataset_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist(addr, 'dir')
   mkdir(addr);
   datasetResize('MTFL/AFLW', addr, rowInter, columnInter);
end
avgRGBFilename = ['avgRGB_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist([avgRGBFilename, '.mat'], 'file')
   avgRGB = calcDatasetAvgRGB(addr);
   save(avgRGBFilename, 'avgRGB');
else
   load(avgRGBFilename);
end
RGBMapFilename = ['RGBMap_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist([RGBMapFilename, '.mat'], 'file')
   [Rmap, Gmap, Bmap] = mapDatasetAvgRGB(avgRGB);
   save(RGBMapFilename, 'Rmap', 'Gmap', 'Bmap');
else
   load(RGBMapFilename);
end
for i = 1 : rowInter : row
   for j = 1 : columnInter : column
       subimg = imread([addr, '/', num2str(chooseImageAvg(img(i, j, 1), img(i, j, 2), img(i, j, 3), Rmap, Gmap, Bmap)), '.jpg']);
       rowBound = i + rowInter - 1;
       columnBound = j + columnInter - 1;
       if rowBound > row
           rowBound = row;
           subimg = subimg(1:rowBound-i+1, :, :);
       end 
       if columnBound > column
           columnBound = column;
           subimg = subimg(:, 1:columnBound-j+1, :);
       end
       img(i:rowBound,j:columnBound,:) = subimg;
   end
   if mod(i, 100) == 0
          disp(i/row);
   end
end

3.3.3 结果

图像重新调整大小行数等于3000，马赛克大小是20*20，图像数据集共有16414张

直观上来看，RGB平均比区间的效果更好。

3.4 直方图平移

当可选的图像数据不够时，填充图像的颜色无法和原图像素颜色很好地匹配，得到的拼图效果也不好。此时，可以通过对RGB的直方图平移一个量，得到的新图像的颜色和原图像素颜色完全一致（RGB平均值完全一致）。

3.4.1 方法

n e w S u b i m g K = s u b i m g K + ( a v g { I K } − a v g { s u b i m g K } ) newSubimg_K = subimg_K + (avg\{I_K\}-avg\{subimg_K\}) newSubimgK​=subimgK​+(avg{IK​}−avg{subimgK​})

3.4.2 直方图平移主程序

% matlab
rowInter = 50;
img = imread('pic.jpg');
[row, ~, ~] = size(img);
block = imresize(img, rowInter/row);
[rowInter, columnInter, ~] = size(block);
img = imresize(img, 3000/1686);
img = mosaic(img, rowInter, columnInter);
[row, column, ~] = size(img);

[blockAvgR, blockAvgG, blockAvgB] = calcAvgRGB(block);

for i = 1 : rowInter : row
    for j = 1 : columnInter : column
        rowBound = i + rowInter - 1;
        columnBound = j + columnInter - 1;
        %[imgAvgR, imgAvgG, imgAvgB] = calcAvgRGB(img(i:rowBound, j:rowBound, :));
        subimg = block;
        subimg(:,:,1) = block(:,:,1) + img(i,j,1) - blockAvgR;
        subimg(:,:,2) = block(:,:,2) + img(i,j,2) - blockAvgG;
        subimg(:,:,3) = block(:,:,3) + img(i,j,3) - blockAvgB;
        if rowBound > row
            rowBound = row;
            subimg = subimg(1:rowBound-i+1, :, :);
        end 
        if columnBound > column
            columnBound = column;
            subimg = subimg(:, 1:columnBound-j+1, :);
        end
        img(i:rowBound,j:columnBound,:) = subimg;
    end
    if mod(i, 100) == 0
           disp(i/row);
    end
end

3.4.3 结果

效果如下：
图像重新调整大小行数等于3000，马赛克大小是50*50，填充图像就是原图经过压缩后端图像。

3.5 平均分类和直方图平移混合

3.5.1 介绍

这个操作是将平均分类和直方图平移结合在一起，既要考虑到尽可能保留填充图像的RGB信息，又要最大程度得满足平均RGB匹配。

3.5.2 程序

代码基本和平均分类的代码一样，只是在循环中多加了一段RGB加减的代码，但是要注意数据类型转换，因为uint8类型的范围是0-255，而img(i,j,1)-subimgAvgR可能小于0。

% matlab
% img = int32(img);
% subimg = int32(subimg);
[subimgAvgR, subimgAvgG, subimgAvgB] = calcAvgRGB(subimg);
subimg(:,:,1) = subimg(:,:,1) + img(i,j,1) - int32(subimgAvgR);
subimg(:,:,2) = subimg(:,:,2) + img(i,j,2) - int32(subimgAvgG);
subimg(:,:,3) = subimg(:,:,3) + img(i,j,3) - int32(subimgAvgB);

完整代码如下：

% matlab
rowInter = 20;
columnInter = 20;
img = imread('pic.jpg');
img = imresize(img, 3000/1686);
img = mosaic(img, rowInter, columnInter);
[row, column, ~] = size(img);

addr = ['image-dataset_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist(addr, 'dir')
    mkdir(addr);
    datasetResize('dataset/MTFL/AFLW', addr, rowInter, columnInter);
end
avgRGBFilename = ['avgRGB_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist([avgRGBFilename, '.mat'], 'file')
    avgRGB = calcDatasetAvgRGB(addr);
    save(avgRGBFilename, 'avgRGB');
else
    load(avgRGBFilename);
end
RGBMapFilename = ['RGBMap_', num2str(rowInter), 'x', num2str(columnInter)];
if ~exist([RGBMapFilename, '.mat'], 'file')
    [Rmap, Gmap, Bmap] = mapDatasetAvgRGB(avgRGB);
    save(RGBMapFilename, 'Rmap', 'Gmap', 'Bmap');
else
    load(RGBMapFilename);
end

img = int32(img);

for i = 1 : rowInter : row
    for j = 1 : columnInter : column
        subimg = imread([addr, '/', num2str(chooseImageAvg(img(i, j, 1), img(i, j, 2), img(i, j, 3), Rmap, Gmap, Bmap)), '.jpg']);
        subimg = int32(subimg);
        rowBound = i + rowInter - 1;
        columnBound = j + columnInter - 1;
        [subimgAvgR, subimgAvgG, subimgAvgB] = calcAvgRGB(subimg);
        subimg(:,:,1) = subimg(:,:,1) + img(i,j,1) - int32(subimgAvgR);
        subimg(:,:,2) = subimg(:,:,2) + img(i,j,2) - int32(subimgAvgG);
        subimg(:,:,3) = subimg(:,:,3) + img(i,j,3) - int32(subimgAvgB);
        if rowBound > row
            rowBound = row;
            subimg = subimg(1:rowBound-i+1, :, :);
        end 
        if columnBound > column
            columnBound = column;
            subimg = subimg(:, 1:columnBound-j+1, :);
        end
        img(i:rowBound,j:columnBound,:) = subimg;
    end
    if mod(i, 100) == 0
           disp(i/row);
    end
end
img = uint8(img);

3.5.3 结果

图像重新调整大小行数等于3000，马赛克大小是20*20，图像数据集共有16414张

3.6 总结

马赛克照片拼图的程序粗略完成，效果勉强能看。仅从马赛克大小20*20的几张图来看，区间分类更擅长还原原图丰富的颜色，平均分类更能反映原图的主体，但是颜色暗淡。所以，两个程序都需要优化，使得生成的图片能尽可能还原原图内容。当把平均分类和直方图平移结合在一起使用时，效果当然是非常不错的。另外，程序只考虑到结果，未考虑时间效率，有多处重复（如多次调用calcAvgRGB、数据类型转换），所以仍然需要修改，使代码更加简洁，逻辑更加清晰。

备注：以上内容仅仅是根据自己的想法写得。