MATLAB学习笔记 学习总结归纳（第二周）

本次学习内容为第6章

这周复习了上一周所学的内容，顽固加深理解以后，学习了第6章前半部分。

MATLAB中的彩色图像的表示（图像格式）

RGB图像

所谓RGB图像，即为R（Red 红） G（Green绿） B（Blue蓝） 的灰度图（三个长宽一样的二维矩阵）组成的图像比如以下这段代码

r = zeros(300, 300);
g = zeros(300, 300);
b = zeros(300, 300);
r(1:100, 1:100) = 1;
g(101:200, 101:200) = 1;
b(201:300, 201:300) = 1;
img = cat(3, r, g, b);
subplot(2, 2, 1), imshow(r), title('R');
subplot(2, 2, 2), imshow(g), title('G');
subplot(2, 2, 3), imshow(b), title('B');
subplot(2, 2, 4), imshow(img), title('RGB');

输入:



输出:



可以看出 其实RGB 三张分别代表红绿蓝的灰度图 叠加在一起的图像。

索引图像

索引图像有两个分量: 一个整数数据矩阵X 和 一个彩色银蛇矩阵map。矩阵map是一个大小为m*3的double类数组，其值是区间[0， 1]上的浮点数。map的长度m等于其定义的颜色的个数。map存的是RGB三个分量。

f  = imread('onion.png');
[X, map] = rgb2ind(f, 2);
[X2, map2] = rgb2ind(f, 256);
[X3, map3] = rgb2ind(f, 65535);
subplot(2, 2, 1), imshow(f), title('原图');
subplot(2, 2, 2), imshow(X, map), title('索引图 2种颜色');
subplot(2, 2, 3), imshow(X2, map2), title('索引图 256种颜色');
subplot(2, 2, 4), imshow(X3, map3), title('索引图 65535种颜色');

输入:



输出:

处理RGB和索引图像的函数

* ind2rgb(X, map) 将索引图变为RGB图*

rgb2ind(image, n) 将RGB转为索引图，索引颜色数规定为n个，可以指定转换时运用什么方式进行颜色处理

有二种方式，一种为dither 抖动，另一种为nodither 不抖动。

抖动是什么意思呢，就是为了让图片看起来更加自然，每个像素的点会中和其邻域的点的颜色，这样就会显得这个点的像素颜色不那么突兀。

f = imread('onion.png');
subplot(2, 3, 1), imshow(f), title('原图');
subplot(2, 3, 4), bar([imhist(f(:, :, 1), 10), imhist(f(:, :, 2), 10), imhist(f(:, :, 3), 10)]), axis tight, title('原图的三色直方图');
[X, map] = rgb2ind(f, 256); % 默认为dither
subplot(2, 3, 2), imshow(ind2rgb(X, map)), title('dither');
f = ind2rgb(X, map);
subplot(2, 3, 5), bar([imhist(f(:, :, 1), 10), imhist(f(:, :, 2), 10), imhist(f(:, :, 3), 10)]), axis tight, title('dither的三色直方图');
[X, map] = rgb2ind(f, 256, 'nodither');
subplot(2, 3, 3), imshow(ind2rgb(X, map)), title('nodither');
f = ind2rgb(X, map);
subplot(2, 3, 6), bar([imhist(f(:, :, 1), 10), imhist(f(:, :, 2), 10), imhist(f(:, :, 3), 10)]), axis tight, title('nodither的三色直方图');

输入:



输出:

dither(image) 该方法即为上面rgb2ind函数的dither方法，该函数处理灰度图效果最明显

f = imread('coins.png');
g = dither(f);
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(g), title('抖动出来的二值图');

输入:



输出:

彩色空间转换

我们之前都是对RGB图像对彩色图像进行操作（间接或直接），很少用到其他的彩色空间，这章我们学习从RGB变换得到的彩色空间（彩色模型）。

这些彩色空间包括NTSC、YCbCr、HSV、CMY、CMYK、HSI。

NTSC彩色空间

NTSC彩色模型为模拟电视， 分为三个分量，YIQ ，分为别亮度，色调， 对比度（就是我们常常调自己家电视的亮度、色调、对比度），这个模型使得同一个信号，可以在彩色电视显示，也可以在黑白电视显示。

rgb2ntsc(image) 将RGB图像变换为NTSC

f = imread('onion.png');
g = rgb2ntsc(f);
subplot(1, 2, 1), imshow(f), title('RGB');
subplot(1, 2, 2), imshow(g), title('NTSC');

输入:



输出:

ntsc2rgb(image) 将NTSC图像变换为RGB图像

