Re: Re[2]: [Dev-C++] resampling

[Mac R. <re...@no...>]

Zodiaq,

Hope you don't mind my sending this to the group, but thought I should
spread it around in case anyone else needed this kind of resampling.

Find a book on how to draw fast lines - The Zen of Graphics Programming, or
almost any other book that discusses the Bresenham line algorithm.  Once
you've got a good line algorithm (Bresenham, or some better algorithm might
work as well), I *think* you could modify it to do the resampling you want.

In most fast line drawing systems, you first decide which of the two
directions (vertical/horizontal) is the "long" way.  You then step along
that "long" axis, incrementing another variable, and whenever that variable
passes a certain threshold, you jump one pixel along the "short" axis.
Basically, you are using the slope and saying "for every pixel I move to
the right (to be concrete for a moment), I need to move 1/3 pixel up".
Rather than working with fractions, it is more common in this case to
simply add 1 to the accumulator at each step, and step up when the
accumulator exceeds 3.  The basic point is that you can scale the issue
such that you do incremental work with integer values, which makes for a
really fast algorithm.

Now, how do we adjust this to your problem?

Instead of the conventional X and Y axes, one of your axes is the new bin
number, and the other axis is in some way related to the energy units in
the old bins (I am thinking out loud here, so we may wander around the
problem a bit before we get to the answer...).  To go with the concrete
example from above, where the line is three times longer horizontally than
vertically, we would replace the vertical Y axis with the new bin number.
You would start in the first bin, and as you distribute energy units you
would add an incremental amount to an accumulator.  When the accumulator
crossed a certain threshold, you would move "up" one unit along the "Y
axis" by incrementing to the next new bin.  Energy units would then be
transferred to this bin.

So now we either have to work out how much to add to the accumulator, or we
have to figure out what the threshold is...  Hmmm...  It seems like we
ought to be able to do something with those integral "number of bins"
values (20 and 13, from your example).  Or maybe we want to look directly
at the energy units that are stored in the bins -- as a histogram, I am
assuming that the number stored in any given bin is an integral value,
namely the number of pixels whose values fell within that bin.

I'm going to stop for a second, because that last statement brings up a
concern.  "Resampling" a histogram is not valid...  Let's keep it simple:
Say I start with a 4 bin histogram.  Now I feed it a picture (10,000 total
pixels - 100x100) where ALL of the pixels are 74% intensity.  That means
that the histogram looks like { 0, 0, 10000, 0 }, because all of the pixels
were in the third bin.  Now if I resample to a 3 bin histogram, those
10,000 pixels are going to get split between the second bin and the third
bin.  Each new bin takes 4/3 (1.33333333333) bins of the old bins' values.
So the new histogram looks like { 0, 6667, 3333 }.  (I had to draw a little
bar chart on a piece of paper and then overlay a 3 way split over it...)
But if we think about the original picture, this is blatantly incorrect --
the 3 bin histogram *should* look like { 0, 0, 10000 }, because all of the
pixels were at 74% intensity, and the threshold between the second and
third bins is at 66.666666...%.

I think that if you want the answers to retain meaning, you are going to
have to re-histogram with the new number of bins.

One method that *might* help is to perform one histogram with the maximum
valid number of bins, and then resample down from there.  For instance, if
you are looking at 8bit grayscale images, the "maximum valid number of
bins" is 256.  You could perform one histogram with 256 bins, and then
resample down from there.  This would avoid the discretization problems
just mentioned.  The 4:3 resampling I did above only introduced errors
because the values had been discretized into 4 values, even though they
originally represented many more (like 256).  If the input histogram fully
represents the range of the input pixels, you no longer have a problem.
You do have to work with a bigger histogram, but it would be correct, and
it would be much faster than re-histogramming the full image -- a maximum
of 256 lookups, compared to 640x480=307,200 lookups for a moderate sized
image.

OK, continuing from where I left off, now that we have a way to keep the
meaning of the numbers -- say you wanted to resample the 256 entry
histogram to a 20 entry histogram.  (I will work in fractions for the
initial layout)  As you stepped past each bin<256>, you would add the value
in that bin<256> into the current bin<20> and add 20.0/256.0 to the
accumulator.  Then, when the accumulator exceeded 1.0, you would step to
the next bin<20> and subtract 1.0 from the accumulator.  Keep that up until
you finish all 256 bins.

But that required 256 floating point additions and compares, and 20
floating point subtracts.  Can we do better?  Well, if we multiply the step
(20.0/256.0) and the threshold (1.0) both by 256, then the step becomes 20
and the threshold becomes 256.  Now for each bin<256> you can add 20 to the
accumulator, and when the accumulator exceeds 256, you move to the next
bin<20> and subtract 256 from the accumulator.  Note that the step is just
the output number of bins, and the threshold is just the input number of
bins...

Or, in code form: (this is short, not error checked, and probably should be
rewritten in terms of STL containers, but we're going for maximal
understanding and shortest typing time for me.) (oh, and I tend to be wordy
in my code, too)

