FFT convolution

[mkiefte]

Joel, would you like me to try to implement FFT convolution in the present zoom() or is that too ambitious?

I found some old, old Fortran code by Danielson and Lanczos that will do the trick. Easy enough to implement in Java.

• M.

[uckelman]

Thus spake “mkiefte”:

Joel, would you like me to try to implement FFT convolution in the present zo
om() or is that too ambitious?

I found some old, old Fortran code by Danielson and Lanczos that will do the
trick. Easy enough to implement in Java.

• M.

Sure, go ahead if you want.

–
J.

---

Messages mailing list
Messages@forums.vassalengine.org
forums.vassalengine.org/mailman/ … engine.org

Post generated using Mail2Forum (mail2forum.com)

[mkiefte]

The two listings on that page aren’t very good. The first is recursive which is not very efficient (e.g., if (n==1) return x;). The second is a verbatim copy of the Numerical Recipes FFT code which has very strict licensing issues that make it worse than useless. Because VASSAL is GPL, I won’t even look at Numerical Recipes lest I be tainted for life. The original FORTRAN code is relatively simple, unencumbered by licensing issues, and very, very fast.

If you’re interested in the long, horrible history of NR licensing, check out:
http://mingus.as.arizona.edu/~bjw/software/boycottnr.html
Real, bloody shame.

• M.

2008/5/9 Jeffrey Brent McBeth <mcbeth@broggs.org (mcbeth@broggs.org)>:

On Fri, May 09, 2008 at 08:27:57PM +0200, Joel Uckelman wrote:

Thus spake “mkiefte”:

Joel, would you like me to try to implement FFT convolution in the present zo
om() or is that too ambitious?

I found some old, old Fortran code by Danielson and Lanczos that will do the
trick. Easy enough to implement in Java.

• M.

Sure, go ahead if you want.

Might I suggest looking at something a little closer to home than Fortran?

http://www.ddj.com/cpp/199500857

–

“The man who does not read good books has no advantage over
the man who cannot read them.”
– Mark Twain

Post generated using Mail2Forum (mail2forum.com)

[uckelman]

Thus spake “Michael Kiefte”:

The two listings on that page aren’t very good. The first is recursive
which is not very efficient (e.g., if (n==1) return x;). The second is a
verbatim copy of the Numerical Recipes FFT code which has very strict
licensing issues that make it worse than useless. Because VASSAL is GPL, I
won’t even look at Numerical Recipes lest I be tainted for life. The
original FORTRAN code is relatively simple, unencumbered by licensing
issues, and very, very fast.

The one completely implemented as templates is a thing of beauty,
though.

J.

---

Messages mailing list
Messages@forums.vassalengine.org
forums.vassalengine.org/mailman/ … engine.org

Post generated using Mail2Forum (mail2forum.com)

[uckelman]

Thus spake Jeffrey Brent McBeth:

And his statement that recusion isn’t as bad as it used to be is spot on. I
only became aware of the NR licensing mess about a year ago. Prior to that,
my primary pet peeve was their coding style of random mixing of 0-based and
1-based arrays, nonsensical variable names, and complete lack of explanation
of what they were doing once they got into it.

1-based arrays? AUGH! MY EYES!!!

–
J.

---

Messages mailing list
Messages@forums.vassalengine.org
forums.vassalengine.org/mailman/ … engine.org

Post generated using Mail2Forum (mail2forum.com)

[mkiefte]

2008/5/9 Joel Uckelman <uckelman@nomic.net (uckelman@nomic.net)>:

Thus spake Jeffrey Brent McBeth:

And his statement that recusion isn’t as bad as it used to be is spot on.

The original algorithm takes advantage of some recurrences in the trigonometric functions as well which you don’t get in the recursive FFT function, so it’s not just the function calls. Either way, it works. And it was written at a time when people who wrote algorithms weren’t jerks.

I

only became aware of the NR licensing mess about a year ago. Prior to that,
my primary pet peeve was their coding style of random mixing of 0-based and
1-based arrays, nonsensical variable names, and complete lack of explanation
of what they were doing once they got into it.

1-based arrays? AUGH! MY EYES!!!

Oh yeah, I remember that. That was a pain in the ass. I recall they made some comment like all you have to do is subtract one from the array pointer to get one-based arrays and that they had never seen a compiler that objected to that. That was so strange. I used to use NR relatively liberally until someone pointed out the licensing restrictions. I’ve avoided it like the plague since.

• M.

