Opened 6 years ago

Closed 5 years ago

#2341 closed enhancement (fixed)

[raster] ST_MapAlgebra mask parameter

Reported by: nclay Owned by: Bborie Park
Priority: medium Milestone: PostGIS 2.2.0
Component: raster Version: trunk
Keywords: history Cc:

Description

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

Change History (14)

comment:1 Changed 6 years ago by Bborie Park

Related tickets are...

#2310 #2311 #2312

comment:2 Changed 6 years ago by Bborie Park

Milestone: PostGIS FuturePostGIS 2.2.0
Status: newassigned
Summary: [raster] st_mapalgebra mask parameter[raster] ST_MapAlgebra mask parameter

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.

comment:3 Changed 6 years ago by 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.

comment:4 Changed 6 years ago by 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?

comment:5 Changed 6 years ago by pracine

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

comment:6 in reply to:  4 ; Changed 6 years ago by 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)

comment:7 in reply to:  5 Changed 6 years ago by 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.

comment:8 Changed 6 years ago by Bborie Park

Status: assignednew

comment:9 in reply to:  6 Changed 6 years ago by 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.

comment:10 Changed 6 years ago by Bborie Park

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

comment:11 Changed 6 years ago by 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?

comment:12 Changed 6 years ago by 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)
}}

comment:13 Changed 5 years ago by 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

comment:14 Changed 5 years ago by Bborie Park

Keywords: history added
Resolution: fixed
Status: newclosed

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.

Note: See TracTickets for help on using tickets.