[GRASS-user] Gauss filter in r.neighbors

GRASS GRASS user
1

Hello!

I am using the the r.neighbors function and I have an idea about what the ‘Gauss’ parameter will do for me, but I can’t find any documentation for it. Experimenting with it, left me clueless.

What I won’t to do:
Apply a Gaussian weighting of cell neighbors and sum the values of the neighbors. If all cells=1, I would get only 1’s back. I have a map of zeros and ones. I tried this:

r.neighbors -c input=‘input_map’ output=‘output_map’ method=sum gauss=10 size=33

The result is a map with only zeros. So, my question is: What does ‘gauss’ do? I’ am guessing I am using it wrong but I would like to know.

Best regards,

Martin

2

Martin Album Ytre-Eide wrote:

I am using the the r.neighbors function and I have an idea about
what the 'Gauss' parameter will do for me, but I can't find any
documentation for it. Experimenting with it, left me clueless.

Using gauss= is equivalent to using weight= with a weights file
corresponding to a Gaussian filter.

The exact weights which are used can be found in the
gaussian_weights() function in raster/r.neighbors/readweights.c. They
are effectively:

  ncb.weights[i][j] = exp(-(x*x+y*y)/(2*sigma2))/(2*M_PI*sigma2);

where x, y is the offset (in cells) from the centre of the
neighbourhood and sigma2 is the square of the value passed to the
gauss= option.

What I won't to do:

Apply a Gaussian weighting of cell neighbors and sum the values of
the neighbors. If all cells=1, I would get only 1's back.

You need to use method=average to make the effective weights sum to
one (when weights are used, method=average divides the weighted sum of
the cell values by the sum of the weights).

IIRC, the weights are such that they would sum to one if you had an
infinitely-large neighbourhood, but the fact that you're cropping the
curve to a finite neighbourhood means that the sum will always be less
than one.

The weights are determined by sigma alone, and are unaffected by the
neighbourhood size.

I have a map of zeros and ones. I tried this:

r.neighbors -c input='input_map' output='output_map' method=sum gauss=10 size=33

The result is a map with only zeros.

It works for me; e.g.

  $ r.mapcalc --o 'foo = rand(0,2)' seed=1
  $ r.neighbors --o -c in=foo out=bar method=average size=33 gauss=10
  $ r.info -r bar
  min=0.442173333781042
  max=0.568925311919939
  $ r.neighbors --o -c in=foo out=bar method=sum size=33 gauss=10
  $ r.info -r bar
  min=0.0985901293262625
  max=0.449425795562072

  $ r.mapcalc --o 'foo = 1'
  $ r.neighbors --o -c in=foo out=bar method=average size=33 gauss=10
  $ r.info -r bar
  min=1
  max=1
  $ r.neighbors --o -c in=foo out=bar method=sum size=33 gauss=10
  $ r.info -r bar
  min=0.221413499330273
  max=0.812157275619236

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

3

Hello and thanks a lot!

I followed you examples and it did not work in grass 6.4.4 so I tried them in grass7.1 and they worked for me as well.

Tried my map - and it worked. I result looks good :slight_smile:

Thanks!

(maybe there is a problem in r.neighbors for grass 6.4.4?)

Martin

________________________________________
Fra: Glynn Clements [glynn@gclements.plus.com]
Sendt: 2. oktober 2014 14:27
Til: Martin Album Ytre-Eide
Kopi: grass-user@lists.osgeo.org
Emne: Re: [GRASS-user] Gauss filter in r.neighbors

Martin Album Ytre-Eide wrote:

I am using the the r.neighbors function and I have an idea about
what the 'Gauss' parameter will do for me, but I can't find any
documentation for it. Experimenting with it, left me clueless.

Using gauss= is equivalent to using weight= with a weights file
corresponding to a Gaussian filter.

The exact weights which are used can be found in the
gaussian_weights() function in raster/r.neighbors/readweights.c. They
are effectively:

        ncb.weights[i][j] = exp(-(x*x+y*y)/(2*sigma2))/(2*M_PI*sigma2);

where x, y is the offset (in cells) from the centre of the
neighbourhood and sigma2 is the square of the value passed to the
gauss= option.

What I won't to do:

Apply a Gaussian weighting of cell neighbors and sum the values of
the neighbors. If all cells=1, I would get only 1's back.

You need to use method=average to make the effective weights sum to
one (when weights are used, method=average divides the weighted sum of
the cell values by the sum of the weights).

IIRC, the weights are such that they would sum to one if you had an
infinitely-large neighbourhood, but the fact that you're cropping the
curve to a finite neighbourhood means that the sum will always be less
than one.

The weights are determined by sigma alone, and are unaffected by the
neighbourhood size.

I have a map of zeros and ones. I tried this:

r.neighbors -c input='input_map' output='output_map' method=sum gauss=10 size=33

The result is a map with only zeros.

It works for me; e.g.

        $ r.mapcalc --o 'foo = rand(0,2)' seed=1
        $ r.neighbors --o -c in=foo out=bar method=average size=33 gauss=10
        $ r.info -r bar
        min=0.442173333781042
        max=0.568925311919939
        $ r.neighbors --o -c in=foo out=bar method=sum size=33 gauss=10
        $ r.info -r bar
        min=0.0985901293262625
        max=0.449425795562072

        $ r.mapcalc --o 'foo = 1'
        $ r.neighbors --o -c in=foo out=bar method=average size=33 gauss=10
        $ r.info -r bar
        min=1
        max=1
        $ r.neighbors --o -c in=foo out=bar method=sum size=33 gauss=10
        $ r.info -r bar
        min=0.221413499330273
        max=0.812157275619236

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

4

Martin Album Ytre-Eide wrote:

(maybe there is a problem in r.neighbors for grass 6.4.4?)

Yes:

Markus Metz wrote:

There is no problem with the Gaussian filter, but a problem with the
output map type of r.neighbours in G64: the output type is in G64 the
same like the input type, whereas in G7 the output type depends on the
chosen function. That means, if you use r.neighbours method=average in
G64 with a CELL map, the output is also CELL.

If you need it to work with 6.x, convert the input to FCELL/DCELL
first (with e.g. r.mapcalc "output = float(input)").

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