[GRASS5] Gaussian filter for GRASS

[Glynn Clements]

Jachym Cepicky wrote:

> > > I have tryed to program the Gaussian filter for raster maps. It
> > > works only with CELL maps for now, but I hope to extend it for DCELL
> > > and FCELL maps too.
> > 
> > The simplest solution is to use DCELL throughout. All CELL and FCELL
> > values can be represented as DCELL without loss of range or precision.
> 
> OH! I didn't know this! I have to read it in the Programers manual - I
> thought, DCELL is just 'some' another type of map :-( 

No, they're just aliases for the standard C types:

	CELL	int
	FCELL	float
	DCELL	double

> > > Because this is my first raster module for GRASS and also my first bigger
> > > C-program, I would be very thankfull, if someone could have a look at the 
> > > code at to tell me, what I could do better/safer.
> > 
> > Use a scrolling window; don't read each row three times.
> > 
> > At the moment, you're reading rows 0/1/2, then 1/2/3, then 2/3/4 and
> > so on. Note that you're reading each row (except the first and last
> > two) three times. Also note that you already have two out of the three
> > rows in memory from the previous pass (again, except for the first and
> > last two passes).
> > 
> > See r.neighbors or r.grow2.
> 
> I looked there, but I did not see it. Thanks again...

Look at how in_rows is used in r.grow2. [OK, reading further, you
already have.]

For r.neighbors, the equivalent is ncb.buf; however, r.neighbors is
split across several files, so the logic is a bit harder to follow.

> > > I could not to solve, how to work on the boundaries of the raster. I had to
> > > skip these cells, like e.g. r.mapcalc does -- so NULL values appear there. 
> > > How to do it some better way?
> > 
> > In general, you can't. Convolution-style filters always result in the
> > map being "shrunk".
> 
> But one could count the value of the central cell only from the values, which
> are laying in the raster
> 
> +-----------+
> | 1| 2| 3| 4|   So for cell number 1, one would take the cells 1,2,5 and 6 
> +-----------+   in the count. For number 2, One would take 1,2,3,5,6,7
> | 5| 6| 7| 8|   first for cell number 6, one would take all 9 values around.
> +-----------+  
> | 9|10|11|12|   But it would be too complicated tough :-/
> +-----------+

For simple aggregates which don't care about the relative positions of
the cells (sum, mean etc), you can get an answer, but it isn't
necessarily meaningful (e.g. for a sum, the values of the edge cells
will tend to be lower because you're summing fewer values).

For convolution filters, silently "ignoring" missing cells tends to
produce the same result as if the cell's value was zero. Usually this
doesn't make sense.

More realistic would be to define separate weight matrices for each
case, but this is usually too involved. Or you could scale the weights
by the total number of cells, then divide by the actual number of
non-NULL cells; that will tend to give reasonable results for typical
smoothing filters.

Whatever approach you use should treat "stored" NULLs the same as
cells which are NULL by virtue of being outside of the current region. 
I.e. if you were to enlarge the region, filling the boundary cells
with NULLs, you should get the same result.

r.grow2 just fills input cells with NULL if they are outside of the
region, so there is guaranteed to be no difference with stored NULLs.

> OK, I will work with DCELL type, than I implement the row swiching from r.grow2
> 
> 		tmp = in_rows[0];
> 		for (i = 0; i < size * 2; i++)
> 			in_rows[i] = in_rows[i + 1];
> 		in_rows[size * 2] = tmp;

Yeah, that's the essence of a scrolling window.

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