f = imread('onion.png');
g = rgb2ntsc(f);
img = ntsc2rgb(g);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('NTSC');
subplot(1, 3, 3), imshow(img), title('RGB');

输入:



输出:

myRgb2ntsc(image) 自制的rgb2ntsc，将RGB变换NTSC格式

f = imread('onion.png');
g = rgb2ntsc(f);
img = myRgb2ntsc(f);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('NTSC rgb2ntsc');
subplot(1, 3, 3), imshow(img), title('NTSC myRgb2ntsc');

输入:



输出:

myNtsc2rgb(image) 自制的ntsc2rgb, 将NTSC变换为RGB

f = imread('onion.png');
g = rgb2ntsc(f);
img = myNtsc2rgb(g);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('NTSC');
subplot(1, 3, 3), imshow(img), title('RGB');

输入:



输出:

YCbCr彩色空间

YCbCr彩色模型用于数字视频中。模型分为3个分量，Y（亮度） Cb（蓝色分量和参数值的差） Cr（红色分量和参考值的差）

rgb2ycbcr(image) 将RGB图像转为YCbCr

f = imread('onion.png');
g = rgb2ycbcr(f);
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(g), title('YCbCr');

输入:



输出:

ycbcr2rgb(image) 将YCbCr转为RGB

f = imread('onion.png');
g = rgb2ycbcr(f);
img = ycbcr2rgb(g);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('YCbCr');
subplot(1, 3, 3), imshow(img), title('RGB');

输入:



输出:

HSV彩色空间

HSV分别代表色调、饱和度、数值，该模型更加贴近人们的经验和描述彩色感觉时所采用的方式，在艺术领域，称为色泽、明暗、调色。

rgb2hsv(image) 将RGB转为HSV

f = imread('onion.png');
g = rgb2hsv(f);
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(g), title('HSV');

输入:



输出:

**hsv2rgb(image) 将HSV转为RGB

f = imread('onion.png');
g = rgb2hsv(f);
img = hsv2rgb(g);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('hsv');
subplot(1, 3, 3), imshow(img), title('RGB');

输入:



输出:

CMY和CMYK彩色空间

该模型常常用于彩色打印。

CMY模型分为三个分量，C 青色（Cyan）、M 深红色（Magenta）、Y 黄色（Yellow）。

CMYK则是在CMY上加一个黑色

imcomplement(image) 获取图片负片, 可以实现CMY模型与RGB互换

f = imread('onion.png');
g = imcomplement(f);
img = imcomplement(g);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('CMY');
subplot(1, 3, 3), imshow(img), title('RGB');

输入:



输出:

HSI彩色空间

HSI模型三个分量为， H 色彩（hue）， S 饱和度（saturation）， I 强度（intensity）

rgb2hsi(image) 将rgb图像转为hsi图像

f = imread('onion.png');
g = rgb2hsi(f);
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(g), title('HSI');

输入:



输出:

hsi2rgb(image) 将hsi图像变换为rgb图像

f = imread('onion.png');
g = rgb2hsi(f);
img = hsi2rgb(g);
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(g), title('HSI');
subplot(1, 3, 3), imshow(img), title('RGB');

输入:



输出:

附录

小实验

HSV与HSI图像显示比较

f = imread('onion.png');
subplot(1, 3, 1), imshow(f), title('原图');
subplot(1, 3, 2), imshow(rgb2hsi(f)), title('HSI');
subplot(1, 3, 3), imshow(rgb2hsv(f)), title('HSV');

输入:



输出:



可以看出HSI与HSV的图像有一些差别，HSI偏深，HSV偏浅，再看一下这二个模型的三个分量

f = imread('onion.png');
subplot(2, 3, 1), imshow(f), title('原图');
subplot(2, 3, 2), imshow(rgb2hsi(f)), title('HSI');
subplot(2, 3, 3), imshow(rgb2hsv(f)), title('HSV');
subplot(2, 3, 4), sBar(f, 10), title('原图'), xlabel('R红 G绿 B蓝');
subplot(2, 3, 5), sBar(rgb2hsi(f), 10), title('HSI'), xlabel('H色彩 S饱和度 I强度');
subplot(2, 3, 6), sBar(rgb2hsv(f), 10), title('HSV'), xlabel('H色调 S饱和度 I数值');

输入:



输出:



主要看第HSI和HSV图，可以看出，色调（H）基本相同，饱和度（S）是有差距的，强度与数值的差距巨大。所以HSI与HSV是不同的，如果只是想要调整色调，则无所谓，但是在调节数值和强度的时候，就需要好好琢磨了。

