[raster] ST_MapAlgebra mask parameter

[nclay]

MapAlgebra should accept a mask to allow for focal statistics and other weighted mask operations.

[Bborie Park]

Related tickets are...

#2310 #2311 #2312

[Bborie Park]

New variants using a mask and weight flag parameters instead of distancex and distancey.

weight parameter is a boolean flag indicating that the values in mask are weights that should be applied to the pixel values.

The mask would be a 2D double precision array of the neighborhood. The length along the X and Y axes must be at least one AND an odd number. The length must be odd as otherwise, where in the mask is the element for the center pixel (pixel of interest)?

A valid mask...

mask = ARRAY[
  [1, 0, 1],
  [0, 1, 0]
  [1, 0, 1]
]

An invalid mask...

mask = ARRAY[
  [1, 0],
  [0, 1],
  [1, 0]
]

In the default case where weight = FALSE, the above valid mask of zeros and ones indicate which pixel values should be passed to the callback function. The values zero or NULL are treated as FALSE.

If weight = TRUE, the values in mask are the weight by which the corresponding pixel values are multiplied. If a weight is NULL, the pixel value is set to NULL.

The basic function should look like...

ST_MapAlgebra(
  rastbandargset rastbandarg[],
  callbackfunc regprocedure,
  pixeltype text=NULL,
  extenttype text=INTERSECTION, customextent raster=NULL,
  mask double precision[][]=NULL, weighted boolean=FALSE,
  VARIADIC userargs text[]=NULL
);

If mask = NULL (or empty array), the center pixel value is passed as-is to the callback function.

[Bborie Park]

Existing distancex and distancey variants will be changed to run on top of the mask variants since the mask ones are more flexible.

[nclay]

That is what I had in mind. I am sill working on it and should have it ready for testing, by Monday. I have the raster iterator done and its pixel fetcher done that considers of a mask. When you say that "Existing distancex and distancey variants will be changed to run on top of the mask variants since the mask ones are more flexible." My current code has it where if the Mask parameter is null then the "old code runs", and if a Mask is provided the 'new code is run'. Are you saying that, the code should be cleaned up so that we create a mask,

mask = ARRAY[
  [1, 1, 1],
  [1, 1, 1]
  [1, 1, 1]
]

To be passed mimicking the old behaviour? Or is the current branching behaviour acceptable?

[pracine]

Why you don't just pass the mask as a userargs argument?

[Bborie Park]

Replying to nclay:

That is what I had in mind. I am sill working on it and should have it ready for testing, by Monday. I have the raster iterator done and its pixel fetcher done that considers of a mask. When you say that "Existing distancex and distancey variants will be changed to run on top of the mask variants since the mask ones are more flexible." My current code has it where if the Mask parameter is null then the "old code runs", and if a Mask is provided the 'new code is run'. Are you saying that, the code should be cleaned up so that we create a mask,

mask = ARRAY[
  [1, 1, 1],
  [1, 1, 1]
  [1, 1, 1]
]

To be passed mimicking the old behaviour? Or is the current branching behaviour acceptable?

How are you defining the SQL-level ST_MapAlgebra() function to do the branching? If you're just appending the mask parameter to existing function signatures, that is messy and redundant.

As part of what I envision, the changes would also require changing RASTER_nMapAlgebra() in rt_pg/rt_pg.c and rt_raster_iterator() in rt_core/rt_api.c to operate upon masks instead of distances. Because of the deep changes, I don't think a separate band pixels by mask getter function is needed. A alternative approach would be a general utility function that operates on a values array based upon a mask array.

So, crudely thought out.

rt_errorstate
rt_util_pixel_mask(rt_pixel *pixels, int npixels, double precision **mask, int dimx, int dimy, int weighted)

[Bborie Park]

Replying to pracine:

Why you don't just pass the mask as a userargs argument?

Actually, I already do this for my code. The problem is that the conversion from text to whatever datatype is a performance hit. For something as generic as a mask (distancex and distancey basically creates a mask), might as well have it in C.

[nclay]

Replying to dustymugs:

Replying to nclay:

That is what I had in mind. I am sill working on it and should have it ready for testing, by Monday. I have the raster iterator done and its pixel fetcher done that considers of a mask. When you say that "Existing distancex and distancey variants will be changed to run on top of the mask variants since the mask ones are more flexible." My current code has it where if the Mask parameter is null then the "old code runs", and if a Mask is provided the 'new code is run'. Are you saying that, the code should be cleaned up so that we create a mask,

mask = ARRAY[
  [1, 1, 1],
  [1, 1, 1]
  [1, 1, 1]
]

To be passed mimicking the old behaviour? Or is the current branching behaviour acceptable?

How are you defining the SQL-level ST_MapAlgebra() function to do the branching? If you're just appending the mask parameter to existing function signatures, that is messy and redundant.

As part of what I envision, the changes would also require changing RASTER_nMapAlgebra() in rt_pg/rt_pg.c and rt_raster_iterator() in rt_core/rt_api.c to operate upon masks instead of distances. Because of the deep changes, I don't think a separate band pixels by mask getter function is needed. A alternative approach would be a general utility function that operates on a values array based upon a mask array.

So, crudely thought out.

rt_errorstate
rt_util_pixel_mask(rt_pixel *pixels, int npixels, double precision **mask, int dimx, int dimy, int weighted)

What about considering the mask in rt_pixel_set_to_array? That would save us some iterations over the pixels.

[Bborie Park]

Yes, that would work. I forgot about rt_pixel_set_to_array...

[nclay]

Nulls?

How should we handle null as a double can not be null in c? In my current code I treat 0 as null for a non-weighted mask, and ignore the possibility for a null in a weighted mask.

             if ( mask == NULL ) {
		  values[_y][_x] = npixel[i].value;
		  nodatas[_y][_x] = 0;
		}else{ 
		  if( !weighted ){
		    if( mask[_y][_x] == 0 ){
		      values[_y][_x] = 0;
		      nodatas[_y][_x] = 1;
		    }else{
		      values[_y][_x] = npixel[i].value;
		      nodatas[_y][_x] = 0;
		    }
		  }else{
		    values[_y][_x] = npixel[i].value * mask[_y][_x];
		    nodatas[_y][_x] = 0;
		  }
		}

What value should be treated as a NULL in the case of a weighted mask? I am thinking 0, but I want some thought on this or should we have a mask nodata array?

[Bborie Park]

You will need an array for mask NODATA.

In your code snippet, if mask is of type double, don't ever evaluate like...

mask[_y][_x] == 0

Instead use...

FLT_EQ(mask[_y][_x], 0)
}}

[nclay]

I have created a pull Request on Git Hub for this ticket #2341 (​http://github.com/radioprop/postgis). Please review.

Thanks,

Nathaniel Hunter Clay

[Bborie Park]

Merged ​https://github.com/postgis/postgis/pull/17 with minor memory leak fix in cunits. r12179 in -trunk.

I will open a new ticket to add docs for the ST_MapAlgebra with mask variant.