/* Last edited on 2025-02-18 19:49:30 by stolfi */ /* ---------------------------------------------------------------------- */ /* pst_integrate.{h,c} */ /* Build the index {ix[0..N-1]} and inverse indices {col[0..N-1],row[0..N-1]}: */ int32_t *ix = talloc(NXY_Z, int32_t); int32_t *col = talloc(NXY_Z, int32_t); int32_t *row = talloc(NXY_Z, int32_t); for (uint32_t xy = 0; xy < NXY_Z; xy++) { ix[xy] = (int32_t)xy; col[xy] = xy % (uint32_t)NX_Z; row[xy] = xy / (uint32_t)NX_Z; } void compute_axial_edge_delta_and_weight(uint32_t axis, int32_t x, int32_t y, int32_t u, double *dP, double *wP) { assert((u == -1) || (u == +1)); (*dP) = 0.0; (*wP) = 0.0; /* If edge does not exist. */ if (axis == 0) { /* X derivative */ int32_t xg; /* Indices of relevant samples in {G,W}. */ if (u < 0) { xg = x - 1; if (xg < 0) { return; } } else { xg = x; if (xg >= NX_G) { return; } } pst_interpolate_four_samples(G,W, 0, xg, y-1, xg, y, dP,wP); (*dP) = (*dP) * u; } else if (axis == 1) { /* Y derivative */ int32_t yg; /* Indices of relevant samples in {G,W}. */ if (u < 0) { yg = y - 1; if (yg < 0) { return; } } else { yg = y; if (yg >= NY_G) { return; } } pst_interpolate_four_samples(G,W, 1, x-1,yg, x,yg, dP,wP); (*dP) = (*dP) * u; } else { demand(FALSE, "invalid axis"); } return; } /* Get the axial deltas and their weights: */ double wmo, wpo, wom, wop; /* Edge weights. */ double dmo, dpo, dom, dop; /* Edge deltas. */ compute_axial_edge_delta_and_weight(G, W, x,y, 0, -1, &dmo,&wmo); compute_axial_edge_delta_and_weight(G, W, x,y, 0, +1, &dpo,&wpo); compute_axial_edge_delta_and_weight(G, W, x,y, 1, -1, &dom,&wom); compute_axial_edge_delta_and_weight(G, W, x,y, 1, +1, &dop,&wop); assert((! isnan(wmo)) && (wmo >= 0)); assert((! isnan(wpo)) && (wpo >= 0)); assert((! isnan(wom)) && (wom >= 0)); assert((! isnan(wop)) && (wop >= 0)); /* Get the skew deltas and weights: */ double wmm, wmp, wpm, wpp; /* Edge weights. */ double dmm, dmp, dpm, dpp; /* Edge deltas. */ compute_skew_edge_delta_and_weight(x,y, -1,-1, &dmm,&wmm); compute_skew_edge_delta_and_weight(x,y, -1,+1, &dmp,&wmp); compute_skew_edge_delta_and_weight(x,y, +1,-1, &dpm,&wpm); compute_skew_edge_delta_and_weight(x,y, +1,+1, &dpp,&wpp); assert((! isnan(wmm)) && (wmm >= 0)); assert((! isnan(wmp)) && (wmp >= 0)); assert((! isnan(wpm)) && (wpm >= 0)); assert((! isnan(wpp)) && (wpp >= 0)); /* Select which mismatch terms to include in the {Q} function: */ bool_t use_mo = FALSE, use_po = FALSE, use_om = FALSE, use_op = FALSE; bool_t use_mm = FALSE, use_mp = FALSE, use_pm = FALSE, use_pp = FALSE; if (wmo > 0) { use_mo = TRUE; } else if ( wmm > 0) { use_mm = TRUE; } else if (wmp > 0) { use_mp = TRUE; } if (wpo > 0) { use_po = TRUE; } else if ( wpm > 0) { use_pm = TRUE; } else if (wpp > 0) { use_pp = TRUE; } if (wom > 0) { use_om = TRUE; } else if ( wmm > 0) { use_mm = TRUE; } else if (wpm > 0) { use_pm = TRUE; } if (wop > 0) { use_op = TRUE; } else if ( wmp > 0) { use_mp = TRUE; } else if (wpp > 0) { use_pp = TRUE; } uint32_t n_use = (uint32_t)use_mo + (uint32_t)use_po + (uint32_t)use_om + (uint32_t)use_op + (uint32_t)use_mm + (uint32_t)use_mp + (uint32_t)use_pm + (uint32_t)use_pp; if (n_use == 0) { /* Vertex is isolated. */ if (fillHoles) { /* Provide axial equations for existing edges. */ double wtiny = 1.0e-50; wmo = (x <= 0 ? 0.0 : wtiny); dmo = 0.0; use_mo = (wmo > 0); wpo = (x >= NX_G ? 0.0 : wtiny); dpo = 0.0; use_po = (wpo > 0); wom = (y <= 0 ? 0.0 : wtiny); dom = 0.0; use_om = (wom > 0); wop = (y >= NY_G ? 0.0 : wtiny); dop = 0.0; use_op = (wop > 0); n_use = use_mo + use_po + use_om + use_op; assert(n_use > 0); } else { /* Leave it isolated: */ } } /* Combine the selected mismatch terms in one normalized equilibrium equation: */ /* Also compute the total vertex weight weight {wtot}: */ pst_imgsys_equation_t *eqk = &(eq[xy]); /* Initialize equation: */ eqk->nt = 0; /* Number of terms in equation. */ eqk->rhs = 0.0; /* Right-hand side. */ eqk->wtot = 0.0; /* Weight of equation. */ /* Set the diagonal term of equation (already normalized): */ eqk->uid[0] = xy; eqk->cf[0] = 1.0; eqk->nt = 1; if (n_use > 0) { assert((1+n_use) <= MAX_COEFFS); /* Collect non-diagonal equation terms of selected mismtch terms with unnormalized coeffs: */ if ( use_mo ) { append_equation_term(eqk, x,y, -1,00, dmo,wmo); } if ( use_po ) { append_equation_term(eqk, x,y, +1,00, dpo,wpo); } if ( use_om ) { append_equation_term(eqk, x,y, 00,-1, dom,wom); } if ( use_op ) { append_equation_term(eqk, x,y, 00,+1, dop,wop); } if ( use_mm ) { append_equation_term(eqk, x,y, -1,-1, dmm,wmm); } if ( use_mp ) { append_equation_term(eqk, x,y, -1,+1, dmp,wmp); } if ( use_pm ) { append_equation_term(eqk, x,y, +1,-1, dpm,wpm); } if ( use_pp ) { append_equation_term(eqk, x,y, +1,+1, dpp,wpp); } assert(eqk->nt == 1 + n_use); /* Compute {wtot} and normalize coeficients: */ assert(eqk->wtot > 0); for (uint32_t j = 1; j < eqk->nt; j++) { eqk->cf[j] /= eqk->wtot; } eqk->rhs /= eqk->wtot; } } int32_t xg, yg; /* Indices of relevant samples in {G,W}. */ (*dP) = 0.0; (*wP) = 0.0; /* If edge does not exist. */ if (ux < 0) { xg = x - 1; if (xg < 0) { return; } } else { xg = x; if (xg >= NX_G) { return; } } if (uy < 0) { yg = y - 1; if (yg < 0) { return; } } else { yg = y; if (yg >= NY_G) { return; } } double dx = float_image_get_sample(G, 0, xg, yg); double dy = float_image_get_sample(G, 1, xg, yg); double d = dx*ux + dy*uy; double w = (W == NULL ? 1.0 : float_image_get_sample(W, 0, xg, yg)); assert((! isnan(w)) && (w >= 0)); (*dP) = d; (*wP) = w; } /* ---------------------------------------------------------------------- */ /* pst_integrate */ auto void compute_axial_edge_delta_and_weight(uint32_t axis, int32_t x, int32_t y, int32_t u, double *dP, double *wP); /* Computes the delta {*dP} and the weight {*wP} of the axial edge that goes from {(x,y)} to {(x+u,y)} if {axis=0} or to {(x,y+u)} if {axis=1}, where {u} is either {+1} or {-1}. Assumes that {x,y} is a vertex of the cell grid. */ auto void compute_skew_edge_delta_and_weight(int32_t x, int32_t y, int32_t ux, int32_t uy, double *dP, double *wP); /* Computes the delta {*dP} and the weight {*wP} of the skew edge that goes from {(x,y)} to {(x+ux,y+uy)}, where {ux} and {uy} are either {+1} or {-1}. Assumes that {x,y} is a vertex of the cell grid. */ void pst_integrate_solve_system ( pst_imgsys_t *S, float_image_t *OZ, uint32_t ord[], uint32_t maxIter, double convTol, bool_t para, bool_t szero, bool_t verbose, uint32_t level, uint32_t reportIter, pst_integrate_report_heights_proc_t *reportHeights ); /* Solves system {S} by Gauss-Seidel iteration, and stores the solution into the height map {OZ}. On input {OZ} should contain the initial guess. The vector {ord[0..S->N-1]} specifies the order in which the equations are to be solved. The image {OZ} must be single-channel. The procedure assumes that the unknown height values of {S} are related to the samples of {OZ}, as defined by the tables {uid[0..NX*NY-1]}, {col[0..S->N-1]}, and {row[0..S->N-1]}. If {reportHeights} is not null and {reportIter} is not zero, will call {reportZ(level,iter,change,final,OZ)} before the first iteration, after the last iteration, and after any number of iterations that is divisible by {reportIter}. */ void pst_integrate_solve_system ( pst_imgsys_t *S, float_image_t *OZ, uint32_t ord[], uint32_t maxIter, double convTol, bool_t para, bool_t szero, bool_t verbose, uint32_t level, uint32_t reportIter, pst_integrate_report_heights_proc_t *reportHeights ) { uint32_t indent = 2*level+2; auto void reportSol(uint32_t iter, double change, bool_t final, uint32_t N, double Z[]); /* A procedure that is called by the Gauss-Seidel solver at each iteration. When appropriate, it copies {Z} into the image {OZ} and calls {reportHeights}. This happens when {final} is TRUE or when {reportIter != 0} and {iter} is a multiple of {reportIter}. */ void reportSol(uint32_t iter, double change, bool_t final, uint32_t N, double Z[]) { bool_t doit = final || ((reportHeights != NULL) && (reportIter != 0) && (iter % reportIter == 0)); if (doit) { pst_imgsys_copy_sol_vec_to_image(S, Z, OZ, 0.0); if (reportHeights != NULL) { reportHeights(level, iter, change, final, OZ); } } } double *Z = rn_alloc(S->N); pst_imgsys_copy_image_to_sol_vec(S, OZ, Z, 0.0); pst_imgsys_solve_iterative(S, Z, ord, maxIter, convTol, para, szero, verbose, indent, reportSol); /* {reportSol} must have copied the final solution into {OZ}. */ free(Z); } void pst_integrate_fill_holes(pst_imgsys_t *S, int32_t NX_Z, int32_t NY_Z); /* Modifies every equation {eqk=S.eq[k]} that is indeterminate {0*Z[x,y]=0} (which should have with zero total weight {eqk.wtot}) so that it will set the corresponding unknown height {h[k]} to the average of the heights of its neighors, with equal weights and an arbitrary positive {wtot}. The equations of those neighbors should not depend on {Z[x,y]}, and are not modified -- unless they too are indeterminate. The effect of this fudging is to fill in any indeterminate regions of the height map with a smooth surface that is attached at the edges to the determinate heights. The indices {S.col[0..S.N-1]} and {S.row[0..S.N-1]} must not be {-1} and must be in the ranges {0..NX_Z-1} and {0..NY_Z-1}, respectively. Assumes that, for every variable/equation {uid} in {0..S.N-1}, the indices {x=S.col[k]} and {y=S.row[k]} are such that {uid = x + y*NX_X}. */ void pst_integrate_fill_holes(pst_imgsys_t *S, int32_t NX_Z, int32_t NY_Z) { uint32_t N = S->N; auto void add_term(pst_imgsys_equation_t *eq, int32_t x1, int32_t y1); /* If the variable {x1,y1} exists, adds to {eq} a term with coefficient 1 that pulls the variable of {eq} towards the mean of its neighbors. */ for (uint32_t k = 0; k < N; k++) { uint32_t uidk = k; pst_imgsys_equation_t *eqk = &(S->eq[k]); if (pst_imgsys_equation_is_null(uidk, eqk, N)) { assert(eqk->nt == 1); assert(eqk->uid[0] == uidk); /* Get neighbors if variable {uidk}: */ int32_t x = S->col[uidk]; int32_t y = S->row[uidk]; assert((x >= 0) && (x < NX_Z)); assert((y >= 0) && (y < NY_Z)); add_term(eqk, x-1,y); add_term(eqk, x+1,y); add_term(eqk, x,y-1); add_term(eqk, x,y+1); } eqk->wtot = 1.0e-5; } return; void add_term(pst_imgsys_equation_t *eq, int32_t x1, int32_t y1) { if ((x1 < 0) || (x1 >= NX_Z)) { return; } if ((y1 < 0) || (y1 >= NY_Z)) { return; } /* Variable {x1,y1} exists: */ int32_t uid1 = x1 + y1*NX_Z; assert((uid1 >= 0) && (uid1 < S->N)); assert((uid1 >= 0) && (uid1 < S->N)); assert((eq->nt >= 1) && (eq->nt < MAX_COEFFS)); int32_t j = (int32_t)eq->nt; eq->uid[j] = (uint32_t)uid1; eq->cf[j] = -1.0; eq->cf[0] += 1.0; (eq->nt)++; } } void pst_integrate_set_equations(pst_imgsys_t *S, float_image_t *G, float_image_t *W, bool_t verbose); /* Builds the equation list {S.eq} for the integration system {S}, given the gradient map {G} and its weight map {W}.*/ /* ---------------------------------------------------------------------- */ /* pst_imgsys.{h,c} */ /* The table {S.uid[0..NX*NY-1]} is the inverse mapping, from the col and row indices {(x,y)} to the corresponding unknown index {k = S.uid[x+y*NX]}; or {-1} if there is no unknown corresponding to that col/row pair. */ /* The array {uid} must have {NX*NY} elements, and must be the inverse of {col,row} as explained under {pst_imgsys_t} */ /* Conversely, for each {x} in {0..NX} and each {y} in {0..NY}, the height map element {Z[x,y]} corresponds to the height value {h[k]} where {k = S.uid[x + y*(NX+1)]}, if {k >= 0}; os is omitted from the system, if {k<0}. */ for (int32_t xy = 0; xy < NX*NY; xy++) { int32_t x = xy % NX; int32_t y = xy / NX; int32_t uid = S->uid[xy]; if (uid != -1) { demand((uid >= 0) && (uid < N), "invalid {uid[k]}"); demand(x == S->col[uid], "inconsistent {S.uid[x,y],S.col[k]"); demand(y == S->row[uid], "inconsistent {S.iy[x,y],S.row[k]"); } } int32_t NXY = S->NX*S->NY; /* Entries in the {S.uid} table. */ /* Beware that there may be equations with {x,y=-1,-1} that are not counted and not in the {ix} table. Thus we must create a separate dictionary {uid_map} to map old unknown/equation indices to new ones. We cannot rely on {S.uid} for that. */ /* Each pixel of the image {U} will be the assumed reliability for the corresponding pixel of the height map, which is just the total weight {eq.wtot} of the corresponding equation. The {S.col} and {S.row} tables are used to map indices of variables and equations to height map pixel indices. */ void pst_imgsys_extract_system_eq_tot_weight_image(pst_imgsys_t *S, float_image_t *U, float vdef); /* Stores into channel 0 of {U} the equation weights {.wtot} of the system {S}. Specifically, fills channel 0 of {U} with {vdef}, then copies each equation weight {S.eq[k].wtot} into {U[0,x,y]}, where {x} and {y} are {S->col[k]} and {S->row[k]}. However, skips the assignment if {x} and {y} are {-1} (meaning that the equation {k} is not associated to an heigh map pixel). */ void pst_imgsys_extract_system_eq_tot_weight_image(pst_imgsys_t *S, float_image_t *U, float vdef) { int32_t NC, NX, NY; float_image_get_size(U, &NC, &NX, &NY); demand(NC == 1, "image should have only one channel"); float_image_fill_channel(U, 0, vdef); for (uint32_t uid = 0; uid < S->N; uid++) { int32_t x = S->col[uid]; int32_t y = S->row[uid]; assert((x == -1) == (y == -1)); if (x != -1) { demand((x >= 0) && (x < NX), "invalid image col in system"); demand((y >= 0) && (y < NY), "invalid image row in system"); double wtot = S->eq[uid].wtot; float_image_set_sample(U, 0, x,y, (float)wtot); } } } void pst_imgsys_remove_holes(pst_imgsys_t *S, int32_t *nuid_from_ouid); /* Removes from the system {S} any variable {Z[k]=Z[x,y} whose equation {S->eq[k]} is indeterminate, namely {0*Z[k]=0}. This will reduce the number {S->N} (but not {S->NX} and {S->NY}), renumber the equations and variables to the new range {0..S->N-1}, and set {S->uid[x + y*NX]} to {-1}. If not {NULL}, the parameter {nuid_from_ouid} must be a table of {N_in} entries where {N_in} is {S.N} upon entry. In that case, the procedure will set {nuid_from_ouid[oi]} to the new index of the variable/equation whose old index was {oi}; or to -1 if that variable/equation was deleted. */ void pst_imgsys_remove_holes(pst_imgsys_t *S, int32_t *nuid_from_ouid) { bool_t debug = TRUE; uint32_t N_in = S->N; /* Provide a local table if the client did not give one: */ int32_t *nfo_old = nuid_from_ouid; if (nuid_from_ouid == NULL) { nuid_from_ouid = talloc(N_in, int32_t); } for (int32_t k = 0; k < N_in; k++) { nuid_from_ouid[k] = -1; } /* Get number {N-ot} of non-null equations and set {nuid_from_ouid}: */ uint32_t N_ot = 0; /* Non-null equations. */ uint32_t N_ex = 0; /* Number of equations excluded. */ for (uint32_t k = 0; k < N_in; k++) { uint32_t uidk_old = k; pst_imgsys_equation_t *eqk_old = &(S->eq[k]); if (pst_imgsys_equation_is_null(uidk_old, eqk_old, N_in)) { /* Will exclude it: */ if (debug) { fprintf(stderr, " excluding null height/equation %d\n", uidk_old); } N_ex++; } else { /* Will keep it: */ uint32_t uidk_new = N_ot; nuid_from_ouid[uidk_old] = (int32_t)uidk_new; N_ot++; } } assert(N_ot + N_ex == N_in); /* Paranoia. */ if (N_ex > 0) { /* Now create the equation, colum, row tables for a new system and copy the non-null equations into them: */ pst_imgsys_equation_t *eq_ot = talloc(N_ot, pst_imgsys_equation_t); int32_t *col_ot = talloc(N_ot, int32_t); int32_t *row_ot = talloc(N_ot, int32_t); uint32_t k_new = 0; for (uint32_t k = 0; k < N_in; k++) { uint32_t uidk_old = k; pst_imgsys_equation_t *eqk = &(S->eq[k]); if (! pst_imgsys_equation_is_null(uidk_old, eqk, N_in)) { /* Will keep it: */ uint32_t uidk_new = k_new; assert(nuid_from_ouid[uidk_old] == uidk_new); eq_ot[k_new] = (*eqk); /* Remap terms; for good measure, discard null ones: */ uint32_t nt = eqk->nt; assert((nt >= 2) && (nt <= MAX_COEFFS)); uint32_t nt_new = 0; for (int32_t j = 0; j < nt; j++) { if ((j == 0) || (fabs(eqk->cf[j]) >= FLUFF)) { uint32_t uidj_old = eqk->uid[j]; int32_t j_new = (int32_t)nt_new; int32_t uidj_new = nuid_from_ouid[uidj_old]; assert(uidj_new != -1); assert((uidj_new >= 0) && (uidj_new < N_ot)); eqk->uid[j_new] = (uint32_t)uidj_new; eqk->cf[j_new] = eqk->cf[j]; nt_new++; } } assert(nt_new >= 2); eqk->nt = nt_new; assert(eqk->uid[0] == uidk_new); /* Remap {col,row,uid}: */ int32_t x = S->col[uidk_old]; int32_t y = S->row[uidk_old]; col_ot[uidk_new] = x; row_ot[uidk_new] = y; /* One more: */ k_new++; } } assert(k_new == N_ot); /* Replace the tables: */ S->N = N_ot; free(S->eq); S->eq = eq_ot; free(S->col); S->col = col_ot; free(S->row); S->row = row_ot; pst_imgsys_check_valid(S, N_ot, -1, -1); } if (nuid_from_ouid != nfo_old) { free(nuid_from_ouid); } return; } /* ---------------------------------------------------------------------- */ /* pst_integrate_iterative.{h,c} */ /* The {OZ} image is returned in {*OZP}; it will be allocated by the procedure and will be bigger than {IG} by one row and one column. If {keepNull} is true, the system is fudged so that any height map pixel that is completely surrounded by zero-weight data (so that its equation has only one term with weight zero) will be ultimately be set to the average of their neighbors. Thus the regions of weight zero in the input data will be filled with a Laplace equilibrium surface. This fudging does not affect the final result for pixels with non-null equations. If {keepNull} is false, those pixels will be excluded from the system, and will be set to zero in the final height map. */ void get_initial_guess(void) { if (verbose) { fprintf(stderr, "%*sobtaining initial guess of heights from {Z} ...\n", indent, ""); } for (int32_t k = 0; k < S->N; k++) { /* Index of equation's main variable: */ int32_t uidk = k; /* Compute indices of pixel corresponding to the variable {Z[k]}: */ int32_t xk = S->col[uidk]; int32_t yk = S->row[uidk]; double wzk = float_image_get_sample(Z, 1, xk, yk); demand(isfinite(wzk) && (wzk >= 0), "invalid {Z} weight"); double zk = (wzk == 0.0 ? 0.0 : float_image_get_sample(Z, 0, xk, yk)); demand(isfinite(zk), "invalid {Z} height value"); double whk = (H == NULL ? 0.0 : float_image_get_sample(H, 1, xk, yk)); demand(isfinite(whk) && (whk >= 0), "invalid {H} weight"); double hk = (whk == 0.0 ? 0.0 : float_image_get_sample(H, 0, xk, yk)); demand(isfinite(hk), "invalid {H} height value"); double wtk = wzk + whk; z[uidk] = (wtk == 0.0 ? 0.0 : (wzk*zk + whk*hk)/wtk); } } /* ---------------------------------------------------------------------- */ /* pst_integrate_recursive */ /* If {keepNull} is true, each system is fudged so that any height map pixel that is completely surrounded by zero-weight data (so that its equation has only one term with weight zero) will be ultimately be set to the average of their neighbors. Thus the regions of weight zero in the input data will be filled with a Laplace equilibrium surface. This fudging does not affect the final result for pixels with non-null equations. If {keepNull} is false, those pixels will be excluded from the system, and will be set to zero in the final height map. If {U} is not null, it must be a single-channel image with the same weight map with the same size as the height map {Z}. Each sample of {U} will be set to the weight {wtot} of the equation that defines the corresponding sample of {Z}. Pixels of {Z} that end up with no equation (if {keepHoles} is false) are set to zero. */ if (trivial) { /* Solution for the 1x1 system: */ /* Get slopes: */ float dZdX = float_image_get_sample(G, 0, 0, 0); float dZdY = float_image_get_sample(G, 1, 0, 0); /* Assume an affine function that is 0 at center, with the given gradient: */ double z10 = 0.5*(dZdX - dZdY); double z11 = 0.5*(dZdX + dZdY); float_image_set_sample(Z, 0, 0, 0, (float)-z11); float_image_set_sample(Z, 0, 0, 1, (float)-z10); float_image_set_sample(Z, 0, 1, 0, (float)+z10); float_image_set_sample(Z, 0, 1, 1, (float)+z11); float_image_fill_channel(Z, 1, 1.0); /* The initial guess is the solution: */ if (reportSys != NULL) { reportSys(level); } if (reportHeights != NULL) { reportHeights(level, 0, 0.0, TRUE, Z); } } /* Build the linear system: */ uint32_t NXY_Z = (uint32_t)(NX_Z*NY_Z); if (verbose) { fprintf(stderr, "%*sbuilding linear system for %d pixels ...\n", indent, "", NXY_Z); } pst_imgsys_t *S = pst_integrate_build_system(G, H, verbose); assert(S->N == NXY_Z); if (verbose) { fprintf(stderr, "%*ssystem has %d equations and %d variables\n", indent, "", S->N, S->N); } if (reportSys != NULL) { reportSys(level, S); } ??? /* Get the initial guess {Z} for the heights: */ double *z = talloc(S->N, double); if (verbose) { fprintf(stderr, "%*staking initial guess of heights ...\n", indent, ""); } pst_imgsys_copy_image_to_sol_vec(S, Z, z, 0.0); auto void reportSol(int32_t level, int32_t iter, double change, bool_t final, uint32_t N, double z[]); if (verbose) { fprintf(stderr, "%*sSolving the system ...\n", indent, ""); } bool_t para = FALSE; /* Should be parameter. 1 means parallel execution, 0 sequential. */ bool_t szero = TRUE; /* Should be parameter. 1 means adjust sum to zero, 0 let it float. */ uint32_t *ord = NULL; if (topoSort) { ord = pst_imgsys_sort_equations(S); } pst_imgsys_solve_iterative ( S, z, ord, maxIter, convTol, para, szero, verbose, level, reportStep, (reportHeights == NULL ? NULL : &reportSol) ); pst_imgsys_copy_sol_vec_to_image(S, z, Z, 0.0); free(z); pst_imgsys_free(S); if (ord != NULL) { free(ord); } void reportSol(int32_t level, int32_t iter, double change, bool_t final, uint32_t N, double z[]) { assert(reportHeights != NULL); pst_imgsys_copy_sol_vec_to_image(S, z, Z, 0.0); reportHeights(level, iter, change, final, Z); }