/* Last edited on 2025-04-26 04:50:39 by stolfi */ /* ---------------------------------------------------------------------- */ /* psgr_read.{h,c} */ fget_skip_spaces_and_match(rd, "at"); r2_t p; fget_skip_spaces_and_test_char(rd, '('); p.c[0] = fget_double(rd); fget_skip_spaces_and_test_char(rd, ','); p.c[1] = fget_double(rd); fget_skip_spaces_and_test_char(rd, ')'); demand(isfinite(p.c[0]) && isfinite(p.c[1]), "invalid vertex coordinates"); fprintf(wr, " at (%f,%f)", vd->coords.c[0], vd->coords.c[1]); int32_t *list_next = talloc(2*NE, int32_t); int32_t *list_long_bit = talloc(2*NE, int32_t); auto psgr_arc_t recover_edge(int32_t index, int32_t long_bit); psgr_arc_t recover_edge(int32_t index, int32_t long_bit) { assert((index >= 0 ) && (index < gr->NE)); psgr_arc_t e = gr->{hedge|dedge}@@[index].arc; assert(e != NULL); if (long_bit == 1) { e = psgr_arc_sym(e); } return e; } /*Hard part*/ for (int32_t i = 0; i < gr->NE;i++) { psgr_arc_t e = gr->{hedge|dedge}@@[i].arc; if (e != NULL) { psgr_arc_t e_next = recover_edge(list_next[2*i],list_long_bit[2*i]); psgr_arc_t e_next_sym = recover_edge(list_next[(2*i)+1],list_long_bit[(2*i)+1]); haf_splice( haf_oprev(e_next),e); haf_splice( haf_oprev(e_next_sym),psgr_arc_sym(e)); } } free(list_next); free(list_long_bit); /* ---------------------------------------------------------------------- */ /* psgr.{h,c} */ void compare_int32(int32_t gval, int32_t hval, char *name1, int32_t ix, char *name2) { if (gval != hval) { blast(name1, ix, name2); fprintf(stderr, " %d %d\n", gval, hval); ok = FALSE; } } void compare_r2(r2_t *gval, r2_t *hval, char *name1, int32_t ix, char *name2) { double d = r2_dist(gval, hval); double tol = 2.0e-6; /* Considering the format of {psgr_write}. */ if (d > tol) { blast(name1, ix, name2); r2_print(stderr, gval); fprintf(stderr, " "); r2_print(stderr, hval); fprintf(stderr, "\n"); ok = FALSE; } } psgr_edge_data_t *psgr_arc_edge_data(psgr_t *gr, psgr_arc_t ai); /* Returns a pointer to the edge data rectord {ed} of the undirected edge of the arc {ai}. */ psgr_edge_data_t *psgr_arc_edge_data(psgr_t *gr, psgr_arc_t ai) { demand(ai != psgr_NONE, "invalid null arc {ai}"); psgr_edge_t ei = psgr_arc_edge(ai); psgr_edge_data_t *ed = &(gr->edata[ei]); return ed; } auto void compare_string(char *gval, char *hval, char *name1, int32_t ix, char *name2); void compare_string(char *gval, char *hval, char *name1, int32_t ix, char *name2) { if (strcmp(gval, hval) != 0) { blast(name1, ix, name2); fprintf(stderr, " «%s» «%s»\n", gval, hval); ok = FALSE; } } void psgr_edge_free(psgr_t *gr, psgr_edge_t e); /* Reclaims the heap data associated with edge {e}, including its label and path, and the {had_edge_rec_t} pointed to by {e}. */ void psgr_edge_free(psgr_t *gr, haf_edge_t e) { if (e == NULL) { return; } psgr_edge_t ei = psgr_arc_edge(a); assert(e == gr->hedge[ei]); /* Trash contents of edge data record: */ psgr_edge_data_t *ed = &(gr->edata[ei]); ed->org[0] = UINT32_MAX; ed->org[1] = UINT32_MAX; ed->weight = 0; ed->delta = 0; free(ed->label); ed->label = NULL; psgr_path_t *p = &(ed->path); if (p->v != NULL) { free(p->v); p->v = NULL; } haf_edge_free(e); } psgr_arc_t psgr_get_arc_id(psgr_arc_t a); /* Returns the ID of the arc {a}, which is twice the ID of the undirected edge {e} plus the longitudinal orientation bit of {a}. If {a} is {NULL}, returns {-1}. */ void psgr_set_arc_weight(psgr_t *gr, psgr_arc_t a, double w); /* Sets the weight of the (undirected) edge underlying arc {a} to {w}. */ void psgr_set_arc_weight(psgr_t *gr, psgr_arc_t a, double w) { psgr_edge_t ei = psgr_edge(a); demand(ei < gr->NE, "invalid edge ID"); assert((w == 0) || (w > 1.0e-100)); psgr_edge_data_t *ed = &(gr->edata[ei]); ed->weight = w; } void psgr_set_arc_delta(psgr_t *gr, psgr_arc_t a, double delta); /* Sets {delta} as the supposed height difference along the (directed) arc {a}. Also sets {-delta} as the difference along {psgr_sym(a)}. */ void psgr_set_arc_delta(psgr_t *gr, psgr_arc_t a, double delta) { psgr_edge_t ei = psgr_edge(a); psgr_dir_t db = psgr_dir(a); demand(ei < gr->NE, "invalid edge ID"); psgr_edge_data_t *ed = &(gr->edata[ei]); ed->delta = (db == 0 ? +1.0 : -1.0)*delta; } void psgr_set_edge_path(psgr_t *gr, psgr_arc_t a, psgr_path_t p) { psgr_edge_t ei = psgr_edge(a); psgr_dir_t db = psgr_dir(a); demand(ei < gr->NE, "invalid edge ID"); psgr_edge_data_t *ed = &(gr->edata[ei]); ed->path = (db == 0 ? p : psgr_path_reverse(p)); } void psgr_set_edge_path(psgr_t *gr,psgr_arc_t a, psgr_path_t p); /* Sets {p} as the plotting path associated with the directed arc {a}. Also creates the reversed path and associates it with the arc {psgr_sym(a)}. */ psgr_path_t psgr_path_create_single(r2_t coords); psgr_path_t psgr_path_create_single(r2_t coords) { r2_t *v = (r2_t*)malloc(sizeof(r2_t)); v[0] = coords; psgr_path_t p = (psgr_path_t){ .n=1, .v=v, .reverse=FALSE }; return p; } void pst_img_graph_vertex_remove ( pst_img_graph_t* gr, uint32_t vi, double* w_i, double wmag, bool_t verbose ); /* ---------------------------------------------------------------------- */ /* pst_img_graph_from_slope_map.{h,c} */ /* The west-pointing arc out of vertex {x} in prev/current row is {west[x]}: */ psgr_arc_t *west = talloc(NX_Z, psgr_arc_t); /* The south-pointing arc out of vertex {x} in prev/current row is {south[x]}: */ psgr_arc_t *south = talloc(NX_Z, psgr_arc_t); for (int32_t x = 0; x < NX_Z; x++) { west[x] = south[x] = NONE; } west[x] = ah; south[x] = av; haf_arc_t a = insert_edge(gr, 0, 0, 1, 0); haf_arc_t edge_y = a; for (int32_t x = 1; x < NX; x++ ) { haf_arc_t e = insert_edge(gr, x, 0, x+1, 0); haf_splice(haf_sym(edge_y),e); edge_y = e; } haf_arc_t e = insert_edge(gr,0,y,0,y+1); haf_splice(a,e); haf_arc_t b = e; haf_arc_t c = a; haf_arc_t ev = pst_img_graph_insert_edge_from_gradient_and_weright_maps(gr,IG,IW,x+1,y+0,Y_AXIS,+1); // haf_arc_t eh = pst_img_graph_insert_edge_from_gradient_and_weright_maps(gr,IG,IW,x+1,y+1,X_AXIS,-1); haf_arc_t eh = insert_edge(gr,x+1,y+1,x,y+1); haf_arc_t t = haf_lnext(c); haf_splice(t,ev); haf_splice(haf_sym(ev),eh); haf_splice(haf_sym(b), haf_sym(eh)); if(o->addDiags) { // haf_arc_t ed = pst_img_graph_insert_diagonal_edge_from_gradient_and_weright_maps(gr,IG,IW,x+1,y,-1,+1); haf_arc_t ed = insert_edge(gr,x+1,y,x,y+1); haf_splice(ed,ev); haf_splice(haf_sym(ed),haf_sym(b)); } c = t; b = ev; } if (!o->addDiags) { a = haf_sym(haf_lprev(haf_lprev(a))); }else{ a = haf_onext(haf_sym(haf_lnext(a))); } } /*now we have to clean the table from the null edges*/ int32_t valid_edges = 0; int32_t i; for (int32_t i = 0; i < gr->m; i++) { if(gr->edge[i].data->weight > 0) { valid_edges++;} } if( valid_edges != gr->m) { fprintf(stderr,"Fixing...\n"); int32_t count_valid = 0; for (int32_t i = 0; i < gr->m; i++) { if(gr->edge[i].data->weight > 0) { count_valid++; }else{ pst_img_graph_edge_remove(gr,gr->edge[i].edge); } } } /* ---------------------------------------------------------------------- */ /* pst_img_graph_from_maps.c */ void get_edge_parameters(int32_t xu, int32_t yu, int32_t xv, int32_t yv, double *d, double *w) { int32_t dx = xv - xu; int32_t dy = yv - yu; if (dx == 0) { assert(abs(dy) == 1); pst_img_graph_get_axial_edge_data_from_maps(IG,IW, xu,yu, Y_AXIS, dy, d,w); } else if (dy == 0) { assert(abs(dx) == 1); pst_img_graph_get_axial_edge_data_from_maps(IG,IW, xu,yu, X_AXIS, dx, d,w); } else { assert(abs(dx) == 1); assert(abs(dy) == 1); pst_img_graph_get_diagonal_edge_data_from_maps(IG,IW, xu,yu, dx,dy, d,w); } /* Is this possible? else{ *d = 0; *w = 1/r2_dist_sqr(&u,&v); } */ return; void get_nearest_grid_corner(r2_t pt, int32_t NX, int32_t NY, int32_t* x, int32_t* y) { assert( (pt.c[0] >= 0) && (pt.c[0] <= NX )); assert( (pt.c[1] >= 0) && (pt.c[1] <= NY )); *x = pt.c[0]; *y = pt.c[1]; } } void pst_img_graph_compute_edge_parameters ( int32_t NX, int32_t NY, int32_t gradientFunction, int32_t weightFunction, r2_t u, r2_t v, double *d, double *w ) { double ww; r2_t p; r2_mix(0.5,&u,0.5,&v,&p); switch(weightFunction) { case WGHT_FUNC_NONE : break; case WGHT_FUNC_CONST : *w = 2.5; break; case WGHT_FUNC_RAMP : *w = (0.25/NX)*p.c[0] + (0.75/NY)*p.c[1] ; break; case WGHT_FUNC_RANDOM : ww = (sin(17*p.c[0]*p.c[0]) + cos(29*p.c[1]*p.c[1]))*47; (*w) = ww - floor(ww); break; default: demand(FALSE,"Invalid Weight function"); }; r2_t uv; r2_sub(&v,&u,&uv); double dx,dy,dd; switch(gradientFunction) { case GRAD_FUNC_NONE : break; case GRAD_FUNC_CONST : *d = 2.5; break; case GRAD_FUNC_RAMP : dx = 2*((0.25/NX)*p.c[0] + (0.75/NY)*p.c[1]) -1 ; dy = 2*((0.75/NX)*p.c[0] + (0.25/NY)*p.c[1]) -1 ; *d = dx*uv.c[0] + dy*uv.c[1]; break; case GRAD_FUNC_RANDOM : dd = (sin(17*p.c[0]*p.c[0]) + cos(29*p.c[1]*p.c[1]))*47; dd = 2*(dd - floor(dd)) -1; dd *= r2_dist(&u,&v); (*d)= dd; break; default: demand(FALSE,"Invalid Grad function"); }; } void pst_img_graph_get_diagonal_edge_data_from_maps( float_image_t* IG, float_image_t* IW, int32_t x, int32_t y, int32_t dirx,int32_t diry, double *d, double*w ) { if ((IG == NULL) && (IW == NULL) ) { *d = 0; *w = 0.5; return; } assert( (abs(dirx) == 1) && (abs(diry) == 1) ); int32_t NX = (IG == NULL ? IW->sz[1] : IG->sz[1]); int32_t NY = (IG == NULL ? IW->sz[2] : IG->sz[2]); int32_t ix = (dirx > 0 ? x:x-1); int32_t iy = (diry > 0 ? y:y-1); assert( (ix >=0) && (ix < NX)); assert( (iy >=0) && (iy < NY)); double dx = dirx*(IG == NULL ? 0 : float_image_get_sample(IG,0,ix,iy)); double dy = diry*(IG == NULL ? 0 : float_image_get_sample(IG,1,ix,iy)); *w = 0.5*(IW == NULL ? 1.0 : float_image_get_sample(IW,0,ix,iy)); *d = dx+dy; if (*w == 0) { *d = 0; } } /* ---------------------------------------------------------------------- */ /* pst_img_graph.{h,c} */ int32_t pst_img_graph_get_vertex_index_from_image_indices(int32_t ix, int32_t iy, int32_t NX, int32_t NY); void pst_img_graph_get_vertex_image_indices(r2_t *p ,int32_t NX, int32_t NY, int32_t *ix, int32_t *iy); void pst_img_graph_connect_vertices(pst_img_graph_t* gr); void pst_img_graph_reorganize(pst_img_graph_t* gr); void pst_img_graph_reorganize(pst_img_graph_t *gr) { demand(FALSE, "pst_img_graph_reorganize not implemented"); } /* ---------------------------------------------------------------------- */ /* pst_img_graph_integrate.{h,c} */ void pst_img_graph_solve_system ( pst_img_graph_t *gr, pst_imgsys_t *S, double *iZ, double *iW, int32_t *ref_tab, uint32_t maxIter, double convTol, bool_t para, bool_t szero, bool_t verbose ); uint32_t *pst_img_graph_sort_equations ( pst_imgsys_t *S, int32_t *ref_tab, double *iW , uint32_t NV ); uint32_t *pst_img_graph_sort_equations ( pst_imgsys_t *S, int32_t *ref_tab, double *iW , uint32_t NV ); void pst_img_graph_integration ( pst_img_graph_t *gr, double *iZ, double *iW, uint32_t maxIter, double convTol, bool_t para, bool_t szero, bool_t verbose, uint32_t level, float_image_t *OZ, /*debug only*/ float_image_t *RZ, char* out_prefix ); void pst_img_graph_put_solution_into_image ( pst_img_graph_t* gr, double *iZ, float_image_t* OZ ); void pst_img_graph_put_error_into_image ( pst_img_graph_t *gr, double *iZ, double *iW, float_image_t *RZ, float_image_t *OZ ); void pst_img_graph_put_solution_into_image ( pst_img_graph_t* gr, double *iZ, float_image_t* OZ ); void pst_img_graph_put_error_into_image ( pst_img_graph_t *gr, double *iZ, double *iW, float_image_t *RZ, float_image_t *OZ ); void pst_img_graph_solve_system ( pst_img_graph_t *gr, pst_imgsys_t *S, double *iZ, double *iW, int32_t *ref_tab, uint32_t maxIter, double convTol, bool_t para, bool_t szero, bool_t verbose ) { if ( gr->NV_valid == S->N ) { fprintf(stderr,"IGUAL\n"); } else { fprintf(stderr,"NAO IGUAL\n"); } double *Z = rn_alloc(S->N); int32_t count_vt = 0; for (uint32_t i = 0; i < gr->NV; i++) { if (ref_tab[i] >= 0) { Z[count_vt] = iZ[i]; count_vt++; } } assert(count_vt == S->N); uint32_t *queue = pst_img_graph_sort_equations(S,ref_tab,iW,gr->NV ); pst_imgsys_solve_iterative(S, Z, queue, maxIter, convTol, para, szero, verbose, 0, NULL); count_vt = 0; for (uint32_t i = 0; i < gr->NV; i++) { if ( ref_tab[i] >= 0) { iZ[i] = Z[count_vt]; count_vt++; } } assert(count_vt == S->N); free(Z); /* free(queue); */ } void pst_img_graph_put_solution_into_image(pst_img_graph_t *gr, double *iZ, float_image_t *OZ) { float_image_fill_channel(OZ, 0, 0.0); for (uint32_t i = 0; i < gr->NV; i++) { pst_img_graph_vertex_data_t *v = &(gr->vertex[i]); if (v->id == -1) continue; double z = iZ[i]; int32_t x,y; pst_img_graph_get_vertex_image_indices(&(v->coords), (int32_t)OZ->sz[1], (int32_t)OZ->sz[2], &x, &y); float_image_set_sample(OZ, 0, x, y, (float)z); } } void pst_img_graph_put_error_into_image(pst_img_graph_t *gr, double *iZ, double*iW, float_image_t *RZ, float_image_t *OZ) { float_image_fill_channel(OZ,0,0.0); double sumWD= 0; double sumW = 1.0e-300; for (uint32_t i = 0; i < gr->NV; i++) { pst_img_graph_vertex_data_t *v = &(gr->vertex[i]); if (v->id == -1) continue; double iz = iZ[i]; double iw = iW[i]; if (iw == 0) continue; int32_t x,y; pst_img_graph_get_vertex_image_indices(&(v->coords), (int32_t)OZ->sz[1], (int32_t)OZ->sz[2], &x, &y); double rz = float_image_get_sample(RZ, 0, x, y); double dz = iz - rz; sumW+= iw; sumWD+= iw*dz; } double avgdz = sumWD/sumW; for (uint32_t i = 0; i < gr->NV; i++) { pst_img_graph_vertex_data_t *v = &(gr->vertex[i]); if (v->id == -1) continue; double iw = iW[i]; if (iw == 0) continue; double iz = iZ[i]; int32_t x,y; pst_img_graph_get_vertex_image_indices(&(v->coords), (int32_t)OZ->sz[1], (int32_t)OZ->sz[2], &x,&y); double rz = float_image_get_sample(RZ,0,x,y); double dz = iz - rz; float_image_set_sample(OZ, 0, x, y, (float)(dz - avgdz)); } } uint32_t *pst_img_graph_sort_equations ( pst_imgsys_t *S, int32_t *ref_tab, double *iW , uint32_t NV ) { #define MAXDEGREE (2*((MAX_COEFFS)-1)) uint32_t N = S->N; assert(MAX_COEFFS <= 20); /* Build the graph of the equations */ uint32_t *graph = talloc(MAXDEGREE*N, uint32_t); /* Value of out_degree[k] is the number of neighbours of pixel[k] with weight less than OW[k] */ uint32_t *out_degree = talloc(N, uint32_t); /* Value of in_degree[k] is the number of unprocessed neighbours of pixel[k] with weight greater than OW[k] */ uint32_t *in_degree = talloc(N, uint32_t); int32_t *inv_ref_tab = (int32_t*)malloc(sizeof(int32_t)*S->N); int32_t count_vt = 0; for (int32_t k = 0; k < NV; k++) { if (ref_tab[k] != -1) { inv_ref_tab[count_vt] = k; count_vt++; } } assert(count_vt == S->N); auto bool_t arrow( uint32_t k1, uint32_t k2); bool_t arrow( uint32_t k1, uint32_t k2) { int32_t refk1 = inv_ref_tab[k1]; int32_t refk2 = inv_ref_tab[k2]; double w1 = iW[refk1]; double w2 = iW[refk2]; return w1 < w2; } for (uint32_t k = 0; k < N; k++) { in_degree[k] = 0; out_degree[k] = 0; } for (uint32_t k = 0; k < N; k++) { pst_imgsys_equation_t *eqk = &(S->eq[k]); assert( eqk->uid[0] == k); for (uint32_t j = 1; j < eqk->nt; j++) { uint32_t i = eqk->uid[j]; if (arrow(k,i)) { in_degree[i]++; graph[MAXDEGREE*k + out_degree[k]] = i; out_degree[k]++;} if (arrow(i,k)) { in_degree[k]++; graph[MAXDEGREE*i + out_degree[i]] = k; out_degree[i]++;} } } uint32_t queue_free = 0; uint32_t queue_start = 0; uint32_t *queue = talloc(N, uint32_t); auto void eq_queue_insert(uint32_t index); void eq_queue_insert(uint32_t index) { assert(queue_free < N); queue[queue_free] = index; queue_free++; } auto uint32_t eq_queue_remove( void ); uint32_t eq_queue_remove( void ) { assert(queue_start < queue_free ); uint32_t i = queue[queue_start]; queue_start++; return i; } for (uint32_t k = 0; k < N; k++) {if ( in_degree[k] == 0) { eq_queue_insert(k); }} while(queue_start < queue_free ) { uint32_t k = eq_queue_remove(); for (uint32_t j = 0; j < out_degree[k]; j++) { uint32_t i = graph[MAXDEGREE*k + j]; assert(in_degree[i] > 0); in_degree[i]--; if (in_degree[i] == 0) { eq_queue_insert(i); } } } assert(queue_free == N); free(graph); free(in_degree); free(out_degree); free(inv_ref_tab); return queue; } /* ---------------------------------------------------------------------- */ /* pst_img_graph_integrate_recursive.{h,c} */ void pst_img_graph_copy_solution_from_shrunk(pst_img_graph_t *jg, double *jZ, pst_img_graph_t *ig, double *iZ) { assert(ig->NV >= jg->NV); uint32_t kvj = 0; for (uint32_t kvi = 0; kvi < ig->NV; kvi++) { pst_img_graph_vertex_data_t *vdi = &(ig->vdata[kvi]); if ( vdi->vmark == DELETED) { iZ[kvi] = 0; continue; } pst_img_graph_vertex_data_t *vdj = NULL; while(kvj < jg->NV) { vdj = &(jg->vdata[kvj]); if ( vdj->vmark != DELETED) break; kvj++; } if ((kvj >= jg->NV) || (vdj->id > vdi->id) ) { iZ[kvi] = 0; vdi->mark = pst_img_graph_mark_REMOVED; } else if (vdj->id == vdi->id) { iZ[kvi] = jZ[kvj]; vdi->mark = pst_img_graph_mark_TO_KEEP; kvj++; }else{ demand(FALSE,"Vertices out of order!"); } } } ] if (debug) { pst_img_graph_put_curl_into_image(gr, OZ); char *filename = jsprintf("%s-%02d-CL.fni",out_prefix,level); FILE *wr_curl = open_write(filename,FALSE); float_image_write(wr_curl,OZ); fclose(wr_curl); free(filename); } if (debug) { pst_img_graph_mark_vertex_removal(gr,NULL); char *filename = jsprintf("%s-%02d-GR.txt",out_prefix,level); FILE *wr_grph = open_write(filename,FALSE); pst_img_graph_write(wr_grph,gr); fclose(wr_grph); free(filename); filename = jsprintf("%s-%02d-GR.grf",out_prefix,level); wr_grph = open_write(filename,FALSE); pst_img_graph_write(wr_grph,gr); fclose(wr_grph); free(filename); } if (debug) { char *filename = jsprintf("%s-%02d-S.txt",out_prefix,level); FILE *wr_dump = open_write(filename,FALSE); pst_imgsys_write(wr_dump,S); fclose(wr_dump); free(filename); pst_img_graph_put_solution_into_image(gr,iZ,OZ); filename = jsprintf("%s-%02d-iZ.fni",out_prefix,level); wr_dump = open_write(filename,FALSE); float_image_write(wr_dump,OZ); fclose(wr_dump); free(filename); pst_img_graph_put_solution_into_image(gr,iW,OZ); filename = jsprintf("%s-%02d-iW.fni",out_prefix,level); wr_dump = open_write(filename,FALSE); float_image_write(wr_dump,OZ); fclose(wr_dump); free(filename); if (RZ != NULL) { pst_img_graph_put_error_into_image(gr,iZ,iW,RZ,OZ); filename = NULL; char *filename = jsprintf("%s-%02d-iE.fni",out_prefix,level); wr_dump = open_write(filename,FALSE); float_image_write(wr_dump,OZ); fclose(wr_dump); free(filename); } } if (debug) { if (RZ != NULL) { pst_img_graph_put_error_into_image(gr,iZ,iW,RZ,OZ); char *filename = jsprintf("%s-%02d-oE.fni",out_prefix,level); FILE *wr_dump = open_write(filename,FALSE); float_image_write(wr_dump,OZ); fclose(wr_dump); free(filename); } } pst_img_graph_put_solution_into_image(gr,iZ,OZ); if (debug) { char *filename = jsprintf("%s-%02d-oZ.fni",out_prefix,level); FILE *wr_dump = open_write(filename,FALSE); float_image_write(wr_dump,OZ); fclose(wr_dump); free(filename); } /* ---------------------------------------------------------------------- */ /* compute_rms.c */ float_image_t* height_map_shrink_one_pixel(float_image_t* Z); float_image_t* height_map_shrink_one_pixel(float_image_t* Z){ int32_t NC = Z->sz[0]; int32_t NX = Z->sz[1]; int32_t NY = Z->sz[2]; float_image_t* IJ = float_image_new(NC,NX-1,NY-1); int32_t c; for( c = 0; c < NC; c++) { int32_t x; for(x = 0; x < NX -1; x++) { int32_t y; for(y = 0; y < NY -1; y++) { double P00 = float_image_get_sample(Z,c,x,y); double P10 = float_image_get_sample(Z,c,x+1,y); double P01 = float_image_get_sample(Z,c,x,y+1); double P11 = float_image_get_sample(Z,c,x+1,y+1); double val = (P00 + P01 + P10 + P11)/4.0; float_image_set_sample(IJ,c,x,y,val); } } } return IJ; } /* ---------------------------------------------------------------------- */ /* psgr_shrink.{h,c} */ /* Disabled additional term for degree 5: */ double B = 0.4425; weight += B*(w_star[k0]*w_star[k3] + w_star[k1]*w_star[k3] ); void psgr_shrink_stacy_trans_gen ( psgr_t* gr, uint32_t n, double d_star[], double w_star[], double wtot, double d_cycle[], double w_cycle[] ) { double lambda = 0.5; auto void compute_wd2l(double w0, double w1, double w2,double *w2l); for (int32_t ke = 0; ke < n; ke++) { int32_t kc = (ke - 2 + n) % n; int32_t kd = (ke - 1 + n) % n; int32_t kf = (ke + 1) % n; int32_t kg = (ke + 2) % n; int32_t kh = (ke + 3) % n; double w0_org = lambda*(w_star[ke]*w_star[kc])/wtot; double w0_dst = lambda*(w_star[kf]*w_star[kg])/wtot; uint32_t kdst = kg; for (uint32_t i = 0; i < n - 3; i++) { uint32_t korg = (kc + i) % n; uint32_t kon = (kd + i) % n; double w2_org = (w_star[ke]*w_star[kon])/wtot; double w1_org = (w_star[korg]*w_star[kon])/wtot; double w2l; compute_wd2l(w0_org,w1_org,w2_org,&w2l); w0_org = lambda*w2_org + w2l; korg = kon; uint32_t kdn = update_ind(kdst,+1); double w2_dst = (w_star[kf]*w_star[kdn])/wtot; double w1_dst = (w_star[kdst]*w_star[kdn])/wtot; compute_wd2l(w0_dst,w1_dst,w2_dst,&w2l); w0_dst = lambda*w2_dst + w2l; kdst = kdn; } w_cycle[ke] = w0_org + w0_dst; d_cycle[ke] = d_star[kf] - d_star[ke]; } return void compute_wd2l(double w0, double w1, double w2,double *w2l) { *w2l = w0*(w1+w2)/w1; } } double *d_default, double *w_default /* deg 4, 5, 6 */ *d_default = delta; *w_default = weight/wtot; /* General deg 3, 7, ... */ *w_default = w0_org + w0_dst ; *d_default = -deltas[ind_e] + deltas[ind_f]; if (psgr_arc_lnext(gr, psgr_arc_lnext(gr, ae)) != psgr_arc_sym(af)) { FILE *wr = open_write("debug.txt", TRUE); psgr_write(wr,gr); fclose(wr); } void psgr_edge_remove(psgr_t *gr, psgr_arc_t ai) { demand(e != ANONE , "invalid null edge {e}"); uint32_t eid = (uint32_t)psgr_arc_edge_id(e); assert(e == gr->hedge[eid]); psgr_edge_data_t *ed = &(gr->edata[eid]); /* Disconnect the edge from he half-edge mesh structure: */ psgr_arc_t a0 = psgr_arc_base_arc(e); psgr_arc_t a1 = psgr_arc_sym(a0); psgr_arc_t oprev_a0 = psgr_arc_oprev(gr, a0); psgr_arc_t oprev_a1 = psgr_arc_oprev(gr, a1); uint32_t org_id = ed->org[0]; uint32_t dst_id = ed->org[1]; if (a0 != oprev_a0) { psgr_arc_splice(a0, oprev_a0); } if (a1 != oprev_a1) { psgr_arc_splice(a1, oprev_a1); } assert(psgr_arc_dprev(gr, a0) == a0); assert(psgr_arc_dnext(gr, a0) == a0); assert(psgr_arc_lnext(gr, a0) == a1); assert(psgr_arc_lprev(gr, a0) == a1); gr->vdata[org_id].aout = (oprev_a0 == a0 ? ANONE : oprev_a0); gr->vdata[dst_id].aout = (oprev_a1 == a1 ? ANONE : oprev_a1); gr->hedge[eid] = NULL; psgr_edge_free(gr, e); } void psgr_vertex_remove(psgr_t *gr, uint32_t vi) { demand(vi < gr->NV, "invalid vertex index {vi}"); psgr_vertex_data_t *vd = &(gr->vdata[vi]); demand(vd->aout == ANONE , "vertex is still connected"); demand(vd->vmark != DELETED, "vertex was already deleted"); vd->vmark = DELETED; } void psgr_choose_vertices_for_removal_no_bucket(psgr_t* gr); void psgr_choose_vertices_for_removal_no_bucket(psgr_t *gr) { for (uint32_t kv = 0; kv < gr->NV; kv++) { psgr_vertex_data_t *vd = &(gr->vdata[kv]); if (kv == -1) continue; if (vd->vmark == psgr_mark_NONE) { vd->vmark = psgr_mark_REMOVED; /* fprintf(stderr,"*"); */ if ( vd->aout.ix != ANONE.ix ) { psgr_arc_t e = vd->aout; do { uint32_t dest = psgr_arc_org(gr,psgr_arc_sym(e)); if (dest == -1) { fprintf(stderr,"ALERT - vertex %d is damaged\n",kv); break; } /* fprintf(stderr,"+"); */ gr->vdata[dest].mark = psgr_mark_TO_KEEP; e = psgr_arc_onext(gr, e); } while(e != vd->aout); } } } /* fprintf(stderr,"\n"); */ } #define DEFAULT_WMAG 1.0 double* w_i, double wmag, double w = w = (w_i != NULL ? w_i[count] : w_default); double d = d_default; fprintf(stderr," WDEF = %12.6lf WMAG = %12.9lf", w_default, w/w_default); fprintf(stderr," WDEF = %12.6lf WMAG = %12.9lf", w_default, w/w_default); if (merge_diagonals) { assert(deg_vi >=4); } /* ---------------------------------------------------------------------- */ /* psgr_shrink_stacy.{h,c} (version RFVS, SAVE/2013-01-14/pst_img_graph.c) */ void pst_img_graph_compute_cycle_delta_weight( double wtot, long int n_neighbours, double* deltas, double* weights, long int ind_e, double *d_cycle, double *w_cycle ){ double lambda = 0.5; auto void compute_wd2l(double w0, double w1, double w2,double* w2l); auto void compute_wd2l(double w0, double w1, double w2,double* w2l){ *w2l = w0*(w1+w2)/w1; } auto long int update_ind(long int i,long int increase); long int update_ind(long int i,long int increase){ return (i + n_neighbours + increase)%n_neighbours; } long int ind_f = update_ind(ind_e,+1); long int ind_org_act = update_ind(ind_e,-2); double w0_org = lambda*(weights[ind_e]*weights[ind_org_act])/wtot; long int ind_dst_act = update_ind(ind_f,+2); double w0_dst = lambda*(weights[ind_f]*weights[ind_dst_act])/wtot; long int i; for(i = 0; i < n_neighbours - 3; i++){ long int ind_org_next = update_ind(ind_org_act,-1); double w2_org = (weights[ind_e]*weights[ind_org_next])/wtot; double w1_org = (weights[ind_org_act]*weights[ind_org_next])/wtot; double w2l; compute_wd2l(w0_org,w1_org,w2_org,&w2l); w0_org = lambda*w2_org + w2l; ind_org_act = ind_org_next; long int ind_dst_next = update_ind(ind_dst_act,+1); double w2_dst = (weights[ind_f]*weights[ind_dst_next])/wtot; double w1_dst = (weights[ind_dst_act]*weights[ind_dst_next])/wtot; compute_wd2l(w0_dst,w1_dst,w2_dst,&w2l); w0_dst = lambda*w2_dst + w2l; ind_dst_act = ind_dst_next; } *w_cycle = w0_org + w0_dst ; *d_cycle = -deltas[ind_e] + deltas[ind_f]; } /* ---------------------------------------------------------------------- */ /* psgr_plot.{h,c} */ void psgr_plot_path(epswr_figure_t *eps, double pixelSize, i2_t *ipo, psgr_path_t P, i2_t *ipd) { bool_t debug = FALSE; if (debug) { fprintf(stderr, "P.reverse = %c\n", "FT"[P.reverse]); } auto void debug_point(int32_t j, r2_t *pj); if (debug) { debug_point(-1, &p); } if (debug) { debug_point(k, &q); } if (debug) { fprintf(stderr, "\n"); } void debug_point(int32_t j, r2_t *pj) { fprintf(stderr, "p[%d] (unscaled) = ", j+1); r2_print(stderr, pj); double d = r2_dist(&po, pj); fprintf(stderr, " d = %12.8f\n", d); } } /* ---------------------------------------------------------------------- */ /* psgr_path.{h,c} */ psgr_path_t psgr_path_throw(i2_t *ip0, uint32_t n, i2_t *ip1) { bool_t debug = FALSE; auto void debug_point(int32_t j, r2_t *pj); if (debug) { fprintf(stderr, "P.reverse = %c\n", "FT"[P.reverse]); } if (debug) { debug_point(-1, &p0); } if (debug) { debug_point(k, &pk); } if (debug) { debug_point((int32_t)P.n, &p1); } if (debug) { fprintf(stderr, "\n"); } void debug_point(int32_t j, r2_t *pj) { fprintf(stderr, "p[%d] = ", j+1); r2_print(stderr, pj); double d = r2_dist(&p0, pj); fprintf(stderr, " d = %12.8f\n", d); }