myHisteq(image) 对彩色图像进行直方图均衡化

f = imread('onion.png');
g = myHisteq(f, 'hsv');
p = myHisteq(f, 'rgb');
subplot(2, 3, 1), imshow(f), title('原图');
subplot(2, 3, 2), imshow(g), title('仅对HSV的V通道进行均衡化后');
subplot(2, 3, 3), imshow(p), title('对RGB通道进行均衡化后');
subplot(2, 3, 4), sBar(f, 25), title('原图');
subplot(2, 3, 5), sBar(g, 25), title('仅对HSV的V通道进行均衡化后');
subplot(2, 3, 6), sBar(p, 25), title('对RGB通道进行均衡化后');

输入:



输出:

clahe(image) 对比度受限的自适应直方图均衡化

本小实验主要看一下自制的clahe、工具箱中的adapthisteq、和利用imadjust进行灰度变换去雾的差别

f = imread('大雾.jpg');
g1 = clahe(f);
g2 = quwu(f);
g3 = cat(3, adapthisteq(f(:, :, 1)), adapthisteq(f(:, :, 2)), adapthisteq(f(:, :, 3)));
subplot(2, 4, 1), imshow(f), title('原图');
subplot(2, 4, 2), imshow(g1), title('CLAHE');
subplot(2, 4, 3), imshow(g2), title('imadjust');
subplot(2, 4, 4), imshow(g3), title('adapthisteq');
subplot(2, 4, 5), sBar(f, 10), title('原图');
subplot(2, 4, 6), sBar(g1, 10), title('CLAHE');
subplot(2, 4, 7), sBar(g2, 10), title('imadjust');
subplot(2, 4, 8), sBar(g3, 10), title('adapthisteq');

输入:



输出:



可以看出去雾效果还是很明显的，imadjust的话，还是有很多雾点没有去除，adapthisteq和clahe的话相对好一些。

代码

myRgb2ntsc.m

function img = myRgb2ntsc(f)
f = im2double(f);
d =  [0.229, 0.587, 0.114; 0.596, -0.274, -0.332; 0.211, -0.523, 0.312] ;
r = f(:, :, 1);
g = f(:, :, 2);
b = f(:, :, 3);
y = r * d(1, 1) + g * d(1, 2) + b * d(1, 3);
i = r * d(2, 1) + g * d(2, 2) + b * d(2, 3);
q = r * d(3, 1) + g * d(3, 2) + b * d(3, 3);
img = cat(3, y, i, q);
end

myNtsc2rgb.m

function img = myNtsc2rgb(f)
f = im2double(f);
d =  [1.0, 0.956, 0.621; 1.0, -0.272, -0.647; 1.0, -1.106, 1.703];
y = f(:, :, 1);
i = f(:, :, 2);
q = f(:, :, 3);
r = y * d(1, 1) + i * d(1, 2) + q * d(1, 3);
g = y * d(2, 1) + i * d(2, 2) + q * d(2, 3);
b = y * d(3, 1) + i * d(3, 2) + q * d(3, 3);
img = cat(3, r, g, b);
end

rgb2hsi.m

function hsi = rgb2hsi(image)
rgb = im2double(image);
r = rgb(:, :, 1);
g = rgb(:, :, 2);
b = rgb(:, :, 3);

num = 0.5 * ((r - g) + (r - b));
den = sqrt((r - g) .^ 2 + (r - b) .* (g - b));
theta = acos(num ./ (den + eps));

H = theta;
H(b > g) = 2 * pi - H(b > g);
H = H / (2 * pi);

num = min(min(r, g), b);
den = r + g + b;
den(den == 0) = eps;
S = 1 - 3 .* num ./ den;
H(S == 0) = 0;
I = (r + g + b) / 3;
hsi = cat(3, H, S, I);
end

hsi2rgb.m

function rgb = hsi2rgb(hsi)
H = hsi(:, :, 1) * 2 * pi;
S = hsi(:, :, 2);
I = hsi(:, :, 3);
width = size(hsi, 1);
height =  size(hsi, 2);
R = zeros(width, height);
G = zeros(width, height);
B = zeros(width, height);

idx = find(0 <= H < 2 * pi / 3);
B(idx) = I(idx) .* (1 - S(idx));
R(idx) = I(idx) .* (1 + S(idx) .* cos(H(idx)) ./ cos(pi / 3 - H(idx)));
G(idx) = 3 * I(idx) - (R(idx) + B(idx));

idx = find((2 * pi / 3 <= H & H < 4 * pi / 3));
R(idx) = I(idx) .* (1 - S(idx));
G(idx) = I(idx) .* (1 + S(idx) .* cos(H(idx) - 2 * pi / 3) ./ cos(pi - H(idx)));
B(idx) = 3 * I(idx) - (R(idx) + G(idx));

