Replaced tcholSolver() with the gmath version

G_math_solver_cholesky_sband(). Removed TcholBand.c. Renamed some
variables to be more descriptive. Added two nc examples from Markus Metz
to the html docs.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@42065 15284696-431f-4ddb-bdfa-cd5b030d7da7

vector/lidar/lidarlib/PolimiFunct.h

@@ -69,11 +69,11 @@
     /*STRUCTS DECLARATION */
     struct Reg_dimens
 {
-    double orlo_h;		/*Horizontal tile edge */
-    double orlo_v;		/*Vertical tile edge */
+    double edge_h;		/*Horizontal tile edge */
+    double edge_v;		/*Vertical tile edge */
     double overlap;		/*Tile's overlapping size */
-    double latoN;		/*South-North side size */
-    double latoE;		/*East-West side size */
+    double sn_size;		/*South-North side size */
+    double ew_size;		/*East-West side size */
 };
 
 struct Point
@@ -102,7 +102,7 @@ int P_set_regions(struct Cell_head *, /**/
 		  struct bound_box *, /**/
 		  struct bound_box *, /**/ struct Reg_dimens, /**/ int /**/);
 
-int P_get_orlo(int, /**/ struct Reg_dimens *, /**/ double, /**/ double /**/);
+int P_get_edge(int, /**/ struct Reg_dimens *, /**/ double, /**/ double /**/);
 
 int P_get_BandWidth(int, /**/ int /**/);
 
