Major G3D library and related modules code review request

[huhabla]

I would like to request a major code review of the g3d library and the related modules to clarify the important question "What kind of cube coordinate system is used in the g3d lib".

While digging in the code of the g3d library i noticed that the function g3d_getValue() uses a right hand side coordinate system (y achses points in north direction) while the left hand coordinate system code was commented out by AV????:

void G3d_getValue(G3D_Map * map, int x, int y, int z, void *value, int type)
{
    double north, east, top;

     /*AV*/
	/* BEGIN OF ORIGINAL CODE */
	/*
	   int row, col, depth;
	 */
	/* END OF ORIGINAL CODE */
	 /*AV*/
	/* BEGIN OF MY CODE */
    double row, col, depth;

    /* END OF MY CODE */

    /* convert (x, y, z) into (north, east, top) */

     /*AV*/
	/* BEGIN OF ORIGINAL CODE */
	/*
	   north = ((double) map->window.rows - y - 0.5) / (double) map->window.rows *
	   (map->window.north - map->window.south) + map->window.south;
	 */
	/* END OF ORIGINAL CODE */
	 /*AV*/
	/* BEGIN OF MY CODE */
	north = ((double)y + 0.5) / (double)map->window.rows *
	(map->window.north - map->window.south) + map->window.south;
    /* END OF MY CODE */

    east = ((double)x + 0.5) / (double)map->window.cols *
	(map->window.east - map->window.west) + map->window.west;
    top = ((double)z + 0.5) / (double)map->window.depths *
	(map->window.top - map->window.bottom) + map->window.bottom;

    /* convert (north, east, top) into (row, col, depth) */

     /*AV*/
	/* BEGIN OF ORIGINAL CODE */
	/*
	   row = map->region.rows -
	   (north - map->region.south) / (map->region.north - map->region.south) *
	   map->region.rows;
	 */
	/* END OF ORIGINAL CODE */
	 /*AV*/
	/* BEGIN OF MY CODE */
	row =
	(north - map->region.south) / (map->region.north -
				       map->region.south) * map->region.rows;
    /* END OF MY CODE */

    col = (east - map->region.west) / (map->region.east - map->region.west) *
	map->region.cols;
    depth =
	(top - map->region.bottom) / (map->region.top -
				      map->region.bottom) *
	map->region.depths;

    /* if (row, col, depth) outside window return NULL value */
    if ((row < 0) || (row >= map->region.rows) ||
	(col < 0) || (col >= map->region.cols) ||
	(depth < 0) || (depth >= map->region.depths)) {
	G3d_setNullValue(value, 1, type);
	return;
    }

    /* get value */
    map->resampleFun(map, (int)row, (int)col, (int)depth, value, type);

}

But the function G3d_getRegionValue() uses a left hand coordinate system (y achses points in south direction):

void
G3d_getRegionValue(G3D_Map * map, double north, double east, double top,
		   void *value, int type)
{
    int row, col, depth;

    /* convert (north, east, top) into (row, col, depth) */

    row = map->region.rows -
	(north - map->region.south) / (map->region.north -
				       map->region.south) * map->region.rows;
    col =
	(east - map->region.west) / (map->region.east -
				     map->region.west) * map->region.cols;
    depth =
	(top - map->region.bottom) / (map->region.top -
				      map->region.bottom) *
	map->region.depths;

    /* if (row, col, depth) outside window return NULL value */
    if ((row < 0) || (row >= map->region.rows) ||
	(col < 0) || (col >= map->region.cols) ||
	(depth < 0) || (depth >= map->region.depths)) {
	G3d_setNullValue(value, 1, type);
	return;
    }

    /* get value */
    map->resampleFun(map, row, col, depth, value, type);
}

Besides of that the manual pages of r3.in.ascii and r3.out.ascii describe opposite positions about the raster ascii import compatibility.

A left hand side cube coordinate system (y-achses points to the South) is compatible with the row approach of the raster library in grass. The conversion of the y coordinate in northing and vis versa is:

Northing = South + (Rows - Y) / Rows * (North - South)
Y = rows - (Northing - South) / (North - South) * Rows

In a right hand side cube coordinate system (y-achses points to the North) the X and Z coordinate computation is equal the left hand side, but the Y coordinate to northing and vis versa is different:

Northing = South + Y / Rows * (North - South)
Y = (Northing - South) / (North - South) * Rows

I fear that a right hand side cube coordinate system is used in g3d lib and all modules converting voxel into raster and vis versa a implemented wrong.

I will implement some tests to get a clear picture of the g3d library behavior.

Any hints and suggestions are very welcome.

Cheers Soeren

[huhabla]

I have updated the g3d library and all related modules to use the same coordinate system approach. The g3d library should now consistently support the raster library compatible north -> south row order. All g3d coordinate transform functions, g3d getValue() functions and all g3d related modules are now using the following apprach:

• Columns (x coordinate) count from west to east
• Rows (y coordinate) count from north to south
• Depths (z coordinate) count from bottom to top

Several g3d coordinate transform function have been reviewed and reimplemented. I have updated the documentation and implemented several tests for many modules and the g3d library to assure correct behavior.

Additionally three new modules have been implemented:

• r3.retile to modify the tile size of g3d maps
• r3.colors to modify g3d map color tables
• r3.colors.out for g3d map color export

Review finished. :) Soeren