基于块的纹理合成 matlab程序,基于块拼接的纹理合成算法

% 设置 MATLAB 初始环境.

clear; clc; randseed(0);

% 设置图像缝接的参数.

sourceTexture = './textures/examples/synthetic.jpg'; % 源纹理

blockSize = [15

15]; % 块大小

blockOverlap =

ceil(blockSize/5); % 块交迭大小

scaleFactor =

3; % scale factor for synthesized texture

preCache =

false; % enable/disable source block caching

errorTol =

0.1; % 容忍的匹配错误值

colorSpace =

'RGB'; % 选择的颜色空间是 RGB或者L*A*B

bndCut =

true; % 是否开启边界最优化

% 源纹理处理之前 (对合成的每个块进行摘录)

% 载入源纹理.

S.image =

im2double(imread(sourceTexture));

S.nrows =

size(S.image,1);

S.ncols =

size(S.image,2);

S.ncolors = size(S.image,3);

S.nblocks = [S.nrows S.ncols]-blockSize+1;

% 转换成 L*A*B* 颜色空间 (如果需要).

if strcmp(colorSpace,'LAB')

S.image =

applycform(S.image,makecform('srgb2lab'));

end

% 显示源纹理.

figure(1); clf;

if strcmp(colorSpace,'RGB')

imagesc(S.image);

else

imagesc(applycform(S.image,makecform('lab2srgb')));

end

set(gcf,'Name','Source Texture');

if S.ncolors == 1

colormap gray;

end

axis image; drawnow;

% 判断块大小和交迭块是否属于正常情况.

if (blockSize(1) > S.nrows) || (blockSize(2) > S.ncols)

error(['The output texture

does not support ',...

int2str(blockSize(1)),'x',int2str(blockSize(2)),' blocks.']);

end

if (blockOverlap(1) >= blockSize(1)) || (blockOverlap(2) >=

blockSize(2))