Post generated using Mail2Forum (mail2forum.com)

[mkiefte]

I wrote the FFT algorithm in GeneralFilter.java only to discover that it’s mostly pointless. The resampled pixels almost never fall on integer indeces, so it would require upsampling followed by downsampling in most cases.

Now I realise why ImageMagick doesn’t use FFT convolution. However, it’s still good for mipmaps, because the resampling is always a factor of two. Here’s the code I wrote in the end based on the original FORTRAN. I’ll put it in the texture demo and see if it works.

Note that Math.sin() is only called in the outermost loop.

• M.

private static final void fft(final Complex[] x, boolean inverse) {
final int length = x.length;
final int nBits = Integer.SIZE - Integer.numberOfLeadingZeros(length);
if (1 << nBits-1 != length) {
throw new IllegalArgumentException(“Input data length must be a power of 2”);
}
else {
final int rshift = Integer.SIZE - nBits;
for (int i = 1; i < length-1; ++i) {
final int j = Integer.reverse(i) >> rshift;
if (j > i) {
final Complex c = x[i];
x[i] = x[j];
x[j] = c;
}
}
}

final Complex wp = new Complex();
final Complex w = new Complex();
final Complex ws = new Complex();
final Complex c = new Complex();

int mmax = 1;
while (length > mmax) {
final int istep = 2 * mmax;
final double theta = (inverse ? -1.0 : 1.0) * Math.PI/mmax;
final double d = Math.sin(0.5theta);
wp.real = (float) (-2.0 * d * d);
wp.imag = (float) Math.sin(theta);
w.real = 1.0f;
w.imag = 0.0f;
for (int m = 0; m < mmax; ++m) {
// ws = w
ws.real = w.real;
ws.imag = w.imag;
for (int i = m; i < length; i += istep) {
final int j = i + mmax;
// c = ws * x
c.real = ws.realx[j].real - ws.imagx[j].imag;
c.imag = ws.realx[j].imag + ws.imagx[j].real;
// x[j] = x[i] - c
x[j].real = x[i].real - c.real;
x[j].imag = x[j].imag - c.imag;
// x[i] = x[i] + c
x[i].real += c.real;
x[i].imag += c.imag;
}
// c = wwp + w
c.real = w.realwp.real - w.imagwp.imag + w.real;
c.imag = w.realwp.imag + w.imagwp.real + w.imag;
// w = c
w.real = c.real;
w.imag = c.imag;
}
mmax = istep;
}
}

Post generated using Mail2Forum (mail2forum.com)

[uckelman]

Thus spake “Michael Kiefte”:

I wrote the FFT algorithm in GeneralFilter.java only to discover that it’s
mostly pointless. The resampled pixels almost never fall on integer
indeces, so it would require upsampling followed by downsampling in most
cases.

Now I realise why ImageMagick doesn’t use FFT convolution. However, it’s
still good for mipmaps, because the resampling is always a factor of two.
Here’s the code I wrote in the end based on the original FORTRAN. I’ll put
it in the texture demo and see if it works.

One thing I hadn’t thought of before is that if you’re downsampling
by a known, fixed factor, then you can use integers instead of floats
for some of the arithmetic.

Note that Math.sin() is only called in the outermost loop.

That looks pretty tight to me.

Let me know when you’ve uploaded this.

–
J.

---

Messages mailing list
Messages@forums.vassalengine.org
forums.vassalengine.org/mailman/ … engine.org

Post generated using Mail2Forum (mail2forum.com)

[uckelman]

Thus spake “Michael Kiefte”:

I wrote the FFT algorithm in GeneralFilter.java only to discover that it’s
mostly pointless. The resampled pixels almost never fall on integer
indeces, so it would require upsampling followed by downsampling in most
cases.

Now I realise why ImageMagick doesn’t use FFT convolution. However, it’s
still good for mipmaps, because the resampling is always a factor of two.
Here’s the code I wrote in the end based on the original FORTRAN. I’ll put
it in the texture demo and see if it works.

Note that Math.sin() is only called in the outermost loop.

Remember to throw the source into the next bundle, too.

–
J.

---

Messages mailing list
Messages@forums.vassalengine.org
forums.vassalengine.org/mailman/ … engine.org

Post generated using Mail2Forum (mail2forum.com)

[mkiefte]

Let me know when you’ve uploaded this.

It won’t be in the next demo. It’s a lot of work, so I’ve switched temporarily to displaying a much larger map.

• M.

Post generated using Mail2Forum (mail2forum.com)