[GRASS-SVN] r61569 - grass/trunk/lib/raster3d
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Aug 8 14:31:07 PDT 2014
Author: annakrat
Date: 2014-08-08 14:31:07 -0700 (Fri, 08 Aug 2014)
New Revision: 61569
Modified:
grass/trunk/lib/raster3d/gradient.c
Log:
r3.gradient: set voxels which have neighbors nulls to zero
Modified: grass/trunk/lib/raster3d/gradient.c
===================================================================
--- grass/trunk/lib/raster3d/gradient.c 2014-08-08 18:44:59 UTC (rev 61568)
+++ grass/trunk/lib/raster3d/gradient.c 2014-08-08 21:31:07 UTC (rev 61569)
@@ -18,7 +18,9 @@
Gradient computation (second order approximation)
using central differencing scheme (plus forward and backward
- difference of second order approx.)
+ difference of second order approximation). When one or more of the cells,
+ from which the gradient for a particular cell is computed, is null,
+ gradient for that particular cell is set to 0.
\param array pointer to RASTER3D_Array with input values
\param step array of x, y, z steps for gradient (resolution values)
@@ -35,115 +37,142 @@
RASTER3D_Array_double *grad_z)
{
int col, row, depth;
- double val;
+ double val0, val1, val2;
for (depth = 0; depth < array->sz; depth++) {
for (row = 0; row < array->sy; row++) {
/* row start */
- val = RASTER3D_ARRAY_ACCESS(array, 0, row, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, 0, row, depth);
+ val1 = RASTER3D_ARRAY_ACCESS(array, 1, row, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, 2, row, depth);
+ if (Rast_is_d_null_value(&val0))
Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
+ RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_x, 0, row, depth) =
- (-3 * val + 4 * RASTER3D_ARRAY_ACCESS(array, 1, row, depth) -
- RASTER3D_ARRAY_ACCESS(array, 2, row, depth)) / (2 * step[0]);
+ (-3 * val0 + 4 * val1 - val2) / (2 * step[0]);
/* row end */
- val = RASTER3D_ARRAY_ACCESS(array, array->sx - 1, row, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, array->sx - 3, row, depth);
+ val1 = RASTER3D_ARRAY_ACCESS(array, array->sx - 2, row, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, array->sx - 1, row, depth);
+ if (Rast_is_d_null_value(&val2))
Rast_set_null_value(
&RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
+ RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_x, array->sx - 1, row, depth) =
- (3 * val - 4 * RASTER3D_ARRAY_ACCESS(array, array->sx - 2, row, depth) +
- RASTER3D_ARRAY_ACCESS(array, array->sx - 3, row, depth)) / (2 * step[0]);
+ (3 * val2 - 4 * val1 + val0) / (2 * step[0]);
/* row */
for (col = 1; col < array->sx - 1; col++) {
- val = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col - 1, row, depth);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col + 1, row, depth);
+ if (Rast_is_d_null_value(&val1))
Rast_set_null_value(
&RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val2))
+ RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_x, col, row, depth) =
- (RASTER3D_ARRAY_ACCESS(array, col + 1, row, depth) -
- RASTER3D_ARRAY_ACCESS(array, col - 1, row, depth)) / (2 * step[0]);
+ (val2 - val0) / (2 * step[0]);
}
}
}
for (depth = 0; depth < array->sz; depth++) {
for (col = 0; col < array->sx; col++) {
/* col start */
- val = RASTER3D_ARRAY_ACCESS(array, col, 0, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col, 0, depth);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, 1, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col, 2, depth);
+ if (Rast_is_d_null_value(&val0))
Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
+ RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_y, col, 0, depth) =
- -(-3 * val + 4 * RASTER3D_ARRAY_ACCESS(array, col, 1, depth) -
- RASTER3D_ARRAY_ACCESS(array, col, 2, depth)) / (2 * step[1]);
+ -(-3 * val0 + 4 * val1 - val2) / (2 * step[1]);
/* col end */
- val = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 1, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 3, depth);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 2, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col, array->sy - 1, depth);
+ if (Rast_is_d_null_value(&val2))
Rast_set_null_value(
&RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
+ RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_y, col, array->sy - 1, depth) =
- -(3 * val - 4 * RASTER3D_ARRAY_ACCESS(array, col, array->sy - 2, depth) +
- RASTER3D_ARRAY_ACCESS(array, col, array->sy - 3, depth)) / (2 * step[1]);
+ -(3 * val2 - 4 * val1 + val0) / (2 * step[1]);
/* col */
for (row = 1; row < array->sy - 1; row++) {
- val = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col, row - 1, depth);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col, row + 1, depth);
+ if (Rast_is_d_null_value(&val1))
Rast_set_null_value(
&RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val2))
+ RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_y, col, row, depth) =
- -(RASTER3D_ARRAY_ACCESS(array, col, row + 1, depth) -
- RASTER3D_ARRAY_ACCESS(array, col, row - 1, depth)) / (2 * step[1]);
+ -(val2 - val0) / (2 * step[1]);
}
}
}
for (row = 0; row < array->sy; row++) {
for (col = 0; col < array->sx; col++) {
/* vertical col start */
- val = RASTER3D_ARRAY_ACCESS(array, col, row, 0);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col, row, 0);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, row, 1);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col, row, 2);
+ if (Rast_is_d_null_value(&val0))
Rast_set_null_value(&RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val1) || Rast_is_d_null_value(&val2))
+ RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_z, col, row, 0) =
- (-3 * val + 4 * RASTER3D_ARRAY_ACCESS(array, col, row, 1) -
- RASTER3D_ARRAY_ACCESS(array, col, row, 2)) / (2 * step[2]);
+ (-3 * val0 + 4 * val1 - val2) / (2 * step[2]);
/* vertical col end */
- val = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 1);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 3);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 2);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 1);
+ if (Rast_is_d_null_value(&val2))
Rast_set_null_value(
&RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val1))
+ RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_z, col, row, array->sz - 1) =
- (3 * val - 4 * RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 2) +
- RASTER3D_ARRAY_ACCESS(array, col, row, array->sz - 3)) / (2 * step[2]);
+ (3 * val2 - 4 * val1 + val0) / (2 * step[2]);
/* vertical col */
for (depth = 1; depth < array->sz - 1; depth++) {
- val = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
- if (Rast_is_d_null_value(&val))
+ val0 = RASTER3D_ARRAY_ACCESS(array, col, row, depth - 1);
+ val1 = RASTER3D_ARRAY_ACCESS(array, col, row, depth);
+ val2 = RASTER3D_ARRAY_ACCESS(array, col, row, depth + 1);
+ if (Rast_is_d_null_value(&val1))
Rast_set_null_value(
&RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth),
1, DCELL_TYPE);
+ else if (Rast_is_d_null_value(&val0) || Rast_is_d_null_value(&val2))
+ RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) = 0;
else
RASTER3D_ARRAY_ACCESS(grad_z, col, row, depth) =
- (RASTER3D_ARRAY_ACCESS(array, col, row, depth + 1) -
- RASTER3D_ARRAY_ACCESS(array, col, row, depth - 1)) / (2 * step[2]);
+ (val2 - val0) / (2 * step[2]);
}
}
}
More information about the grass-commit
mailing list