Drawing with Fourier descriptors?

November 5, 2006 2:55 PM

I need help finding code to draw (not extract) fourier descriptors.

I need to draw some simple FD-specified shapes (basically little blobs - power at only a few frequencies). I'm having a lot of trouble doing it myself, and it seems that this isn't something that anyone else ever does - there's lots of code out there for taking an existing contour and fitting it, but I need to go in the other direction. Anyone have some code or know where I can find it? I'm language-agnostic; I'm happy to translate from whatever I get into the language I'm working in.

posted by dmd to science & nature (15 comments total)

Based on this page, it looks like there are at least two kinds of FDs. One describing angle delta, and one using contour points directly via complex numbers.
The second page (linked from the bottom of the page linked above) has a formula (using sigma notation) for converting back and forth between the latter (complex) ones. Looks really simple to translate into code.
posted by blenderfish at 3:33 PM on November 5, 2006

I may be confused about what you want, but it seems like you should be able to do what you want with Excel, or any other plotting spreadsheet program.

Create three columns. The first column is just rising numbers (integers, or fractional steps of integers). The second and third columns are functions based on the first column.

Then do an X-Y plot of the second and third columns against one another.
posted by Steven C. Den Beste at 3:56 PM on November 5, 2006

Ok, let me be a little more specific. I guess I'm a little unclear on the math, and I should probably become more clear on it.

I'm looking at that and I'm not entirely sure how to relate it to the inputs that I have. I want to specify a list of frequencies (given in cycles per perimeter - say, 2, 4, and 6 cycles per perimeter), and for each frequency, I have a phase and an amplitude (amplitude is given in radians, and represents the orientation of a vector tangent to the countour of the shape - cumulative angular bend).

Does that make any sense?

If it helps any, I'm trying to recreate the shapes used in this paper: Perceptual Similarity of Shapes (see page 3).
posted by dmd at 4:09 PM on November 5, 2006

I don't know if this helps, but an outline of a naive (that is, sub-optimal in speed, memory, etc etc etc.) algorithm would be to take the inverse fourer series of the data that you have, to recover an approximation of the original function that describes the angular relationship between the tangent at the origin, and the tangent s units along the curve.

Once you have this function, it shouldn't be too difficult to plot the curve.

If you're unfamiliar with fourier analysis, there are lots of resources on the net, and I heartily recommend reading a few pages in Oppenheim and Shafer, Signals and Systems, the cannoical reference for learning about Fourier analysis.
posted by scalespace at 4:48 PM on November 5, 2006

Did you look up the reference?

Zahn, C. T., & Roskies, R. Z. (1972). Fourier descriptors for plane closed curves. IEEE Transactions on Computers, C-21, 269281.
posted by nonmyopicdave at 8:13 PM on November 5, 2006

I'm still waiting for my library to retrieve the reference from storage. Apparently it's too old to have in general stacks.

But I think taking the inverse fourier is the 'aha' that I needed. Matlab should be able to do that, and whatever that gives me I bet will be amenable to my further tweakings.
posted by dmd at 5:10 AM on November 6, 2006

You know what? Looking at the shapes and descriptors (and not looking at the reference, which probably has the answer) I don't think the inverse fourier is it. I think it may be much simpler than you think. I think your cumulative angular bend function f(t, amp, phase) is just given by:

f(t, amp, phase) = 0.5(sin(2t) + cos(2t) + sin(4t) + cos(4t)) + amp*(sin((6+phase)*t) + cos((6+phase)*t))

t is your parameter that runs from 0 to 2*pi.
The nine shapes come from amp = 0.5, 0.75, 1.00 and phase = 0, 1/8, 1/4.

This is the cumulative angular bend function shown in figure one.

Once you have that, you step t from 0 to 2*pi and write another function that spits out coordinates based on the location and change in angle from the last point.

Again, that's just a guess based on what the shapes look like. I think you'll find the definitive answer in that paper, though.
posted by nonmyopicdave at 12:23 PM on November 6, 2006

Actually, I realized I gave the wrong pair of authors for Signals and Systems; it's actually Oppenheim and Willsky.
posted by scalespace at 12:44 PM on November 6, 2006

I'm not sure about the values for amp.

fracphase = phase/360;
theta = 0.5*(sin(2*t)+cos(2*t)) + 0.5*(sin(4*t)+cos(4*t)) + amp*(sin((6+fracphase)*t)+cos((6+fracphase)*t));


Putting in 0.5 and 45 gives this plot, which doesn't look at all right to me.
posted by dmd at 9:17 AM on November 7, 2006

never mind - I got it!

It also needed a -t term.

Thank you! it's working!
posted by dmd at 11:00 AM on November 7, 2006

Solved really now:

function [ theta ] = cumubend( t, FD )
%cumubend cumulative angular bend function. 
%  cumubend(t,FD) is the theta value of the 
%      cumulative angular bend function at t

theta = -t;

for freq = 1:length(FD(:,1)),
    amp = FD(freq,1);
    phase = FD(freq,2) / 360 * 2* pi;
    theta = theta + amp * cos(freq * t - phase);
end



function [ points ] = cumubend2points(fd,steps)

points = [0];
thispoint = [0];

for t = 0:2*pi/steps:2*pi,
    bendangle = cumubend(t,fd);
    nextpoint = cos(bendangle) + i * sin(bendangle) + thispoint;
    points = [points nextpoint];
    thispoint = nextpoint;
end

posted by dmd at 3:35 PM on November 8, 2006

Sweet! Can you post an image of the shapes you got?
posted by nonmyopicdave at 12:17 AM on November 11, 2006

some images
posted by dmd at 9:07 AM on November 11, 2006

Excellent! Well done.

I'm embarrassed I wrote (6+phase)*t instead of (6*t+phase) but hopefully it still helped!
posted by nonmyopicdave at 1:09 PM on November 11, 2006

looks awesome!
posted by scalespace at 9:36 PM on November 11, 2006