Deformable registration in D

[Don Clugston]

Robert Jacques wrote:
> I recently presented at the 49th Annual Meeting of the American 
> Association of Physicists in Medicine an application of a deformable 
> registration algorithm during the general poster discussion session. Of 
> course, the reason I?m posting this here is that I wrote it in D. Or 
> more precisely, I refactored it to D from C++, with some surprising 
> results. First, the main algorithm dropped to half the lines of code and 
> I was able to completely remove some supporting routines. Second, the 
> numerical stability of the algorithm increased, helping reduce a major 
> fundamental issue/flaw in algorithm (Ah, the joys of a continuous 
> algorithm working on discrete data). And lastly, it ran three to four 
> times slower. The culprit turned out to be the exponential function, 
> which is the crux of the algorithm and also what I assume was 
> responsible for the increased stability. I switched to using an 
> exponential with finite support and regained the speed while keeping the 
> stability. For reference, I was using Microsoft Visual Studio 2003 and 
> D1.x/Phobos on Intel/XP box.
> The code is now open source and available at:
> http://dailabs.duhs.duke.edu/imagequality.html
> And the associated abstract is published in Medical Physics:
> http://www.medphys.org/
> 
> P.S. Thank you Walter and the community for this great language.

Awesome!
A few comments:
(1) D's exp() shouldn't be too much slower than C's. Certainly not a factor of 
four. (Unless you were creating a lot of denormalised numbers). So your result 
is interesting. Perhaps there is a performance bug in exp().
(2) DMD does essentially no floating-point optimisation. The fact that it can 
compete directly with Visual C++ is impressive, since there is enormous 
potential for improvement (a factor of 2 at the very least).
(3) The increased stability could also be due to the use of 80-bit reals instead 
of 64 bits. You could easily test this by putting

import std.c.fenv;

fesetprec(FE_DBLPREC);

into main(), and seeing if the stability disappears.

(4) My experience has been that D is a superb language for developing 
floating-point algorithms. It's great to see further confirmation of this.

I presume that you were using DMD, not GDC. If not, it would explain (1) and (2).