Some simplifications and some final insights on the stretching algorithm.
We started with;
Code:
void stretch_line_bhm(img *imgSrc, img *imgDst, int widthSrc, int widthDst)
{
int x = 0;
int y = 0;
int dx = widthDst;
int dy = widthSrc;
int e = 0;
// determine iteration direction
if(dy < dx)
{
e = 2*dy - dx;
// iterate over x (scale up, widthSrc < widthDst)
while(x < dx)
{
// setpixel(x,y)
imgDst[x] = imgSrc[y];
if(e > 0)
{
y++;
e -= 2*dx;
}
x++;
e += 2*dy;
}
}
else
{
e = dy - 2*dx;
// iterate over y (scale down, widthSrc >= widthDst)
while(y < dy)
{
// setpixel(x,y)
imgDst[x] = imgSrc[y];
if(e < 0)
{
x++;
e += 2*dy;
}
y++;
e -= 2*dx;
}
}
}
We can modify the algorithm such that we don't need to run through the else-
part (slope dy/dx >= 1). Well, if dy/dx is greater than 1, we will actually
scale the source array down. Hence, some pixels need to be skipped. So let us
see what happens if we leave the else-part out;
Code:
// does not work
void stretch_line_bhm(img *imgSrc, img *imgDst, int widthSrc, int widthDst)
{
int x = 0;
int y = 0;
int dx = widthDst;
int dy = widthSrc;
int e = 2*dy - dx;
// iterate over x
while(x < dx)
{
// setpixel(x,y)
imgDst[x] = imgSrc[y];
if(e > 0)
{
y++;
e -= 2*dx;
}
x++;
e += 2*dy;
}
}
With dy > dx the error term is > 0 and remains > 0 after the if-statement.
Hence, in each iteration y will be incremented by one which is not what we
want, because incrementing y by one in each iteration would imply no down
scaling at all. But since dy > dx, we need to skip some of those y's somehow.
In the standard Bresenham integer line algorithm the error term e in the
if(dy < dx) { ... if(e > 0){... e -= 2*dx; } ...} section will be
reinitialized when we move one step in y by subtracting 2*dx from e making it
<= 0. This is based on the assumption that dy < dx. If dy is greater than dx,
we would need to subtract 2*dx multiple times from e, because for each step in
x we would make more than one step in y. Doing so, and with respect to drawing
a line, would introduce gaps along the line, because we would jump more than
one pixel in y per step in x. However, since we don't want to draw the line
at all, and because the y in our algorithm indexes our source array, we will
arrive at the correct y value when subtracting 2*dx from e multiple times
until e <= 0 -- while incrementing y in the process skipping those y's we
can't account for due to down scaling. This effect will become much more clear
if we study some further simplifications. For the time being we will replace
if(e > 0) with while(e > 0) and the algorithm will work as expected;
Code:
void stretch_line_bhm(img *imgSrc, img *imgDst, int widthSrc, int widthDst)
{
int x = 0;
int y = 0;
int dx = widthDst;
int dy = widthSrc;
int e = 2*dy - dx;
// iterate over x
while(x < dx)
{
imgDst[x] = imgSrc[y];
while(e > 0)
{
y++;
e -= 2*dx;
}
x++;
e += 2*dy;
}
}
Now let us remove Bresenham's error-sauce. The error term can be derived from
the condition; dy/dx <= 1/2 (midpoint rule) <=> e := 2*dy - dx <= 0. Well,
lets modify this condition somewhat, i.e. to dy/dx <= 1 to show something
much more interesting, and let us define a similar error term (named e again)
by e := dy <= dx. For each step in x we keep track of the differential
increment in y, i.e. dy, by adding dy to e in each iteration on x. Now it so
happens that after some iterations (depending on dy/dx) the error term e will
be greater than dx which means; we now have to increment y to follow the line,
or to get a new index into our source array, while subtracting dx from e to
reinitialize the error term again. Ok, lets do this;
Code:
void stretch_line_bhm(img *imgSrc, img *imgDst, int widthSrc, int widthDst)
{
int x = 0;
int y = 0;
int dx = widthDst;
int dy = widthSrc;
int e = dy;
// iterate over x
while(x < dx)
{
imgDst[x] = imgSrc[y];
while(e > dx)
{
y++;
e -= dx;
}
x++;
e += dy;
}
}
Fine.
Now lets apply some magic. Lets relabel some of the variables and let us see
if we won't suddenly get enlighted;
Code:
void stretch_line_bhm(img *imgSrc, img *imgDst, int widthSrc, int widthDst)
{
int x = 0;
int y = 0;
int numerator = widthSrc; // dx
int denominator = widthDst; // dy
int reminder = numerator; // e
// iterate over x
while(x < widthDst)
{
imgDst[x] = imgSrc[y];
while(reminder > denominator)
{
y++;
reminder = reminder - denominator;
}
x++;
reminder = reminder + numerator;
}
}
Now what does it look like? It's division by repeated subtraction! We want to
compute dy/dx = widthSrc/widthDst = numerator/denominator, yet we won't use a
floating-point divide operation, we instead emulate it by doing an integer
divide using repeated subtraction, and we also compute the reminder to keep
track of the fractional part. We initialize the reminder with the numerator
and subtract the denominator (divisor) from it as along as the reminder
remains larger than the denominator. Hence, if (for example) widthSrc >
widthDst (scaling down), then the subtracting-loop skips whole pixels, which
is what we want. Now after subtracting multiples of the denominator out of the
numerator, we will be left with a reminder > 0 if the numerator and
denominator won't divide evenly. The reminder actually represent the
fractional part/pixel we cannot address because we can only work with integer
pixels. The actual fraction is computed as reminder/denominator.
Lets say we have widthSrc = 64 and widthDst = 40 (scale down, 64/40 = 1.6).
Then we will be left with a reminder of 24 (fraction = 24/40 = 0.6). Now,
ideally, we would want to step 1.6 pixels in the source array for each setp in
the destination array. But we need integers here. If we would just round down
to 1, we would just copy every pixel out of the source array (leading to no
scaling at all for this example) up to 1*40 = 40 pixels with 24 pixels of the
source array not considered. On the other hand, if we round-up to 2, we would
scale down by 2 (0.4 more than what was asked for). Additionally, we would
over run our source array (2*40 = 80) by 16 pixels. Those issues are resolved
by keeping the fractional part (the reminder) for the next iteration, such
that we increment (in this example) our source array index y by 1, and at
times by 2 if the reminder has filled up sufficiently. Looking at the code
again, we see; after one iteration of the outer loop, the reminder equals 24 +
64 = 80. Hence, in the next iteration y will be incremented two times, yet in
the next² iteration only once. And so on. You may also go through the
algorithm for the case of scaling up with widthSrc = 40 and widthDst = 64 for
example.
Finally, lets rephrase our algorithm somewhat;
Code:
void stretch_line_bhm(img *imgSrc, img *imgDst, int widthSrc, int widthDst)
{
int x = 0;
int y = 0;
int rem = widthSrc;
// iterate over x
while(x < widthDst)
{
imgDst[x] = imgSrc[y];
while(rem > widthDst)
{
y++;
rem -= widthDst;
}
x++;
rem += widthSrc;
}
}
That's it! Sweet!
We could also now start to modify the integer scaler to account for the
skipped pixels when scaling down such that every destination pixel will be
an average of widthSrc/widthDst (= 64/40 = 1.6) source pixels making the
scaler very smooth. But this is way off of the original intention I had for
writing this article which follows below. Yet I implemented it nevertheless.
Having said all that, we now come to a very interesting conclusion; The
strength of Bresenham's integer line algorithm is not so much the optimal
error term, or of being all integer, which is cool nevertheless, mind you, but
the real insight into his algorithm is that he implicitly converted the
computation of the line (its slope) into a division algorithm based on
repeated subtraction!