@@ -159,16 +159,6 @@ void P_Aux_to_Vector(struct Map_info *, /**/
 
 double **P_Null_Matrix(double ** /**/);
 
-/*----------------------------------------------------------------------------------------------------------*/
-/*tcholBand */
-void tcholDec(double **N, double **T, int n, int BW);
-void tcholSolve(double **N, double *TN, double *parVect, int n, int BW);
-void tcholSolve2(double **N, double *TN, double **T, double *parVect, int n,
-		 int BW);
-void tcholInv(double **N, double *invNdiag, int n, int BW);
-void tcholSolveInv(double **N, double *TN, double *invNdiag, double *parVect,
-		   int n, int BW);
-
 /*---------------------------------------------------------------------------------------*/
 /*interpSpline */
 void normalDefBicubic(double **N, double *TN, double *Q, double **obsVect,

vector/lidar/lidarlib/TcholBand.c

@@ -1,236 +0,0 @@
-#include <stdlib.h>		/* imported libraries */
-#include <stdio.h>
-#include <math.h>
-#include <grass/gis.h>
-#include <grass/PolimiFunct.h>
-
-/*--------------------------------------------------------------------------------------*/
-/* Tcholetsky decomposition -> T= Lower Triangular Matrix */
-
-void tcholDec(double **N, double **T, int n, int BW)
-{
-    int i, j, k, end;
-    double somma;
-
-    G_debug(2, "tcholDec(): n=%d  BW=%d", n, BW);
-
-    for (i = 0; i < n; i++) {
-	G_percent(i, n, 2);
-	for (j = 0; j < BW; j++) {
-	    somma = N[i][j];
-	    /* start = 1 */
-	    /* end = BW - j or i + 1 */
-	    end = ((BW - j) < (i + 1) ? (BW - j) : (i + 1));
-	    for (k = 1; k < end; k++)
-		somma -= T[i - k][k] * T[i - k][j + k];
-	    if (j == 0) {
-		if (somma <= 0.0)
-		    G_fatal_error(_("Decomposition failed"));
-		T[i][0] = sqrt(somma);
-	    }
-	    else
-		T[i][j] = somma / T[i][0];
-	}
-    }
-
-    G_percent(i, n, 2);
-    return;
-}
-
-/*--------------------------------------------------------------------------------------*/
-/* Tcholetsky matrix solution */
-
-void tcholSolve(double **N, double *TN, double *parVect, int n, int BW)
-{
-
-    double **T;
-    int i, j, start, end;
-
-	/*--------------------------------------*/
-    T = G_alloc_matrix(n, BW);
-
-	/*--------------------------------------*/
-    tcholDec(N, T, n, BW);	/* T computation                */
-
-    /* Forward substitution */
-    parVect[0] = TN[0] / T[0][0];
-    for (i = 1; i < n; i++) {
-	parVect[i] = TN[i];
-	/* start = 0 or i - BW + 1 */
-	start = ((i - BW + 1) < 0 ? 0 : (i - BW + 1));
-	/* end = i */
-	for (j = start; j < i; j++)
-	    parVect[i] -= T[j][i - j] * parVect[j];
-	parVect[i] = parVect[i] / T[i][0];
-    }
-
-    /* Backward substitution */
-    parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
-    for (i = n - 2; i >= 0; i--) {
-	/* start = i + 1 */
-	/* end = n or BW + i */
-	end = (n < (BW + i) ? n : (BW + i));
-	for (j = i + 1; j < end; j++)
-	    parVect[i] -= T[i][j - i] * parVect[j];
-	parVect[i] = parVect[i] / T[i][0];
-    }
-
-	/*--------------------------------------*/
-    G_free_matrix(T);
-
-    return;
-}
-
-
-/*--------------------------------------------------------------------------------------*/
-/* Soluzione con Tcholetsky -> la matrice T triangolare viene passata come paramtero e 
-   non calcolata internamente alla procedura -> T = dmatrix (0, n-1, 0, BW-1) */
-
-void tcholSolve2(double **N, double *TN, double **T, double *parVect, int n,
-		 int BW)
-{
-
-    int i, j, start, end;
-
-    /* Forward substitution */
-    parVect[0] = TN[0] / T[0][0];
-    for (i = 1; i < n; i++) {
-	parVect[i] = TN[i];
-	/* start = 0 or i - BW + 1 */
-	start = ((i - BW + 1) < 0 ? 0 : (i - BW + 1));
-	/* end = i */
-	for (j = start; j < i; j++)
-	    parVect[i] -= T[j][i - j] * parVect[j];
-	parVect[i] = parVect[i] / T[i][0];
-    }
-
-    /* Backward substitution */
-    parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
-    for (i = n - 2; i >= 0; i--) {
-	/* start = i + 1 */
-	/* end = n or BW + i */
-	end = (n < (BW + i) ? n : (BW + i));
-	for (j = i + 1; j < end; j++)
-	    parVect[i] -= T[i][j - i] * parVect[j];
-	parVect[i] = parVect[i] / T[i][0];
-    }
-
-    return;
-}
-
-/*--------------------------------------------------------------------------------------*/
-/* Tcholetsky matrix invertion */
-
-void tcholInv(double **N, double *invNdiag, int n, int BW)
-{
-    double **T = NULL;
-    double *vect = NULL;
-    int i, j, k, start;
-    double somma;
-
-	/*--------------------------------------*/
-    T = G_alloc_matrix(n, BW);
-    vect = G_alloc_vector(n);
-
-    /* T computation                */
-    tcholDec(N, T, n, BW);
-
-    /* T Diagonal invertion */
-    for (i = 0; i < n; i++) {
-	T[i][0] = 1.0 / T[i][0];
-    }
-
-    /* N Diagonal invertion */
-    for (i = 0; i < n; i++) {
-	vect[0] = T[i][0];
-	invNdiag[i] = vect[0] * vect[0];
-	for (j = i + 1; j < n; j++) {
-	    somma = 0.0;
-	    /* start = i or j - BW + 1 */
-	    start = ((j - BW + 1) < i ? i : (j - BW + 1));
-	    /* end = j */
-	    for (k = start; k < j; k++) {
-		somma -= vect[k - i] * T[k][j - k];
-	    }
-	    vect[j - i] = somma * T[j][0];
-	    invNdiag[i] += vect[j - i] * vect[j - i];
-	}
-    }
-
-	/*--------------------------------------*/
-    G_free_matrix(T);
-    G_free_vector(vect);
-
-    return;
-}
-
-/*--------------------------------------------------------------------------------------*/
-/* Tcholetsky matrix solution and invertion */
-
-void tcholSolveInv(double **N, double *TN, double *invNdiag, double *parVect,
-		   int n, int BW)
-{
-
-    double **T = NULL;
-    double *vect = NULL;
-    int i, j, k, start, end;
-    double somma;
-
-	/*--------------------------------------*/
-    T = G_alloc_matrix(n, BW);
-    vect = G_alloc_vector(n);
-
-    /* T computation                */
-    tcholDec(N, T, n, BW);
-
-    /* Forward substitution */
-    parVect[0] = TN[0] / T[0][0];
-    for (i = 1; i < n; i++) {
-	parVect[i] = TN[i];
-	/* start = 0 or i - BW + 1 */
-	start = ((i - BW + 1) < 0 ? 0 : (i - BW + 1));
-	/* end = i */
-	for (j = start; j < i; j++)
-	    parVect[i] -= T[j][i - j] * parVect[j];
-	parVect[i] = parVect[i] / T[i][0];
-    }
-
-    /* Backward substitution */
-    parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
-    for (i = n - 2; i >= 0; i--) {
-	/* start = i + 1 */
-	/* end = n or BW + i */
-	end = (n < (BW + i) ? n : (BW + i));
-	for (j = i + 1; j < end; j++)
-	    parVect[i] -= T[i][j - i] * parVect[j];
-	parVect[i] = parVect[i] / T[i][0];
-    }
-
-    /* T Diagonal invertion */
-    for (i = 0; i < n; i++) {
-	T[i][0] = 1.0 / T[i][0];
-    }
-
-    /* N Diagonal invertion */
-    for (i = 0; i < n; i++) {
-	vect[0] = T[i][0];
-	invNdiag[i] = vect[0] * vect[0];
-	for (j = i + 1; j < n; j++) {
-	    somma = 0.0;
-	    /* start = i or j - BW + 1 */
-	    start = ((j - BW + 1) < i ? i : (j - BW + 1));
-	    /* end = j */
-	    for (k = start; k < j; k++) {
-		somma -= vect[k - i] * T[k][j - k];
-	    }
-	    vect[j - i] = somma * T[j][0];
-	    invNdiag[i] += vect[j - i] * vect[j - i];
-	}
-    }
-
-	/*--------------------------------------*/
-    G_free_matrix(T);
-    G_free_vector(vect);
-
-    return;
-}

vector/lidar/lidarlib/zones.c

@@ -14,11 +14,11 @@
 /*----------------------------------------------------------------------------------------*/
 void P_zero_dim(struct Reg_dimens *dim)
 {
-    dim->orlo_h = 0.0;
-    dim->orlo_v = 0.0;
+    dim->edge_h = 0.0;
+    dim->edge_v = 0.0;
     dim->overlap = 0.0;
-    dim->latoN = 0.0;
-    dim->latoE = 0.0;
+    dim->sn_size = 0.0;
+    dim->ew_size = 0.0;
     return;
 }
 
@@ -69,51 +69,51 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
     switch (type) {
     case GENERAL_ROW:		/* General case N-S direction */
 	Elaboration->north =
-	    Elaboration->south + dim.overlap + (2 * dim.orlo_h);
-	Elaboration->south = Elaboration->north - dim.latoN;
-	General->N = Elaboration->north - dim.orlo_h;
-	General->S = Elaboration->south + dim.orlo_h;
+	    Elaboration->south + dim.overlap + (2 * dim.edge_h);
+	Elaboration->south = Elaboration->north - dim.sn_size;
+	General->N = Elaboration->north - dim.edge_h;
+	General->S = Elaboration->south + dim.edge_h;
 	Overlap->N = General->N - dim.overlap;
 	Overlap->S = General->S + dim.overlap;
 	return 0;
 
     case GENERAL_COLUMN:	/* General case E-W direction */
 	Elaboration->west =
-	    Elaboration->east - dim.overlap - (2 * dim.orlo_v);
-	Elaboration->east = Elaboration->west + dim.latoE;
-	General->W = Elaboration->west + dim.orlo_v;
-	General->E = Elaboration->east - dim.orlo_v;
+	    Elaboration->east - dim.overlap - (2 * dim.edge_v);
+	Elaboration->east = Elaboration->west + dim.ew_size;
+	General->W = Elaboration->west + dim.edge_v;
+	General->E = Elaboration->east - dim.edge_v;
 	Overlap->W = General->W + dim.overlap;
 	Overlap->E = General->E - dim.overlap;
 	return 0;
 
     case FIRST_ROW:		/* Just started with first row */
-	Elaboration->north = orig.north + 2 * dim.orlo_h;
-	Elaboration->south = Elaboration->north - dim.latoN;
-	General->N = Elaboration->north - 2 * dim.orlo_h;
-	General->S = Elaboration->south + dim.orlo_h;
+	Elaboration->north = orig.north + 2 * dim.edge_h;
+	Elaboration->south = Elaboration->north - dim.sn_size;
+	General->N = Elaboration->north - 2 * dim.edge_h;
+	General->S = Elaboration->south + dim.edge_h;
 	Overlap->N = General->N;
 	Overlap->S = General->S + dim.overlap;
 	return 0;
 
     case LAST_ROW:		/* Reached last row */
-	Elaboration->south = orig.south - 2 * dim.orlo_h;
-	General->S = Elaboration->south + 2 * dim.orlo_h;
+	Elaboration->south = orig.south - 2 * dim.edge_h;
+	General->S = Elaboration->south + 2 * dim.edge_h;
 	Overlap->S = General->S;
 	return 0;
 
     case FIRST_COLUMN:		/* Just started with first column */
-	Elaboration->west = orig.west - 2 * dim.orlo_v;
-	Elaboration->east = Elaboration->west + dim.latoE;
-	General->W = Elaboration->west + 2 * dim.orlo_v;
-	General->E = Elaboration->east - dim.orlo_v;
+	Elaboration->west = orig.west - 2 * dim.edge_v;
+	Elaboration->east = Elaboration->west + dim.ew_size;
+	General->W = Elaboration->west + 2 * dim.edge_v;
+	General->E = Elaboration->east - dim.edge_v;
 	Overlap->W = General->W;
 	Overlap->E = General->E - dim.overlap;
 	return 0;
 
     case LAST_COLUMN:		/* Reached last column */
-	Elaboration->east = orig.east + 2 * dim.orlo_v;
-	General->E = Elaboration->east - 2 * dim.orlo_v;
+	Elaboration->east = orig.east + 2 * dim.edge_v;
+	General->E = Elaboration->east - 2 * dim.edge_v;
 	Overlap->E = General->E;
 	return 0;
     }
@@ -124,8 +124,8 @@ P_set_regions(struct Cell_head *Elaboration, struct bound_box * General,
 /*----------------------------------------------------------------------------------------*/
 int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsply)
 {
-    int total_splines, orlo_splines, n_windows, minsplines;
-    double E_extension, N_extension, orloE, orloN;
+    int total_splines, edge_splines, n_windows, minsplines;
+    double E_extension, N_extension, edgeE, edgeN;
     struct Cell_head orig;
     int ret = 0;
 
@@ -133,67 +133,67 @@ int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsp
 
     E_extension = orig.east - orig.west;
     N_extension = orig.north - orig.south;
-    dim->latoE = *nsplx * pe;
-    dim->latoN = *nsply * pn;
-    orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
-    orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
-
-    /* number of moving windows: E_extension / orloE */
-    /* remaining steps: total steps - (floor(E_extension / orloE) * E_extension) / passoE */
-    /* remaining steps must be larger than orlo_v + overlap + half of overlap window */
+    dim->ew_size = *nsplx * pe;
+    dim->sn_size = *nsply * pn;
+    edgeE = dim->ew_size - dim->overlap - 2 * dim->edge_v;
+    edgeN = dim->sn_size - dim->overlap - 2 * dim->edge_h;
+
+    /* number of moving windows: E_extension / edgeE */
+    /* remaining steps: total steps - (floor(E_extension / edgeE) * E_extension) / passoE */
+    /* remaining steps must be larger than edge_v + overlap + half of overlap window */
     total_splines = ceil(E_extension / pe);
-    orlo_splines = orloE / pe;
-    n_windows = floor(E_extension / orloE); /* without last one */
+    edge_splines = edgeE / pe;
+    n_windows = floor(E_extension / edgeE); /* without last one */
     if (n_windows > 0) {
-	minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
-	while (total_splines - orlo_splines * n_windows < minsplines) {
+	minsplines = ceil((double)(dim->ew_size / 2.0 - dim->edge_v - dim->overlap) / pe);
+	while (total_splines - edge_splines * n_windows < minsplines) {
 	    *nsplx -= 1;
-	    dim->latoE = *nsplx * pe;
-	    orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+	    dim->ew_size = *nsplx * pe;
+	    edgeE = dim->ew_size - dim->overlap - 2 * dim->edge_v;
 
-	    orlo_splines = orloE / pe;
-	    n_windows = floor(E_extension / orloE); /* without last one */
-	    minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+	    edge_splines = edgeE / pe;
+	    n_windows = floor(E_extension / edgeE); /* without last one */
+	    minsplines = ceil((double)(dim->ew_size / 2.0 - dim->edge_v - dim->overlap) / pe);
 	    if (ret == 0)
 		ret = 1;
 	}
-	while (total_splines - orlo_splines * n_windows < minsplines * 2 && minsplines > 30) {
+	while (total_splines - edge_splines * n_windows < minsplines * 2 && minsplines > 30) {
 	    *nsplx -= 1;
-	    dim->latoE = *nsplx * pe;
-	    orloE = dim->latoE - dim->overlap - 2 * dim->orlo_v;
+	    dim->ew_size = *nsplx * pe;
+	    edgeE = dim->ew_size - dim->overlap - 2 * dim->edge_v;
 
-	    orlo_splines = orloE / pe;
-	    n_windows = floor(E_extension / orloE); /* without last one */
-	    minsplines = ceil((double)(dim->latoE / 2.0 - dim->orlo_v - dim->overlap) / pe);
+	    edge_splines = edgeE / pe;
+	    n_windows = floor(E_extension / edgeE); /* without last one */
+	    minsplines = ceil((double)(dim->ew_size / 2.0 - dim->edge_v - dim->overlap) / pe);
 	    if (ret == 0)
 		ret = 1;
 	}
     }
 
     total_splines = ceil(N_extension / pn);
-    orlo_splines = orloN / pn;
-    n_windows = floor(N_extension / orloN); /* without last one */
+    edge_splines = edgeN / pn;
+    n_windows = floor(N_extension / edgeN); /* without last one */
     if (n_windows > 0) {
-	minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
-	while (total_splines - orlo_splines * n_windows < minsplines) {
+	minsplines = ceil((double)(dim->sn_size / 2.0 - dim->edge_h - dim->overlap) / pn);
+	while (total_splines - edge_splines * n_windows < minsplines) {
 	    *nsply -= 1;
-	    dim->latoN = *nsply * pn;
-	    orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+	    dim->sn_size = *nsply * pn;
+	    edgeN = dim->sn_size - dim->overlap - 2 * dim->edge_h;
 
-	    orlo_splines = orloN / pn;
-	    n_windows = floor(N_extension / orloN); /* without last one */
-	    minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+	    edge_splines = edgeN / pn;
+	    n_windows = floor(N_extension / edgeN); /* without last one */
+	    minsplines = ceil((double)(dim->sn_size / 2.0 - dim->edge_h - dim->overlap) / pn);
 	    if (ret < 2)
 		ret += 2;
 	}
-	while (total_splines - orlo_splines * n_windows < minsplines * 2 && minsplines > 30) {
+	while (total_splines - edge_splines * n_windows < minsplines * 2 && minsplines > 30) {
 	    *nsply -= 1;
-	    dim->latoN = *nsply * pn;
-	    orloN = dim->latoN - dim->overlap - 2 * dim->orlo_h;
+	    dim->sn_size = *nsply * pn;
+	    edgeN = dim->sn_size - dim->overlap - 2 * dim->edge_h;
 
-	    orlo_splines = orloN / pn;
-	    n_windows = floor(N_extension / orloN); /* without last one */
-	    minsplines = ceil((double)(dim->latoN / 2.0 - dim->orlo_h - dim->overlap) / pn);
+	    edge_splines = edgeN / pn;
+	    n_windows = floor(N_extension / edgeN); /* without last one */
+	    minsplines = ceil((double)(dim->sn_size / 2.0 - dim->edge_h - dim->overlap) / pn);
 	    if (ret < 2)
 		ret += 2;
 	}
@@ -202,20 +202,20 @@ int P_set_dim(struct Reg_dimens *dim, double pe, double pn, int *nsplx, int *nsp
 }
 
 /*----------------------------------------------------------------------------------------*/
-int P_get_orlo(int interpolator, struct Reg_dimens *dim, double pe, double pn)
+int P_get_edge(int interpolator, struct Reg_dimens *dim, double pe, double pn)
 {
-    /* Set the orlo regions dimension
+    /* Set the edge regions dimension
      * Returns 1 on success of bilinear; 2 on success of bicubic, 0 on failure */
     if (interpolator == P_BILINEAR) {
        	/* in case of edge artefacts, increase as multiples of 3 */
-	dim->orlo_v = 9 * pe;
-	dim->orlo_h = 9 * pn;
+	dim->edge_v = 9 * pe;
+	dim->edge_h = 9 * pn;
 	return 1;
     }
     else if (interpolator == P_BICUBIC) {
        	/* in case of edge artefacts, increase as multiples of 4 */
-	dim->orlo_v = 12 * pe;	/*3 */
-	dim->orlo_h = 12 * pn;
+	dim->edge_v = 12 * pe;	/*3 */
+	dim->edge_h = 12 * pn;
 	return 2;
     }
     else
@@ -571,27 +571,6 @@ P_Aux_to_Vector(struct Map_info *Map, struct Map_info *Out, dbDriver * driver,
     return;
 }
 
-/*------------------------------------------------------------------------------------------------*/
-#ifdef notdef
-double **P_Null_Matrix(double **matrix)
-{
-    int nrows, row, ncols, col;
-    struct Cell_head Original;
-
-
-    G_get_window(&Original);
-    Rast_set_window(&Original);
-    nrows = Rast_window_rows();
-    ncols = Rast_window_cols();
-
-    for (row = 0; row < nrows; row++) {
-	for (col = 0; col < ncols; col++) {
-	    matrix[row][col] = NULL;
-	}
-    }
-}
-#endif
-
 /*! DEFINITION OF THE SUBZONES
 
   5: inside Overlap region

vector/lidar/v.lidar.correction/correction.c

@@ -32,7 +32,7 @@ P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 		    struct Map_info *Terrain, struct Cell_head *Elaboration,
 		    struct bound_box General, struct bound_box Overlap,
 		    double **obs, struct lidar_cat *lcat, double *param,
-		    int *line_num, double passoN, double passoE,
+		    int *line_num, double stepN, double stepE,
 		    double overlap, double HighThresh, double LowThresh,
 		    int nsplx, int nsply, int num_points,
 		    dbDriver * driver, double mean, char *tab_name)
@@ -58,7 +58,7 @@ P_Sparse_Correction(struct Map_info *In, struct Map_info *Out,
 	
 	if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {
 	    interpolation =
-		dataInterpolateBilin(obs[i][0], obs[i][1], passoE, passoN,
+		dataInterpolateBilin(obs[i][0], obs[i][1], stepE, stepN,
 				     nsplx, nsply, Elaboration->west,
 				     Elaboration->south, param);
 	    interpolation += mean;

vector/lidar/v.lidar.correction/main.c

@@ -41,9 +41,9 @@ int main(int argc, char *argv[])
     const char *dvr, *db, *mapset;
     char table_name[GNAME_MAX];
     char xname[GNAME_MAX], xmapset[GMAPSET_MAX];
-    double lambda, ew_resol, ns_resol, mean, passoN, passoE, HighThresh,
+    double lambda, ew_resol, ns_resol, mean, stepN, stepE, HighThresh,
 	LowThresh;
-    double N_extension, E_extension, orloE, orloN;
+    double N_extension, E_extension, edgeE, edgeN;
 
     int i, nterrain, count_terrain;
 
@@ -54,8 +54,8 @@ int main(int argc, char *argv[])
     double **N, **obsVect, **obsVect_all;	/* Interpolation and least-square matrix */
 
     struct Map_info In, Out, Terrain;
-    struct Option *in_opt, *out_opt, *out_terrain_opt, *passoE_opt,
-	*passoN_opt, *lambda_f_opt, *Thresh_A_opt, *Thresh_B_opt;
+    struct Option *in_opt, *out_opt, *out_terrain_opt, *stepE_opt,
+	*stepN_opt, *lambda_f_opt, *Thresh_A_opt, *Thresh_B_opt;
     struct Flag *spline_step_flag;
     struct GModule *module;
 
@@ -98,20 +98,20 @@ int main(int argc, char *argv[])
     out_terrain_opt->description =
 	_("Only 'terrain' points output vector map");
 
-    passoE_opt = G_define_option();
-    passoE_opt->key = "sce";
-    passoE_opt->type = TYPE_DOUBLE;
-    passoE_opt->required = NO;
-    passoE_opt->answer = "25";
-    passoE_opt->description =
+    stepE_opt = G_define_option();
+    stepE_opt->key = "sce";
+    stepE_opt->type = TYPE_DOUBLE;
+    stepE_opt->required = NO;
+    stepE_opt->answer = "25";
+    stepE_opt->description =
 	_("Interpolation spline step value in east direction");
 
-    passoN_opt = G_define_option();
-    passoN_opt->key = "scn";
-    passoN_opt->type = TYPE_DOUBLE;
-    passoN_opt->required = NO;
-    passoN_opt->answer = "25";
-    passoN_opt->description =
+    stepN_opt = G_define_option();
+    stepN_opt->key = "scn";
+    stepN_opt->type = TYPE_DOUBLE;
+    stepN_opt->required = NO;
+    stepN_opt->answer = "25";
+    stepN_opt->description =
 	_("Interpolation spline step value in north direction");
 
     lambda_f_opt = G_define_option();
@@ -144,8 +144,8 @@ int main(int argc, char *argv[])
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-    passoN = atof(passoN_opt->answer);
-    passoE = atof(passoE_opt->answer);
+    stepN = atof(stepN_opt->answer);
+    stepE = atof(stepE_opt->answer);
     lambda = atof(lambda_f_opt->answer);
     HighThresh = atof(Thresh_A_opt->answer);
     LowThresh = atof(Thresh_B_opt->answer);
@@ -262,7 +262,7 @@ int main(int argc, char *argv[])
       | Each original region will be divided into several subregions. 
       | Each one will be overlaped by its neighbouring subregions. 
       | The overlapping is calculated as a fixed OVERLAP_SIZE times
-      | the largest spline step plus 2 * orlo
+      | the largest spline step plus 2 * edge
       ----------------------------------------------------------------*/
 
     /* Fixing parameters of the elaboration region */
@@ -270,25 +270,25 @@ int main(int argc, char *argv[])
 
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    if (passoN > passoE)
-	dims.overlap = OVERLAP_SIZE * passoN;
+    if (stepN > stepE)
+	dims.overlap = OVERLAP_SIZE * stepN;
     else
-	dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
-    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+	dims.overlap = OVERLAP_SIZE * stepE;
+    P_get_edge(P_BILINEAR, &dims, stepE, stepN);
+    P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
     G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
 
     /* calculate number of subregions */
-    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
-    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+    edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
+    edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
 
     N_extension = original_reg.north - original_reg.south;
     E_extension = original_reg.east - original_reg.west;
 
-    nsubregion_col = ceil(E_extension / orloE) + 0.5;
-    nsubregion_row = ceil(N_extension / orloN) + 0.5;
+    nsubregion_col = ceil(E_extension / edgeE) + 0.5;
+    nsubregion_row = ceil(N_extension / edgeN) + 0.5;
 
     if (nsubregion_col < 0)
 	nsubregion_col = 0;
@@ -318,7 +318,7 @@ int main(int argc, char *argv[])
 
 	nsply =
 	    ceil((elaboration_reg.north -
-		  elaboration_reg.south) / passoN) + 0.5;
+		  elaboration_reg.south) / stepN) + 0.5;
 	/*
 	if (nsply > NSPLY_MAX) {
 	    nsply = NSPLY_MAX;
@@ -350,7 +350,7 @@ int main(int argc, char *argv[])
 	    }
 
 	    nsplx =
-		ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
+		ceil((elaboration_reg.east - elaboration_reg.west) / stepE) +
 		0.5;
 	    /*
 	    if (nsplx > NSPLX_MAX) {
@@ -403,12 +403,12 @@ int main(int argc, char *argv[])
 		G_free(observ);
 
 		G_verbose_message(_("Bilinear interpolation"));
-		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
+		normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
 			       nsply, elaboration_reg.west,
 			       elaboration_reg.south, nterrain, nparameters,
 			       BW);
-		nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
-		tcholSolve(N, TN, parVect, nparameters, BW);
+		nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
+		G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
 
 		G_free_matrix(N);
 		G_free_vector(TN);
@@ -418,7 +418,7 @@ int main(int argc, char *argv[])
 		G_verbose_message( _("Correction and creation of terrain vector"));
 		P_Sparse_Correction(&In, &Out, &Terrain, &elaboration_reg,
 				    general_box, overlap_box, obsVect_all, lcat,
-				    parVect, lineVect, passoN, passoE,
+				    parVect, lineVect, stepN, stepE,
 				    dims.overlap, HighThresh, LowThresh,
 				    nsplx, nsply, npoints, driver, mean, table_name);
 

vector/lidar/v.lidar.edgedetection/edgedetection.c

@@ -45,7 +45,7 @@ int edge_detection(struct Cell_head elaboration_reg, struct bound_box Overlap_Bo
 
     int c1, c2;
     double g[9][2], gradient[2], gradPto, dirPto;
-    extern double passoE, passoN;
+    extern double stepE, stepN;
     static struct Cell_head Elaboration;
 
     g[0][0] = partial[0];
@@ -62,13 +62,13 @@ int edge_detection(struct Cell_head elaboration_reg, struct bound_box Overlap_Bo
     else if ((gradPto > gradLow) && (residual > 0)) {	/* Soft condition for 'edge' points */
 
 	if (Vect_point_in_box(obsX, obsY, 0.0, &Overlap_Box)) {
-	    Get_Gradient(Elaboration, obsX + passoE * cos(dirPto),
-			     obsY + passoN * sin(dirPto), parBilin, gradient);
+	    Get_Gradient(Elaboration, obsX + stepE * cos(dirPto),
+			     obsY + stepN * sin(dirPto), parBilin, gradient);
 	    g[2][0] = gradient[0];
 	    g[2][1] = gradient[1];
 
-	    Get_Gradient(Elaboration, obsX + passoE * cos(dirPto + M_PI),
-			     obsY + passoN * sin(dirPto + M_PI), parBilin, gradient);
+	    Get_Gradient(Elaboration, obsX + stepE * cos(dirPto + M_PI),
+			     obsY + stepN * sin(dirPto + M_PI), parBilin, gradient);
 	    g[7][0] = gradient[0];
 	    g[7][1] = gradient[1];
 
@@ -76,43 +76,43 @@ int edge_detection(struct Cell_head elaboration_reg, struct bound_box Overlap_Bo
 		(fabs(atan(g[7][1] / g[7][0]) + M_PI / 2 - dirPto) < alpha)) {
 
 		Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto + M_PI / 4),
-				 obsY + passoN * sin(dirPto + M_PI / 4),
+				 obsX + stepE * cos(dirPto + M_PI / 4),
+				 obsY + stepN * sin(dirPto + M_PI / 4),
 				 parBilin, gradient);
 		g[1][0] = gradient[0];
 		g[1][1] = gradient[1];
 
 		Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto - M_PI / 4),
-				 obsY + passoN * sin(dirPto - M_PI / 4),
+				 obsX + stepE * cos(dirPto - M_PI / 4),
+				 obsY + stepN * sin(dirPto - M_PI / 4),
 				 parBilin, gradient);
 		g[3][0] = gradient[0];
 		g[3][1] = gradient[1];
 
 		Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto + M_PI / 2),
-				 obsY + passoN * sin(dirPto + M_PI / 2),
+				 obsX + stepE * cos(dirPto + M_PI / 2),
+				 obsY + stepN * sin(dirPto + M_PI / 2),
 				 parBilin, gradient);
 		g[4][0] = gradient[0];
 		g[4][1] = gradient[1];
 
 		Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto - M_PI / 2),
-				 obsY + passoN * sin(dirPto - M_PI / 2),
+				 obsX + stepE * cos(dirPto - M_PI / 2),
+				 obsY + stepN * sin(dirPto - M_PI / 2),
 				 parBilin, gradient);
 		g[5][0] = gradient[0];
 		g[5][1] = gradient[1];
 
 		Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto + M_PI * 3 / 4),