void Resample(long InputHist[], long OutputHist[], 
              long InputSize, long OutputSize)
{
	long Accumulator = 0;
	long OutputBin   = 0;
	long InputBin    = 0;

	// Clear out the output histogram.
	memset(&(OutputHist[0]), 0, sizeof(long)*OutputSize);

	for (InputBin = 0; InputBin < InputSize; ++InputBin)
	{
		OutputHist[OutputBin] += InputHist[InputBin];
		Accumulator += OutputSize;
		if (Accumulator > InputSize)
		{
			++OutputBin;
			Accumulator -= InputSize;
		}
	}
}

Note that, at least in theory, you could apply this to the 20:13
resampling, but I don't think the result would mean exactly what you
wanted.  I would stick with the 256:20 and then 256:13 resampling...

There is still an error, but there is no "correct" solution -- a 20 bin
histogram places bin boundaries at fractional intensity values.  Since your
original input was already discretized in the process of capture, you have
no way of knowing which side of the fractional intensity boundary the pixel
belongs on.  I believe the code above will always place pixels below
fractional boundaries (all 173 brightness pixels will be placed in the bin
below the 173.259382... boundary).  I have not spent much time thinking
about this to verify it, and I don't see any way to fix it that is any less
arbitrary than the above.  Also, any attempt to "fix" it will add
computational complexity, and if it doesn't increase accuracy, why do it?

I should point out that I have not compiled or tested this in ANY way
whatsoever.

I should also explain why I used "++InputBin" and "++OutputBin" instead of
the more common "InputBin++" and "OutputBin++".  A post increment (var++)
forces the compiler to store the value of the variable, increment the
variable, and then return the stored value (the value before the
increment).  A pre increment (++var) allows the compiler to increment the
variable and then return the variable, with no additional storage.  In both
of the above cases, it *should* be irrelevant whether the addition occurred
before or after the value was considered (the increments are in standalone
statements).  I would expect any decent optimizing compiler to laugh and
just do the increment without worrying about storage and return values.
But if this really is speed critical code, it never hurts to think about
such things...

Oh, and if you are going to run this on a x86 CPU (386, 486, Pentium, K6,
Athlon, etc), avoid "short integer" like the plague!  x86 architectures can
manipulate chars (8 bits) and longs (32 bits) normally, and at full speed.
But shorts (16 bits) require adding a special addressing byte onto the
opcode, and that extra addressing mode can slow things down by a factor of
50%-100%.  I ran into that on a video processing app I was writing, where I
used loop indices that were short integers.  It got tremendously faster
when I switched to longs, or just ints, but I try to always specify what I
mean.  I suppose in this case, "int" would have specified what I meant --
"int" generally means "the fastest integral type for the given hardware".

Just as a side note - I'm kinda of a video nut.  What are you doing with
pictures that you are equalizing 200,000 times?  Is it a project that I
could go poke around on a web site and learn more about?

Good luck, and I hope that the above works out for you.
Mac

At 06:38 PM 7/14/02 +0200, Zodiaq wrote:
>Hello Mac and all Dev-C++ users,
>
>Sunday, July 14, 2002, 5:48:45 PM, Mac wrote as follows:
>
>MR> Will the new array always be larger, always be smaller, or be larger or
>MR> smaller at different times during the run?
>
>Smaller and in special cases the same size... but with last option of
>course no problem...
>
>MR> If the new array is always smaller, a relatively simple folding can be
>MR> done.  Depending on the meaning of your arrays, you may need to split
>MR> entries and fold part of the value down and part of the value up.
>
>Yes I though about that... I think also that best way is to calculate
>area under the graph (between two points) , and then to get final
>result and devide it by number of area in x taken... heh, for example:
>
>first array 20 points...
>second that I want 13 points...
>20 / 13 = 1.538...
>so one point in second arrays need 1.538 of point in first... so to
>get it I need to calculate area between point 1 and 2 and add area
>between 2 and 2.538... but this makes a lot of calculations so I
>wonder is there fastest way, and more important stable...
>
>MR> What do the values in these arrays represent?
>
>These values represent points of histogram made from pictures... I
>need to scale these histograms and then make some more calculations,
>and then save data to make even more calculations later... so values
>in first array are in range <0, ~500>
>
>And one more thing - this function is time critical to my program, it
>will be executed 200 000 times... so I'm looking for best way...
>
>MR> Mac
>MR> -
>MR>  _____________________________     /"\
>MR>  Mac Reiter                        \ /    ASCII Ribbon Campaign
>MR>  Nomadics, Inc.                     X     Against HTML Mail
>MR>  re...@no...               / \   (To join the campaign, simply
use
>MR>                                           this in your signature.)
>
>
>
>-- 
>regards, Zodiaq
>
>
-
 _____________________________     /"\
 Mac Reiter                        \ /    ASCII Ribbon Campaign
 Nomadics, Inc.                     X     Against HTML Mail
 re...@no...               / \   (To join the campaign, simply use
                                          this in your signature.)