title: Digital Image Processing Assignment IV tags:
We modify our gen_color function in order to generate a blank image with given scale.
// generate a blank color image
void gen_color(BITMAP *bmImg, BITMAP *bmColor, uint32_t biHeight, uint32_t biWidth) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
bmColor->bmHeader = bmImg->bmHeader;
bmColor->bmInfoSize = bmImg->bmInfoSize;
bmColor->bmInfo = (BITMAPINFO *) malloc(bmColor->bmInfoSize);
bmColor->bmInfo->bmiHeader = *bmiHeader;
bmColor->bmInfo->bmiHeader.biWidth = biWidth;
bmColor->bmInfo->bmiHeader.biHeight = biHeight;
init_bmp(bmColor);
}
Normally, we enumerate each pixel in the new image, map it back to the original image, and use the neighboring pixels to interpolate. However, for translating, mirroring and shearing, there is a one-to-one match between the pixels in two images, so interpolation is unnecessary.
void translate(BITMAP *bmImg, BITMAP *bmTranslate, uint32_t disX, uint32_t disY) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
uint32_t height = bmiHeader->biHeight + disX;
uint32_t width = bmiHeader->biWidth + disY;
gen_color(bmImg, bmTranslate, height, width);
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
int32_t x = h;
int32_t y = w - disY;
uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
uint32_t _pos = h * bmTranslate->bmBytesPerRow + w * bmTranslate->bmBytesPerPel;
if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
for (uint8_t k = 0; k < 3; ++k) {
bmTranslate->bmData[_pos + k] = bmImg->bmData[pos + k];
}
}
}
}
}
void mirror(BITMAP *bmImg, BITMAP *bmMirror) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
uint32_t height = bmiHeader->biHeight << 1;
uint32_t width = bmiHeader->biWidth;
gen_color(bmImg, bmMirror, height, width);
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
int32_t x = bmiHeader->biHeight - 1 - h;
int32_t y = w;
uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
uint32_t _pos = h * bmMirror->bmBytesPerRow + w * bmMirror->bmBytesPerPel;
if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
for (uint8_t k = 0; k < 3; ++k) {
bmMirror->bmData[_pos + k] = bmImg->bmData[pos + k];
}
}
}
}
}
void shear(BITMAP *bmImg, BITMAP *bmShear, uint32_t disX) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
uint32_t height = bmiHeader->biHeight + disX;
uint32_t width = bmiHeader->biWidth;
gen_color(bmImg, bmShear, height, width);
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
int32_t x = h + disX * w / (width - 1) - disX;
int32_t y = w;
uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
uint32_t _pos = h * bmShear->bmBytesPerRow + w * bmShear->bmBytesPerPel;
if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
for (uint8_t k = 0; k < 3; ++k) {
bmShear->bmData[_pos + k] = bmImg->bmData[pos + k];
}
}
}
}
}
There are several methods of interpolation, such as nearest-neighbor interpolation, bilinear interpolation, and bicubic interpolation. Among these, nearest-neighbor interpolation has the highest speed, bicubic interpolation has the highest quality.
The bicubic interpolation can be summarized as solving a system of linear equations with $16$ variables, i.e., find $a{ij}$, s.t. $p(x, y)=\sum{i=0}^3 \sum{j=0}^3 a{ij} x^i y^j$. An interpolator with similar properties can be obtained by applying a convolution with the following kernel in both dimensions,
$$ \begin{equation}
W(x)=\begin{cases}
(a+2)|x|^3-(a+3)|x|^2+1 & ,\ |x|\le 1 \\\\
a|x|^3-5a|x|^2+8a|x|-4a & ,\ 1<|x|\le 2 \\\\
0 & ,\ |x|>2
\end{cases}
\end{equation} $$
where $a$ is usually set to $-0.5$. At this time, the equation can be expressed in a more friendly manner,
$$ \begin{equation}
p(t)=\frac{1}{2}
\left[\begin{matrix}
1 & t & t^2 & t^3
\end{matrix}\right]
\left[\begin{matrix}
0 & 2 & 0 & 0 \\\\
-1 & 0 & 1 & 0 \\\\
2 & -5 & 4 & -1 \\\\
-1 & 3 & -3 & 1
\end{matrix}\right]
\left[\begin{matrix}
f_{-1} \\\\ f_0 \\\\ f_1 \\\\ f_2
\end{matrix}\right]
\end{equation} $$
for $t$ between $0$ and $1$ for one dimension. Note that for $1$-dimensional cubic convolution interpolation $4$ sample points are required. For each inquiry two samples are located on its left and two samples on the right. These points are indexed from $−1$ to $2$ in this text. The distance from the point indexed with $0$ to the inquiry point is denoted by $t$ here. For two dimensions first applied once in $y$ and again in $x$.
double interpolate(double t1, double f0, double f1, double f2, double f3) {
double t2 = t1 * t1, t3 = t2 * t1;
double ret = (-t1 + 2 * t2 - t3) * f0;
ret += (2 - 5 * t2 + 3 * t3) * f1;
ret += (t1 + 4 * t2 - 3 * t3) * f2;
ret += (-t2 + t3) * f3;
return ret / 2;
}
First we need to calculate the size of the new image. Then we use nearest-neighbor interpolation and bicubic interpolation separately.
// nearest-neighbor interpolation
void simple_rotate(BITMAP *bmImg, BITMAP *bmRotate, double theta) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
double minX = fmin(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
double maxX = fmax(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
double minY = fmin(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
double maxY = fmax(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
minX = fmin(minX, fmin(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
maxX = fmax(maxX, fmax(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
minY = fmin(minY, fmin(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
maxY = fmax(maxY, fmax(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
uint32_t height = ceil(maxX) - floor(minX);
uint32_t width = ceil(maxY) - floor(minY);
gen_color(bmImg, bmRotate, height, width);
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
int32_t x = round((h + minX) * cos(theta) + (w + minY) * sin(theta));
int32_t y = round((h + minX) * -sin(theta) + (w + minY) * cos(theta));
uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
uint32_t _pos = h * bmRotate->bmBytesPerRow + w * bmRotate->bmBytesPerPel;
if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
for (uint8_t k = 0; k < 3; ++k) {
bmRotate->bmData[_pos + k] = bmImg->bmData[pos + k];
}
}
}
}
}
// bicubic interpolation
void rotate(BITMAP *bmImg, BITMAP *bmRotate, double theta) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
double minX = fmin(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
double maxX = fmax(bmiHeader->biHeight * cos(theta), bmiHeader->biWidth * -sin(theta));
double minY = fmin(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
double maxY = fmax(bmiHeader->biHeight * sin(theta), bmiHeader->biWidth * cos(theta));
minX = fmin(minX, fmin(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
maxX = fmax(maxX, fmax(0, bmiHeader->biHeight * cos(theta) - bmiHeader->biWidth * sin(theta)));
minY = fmin(minY, fmin(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
maxY = fmax(maxY, fmax(0, bmiHeader->biHeight * sin(theta) + bmiHeader->biWidth * cos(theta)));
uint32_t height = ceil(maxX) - floor(minX);
uint32_t width = ceil(maxY) - floor(minY);
gen_color(bmImg, bmRotate, height, width);
for (uint8_t k = 0; k < 3; ++k) {
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
double x = (h + minX) * cos(theta) + (w + minY) * sin(theta);
double y = (h + minX) * -sin(theta) + (w + minY) * cos(theta);
uint32_t _pos = h * bmRotate->bmBytesPerRow + w * bmRotate->bmBytesPerPel;
if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
double g[4];
for (int8_t u = -1; u < 3; ++u) {
double f[4];
int32_t _x = max(0, min(bmiHeader->biHeight - 1, (int32_t) x + u));
for (int8_t v = -1; v < 3; ++v) {
int32_t _y = max(0, min(bmiHeader->biWidth - 1, (int32_t) y + v));
f[v + 1] = bmImg->bmData[_x * bmImg->bmBytesPerRow + _y * bmImg->bmBytesPerPel + k];
}
// interpolate horizontally
g[u + 1] = interpolate(y - (int32_t) y, f[0], f[1], f[2], f[3]);
}
// interpolate vertically
bmRotate->bmData[_pos + k] = adjust(interpolate(x - (int32_t) x, g[0], g[1], g[2], g[3]));
}
}
}
}
}
Similar to rotating, interpolation is vital for scaling. A slight difference is that, we can storage the results of all horizontal interpolations before performing the vertical, since they will be reused often.
// nearest-neighbor interpolation
void simple_scale(BITMAP *bmImg, BITMAP *bmScale, double ratioH, double ratioW) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
uint32_t height = ceil(ratioH * bmiHeader->biHeight);
uint32_t width = ceil(ratioW * bmiHeader->biWidth);
gen_color(bmImg, bmScale, height, width);
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
int32_t x = round(h / ratioH);
int32_t y = round(w / ratioW);
x = max(0, min(bmiHeader->biHeight - 1, x));
y = max(0, min(bmiHeader->biWidth - 1, y));
uint32_t pos = x * bmImg->bmBytesPerRow + y * bmImg->bmBytesPerPel;
uint32_t _pos = h * bmScale->bmBytesPerRow + w * bmScale->bmBytesPerPel;
if (x >= 0 && x < bmiHeader->biHeight && y >= 0 && y < bmiHeader->biWidth) {
for (uint8_t k = 0; k < 3; ++k) {
bmScale->bmData[_pos + k] = bmImg->bmData[pos + k];
}
}
}
}
}
// bicubic interpolation
void scale(BITMAP *bmImg, BITMAP *bmScale, double ratioH, double ratioW) {
BITMAPINFOHEADER *bmiHeader = &bmImg->bmInfo->bmiHeader;
uint32_t height = ceil(ratioH * bmiHeader->biHeight);
uint32_t width = ceil(ratioW * bmiHeader->biWidth);
gen_color(bmImg, bmScale, height, width);
double *p = (double *) malloc(bmiHeader->biHeight * width << 3);
for (uint8_t k = 0; k < 3; ++k) {
for (uint32_t h = 0; h < bmiHeader->biHeight; ++h) {
for (uint32_t w = 0; w < width; ++w) {
double y = w / ratioW;
double f[4];
for (int8_t u = -1; u < 3; ++u) {
int32_t _y = max(0, min(bmiHeader->biWidth - 1, (int32_t) y + u));
f[u + 1] = bmImg->bmData[h * bmImg->bmBytesPerRow + _y * bmImg->bmBytesPerPel + k];
}
// interpolate horizontally
p[h * width + w] = interpolate(y - (int32_t) y, f[0], f[1], f[2], f[3]);
}
}
for (uint32_t h = 0; h < height; ++h) {
for (uint32_t w = 0; w < width; ++w) {
double x = h / ratioH;
uint32_t _pos = h * bmScale->bmBytesPerRow + w * bmScale->bmBytesPerPel;
double f[4];
for (int8_t u = -1; u < 3; ++u) {
int32_t _x = max(0, min(bmiHeader->biHeight - 1, (int32_t) x + u));
f[u + 1] = p[_x * width + w];
}
// interpolate vertically
bmScale->bmData[_pos + k] = adjust(interpolate(x - (int32_t) x, f[0], f[1], f[2], f[3]));
}
}
}
free(p);
}