[GRASSLIST:9635] Re: filtering with FFT

[Glynn Clements]

Markus Neteler wrote:

> I have made these steps
> 
>  g.region rast=tm1
>  i.fft in=tm1 real=tm1_r imag=tm1_i --o
> # g.region rast=tm1_r
> # r.circle -b out=ff coord=-769058,-988168 max=1550 --o
>  r.circle -b out=ff coord=-769058,-988168 max=10550 --o
>  d.rast tm1_r
>  d.rast -o ff
>  
>  r.mapcalc 'mymask=if(isnull(ff),1,null())'
>  g.rename rast=mymask,MASK
>  i.ifft real=tm1_r imag=tm1_i out=tm1_hp --o
>  g.remove rast=MASK
> # g.region rast=tm1
>  r.colors map=tm1_hp col=grey.eq
>  d.rast tm1
>  d.rast tm1_hp
> 
> and get the same ugly result.
> 
> I suspect that 'krovak' with it's negative coordinates may be
> the case to let misbehave i.[i]fft. 
> 
> You may reproject to an UTM location and try again.
> If this works, we have to understand what i.[i]fft doesn't
> work with negative coordinates. If it doesn't work, the
> problem is more generic.

It might be useful to mention that i.ifft doesn't do what one might
normally assume, in that it doesn't use the actual maps which are
specified by the real= and imag= options in the way that any other
module would.

Instead of reading the maps using the normal libgis functions
(G_open_cell_old, G_get_*_row etc), it reads the data from the
fftreal/fftimag files in the maps' cell_misc subdirectories, and
handles the region and mask itself (libgis does this part
automatically when reading raster maps).

It isn't entirely clear how the region or mask are meant to be
interpreted, given that the fftreal/fftimag files aren't in any
cartographic coordinate system.

Arguably, the results of i.fft (and the inputs to i.ifft) should live
in a completely separate location, typically an X-Y location with a
region centred around the origin of the complex plane (although the
units would still be metres or feet).

I suspect the main reason it doesn't work that way is that it is
awkward to use multiple locations concurrently (the situation has been
improved recently, but you still can't e.g. manage data from multiple
locations in a single d.m instance). Also, 4.x didn't support FP maps,
so the FFT coefficient arrays couldn't realistically have been
represented as maps back then.

-- 
Glynn Clements <glynn at gclements.plus.com>