error(['A

',int2str(blockSize(1)),'x',int2str(blockSize(2)),...

' block does not support ',...

int2str(blockOverlap(1)),'x',int2str(blockOverlap(2)),'

overlaps.']);

end

%取源块的决定性区域.

[C,R] = meshgrid(1:S.nblocks(2),1:S.nblocks(1));

R = R'; C = C'; S.blockIndex = [R(:) C(:)];

blockRows = 0:blockSize(1)-1;

blockCols = 0:blockSize(2)-1;

% 连接源块到单个矩阵.

if preCache

S.blocks =

zeros(S.ncolors*prod(blockSize),prod(S.nblocks));

for i = 1:S.nblocks(1)

for j =1:S.nblocks(2)

S.blocks(:,(i-1)*S.nblocks(2)+j) = ...

reshape(S.image(i+blockRows,j+blockCols,:),[],1);

end

end

end

% 用图像缝合(IQ)产生一个合成纹理.

% Allocate storage for the synthesized output image/mask.

% 分配存储

% Note: Increase size to prevent clipped blocks at edges.

% 注意: 为保证省略块边缘增加大小,合成后修剪输出图像

% Remember to crop output image after synthesis.

% 合成后保存输出图像

overlapSize = blockSize-blockOverlap;

T.nrows =

overlapSize(1)*(ceil(scaleFactor*S.nrows/overlapSize(1))-1)+blockSize(1);

T.ncols =

overlapSize(2)*(ceil(scaleFactor*S.ncols/overlapSize(2))-1)+blockSize(2);

T.ncolors = S.ncolors;

T.nblocks = floor([T.nrows T.ncols]./overlapSize);

T.image =

zeros(T.nrows,T.ncols,T.ncolors);

T.mask =

false(T.nrows,T.ncols,T.ncolors);

T.cut = false(T.nrows,T.ncols);

% Determine locations of overlapping output blocks.

% 确定重叠输出块的位置

[C,R] = meshgrid(overlapSize(2)*(0:T.nblocks(2)-1)+1,...

overlapSize(1)*(0:T.nblocks(1)-1)+1);

R = R'; C = C'; T.blockIndex = [R(:) C(:)];

% Seed the synthesized texture with a random block.

i = ceil(prod(S.nblocks)*rand);

B = S.image(S.blockIndex(i,1)+blockRows,...

S.blockIndex(i,2)+blockCols,:);

T.image(blockRows+1,blockCols+1,:) = B;

T.mask(blockRows+1,blockCols+1,:) = 1;

% Display synthesized texture.

figure(2); clf;

if strcmp(colorSpace,'RGB')

imagesc(T.image);

else

imagesc(applycform(T.image,makecform('lab2srgb')));

end

set(gcf,'Name','Synthesized Texture');

if S.ncolors == 1

colormap gray;

end

axis image; drawnow;

% Synthesize the remaining pixels using IQ.

for i = 2:prod(T.nblocks)

% Extract current value

(i.e., the overlapping region).

B1 =

T.image(T.blockIndex(i,1)+blockRows,...

T.blockIndex(i,2)+blockCols,:);

M =

T.mask(T.blockIndex(i,1)+blockRows,...

T.blockIndex(i,2)+blockCols,:);

% 在源纹理中寻找相似块.

% Note: Either use pre-cached

blocks or draw from input texture.

B2 =

findmatch(S,B1,M,preCache,errorTol);

% Evaluate minimum boundary

cut.

if bndCut

[B3,C] = mincut(B1,B2,M);

T.cut(T.blockIndex(i,1)+blockRows,...

T.blockIndex(i,2)+blockCols,:) = ...

T.cut(T.blockIndex(i,1)+blockRows,...

T.blockIndex(i,2)+blockCols,:) | C;

else

B3 = B2;

end

% Update output/mask values

for this block.

T.image(T.blockIndex(i,1)+blockRows,...

T.blockIndex(i,2)+blockCols,:) = B3;

T.mask(T.blockIndex(i,1)+blockRows,...

T.blockIndex(i,2)+blockCols,:) = 1;

% Display synthesis

progress.

figure(2);

clf; % 关闭进程体

if

strcmp(colorSpace,'RGB')

imagesc(T.image);

else

imagesc(applycform(T.image,makecform('lab2srgb')));

end

axis image; drawnow;

end % End of IQ synthesis.

% 显示结果.

% Display synthesized texture.

figure(2); clf; %

关闭进程体

if strcmp(colorSpace,'RGB')

imagesc(T.image);

else

imagesc(applycform(T.image,makecform('lab2srgb')));

end

axis image; set(gcf,'Name','Synthesized Texture');

if S.ncolors == 1

colormap gray;

end

% 显示合成纹理的分割分界线.

if bndCut

figure(3); clf;

if

strcmp(colorSpace,'RGB')

I = T.image;

else

I = applycform(T.image,makecform('lab2srgb'));

end

if S.ncolors ==1

I = repmat(I,[1 1 3]);

end

R = I(:,:,1); G = I(:,:,2); B

= I(:,:,3);

idx = find(T.cut); R(idx) =

1*T.cut(idx) + 0*R(idx);

G(idx) = 0; B(idx) = 0;

I = cat(3,R,G,B); imagesc(I);

axis image;

axis image;

set(gcf,'Name','Synthesized Texture and Boundary Cuts');

end

findmatch.m：

function B =

findmatch(S,B1,M,preCache,errorTol)

% Vectorize search block (using input

mask).

blockRows = 0:size(B1,1)-1;

blockCols = 0:size(B1,2)-1;

B1 = B1(M);

% Perform exhaustive search to determine

similar blocks.

E = zeros(1,prod(S.nblocks));

if preCache

E =

sum((repmat(B1,1,prod(S.nblocks))-S.blocks(M,:)).^2,1);

else

for i =

1:prod(S.nblocks)

B2 =

S.image(S.blockIndex(i,1)+blockRows,S.blockIndex(i,2)+blockCols,:);

E(i) = sum((B1-B2(M)).^2);

end

end

% Randomly select output block among

matches.

% 随机选择一个匹配输出块

match_index = find(E <= (1+errorTol^2)*min(E));

i = ceil(length(match_index)*rand);

B = S.image(S.blockIndex(match_index(i),1)+blockRows,...

S.blockIndex(match_index(i),2)+blockCols,:);

mincut.m：

function [B,C] = mincut(B1,B2,M)

% Part I: Evaluate overlap error and extract horizontal/vertical

regions.

% 第一部分:

% Evaluate overlap error (using L2 norm).

E = M(:,:,1).*sum((B1-B2).^2,3);

% Extract horizontal and vertical overlap regions.

if M(1,end,1) == 1

blockOverlap(1) =

find(M(:,end,1),1,'last');

H =

E(1:blockOverlap(1),:);

else

H = [];

end

if M(end,1,1) == 1

blockOverlap(2) =

find(M(end,:,1),1,'last');

V =

E(:,1:blockOverlap(2));

else

V = [];

end

% Part II: minimum-error horizontal boundary cut (if

necessary).

% 第二部分:

% Evaluate minimum-error horizontal boundary cut.

if ~isempty(H)

% Recursively compute path

costs.

% 递归计算路径开销

CH = zeros(size(H));

CH(:,end) = H(:,end);

for j =

(size(H,2)-1):-1:1

for i = 1:size(H,1)

if i == 1

CH(i,j) = H(i,j) + min(CH(i+(0:1),j+1));

elseif i == size(H,1)

CH(i,j) = H(i,j) + min(CH(i+(-1:0),j+1));

else

CH(i,j) = H(i,j) + min(CH(i+(-1:1),j+1));

end

end

end

end % End of horizontal boundary cut.

% Part III: minimum-error vertical boundary cut (if

necessary).

% 第三部分:

% Evaluate minimum-error vertical boundary cut.

if ~isempty(V)

% Recursively compute path

costs.

% 递归计算路径开销

CV = zeros(size(V));

CV(end,:) = V(end,:);

for i =

(size(V,1)-1):-1:1

for j = 1:size(V,2)

if j == 1

CV(i,j) = V(i,j) + min(CV(i+1,j+(0:1)));

elseif j == size(V,2)

CV(i,j) = V(i,j) + min(CV(i+1,j+(-1:0)));

else

CV(i,j) = V(i,j) + min(CV(i+1,j+(-1:1)));

end

end

end

end % End of vertical boundary cut.

% Part IV: Evaluate total minimum-error boundary cut (if

necessary).

% 第四部分:

% Allocate storage for image mask (i.e., filled

boundaries).

C = false(size(E));

M = false(size(B1));

% Determine total minimum-error boundary cut.

% 确定总最小误差分割分界线

% Case 1: Only horizonal overlap (e.g., along first scan

column).

if (~isempty(H) && isempty(V))

% Trace cost matrix to

determine boundary cut.

[ignore,i] =

min(CH(:,1));

C(i,1) = 1; M(i:end,1,:) =

1;

for j = 2:size(H,2)

if i == 1

[ignore,i] = min(CH(i+(0:1),j));

elseif i == size(H,1)

[ignore,di] = min(CH(i+(-1:0),j));

i = i+(di-2);

else

[ignore,di] = min(CH(i+(-1:1),j));

i = i+(di-2);

end

C(i,j) = 1; M(i:end,j,:) = 1;

end

% Case 2: Only vertical overlap (e.g., along first scan row).

elseif (isempty(H) && ~isempty(V))

% Trace cost matrix to

determine boundary cut.

[ignore,j] =

min(CV(1,:));

C(1,j) = 1; M(1,j:end,:) =

1;

for i = 2:size(V,1)

if j == 1

[ignore,j] = min(CV(i,j+(0:1)));

elseif j == size(V,2)

[ignore,dj] = min(CV(i,j+(-1:0)));

j = j+(dj-2);

else

[ignore,dj] = min(CV(i,j+(-1:1)));

j = j+(dj-2);

end

C(i,j) = 1; M(i,j:end,:) = 1;

end

% Case 3: Both vertical and horizontal overlaps.

else

% Determine pixel where

vertical/horizontal cuts intersect.

CC =

CH(1:blockOverlap(1),1:blockOverlap(2)) + ...

CV(1:blockOverlap(1),1:blockOverlap(2)) - ...

E(1:blockOverlap(1),1:blockOverlap(2));

% Invert pixel assignment

mask.

M = ~M;

% Trace back from intersection

to determine horizontal boundary cut.

[ignore,p] = min(CC(:)); [i,j]

= ind2sub(blockOverlap,p);

C(i,j) = 1; M(1:i,j,:) =

0;

for j = (j+1):size(H,2)

if i == 1

[ignore,i] = min(CH(i+(0:1),j));

elseif i == size(H,1)

[ignore,di] = min(CH(i+(-1:0),j));

i = i+(di-2);

else

[ignore,di] = min(CH(i+(-1:1),j));

i = i+(di-2);

end

C(i,j) = 1; M(1:i,j,:) = 0;

end

% Trace back from intersection

to determine vertical boundary cut.

[ignore,p] = min(CC(:)); [i,j]

= ind2sub(blockOverlap,p);

C(i,j) = 1; M(i,1:j,:) =

0;

for i = (i+1):size(V,1)

if j == 1

[ignore,j] = min(CV(i,j+(0:1)));

elseif j == size(V,2)

[ignore,dj] = min(CV(i,j+(-1:0)));

j = j+(dj-2);

else

[ignore,dj] = min(CV(i,j+(-1:1)));