idx = find((4 * pi / 3 <= H) & (H <= 2 * pi));
G(idx) = I(idx) .* (1 - S(idx));
B(idx) = I(idx) .* (1 + S(idx) .* cos(H(idx) - 4 * pi / 3) ./ cos(5 * pi / 3 - H(idx)));
R(idx) = 3 * I(idx) - (G(idx) + B(idx));

rgb = cat(3, R, G, B);
rgb = max(min(rgb, 1), 0);
end

sBar.m

function sBar(data, varargin)
if ndims(data) == 3
if nargin() == 1
b = bar([imhist(data(:, :, 1)) imhist(data(:, :, 2)) imhist(data(:, :, 3))]);
elseif nargin() == 2
b = bar([imhist(data(:, :, 1), varargin{1}) imhist(data(:, :, 2), varargin{1}) imhist(data(:, :, 3), varargin{1})]);
end
elseif ndims(data) == 2
b = bar(imhist(data));
elseif ndims(data) == 1
b = bar(data)
else
error '错误的输入';
end
axis tight;
color=['r', 'g', 'b'];
for i=1:3
set(b(i),'FaceColor',color(i));
end
end

myHisteq.m

function img = myHisteq(f, varargin)
% 对彩色图像进行均衡化
% 可以传入的图片为灰度图，RGB，索引图
% g = myHisteq(img);
% g = myHisteq(X, map);

if numel(varargin) == 1
type = getImageType(f, varargin);
else
type = getImageType(f);
end
if strcmp(type,'RGB图') == 1
p = rgb2hsv(f); % 对HSV 的 V 通道进行处理（原图不会失真）
img = cat(3, p(:,:,1) .* 256, p(:,:,2) .* 256, myHisteq(uint8(p(:,:,3) .* 256)));
img = hsv2rgb(tofloat(img));
elseif strcmp(type,'索引图') == 1
img = myHisteq(ind2rgb(f, varargin));
else
[height, width] = size(f);
numPixel = zeros(1, 256);
for i=1:height
for j=1:width
numPixel(f(i, j) + 1) = numPixel(f(i, j) + 1) + 1;
end
end
probPixel = numPixel ./ (height * width);
cumuPixel = cumsum(probPixel);
cumuPixel = uint8(256 .* (cumuPixel + 0.1));
img = cumuPixel(f(:,:) + 1);
end
end

quwu.m

function quwu(f, varargin)
% QUWU
% 使用imadjust 和 stretchlim 进行简单的去雾
% quwu(f)
% quwu(f, 0.1)
% quwu(f, 0.95);
% 第二，第三个参数分别为stretchlim的TOL的Low High

r = f(:,:,1);
g = f(:,:,2);
b = f(:,:,3);
if nargin == 1
low = 0;
high = 1;
elseif nargin == 2
if isfloat(varargin{1}) && varargin{1} >= 0 && varargin{1} <= 1
low = varargin{1};
else
error '第二个参数必须为小数([0, 1])';
end
high = 1;
elseif nargin == 3
if isfloat(varargin{1}) && varargin{1} >= 0 && varargin{1} <= 1
low = varargin{1};
else
error '第二个参数必须为小数([0, 1])';
end
if isfloat(varargin{2}) && varargin{2} >= 0 && varargin{2} <= 1
high = varargin{2};
else
error '第三个参数必须为小数([0, 1])';
end
else
error '传入参数个数错误';
end
Low_High = stretchlim(r, [low high]);
rr = imadjust(r, Low_High, []);
Low_High = stretchlim(g, [low high]);
gg = imadjust(g, Low_High, []);
Low_High = stretchlim(b, [low high]);
bb = imadjust(b, Low_High, []);
l = cat(3, rr, gg, bb);
imshow(l);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(l);
end

CLAHE

以下函数除了clahe, 其他都为网上的资源

[b]clipHistogram.m[/b]

function [Hist] = clipHistogram(Hist,NrBins,ClipLimit,NrX,NrY)
%  This function performs clipping of the histogram and redistribution of bins.
%  The histogram is clipped and the number of excess pixels is counted. Afterwards
%  the excess pixels are equally redistributed across the whole histogram (providing
%  the bin count is smaller than the cliplimit).

for i = 1:NrX
for j = 1:NrY
%   Calculate the total number of excess pixels.
NrExcess = 0;
for nr = 1:NrBins
excess=Hist(i,j,nr) - ClipLimit;
if excess > 0
NrExcess = NrExcess + excess;
end
end

