Opened 10 years ago

Closed 9 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: master
Keywords: history Cc:

Description

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

Change History (14)

comment:1 by Bborie Park, 10 years ago

Related tickets are…

#2310 #2311 #2312

comment:2 by Bborie Park, 10 years ago

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 by Bborie Park, 10 years ago

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 by nclay, 10 years ago

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 by pracine, 10 years ago

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

in reply to:  4 ; comment:6 by Bborie Park, 10 years ago

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)

in reply to:  5 comment:7 by Bborie Park, 10 years ago

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 by Bborie Park, 10 years ago

Status: assignednew

in reply to:  6 comment:9 by nclay, 9 years ago

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 by Bborie Park, 9 years ago

Yes, that would work. I forgot about rt_pixel_set_to_array…

comment:11 by nclay, 9 years ago

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 by Bborie Park, 9 years ago

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 by nclay, 9 years ago

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 by Bborie Park, 9 years ago

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.