-				 obsY + passoN * sin(dirPto + M_PI * 3 / 4),
+				 obsX + stepE * cos(dirPto + M_PI * 3 / 4),
+				 obsY + stepN * sin(dirPto + M_PI * 3 / 4),
 				 parBilin, gradient);
 		g[6][0] = gradient[0];
 		g[6][1] = gradient[1];
 
 		Get_Gradient(Elaboration,
-				 obsX + passoE * cos(dirPto - M_PI * 3 / 4),
-				 obsY + passoN * sin(dirPto - M_PI * 3 / 4),
+				 obsX + stepE * cos(dirPto - M_PI * 3 / 4),
+				 obsY + stepN * sin(dirPto - M_PI * 3 / 4),
 				 parBilin, gradient);
 		g[8][0] = gradient[0];
 		g[8][1] = gradient[1];
@@ -145,13 +145,13 @@ int Get_Gradient(struct Cell_head Elaboration, double X, double Y,
     double csi, eta, d, b, a, c;
 
     extern int nsply;
-    extern double passoN, passoE;
+    extern double stepN, stepE;
 
-    row = (int)((Y - Elaboration.south) / passoN);
-    col = (int)((X - Elaboration.west) / passoE);
+    row = (int)((Y - Elaboration.south) / stepN);
+    col = (int)((X - Elaboration.west) / stepE);
     N = nsply * col + row;
-    eta = X - (Elaboration.west + (col * passoE));
-    csi = Y - (Elaboration.south + (row * passoN));
+    eta = X - (Elaboration.west + (col * stepE));
+    csi = Y - (Elaboration.south + (row * stepN));
     d = parVect[N];
     b = parVect[N + 1] - d;
     a = parVect[N + nsply] - d;
@@ -172,7 +172,7 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
     double interpolation, weight, residual, eta, csi, gradient[2];
 
     extern int nsplx, nsply, line_out_counter;
-    extern double passoN, passoE;
+    extern double stepN, stepE;
 
     struct line_pnts *point;
     struct line_cats *categories;
@@ -190,7 +190,7 @@ void classification(struct Map_info *Out, struct Cell_head Elaboration,
 
 	if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {
 	    interpolation =
-		dataInterpolateBicubic(obs[i][0], obs[i][1], passoE, passoN,
+		dataInterpolateBicubic(obs[i][0], obs[i][1], stepE, stepN,
 				       nsplx, nsply, Elaboration.west,
 				       Elaboration.south, parBicub);
 	    interpolation += mean;

vector/lidar/v.lidar.edgedetection/main.c

@@ -35,7 +35,7 @@
 #include "edgedetection.h"
 
 int nsply, nsplx, line_out_counter;
-double passoN, passoE;
+double stepN, stepE;
 
 /**************************************************************************************
 **************************************************************************************/
@@ -44,7 +44,7 @@ int main(int argc, char *argv[])
     /* Variables' declarations */
     int nsplx_adj, nsply_adj;
     int nsubregion_col, nsubregion_row, subregion = 0, nsubregions = 0;
-    double N_extension, E_extension, orloE, orloN;
+    double N_extension, E_extension, edgeE, edgeN;
     int dim_vect, nparameters, BW, npoints;
     double lambda_B, lambda_F, grad_H, grad_L, alpha, mean;
     const char *dvr, *db, *mapset;
@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
 
     /* Structs' declarations */
     struct Map_info In, Out;
-    struct Option *in_opt, *out_opt, *passoE_opt, *passoN_opt,
+    struct Option *in_opt, *out_opt, *stepE_opt, *stepN_opt,
 	*lambdaF_opt, *lambdaB_opt, *gradH_opt, *gradL_opt, *alfa_opt;
     struct Flag *spline_step_flag;
     struct GModule *module;
@@ -91,23 +91,23 @@ int main(int argc, char *argv[])
 
     out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
 
-    passoE_opt = G_define_option();
-    passoE_opt->key = "see";
-    passoE_opt->type = TYPE_DOUBLE;
-    passoE_opt->required = NO;
-    passoE_opt->answer = "4";
-    passoE_opt->description =
+    stepE_opt = G_define_option();
+    stepE_opt->key = "see";
+    stepE_opt->type = TYPE_DOUBLE;
+    stepE_opt->required = NO;
+    stepE_opt->answer = "4";
+    stepE_opt->description =
 	_("Interpolation spline step value in east direction");
-    passoE_opt->guisection = _("Settings");
-
-    passoN_opt = G_define_option();
-    passoN_opt->key = "sen";
-    passoN_opt->type = TYPE_DOUBLE;
-    passoN_opt->required = NO;
-    passoN_opt->answer = "4";
-    passoN_opt->description =
+    stepE_opt->guisection = _("Settings");
+
+    stepN_opt = G_define_option();
+    stepN_opt->key = "sen";
+    stepN_opt->type = TYPE_DOUBLE;
+    stepN_opt->required = NO;
+    stepN_opt->answer = "4";
+    stepN_opt->description =
 	_("Interpolation spline step value in north direction");
-    passoN_opt->guisection = _("Settings");
+    stepN_opt->guisection = _("Settings");
 
     lambdaB_opt = G_define_option();
     lambdaB_opt->key = "lambda_g";
@@ -160,8 +160,8 @@ int main(int argc, char *argv[])
 	exit(EXIT_FAILURE);
 
     line_out_counter = 1;
-    passoN = atof(passoN_opt->answer);
-    passoE = atof(passoE_opt->answer);
+    stepN = atof(stepN_opt->answer);
+    stepE = atof(stepE_opt->answer);
     lambda_F = atof(lambdaF_opt->answer);
     lambda_B = atof(lambdaB_opt->answer);
     grad_H = atof(gradH_opt->answer);
@@ -282,7 +282,7 @@ int main(int argc, char *argv[])
       | Each original region will be divided into several subregions. 
       | Each one will be overlaped by its neighbouring subregions. 
       | The overlapping is calculated as a fixed OVERLAP_SIZE times
-      | the largest spline step plus 2 * orlo
+      | the largest spline step plus 2 * edge
       ----------------------------------------------------------------*/
 
     /* Fixing parameters of the elaboration region */
@@ -290,25 +290,25 @@ int main(int argc, char *argv[])
 
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    if (passoN > passoE)
-	dims.overlap = OVERLAP_SIZE * passoN;
+    if (stepN > stepE)
+	dims.overlap = OVERLAP_SIZE * stepN;
     else
-	dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(P_BICUBIC, &dims, passoE, passoN);
-    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+	dims.overlap = OVERLAP_SIZE * stepE;
+    P_get_edge(P_BICUBIC, &dims, stepE, stepN);
+    P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
     G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
 
     /* calculate number of subregions */
-    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
-    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+    edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
+    edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
 
     N_extension = original_reg.north - original_reg.south;
     E_extension = original_reg.east - original_reg.west;
 
-    nsubregion_col = ceil(E_extension / orloE) + 0.5;
-    nsubregion_row = ceil(N_extension / orloN) + 0.5;
+    nsubregion_col = ceil(E_extension / edgeE) + 0.5;
+    nsubregion_row = ceil(N_extension / edgeN) + 0.5;
 
     if (nsubregion_col < 0)
 	nsubregion_col = 0;
@@ -337,7 +337,7 @@ int main(int argc, char *argv[])
 	}
 
 	nsply =
-	    ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
+	    ceil((elaboration_reg.north - elaboration_reg.south) / stepN) +
 	    0.5;
 	/*
 	if (nsply > NSPLY_MAX) {
@@ -370,7 +370,7 @@ int main(int argc, char *argv[])
 	    }
 
 	    nsplx =
-		ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
+		ceil((elaboration_reg.east - elaboration_reg.west) / stepE) +
 		0.5;
 	    /*
 	    if (nsplx > NSPLX_MAX) {
@@ -417,12 +417,12 @@ int main(int argc, char *argv[])
 		G_free(observ);
 
 		G_verbose_message(_("Bilinear interpolation"));
-		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
+		normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
 			       nsply, elaboration_reg.west,
 			       elaboration_reg.south, npoints, nparameters,
 			       BW);
-		nCorrectGrad(N, lambda_B, nsplx, nsply, passoE, passoN);
-		tcholSolve(N, TN, parVect_bilin, nparameters, BW);
+		nCorrectGrad(N, lambda_B, nsplx, nsply, stepE, stepN);
+		G_math_solver_cholesky_sband(N, parVect_bilin, TN, nparameters, BW);
 
 		G_free_matrix(N);
 		for (tn = 0; tn < nparameters; tn++)
@@ -434,12 +434,12 @@ int main(int argc, char *argv[])
 		parVect_bicub = G_alloc_vector(nparameters);	/* Bicubic parameters vector */
 
 		G_verbose_message(_("Bicubic interpolation"));
-		normalDefBicubic(N, TN, Q, obsVect, passoE, passoN, nsplx,
+		normalDefBicubic(N, TN, Q, obsVect, stepE, stepN, nsplx,
 				 nsply, elaboration_reg.west,
 				 elaboration_reg.south, npoints, nparameters,
 				 BW);
-		nCorrectLapl(N, lambda_F, nsplx, nsply, passoE, passoN);
-		tcholSolve(N, TN, parVect_bicub, nparameters, BW);
+		nCorrectLapl(N, lambda_F, nsplx, nsply, stepE, stepN);
+		G_math_solver_cholesky_sband(N, parVect_bicub, TN, nparameters, BW);
 
 		G_free_matrix(N);
 		G_free_vector(TN);

vector/lidar/v.lidar.growing/main.c

@@ -32,7 +32,7 @@
 #include "growing.h"
     /* GLOBAL DEFINITIONS */
 int nsply, nsplx, count_obj;
-double passoN, passoE;
+double stepN, stepE;
 
 /*--------------------------------------------------------------------------------*/
 int main(int argc, char *argv[])

vector/lidar/v.outlier/main.c

@@ -30,7 +30,7 @@
 #include "outlier.h"
     /* GLOBAL VARIABLES DEFINITIONS */
 int nsply, nsplx;
-double passoN, passoE, Thres_Outlier;
+double stepN, stepE, Thres_Outlier;
 
 /*--------------------------------------------------------------------------------------*/
 int main(int argc, char *argv[])
@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
     int nsplx_adj, nsply_adj;
     int nsubregion_col, nsubregion_row;
     int subregion = 0, nsubregions = 0;
-    double N_extension, E_extension, orloE, orloN;
+    double N_extension, E_extension, edgeE, edgeN;
     int dim_vect, nparameters, BW, npoints;
     double mean, lambda;
     const char *dvr, *db, *mapset;
@@ -54,8 +54,8 @@ int main(int argc, char *argv[])
 
     /* Structs' declarations */
     struct Map_info In, Out, Outlier, Qgis;
-    struct Option *in_opt, *out_opt, *outlier_opt, *qgis_opt, *passoE_opt,
-	*passoN_opt, *lambda_f_opt, *Thres_O_opt;
+    struct Option *in_opt, *out_opt, *outlier_opt, *qgis_opt, *stepE_opt,
+	*stepN_opt, *lambda_f_opt, *Thres_O_opt;
     struct Flag *spline_step_flag;
     struct GModule *module;
 
@@ -100,20 +100,20 @@ int main(int argc, char *argv[])
     qgis_opt->gisprompt = "new,vector,vector";
     qgis_opt->description = _("Name of vector map for visualization in QGIS");
 
-    passoE_opt = G_define_option();
-    passoE_opt->key = "soe";
-    passoE_opt->type = TYPE_DOUBLE;
-    passoE_opt->required = NO;
-    passoE_opt->answer = "10";
-    passoE_opt->description =
+    stepE_opt = G_define_option();
+    stepE_opt->key = "soe";
+    stepE_opt->type = TYPE_DOUBLE;
+    stepE_opt->required = NO;
+    stepE_opt->answer = "10";
+    stepE_opt->description =
 	_("Interpolation spline step value in east direction");
 
-    passoN_opt = G_define_option();
-    passoN_opt->key = "son";
-    passoN_opt->type = TYPE_DOUBLE;
-    passoN_opt->required = NO;
-    passoN_opt->answer = "10";
-    passoN_opt->description =
+    stepN_opt = G_define_option();
+    stepN_opt->key = "son";
+    stepN_opt->type = TYPE_DOUBLE;
+    stepN_opt->required = NO;
+    stepN_opt->answer = "10";
+    stepN_opt->description =
 	_("Interpolation spline step value in north direction");
 
     lambda_f_opt = G_define_option();
@@ -142,8 +142,8 @@ int main(int argc, char *argv[])
     if (!(dvr = G__getenv2("DB_DRIVER", G_VAR_MAPSET)))
 	G_fatal_error(_("Unable to read name of driver"));
 
-    passoN = atof(passoN_opt->answer);
-    passoE = atof(passoE_opt->answer);
+    stepN = atof(stepN_opt->answer);
+    stepE = atof(stepE_opt->answer);
     lambda = atof(lambda_f_opt->answer);
     Thres_Outlier = atof(Thres_O_opt->answer);
 
@@ -259,7 +259,7 @@ int main(int argc, char *argv[])
       | Each original region will be divided into several subregions. 
       | Each one will be overlaped by its neighbouring subregions. 
       | The overlapping is calculated as a fixed OVERLAP_SIZE times
-      | the largest spline step plus 2 * orlo
+      | the largest spline step plus 2 * edge
       ----------------------------------------------------------------*/
 
     /* Fixing parameters of the elaboration region */
@@ -267,25 +267,25 @@ int main(int argc, char *argv[])
 
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    if (passoN > passoE)
-	dims.overlap = OVERLAP_SIZE * passoN;
+    if (stepN > stepE)
+	dims.overlap = OVERLAP_SIZE * stepN;
     else
-	dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(P_BILINEAR, &dims, passoE, passoN);
-    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+	dims.overlap = OVERLAP_SIZE * stepE;
+    P_get_edge(P_BILINEAR, &dims, stepE, stepN);
+    P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
     G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
 
     /* calculate number of subregions */
-    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
-    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+    edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
+    edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
 
     N_extension = original_reg.north - original_reg.south;
     E_extension = original_reg.east - original_reg.west;
 
-    nsubregion_col = ceil(E_extension / orloE) + 0.5;
-    nsubregion_row = ceil(N_extension / orloN) + 0.5;
+    nsubregion_col = ceil(E_extension / edgeE) + 0.5;
+    nsubregion_row = ceil(N_extension / edgeN) + 0.5;
 
     if (nsubregion_col < 0)
 	nsubregion_col = 0;
@@ -314,7 +314,7 @@ int main(int argc, char *argv[])
 	}
 
 	nsply =
-	    ceil((elaboration_reg.north - elaboration_reg.south) / passoN) +
+	    ceil((elaboration_reg.north - elaboration_reg.south) / stepN) +
 	    0.5;
 	/*
 	if (nsply > NSPLY_MAX) {
@@ -347,7 +347,7 @@ int main(int argc, char *argv[])
 	    }
 
 	    nsplx =
-		ceil((elaboration_reg.east - elaboration_reg.west) / passoE) +
+		ceil((elaboration_reg.east - elaboration_reg.west) / stepE) +
 		0.5;
 	    /*
 	    if (nsplx > NSPLX_MAX) {
@@ -392,12 +392,12 @@ int main(int argc, char *argv[])
 		G_free(observ);
 
 		G_verbose_message(_("Bilinear interpolation"));
-		normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
+		normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
 			       nsply, elaboration_reg.west,
 			       elaboration_reg.south, npoints, nparameters,
 			       BW);
-		nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
-		tcholSolve(N, TN, parVect, nparameters, BW);
+		nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
+		G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
 
 		G_free_matrix(N);
 		G_free_vector(TN);

vector/lidar/v.outlier/outlier.c

@@ -20,7 +20,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
     int i;
     double interpolation, weight, residual, eta, csi;
     extern int nsplx, nsply;
-    extern double passoN, passoE;
+    extern double stepN, stepE;
     struct line_pnts *point;
     struct line_cats *categories;
 
@@ -36,7 +36,7 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 
 	if (Vect_point_in_box(obs[i][0], obs[i][1], mean, &General)) {
 	    interpolation =
-		dataInterpolateBicubic(obs[i][0], obs[i][1], passoE, passoN,
+		dataInterpolateBicubic(obs[i][0], obs[i][1], stepE, stepN,
 				       nsplx, nsply, Elaboration.west,
 				       Elaboration.south, parBilin);
 
@@ -123,7 +123,6 @@ void P_Outlier(struct Map_info *Out, struct Map_info *Outlier,
 			else {
 			    Vect_write_line(Outlier, GV_POINT, point,
 					    categories);
-			    G_message("here we are");
 			}
 		    }
 		    else if ((*point->y < Overlap.S) && (*point->y > General.S)) {	/*(2) */

vector/lidar/v.outlier/v.outlier.html

@@ -32,6 +32,13 @@ Now, the outlier removal uses the default threshold and there is also
 an output vector available for visualizaton in QGIS
  (<a href="http://www.qgis.org">http://www.qgis.org</a>).
 
+<h4>North carolina location example</h4>
+
+<div class="code"><pre>
+v.outlier --o --v input=elev_lid792_bepts@PERMANENT output=elev_lid792_bepts_nooutliers outlier=elev_lid792_bepts_outliers soe=5 son=5 thres_o=0.1
+</pre></div>
+
+
 <h2>NOTES</h2>
 
 This module is designed to work with LIDAR data, so not topology is

vector/lidar/v.surf.bspline/crosscorr.c

@@ -47,7 +47,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
 {
     int bilin = TRUE;		/*booleans */
     int nsplx, nsply, nparam_spl, ndata;
-    double *mean, *rms, *stdev, rms_min, stdev_min;
+    double *mean, *rms, *stdev;
 
     /* double lambda[PARAM_LAMBDA] = { 0.0001, 0.001, 0.01, 0.1, 1.0, 10.0 }; */	/* Fixed values (by the moment) */
     double lambda[PARAM_LAMBDA] = { 0.01, 0.05, 0.1, 0.2, 0.3, 0.4 };	/* Fixed values (by the moment) */
@@ -94,7 +94,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
 
     if (ndata > 0) {		/* If at least one point is in the region */
 	int i, j, lbd;		/* lbd: lambda index */
-	int BW, lbd_min;	/* lbd_min: index where minimun is found */
+	int BW;	
 	double mean_reg, *obs_mean;
 
 	int nrec, ctype = 0;
@@ -172,7 +172,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
 	    /*
 	       How cross correlation algorithm is done:
 	       For each cicle, only the first ndata-1 "observ" elements are considered for the 
-	       interpolation. Within every interpolation mean is calculated to lowering border 
+	       interpolation. Within every interpolation mean is calculated to lowering edge 
 	       errors. The point let out will be used for an estimation. The error between the 
 	       estimation and the observation is recorded for further statistics.
 	       At the end of the cicle, the last point, that is, the ndata-1 index, and the point 
@@ -257,7 +257,7 @@ int cross_correlation(struct Map_info *Map, double passWE, double passNS)
 		   if (bilin) interpolation (&interp, P_BILINEAR);
 		   else interpolation (&interp, P_BICUBIC);
 		 */
-		tcholSolve(N, TN, parVect, nparam_spl, BW);
+		G_math_solver_cholesky_sband(N, parVect, TN, nparam_spl, BW);
 
 		/* Estimation of j-point */
 		if (bilin)
@@ -335,23 +335,23 @@ void interpolation(struct ParamInterp *interp, boolean bilin)
 {
     if (bilin == P_BILINEAR) {	/* Bilinear interpolation */
 	normalDefBilin(interp->N, interp->TN, interp->Q, interp->obsVect,
-		       interp->passoE, interp->passoN, interp->nsplx,
+		       interp->stepE, interp->stepN, interp->nsplx,
 		       interp->nsply, interp->region.west,
 		       interp->region.south, interp->ndata,
 		       interp->nparam_spl, interp->BW);
 
 	nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
-		     interp->nsply, interp->passoE, interp->passoN);
+		     interp->nsply, interp->stepE, interp->stepN);
     }
     else {			/* Bicubic interpolation */
 	normalDefBicubic(interp->N, interp->TN, interp->Q, interp->obsVect,
-			 interp->passoE, interp->passoN, interp->nsplx,
+			 interp->stepE, interp->stepN, interp->nsplx,
 			 interp->nsply, interp->region.west,
 			 interp->region.south, interp->ndata,
 			 interp->nparam_spl, interp->BW);
 
 	nCorrectGrad(interp->N, interp->lambda[lbd], interp->nsplx,
-		     interp->nsply, interp->passoE, interp->passoN);
+		     interp->nsply, interp->stepE, interp->stepN);
     }
     return TRUE;
 }

vector/lidar/v.surf.bspline/main.c

@@ -43,8 +43,8 @@ int main(int argc, char *argv[])
     int nsubregion_col, nsubregion_row, subregion_row, subregion_col;
     int subregion = 0, nsubregions = 0;
     int last_row, last_column, grid, bilin, ext, flag_auxiliar, cross;	/* booleans */
-    double passoN, passoE, lambda, mean;
-    double N_extension, E_extension, orloE, orloN;
+    double stepN, stepE, lambda, mean;
+    double N_extension, E_extension, edgeE, edgeN;
 
     const char *mapset, *dvr, *db, *vector, *map;
     char table_name[GNAME_MAX], title[64];
@@ -62,8 +62,8 @@ int main(int argc, char *argv[])
     struct History history;
 
     struct GModule *module;
-    struct Option *in_opt, *in_ext_opt, *out_opt, *out_map_opt, *passoE_opt,
-	*passoN_opt, *lambda_f_opt, *type_opt, *dfield_opt, *col_opt;
+    struct Option *in_opt, *in_ext_opt, *out_opt, *out_map_opt, *stepE_opt,
+	*stepN_opt, *lambda_f_opt, *type_opt, *dfield_opt, *col_opt;
     struct Flag *cross_corr_flag, *spline_step_flag, *withz_flag;
 
     struct Reg_dimens dims;
@@ -120,23 +120,23 @@ int main(int argc, char *argv[])
     out_map_opt->key = "raster";
     out_map_opt->required = NO;
 
-    passoE_opt = G_define_option();
-    passoE_opt->key = "sie";
-    passoE_opt->type = TYPE_DOUBLE;
-    passoE_opt->required = NO;
-    passoE_opt->answer = "4";
-    passoE_opt->description =
+    stepE_opt = G_define_option();
+    stepE_opt->key = "sie";
+    stepE_opt->type = TYPE_DOUBLE;
+    stepE_opt->required = NO;
+    stepE_opt->answer = "4";
+    stepE_opt->description =
 	_("Length of each spline step in the east-west direction");
-    passoE_opt->guisection = _("Settings");
-
-    passoN_opt = G_define_option();
-    passoN_opt->key = "sin";
-    passoN_opt->type = TYPE_DOUBLE;
-    passoN_opt->required = NO;
-    passoN_opt->answer = "4";
-    passoN_opt->description =
+    stepE_opt->guisection = _("Settings");
+
+    stepN_opt = G_define_option();
+    stepN_opt->key = "sin";
+    stepN_opt->type = TYPE_DOUBLE;
+    stepN_opt->required = NO;
+    stepN_opt->answer = "4";
+    stepN_opt->description =
 	_("Length of each spline step in the north-south direction");
-    passoN_opt->guisection = _("Settings");
+    stepN_opt->guisection = _("Settings");
 
     type_opt = G_define_option();
     type_opt->key = "method";
@@ -186,8 +186,8 @@ int main(int argc, char *argv[])
     else
 	bilin = P_BICUBIC;
 
-    passoN = atof(passoN_opt->answer);
-    passoE = atof(passoE_opt->answer);
+    stepN = atof(stepN_opt->answer);
+    stepE = atof(stepE_opt->answer);
     lambda = atof(lambda_f_opt->answer);
 
     flag_auxiliar = FALSE;
@@ -259,7 +259,7 @@ int main(int argc, char *argv[])
     /* Cross-correlation begins */
     if (cross_corr_flag->answer) {
 	G_debug(1, "CrossCorrelation()");
-	cross = cross_correlation(&In, passoE, passoN);
+	cross = cross_correlation(&In, stepE, stepN);
 
 	if (cross != TRUE)
 	    G_fatal_error(_("Cross validation didn't finish correctly"));
@@ -268,7 +268,7 @@ int main(int argc, char *argv[])
 
 	    Vect_close(&In);
 
-	    G_done_msg(_("Cross validation finished for sie = %f and sin = %f"), passoE, passoN);
+	    G_done_msg(_("Cross validation finished for sie = %f and sin = %f"), stepE, stepN);
 	    exit(EXIT_SUCCESS);
 	}
     }
@@ -406,7 +406,7 @@ int main(int argc, char *argv[])
       | Each original region will be divided into several subregions. 
       | Each one will be overlaped by its neighbouring subregions. 
       | The overlapping is calculated as a fixed OVERLAP_SIZE times
-      | the largest spline step plus 2 * orlo
+      | the largest spline step plus 2 * edge
       ----------------------------------------------------------------*/
 
     /* Fixing parameters of the elaboration region */
@@ -414,25 +414,25 @@ int main(int argc, char *argv[])
 
     nsplx_adj = NSPLX_MAX;
     nsply_adj = NSPLY_MAX;
-    if (passoN > passoE)
-	dims.overlap = OVERLAP_SIZE * passoN;
+    if (stepN > stepE)
+	dims.overlap = OVERLAP_SIZE * stepN;
     else
-	dims.overlap = OVERLAP_SIZE * passoE;
-    P_get_orlo(bilin, &dims, passoE, passoN);
-    P_set_dim(&dims, passoE, passoN, &nsplx_adj, &nsply_adj);
+	dims.overlap = OVERLAP_SIZE * stepE;
+    P_get_edge(bilin, &dims, stepE, stepN);
+    P_set_dim(&dims, stepE, stepN, &nsplx_adj, &nsply_adj);
 
     G_verbose_message(_("adjusted EW splines %d"), nsplx_adj);
     G_verbose_message(_("adjusted NS splines %d"), nsply_adj);
 
     /* calculate number of subregions */
-    orloE = dims.latoE - dims.overlap - 2 * dims.orlo_v;
-    orloN = dims.latoN - dims.overlap - 2 * dims.orlo_h;
+    edgeE = dims.ew_size - dims.overlap - 2 * dims.edge_v;
+    edgeN = dims.sn_size - dims.overlap - 2 * dims.edge_h;
 
     N_extension = original_reg.north - original_reg.south;
     E_extension = original_reg.east - original_reg.west;
 
-    nsubregion_col = ceil(E_extension / orloE) + 0.5;
-    nsubregion_row = ceil(N_extension / orloN) + 0.5;
+    nsubregion_col = ceil(E_extension / edgeE) + 0.5;
+    nsubregion_row = ceil(N_extension / edgeN) + 0.5;
 
     if (nsubregion_col < 0)
 	nsubregion_col = 0;
@@ -470,7 +470,7 @@ int main(int argc, char *argv[])
 
 	nsply =
 	    ceil((elaboration_reg.north -
-		  elaboration_reg.south) / passoN) + 0.5;
+		  elaboration_reg.south) / stepN) + 0.5;
 	G_debug(1, "Interpolation: nsply = %d", nsply);
 	/*
 	if (nsply > NSPLY_MAX)
@@ -505,7 +505,7 @@ int main(int argc, char *argv[])
 	    }
 	    nsplx =
 		ceil((elaboration_reg.east -
-		      elaboration_reg.west) / passoE) + 0.5;
+		      elaboration_reg.west) / stepE) + 0.5;
 	    G_debug(1, "Interpolation: nsplx = %d", nsplx);
 	    /*
 	    if (nsplx > NSPLX_MAX)
@@ -601,25 +601,25 @@ int main(int argc, char *argv[])
 		    G_debug(1,
 			    "Interpolation: (%d,%d): Bilinear interpolation...",
 			    subregion_row, subregion_col);
-		    normalDefBilin(N, TN, Q, obsVect, passoE, passoN, nsplx,
+		    normalDefBilin(N, TN, Q, obsVect, stepE, stepN, nsplx,
 				   nsply, elaboration_reg.west,
 				   elaboration_reg.south, npoints,
 				   nparameters, BW);
-		    nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
+		    nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
 		}
 		/* Bicubic interpolation */
 		else {
 		    G_debug(1,
 			    "Interpolation: (%d,%d): Bicubic interpolation...",
 			    subregion_row, subregion_col);
-		    normalDefBicubic(N, TN, Q, obsVect, passoE, passoN, nsplx,
+		    normalDefBicubic(N, TN, Q, obsVect, stepE, stepN, nsplx,
 				     nsply, elaboration_reg.west,
 				     elaboration_reg.south, npoints,
 				     nparameters, BW);
-		    nCorrectGrad(N, lambda, nsplx, nsply, passoE, passoN);
+		    nCorrectGrad(N, lambda, nsplx, nsply, stepE, stepN);
 		}
 
-		tcholSolve(N, TN, parVect, nparameters, BW);
+		G_math_solver_cholesky_sband(N, parVect, TN, nparameters, BW);
 
 		G_free_matrix(N);
 		G_free_vector(TN);
@@ -631,7 +631,7 @@ int main(int argc, char *argv[])
 		    raster_matrix =
 			P_Regular_Points(&elaboration_reg, general_box,
 					 overlap_box, raster_matrix, parVect,
-					 passoN, passoE, dims.overlap, mean,
+					 stepN, stepE, dims.overlap, mean,
 					 nsplx, nsply, nrows, ncols, bilin);
 		}
 		else {		/* OBSERVATION POINTS INTERPOLATION */
@@ -640,7 +640,7 @@ int main(int argc, char *argv[])
 				subregion_row, subregion_col);
 			P_Sparse_Points(&Out, &elaboration_reg, general_box,
 					overlap_box, obsVect, parVect,
-					lineVect, passoE, passoN,
+					lineVect, stepE, stepN,
 					dims.overlap, nsplx, nsply, npoints,
 					bilin, Cats, driver, mean,
 					table_name);
@@ -672,7 +672,7 @@ int main(int argc, char *argv[])
 				subregion_row, subregion_col);
 			P_Sparse_Points(&Out, &elaboration_reg, general_box,
 					overlap_box, obsVect_ext, parVect,
-					lineVect_ext, passoE, passoN,
+					lineVect_ext, stepE, stepN,
 					dims.overlap, nsplx, nsply,
 					npoints_ext, bilin, Cats, driver,
 					mean, table_name);

vector/lidar/v.surf.bspline/v.surf.bspline.html

@@ -128,6 +128,12 @@ v.surf.bspline input=point_vector raster=interpolate_surface layer=1 column=attr
 The interpolation will be done using the values in attrib_column, in the
 table associated with layer 1.
 
+<h4>North carolina location example using Z-coordinates for interpolation</h4>
+<div class="code"><pre>
+v.surf.bspline --o --v input=elev_lid792_bepts@PERMANENT raster=elev_lid792_rast sie=5 sin=5 method=bicubic lambda_i=0.1 -z
+</pre></div>
+
+
 <h2>BUGS</h2>
 Known issues:
 <p>