%  Clip histogram and redistribute excess pixels in each bin
binIncr = NrExcess / NrBins;
upper = ClipLimit - binIncr;
for nr = 1:NrBins
if Hist(i,j,nr) > ClipLimit
Hist(i,j,nr) = ClipLimit;
else
if Hist(i,j,nr) > upper
NrExcess = NrExcess + upper - Hist(i,j,nr);
Hist(i,j,nr) = ClipLimit;
else
NrExcess = NrExcess - binIncr;
Hist(i,j,nr) = Hist(i,j,nr) + binIncr;
end
end
end

if NrExcess > 0
stepSize = max(1,fix(1+NrExcess/NrBins));
for nr = 1:NrBins
NrExcess = NrExcess - stepSize;
Hist(i,j,nr) = Hist(i,j,nr) + stepSize;
if NrExcess < 1
break;
end
end
end

end
end

[b]interpolate.m[/b]

function [subImage] = interpolate(subBin,LU,RU,LB,RB,XSize,YSize)
%  pImage      - pointer to input/output image
%  uiXRes      - resolution of image in x-direction
%  pulMap*     - mappings of greylevels from histograms
%  uiXSize     - uiXSize of image submatrix
%  uiYSize     - uiYSize of image submatrix
%  pLUT        - lookup table containing mapping greyvalues to bins
%  This function calculates the new greylevel assignments of pixels within a submatrix
%  of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation
%  between four different mappings in order to eliminate boundary artifacts.
%  It uses a division; since division is often an expensive operation, I added code to
%  perform a logical shift instead when feasible.
%

subImage = zeros(size(subBin));
num = XSize * YSize;
for i = 0:XSize-1
inverseI = XSize - i;
for j = 0:YSize-1
inverseJ = YSize - j;
val = subBin(i+1,j+1);
subImage(i+1,j+1) = fix((inverseI*(inverseJ*LU(val) + j*RU(val)) ...
+ i*(inverseJ*LB(val) + j*RB(val)))/num);
end
end

[b]mapHistogram.m[/b]

function [Map] = mapHistogram(Hist,Min,Max,NrBins,NrPixels,NrX,NrY)
%  This function calculates the equalized lookup table (mapping) by
%  cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max].

Map=zeros(NrX,NrY,NrBins);

Scale = (Max - Min)/NrPixels;

for i = 1:NrX
for j = 1:NrY

Sum = 0;
for nr = 1:NrBins
Sum = Sum + Hist(i,j,nr);
Map(i,j,nr) = fix(min(Min + Sum*Scale,Max));
end

end
end

[b]makeLUT.m[/b]

function [LUT] = makeLUT(Min,Max,NrBins)
%  To speed up histogram clipping, the input image [Min,Max] is scaled down to
%  [0,uiNrBins-1]. This function calculates the LUT.

Max1 = Max + max(1,Min) - Min;
Min1 = max(1,Min);

BinSize = fix(1 + (Max - Min)/NrBins);
LUT = zeros(fix(Max - Min),1);
for i=Min1:Max1
LUT(i) = fix((i - Min1)/BinSize);
end

[b]makeHistogram.m[/b]

%  This function classifies the greylevels present in the array image into
%  a greylevel histogram. The pLookupTable specifies the relationship
%  between the greyvalue of the pixel (typically between 0 and 4095) and
%  the corresponding bin in the histogram (usually containing only 128 bins).

Hist=zeros(NrX,NrY,NrBins);

for i=1:NrX
for j=1:NrY
bin=Bin(1+(i-1)*XSize:i*XSize,1+(j-1)*YSize:j*YSize);
for i1=1:XSize
for j1=1:YSize
Hist(i,j,bin(i1,j1)) = Hist(i,j,bin(i1,j1)) + 1;
end
end
end
end

[b]clahe.m[/b]

function output = clahe(Img, varargin)
ClipLimit = 1.75;
if numel(varargin) == 1
ClipLimit = varargin{1};
end
type = getImageType(Img);
if strcmp(type,'RGB图') == 1
output = cat(3, clahe(Img(:,:,1)), clahe(Img(:,:,2)), clahe(Img(:,:,3)));
else
[h,w] = size(Img);
output = zeros(h, w);
output = padarray(output, [5, 5], 'both');
minV = double(min(min(Img)));
maxV = double(max(max(Img)));
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 = 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));
[Hist(i, j, :), x] = imhist(tmp, 256);
end
end

Hist = circshift(Hist,[0, 0, -1]);
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 = uint8(output(1:h, 1:w));
end
end