From c2b2bb16612fd68d05ac50665ce074fcaa67bd77 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 6 Nov 2018 17:00:00 +0100 Subject: [PATCH 01/53] fix pip-grid bug of double* and pPipoint --- thirdparty/ptinpoly.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/thirdparty/ptinpoly.c b/thirdparty/ptinpoly.c index 7c833924..383e7dc2 100644 --- a/thirdparty/ptinpoly.c +++ b/thirdparty/ptinpoly.c @@ -138,9 +138,11 @@ int y_flag, io_state ; } /* loop through edges and insert into grid structure */ - vtx0 = pgon[numverts-1] ; + double pt[2] = { pgon[numverts - 1]->x,pgon[numverts - 1]->y } ; + vtx0 = &pt[0] ; for ( i = 0 ; i < numverts ; i++ ) { - vtx1 = pgon[i] ; + double pt1[2] = { pgon[i]->x, pgon[i]->y } ; + vtx1 = &pt1[0] ; if ( vtx0[Y] < vtx1[Y] ) { vtxa = vtx0 ; From 39ff0788d8a72d969adc78c10d91d616016d5c4e Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 6 Nov 2018 17:03:25 +0100 Subject: [PATCH 02/53] pretty formatting of ptinpoly code --- thirdparty/ptinpoly.c | 1160 +++++++++++++++++++++-------------------- thirdparty/ptinpoly.h | 76 +-- 2 files changed, 629 insertions(+), 607 deletions(-) diff --git a/thirdparty/ptinpoly.c b/thirdparty/ptinpoly.c index 383e7dc2..c2e097fe 100644 --- a/thirdparty/ptinpoly.c +++ b/thirdparty/ptinpoly.c @@ -3,7 +3,7 @@ by Eric Haines, 3D/Eye Inc, erich@eye.com This code contains the following algorithms: - grid testing - grid imposed on polygon + grid testing - grid imposed on polygon */ #include @@ -47,389 +47,406 @@ * Call cleanup with pointer to grid structure to free space. */ -/* Strategy for setup: - * Get bounds of polygon, allocate grid. - * "Walk" each edge of the polygon and note which edges have been crossed - * and what cells are entered (points on a grid edge are always considered - * to be above that edge). Keep a record of the edges overlapping a cell. - * For cells with edges, determine if any cell border has no edges passing - * through it and so can be used for shooting a test ray. - * Keep track of the parity of the x (horizontal) grid cell borders for - * use in determining whether the grid corners are inside or outside. - */ -void GridSetup(pPipoint pgon[], int numverts, int resolution, pGridSet p_gs) -{ -double *vtx0, *vtx1, *vtxa, *vtxb, *p_gl ; -int i, j, gc_clear_flags ; -double vx0, vx1, vy0, vy1, gxdiff, gydiff, eps ; -pGridCell p_gc, p_ngc ; -double xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty ; -double tgcx, tgcy ; -int gcx, gcy, sign_x ; -int y_flag, io_state ; - - p_gs->xres = p_gs->yres = resolution ; - p_gs->tot_cells = p_gs->xres * p_gs->yres ; - p_gs->glx = (double *)malloc( (p_gs->xres+1) * sizeof(double)); - MALLOC_CHECK( p_gs->glx ) ; - p_gs->gly = (double *)malloc( (p_gs->yres+1) * sizeof(double)); - MALLOC_CHECK( p_gs->gly ) ; - p_gs->gc = (pGridCell)malloc( p_gs->tot_cells * sizeof(GridCell)); - MALLOC_CHECK( p_gs->gc ) ; - - p_gs->minx = - p_gs->maxx = pgon[0]->x ; - p_gs->miny = - p_gs->maxy = pgon[0]->y ; - - /* find bounds of polygon */ - for ( i = 1 ; i < numverts ; i++ ) { - vx0 = pgon[i]->x ; - if ( p_gs->minx > vx0 ) { - p_gs->minx = vx0 ; - } else if ( p_gs->maxx < vx0 ) { - p_gs->maxx = vx0 ; - } - - vy0 = pgon[i]->y ; - if ( p_gs->miny > vy0 ) { - p_gs->miny = vy0 ; - } else if ( p_gs->maxy < vy0 ) { - p_gs->maxy = vy0 ; - } + /* Strategy for setup: + * Get bounds of polygon, allocate grid. + * "Walk" each edge of the polygon and note which edges have been crossed + * and what cells are entered (points on a grid edge are always considered + * to be above that edge). Keep a record of the edges overlapping a cell. + * For cells with edges, determine if any cell border has no edges passing + * through it and so can be used for shooting a test ray. + * Keep track of the parity of the x (horizontal) grid cell borders for + * use in determining whether the grid corners are inside or outside. + */ +void GridSetup(pPipoint pgon[], int numverts, int resolution, pGridSet p_gs) { + double *vtx0, *vtx1, *vtxa, *vtxb, *p_gl; + int i, j, gc_clear_flags; + double vx0, vx1, vy0, vy1, gxdiff, gydiff, eps; + pGridCell p_gc, p_ngc; + double xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty; + double tgcx, tgcy; + int gcx, gcy, sign_x; + int y_flag, io_state; + + p_gs->xres = p_gs->yres = resolution; + p_gs->tot_cells = p_gs->xres * p_gs->yres; + p_gs->glx = (double *)malloc((p_gs->xres + 1) * sizeof(double)); + MALLOC_CHECK(p_gs->glx); + p_gs->gly = (double *)malloc((p_gs->yres + 1) * sizeof(double)); + MALLOC_CHECK(p_gs->gly); + p_gs->gc = (pGridCell)malloc(p_gs->tot_cells * sizeof(GridCell)); + MALLOC_CHECK(p_gs->gc); + + p_gs->minx = p_gs->maxx = pgon[0]->x; + p_gs->miny = p_gs->maxy = pgon[0]->y; + + /* find bounds of polygon */ + for (i = 1; i < numverts; i++) { + vx0 = pgon[i]->x; + if (p_gs->minx > vx0) { + p_gs->minx = vx0; } - - /* add a little to the bounds to ensure everything falls inside area */ - gxdiff = p_gs->maxx - p_gs->minx ; - gydiff = p_gs->maxy - p_gs->miny ; - p_gs->minx -= EPSILON * gxdiff ; - p_gs->maxx += EPSILON * gxdiff ; - p_gs->miny -= EPSILON * gydiff ; - p_gs->maxy += EPSILON * gydiff ; - - /* avoid roundoff problems near corners by not getting too close to them */ - eps = 1e-9 * ( gxdiff + gydiff ) ; - - /* use the new bounds to compute cell widths */ - TryAgain: - p_gs->xdelta = - (p_gs->maxx-p_gs->minx) / (double)p_gs->xres ; - p_gs->inv_xdelta = 1.0 / p_gs->xdelta ; - - p_gs->ydelta = - (p_gs->maxy-p_gs->miny) / (double)p_gs->yres ; - p_gs->inv_ydelta = 1.0 / p_gs->ydelta ; - - for ( i = 0, p_gl = p_gs->glx ; i < p_gs->xres ; i++ ) { - *p_gl++ = p_gs->minx + i * p_gs->xdelta ; + else if (p_gs->maxx < vx0) { + p_gs->maxx = vx0; } - /* make last grid corner precisely correct */ - *p_gl = p_gs->maxx ; - for ( i = 0, p_gl = p_gs->gly ; i < p_gs->yres ; i++ ) { - *p_gl++ = p_gs->miny + i * p_gs->ydelta ; + vy0 = pgon[i]->y; + if (p_gs->miny > vy0) { + p_gs->miny = vy0; + } + else if (p_gs->maxy < vy0) { + p_gs->maxy = vy0; } - *p_gl = p_gs->maxy ; + } + + /* add a little to the bounds to ensure everything falls inside area */ + gxdiff = p_gs->maxx - p_gs->minx; + gydiff = p_gs->maxy - p_gs->miny; + p_gs->minx -= EPSILON * gxdiff; + p_gs->maxx += EPSILON * gxdiff; + p_gs->miny -= EPSILON * gydiff; + p_gs->maxy += EPSILON * gydiff; + + /* avoid roundoff problems near corners by not getting too close to them */ + eps = 1e-9 * (gxdiff + gydiff); + + /* use the new bounds to compute cell widths */ +TryAgain: + p_gs->xdelta = + (p_gs->maxx - p_gs->minx) / (double)p_gs->xres; + p_gs->inv_xdelta = 1.0 / p_gs->xdelta; + + p_gs->ydelta = + (p_gs->maxy - p_gs->miny) / (double)p_gs->yres; + p_gs->inv_ydelta = 1.0 / p_gs->ydelta; + + for (i = 0, p_gl = p_gs->glx; i < p_gs->xres; i++) { + *p_gl++ = p_gs->minx + i * p_gs->xdelta; + } + /* make last grid corner precisely correct */ + *p_gl = p_gs->maxx; + + for (i = 0, p_gl = p_gs->gly; i < p_gs->yres; i++) { + *p_gl++ = p_gs->miny + i * p_gs->ydelta; + } + *p_gl = p_gs->maxy; + + for (i = 0, p_gc = p_gs->gc; i < p_gs->tot_cells; i++, p_gc++) { + p_gc->tot_edges = 0; + p_gc->gc_flags = 0x0; + p_gc->gr = NULL; + } + + /* loop through edges and insert into grid structure */ + double pt[2] = { pgon[numverts - 1]->x,pgon[numverts - 1]->y }; + vtx0 = &pt[0]; + for (i = 0; i < numverts; i++) { + double pt1[2] = { pgon[i]->x, pgon[i]->y }; + vtx1 = &pt1[0]; + + if (vtx0[Y] < vtx1[Y]) { + vtxa = vtx0; + vtxb = vtx1; + } + else { + vtxa = vtx1; + vtxb = vtx0; + } + + /* Set x variable for the direction of the ray */ + xdiff = vtxb[X] - vtxa[X]; + ydiff = vtxb[Y] - vtxa[Y]; + tmax = sqrt(xdiff * xdiff + ydiff * ydiff); - for ( i = 0, p_gc = p_gs->gc ; i < p_gs->tot_cells ; i++, p_gc++ ) { - p_gc->tot_edges = 0 ; - p_gc->gc_flags = 0x0 ; - p_gc->gr = NULL ; + /* if edge is of 0 length, ignore it (useless edge) */ + if (tmax == 0.0) goto NextEdge; + + xdir = xdiff / tmax; + ydir = ydiff / tmax; + + gcx = (int)((vtxa[X] - p_gs->minx) * p_gs->inv_xdelta); + gcy = (int)((vtxa[Y] - p_gs->miny) * p_gs->inv_ydelta); + + /* get information about slopes of edge, etc */ + if (vtxa[X] == vtxb[X]) { + sign_x = 0; + tx = HUGE; + } + else { + inv_x = tmax / xdiff; + tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X]; + if (vtxa[X] < vtxb[X]) { + sign_x = 1; + tx += p_gs->xdelta; + tgcx = p_gs->xdelta * inv_x; + } + else { + sign_x = -1; + tgcx = -p_gs->xdelta * inv_x; + } + tx *= inv_x; } - /* loop through edges and insert into grid structure */ - double pt[2] = { pgon[numverts - 1]->x,pgon[numverts - 1]->y } ; - vtx0 = &pt[0] ; - for ( i = 0 ; i < numverts ; i++ ) { - double pt1[2] = { pgon[i]->x, pgon[i]->y } ; - vtx1 = &pt1[0] ; - - if ( vtx0[Y] < vtx1[Y] ) { - vtxa = vtx0 ; - vtxb = vtx1 ; - } else { - vtxa = vtx1 ; - vtxb = vtx0 ; - } - - /* Set x variable for the direction of the ray */ - xdiff = vtxb[X] - vtxa[X] ; - ydiff = vtxb[Y] - vtxa[Y] ; - tmax = sqrt( xdiff * xdiff + ydiff * ydiff ) ; - - /* if edge is of 0 length, ignore it (useless edge) */ - if ( tmax == 0.0 ) goto NextEdge ; - - xdir = xdiff / tmax ; - ydir = ydiff / tmax ; - - gcx = (int)(( vtxa[X] - p_gs->minx ) * p_gs->inv_xdelta) ; - gcy = (int)(( vtxa[Y] - p_gs->miny ) * p_gs->inv_ydelta) ; - - /* get information about slopes of edge, etc */ - if ( vtxa[X] == vtxb[X] ) { - sign_x = 0 ; - tx = HUGE ; - } else { - inv_x = tmax / xdiff ; - tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X] ; - if ( vtxa[X] < vtxb[X] ) { - sign_x = 1 ; - tx += p_gs->xdelta ; - tgcx = p_gs->xdelta * inv_x ; - } else { - sign_x = -1 ; - tgcx = -p_gs->xdelta * inv_x ; - } - tx *= inv_x ; - } - - if ( vtxa[Y] == vtxb[Y] ) { - ty = HUGE ; - } else { - inv_y = tmax / ydiff ; - ty = (p_gs->ydelta * (double)(gcy+1) + p_gs->miny - vtxa[Y]) - * inv_y ; - tgcy = p_gs->ydelta * inv_y ; - } - - p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ; - - vx0 = vtxa[X] ; - vy0 = vtxa[Y] ; - - t_near = 0.0 ; - - do { - /* choose the next boundary, but don't move yet */ - if ( tx <= ty ) { - gcx += sign_x ; - - ty -= tx ; - t_near += tx ; - tx = tgcx ; - - /* note which edge is hit when leaving this cell */ - if ( t_near < tmax ) { - if ( sign_x > 0 ) { - p_gc->gc_flags |= GC_R_EDGE_HIT ; - vx1 = p_gs->glx[gcx] ; - } else { - p_gc->gc_flags |= GC_L_EDGE_HIT ; - vx1 = p_gs->glx[gcx+1] ; - } - - /* get new location */ - vy1 = t_near * ydir + vtxa[Y] ; - } else { - /* end of edge, so get exact value */ - vx1 = vtxb[X] ; - vy1 = vtxb[Y] ; - } - - y_flag = FALSE ; - - } else { - - gcy++ ; - - tx -= ty ; - t_near += ty ; - ty = tgcy ; - - /* note top edge is hit when leaving this cell */ - if ( t_near < tmax ) { - p_gc->gc_flags |= GC_T_EDGE_HIT ; - /* this toggles the parity bit */ - p_gc->gc_flags ^= GC_T_EDGE_PARITY ; - - /* get new location */ - vx1 = t_near * xdir + vtxa[X] ; - vy1 = p_gs->gly[gcy] ; - } else { - /* end of edge, so get exact value */ - vx1 = vtxb[X] ; - vy1 = vtxb[Y] ; - } - - y_flag = TRUE ; - } - - /* check for corner crossing, then mark the cell we're in */ - if ( !AddGridRecAlloc( p_gc, vx0, vy0, vx1, vy1, eps ) ) { - /* warning, danger - we have just crossed a corner. - * There are all kinds of topological messiness we could - * do to get around this case, but they're a headache. - * The simplest recovery is just to change the extents a bit - * and redo the meshing, so that hopefully no edges will - * perfectly cross a corner. Since it's a preprocess, we - * don't care too much about the time to do it. - */ - - /* clean out all grid records */ - for ( i = 0, p_gc = p_gs->gc - ; i < p_gs->tot_cells - ; i++, p_gc++ ) { - - if ( p_gc->gr ) { - free( p_gc->gr ) ; - } - } - - /* make the bounding box ever so slightly larger, hopefully - * changing the alignment of the corners. - */ - p_gs->minx -= EPSILON * gxdiff * 0.24 ; - p_gs->miny -= EPSILON * gydiff * 0.10 ; - - /* yes, it's the dreaded goto - run in fear for your lives! */ - goto TryAgain ; - } - - if ( t_near < tmax ) { - /* note how we're entering the next cell */ - /* TBD: could be done faster by incrementing index in the - * incrementing code, above */ - p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ; - - if ( y_flag ) { - p_gc->gc_flags |= GC_B_EDGE_HIT ; - /* this toggles the parity bit */ - p_gc->gc_flags ^= GC_B_EDGE_PARITY ; - } else { - p_gc->gc_flags |= - ( sign_x > 0 ) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT ; - } - } - - vx0 = vx1 ; - vy0 = vy1 ; - } - /* have we gone further than the end of the edge? */ - while ( t_near < tmax ) ; - - NextEdge: - vtx0 = vtx1 ; + if (vtxa[Y] == vtxb[Y]) { + ty = HUGE; + } + else { + inv_y = tmax / ydiff; + ty = (p_gs->ydelta * (double)(gcy + 1) + p_gs->miny - vtxa[Y]) + * inv_y; + tgcy = p_gs->ydelta * inv_y; } - /* the grid is all set up, now set up the inside/outside value of each - * corner. - */ - p_gc = p_gs->gc ; - p_ngc = &p_gs->gc[p_gs->xres] ; - - /* we know the bottom and top rows are all outside, so no flag is set */ - for ( i = 1; i < p_gs->yres ; i++ ) { - /* start outside */ - io_state = 0x0 ; - - for ( j = 0; j < p_gs->xres ; j++ ) { - - if ( io_state ) { - /* change cell left corners to inside */ - p_gc->gc_flags |= GC_TL_IN ; - p_ngc->gc_flags |= GC_BL_IN ; - } - - if ( p_gc->gc_flags & GC_T_EDGE_PARITY ) { - io_state = !io_state ; - } - - if ( io_state ) { - /* change cell right corners to inside */ - p_gc->gc_flags |= GC_TR_IN ; - p_ngc->gc_flags |= GC_BR_IN ; - } - - p_gc++ ; - p_ngc++ ; - } + p_gc = &p_gs->gc[gcy*p_gs->xres + gcx]; + + vx0 = vtxa[X]; + vy0 = vtxa[Y]; + + t_near = 0.0; + + do { + /* choose the next boundary, but don't move yet */ + if (tx <= ty) { + gcx += sign_x; + + ty -= tx; + t_near += tx; + tx = tgcx; + + /* note which edge is hit when leaving this cell */ + if (t_near < tmax) { + if (sign_x > 0) { + p_gc->gc_flags |= GC_R_EDGE_HIT; + vx1 = p_gs->glx[gcx]; + } + else { + p_gc->gc_flags |= GC_L_EDGE_HIT; + vx1 = p_gs->glx[gcx + 1]; + } + + /* get new location */ + vy1 = t_near * ydir + vtxa[Y]; + } + else { + /* end of edge, so get exact value */ + vx1 = vtxb[X]; + vy1 = vtxb[Y]; + } + + y_flag = FALSE; + + } + else { + + gcy++; + + tx -= ty; + t_near += ty; + ty = tgcy; + + /* note top edge is hit when leaving this cell */ + if (t_near < tmax) { + p_gc->gc_flags |= GC_T_EDGE_HIT; + /* this toggles the parity bit */ + p_gc->gc_flags ^= GC_T_EDGE_PARITY; + + /* get new location */ + vx1 = t_near * xdir + vtxa[X]; + vy1 = p_gs->gly[gcy]; + } + else { + /* end of edge, so get exact value */ + vx1 = vtxb[X]; + vy1 = vtxb[Y]; + } + + y_flag = TRUE; + } + + /* check for corner crossing, then mark the cell we're in */ + if (!AddGridRecAlloc(p_gc, vx0, vy0, vx1, vy1, eps)) { + /* warning, danger - we have just crossed a corner. + * There are all kinds of topological messiness we could + * do to get around this case, but they're a headache. + * The simplest recovery is just to change the extents a bit + * and redo the meshing, so that hopefully no edges will + * perfectly cross a corner. Since it's a preprocess, we + * don't care too much about the time to do it. + */ + + /* clean out all grid records */ + for (i = 0, p_gc = p_gs->gc + ; i < p_gs->tot_cells + ; i++, p_gc++) { + + if (p_gc->gr) { + free(p_gc->gr); + } + } + + /* make the bounding box ever so slightly larger, hopefully + * changing the alignment of the corners. + */ + p_gs->minx -= EPSILON * gxdiff * 0.24; + p_gs->miny -= EPSILON * gydiff * 0.10; + + /* yes, it's the dreaded goto - run in fear for your lives! */ + goto TryAgain; + } + + if (t_near < tmax) { + /* note how we're entering the next cell */ + /* TBD: could be done faster by incrementing index in the + * incrementing code, above */ + p_gc = &p_gs->gc[gcy*p_gs->xres + gcx]; + + if (y_flag) { + p_gc->gc_flags |= GC_B_EDGE_HIT; + /* this toggles the parity bit */ + p_gc->gc_flags ^= GC_B_EDGE_PARITY; + } + else { + p_gc->gc_flags |= + (sign_x > 0) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT; + } + } + + vx0 = vx1; + vy0 = vy1; } + /* have we gone further than the end of the edge? */ + while (t_near < tmax); + + NextEdge: + vtx0 = vtx1; + } + + /* the grid is all set up, now set up the inside/outside value of each + * corner. + */ + p_gc = p_gs->gc; + p_ngc = &p_gs->gc[p_gs->xres]; + + /* we know the bottom and top rows are all outside, so no flag is set */ + for (i = 1; i < p_gs->yres; i++) { + /* start outside */ + io_state = 0x0; + + for (j = 0; j < p_gs->xres; j++) { + + if (io_state) { + /* change cell left corners to inside */ + p_gc->gc_flags |= GC_TL_IN; + p_ngc->gc_flags |= GC_BL_IN; + } + + if (p_gc->gc_flags & GC_T_EDGE_PARITY) { + io_state = !io_state; + } + + if (io_state) { + /* change cell right corners to inside */ + p_gc->gc_flags |= GC_TR_IN; + p_ngc->gc_flags |= GC_BR_IN; + } + + p_gc++; + p_ngc++; + } + } + + p_gc = p_gs->gc; + for (i = 0; i < p_gs->tot_cells; i++) { - p_gc = p_gs->gc ; - for ( i = 0; i < p_gs->tot_cells ; i++ ) { - - /* reverse parity of edge clear (1==edge clear) */ - gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR ; - if ( gc_clear_flags & GC_L_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_L ; - } else - if ( gc_clear_flags & GC_B_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_B ; - } else - if ( gc_clear_flags & GC_R_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_R ; - } else - if ( gc_clear_flags & GC_T_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_T ; - } else { - /* all edges have something on them, do full test */ - p_gc->gc_flags |= GC_AIM_C ; - } - p_gc++ ; + /* reverse parity of edge clear (1==edge clear) */ + gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR; + if (gc_clear_flags & GC_L_EDGE_CLEAR) { + p_gc->gc_flags |= GC_AIM_L; } + else + if (gc_clear_flags & GC_B_EDGE_CLEAR) { + p_gc->gc_flags |= GC_AIM_B; + } + else + if (gc_clear_flags & GC_R_EDGE_CLEAR) { + p_gc->gc_flags |= GC_AIM_R; + } + else + if (gc_clear_flags & GC_T_EDGE_CLEAR) { + p_gc->gc_flags |= GC_AIM_T; + } + else { + /* all edges have something on them, do full test */ + p_gc->gc_flags |= GC_AIM_C; + } + p_gc++; + } } -int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps ) -{ -pGridRec p_gr ; -double slope, inv_slope ; - - if ( Near(ya, yb, eps) ) { - if ( Near(xa, xb, eps) ) { - /* edge is 0 length, so get rid of it */ - return( FALSE ) ; - } else { - /* horizontal line */ - slope = HUGE ; - inv_slope = 0.0 ; - } - } else { - if ( Near(xa, xb, eps) ) { - /* vertical line */ - slope = 0.0 ; - inv_slope = HUGE ; - } else { - slope = (xb-xa)/(yb-ya) ; - inv_slope = (yb-ya)/(xb-xa) ; - } - } +int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps) { + pGridRec p_gr; + double slope, inv_slope; - p_gc->tot_edges++ ; - if ( p_gc->tot_edges <= 1 ) { - p_gc->gr = (pGridRec)malloc( sizeof(GridRec) ) ; - } else { - p_gc->gr = (pGridRec)realloc( p_gc->gr, - p_gc->tot_edges * sizeof(GridRec) ) ; + if (Near(ya, yb, eps)) { + if (Near(xa, xb, eps)) { + /* edge is 0 length, so get rid of it */ + return(FALSE); } - MALLOC_CHECK( p_gc->gr ) ; - p_gr = &p_gc->gr[p_gc->tot_edges-1] ; - - p_gr->slope = slope ; - p_gr->inv_slope = inv_slope ; - - p_gr->xa = xa ; - p_gr->ya = ya ; - if ( xa <= xb ) { - p_gr->minx = xa ; - p_gr->maxx = xb ; - } else { - p_gr->minx = xb ; - p_gr->maxx = xa ; + else { + /* horizontal line */ + slope = HUGE; + inv_slope = 0.0; } - if ( ya <= yb ) { - p_gr->miny = ya ; - p_gr->maxy = yb ; - } else { - p_gr->miny = yb ; - p_gr->maxy = ya ; + } + else { + if (Near(xa, xb, eps)) { + /* vertical line */ + slope = 0.0; + inv_slope = HUGE; } - - /* P2 - P1 */ - p_gr->ax = xb - xa ; - p_gr->ay = yb - ya ; - - return( TRUE ) ; + else { + slope = (xb - xa) / (yb - ya); + inv_slope = (yb - ya) / (xb - xa); + } + } + + p_gc->tot_edges++; + if (p_gc->tot_edges <= 1) { + p_gc->gr = (pGridRec)malloc(sizeof(GridRec)); + } + else { + p_gc->gr = (pGridRec)realloc(p_gc->gr, + p_gc->tot_edges * sizeof(GridRec)); + } + MALLOC_CHECK(p_gc->gr); + p_gr = &p_gc->gr[p_gc->tot_edges - 1]; + + p_gr->slope = slope; + p_gr->inv_slope = inv_slope; + + p_gr->xa = xa; + p_gr->ya = ya; + if (xa <= xb) { + p_gr->minx = xa; + p_gr->maxx = xb; + } + else { + p_gr->minx = xb; + p_gr->maxx = xa; + } + if (ya <= yb) { + p_gr->miny = ya; + p_gr->maxy = yb; + } + else { + p_gr->miny = yb; + p_gr->maxy = ya; + } + + /* P2 - P1 */ + p_gr->ax = xb - xa; + p_gr->ay = yb - ya; + + return(TRUE); } /* Test point against grid and edges in the cell (if any). Algorithm: @@ -440,214 +457,219 @@ double slope, inv_slope ; * state of the edge or corner the ray went to and so determine the * state of the point (inside or outside). */ -int GridTest(pGridSet p_gs, pPipoint point) -{ -int j, count, init_flag ; -pGridCell p_gc ; -pGridRec p_gr ; -double tx, ty, xcell, ycell, bx,by,cx,cy, cornerx, cornery ; -double alpha, beta, denom ; -unsigned short gc_flags ; -int inside_flag ; - - /* first, is point inside bounding rectangle? */ - if ( ( ty = point->y ) < p_gs->miny || - ty >= p_gs->maxy || - ( tx = point->x ) < p_gs->minx || - tx >= p_gs->maxx ) { - - /* outside of box */ - inside_flag = FALSE ; - } else { - - /* what cell are we in? */ - ycell = ( ty - p_gs->miny ) * p_gs->inv_ydelta ; - xcell = ( tx - p_gs->minx ) * p_gs->inv_xdelta ; - p_gc = &p_gs->gc[((int)ycell)*p_gs->xres + (int)xcell] ; - - /* is cell simple? */ - count = p_gc->tot_edges ; - if ( count ) { - /* no, so find an edge which is free. */ - gc_flags = p_gc->gc_flags ; - p_gr = p_gc->gr ; - switch( gc_flags & GC_AIM ) { - case GC_AIM_L: - /* left edge is clear, shoot X- ray */ - /* note - this next statement requires that GC_BL_IN is 1 */ - inside_flag = gc_flags & GC_BL_IN ; - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if y is between edges */ - if ( ty >= p_gr->miny && ty < p_gr->maxy ) { - if ( tx > p_gr->maxx ) { - inside_flag = !inside_flag ; - } else if ( tx > p_gr->minx ) { - /* full computation */ - if ( ( p_gr->xa - - ( p_gr->ya - ty ) * p_gr->slope ) < tx ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_B: - /* bottom edge is clear, shoot Y+ ray */ - /* note - this next statement requires that GC_BL_IN is 1 */ - inside_flag = gc_flags & GC_BL_IN ; - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if x is between edges */ - if ( tx >= p_gr->minx && tx < p_gr->maxx ) { - if ( ty > p_gr->maxy ) { - inside_flag = !inside_flag ; - } else if ( ty > p_gr->miny ) { - /* full computation */ - if ( ( p_gr->ya - ( p_gr->xa - tx ) * - p_gr->inv_slope ) < ty ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_R: - /* right edge is clear, shoot X+ ray */ - inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ; - - /* TBD: Note, we could have sorted the edges to be tested - * by miny or somesuch, and so be able to cut testing - * short when the list's miny > point.y . - */ - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if y is between edges */ - if ( ty >= p_gr->miny && ty < p_gr->maxy ) { - if ( tx <= p_gr->minx ) { - inside_flag = !inside_flag ; - } else if ( tx <= p_gr->maxx ) { - /* full computation */ - if ( ( p_gr->xa - - ( p_gr->ya - ty ) * p_gr->slope ) >= tx ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_T: - /* top edge is clear, shoot Y+ ray */ - inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ; - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if x is between edges */ - if ( tx >= p_gr->minx && tx < p_gr->maxx ) { - if ( ty <= p_gr->miny ) { - inside_flag = !inside_flag ; - } else if ( ty <= p_gr->maxy ) { - /* full computation */ - if ( ( p_gr->ya - ( p_gr->xa - tx ) * - p_gr->inv_slope ) >= ty ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_C: - /* no edge is clear, bite the bullet and test - * against the bottom left corner. - * We use Franklin Antonio's algorithm (Graphics Gems III). - */ - /* TBD: Faster yet might be to test against the closest - * corner to the cell location, but our hope is that we - * rarely need to do this testing at all. - */ - inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN) ; - init_flag = TRUE ; - - /* get lower left corner coordinate */ - cornerx = p_gs->glx[(int)xcell] ; - cornery = p_gs->gly[(int)ycell] ; - for ( j = count+1 ; --j ; p_gr++ ) { - - /* quick out test: if test point is - * less than minx & miny, edge cannot overlap. - */ - if ( tx >= p_gr->minx && ty >= p_gr->miny ) { - - /* quick test failed, now check if test point and - * corner are on different sides of edge. - */ - if ( init_flag ) { - /* Compute these at most once for test */ - /* P3 - P4 */ - bx = tx - cornerx ; - by = ty - cornery ; - init_flag = FALSE ; - } - denom = p_gr->ay * bx - p_gr->ax * by ; - if ( denom != 0.0 ) { - /* lines are not collinear, so continue */ - /* P1 - P3 */ - cx = p_gr->xa - tx ; - cy = p_gr->ya - ty ; - alpha = by * cx - bx * cy ; - if ( denom > 0.0 ) { - if ( alpha < 0.0 || alpha >= denom ) { - /* test edge not hit */ - goto NextEdge ; - } - beta = p_gr->ax * cy - p_gr->ay * cx ; - if ( beta < 0.0 || beta >= denom ) { - /* polygon edge not hit */ - goto NextEdge ; - } - } else { - if ( alpha > 0.0 || alpha <= denom ) { - /* test edge not hit */ - goto NextEdge ; - } - beta = p_gr->ax * cy - p_gr->ay * cx ; - if ( beta > 0.0 || beta <= denom ) { - /* polygon edge not hit */ - goto NextEdge ; - } - } - inside_flag = !inside_flag ; - } - - } - NextEdge: ; - } - break ; - } - } else { - /* simple cell, so if lower left corner is in, - * then cell is inside. - */ - inside_flag = p_gc->gc_flags & GC_BL_IN ; - } +int GridTest(pGridSet p_gs, pPipoint point) { + int j, count, init_flag; + pGridCell p_gc; + pGridRec p_gr; + double tx, ty, xcell, ycell, bx, by, cx, cy, cornerx, cornery; + double alpha, beta, denom; + unsigned short gc_flags; + int inside_flag; + + /* first, is point inside bounding rectangle? */ + if ((ty = point->y) < p_gs->miny || + ty >= p_gs->maxy || + (tx = point->x) < p_gs->minx || + tx >= p_gs->maxx) { + + /* outside of box */ + inside_flag = FALSE; + } + else { + + /* what cell are we in? */ + ycell = (ty - p_gs->miny) * p_gs->inv_ydelta; + xcell = (tx - p_gs->minx) * p_gs->inv_xdelta; + p_gc = &p_gs->gc[((int)ycell)*p_gs->xres + (int)xcell]; + + /* is cell simple? */ + count = p_gc->tot_edges; + if (count) { + /* no, so find an edge which is free. */ + gc_flags = p_gc->gc_flags; + p_gr = p_gc->gr; + switch (gc_flags & GC_AIM) { + case GC_AIM_L: + /* left edge is clear, shoot X- ray */ + /* note - this next statement requires that GC_BL_IN is 1 */ + inside_flag = gc_flags & GC_BL_IN; + for (j = count + 1; --j; p_gr++) { + /* test if y is between edges */ + if (ty >= p_gr->miny && ty < p_gr->maxy) { + if (tx > p_gr->maxx) { + inside_flag = !inside_flag; + } + else if (tx > p_gr->minx) { + /* full computation */ + if ((p_gr->xa - + (p_gr->ya - ty) * p_gr->slope) < tx) { + inside_flag = !inside_flag; + } + } + } + } + break; + + case GC_AIM_B: + /* bottom edge is clear, shoot Y+ ray */ + /* note - this next statement requires that GC_BL_IN is 1 */ + inside_flag = gc_flags & GC_BL_IN; + for (j = count + 1; --j; p_gr++) { + /* test if x is between edges */ + if (tx >= p_gr->minx && tx < p_gr->maxx) { + if (ty > p_gr->maxy) { + inside_flag = !inside_flag; + } + else if (ty > p_gr->miny) { + /* full computation */ + if ((p_gr->ya - (p_gr->xa - tx) * + p_gr->inv_slope) < ty) { + inside_flag = !inside_flag; + } + } + } + } + break; + + case GC_AIM_R: + /* right edge is clear, shoot X+ ray */ + inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0; + + /* TBD: Note, we could have sorted the edges to be tested + * by miny or somesuch, and so be able to cut testing + * short when the list's miny > point.y . + */ + for (j = count + 1; --j; p_gr++) { + /* test if y is between edges */ + if (ty >= p_gr->miny && ty < p_gr->maxy) { + if (tx <= p_gr->minx) { + inside_flag = !inside_flag; + } + else if (tx <= p_gr->maxx) { + /* full computation */ + if ((p_gr->xa - + (p_gr->ya - ty) * p_gr->slope) >= tx) { + inside_flag = !inside_flag; + } + } + } + } + break; + + case GC_AIM_T: + /* top edge is clear, shoot Y+ ray */ + inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0; + for (j = count + 1; --j; p_gr++) { + /* test if x is between edges */ + if (tx >= p_gr->minx && tx < p_gr->maxx) { + if (ty <= p_gr->miny) { + inside_flag = !inside_flag; + } + else if (ty <= p_gr->maxy) { + /* full computation */ + if ((p_gr->ya - (p_gr->xa - tx) * + p_gr->inv_slope) >= ty) { + inside_flag = !inside_flag; + } + } + } + } + break; + + case GC_AIM_C: + /* no edge is clear, bite the bullet and test + * against the bottom left corner. + * We use Franklin Antonio's algorithm (Graphics Gems III). + */ + /* TBD: Faster yet might be to test against the closest + * corner to the cell location, but our hope is that we + * rarely need to do this testing at all. + */ + inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN); + init_flag = TRUE; + + /* get lower left corner coordinate */ + cornerx = p_gs->glx[(int)xcell]; + cornery = p_gs->gly[(int)ycell]; + for (j = count + 1; --j; p_gr++) { + + /* quick out test: if test point is + * less than minx & miny, edge cannot overlap. + */ + if (tx >= p_gr->minx && ty >= p_gr->miny) { + + /* quick test failed, now check if test point and + * corner are on different sides of edge. + */ + if (init_flag) { + /* Compute these at most once for test */ + /* P3 - P4 */ + bx = tx - cornerx; + by = ty - cornery; + init_flag = FALSE; + } + denom = p_gr->ay * bx - p_gr->ax * by; + if (denom != 0.0) { + /* lines are not collinear, so continue */ + /* P1 - P3 */ + cx = p_gr->xa - tx; + cy = p_gr->ya - ty; + alpha = by * cx - bx * cy; + if (denom > 0.0) { + if (alpha < 0.0 || alpha >= denom) { + /* test edge not hit */ + goto NextEdge; + } + beta = p_gr->ax * cy - p_gr->ay * cx; + if (beta < 0.0 || beta >= denom) { + /* polygon edge not hit */ + goto NextEdge; + } + } + else { + if (alpha > 0.0 || alpha <= denom) { + /* test edge not hit */ + goto NextEdge; + } + beta = p_gr->ax * cy - p_gr->ay * cx; + if (beta > 0.0 || beta <= denom) { + /* polygon edge not hit */ + goto NextEdge; + } + } + inside_flag = !inside_flag; + } + + } + NextEdge:; + } + break; + } + } + else { + /* simple cell, so if lower left corner is in, + * then cell is inside. + */ + inside_flag = p_gc->gc_flags & GC_BL_IN; } + } - return( inside_flag ) ; + return(inside_flag); } -void GridCleanup(pGridSet p_gs) -{ -int i ; -pGridCell p_gc ; +void GridCleanup(pGridSet p_gs) { + int i; + pGridCell p_gc; - for ( i = 0, p_gc = p_gs->gc - ; i < p_gs->tot_cells - ; i++, p_gc++ ) { + for (i = 0, p_gc = p_gs->gc + ; i < p_gs->tot_cells + ; i++, p_gc++) { - if ( p_gc->gr ) { - free( p_gc->gr ) ; - } + if (p_gc->gr) { + free(p_gc->gr); } - free( p_gs->glx ) ; - free( p_gs->gly ) ; - free( p_gs->gc ) ; + } + free(p_gs->glx); + free(p_gs->gly); + free(p_gs->gc); } diff --git a/thirdparty/ptinpoly.h b/thirdparty/ptinpoly.h index 7a0d26be..8ec7a0cf 100644 --- a/thirdparty/ptinpoly.h +++ b/thirdparty/ptinpoly.h @@ -7,8 +7,8 @@ extern "C" { #endif -/* =========================== System Related ============================= */ -/* SRAN initializes random number generator, if needed */ + /* =========================== System Related ============================= */ + /* SRAN initializes random number generator, if needed */ #define SRAN() srand(1) /* RAN01 returns a double from [0..1) */ #define RAN01() ((double)rand() / 32768.0) @@ -18,12 +18,12 @@ extern "C" #define GR_FULL_VERT 0x01 /* line crosses vertically */ #define GR_FULL_HORZ 0x02 /* line crosses horizontally */ -typedef struct { - double xa,ya ; - double minx, maxx, miny, maxy ; - double ax, ay ; - double slope, inv_slope ; -} GridRec, *pGridRec; + typedef struct { + double xa, ya; + double minx, maxx, miny, maxy; + double ax, ay; + double slope, inv_slope; + } GridRec, *pGridRec; #define GC_BL_IN 0x0001 /* bottom left corner is in (else out) */ #define GC_BR_IN 0x0002 /* bottom right corner is in (else out) */ @@ -52,42 +52,42 @@ typedef struct { GC_B_EDGE_HIT | \ GC_T_EDGE_HIT ) -typedef struct { - short tot_edges ; - unsigned short gc_flags ; - GridRec *gr ; -} GridCell, *pGridCell; + typedef struct { + short tot_edges; + unsigned short gc_flags; + GridRec *gr; + } GridCell, *pGridCell; -typedef struct { - int xres, yres ; /* grid size */ - int tot_cells ; /* xres * yres */ - double minx, maxx, miny, maxy ; /* bounding box */ - double xdelta, ydelta ; - double inv_xdelta, inv_ydelta ; - double *glx, *gly ; - GridCell *gc ; -} GridSet, *pGridSet ; + typedef struct { + int xres, yres; /* grid size */ + int tot_cells; /* xres * yres */ + double minx, maxx, miny, maxy; /* bounding box */ + double xdelta, ydelta; + double inv_xdelta, inv_ydelta; + double *glx, *gly; + GridCell *gc; + } GridSet, *pGridSet; -typedef struct { - double x, y; -} Pipoint, *pPipoint; + typedef struct { + double x, y; + } Pipoint, *pPipoint; -/* add a little to the limits of the polygon bounding box to avoid precision -* problems. -*/ + /* add a little to the limits of the polygon bounding box to avoid precision + * problems. + */ #define EPSILON 0.00001 -/* The following structure is associated with a polygon */ -typedef struct { - int id ; /* vertex number of edge */ - int full_cross ; /* 1 if extends from top to bottom */ - double minx, maxx ; /* X bounds for bin */ -} Edge, *pEdge ; + /* The following structure is associated with a polygon */ + typedef struct { + int id; /* vertex number of edge */ + int full_cross; /* 1 if extends from top to bottom */ + double minx, maxx; /* X bounds for bin */ + } Edge, *pEdge; -void GridSetup(pPipoint pgon[], int numverts, int resolution, pGridSet p_gs); -int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps); -void GridCleanup(pGridSet p_gs); -int GridTest(pGridSet p_gs, pPipoint point); + void GridSetup(pPipoint pgon[], int numverts, int resolution, pGridSet p_gs); + int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps); + void GridCleanup(pGridSet p_gs); + int GridTest(pGridSet p_gs, pPipoint point); #ifdef __cplusplus } From 4f9955280fdd5f51e00b4f276da74b242de82476 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 19 Nov 2018 16:42:09 +0100 Subject: [PATCH 03/53] correct citygml schema references --- src/Water.cpp | 8 ++++---- src/io.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Water.cpp b/src/Water.cpp index f5941610..41c82d12 100644 --- a/src/Water.cpp +++ b/src/Water.cpp @@ -108,10 +108,10 @@ void Water::get_citygml_imgeo(std::wostream& of) { std::string attribute; if (ondersteunend) { if (get_attribute("bgt-type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } else if (get_attribute("bgt_type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } if (get_attribute("plus-type", attribute)) { of << "" << attribute << ""; @@ -123,10 +123,10 @@ void Water::get_citygml_imgeo(std::wostream& of) { } else { if (get_attribute("bgt-type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } else if (get_attribute("bgt_type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } if (get_attribute("plus-type", attribute)) { of << "" << attribute << ""; diff --git a/src/io.cpp b/src/io.cpp index 87ea2c85..afa64958 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -67,7 +67,7 @@ void get_citygml_namespaces(std::wostream& of) { of << "xmlns:app=\"http://www.opengis.net/citygml/appearance/2.0\"\n"; of << "xmlns:tun=\"http://www.opengis.net/citygml/tunnel/2.0\"\n"; of << "xmlns:cif=\"http://www.opengis.net/citygml/cityfurniture/2.0\"\n"; - of << "xsi:schemaLocation=\"http://www.opengis.net/citygml/2.0 ./CityGML_2.0/CityGML.xsd\">\n"; + of << "xsi:schemaLocation=\"http://www.opengis.net/citygml/2.0 http://schemas.opengis.net/citygml/profiles/base/2.0/CityGML.xsd\">\n"; } void get_citygml_imgeo_namespaces(std::wostream& of) { From d2c56e280f99f89b22e738893296753819f287c2 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 20 Nov 2018 17:55:06 +0100 Subject: [PATCH 04/53] fix liblas version error catching --- src/Map3d.cpp | 90 ++++++++++++++++++++++++++------------------------- 1 file changed, 46 insertions(+), 44 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 4cfeb31c..c39e866b 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -1143,45 +1143,47 @@ void Map3d::extract_feature(OGRFeature *f, std::string layername, const char *id bool Map3d::add_las_file(PointFile pointFile) { std::clog << "Reading LAS/LAZ file: " << pointFile.filename << std::endl; std::ifstream ifs; - ifs.open(pointFile.filename.c_str(), std::ios::in | std::ios::binary); - if (ifs.is_open() == false) { - std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl; - return false; - } - //-- LAS classes to omit - std::vector liblasomits; - for (int i : pointFile.lasomits) { - liblasomits.push_back(liblas::Classification(i)); - } - - //-- read each point 1-by-1 - liblas::ReaderFactory f; - liblas::Reader reader = f.CreateWithStream(ifs); - liblas::Header const& header = reader.GetHeader(); - - //-- check if the file overlaps the polygons - liblas::Bounds bounds = header.GetExtent(); - liblas::Bounds polygonBounds = get_bounds(); - uint32_t pointCount = header.GetPointRecordsCount(); - if (polygonBounds.intersects(bounds)) { - std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n"; - if ((pointFile.thinning > 1)) { - std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus "; - std::clog << boost::locale::as::number << (pointCount / pointFile.thinning) << " are used)\n"; + + try { + ifs.open(pointFile.filename.c_str(), std::ios::in | std::ios::binary); + if (ifs.is_open() == false) { + std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl; + return false; } - else - std::clog << "\t(all points used, no skipping)\n"; - - if (pointFile.lasomits.empty() == false) { - std::clog << "\t(omitting LAS classes: "; - for (int i : pointFile.lasomits) - std::clog << i << " "; - std::clog << ")\n"; + //-- LAS classes to omit + std::vector liblasomits; + for (int i : pointFile.lasomits) { + liblasomits.push_back(liblas::Classification(i)); } - printProgressBar(0); - int i = 0; - - try { + + //-- read each point 1-by-1 + liblas::ReaderFactory f; + liblas::Reader reader = f.CreateWithStream(ifs); + liblas::Header const& header = reader.GetHeader(); + + //-- check if the file overlaps the polygons + liblas::Bounds bounds = header.GetExtent(); + liblas::Bounds polygonBounds = get_bounds(); + uint32_t pointCount = header.GetPointRecordsCount(); + + if (polygonBounds.intersects(bounds)) { + std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n"; + if ((pointFile.thinning > 1)) { + std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus "; + std::clog << boost::locale::as::number << (pointCount / pointFile.thinning) << " are used)\n"; + } + else + std::clog << "\t(all points used, no skipping)\n"; + + if (pointFile.lasomits.empty() == false) { + std::clog << "\t(omitting LAS classes: "; + for (int i : pointFile.lasomits) + std::clog << i << " "; + std::clog << ")\n"; + } + printProgressBar(0); + int i = 0; + while (reader.ReadNextPoint()) { liblas::Point const& p = reader.GetPoint(); //-- set the thinning filter @@ -1201,16 +1203,16 @@ bool Map3d::add_las_file(PointFile pointFile) { printProgressBar(100); std::clog << std::endl; } - catch (std::exception e) { - std::cerr << std::endl << e.what() << std::endl; - ifs.close(); - return false; + else { + std::clog << "\tskipping file, bounds do not intersect polygon extent\n"; } + ifs.close(); } - else { - std::clog << "\tskipping file, bounds do not intersect polygon extent\n"; + catch (std::exception e) { + std::cerr << std::endl << e.what() << std::endl; + ifs.close(); + return false; } - ifs.close(); return true; } From 5fcd2f7e7be7b8be0a8813ba962cadacff67d3ba Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 26 Nov 2018 12:33:35 +0100 Subject: [PATCH 05/53] set correct number output for xml streams --- src/Map3d.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index c39e866b..8dc232f0 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -454,6 +454,7 @@ bool Map3d::get_pdok_output(std::string filename) { //Add additional attribute describing CityGML of feature std::wstring_convert>converter; std::wstringstream ss; + ss << std::fixed << std::setprecision(3); f->get_citygml_imgeo(ss); std::string gmlAttribute = converter.to_bytes(ss.str()); ss.clear(); @@ -538,6 +539,7 @@ bool Map3d::get_pdok_citygml_output(std::string filename) { //Add additional attribute describing CityGML of feature std::wstring_convert>converter; std::wstringstream ss; + ss << std::fixed << std::setprecision(3); f->get_citygml(ss); std::string gmlAttribute = converter.to_bytes(ss.str()); ss.clear(); From c3284966eb14199a0d3f5ee80949fd6a0c8859db Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 28 Nov 2018 09:31:42 +0100 Subject: [PATCH 06/53] set correct z-coordinate precision --- src/io.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/io.cpp b/src/io.cpp index afa64958..00418927 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -132,11 +132,11 @@ void get_extruded_line_gml(std::wostream& of, Point2* a, Point2* b, double high, of << ""; of << ""; of << ""; - of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << low << " "; - of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << low << " "; - of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << high << " "; - of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << high << " "; - of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << low; + of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3) << " "; + of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << low << std::setprecision(3) << " "; + of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << high << std::setprecision(3) << " "; + of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << high << std::setprecision(3) << " "; + of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3); of << ""; of << ""; of << ""; From 9f123ef9dc6c9231a31365a305a55ca145fdb339 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 5 Dec 2018 15:52:13 +0100 Subject: [PATCH 07/53] change to lastools --- src/Map3d.cpp | 155 ++++++++++++++++++++++------------------ src/Map3d.h | 8 ++- src/Water.cpp | 8 +-- src/definitions.h | 2 +- src/io.cpp | 4 +- src/main.cpp | 4 +- vs_build/3dfier.vcxproj | 12 ++-- 7 files changed, 105 insertions(+), 88 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index f6740dda..ae8cc08a 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -53,6 +53,10 @@ Map3d::Map3d() { _threshold_jump_edges = 50; _requestedExtent = Box2(Point2(0, 0), Point2(0, 0)); _bbox = Box2(Point2(999999, 999999), Point2(-999999, -999999)); + _minxradius = 999999; + _maxxradius = 999999; + _minyradius = -999999; + _maxyradius = -999999; } Map3d::~Map3d() { @@ -155,10 +159,12 @@ Box2 Map3d::get_bbox() { return _bbox; } -liblas::Bounds Map3d::get_bounds() { - double radius = std::max(_radius_vertex_elevation, _building_radius_vertex_elevation); - liblas::Bounds bounds(_bbox.min_corner().x() - radius, _bbox.min_corner().y() - radius, _bbox.max_corner().x() + radius, _bbox.max_corner().y() + radius); - return bounds; +bool Map3d::check_bounds(const double xmin, const double xmax, const double ymin, const double ymax) { + if ((xmin < _maxxradius || xmax > _minxradius) && + (ymin < _maxyradius || ymax > _minyradius)) { + return true; + } + return false; } bool Map3d::get_cityjson(std::string filename) { @@ -189,8 +195,7 @@ bool Map3d::get_cityjson(std::string filename) { boost::split(c, p, boost::is_any_of(" ")); j["vertices"].push_back({std::stod(c[0], NULL), std::stod(c[1], NULL), std::stod(c[2], NULL) }); } - std::ofstream o(filename); - // o << j.dump(2) << std::endl; + std::ofstream o(filename); o << j.dump() << std::endl; return true; } @@ -736,27 +741,30 @@ const std::vector& Map3d::get_polygons3d() { return _lsFeatures; } -void Map3d::add_elevation_point(liblas::Point const& laspt) { +void Map3d::add_elevation_point(LASpoint const& laspt) { + //-- only process last returns; + //-- although perhaps not smart for vegetation/forest in the future + //-- TODO: always ignore the non-last-return points? + if (laspt.return_number != laspt.number_of_returns) + return; + std::vector re; - Point2 minp(laspt.GetX() - _radius_vertex_elevation, laspt.GetY() - _radius_vertex_elevation); - Point2 maxp(laspt.GetX() + _radius_vertex_elevation, laspt.GetY() + _radius_vertex_elevation); + float x = laspt.get_x(); + float y = laspt.get_y(); + Point2 minp(x - _radius_vertex_elevation, y - _radius_vertex_elevation); + Point2 maxp(x + _radius_vertex_elevation, y + _radius_vertex_elevation); Box2 querybox(minp, maxp); _rtree.query(bgi::intersects(querybox), std::back_inserter(re)); - minp = Point2(laspt.GetX() - _building_radius_vertex_elevation, laspt.GetY() - _building_radius_vertex_elevation); - maxp = Point2(laspt.GetX() + _building_radius_vertex_elevation, laspt.GetY() + _building_radius_vertex_elevation); + minp = Point2(x - _building_radius_vertex_elevation, y - _building_radius_vertex_elevation); + maxp = Point2(x + _building_radius_vertex_elevation, y + _building_radius_vertex_elevation); querybox = Box2(minp, maxp); _rtree_buildings.query(bgi::intersects(querybox), std::back_inserter(re)); for (auto& v : re) { TopoFeature* f = v.second; float radius = _radius_vertex_elevation; - //-- only process last returns; - //-- although perhaps not smart for vegetation/forest in the future - //-- TODO: always ignore the non-last-return points? - if (laspt.GetReturnNumber() != laspt.GetNumberOfReturns()) - continue; - int c = laspt.GetClassification().GetClass(); + int c = (int)laspt.classification; bool bInsert = false; bool bWithin = false; if (f->get_class() == BUILDING) { @@ -817,10 +825,10 @@ void Map3d::add_elevation_point(liblas::Point const& laspt) { bWithin = true; } } - + if (bInsert == true) { //-- only insert if in the allowed LAS classes - Point2 p(laspt.GetX(), laspt.GetY()); - f->add_elevation_point(p, laspt.GetZ(), radius, c, bWithin); + Point2 p(x, y); + f->add_elevation_point(p, laspt.get_z(), radius, c, bWithin); } } } @@ -935,6 +943,11 @@ bool Map3d::construct_rtree() { std::min(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds()))), Point2(std::max(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())), std::max(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())))); + double radius = std::max(_radius_vertex_elevation, _building_radius_vertex_elevation); + _minxradius = std::min(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())) - radius; + _maxxradius = std::min(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())) + radius; + _minyradius = std::max(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())) - radius; + _maxyradius = std::max(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())) + radius; return true; } @@ -1135,57 +1148,57 @@ void Map3d::extract_feature(OGRFeature *f, std::string layername, const char *id } } -//-- http://www.liblas.org/tutorial/cpp.html#applying-filters-to-a-reader-to-extract-specified-classes bool Map3d::add_las_file(PointFile pointFile) { std::clog << "Reading LAS/LAZ file: " << pointFile.filename << std::endl; - std::ifstream ifs; - ifs.open(pointFile.filename.c_str(), std::ios::in | std::ios::binary); - if (ifs.is_open() == false) { - std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl; - return false; - } - //-- LAS classes to omit - std::vector liblasomits; - for (int i : pointFile.lasomits) { - liblasomits.push_back(liblas::Classification(i)); - } - - //-- read each point 1-by-1 - liblas::ReaderFactory f; - liblas::Reader reader = f.CreateWithStream(ifs); - liblas::Header const& header = reader.GetHeader(); - - //-- check if the file overlaps the polygons - liblas::Bounds bounds = header.GetExtent(); - liblas::Bounds polygonBounds = get_bounds(); - uint32_t pointCount = header.GetPointRecordsCount(); - if (polygonBounds.intersects(bounds)) { - std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n"; - if ((pointFile.thinning > 1)) { - std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus "; - std::clog << boost::locale::as::number << (pointCount / pointFile.thinning) << " are used)\n"; - } - else - std::clog << "\t(all points used, no skipping)\n"; - - if (pointFile.lasomits.empty() == false) { - std::clog << "\t(omitting LAS classes: "; - for (int i : pointFile.lasomits) - std::clog << i << " "; - std::clog << ")\n"; + + LASreadOpener lasreadopener; + lasreadopener.set_file_name(pointFile.filename.c_str()); + //-- set to compute bounding box + lasreadopener.set_populate_header(true); + LASreader* lasreader = lasreadopener.open(); + + try { + //-- check if file is open + if (lasreader == 0) { + std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl; + return false; } - printProgressBar(0); - int i = 0; - - try { - while (reader.ReadNextPoint()) { - liblas::Point const& p = reader.GetPoint(); + LASheader header = lasreader->header; + + if (check_bounds(header.min_x, header.max_x, header.min_y, header.max_y)) { + //-- LAS classes to omit + std::vector lasomits; + for (int i : pointFile.lasomits) { + lasomits.push_back(i); + } + + //-- read each point 1-by-1 + uint32_t pointCount = header.number_of_point_records; + + std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n"; + if ((pointFile.thinning > 1)) { + std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus "; + std::clog << boost::locale::as::number << (pointCount / pointFile.thinning) << " are used)\n"; + } + else + std::clog << "\t(all points used, no skipping)\n"; + + if (pointFile.lasomits.empty() == false) { + std::clog << "\t(omitting LAS classes: "; + for (int i : pointFile.lasomits) + std::clog << i << " "; + std::clog << ")\n"; + } + printProgressBar(0); + int i = 0; + while (lasreader->read_point()) { + LASpoint const& p = lasreader->point; //-- set the thinning filter if (i % pointFile.thinning == 0) { //-- set the classification filter - if (std::find(liblasomits.begin(), liblasomits.end(), p.GetClassification()) == liblasomits.end()) { + if (std::find(lasomits.begin(), lasomits.end(), (int)p.classification) == lasomits.end()) { //-- set the bounds filter - if (polygonBounds.contains(p)) { + if (check_bounds(p.X, p.X, p.Y, p.Y)) { this->add_elevation_point(p); } } @@ -1197,16 +1210,16 @@ bool Map3d::add_las_file(PointFile pointFile) { printProgressBar(100); std::clog << std::endl; } - catch (std::exception e) { - std::cerr << std::endl << e.what() << std::endl; - ifs.close(); - return false; + else { + std::clog << "\tskipping file, bounds do not intersect polygon extent\n"; } + lasreader->close(); } - else { - std::clog << "\tskipping file, bounds do not intersect polygon extent\n"; + catch (std::exception e) { + std::cerr << std::endl << e.what() << std::endl; + lasreader->close(); + return false; } - ifs.close(); return true; } diff --git a/src/Map3d.h b/src/Map3d.h index 931cb017..ab969720 100644 --- a/src/Map3d.h +++ b/src/Map3d.h @@ -56,13 +56,13 @@ class Map3d { bool construct_rtree(); bool threeDfy(bool stitching = true); bool construct_CDT(); - void add_elevation_point(liblas::Point const& laspt); + void add_elevation_point(LASpoint const& laspt); void cleanup_elevations(); unsigned long get_num_polygons(); const std::vector& get_polygons3d(); Box2 get_bbox(); - liblas::Bounds get_bounds(); + bool check_bounds(const double xmin, const double xmax, const double ymin, const double ymax); void get_citygml(std::wostream& of); void get_citygml_multifile(std::string); @@ -134,6 +134,10 @@ class Map3d { float _building_radius_vertex_elevation; int _threshold_jump_edges; //-- in cm/integer Box2 _bbox; + double _minxradius; + double _maxxradius; + double _minyradius; + double _maxyradius; Box2 _requestedExtent; //-- storing the LAS allowed for each TopoFeature diff --git a/src/Water.cpp b/src/Water.cpp index f5941610..4668ad68 100644 --- a/src/Water.cpp +++ b/src/Water.cpp @@ -108,10 +108,10 @@ void Water::get_citygml_imgeo(std::wostream& of) { std::string attribute; if (ondersteunend) { if (get_attribute("bgt-type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } else if (get_attribute("bgt_type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } if (get_attribute("plus-type", attribute)) { of << "" << attribute << ""; @@ -123,10 +123,10 @@ void Water::get_citygml_imgeo(std::wostream& of) { } else { if (get_attribute("bgt-type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } else if (get_attribute("bgt_type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } if (get_attribute("plus-type", attribute)) { of << "" << attribute << ""; diff --git a/src/definitions.h b/src/definitions.h index bd0cf00d..8ea41ac1 100644 --- a/src/definitions.h +++ b/src/definitions.h @@ -6,7 +6,7 @@ #include #include -#include +#include #include #include diff --git a/src/io.cpp b/src/io.cpp index 28bccd41..ae608d79 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -67,7 +67,7 @@ void get_citygml_namespaces(std::wostream& of) { of << "xmlns:app=\"http://www.opengis.net/citygml/appearance/2.0\"\n"; of << "xmlns:tun=\"http://www.opengis.net/citygml/tunnel/2.0\"\n"; of << "xmlns:cif=\"http://www.opengis.net/citygml/cityfurniture/2.0\"\n"; - of << "xsi:schemaLocation=\"http://www.opengis.net/citygml/2.0 ./CityGML_2.0/CityGML.xsd\">\n"; + of << "xsi:schemaLocation=\"http://www.opengis.net/citygml/2.0 http://schemas.opengis.net/citygml/profiles/base/2.0/CityGML.xsd\">\n"; } void get_citygml_imgeo_namespaces(std::wostream& of) { @@ -136,7 +136,7 @@ void get_extruded_line_gml(std::wostream& of, Point2* a, Point2* b, double high, of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << low << std::setprecision(3) << " "; of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << high << std::setprecision(3) << " "; of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << high << std::setprecision(3) << " "; - of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3); + of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3); of << ""; of << ""; of << ""; diff --git a/src/main.cpp b/src/main.cpp index eb518d73..3134bc2a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -36,7 +36,7 @@ #include #include -std::string VERSION = "1.1"; +std::string VERSION = "1.1.1"; bool validate_yaml(const char* arg, std::set& allowedFeatures); int main(int argc, const char * argv[]); @@ -129,7 +129,7 @@ int main(int argc, const char * argv[]) { } if (vm.count("version")) { std::cout << "3dfier " << VERSION << std::endl; - std::cout << liblas::GetFullVersion() << std::endl; + std::cout << "LASlib " << LAS_TOOLS_VERSION << std::endl; std::cout << "GDAL " << GDALVersionInfo("--version") << std::endl; //-- TODO : put here the date and/or the git-commit? return EXIT_SUCCESS; diff --git a/vs_build/3dfier.vcxproj b/vs_build/3dfier.vcxproj index e341354a..fcdebd6a 100644 --- a/vs_build/3dfier.vcxproj +++ b/vs_build/3dfier.vcxproj @@ -58,7 +58,7 @@ - $(LIBLAS_ROOT)\include;$(LASZIP_ROOT)\include;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories) + $(LASLIB_ROOT)\inc;$(LASZIP_ROOT)\src;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories) $(IntDir) Default Sync @@ -89,8 +89,8 @@ /machine:x64 %(AdditionalOptions) - gdal_i.lib;liblas.lib;laszip.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib - $(LIBLAS_ROOT)\bin;$(LIBLAS_ROOT)\lib;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) + gdal_i.lib;LASlib.lib;CGAL_Core-vc140-mt-4.12.lib;CGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib + $(LASLIB_ROOT)\lib\Release;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) @@ -106,7 +106,7 @@ - $(LIBLAS_ROOT)\include;$(LASZIP_ROOT)\include;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories) + $(LASLIB_ROOT)\inc;$(LASZIP_ROOT)\src;$(BOOST_ROOT);$(OSGEO4W_ROOT)\include;$(YAML-CPP_DIR)\include;$(CGAL_DIR)\include;$(CGAL_DIR)\auxiliary\gmp\include;..\thirdparty;%(AdditionalIncludeDirectories) $(IntDir) Default Sync @@ -137,8 +137,8 @@ %(Filename)_p.c - gdal_i.lib;liblas.lib;laszip.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib - $(LIBLAS_ROOT)\bin;$(LIBLAS_ROOT)\lib;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) + gdal_i.lib;LASlib.lib;CGAL_Core-vc140-mt-4.12.lib;CGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib + $(LASLIB_ROOT)\lib\Release;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) From 0857b52d83eadc91ed7237afea2378d54a1943d8 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 5 Dec 2018 15:56:09 +0100 Subject: [PATCH 08/53] remove tab --- src/Map3d.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index ae8cc08a..dc4cc8ac 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -195,8 +195,8 @@ bool Map3d::get_cityjson(std::string filename) { boost::split(c, p, boost::is_any_of(" ")); j["vertices"].push_back({std::stod(c[0], NULL), std::stod(c[1], NULL), std::stod(c[2], NULL) }); } - std::ofstream o(filename); - o << j.dump() << std::endl; + std::ofstream o(filename); + o << j.dump() << std::endl; return true; } From 86a915816ff9468f24b447b0e797b9e125182085 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 5 Dec 2018 16:11:56 +0100 Subject: [PATCH 09/53] styling --- src/Map3d.cpp | 2 +- src/Map3d.h | 3 +-- vs_build/3dfier.vcxproj | 6 +++--- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index dc4cc8ac..5b7bd9e7 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -825,7 +825,6 @@ void Map3d::add_elevation_point(LASpoint const& laspt) { bWithin = true; } } - if (bInsert == true) { //-- only insert if in the allowed LAS classes Point2 p(x, y); f->add_elevation_point(p, laspt.get_z(), radius, c, bWithin); @@ -943,6 +942,7 @@ bool Map3d::construct_rtree() { std::min(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds()))), Point2(std::max(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())), std::max(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())))); + double radius = std::max(_radius_vertex_elevation, _building_radius_vertex_elevation); _minxradius = std::min(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())) - radius; _maxxradius = std::min(bg::get(_rtree.bounds()), bg::get(_rtree_buildings.bounds())) + radius; diff --git a/src/Map3d.h b/src/Map3d.h index ab969720..006cd2b1 100644 --- a/src/Map3d.h +++ b/src/Map3d.h @@ -56,7 +56,7 @@ class Map3d { bool construct_rtree(); bool threeDfy(bool stitching = true); bool construct_CDT(); - void add_elevation_point(LASpoint const& laspt); + void add_elevation_point(LASpoint const& laspt); void cleanup_elevations(); unsigned long get_num_polygons(); @@ -148,7 +148,6 @@ class Map3d { NodeColumn _nc_building_walls; std::unordered_map _bridge_stitches; std::vector _lsFeatures; - std::vector _allowed_layers; bgi::rtree< PairIndexed, bgi::rstar<16> > _rtree; bgi::rtree< PairIndexed, bgi::rstar<16> > _rtree_buildings; diff --git a/vs_build/3dfier.vcxproj b/vs_build/3dfier.vcxproj index fcdebd6a..c325f51f 100644 --- a/vs_build/3dfier.vcxproj +++ b/vs_build/3dfier.vcxproj @@ -1,4 +1,4 @@ - + @@ -90,7 +90,7 @@ /machine:x64 %(AdditionalOptions) gdal_i.lib;LASlib.lib;CGAL_Core-vc140-mt-4.12.lib;CGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib - $(LASLIB_ROOT)\lib\Release;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) + $(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) @@ -138,7 +138,7 @@ gdal_i.lib;LASlib.lib;CGAL_Core-vc140-mt-4.12.lib;CGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib - $(LASLIB_ROOT)\lib\Release;$(LASZIP_ROOT)\lib;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) + $(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) From 3eb0bd5a68bf523b9d46e28086fcbea002cceba9 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 5 Dec 2018 16:27:33 +0100 Subject: [PATCH 10/53] fix fstream import --- src/definitions.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/definitions.h b/src/definitions.h index 8ea41ac1..866e1cf7 100644 --- a/src/definitions.h +++ b/src/definitions.h @@ -2,6 +2,7 @@ #define __3DFIER__Definitions__ #include +#include #include #include #include From 139b419f3fd836f2a995ddd963d5711f7e5efeaf Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 5 Dec 2018 16:30:17 +0100 Subject: [PATCH 11/53] reset cgal to 4.8 --- vs_build/3dfier.vcxproj | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vs_build/3dfier.vcxproj b/vs_build/3dfier.vcxproj index c325f51f..4dec2b74 100644 --- a/vs_build/3dfier.vcxproj +++ b/vs_build/3dfier.vcxproj @@ -89,7 +89,7 @@ /machine:x64 %(AdditionalOptions) - gdal_i.lib;LASlib.lib;CGAL_Core-vc140-mt-4.12.lib;CGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib + gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib $(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) @@ -137,7 +137,7 @@ %(Filename)_p.c - gdal_i.lib;LASlib.lib;CGAL_Core-vc140-mt-4.12.lib;CGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib + gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib $(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) From 3cb3356be59adb0b2ca2575ed2cfd4adff139155 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 27 Dec 2018 12:29:44 +0100 Subject: [PATCH 12/53] remove ptpoly files --- src/TopoFeature.h | 1 - thirdparty/ptinpoly.c | 651 ------------------------------------------ thirdparty/ptinpoly.h | 94 ------ 3 files changed, 746 deletions(-) delete mode 100644 thirdparty/ptinpoly.c delete mode 100644 thirdparty/ptinpoly.h diff --git a/src/TopoFeature.h b/src/TopoFeature.h index 756eed7a..83ebb936 100644 --- a/src/TopoFeature.h +++ b/src/TopoFeature.h @@ -34,7 +34,6 @@ #include "io.h" #include "polyfit.hpp" #include "nlohmann-json/json.hpp" -#include "ptinpoly.h" class TopoFeature { public: diff --git a/thirdparty/ptinpoly.c b/thirdparty/ptinpoly.c deleted file mode 100644 index 7c833924..00000000 --- a/thirdparty/ptinpoly.c +++ /dev/null @@ -1,651 +0,0 @@ -/* ptinpoly.c - point in polygon inside/outside code. - - by Eric Haines, 3D/Eye Inc, erich@eye.com - - This code contains the following algorithms: - grid testing - grid imposed on polygon -*/ - -#include -#include -#include -#include "ptinpoly.h" - -#define X 0 -#define Y 1 - -#ifndef TRUE -#define TRUE 1 -#define FALSE 0 -#endif - -#ifndef HUGE -#define HUGE 1.797693134862315e+308 -#endif - -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif - -/* test if a & b are within epsilon. Favors cases where a < b */ -#define Near(a,b,eps) ( ((b)-(eps)<(a)) && ((a)-(eps)<(b)) ) - -#define MALLOC_CHECK( a ) if ( !(a) ) { \ - fprintf( stderr, "out of memory\n" ) ; \ - exit(1) ; \ - } - -/* ======= Grid algorithm ================================================= */ - -/* Impose a grid upon the polygon and test only the local edges against the - * point. - * - * Call setup with 2D polygon _pgon_ with _numverts_ number of vertices, - * grid resolution _resolution_ and a pointer to a grid structure _p_gs_. - * Call testing procedure with a pointer to this array and test point _point_, - * returns 1 if inside, 0 if outside. - * Call cleanup with pointer to grid structure to free space. - */ - -/* Strategy for setup: - * Get bounds of polygon, allocate grid. - * "Walk" each edge of the polygon and note which edges have been crossed - * and what cells are entered (points on a grid edge are always considered - * to be above that edge). Keep a record of the edges overlapping a cell. - * For cells with edges, determine if any cell border has no edges passing - * through it and so can be used for shooting a test ray. - * Keep track of the parity of the x (horizontal) grid cell borders for - * use in determining whether the grid corners are inside or outside. - */ -void GridSetup(pPipoint pgon[], int numverts, int resolution, pGridSet p_gs) -{ -double *vtx0, *vtx1, *vtxa, *vtxb, *p_gl ; -int i, j, gc_clear_flags ; -double vx0, vx1, vy0, vy1, gxdiff, gydiff, eps ; -pGridCell p_gc, p_ngc ; -double xdiff, ydiff, tmax, inv_x, inv_y, xdir, ydir, t_near, tx, ty ; -double tgcx, tgcy ; -int gcx, gcy, sign_x ; -int y_flag, io_state ; - - p_gs->xres = p_gs->yres = resolution ; - p_gs->tot_cells = p_gs->xres * p_gs->yres ; - p_gs->glx = (double *)malloc( (p_gs->xres+1) * sizeof(double)); - MALLOC_CHECK( p_gs->glx ) ; - p_gs->gly = (double *)malloc( (p_gs->yres+1) * sizeof(double)); - MALLOC_CHECK( p_gs->gly ) ; - p_gs->gc = (pGridCell)malloc( p_gs->tot_cells * sizeof(GridCell)); - MALLOC_CHECK( p_gs->gc ) ; - - p_gs->minx = - p_gs->maxx = pgon[0]->x ; - p_gs->miny = - p_gs->maxy = pgon[0]->y ; - - /* find bounds of polygon */ - for ( i = 1 ; i < numverts ; i++ ) { - vx0 = pgon[i]->x ; - if ( p_gs->minx > vx0 ) { - p_gs->minx = vx0 ; - } else if ( p_gs->maxx < vx0 ) { - p_gs->maxx = vx0 ; - } - - vy0 = pgon[i]->y ; - if ( p_gs->miny > vy0 ) { - p_gs->miny = vy0 ; - } else if ( p_gs->maxy < vy0 ) { - p_gs->maxy = vy0 ; - } - } - - /* add a little to the bounds to ensure everything falls inside area */ - gxdiff = p_gs->maxx - p_gs->minx ; - gydiff = p_gs->maxy - p_gs->miny ; - p_gs->minx -= EPSILON * gxdiff ; - p_gs->maxx += EPSILON * gxdiff ; - p_gs->miny -= EPSILON * gydiff ; - p_gs->maxy += EPSILON * gydiff ; - - /* avoid roundoff problems near corners by not getting too close to them */ - eps = 1e-9 * ( gxdiff + gydiff ) ; - - /* use the new bounds to compute cell widths */ - TryAgain: - p_gs->xdelta = - (p_gs->maxx-p_gs->minx) / (double)p_gs->xres ; - p_gs->inv_xdelta = 1.0 / p_gs->xdelta ; - - p_gs->ydelta = - (p_gs->maxy-p_gs->miny) / (double)p_gs->yres ; - p_gs->inv_ydelta = 1.0 / p_gs->ydelta ; - - for ( i = 0, p_gl = p_gs->glx ; i < p_gs->xres ; i++ ) { - *p_gl++ = p_gs->minx + i * p_gs->xdelta ; - } - /* make last grid corner precisely correct */ - *p_gl = p_gs->maxx ; - - for ( i = 0, p_gl = p_gs->gly ; i < p_gs->yres ; i++ ) { - *p_gl++ = p_gs->miny + i * p_gs->ydelta ; - } - *p_gl = p_gs->maxy ; - - for ( i = 0, p_gc = p_gs->gc ; i < p_gs->tot_cells ; i++, p_gc++ ) { - p_gc->tot_edges = 0 ; - p_gc->gc_flags = 0x0 ; - p_gc->gr = NULL ; - } - - /* loop through edges and insert into grid structure */ - vtx0 = pgon[numverts-1] ; - for ( i = 0 ; i < numverts ; i++ ) { - vtx1 = pgon[i] ; - - if ( vtx0[Y] < vtx1[Y] ) { - vtxa = vtx0 ; - vtxb = vtx1 ; - } else { - vtxa = vtx1 ; - vtxb = vtx0 ; - } - - /* Set x variable for the direction of the ray */ - xdiff = vtxb[X] - vtxa[X] ; - ydiff = vtxb[Y] - vtxa[Y] ; - tmax = sqrt( xdiff * xdiff + ydiff * ydiff ) ; - - /* if edge is of 0 length, ignore it (useless edge) */ - if ( tmax == 0.0 ) goto NextEdge ; - - xdir = xdiff / tmax ; - ydir = ydiff / tmax ; - - gcx = (int)(( vtxa[X] - p_gs->minx ) * p_gs->inv_xdelta) ; - gcy = (int)(( vtxa[Y] - p_gs->miny ) * p_gs->inv_ydelta) ; - - /* get information about slopes of edge, etc */ - if ( vtxa[X] == vtxb[X] ) { - sign_x = 0 ; - tx = HUGE ; - } else { - inv_x = tmax / xdiff ; - tx = p_gs->xdelta * (double)gcx + p_gs->minx - vtxa[X] ; - if ( vtxa[X] < vtxb[X] ) { - sign_x = 1 ; - tx += p_gs->xdelta ; - tgcx = p_gs->xdelta * inv_x ; - } else { - sign_x = -1 ; - tgcx = -p_gs->xdelta * inv_x ; - } - tx *= inv_x ; - } - - if ( vtxa[Y] == vtxb[Y] ) { - ty = HUGE ; - } else { - inv_y = tmax / ydiff ; - ty = (p_gs->ydelta * (double)(gcy+1) + p_gs->miny - vtxa[Y]) - * inv_y ; - tgcy = p_gs->ydelta * inv_y ; - } - - p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ; - - vx0 = vtxa[X] ; - vy0 = vtxa[Y] ; - - t_near = 0.0 ; - - do { - /* choose the next boundary, but don't move yet */ - if ( tx <= ty ) { - gcx += sign_x ; - - ty -= tx ; - t_near += tx ; - tx = tgcx ; - - /* note which edge is hit when leaving this cell */ - if ( t_near < tmax ) { - if ( sign_x > 0 ) { - p_gc->gc_flags |= GC_R_EDGE_HIT ; - vx1 = p_gs->glx[gcx] ; - } else { - p_gc->gc_flags |= GC_L_EDGE_HIT ; - vx1 = p_gs->glx[gcx+1] ; - } - - /* get new location */ - vy1 = t_near * ydir + vtxa[Y] ; - } else { - /* end of edge, so get exact value */ - vx1 = vtxb[X] ; - vy1 = vtxb[Y] ; - } - - y_flag = FALSE ; - - } else { - - gcy++ ; - - tx -= ty ; - t_near += ty ; - ty = tgcy ; - - /* note top edge is hit when leaving this cell */ - if ( t_near < tmax ) { - p_gc->gc_flags |= GC_T_EDGE_HIT ; - /* this toggles the parity bit */ - p_gc->gc_flags ^= GC_T_EDGE_PARITY ; - - /* get new location */ - vx1 = t_near * xdir + vtxa[X] ; - vy1 = p_gs->gly[gcy] ; - } else { - /* end of edge, so get exact value */ - vx1 = vtxb[X] ; - vy1 = vtxb[Y] ; - } - - y_flag = TRUE ; - } - - /* check for corner crossing, then mark the cell we're in */ - if ( !AddGridRecAlloc( p_gc, vx0, vy0, vx1, vy1, eps ) ) { - /* warning, danger - we have just crossed a corner. - * There are all kinds of topological messiness we could - * do to get around this case, but they're a headache. - * The simplest recovery is just to change the extents a bit - * and redo the meshing, so that hopefully no edges will - * perfectly cross a corner. Since it's a preprocess, we - * don't care too much about the time to do it. - */ - - /* clean out all grid records */ - for ( i = 0, p_gc = p_gs->gc - ; i < p_gs->tot_cells - ; i++, p_gc++ ) { - - if ( p_gc->gr ) { - free( p_gc->gr ) ; - } - } - - /* make the bounding box ever so slightly larger, hopefully - * changing the alignment of the corners. - */ - p_gs->minx -= EPSILON * gxdiff * 0.24 ; - p_gs->miny -= EPSILON * gydiff * 0.10 ; - - /* yes, it's the dreaded goto - run in fear for your lives! */ - goto TryAgain ; - } - - if ( t_near < tmax ) { - /* note how we're entering the next cell */ - /* TBD: could be done faster by incrementing index in the - * incrementing code, above */ - p_gc = &p_gs->gc[gcy*p_gs->xres+gcx] ; - - if ( y_flag ) { - p_gc->gc_flags |= GC_B_EDGE_HIT ; - /* this toggles the parity bit */ - p_gc->gc_flags ^= GC_B_EDGE_PARITY ; - } else { - p_gc->gc_flags |= - ( sign_x > 0 ) ? GC_L_EDGE_HIT : GC_R_EDGE_HIT ; - } - } - - vx0 = vx1 ; - vy0 = vy1 ; - } - /* have we gone further than the end of the edge? */ - while ( t_near < tmax ) ; - - NextEdge: - vtx0 = vtx1 ; - } - - /* the grid is all set up, now set up the inside/outside value of each - * corner. - */ - p_gc = p_gs->gc ; - p_ngc = &p_gs->gc[p_gs->xres] ; - - /* we know the bottom and top rows are all outside, so no flag is set */ - for ( i = 1; i < p_gs->yres ; i++ ) { - /* start outside */ - io_state = 0x0 ; - - for ( j = 0; j < p_gs->xres ; j++ ) { - - if ( io_state ) { - /* change cell left corners to inside */ - p_gc->gc_flags |= GC_TL_IN ; - p_ngc->gc_flags |= GC_BL_IN ; - } - - if ( p_gc->gc_flags & GC_T_EDGE_PARITY ) { - io_state = !io_state ; - } - - if ( io_state ) { - /* change cell right corners to inside */ - p_gc->gc_flags |= GC_TR_IN ; - p_ngc->gc_flags |= GC_BR_IN ; - } - - p_gc++ ; - p_ngc++ ; - } - } - - p_gc = p_gs->gc ; - for ( i = 0; i < p_gs->tot_cells ; i++ ) { - - /* reverse parity of edge clear (1==edge clear) */ - gc_clear_flags = p_gc->gc_flags ^ GC_ALL_EDGE_CLEAR ; - if ( gc_clear_flags & GC_L_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_L ; - } else - if ( gc_clear_flags & GC_B_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_B ; - } else - if ( gc_clear_flags & GC_R_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_R ; - } else - if ( gc_clear_flags & GC_T_EDGE_CLEAR ) { - p_gc->gc_flags |= GC_AIM_T ; - } else { - /* all edges have something on them, do full test */ - p_gc->gc_flags |= GC_AIM_C ; - } - p_gc++ ; - } -} - -int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps ) -{ -pGridRec p_gr ; -double slope, inv_slope ; - - if ( Near(ya, yb, eps) ) { - if ( Near(xa, xb, eps) ) { - /* edge is 0 length, so get rid of it */ - return( FALSE ) ; - } else { - /* horizontal line */ - slope = HUGE ; - inv_slope = 0.0 ; - } - } else { - if ( Near(xa, xb, eps) ) { - /* vertical line */ - slope = 0.0 ; - inv_slope = HUGE ; - } else { - slope = (xb-xa)/(yb-ya) ; - inv_slope = (yb-ya)/(xb-xa) ; - } - } - - p_gc->tot_edges++ ; - if ( p_gc->tot_edges <= 1 ) { - p_gc->gr = (pGridRec)malloc( sizeof(GridRec) ) ; - } else { - p_gc->gr = (pGridRec)realloc( p_gc->gr, - p_gc->tot_edges * sizeof(GridRec) ) ; - } - MALLOC_CHECK( p_gc->gr ) ; - p_gr = &p_gc->gr[p_gc->tot_edges-1] ; - - p_gr->slope = slope ; - p_gr->inv_slope = inv_slope ; - - p_gr->xa = xa ; - p_gr->ya = ya ; - if ( xa <= xb ) { - p_gr->minx = xa ; - p_gr->maxx = xb ; - } else { - p_gr->minx = xb ; - p_gr->maxx = xa ; - } - if ( ya <= yb ) { - p_gr->miny = ya ; - p_gr->maxy = yb ; - } else { - p_gr->miny = yb ; - p_gr->maxy = ya ; - } - - /* P2 - P1 */ - p_gr->ax = xb - xa ; - p_gr->ay = yb - ya ; - - return( TRUE ) ; -} - -/* Test point against grid and edges in the cell (if any). Algorithm: - * Check bounding box; if outside then return. - * Check cell point is inside; if simple inside or outside then return. - * Find which edge or corner is considered to be the best for testing and - * send a test ray towards it, counting the crossings. Add in the - * state of the edge or corner the ray went to and so determine the - * state of the point (inside or outside). - */ -int GridTest(pGridSet p_gs, pPipoint point) -{ -int j, count, init_flag ; -pGridCell p_gc ; -pGridRec p_gr ; -double tx, ty, xcell, ycell, bx,by,cx,cy, cornerx, cornery ; -double alpha, beta, denom ; -unsigned short gc_flags ; -int inside_flag ; - - /* first, is point inside bounding rectangle? */ - if ( ( ty = point->y ) < p_gs->miny || - ty >= p_gs->maxy || - ( tx = point->x ) < p_gs->minx || - tx >= p_gs->maxx ) { - - /* outside of box */ - inside_flag = FALSE ; - } else { - - /* what cell are we in? */ - ycell = ( ty - p_gs->miny ) * p_gs->inv_ydelta ; - xcell = ( tx - p_gs->minx ) * p_gs->inv_xdelta ; - p_gc = &p_gs->gc[((int)ycell)*p_gs->xres + (int)xcell] ; - - /* is cell simple? */ - count = p_gc->tot_edges ; - if ( count ) { - /* no, so find an edge which is free. */ - gc_flags = p_gc->gc_flags ; - p_gr = p_gc->gr ; - switch( gc_flags & GC_AIM ) { - case GC_AIM_L: - /* left edge is clear, shoot X- ray */ - /* note - this next statement requires that GC_BL_IN is 1 */ - inside_flag = gc_flags & GC_BL_IN ; - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if y is between edges */ - if ( ty >= p_gr->miny && ty < p_gr->maxy ) { - if ( tx > p_gr->maxx ) { - inside_flag = !inside_flag ; - } else if ( tx > p_gr->minx ) { - /* full computation */ - if ( ( p_gr->xa - - ( p_gr->ya - ty ) * p_gr->slope ) < tx ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_B: - /* bottom edge is clear, shoot Y+ ray */ - /* note - this next statement requires that GC_BL_IN is 1 */ - inside_flag = gc_flags & GC_BL_IN ; - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if x is between edges */ - if ( tx >= p_gr->minx && tx < p_gr->maxx ) { - if ( ty > p_gr->maxy ) { - inside_flag = !inside_flag ; - } else if ( ty > p_gr->miny ) { - /* full computation */ - if ( ( p_gr->ya - ( p_gr->xa - tx ) * - p_gr->inv_slope ) < ty ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_R: - /* right edge is clear, shoot X+ ray */ - inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ; - - /* TBD: Note, we could have sorted the edges to be tested - * by miny or somesuch, and so be able to cut testing - * short when the list's miny > point.y . - */ - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if y is between edges */ - if ( ty >= p_gr->miny && ty < p_gr->maxy ) { - if ( tx <= p_gr->minx ) { - inside_flag = !inside_flag ; - } else if ( tx <= p_gr->maxx ) { - /* full computation */ - if ( ( p_gr->xa - - ( p_gr->ya - ty ) * p_gr->slope ) >= tx ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_T: - /* top edge is clear, shoot Y+ ray */ - inside_flag = (gc_flags & GC_TR_IN) ? 1 : 0 ; - for ( j = count+1 ; --j ; p_gr++ ) { - /* test if x is between edges */ - if ( tx >= p_gr->minx && tx < p_gr->maxx ) { - if ( ty <= p_gr->miny ) { - inside_flag = !inside_flag ; - } else if ( ty <= p_gr->maxy ) { - /* full computation */ - if ( ( p_gr->ya - ( p_gr->xa - tx ) * - p_gr->inv_slope ) >= ty ) { - inside_flag = !inside_flag ; - } - } - } - } - break ; - - case GC_AIM_C: - /* no edge is clear, bite the bullet and test - * against the bottom left corner. - * We use Franklin Antonio's algorithm (Graphics Gems III). - */ - /* TBD: Faster yet might be to test against the closest - * corner to the cell location, but our hope is that we - * rarely need to do this testing at all. - */ - inside_flag = ((gc_flags & GC_BL_IN) == GC_BL_IN) ; - init_flag = TRUE ; - - /* get lower left corner coordinate */ - cornerx = p_gs->glx[(int)xcell] ; - cornery = p_gs->gly[(int)ycell] ; - for ( j = count+1 ; --j ; p_gr++ ) { - - /* quick out test: if test point is - * less than minx & miny, edge cannot overlap. - */ - if ( tx >= p_gr->minx && ty >= p_gr->miny ) { - - /* quick test failed, now check if test point and - * corner are on different sides of edge. - */ - if ( init_flag ) { - /* Compute these at most once for test */ - /* P3 - P4 */ - bx = tx - cornerx ; - by = ty - cornery ; - init_flag = FALSE ; - } - denom = p_gr->ay * bx - p_gr->ax * by ; - if ( denom != 0.0 ) { - /* lines are not collinear, so continue */ - /* P1 - P3 */ - cx = p_gr->xa - tx ; - cy = p_gr->ya - ty ; - alpha = by * cx - bx * cy ; - if ( denom > 0.0 ) { - if ( alpha < 0.0 || alpha >= denom ) { - /* test edge not hit */ - goto NextEdge ; - } - beta = p_gr->ax * cy - p_gr->ay * cx ; - if ( beta < 0.0 || beta >= denom ) { - /* polygon edge not hit */ - goto NextEdge ; - } - } else { - if ( alpha > 0.0 || alpha <= denom ) { - /* test edge not hit */ - goto NextEdge ; - } - beta = p_gr->ax * cy - p_gr->ay * cx ; - if ( beta > 0.0 || beta <= denom ) { - /* polygon edge not hit */ - goto NextEdge ; - } - } - inside_flag = !inside_flag ; - } - - } - NextEdge: ; - } - break ; - } - } else { - /* simple cell, so if lower left corner is in, - * then cell is inside. - */ - inside_flag = p_gc->gc_flags & GC_BL_IN ; - } - } - - return( inside_flag ) ; -} - -void GridCleanup(pGridSet p_gs) -{ -int i ; -pGridCell p_gc ; - - for ( i = 0, p_gc = p_gs->gc - ; i < p_gs->tot_cells - ; i++, p_gc++ ) { - - if ( p_gc->gr ) { - free( p_gc->gr ) ; - } - } - free( p_gs->glx ) ; - free( p_gs->gly ) ; - free( p_gs->gc ) ; -} diff --git a/thirdparty/ptinpoly.h b/thirdparty/ptinpoly.h deleted file mode 100644 index 7a0d26be..00000000 --- a/thirdparty/ptinpoly.h +++ /dev/null @@ -1,94 +0,0 @@ -/* ptinpoly.h - point in polygon inside/outside algorithms header file. - * - * by Eric Haines, 3D/Eye Inc, erich@eye.com - */ - -#ifdef __cplusplus -extern "C" -{ -#endif -/* =========================== System Related ============================= */ -/* SRAN initializes random number generator, if needed */ -#define SRAN() srand(1) -/* RAN01 returns a double from [0..1) */ -#define RAN01() ((double)rand() / 32768.0) - -/* =========== Grid stuff ================================================= */ - -#define GR_FULL_VERT 0x01 /* line crosses vertically */ -#define GR_FULL_HORZ 0x02 /* line crosses horizontally */ - -typedef struct { - double xa,ya ; - double minx, maxx, miny, maxy ; - double ax, ay ; - double slope, inv_slope ; -} GridRec, *pGridRec; - -#define GC_BL_IN 0x0001 /* bottom left corner is in (else out) */ -#define GC_BR_IN 0x0002 /* bottom right corner is in (else out) */ -#define GC_TL_IN 0x0004 /* top left corner is in (else out) */ -#define GC_TR_IN 0x0008 /* top right corner is in (else out) */ -#define GC_L_EDGE_HIT 0x0010 /* left edge is crossed */ -#define GC_R_EDGE_HIT 0x0020 /* right edge is crossed */ -#define GC_B_EDGE_HIT 0x0040 /* bottom edge is crossed */ -#define GC_T_EDGE_HIT 0x0080 /* top edge is crossed */ -#define GC_B_EDGE_PARITY 0x0100 /* bottom edge parity */ -#define GC_T_EDGE_PARITY 0x0200 /* top edge parity */ -#define GC_AIM_L (0<<10) /* aim towards left edge */ -#define GC_AIM_B (1<<10) /* aim towards bottom edge */ -#define GC_AIM_R (2<<10) /* aim towards right edge */ -#define GC_AIM_T (3<<10) /* aim towards top edge */ -#define GC_AIM_C (4<<10) /* aim towards a corner */ -#define GC_AIM 0x1c00 - -#define GC_L_EDGE_CLEAR GC_L_EDGE_HIT -#define GC_R_EDGE_CLEAR GC_R_EDGE_HIT -#define GC_B_EDGE_CLEAR GC_B_EDGE_HIT -#define GC_T_EDGE_CLEAR GC_T_EDGE_HIT - -#define GC_ALL_EDGE_CLEAR (GC_L_EDGE_HIT | \ - GC_R_EDGE_HIT | \ - GC_B_EDGE_HIT | \ - GC_T_EDGE_HIT ) - -typedef struct { - short tot_edges ; - unsigned short gc_flags ; - GridRec *gr ; -} GridCell, *pGridCell; - -typedef struct { - int xres, yres ; /* grid size */ - int tot_cells ; /* xres * yres */ - double minx, maxx, miny, maxy ; /* bounding box */ - double xdelta, ydelta ; - double inv_xdelta, inv_ydelta ; - double *glx, *gly ; - GridCell *gc ; -} GridSet, *pGridSet ; - -typedef struct { - double x, y; -} Pipoint, *pPipoint; - -/* add a little to the limits of the polygon bounding box to avoid precision -* problems. -*/ -#define EPSILON 0.00001 - -/* The following structure is associated with a polygon */ -typedef struct { - int id ; /* vertex number of edge */ - int full_cross ; /* 1 if extends from top to bottom */ - double minx, maxx ; /* X bounds for bin */ -} Edge, *pEdge ; - -void GridSetup(pPipoint pgon[], int numverts, int resolution, pGridSet p_gs); -int AddGridRecAlloc(pGridCell p_gc, double xa, double ya, double xb, double yb, double eps); -void GridCleanup(pGridSet p_gs); -int GridTest(pGridSet p_gs, pPipoint point); - -#ifdef __cplusplus -} -#endif \ No newline at end of file From 0b582273b5d3b16b1e4ef5c4c83bbecfbe546697 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 27 Dec 2018 14:13:09 +0100 Subject: [PATCH 13/53] merge all postgis (pdok + pdok citygml) into single function --- src/Map3d.cpp | 114 ++++++++------------------------------------------ src/Map3d.h | 3 +- src/main.cpp | 14 ++----- 3 files changed, 23 insertions(+), 108 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 5b7bd9e7..17215ac6 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -406,7 +406,7 @@ void Map3d::get_obj_per_class(std::wostream& of) { of << fs << std::endl; } -bool Map3d::get_pdok_output(std::string filename) { +bool Map3d::get_postgis_output(std::string connstr, bool pdok = false, bool citygml = false) { #if GDAL_VERSION_MAJOR < 2 std::cerr << "ERROR: cannot write MultiPolygonZ files with GDAL < 2.0.\n"; return false; @@ -414,7 +414,7 @@ bool Map3d::get_pdok_output(std::string filename) { if (GDALGetDriverCount() == 0) GDALAllRegister(); GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("PostgreSQL"); - GDALDataset* dataSource = driver->Create(filename.c_str(), 0, 0, 0, GDT_Unknown, NULL); + GDALDataset* dataSource = driver->Create(connstr.c_str(), 0, 0, 0, GDT_Unknown, NULL); if (dataSource == NULL) { std::cerr << "Starting database connection failed.\n"; return false; @@ -432,9 +432,9 @@ bool Map3d::get_pdok_output(std::string filename) { AttributeMap &attributes = f->get_attributes(); //Add additional attribute to list for layer creation attributes["xml"] = std::make_pair(OFTString, ""); - OGRLayer *layer = create_gdal_layer(driver, dataSource, filename, layername, attributes, f->get_class() == BUILDING); + OGRLayer *layer = create_gdal_layer(driver, dataSource, connstr, layername, attributes, f->get_class() == BUILDING); if (layer == NULL) { - std::cerr << "ERROR: Cannot open database '" + filename + "' for writing" << std::endl; + std::cerr << "ERROR: Cannot open database '" + connstr + "' for writing" << std::endl; dataSource->RollbackTransaction(); GDALClose(dataSource); GDALClose(driver); @@ -456,101 +456,23 @@ bool Map3d::get_pdok_output(std::string filename) { int i = 1; for (auto& f : _lsFeatures) { std::string layername = f->get_layername(); - //Add additional attribute describing CityGML of feature - std::wstring_convert>converter; - std::wstringstream ss; - ss << std::fixed << std::setprecision(3); - f->get_citygml_imgeo(ss); - std::string gmlAttribute = converter.to_bytes(ss.str()); - ss.clear(); AttributeMap extraAttribute = AttributeMap(); - extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute); - - if (!f->get_shape(layers[layername], true, extraAttribute)) { - return false; - } - if (i % 1000 == 0) { - if (dataSource->CommitTransaction() != OGRERR_NONE) { - std::cerr << "Writing to database failed.\n"; - return false; - } - if (dataSource->StartTransaction() != OGRERR_NONE) { - std::cerr << "Starting database transaction failed.\n"; - return false; + + if (pdok) { + //Add additional attribute describing CityGML of feature + std::wstring_convert>converter; + std::wstringstream ss; + ss << std::fixed << std::setprecision(3); + if (citygml) { + f->get_citygml(ss); } - } - i++; - } - if (dataSource->CommitTransaction() != OGRERR_NONE) { - std::cerr << "Writing to database failed.\n"; - return false; - } - GDALClose(dataSource); - GDALClose(driver); - return true; -#endif -} - -bool Map3d::get_pdok_citygml_output(std::string filename) { -#if GDAL_VERSION_MAJOR < 2 - std::cerr << "ERROR: cannot write MultiPolygonZ files with GDAL < 2.0.\n"; - return false; -#else - if (GDALGetDriverCount() == 0) - GDALAllRegister(); - GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("PostgreSQL"); - GDALDataset* dataSource = driver->Create(filename.c_str(), 0, 0, 0, GDT_Unknown, NULL); - if (dataSource == NULL) { - std::cerr << "Starting database connection failed.\n"; - return false; - } - if (dataSource->StartTransaction() != OGRERR_NONE) { - std::cerr << "Starting database transaction failed.\n"; - return false; - } - - std::unordered_map layers; - // create and write layers first - for (auto& f : _lsFeatures) { - std::string layername = f->get_layername(); - if (layers.find(layername) == layers.end()) { - AttributeMap &attributes = f->get_attributes(); - //Add additional attribute to list for layer creation - attributes["xml"] = std::make_pair(OFTString, ""); - OGRLayer *layer = create_gdal_layer(driver, dataSource, filename, layername, attributes, f->get_class() == BUILDING); - if (layer == NULL) { - std::cerr << "ERROR: Cannot open database '" + filename + "' for writing" << std::endl; - dataSource->RollbackTransaction(); - GDALClose(dataSource); - GDALClose(driver); - return false; + else { + f->get_citygml_imgeo(ss); } - layers.emplace(layername, layer); + std::string gmlAttribute = converter.to_bytes(ss.str()); + ss.clear(); + extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute); } - } - if (dataSource->CommitTransaction() != OGRERR_NONE) { - std::cerr << "Writing to database failed.\n"; - return false; - } - if (dataSource->StartTransaction() != OGRERR_NONE) { - std::cerr << "Starting database transaction failed.\n"; - return false; - } - - // create and write features to layers - int i = 1; - for (auto& f : _lsFeatures) { - std::string layername = f->get_layername(); - //Add additional attribute describing CityGML of feature - std::wstring_convert>converter; - std::wstringstream ss; - ss << std::fixed << std::setprecision(3); - f->get_citygml(ss); - std::string gmlAttribute = converter.to_bytes(ss.str()); - ss.clear(); - AttributeMap extraAttribute = AttributeMap(); - extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute); - if (!f->get_shape(layers[layername], true, extraAttribute)) { return false; } @@ -558,7 +480,7 @@ bool Map3d::get_pdok_citygml_output(std::string filename) { if (dataSource->CommitTransaction() != OGRERR_NONE) { std::cerr << "Writing to database failed.\n"; return false; - } + } if (dataSource->StartTransaction() != OGRERR_NONE) { std::cerr << "Starting database transaction failed.\n"; return false; diff --git a/src/Map3d.h b/src/Map3d.h index 006cd2b1..08d31a88 100644 --- a/src/Map3d.h +++ b/src/Map3d.h @@ -71,8 +71,7 @@ class Map3d { bool get_cityjson(std::string filename); void get_citygml_imgeo_multifile(std::string ofname); void create_citygml_imgeo_header(std::wostream& of); - bool get_pdok_output(std::string filename); - bool get_pdok_citygml_output(std::string filename); + bool get_postgis_output(std::string filename, bool pdok = false, bool citygml = false); bool get_gdal_output(std::string filename, std::string drivername, bool multi); void get_csv_buildings(std::wostream& of); void get_csv_buildings_multiple_heights(std::wostream& of); diff --git a/src/main.cpp b/src/main.cpp index 3134bc2a..5e7a3446 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -71,7 +71,6 @@ int main(int argc, const char * argv[]) { outputs["Shapefile"] = ""; outputs["Shapefile-Multifile"] = ""; outputs["PostGIS"] = ""; - outputs["PostGIS-Multi"] = ""; outputs["PostGIS-PDOK"] = ""; outputs["PostGIS-PDOK-CityGML"] = ""; outputs["GDAL"] = ""; @@ -96,7 +95,6 @@ int main(int argc, const char * argv[]) { ("Shapefile", po::value(&outputs["Shapefile"]), "Output ") ("Shapefile-Multifile", po::value(&outputs["Shapefile-Multifile"]), "Output ") ("PostGIS", po::value(&outputs["PostGIS"]), "Output ") - ("PostGIS-Multi", po::value(&outputs["PostGIS-Multi"]), "Output ") ("PostGIS-PDOK", po::value(&outputs["PostGIS-PDOK"]), "Output ") ("PostGIS-PDOK-CityGML", po::value(&outputs["PostGIS-PDOK-CityGML"]), "Output ") ("GDAL", po::value(&outputs["GDAL"]), "Output ") @@ -609,7 +607,7 @@ int main(int argc, const char * argv[]) { std::string ofname = output.second; if (format != "CityGML-Multifile" && format != "CityGML-IMGeo-Multifile" && format != "CityJSON" && format != "Shapefile" && format != "Shapefile-Multifile" && - format != "PostGIS" && format != "PostGIS-Multi" && format != "PostGIS-PDOK" && format != "PostGIS-PDOK-CityGML" && + format != "PostGIS" && format != "PostGIS-PDOK" && format != "PostGIS-PDOK-CityGML" && format != "GDAL") { of.open(ofname); } @@ -663,19 +661,15 @@ int main(int argc, const char * argv[]) { } else if (format == "PostGIS") { std::clog << "PostGIS output\n"; - fileWritten = map3d.get_gdal_output(ofname, "PostgreSQL", false); - } - else if (format == "PostGIS-Multi") { - std::clog << "PostGIS multiple table output\n"; - fileWritten = map3d.get_gdal_output(ofname, "PostgreSQL", true); + fileWritten = map3d.get_postgis_output(ofname, false, false); } else if (format == "PostGIS-PDOK") { std::clog << "PostGIS with IMGeo GML string output\n"; - fileWritten = map3d.get_pdok_output(ofname); + fileWritten = map3d.get_postgis_output(ofname, true, false); } else if (format == "PostGIS-PDOK-CityGML") { std::clog << "PostGIS with CityGML string output\n"; - fileWritten = map3d.get_pdok_citygml_output(ofname); + fileWritten = map3d.get_postgis_output(ofname, true, true); } else if (format == "GDAL") { //-- TODO: what is this? a path? how to use? if (nodes["output"] && nodes["output"]["gdal_driver"]) { From 7a7d586411f39e431b4f2d19d14726a678f12a1a Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 27 Dec 2018 14:15:33 +0100 Subject: [PATCH 14/53] remove redefinition of defaults --- src/Map3d.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 17215ac6..635edf1c 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -406,7 +406,7 @@ void Map3d::get_obj_per_class(std::wostream& of) { of << fs << std::endl; } -bool Map3d::get_postgis_output(std::string connstr, bool pdok = false, bool citygml = false) { +bool Map3d::get_postgis_output(std::string connstr, bool pdok, bool citygml) { #if GDAL_VERSION_MAJOR < 2 std::cerr << "ERROR: cannot write MultiPolygonZ files with GDAL < 2.0.\n"; return false; @@ -1733,4 +1733,3 @@ void Map3d::add_allowed_las_class(AllowedLASTopo c, int i) { void Map3d::add_allowed_las_class_within(AllowedLASTopo c, int i) { _las_classes_allowed_within[c].insert(i); } - From 7909882679a24f89b8d0dd02b03dc843c8d72245 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 27 Dec 2018 14:28:50 +0100 Subject: [PATCH 15/53] remove empty xml from attributes --- src/Map3d.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 635edf1c..49df33d2 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -430,8 +430,10 @@ bool Map3d::get_postgis_output(std::string connstr, bool pdok, bool citygml) { std::string layername = f->get_layername(); if (layers.find(layername) == layers.end()) { AttributeMap &attributes = f->get_attributes(); - //Add additional attribute to list for layer creation - attributes["xml"] = std::make_pair(OFTString, ""); + if (pdok) { + //Add additional attribute to list for layer creation + attributes["xml"] = std::make_pair(OFTString, ""); + } OGRLayer *layer = create_gdal_layer(driver, dataSource, connstr, layername, attributes, f->get_class() == BUILDING); if (layer == NULL) { std::cerr << "ERROR: Cannot open database '" + connstr + "' for writing" << std::endl; From 1cff0ec45a92a993e10a2851b22212ce597b96d6 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 27 Dec 2018 15:30:43 +0100 Subject: [PATCH 16/53] do use gdal_output for simple postgis writing for correct attribute names with dashes --- src/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.cpp b/src/main.cpp index 5e7a3446..ec8d153a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -661,7 +661,7 @@ int main(int argc, const char * argv[]) { } else if (format == "PostGIS") { std::clog << "PostGIS output\n"; - fileWritten = map3d.get_postgis_output(ofname, false, false); + fileWritten = map3d.get_gdal_output(ofname, "PostgreSQL", true); } else if (format == "PostGIS-PDOK") { std::clog << "PostGIS with IMGeo GML string output\n"; From 0c65b9df3d76e2364068bfe62073aa7e17cd1b94 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 27 Dec 2018 16:24:27 +0100 Subject: [PATCH 17/53] make writeAttribute check attribute by replacing '-' with '_' --- src/Map3d.cpp | 4 +++- src/TopoFeature.cpp | 9 +++++++-- src/main.cpp | 2 +- 3 files changed, 11 insertions(+), 4 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 49df33d2..57d0bb08 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -540,7 +540,9 @@ bool Map3d::get_gdal_output(std::string filename, std::string drivername, bool m } layers.emplace(layername, layer); } - f->get_shape(layers[layername], true); + if (!f->get_shape(layers[layername], true)) { + return false; + } } close_gdal_resources(driver, layers); } diff --git a/src/TopoFeature.cpp b/src/TopoFeature.cpp index 85db134d..164b2fa3 100644 --- a/src/TopoFeature.cpp +++ b/src/TopoFeature.cpp @@ -373,8 +373,13 @@ bool TopoFeature::get_multipolygon_features(OGRLayer* layer, std::string classNa bool TopoFeature::writeAttribute(OGRFeature* feature, OGRFeatureDefn* featureDefn, std::string name, std::string value) { int fi = featureDefn->GetFieldIndex(name.c_str()); if (fi == -1) { - std::cerr << "Failed to write attribute " << name << ".\n"; - return false; + // try replace '-' with '_' for postgresql column names + std::replace(name.begin(), name.end(), '-', '_'); + fi = featureDefn->GetFieldIndex(name.c_str()); + if (fi == -1) { + std::cerr << "Failed to write attribute " << name << ".\n"; + return false; + } } // perform extra character encoding for gdal. char* attrcpl = CPLRecode(value.c_str(), "", CPL_ENC_UTF8); diff --git a/src/main.cpp b/src/main.cpp index ec8d153a..5e7a3446 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -661,7 +661,7 @@ int main(int argc, const char * argv[]) { } else if (format == "PostGIS") { std::clog << "PostGIS output\n"; - fileWritten = map3d.get_gdal_output(ofname, "PostgreSQL", true); + fileWritten = map3d.get_postgis_output(ofname, false, false); } else if (format == "PostGIS-PDOK") { std::clog << "PostGIS with IMGeo GML string output\n"; From f5e49f39e4febfc0a08e6b9eed9191207a10a318 Mon Sep 17 00:00:00 2001 From: Ylannl Date: Tue, 8 Jan 2019 17:11:12 +0100 Subject: [PATCH 18/53] attempt to fix tinsimp crash with duplicate point --- src/geomtools.cpp | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index aee5fd99..23ec78c2 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -270,7 +270,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { for(int i=0; i &pts, double threshold) { std::vector faces; T.get_conflicts ( max_p, std::back_inserter(faces) ); + // handle case where max_p somehow coincides with a polygon vertex + if (faces.size()==0) { + // remove this points from the heap + heap.pop(); + + // this should return the already existing vertex + auto v = T.insert(max_p); + + // check the incident faces and erase any references to max_p + auto face_circulator = T.incident_faces(v); + auto start = face_circulator; + do { + auto face = *face_circulator; + + if (face->info().points_inside) { + for (auto it =*face->info().points_inside.begin(); it != *face->info().points_inside.end(); ){ + auto h = *it; + if( maxelement.index == (*h).index) + *face->info().points_inside.erase(it); + } + face->info().points_inside->clear(); + } + face_circulator++; + } while (face_circulator != start) + + // skip to next iteration + continue; + } + // insert max_p in triangulation auto face_hint = faces[0]; + auto v = T.insert(max_p, face_hint); face_hint = v->face(); From 945eddcfddfd85191392e1ef661ee7efc04a432d Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 10 Jan 2019 11:57:17 +0100 Subject: [PATCH 19/53] greedy_insert move heap.pop to end --- src/geomtools.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 23ec78c2..84e0a8be 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -292,9 +292,6 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // handle case where max_p somehow coincides with a polygon vertex if (faces.size()==0) { - // remove this points from the heap - heap.pop(); - // this should return the already existing vertex auto v = T.insert(max_p); @@ -304,17 +301,19 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { do { auto face = *face_circulator; - if (face->info().points_inside) { - for (auto it =*face->info().points_inside.begin(); it != *face->info().points_inside.end(); ){ + if (face.info().points_inside) { + for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); it++){ auto h = *it; if( maxelement.index == (*h).index) - *face->info().points_inside.erase(it); + face.info().points_inside->erase(it); } - face->info().points_inside->clear(); } face_circulator++; - } while (face_circulator != start) - + } while (face_circulator != start); + + // remove this points from the heap + heap.pop(); + // skip to next iteration continue; } From 9a2761f6f2f03199ffcfc2607fc483455b6d02df Mon Sep 17 00:00:00 2001 From: Ylannl Date: Thu, 10 Jan 2019 15:59:53 +0100 Subject: [PATCH 20/53] fix nasty duplicate point bug in tinsimp --- src/geomtools.cpp | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 84e0a8be..2b4f0303 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -52,6 +52,7 @@ struct point_error { }; typedef boost::heap::fibonacci_heap Heap; typedef Heap::handle_type heap_handle; +typedef std::vector heap_handle_vec; typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Projection_traits_xy_3 Gt; @@ -63,7 +64,7 @@ struct FaceInfo2 bool in_domain() { return nesting_level % 2 == 1; } - std::vector* points_inside = nullptr; + heap_handle_vec* points_inside = nullptr; CGAL::Plane_3* plane = nullptr; }; typedef CGAL::Triangulation_face_base_with_info_2 Fbb; @@ -245,7 +246,7 @@ inline double compute_error(Point &p, CDT::Face_handle &face) { face->vertex(1)->point(), face->vertex(2)->point()); if(!face->info().points_inside) - face->info().points_inside = new std::vector(); + face->info().points_inside = new heap_handle_vec(); auto plane = face->info().plane; auto interpolate = - plane->a()/plane->c() * p.x() - plane->b()/plane->c()*p.y() - plane->d()/plane->c(); @@ -296,20 +297,24 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { auto v = T.insert(max_p); // check the incident faces and erase any references to max_p - auto face_circulator = T.incident_faces(v); - auto start = face_circulator; + CDT::Face_circulator fcirculator = T.incident_faces(v), done(fcirculator); + auto start = fcirculator; do { - auto face = *face_circulator; - + auto face = *fcirculator; if (face.info().points_inside) { - for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); it++){ - auto h = *it; - if( maxelement.index == (*h).index) - face.info().points_inside->erase(it); + // collect the heap_handles that need to be removed + std::vector to_erase; + for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); ++it){ + if( maxelement.index == (**it).index) { + to_erase.push_back(it); + } + } + // remove the collected heap_handles + for (auto it : to_erase) { + face.info().points_inside->erase(it); } } - face_circulator++; - } while (face_circulator != start); + } while (++fcirculator != done); // remove this points from the heap heap.pop(); @@ -325,7 +330,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { face_hint = v->face(); // update clear info of triangles that just changed, collect points that were inside these triangles - std::vector points_to_update; + heap_handle_vec points_to_update; for (auto face : faces) { if (face->info().plane){ delete face->info().plane; From a4421ed249a32b344bb370a828fc368b675e88de Mon Sep 17 00:00:00 2001 From: Ylannl Date: Thu, 10 Jan 2019 16:28:28 +0100 Subject: [PATCH 21/53] make cmake attempt to find LASlib automatically --- CMakeLists.txt | 23 ++++++------ cmake/Modules/FindLASlib.cmake | 45 ++++++++++++++++++++++++ cmake/Modules/FindLASzip.cmake | 56 ----------------------------- cmake/Modules/FindlibLAS.cmake | 64 ---------------------------------- 4 files changed, 57 insertions(+), 131 deletions(-) create mode 100644 cmake/Modules/FindLASlib.cmake delete mode 100644 cmake/Modules/FindLASzip.cmake delete mode 100644 cmake/Modules/FindlibLAS.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index dd76ecaa..207e38d9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,9 +34,10 @@ if ( NOT GDAL_FOUND ) message(SEND_ERROR "3dfier requires the GDAL library") endif() -# LIBLAS -find_package(libLAS REQUIRED) -find_package (LASzip QUIET) +# LASLIB +# find_package(libLAS REQUIRED) +# find_package (LASzip QUIET) +find_package (LASlib REQUIRED) # YamlCpp find_package(YamlCpp REQUIRED) @@ -53,15 +54,15 @@ include( ${CGAL_USE_FILE} ) include_directories( ${CMAKE_SOURCE_DIR}/thirdparty ) -include_directories( ${CGAL_INCLUDE_DIR} ${CGAL_3RD_PARTY_INCLUDE_DIR} ${Boost_INCLUDE_DIRS} ${LIBLAS_INCLUDE_DIR} ${LASZIP_INCLUDE_DIR} ${YAMLCPP_INCLUDE_DIR} ${GDAL_INCLUDE_DIR}) +include_directories( ${CGAL_INCLUDE_DIR} ${CGAL_3RD_PARTY_INCLUDE_DIR} ${Boost_INCLUDE_DIRS} ${YAMLCPP_INCLUDE_DIR} ${GDAL_INCLUDE_DIR} ${LASLIB_INCLUDE_DIR}) link_directories(${YamlCpp_LIBRARY_DIRS}) -# build ptinpoly as static lib -add_library(ptinpoly STATIC thirdparty/ptinpoly.c) -set_target_properties( - ptinpoly - PROPERTIES C_STANDARD 11 -) +# # build ptinpoly as static lib +# add_library(ptinpoly STATIC thirdparty/ptinpoly.c) +# set_target_properties( +# ptinpoly +# PROPERTIES C_STANDARD 11 +# ) # Creating entries for target: 3dfier FILE(GLOB SRC_FILES src/*.cpp) @@ -71,6 +72,6 @@ set_target_properties( PROPERTIES CXX_STANDARD 11 ) -target_link_libraries( 3dfier ${CGAL_LIBRARIES} ${CGAL_3RD_PARTY_LIBRARIES} ${GDAL_LIBRARY} ${LIBLAS_LIBRARY} ${LASZIP_LIBRARY} ${YAMLCPP_LIBRARY} Boost::program_options Boost::filesystem Boost::locale ptinpoly) +target_link_libraries( 3dfier ${CGAL_LIBRARIES} ${CGAL_3RD_PARTY_LIBRARIES} ${GDAL_LIBRARY} ${YAMLCPP_LIBRARY} Boost::program_options Boost::filesystem Boost::locale ${LASLIB_LIBRARY}) install(TARGETS 3dfier DESTINATION bin) \ No newline at end of file diff --git a/cmake/Modules/FindLASlib.cmake b/cmake/Modules/FindLASlib.cmake new file mode 100644 index 00000000..fd9d22a1 --- /dev/null +++ b/cmake/Modules/FindLASlib.cmake @@ -0,0 +1,45 @@ +############################################################################### +# +# CMake module to search for LASlib library +# +# On success, the macro sets the following variables: +# LASLIB_FOUND = if the library found +# LASLIB_LIBRARY = full path to the library +# LASLIB_INCLUDE_DIR = where to find the library headers also defined, +# but not for general use are +# LASLIB_LIBRARY = where to find the library. +# LASLIB_VERSION = version of library which was found, e.g. "1.2.5" +# +############################################################################### +MESSAGE(STATUS "Searching for LASlib library") + +IF(LASLIB_INCLUDE_DIR) + # Already in cache, be silent + SET(LASLIB_FIND_QUIETLY TRUE) +ENDIF() + +# IF(APPLE) +# SET(LASLIB_ROOT_DIR /usr/local/) +# MESSAGE(STATUS "Trying LasZip using default location LASLIB_ROOT=c:/laszip") +# ENDIF() + +FIND_PATH(LASLIB_INCLUDE_DIR laszip.hpp + PATH_PREFIXES laszip + PATHS + /usr/include + /usr/local/include + /usr/local/include/LASlib + # ${LASLIB_ROOT_DIR}/include + NO_DEFAULT_PATH) + +SET(LASLIB_NAMES ${LASLIB_LIBRARY} LASlib) + +FIND_LIBRARY(LASLIB_LIBRARY + NAMES ${LASLIB_NAMES} + PATHS + /usr/lib + /usr/local/lib + /usr/local/lib/LASlib + ${LASLIB_ROOT_DIR}/lib) + + diff --git a/cmake/Modules/FindLASzip.cmake b/cmake/Modules/FindLASzip.cmake deleted file mode 100644 index 70d1a4be..00000000 --- a/cmake/Modules/FindLASzip.cmake +++ /dev/null @@ -1,56 +0,0 @@ -############################################################################### -# -# CMake module to search for LASzip library -# -# On success, the macro sets the following variables: -# LASZIP_FOUND = if the library found -# LASZIP_LIBRARIES = full path to the library -# LASZIP_INCLUDE_DIR = where to find the library headers also defined, -# but not for general use are -# LASZIP_LIBRARY = where to find the PROJ.4 library. -# LASZIP_VERSION = version of library which was found, e.g. "1.2.5" -# -# Copyright (c) 2009 Mateusz Loskot -# -# Module source: http://github.com/mloskot/workshop/tree/master/cmake/ -# -# Redistribution and use is allowed according to the terms of the BSD license. -# For details see the accompanying COPYING-CMAKE-SCRIPTS file. -# -############################################################################### -MESSAGE(STATUS "Searching for LASzip ${LASzip_FIND_VERSION}+ library") - -IF(LASZIP_INCLUDE_DIR) - # Already in cache, be silent - SET(LASZIP_FIND_QUIETLY TRUE) -ENDIF() - -IF(WIN32) - IF(DEFINED ENV{LASZIP_ROOT}) - SET(LASZIP_ROOT_DIR $ENV{LASZIP_ROOT}) - MESSAGE(STATUS "Trying LasZip using environment variable LASZIP_ROOT=$ENV{LASZIP_ROOT}") - ELSE() - SET(LASZIP_ROOT_DIR c:/laszip) - MESSAGE(STATUS "Trying LasZip using default location LASZIP_ROOT=c:/laszip") - ENDIF() -ENDIF() - - -FIND_PATH(LASZIP_INCLUDE_DIR laszip.hpp - PATH_PREFIXES laszip - PATHS - /usr/include - /usr/local/include - ${LASZIP_ROOT_DIR}/include - NO_DEFAULT_PATH) - -SET(LASZIP_NAMES ${LASZIP_LIBRARY} laszip) - -FIND_LIBRARY(LASZIP_LIBRARY - NAMES ${LASZIP_NAMES} - PATHS - /usr/lib - /usr/local/lib - ${LASZIP_ROOT_DIR}/lib) - - diff --git a/cmake/Modules/FindlibLAS.cmake b/cmake/Modules/FindlibLAS.cmake deleted file mode 100644 index 7d6ce712..00000000 --- a/cmake/Modules/FindlibLAS.cmake +++ /dev/null @@ -1,64 +0,0 @@ -############################################################################### -# -# CMake module to search for libLAS library -# -# On success, the macro sets the following variables: -# LIBLAS_FOUND = if the library found -# LIBLAS_LIBRARIES = full path to the library -# LIBLAS_INCLUDE_DIR = where to find the library headers also defined, -# but not for general use are -# LIBLAS_LIBRARY = where to find the PROJ.4 library. -# LIBLAS_VERSION = version of library which was found, e.g. "1.2.5" -# -# Copyright (c) 2009 Mateusz Loskot -# -# Module source: http://github.com/mloskot/workshop/tree/master/cmake/ -# -# Redistribution and use is allowed according to the terms of the BSD license. -# For details see the accompanying COPYING-CMAKE-SCRIPTS file. -# -############################################################################### -MESSAGE(STATUS "Searching for LibLAS ${LibLAS_FIND_VERSION}+ library") - -IF(LIBLAS_INCLUDE_DIR) - # Already in cache, be silent - SET(LIBLAS_FIND_QUIETLY TRUE) -ENDIF() - -IF(WIN32) - IF(DEFINED ENV{LIBLAS_ROOT}) - SET(LIBLAS_ROOT_DIR $ENV{LIBLAS_ROOT}) - #MESSAGE(STATUS " FindLibLAS: trying LIBLAS using environment variable LIBLAS_ROOT=$ENV{LIBLAS_ROOT}") - ELSE() - SET(LIBLAS_ROOT_DIR c:/liblas) - #MESSAGE(STATUS " FindLibLAS: trying LIBLAS using default location LIBLAS_ROOT=c:/liblas") - ENDIF() -ENDIF() - - -FIND_PATH(LIBLAS_INCLUDE_DIR - liblas.hpp - PATH_PREFIXES liblas - PATHS - /usr/include - /usr/local/include - /tmp/lasjunk/include - ${LIBLAS_ROOT_DIR}/include) - -find_library(LIBLAS_LIBRARY - NAMES liblas.dylib liblas - PATHS - /usr/lib - /usr/local/lib - /tmp/lasjunk/lib - ${LIBLAS_ROOT_DIR}/lib -) - - - - - -# Handle the QUIETLY and REQUIRED arguments and set LIBLAS_FOUND to TRUE -# if all listed variables are TRUE -#INCLUDE(FindPackageHandleStandardArgs) -#FIND_PACKAGE_HANDLE_STANDARD_ARGS(libLAS DEFAULT_MSG LIBLAS_LIBRARY LIBLAS_INCLUDE_DIR) \ No newline at end of file From a031fd2c5a26fd0c2d4e5456c3c60017e9f2daed Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 11:04:16 +0100 Subject: [PATCH 22/53] update .gitignore --- .gitignore | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index f1ebd88c..57b65da9 100644 --- a/.gitignore +++ b/.gitignore @@ -41,8 +41,13 @@ vs_build/.vs *.pdb *.vcxproj.user *.lib -x64/ -x86/ + +# Folders +x64 +x86 +vendor +_site +.sass-cache /resources/convert-to-simulation/simulatorready /resources/convert-to-simulation/out.off /resources/convert-to-simulation/bin/ From e3e2668d07b7070376921ebc2e8d185e1114762d Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 16:10:51 +0100 Subject: [PATCH 23/53] update to cgal 4.12 --- vs_build/3dfier.vcxproj | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vs_build/3dfier.vcxproj b/vs_build/3dfier.vcxproj index 4dec2b74..61fcf733 100644 --- a/vs_build/3dfier.vcxproj +++ b/vs_build/3dfier.vcxproj @@ -89,7 +89,7 @@ /machine:x64 %(AdditionalOptions) - gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib + gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.12.lib;libCGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib $(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) @@ -137,7 +137,7 @@ %(Filename)_p.c - gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.8.lib;libCGAL-vc140-mt-4.8.lib;libgmp-10.lib;libyaml-cppmd.lib + gdal_i.lib;LASlib.lib;libCGAL_Core-vc140-mt-4.12.lib;libCGAL-vc140-mt-4.12.lib;libgmp-10.lib;libyaml-cppmd.lib $(LASLIB_ROOT)\lib\Release;$(BOOST_LIBRARYDIR);$(OSGEO4W_ROOT)\lib;$(CGAL_DIR)\lib\;$(CGAL_DIR)\auxiliary\gmp\lib\;$(YAML-CPP_DIR)\vs_build\lib\;%(AdditionalLibraryDirectories) true %(IgnoreSpecificDefaultLibraries) From ea19be8a40cece923818d8c63af589615e465c37 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 16:11:04 +0100 Subject: [PATCH 24/53] fix tinsimp point insert on vertex --- src/geomtools.cpp | 71 ++++++++++++++++++++++++++++------------------- 1 file changed, 43 insertions(+), 28 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 2b4f0303..d6fcf245 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -75,21 +75,6 @@ typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; typedef CDT::Point Point; typedef CGAL::Polygon_2 Polygon_2; -struct PointXYHash { - std::size_t operator()(Point const& p) const noexcept { - std::size_t h1 = std::hash{}(p.x()); - std::size_t h2 = std::hash{}(p.y()); - return h1 ^ (h2 << 1); - } -}; -struct PointXYEqual { - bool operator()(Point const& p1, Point const& p2) const noexcept { - auto ex = p1.x() == p2.x(); - auto ey = p1.y() == p2.y(); - return ex && ey; - } -}; - inline double compute_error(Point &p, CDT::Face_handle &face); void greedy_insert(CDT &T, const std::vector &pts, double threshold); @@ -267,17 +252,29 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // compute initial point errors, build heap, store point indices in triangles { - std::unordered_set set; - for(int i=0; iinfo().points_inside->push_back(handle); } + else { + std::cout << "CDT insert; point location not in face but "; + if (lt == CDT::VERTEX) { + std::cout << "on vertex."; + } + else if (lt == CDT::OUTSIDE_CONVEX_HULL) { + std::cout << "outside convex hull."; + } + else if (lt == CDT::OUTSIDE_AFFINE_HULL) { + std::cout << "outside affine hull."; + } + std::cout << " Point; " << std::fixed << std::setprecision(3) << p << std::endl; + } } } @@ -289,7 +286,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // get triangles that will change after inserting this max_p std::vector faces; - T.get_conflicts ( max_p, std::back_inserter(faces) ); + T.get_conflicts(max_p, std::back_inserter(faces)); // handle case where max_p somehow coincides with a polygon vertex if (faces.size()==0) { @@ -351,11 +348,29 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // update the errors of affected elevation points for (auto curelement : points_to_update){ auto p = cpts[(*curelement).index]; - auto containing_face = T.locate(p, face_hint); - const double e = compute_error(p, containing_face); - const point_error new_pe = point_error((*curelement).index, e); - heap.update(curelement, new_pe); - containing_face->info().points_inside->push_back(curelement); + //auto containing_face = T.locate(p, face_hint); + CDT::Locate_type lt; + int li; + CDT::Face_handle containing_face = T.locate(p, lt, li, face_hint); + if (lt == CDT::EDGE || lt == CDT::FACE) { + const double e = compute_error(p, containing_face); + const point_error new_pe = point_error((*curelement).index, e); + heap.update(curelement, new_pe); + containing_face->info().points_inside->push_back(curelement); + } + else { + std::cout << "CDT update; point location not in face but "; + if (lt == CDT::VERTEX) { + std::cout << "on vertex."; + } + else if (lt == CDT::OUTSIDE_CONVEX_HULL) { + std::cout << "outside convex hull."; + } + else if (lt == CDT::OUTSIDE_AFFINE_HULL) { + std::cout << "outside affine hull."; + } + std::cout << " Point; " << std::fixed << std::setprecision(3) << p << std::endl; + } } } From 702e8aad81850b4a9a178b690040c720086a795b Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 16:29:30 +0100 Subject: [PATCH 25/53] undo merge --- src/Map3d.cpp | 74 +++++++++++++++------------------------------------ src/Water.cpp | 8 +++--- src/io.cpp | 2 +- 3 files changed, 26 insertions(+), 58 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index baafff2f..57d0bb08 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -458,13 +458,6 @@ bool Map3d::get_postgis_output(std::string connstr, bool pdok, bool citygml) { int i = 1; for (auto& f : _lsFeatures) { std::string layername = f->get_layername(); - //Add additional attribute describing CityGML of feature - std::wstring_convert>converter; - std::wstringstream ss; - ss << std::fixed << std::setprecision(3); - f->get_citygml_imgeo(ss); - std::string gmlAttribute = converter.to_bytes(ss.str()); - ss.clear(); AttributeMap extraAttribute = AttributeMap(); if (pdok) { @@ -482,30 +475,6 @@ bool Map3d::get_postgis_output(std::string connstr, bool pdok, bool citygml) { ss.clear(); extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute); } - } - if (dataSource->CommitTransaction() != OGRERR_NONE) { - std::cerr << "Writing to database failed.\n"; - return false; - } - if (dataSource->StartTransaction() != OGRERR_NONE) { - std::cerr << "Starting database transaction failed.\n"; - return false; - } - - // create and write features to layers - int i = 1; - for (auto& f : _lsFeatures) { - std::string layername = f->get_layername(); - //Add additional attribute describing CityGML of feature - std::wstring_convert>converter; - std::wstringstream ss; - ss << std::fixed << std::setprecision(3); - f->get_citygml(ss); - std::string gmlAttribute = converter.to_bytes(ss.str()); - ss.clear(); - AttributeMap extraAttribute = AttributeMap(); - extraAttribute["xml"] = std::make_pair(OFTString, gmlAttribute); - if (!f->get_shape(layers[layername], true, extraAttribute)) { return false; } @@ -1107,31 +1076,31 @@ void Map3d::extract_feature(OGRFeature *f, std::string layername, const char *id bool Map3d::add_las_file(PointFile pointFile) { std::clog << "Reading LAS/LAZ file: " << pointFile.filename << std::endl; - std::ifstream ifs; + + LASreadOpener lasreadopener; + lasreadopener.set_file_name(pointFile.filename.c_str()); + //-- set to compute bounding box + lasreadopener.set_populate_header(true); + LASreader* lasreader = lasreadopener.open(); try { - ifs.open(pointFile.filename.c_str(), std::ios::in | std::ios::binary); - if (ifs.is_open() == false) { + //-- check if file is open + if (lasreader == 0) { std::cerr << "\tERROR: could not open file: " << pointFile.filename << std::endl; return false; } - //-- LAS classes to omit - std::vector liblasomits; - for (int i : pointFile.lasomits) { - liblasomits.push_back(liblas::Classification(i)); - } + LASheader header = lasreader->header; - //-- read each point 1-by-1 - liblas::ReaderFactory f; - liblas::Reader reader = f.CreateWithStream(ifs); - liblas::Header const& header = reader.GetHeader(); + if (check_bounds(header.min_x, header.max_x, header.min_y, header.max_y)) { + //-- LAS classes to omit + std::vector lasomits; + for (int i : pointFile.lasomits) { + lasomits.push_back(i); + } - //-- check if the file overlaps the polygons - liblas::Bounds bounds = header.GetExtent(); - liblas::Bounds polygonBounds = get_bounds(); - uint32_t pointCount = header.GetPointRecordsCount(); + //-- read each point 1-by-1 + uint32_t pointCount = header.number_of_point_records; - if (polygonBounds.intersects(bounds)) { std::clog << "\t(" << boost::locale::as::number << pointCount << " points in the file)\n"; if ((pointFile.thinning > 1)) { std::clog << "\t(skipping every " << pointFile.thinning << "th points, thus "; @@ -1148,9 +1117,8 @@ bool Map3d::add_las_file(PointFile pointFile) { } printProgressBar(0); int i = 0; - - while (reader.ReadNextPoint()) { - liblas::Point const& p = reader.GetPoint(); + while (lasreader->read_point()) { + LASpoint const& p = lasreader->point; //-- set the thinning filter if (i % pointFile.thinning == 0) { //-- set the classification filter @@ -1171,11 +1139,11 @@ bool Map3d::add_las_file(PointFile pointFile) { else { std::clog << "\tskipping file, bounds do not intersect polygon extent\n"; } - ifs.close(); + lasreader->close(); } catch (std::exception e) { std::cerr << std::endl << e.what() << std::endl; - ifs.close(); + lasreader->close(); return false; } return true; diff --git a/src/Water.cpp b/src/Water.cpp index 41c82d12..4668ad68 100644 --- a/src/Water.cpp +++ b/src/Water.cpp @@ -108,10 +108,10 @@ void Water::get_citygml_imgeo(std::wostream& of) { std::string attribute; if (ondersteunend) { if (get_attribute("bgt-type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } else if (get_attribute("bgt_type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } if (get_attribute("plus-type", attribute)) { of << "" << attribute << ""; @@ -123,10 +123,10 @@ void Water::get_citygml_imgeo(std::wostream& of) { } else { if (get_attribute("bgt-type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } else if (get_attribute("bgt_type", attribute)) { - of << "" << attribute << ""; + of << "" << attribute << ""; } if (get_attribute("plus-type", attribute)) { of << "" << attribute << ""; diff --git a/src/io.cpp b/src/io.cpp index 00418927..ae608d79 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -136,7 +136,7 @@ void get_extruded_line_gml(std::wostream& of, Point2* a, Point2* b, double high, of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << low << std::setprecision(3) << " "; of << bg::get<0>(a) << " " << bg::get<1>(a) << " " << std::setprecision(2) << high << std::setprecision(3) << " "; of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << high << std::setprecision(3) << " "; - of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3); + of << bg::get<0>(b) << " " << bg::get<1>(b) << " " << std::setprecision(2) << low << std::setprecision(3); of << ""; of << ""; of << ""; From 2531ca40677211d397abe1d1caa0bd53b52858a1 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 16:40:33 +0100 Subject: [PATCH 26/53] code outlining --- src/geomtools.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index d6fcf245..a9e5f30e 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -67,13 +67,13 @@ struct FaceInfo2 heap_handle_vec* points_inside = nullptr; CGAL::Plane_3* plane = nullptr; }; -typedef CGAL::Triangulation_face_base_with_info_2 Fbb; -typedef CGAL::Constrained_triangulation_face_base_2 Fb; -typedef CGAL::Triangulation_data_structure_2 Tds; -typedef CGAL::Exact_predicates_tag Itag; -typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; -typedef CDT::Point Point; -typedef CGAL::Polygon_2 Polygon_2; +typedef CGAL::Triangulation_face_base_with_info_2 Fbb; +typedef CGAL::Constrained_triangulation_face_base_2 Fb; +typedef CGAL::Triangulation_data_structure_2 Tds; +typedef CGAL::Exact_predicates_tag Itag; +typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; +typedef CDT::Point Point; +typedef CGAL::Polygon_2 Polygon_2; inline double compute_error(Point &p, CDT::Face_handle &face); void greedy_insert(CDT &T, const std::vector &pts, double threshold); From 970694c21687720b8866e4196189cd3d58edc045 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 17:17:30 +0100 Subject: [PATCH 27/53] Update gitignore before merge --- .gitignore | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index f1ebd88c..57b65da9 100644 --- a/.gitignore +++ b/.gitignore @@ -41,8 +41,13 @@ vs_build/.vs *.pdb *.vcxproj.user *.lib -x64/ -x86/ + +# Folders +x64 +x86 +vendor +_site +.sass-cache /resources/convert-to-simulation/simulatorready /resources/convert-to-simulation/out.off /resources/convert-to-simulation/bin/ From f1607bec3966551d349c2e4d70de6a17ba01ff45 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 4 Mar 2019 17:27:02 +0100 Subject: [PATCH 28/53] update header text to year 2019 --- src/Bridge.cpp | 2 +- src/Bridge.h | 2 +- src/Building.cpp | 2 +- src/Building.h | 2 +- src/Forest.cpp | 2 +- src/Forest.h | 2 +- src/Map3d.cpp | 2 +- src/Map3d.h | 2 +- src/Road.cpp | 2 +- src/Road.h | 2 +- src/Separation.cpp | 2 +- src/Separation.h | 2 +- src/Terrain.cpp | 2 +- src/Terrain.h | 2 +- src/TopoFeature.cpp | 2 +- src/TopoFeature.h | 2 +- src/Water.cpp | 2 +- src/Water.h | 2 +- src/geomtools.cpp | 2 +- src/geomtools.h | 2 +- src/io.cpp | 2 +- src/io.h | 2 +- src/main.cpp | 8 ++++---- 23 files changed, 26 insertions(+), 26 deletions(-) diff --git a/src/Bridge.cpp b/src/Bridge.cpp index 5baef696..5d81e76d 100644 --- a/src/Bridge.cpp +++ b/src/Bridge.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Bridge.h b/src/Bridge.h index 1d0cc388..cf683b64 100644 --- a/src/Bridge.h +++ b/src/Bridge.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Building.cpp b/src/Building.cpp index c38eefed..d5a91dfa 100644 --- a/src/Building.cpp +++ b/src/Building.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Building.h b/src/Building.h index dc3215fd..159fb9dd 100644 --- a/src/Building.h +++ b/src/Building.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Forest.cpp b/src/Forest.cpp index 4c0e39b1..cff77662 100644 --- a/src/Forest.cpp +++ b/src/Forest.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Forest.h b/src/Forest.h index c031fe51..1cc3558f 100644 --- a/src/Forest.h +++ b/src/Forest.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 57d0bb08..0a557379 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Map3d.h b/src/Map3d.h index 08d31a88..1e5a03fb 100644 --- a/src/Map3d.h +++ b/src/Map3d.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Road.cpp b/src/Road.cpp index 1dcbfcbc..6a5e6c38 100644 --- a/src/Road.cpp +++ b/src/Road.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Road.h b/src/Road.h index 75ee95de..ec23520d 100644 --- a/src/Road.h +++ b/src/Road.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Separation.cpp b/src/Separation.cpp index 2f11428e..c2260b14 100644 --- a/src/Separation.cpp +++ b/src/Separation.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Separation.h b/src/Separation.h index a7f0b045..04e6e59f 100644 --- a/src/Separation.h +++ b/src/Separation.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Terrain.cpp b/src/Terrain.cpp index 264c3859..16dfb65c 100644 --- a/src/Terrain.cpp +++ b/src/Terrain.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Terrain.h b/src/Terrain.h index c04d6da6..f8483d9d 100644 --- a/src/Terrain.h +++ b/src/Terrain.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/TopoFeature.cpp b/src/TopoFeature.cpp index 164b2fa3..d8d1c464 100644 --- a/src/TopoFeature.cpp +++ b/src/TopoFeature.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/TopoFeature.h b/src/TopoFeature.h index 83ebb936..73a3aed9 100644 --- a/src/TopoFeature.h +++ b/src/TopoFeature.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Water.cpp b/src/Water.cpp index 4668ad68..c0e86070 100644 --- a/src/Water.cpp +++ b/src/Water.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/Water.h b/src/Water.h index 0edeb611..820c23b7 100644 --- a/src/Water.h +++ b/src/Water.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/geomtools.cpp b/src/geomtools.cpp index a9e5f30e..1ab2c932 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/geomtools.h b/src/geomtools.h index 7453ef24..a15fdc9f 100644 --- a/src/geomtools.h +++ b/src/geomtools.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/io.cpp b/src/io.cpp index ae608d79..9ae8580f 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/io.h b/src/io.h index b06bb92d..a66e54d7 100644 --- a/src/io.h +++ b/src/io.h @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. diff --git a/src/main.cpp b/src/main.cpp index 5e7a3446..0ae47216 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,7 +1,7 @@ /* 3dfier: takes 2D GIS datasets and "3dfies" to create 3D city models. - Copyright (C) 2015-2018 3D geoinformation research group, TU Delft + Copyright (C) 2015-2019 3D geoinformation research group, TU Delft This file is part of 3dfier. @@ -36,7 +36,7 @@ #include #include -std::string VERSION = "1.1.1"; +std::string VERSION = "1.1.2"; bool validate_yaml(const char* arg, std::set& allowedFeatures); int main(int argc, const char * argv[]); @@ -52,7 +52,7 @@ int main(int argc, const char * argv[]) { std::cout.imbue(loc); std::string licensewarning = - "3dfier Copyright (C) 2015-2018 3D geoinformation research group, TU Delft\n" + "3dfier Copyright (C) 2015-2019 3D geoinformation research group, TU Delft\n" "This program comes with ABSOLUTELY NO WARRANTY.\n" "This is free software, and you are welcome to redistribute it\n" "under certain conditions; for details run 3dfier with the '--license' option.\n"; @@ -698,7 +698,7 @@ std::string print_license() { std::string thelicense = "======================================================================\n" "\n3dfier: takes 2D GIS datasets and '3dfies' to create 3D city models.\n\n" - "Copyright (C) 2015-2018 3D geoinformation research group, TU Delft\n\n" + "Copyright (C) 2015-2019 3D geoinformation research group, TU Delft\n\n" "3dfier is free software: you can redistribute it and/or modify\n" "it under the terms of the GNU General Public License as published by\n" "the Free Software Foundation, either version 3 of the License, or\n" From 2eba67f7836a39b180c809f19a841ea6f9d1d139 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 5 Mar 2019 17:12:59 +0100 Subject: [PATCH 29/53] greedy_insert change iteration to find, only single element to remove --- src/geomtools.cpp | 14 +++----------- 1 file changed, 3 insertions(+), 11 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 1ab2c932..4e7a4742 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -299,17 +299,9 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { do { auto face = *fcirculator; if (face.info().points_inside) { - // collect the heap_handles that need to be removed - std::vector to_erase; - for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); ++it){ - if( maxelement.index == (**it).index) { - to_erase.push_back(it); - } - } - // remove the collected heap_handles - for (auto it : to_erase) { - face.info().points_inside->erase(it); - } + // remove the maxelement heap_handles + auto it = std::find(face.info().points_inside->begin(), face.info().points_inside->end(), maxelement); + face.info().points_inside->erase(it); } } while (++fcirculator != done); From 447018ee667d4ddafe8a3314c5ba9e40a5415bd7 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 5 Mar 2019 17:49:12 +0100 Subject: [PATCH 30/53] fix greedy insertion std::find and clear vector properly --- src/geomtools.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 4e7a4742..4f351c41 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -299,9 +299,17 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { do { auto face = *fcirculator; if (face.info().points_inside) { - // remove the maxelement heap_handles - auto it = std::find(face.info().points_inside->begin(), face.info().points_inside->end(), maxelement); - face.info().points_inside->erase(it); + // collect the heap_handles that need to be removed + std::vector to_erase; + for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); ++it) { + if (maxelement.index == (**it).index) { + to_erase.push_back(it); + } + } + // remove the collected heap_handles + for (auto it : to_erase) { + face.info().points_inside->erase(it); + } } } while (++fcirculator != done); @@ -330,7 +338,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { if( maxelement.index != (*h).index) points_to_update.push_back(h); } - face->info().points_inside->clear(); + heap_handle_vec().swap((*face->info().points_inside)); } } From 284656d97b19e15c4db7d0dc948094c310782d92 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 27 Mar 2019 11:50:26 +0100 Subject: [PATCH 31/53] introduce threshold_bridge_jump_edges --- resources/config_files/myconfig.yml | 1 + resources/config_files/myconfig_DEFAULTS.yml | 3 ++- resources/config_files/myconfig_README.yml | 1 + src/Map3d.cpp | 9 +++++++-- src/Map3d.h | 2 ++ src/main.cpp | 14 ++++++++++++++ 6 files changed, 27 insertions(+), 3 deletions(-) diff --git a/resources/config_files/myconfig.yml b/resources/config_files/myconfig.yml index 6f8de9eb..8eb99a64 100644 --- a/resources/config_files/myconfig.yml +++ b/resources/config_files/myconfig.yml @@ -96,3 +96,4 @@ options: building_radius_vertex_elevation: 3.0 radius_vertex_elevation: 1.0 threshold_jump_edges: 0.5 + threshold_bridge_jump_edges: 0.5 diff --git a/resources/config_files/myconfig_DEFAULTS.yml b/resources/config_files/myconfig_DEFAULTS.yml index d1704132..9fd1c091 100644 --- a/resources/config_files/myconfig_DEFAULTS.yml +++ b/resources/config_files/myconfig_DEFAULTS.yml @@ -36,4 +36,5 @@ input_elevation: options: building_radius_vertex_elevation: 3.0 radius_vertex_elevation: 1.0 - threshold_jump_edges: 0.5 \ No newline at end of file + threshold_jump_edges: 0.5 + threshold_bridge_jump_edges: 0.5 \ No newline at end of file diff --git a/resources/config_files/myconfig_README.yml b/resources/config_files/myconfig_README.yml index d7fa1ff1..f3dcc8df 100644 --- a/resources/config_files/myconfig_README.yml +++ b/resources/config_files/myconfig_README.yml @@ -90,4 +90,5 @@ options: # Global options building_radius_vertex_elevation: 3.0 # Radius in meters used for point-vertex distance between 3D points and vertices of building polygons, radius_vertex_elevation used when not specified radius_vertex_elevation: 1.0 # Radius in meters used for point-vertex distance between 3D points and vertices of polygons threshold_jump_edges: 0.5 # Threshold in meters for stitching adjacent objects, when the height difference is larger then the threshold a vertical wall is created + threshold_bridge_jump_edges: 0.5 # Threshold in meters for stitching bridges to adjacent objects, if not specified it falls back to threshold_jump_edges extent: xmin, ymin, xmax, ymax # Filter the input polygons to this extent diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 0a557379..f41aed65 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -51,6 +51,7 @@ Map3d::Map3d() { _radius_vertex_elevation = 1.0; _building_radius_vertex_elevation = 3.0; _threshold_jump_edges = 50; + _threshold_bridge_jump_edges = 50; _requestedExtent = Box2(Point2(0, 0), Point2(0, 0)); _bbox = Box2(Point2(999999, 999999), Point2(-999999, -999999)); _minxradius = 999999; @@ -83,6 +84,10 @@ void Map3d::set_threshold_jump_edges(float threshold) { _threshold_jump_edges = int(threshold * 100); } +void Map3d::set_threshold_bridge_jump_edges(float threshold) { + _threshold_bridge_jump_edges = int(threshold * 100); +} + void Map3d::set_building_include_floor(bool include) { _building_include_floor = include; } @@ -1533,7 +1538,7 @@ void Map3d::stitch_bridges() { pis.clear(); if (!(fadj->get_class() == BRIDGE && fadj->get_top_level()) && fadj->has_point2(ring[i], ringis, pis)) { int z = fadj->get_vertex_elevation(ringis[0], pis[0]); - if (abs(f->get_vertex_elevation(ringi, i) - z) < _threshold_jump_edges) { + if (abs(f->get_vertex_elevation(ringi, i) - z) < _threshold_bridge_jump_edges) { f->set_vertex_elevation(ringi, i, z); if (!(fadj->get_class() == BRIDGE && fadj->get_top_level() == f->get_top_level())) { // Add height to NC @@ -1705,7 +1710,7 @@ void Map3d::stitch_bridges() { int interz = interpolate_height(f, p, ringi, previ, ringi, endCorner.first); int prevz = f->get_vertex_elevation(ringi, previ); //Allways stich to lower object or if interpolated between corners within threshold or previous within threshold - if (stitchz < interz || abs(stitchz - interz) < _threshold_jump_edges || abs(stitchz - prevz) < _threshold_jump_edges) { + if (stitchz < interz || abs(stitchz - interz) < _threshold_bridge_jump_edges || abs(stitchz - prevz) < _threshold_bridge_jump_edges) { f->set_vertex_elevation(ringi, pi, stitchz); } else { diff --git a/src/Map3d.h b/src/Map3d.h index 1e5a03fb..1c720596 100644 --- a/src/Map3d.h +++ b/src/Map3d.h @@ -102,6 +102,7 @@ class Map3d { void set_radius_vertex_elevation(float radius); void set_building_radius_vertex_elevation(float radius); void set_threshold_jump_edges(float threshold); + void set_threshold_bridge_jump_edges(float threshold); void set_requested_extent(double xmin, double ymin, double xmax, double ymax); void add_allowed_las_class(AllowedLASTopo c, int i); @@ -132,6 +133,7 @@ class Map3d { float _radius_vertex_elevation; float _building_radius_vertex_elevation; int _threshold_jump_edges; //-- in cm/integer + int _threshold_bridge_jump_edges; //-- in cm/integer Box2 _bbox; double _minxradius; double _maxxradius; diff --git a/src/main.cpp b/src/main.cpp index 0ae47216..438ff26b 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -351,6 +351,11 @@ int main(int argc, const char * argv[]) { map3d.set_building_radius_vertex_elevation(n["building_radius_vertex_elevation"].as()); if (n["threshold_jump_edges"]) map3d.set_threshold_jump_edges(n["threshold_jump_edges"].as()); + if (n["threshold_bridge_jump_edges"]) + map3d.set_threshold_bridge_jump_edges(n["threshold_bridge_jump_edges"].as()); + else if (n["threshold_jump_edges"]) // set threshold_jump_edges same for bridge + map3d.set_threshold_bridge_jump_edges(n["threshold_jump_edges"].as()); + if (n["stitching"] && n["stitching"].as() == "false") bStitching = false; @@ -1143,6 +1148,15 @@ bool validate_yaml(const char* arg, std::set& allowedFeatures) { std::cerr << "\tOption 'options.threshold_jump_edges' invalid.\n"; } } + if (n["threshold_bridge_jump_edges"]) { + try { + boost::lexical_cast(n["threshold_bridge_jump_edges"].as()); + } + catch (boost::bad_lexical_cast& e) { + wentgood = false; + std::cerr << "\tOption 'options.threshold_bridge_jump_edges' invalid.\n"; + } + } if (n["stitching"]) { std::string s = n["stitching"].as(); if ((s != "true") && (s != "false")) { From 02e46f0a93f73c68c36966cae1a32ed115b6ad5c Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 27 Mar 2019 11:58:03 +0100 Subject: [PATCH 32/53] up version to v1.2 --- src/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.cpp b/src/main.cpp index 438ff26b..ba540a60 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -36,7 +36,7 @@ #include #include -std::string VERSION = "1.1.2"; +std::string VERSION = "1.2"; bool validate_yaml(const char* arg, std::set& allowedFeatures); int main(int argc, const char * argv[]); From 6fa54f1201eaab333ae4b80594c94e3eba545edd Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 11 Apr 2019 10:33:17 +0200 Subject: [PATCH 33/53] add setting for threshold_bridge_jump_edges --- src/Map3d.cpp | 9 +++++++-- src/Map3d.h | 2 ++ src/main.cpp | 4 ++++ 3 files changed, 13 insertions(+), 2 deletions(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index 0a557379..f41aed65 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -51,6 +51,7 @@ Map3d::Map3d() { _radius_vertex_elevation = 1.0; _building_radius_vertex_elevation = 3.0; _threshold_jump_edges = 50; + _threshold_bridge_jump_edges = 50; _requestedExtent = Box2(Point2(0, 0), Point2(0, 0)); _bbox = Box2(Point2(999999, 999999), Point2(-999999, -999999)); _minxradius = 999999; @@ -83,6 +84,10 @@ void Map3d::set_threshold_jump_edges(float threshold) { _threshold_jump_edges = int(threshold * 100); } +void Map3d::set_threshold_bridge_jump_edges(float threshold) { + _threshold_bridge_jump_edges = int(threshold * 100); +} + void Map3d::set_building_include_floor(bool include) { _building_include_floor = include; } @@ -1533,7 +1538,7 @@ void Map3d::stitch_bridges() { pis.clear(); if (!(fadj->get_class() == BRIDGE && fadj->get_top_level()) && fadj->has_point2(ring[i], ringis, pis)) { int z = fadj->get_vertex_elevation(ringis[0], pis[0]); - if (abs(f->get_vertex_elevation(ringi, i) - z) < _threshold_jump_edges) { + if (abs(f->get_vertex_elevation(ringi, i) - z) < _threshold_bridge_jump_edges) { f->set_vertex_elevation(ringi, i, z); if (!(fadj->get_class() == BRIDGE && fadj->get_top_level() == f->get_top_level())) { // Add height to NC @@ -1705,7 +1710,7 @@ void Map3d::stitch_bridges() { int interz = interpolate_height(f, p, ringi, previ, ringi, endCorner.first); int prevz = f->get_vertex_elevation(ringi, previ); //Allways stich to lower object or if interpolated between corners within threshold or previous within threshold - if (stitchz < interz || abs(stitchz - interz) < _threshold_jump_edges || abs(stitchz - prevz) < _threshold_jump_edges) { + if (stitchz < interz || abs(stitchz - interz) < _threshold_bridge_jump_edges || abs(stitchz - prevz) < _threshold_bridge_jump_edges) { f->set_vertex_elevation(ringi, pi, stitchz); } else { diff --git a/src/Map3d.h b/src/Map3d.h index 1e5a03fb..1c720596 100644 --- a/src/Map3d.h +++ b/src/Map3d.h @@ -102,6 +102,7 @@ class Map3d { void set_radius_vertex_elevation(float radius); void set_building_radius_vertex_elevation(float radius); void set_threshold_jump_edges(float threshold); + void set_threshold_bridge_jump_edges(float threshold); void set_requested_extent(double xmin, double ymin, double xmax, double ymax); void add_allowed_las_class(AllowedLASTopo c, int i); @@ -132,6 +133,7 @@ class Map3d { float _radius_vertex_elevation; float _building_radius_vertex_elevation; int _threshold_jump_edges; //-- in cm/integer + int _threshold_bridge_jump_edges; //-- in cm/integer Box2 _bbox; double _minxradius; double _maxxradius; diff --git a/src/main.cpp b/src/main.cpp index 0ae47216..d987168a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -351,6 +351,10 @@ int main(int argc, const char * argv[]) { map3d.set_building_radius_vertex_elevation(n["building_radius_vertex_elevation"].as()); if (n["threshold_jump_edges"]) map3d.set_threshold_jump_edges(n["threshold_jump_edges"].as()); + if (n["threshold_bridge_jump_edges"]) + map3d.set_threshold_bridge_jump_edges(n["threshold_bridge_jump_edges"].as()); + else if (n["threshold_jump_edges"]) + map3d.set_threshold_bridge_jump_edges(n["threshold_jump_edges"].as()); if (n["stitching"] && n["stitching"].as() == "false") bStitching = false; From 0b750c12059a1451791f55d699a303434ac15737 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 11 Apr 2019 10:50:50 +0200 Subject: [PATCH 34/53] remove newline --- src/main.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/main.cpp b/src/main.cpp index ba540a60..a22a249e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -355,7 +355,6 @@ int main(int argc, const char * argv[]) { map3d.set_threshold_bridge_jump_edges(n["threshold_bridge_jump_edges"].as()); else if (n["threshold_jump_edges"]) // set threshold_jump_edges same for bridge map3d.set_threshold_bridge_jump_edges(n["threshold_jump_edges"].as()); - if (n["stitching"] && n["stitching"].as() == "false") bStitching = false; From 0d50e582b594322b72fb8d599aab4e2f7da3b780 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 11 Apr 2019 16:48:56 +0200 Subject: [PATCH 35/53] add fibheap header from github --- thirdparty/FibHeap.hpp | 361 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 361 insertions(+) create mode 100644 thirdparty/FibHeap.hpp diff --git a/thirdparty/FibHeap.hpp b/thirdparty/FibHeap.hpp new file mode 100644 index 00000000..be5beb85 --- /dev/null +++ b/thirdparty/FibHeap.hpp @@ -0,0 +1,361 @@ +// +// FibHeap.hpp +// +// Created by Gabe Montague on 4/16/17. +// +// + +#ifndef FibHeap_hpp +#define FibHeap_hpp + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace FH { + + using std::logic_error; + using std::cout; + using std::vector; + using std::log; + + static const float PHI = 1.6180339887498948482; + + template + class FibHeap { + public: + + typedef std::size_t SizeType; + + // Define a type for a key-value pair + struct KVPair { + KeyType key; + ValueType value; + }; + + FibHeap() : _n(0), _max(nullptr) {} + + // TODO: copying + FibHeap(const FibHeap& that) = delete; + FibHeap& operator=(const FibHeap&) = delete; + + // Gets the number of items + SizeType getSize() const { + return _n; + } + + // Adds a new value to the heap + void insert(const KeyType key, const ValueType value) { + + Node* created = new Node(); + created->key = key; + created->value = value; + + if (_max == nullptr) { + + // Create a new root list + _max = created; + created->next = created; + created->prev = created; + + } else { + + // Insert into an existing root list + created->next = _max->next; + created->prev = _max; + _max->next->prev = created; + _max->next = created; + + // Check if we should change the min + if (created->key > _max->key) { + _max = created; + } + } + + // Increase the total number of elements + _n += 1; + } + + // Merges another Fibonacci heap to the current one, clearing the other heap + void combineWith(FibHeap & other) { + + _n += other._n; + + // Combine the root linked lists + Node* thisStart = _max; + Node* thisEnd = _max->prev; + Node* otherStart = other._max; + Node* otherEnd = other._max->prev; + + thisEnd->next = otherStart; + otherStart->prev = thisEnd; + thisStart->prev = otherEnd; + otherEnd->next = thisStart; + + // Update the minimum pointer + if (_max == nullptr) { + _max = other._max; // whether or not it is nullptr, other must be the min + } else if (other._max != nullptr && other._max->key > _max->key) { + _max = other._max; + } + + // Destroy the other heap + other._max = nullptr; + other._n = 0; + } + + // Peeks at the new value + KVPair peekMax() const { + return _max->getKV(); + } + + // Pops the value from the heap, removing and returning it + KVPair popMax() { + + if (_max == nullptr) { + throw logic_error("Attempted to pop item from empty heap"); + } + + const KVPair result = _max->getKV(); + const bool isLastElement = _max == _max->next && _max->child == nullptr; + + // Handle case of one element + if (isLastElement) { + delete _max; + _max = nullptr; + } else { + + // Add children of _max to the root list + if (_max->child != nullptr) { + Node* firstChild = _max->child; + Node* lastChild = _max->child->prev; + Node* insertionLower = _max; + Node* insertionUpper = _max->next; + + insertionLower->next = firstChild; + insertionUpper->prev = lastChild; + firstChild->prev = insertionLower; + lastChild->next = insertionUpper; + + // Iterate over added children and set parents to null + firstChild->parent = nullptr; + Node* iteratedNode = firstChild->next; + while (iteratedNode != firstChild) { + iteratedNode->parent = nullptr; + iteratedNode = iteratedNode->next; + } + } + + // Now remove the old min from the root list, and set the new min to the first child added (temporarily). + _max->prev->next = _max->next; + _max->next->prev = _max->prev; + const Node* const toDelete = _max; + _max = _max->next; + delete toDelete; + + //cout << "\nPRECONSOLIDATE\n"; + //this->print(); + + // Consolidate nodes + _consolidate(); + } + + _n -= 1; + + return result; + } + + // Run tests + static void test(); + + // Print + void print() const { + + cout << "\nFibHeap of " << _n << " nodes:\n"; + _max->print(); + cout << "\n"; + } + + // Destructor + virtual ~FibHeap() { + if (_max != nullptr) { + _max->deleteTree(); + } + } + + private: + + // Number of elements + SizeType _n; + + class Node { + public: + KeyType key; + ValueType value; + SizeType degree; + + Node* next; + Node* prev; + Node* parent; + Node* child; + + Node() : degree(0), next(nullptr), prev(nullptr), parent(nullptr), child(nullptr) {} + + KVPair getKV() const { + KVPair kv; + kv.key = key; + kv.value = value; + return kv; + } + + void print(const Node* start=nullptr) const { + + // Base case: no more siblings + if (this == start) { + return; + } + + // Print key and children + cout << "(" << key << ",{"; + if (child != nullptr) { + child->print(); + } + cout << "}),"; + + // Print siblings + next->print(start == nullptr ? this : start); + } + + void deleteTree(const Node* start=nullptr) { + + // Base case: no more siblings - delete self + if (this == start) { + delete this; + return; + } + + // Delete children + if (child != nullptr) { + child->deleteTree(); + } + + // Delete siblings + next->deleteTree(start == nullptr ? this : start); + } + }; + + // Pointer to the minimum element. Don't throw away the value or you leak the entire heap. + Node* _max; + + void _consolidate() { + + // The non-inclusive upper bound of the degree of any tree + const SizeType maxDegree = static_cast(floor(log(_n) / log(PHI))) + 1; + + // Allocate a vector with space maxDegree - all null. + auto a = vector(maxDegree, nullptr); + + // Iterate over the root list + const Node* const last = _max->prev; + Node* iteratedNode = _max; + + while (true) { + + Node* const iteratedNodeNext = iteratedNode->next; + + SizeType degree = iteratedNode->degree; + + // See if there is an existing tree we have already iterated over with matching degree + Node* smaller = iteratedNode; + while (a[degree] != nullptr) { + + // This will have the same degree as the iterated node + Node* otherNode = a[degree]; // Called y in notes + + const bool iteratedIsSmaller = smaller->key >= otherNode->key; + Node* const smallerSave = smaller; + smaller = iteratedIsSmaller ? smaller : otherNode; + Node* larger = iteratedIsSmaller ? otherNode : smallerSave; + + // Parent the smaller under the larger: first remove larger from the root list + larger->prev->next = larger->next; + larger->next->prev = larger->prev; + + // Now parent the larger under the smaller + larger->parent = smaller; + if (smaller->child == nullptr) { + larger->next = larger; + larger->prev = larger; + smaller->child = larger; + } else { + Node* sibling = smaller->child; + larger->next = sibling; + larger->prev = sibling->prev; + sibling->prev->next = larger; + sibling->prev = larger; + smaller->child = larger; + } + smaller->degree++; + + a[degree] = nullptr; + + // Now set the while loop to search for trees of the new degree + degree++; + } + + // Record this tree in the array as having degree degree + a[degree] = smaller; + + // Break out of the loop if we have exhausted nodes in our original root list + if (iteratedNode == last) { + break; + } + + iteratedNode = iteratedNodeNext; + } + + _max = nullptr; // We are safe to do this because currently all nodes are stored in the degree array + + // Recalculate the min, and reform connections (necessary?) + for (int i = 0; i < maxDegree; i++) { + + Node* node = a[i]; + + // Skip null entries in the degree array + if (node == nullptr) { + continue; + } + + // Create a new root list from a single node + if (_max == nullptr) { + node->parent = nullptr; + node->next = node; + node->prev = node; + _max = node; + + // ...or grow the new root list if it's been made + } else { + + // Insert node + node->next = _max->next; + node->prev = _max; + _max->next->prev = node; + _max->next = node; + + // Update the min if necessary + if (node->key > _max->key) { + _max = node; + } + } + } + } + }; +} + +#endif /* FibHeap_hpp */ \ No newline at end of file From 07ed1c6e392b5ac9400ac5a5e04241d8549cbe3f Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 11 Apr 2019 17:42:20 +0200 Subject: [PATCH 36/53] this FibHeap doesn't have update function --- thirdparty/FibHeap.hpp | 361 ----------------------------------------- 1 file changed, 361 deletions(-) delete mode 100644 thirdparty/FibHeap.hpp diff --git a/thirdparty/FibHeap.hpp b/thirdparty/FibHeap.hpp deleted file mode 100644 index be5beb85..00000000 --- a/thirdparty/FibHeap.hpp +++ /dev/null @@ -1,361 +0,0 @@ -// -// FibHeap.hpp -// -// Created by Gabe Montague on 4/16/17. -// -// - -#ifndef FibHeap_hpp -#define FibHeap_hpp - -#include - -#include -#include -#include -#include -#include -#include -#include - -namespace FH { - - using std::logic_error; - using std::cout; - using std::vector; - using std::log; - - static const float PHI = 1.6180339887498948482; - - template - class FibHeap { - public: - - typedef std::size_t SizeType; - - // Define a type for a key-value pair - struct KVPair { - KeyType key; - ValueType value; - }; - - FibHeap() : _n(0), _max(nullptr) {} - - // TODO: copying - FibHeap(const FibHeap& that) = delete; - FibHeap& operator=(const FibHeap&) = delete; - - // Gets the number of items - SizeType getSize() const { - return _n; - } - - // Adds a new value to the heap - void insert(const KeyType key, const ValueType value) { - - Node* created = new Node(); - created->key = key; - created->value = value; - - if (_max == nullptr) { - - // Create a new root list - _max = created; - created->next = created; - created->prev = created; - - } else { - - // Insert into an existing root list - created->next = _max->next; - created->prev = _max; - _max->next->prev = created; - _max->next = created; - - // Check if we should change the min - if (created->key > _max->key) { - _max = created; - } - } - - // Increase the total number of elements - _n += 1; - } - - // Merges another Fibonacci heap to the current one, clearing the other heap - void combineWith(FibHeap & other) { - - _n += other._n; - - // Combine the root linked lists - Node* thisStart = _max; - Node* thisEnd = _max->prev; - Node* otherStart = other._max; - Node* otherEnd = other._max->prev; - - thisEnd->next = otherStart; - otherStart->prev = thisEnd; - thisStart->prev = otherEnd; - otherEnd->next = thisStart; - - // Update the minimum pointer - if (_max == nullptr) { - _max = other._max; // whether or not it is nullptr, other must be the min - } else if (other._max != nullptr && other._max->key > _max->key) { - _max = other._max; - } - - // Destroy the other heap - other._max = nullptr; - other._n = 0; - } - - // Peeks at the new value - KVPair peekMax() const { - return _max->getKV(); - } - - // Pops the value from the heap, removing and returning it - KVPair popMax() { - - if (_max == nullptr) { - throw logic_error("Attempted to pop item from empty heap"); - } - - const KVPair result = _max->getKV(); - const bool isLastElement = _max == _max->next && _max->child == nullptr; - - // Handle case of one element - if (isLastElement) { - delete _max; - _max = nullptr; - } else { - - // Add children of _max to the root list - if (_max->child != nullptr) { - Node* firstChild = _max->child; - Node* lastChild = _max->child->prev; - Node* insertionLower = _max; - Node* insertionUpper = _max->next; - - insertionLower->next = firstChild; - insertionUpper->prev = lastChild; - firstChild->prev = insertionLower; - lastChild->next = insertionUpper; - - // Iterate over added children and set parents to null - firstChild->parent = nullptr; - Node* iteratedNode = firstChild->next; - while (iteratedNode != firstChild) { - iteratedNode->parent = nullptr; - iteratedNode = iteratedNode->next; - } - } - - // Now remove the old min from the root list, and set the new min to the first child added (temporarily). - _max->prev->next = _max->next; - _max->next->prev = _max->prev; - const Node* const toDelete = _max; - _max = _max->next; - delete toDelete; - - //cout << "\nPRECONSOLIDATE\n"; - //this->print(); - - // Consolidate nodes - _consolidate(); - } - - _n -= 1; - - return result; - } - - // Run tests - static void test(); - - // Print - void print() const { - - cout << "\nFibHeap of " << _n << " nodes:\n"; - _max->print(); - cout << "\n"; - } - - // Destructor - virtual ~FibHeap() { - if (_max != nullptr) { - _max->deleteTree(); - } - } - - private: - - // Number of elements - SizeType _n; - - class Node { - public: - KeyType key; - ValueType value; - SizeType degree; - - Node* next; - Node* prev; - Node* parent; - Node* child; - - Node() : degree(0), next(nullptr), prev(nullptr), parent(nullptr), child(nullptr) {} - - KVPair getKV() const { - KVPair kv; - kv.key = key; - kv.value = value; - return kv; - } - - void print(const Node* start=nullptr) const { - - // Base case: no more siblings - if (this == start) { - return; - } - - // Print key and children - cout << "(" << key << ",{"; - if (child != nullptr) { - child->print(); - } - cout << "}),"; - - // Print siblings - next->print(start == nullptr ? this : start); - } - - void deleteTree(const Node* start=nullptr) { - - // Base case: no more siblings - delete self - if (this == start) { - delete this; - return; - } - - // Delete children - if (child != nullptr) { - child->deleteTree(); - } - - // Delete siblings - next->deleteTree(start == nullptr ? this : start); - } - }; - - // Pointer to the minimum element. Don't throw away the value or you leak the entire heap. - Node* _max; - - void _consolidate() { - - // The non-inclusive upper bound of the degree of any tree - const SizeType maxDegree = static_cast(floor(log(_n) / log(PHI))) + 1; - - // Allocate a vector with space maxDegree - all null. - auto a = vector(maxDegree, nullptr); - - // Iterate over the root list - const Node* const last = _max->prev; - Node* iteratedNode = _max; - - while (true) { - - Node* const iteratedNodeNext = iteratedNode->next; - - SizeType degree = iteratedNode->degree; - - // See if there is an existing tree we have already iterated over with matching degree - Node* smaller = iteratedNode; - while (a[degree] != nullptr) { - - // This will have the same degree as the iterated node - Node* otherNode = a[degree]; // Called y in notes - - const bool iteratedIsSmaller = smaller->key >= otherNode->key; - Node* const smallerSave = smaller; - smaller = iteratedIsSmaller ? smaller : otherNode; - Node* larger = iteratedIsSmaller ? otherNode : smallerSave; - - // Parent the smaller under the larger: first remove larger from the root list - larger->prev->next = larger->next; - larger->next->prev = larger->prev; - - // Now parent the larger under the smaller - larger->parent = smaller; - if (smaller->child == nullptr) { - larger->next = larger; - larger->prev = larger; - smaller->child = larger; - } else { - Node* sibling = smaller->child; - larger->next = sibling; - larger->prev = sibling->prev; - sibling->prev->next = larger; - sibling->prev = larger; - smaller->child = larger; - } - smaller->degree++; - - a[degree] = nullptr; - - // Now set the while loop to search for trees of the new degree - degree++; - } - - // Record this tree in the array as having degree degree - a[degree] = smaller; - - // Break out of the loop if we have exhausted nodes in our original root list - if (iteratedNode == last) { - break; - } - - iteratedNode = iteratedNodeNext; - } - - _max = nullptr; // We are safe to do this because currently all nodes are stored in the degree array - - // Recalculate the min, and reform connections (necessary?) - for (int i = 0; i < maxDegree; i++) { - - Node* node = a[i]; - - // Skip null entries in the degree array - if (node == nullptr) { - continue; - } - - // Create a new root list from a single node - if (_max == nullptr) { - node->parent = nullptr; - node->next = node; - node->prev = node; - _max = node; - - // ...or grow the new root list if it's been made - } else { - - // Insert node - node->next = _max->next; - node->prev = _max; - _max->next->prev = node; - _max->next = node; - - // Update the min if necessary - if (node->key > _max->key) { - _max = node; - } - } - } - } - }; -} - -#endif /* FibHeap_hpp */ \ No newline at end of file From 92225d325bbe62f7f443d2ca95ca8c99ab711e46 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 18 Apr 2019 13:49:36 +0200 Subject: [PATCH 37/53] working with debug text --- src/geomtools.cpp | 89 ++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 81 insertions(+), 8 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 4f351c41..d077f65e 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -75,6 +75,21 @@ typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; typedef CDT::Point Point; typedef CGAL::Polygon_2 Polygon_2; +struct PointXYHash { + std::size_t operator()(Point const& p) const noexcept { + std::size_t h1 = std::hash{}(p.x()); + std::size_t h2 = std::hash{}(p.y()); + return h1 ^ (h2 << 1); + } +}; +struct PointXYEqual { + bool operator()(Point const& p1, Point const& p2) const noexcept { + auto ex = p1.x() == p2.x(); + auto ey = p1.y() == p2.y(); + return ex && ey; + } +}; + inline double compute_error(Point &p, CDT::Face_handle &face); void greedy_insert(CDT &T, const std::vector &pts, double threshold); @@ -244,11 +259,14 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { Heap heap; // Convert all elevation points to CGAL points + std::unordered_set cpts_set; std::vector cpts; cpts.reserve(pts.size()); for (auto& p : pts) { - cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); + if(cpts_set.insert(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))).second); } + cpts.assign(cpts_set.begin(), cpts_set.end()); + std::unordered_set().swap(cpts_set); // compute initial point errors, build heap, store point indices in triangles { @@ -325,36 +343,91 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { auto v = T.insert(max_p, face_hint); face_hint = v->face(); + bool maxe = false; // update clear info of triangles that just changed, collect points that were inside these triangles heap_handle_vec points_to_update; for (auto face : faces) { - if (face->info().plane){ + if (face->info().points_inside) { + for (auto h : *face->info().points_inside) { + if (maxelement.index == (*h).index) { + std::cout << "maxelement " << maxelement.index << " found " << std::endl; + std::cout << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; + maxe = true; + } + } + } + } + if (!maxe) { + std::cout << "maxelement " << maxelement.index << " not found " << std::endl; + } + + for (auto face : faces) { + if (face->info().plane) { delete face->info().plane; face->info().plane = nullptr; } if (face->info().points_inside) { - for (auto h :*face->info().points_inside){ - if( maxelement.index != (*h).index) + for (auto h : *face->info().points_inside) { + if (maxelement.index != (*h).index) { points_to_update.push_back(h); + + + //TESTING + std::cout << "pushed: " << (*h).index << std::endl; + + int index = (*h).index; + if (index < 0 || index >= cpts.size()) { + std::cout << "looking for index: " << maxelement.index << std::endl; + std::cout << "index out of range: " << index << std::endl; + if (!maxe) { + // Locate the face where + CDT::Locate_type lt; + int li; + CDT::Face_handle containing_face = T.locate(cpts[maxelement.index], lt, li, face_hint); + std::cout << "CDT update; point location: "; + if (lt == CDT::EDGE) { + std::cout << "on edge."; + } + else if (lt == CDT::FACE) { + std::cout << "in face."; + } + else if (lt == CDT::VERTEX) { + std::cout << "on vertex."; + } + else if (lt == CDT::OUTSIDE_CONVEX_HULL) { + std::cout << "outside convex hull."; + } + else if (lt == CDT::OUTSIDE_AFFINE_HULL) { + std::cout << "outside affine hull."; + } + std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; + } + } + } } heap_handle_vec().swap((*face->info().points_inside)); } } - - // remove the point we just inserted in the triangulation from the heap + + // remove this points from the heap heap.pop(); // update the errors of affected elevation points for (auto curelement : points_to_update){ - auto p = cpts[(*curelement).index]; + int index = (*curelement).index; + std::cout << "pulled: " << (*curelement).index << std::endl; + if (index > cpts.size()) { + std::cout << "index out of range: " << index << std::endl; + } + auto p = cpts[index]; //auto containing_face = T.locate(p, face_hint); CDT::Locate_type lt; int li; CDT::Face_handle containing_face = T.locate(p, lt, li, face_hint); if (lt == CDT::EDGE || lt == CDT::FACE) { const double e = compute_error(p, containing_face); - const point_error new_pe = point_error((*curelement).index, e); + const point_error new_pe = point_error(index, e); heap.update(curelement, new_pe); containing_face->info().points_inside->push_back(curelement); } From 2efef6d70958d505f2e7e42110a6f5e3bfa10264 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Thu, 18 Apr 2019 14:26:38 +0200 Subject: [PATCH 38/53] wip, still no solution --- src/geomtools.cpp | 35 +++++++++++++++++++++++++++++------ 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 4f351c41..a85a5116 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -75,6 +75,21 @@ typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; typedef CDT::Point Point; typedef CGAL::Polygon_2 Polygon_2; +struct PointXYHash { + std::size_t operator()(Point const& p) const noexcept { + std::size_t h1 = std::hash{}(p.x()); + std::size_t h2 = std::hash{}(p.y()); + return h1 ^ (h2 << 1); + } +}; +struct PointXYEqual { + bool operator()(Point const& p1, Point const& p2) const noexcept { + auto ex = p1.x() == p2.x(); + auto ey = p1.y() == p2.y(); + return ex && ey; + } +}; + inline double compute_error(Point &p, CDT::Face_handle &face); void greedy_insert(CDT &T, const std::vector &pts, double threshold); @@ -243,15 +258,22 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // assumes all lidar points are inside a triangle Heap heap; - // Convert all elevation points to CGAL points std::vector cpts; cpts.reserve(pts.size()); - for (auto& p : pts) { - cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); - } // compute initial point errors, build heap, store point indices in triangles { + // Convert all elevation points to CGAL points + std::unordered_set cpts_set; + + for (auto& p : pts) { + cpts_set.insert(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); + } + cpts.assign(cpts_set.begin(), cpts_set.end()); + if (cpts_set.size() != cpts.size()) { + std::cout << "size differs; " << cpts_set.size() << "&" << cpts.size(); + } + for (int i = 0; i < cpts.size(); i++) { auto p = cpts[i]; CDT::Locate_type lt; @@ -347,14 +369,15 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // update the errors of affected elevation points for (auto curelement : points_to_update){ - auto p = cpts[(*curelement).index]; + int idx = (*curelement).index; + auto p = cpts[idx]; //auto containing_face = T.locate(p, face_hint); CDT::Locate_type lt; int li; CDT::Face_handle containing_face = T.locate(p, lt, li, face_hint); if (lt == CDT::EDGE || lt == CDT::FACE) { const double e = compute_error(p, containing_face); - const point_error new_pe = point_error((*curelement).index, e); + const point_error new_pe = point_error(idx, e); heap.update(curelement, new_pe); containing_face->info().points_inside->push_back(curelement); } From c16e37725504ebfc8ce97e076624619f77e249bd Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 24 Apr 2019 15:13:00 +0200 Subject: [PATCH 39/53] swap points_to_update loop and insert statement --- src/geomtools.cpp | 112 +++++++++++++++++++++++++--------------------- 1 file changed, 62 insertions(+), 50 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index d077f65e..e7216073 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -84,9 +84,11 @@ struct PointXYHash { }; struct PointXYEqual { bool operator()(Point const& p1, Point const& p2) const noexcept { - auto ex = p1.x() == p2.x(); - auto ey = p1.y() == p2.y(); - return ex && ey; + return CGAL::compare_xy(p1, p2) == CGAL::Comparison_result::EQUAL; + + //auto ex = p1.x() == p2.x(); + //auto ey = p1.y() == p2.y(); + //return ex && ey; } }; @@ -263,7 +265,10 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { std::vector cpts; cpts.reserve(pts.size()); for (auto& p : pts) { - if(cpts_set.insert(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))).second); + auto r = cpts_set.insert(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); + if (!r.second) { + std::cout << "Duplicate point: " << *(r.first) << " && " << Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p)) << std::endl; + } } cpts.assign(cpts_set.begin(), cpts_set.end()); std::unordered_set().swap(cpts_set); @@ -338,30 +343,11 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { continue; } - // insert max_p in triangulation - auto face_hint = faces[0]; - - auto v = T.insert(max_p, face_hint); - face_hint = v->face(); - bool maxe = false; - // update clear info of triangles that just changed, collect points that were inside these triangles heap_handle_vec points_to_update; - for (auto face : faces) { - if (face->info().points_inside) { - for (auto h : *face->info().points_inside) { - if (maxelement.index == (*h).index) { - std::cout << "maxelement " << maxelement.index << " found " << std::endl; - std::cout << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; - maxe = true; - } - } - } - } - if (!maxe) { - std::cout << "maxelement " << maxelement.index << " not found " << std::endl; - } + bool maxe = false; + int pointless = 0; for (auto face : faces) { if (face->info().plane) { delete face->info().plane; @@ -374,41 +360,67 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { //TESTING - std::cout << "pushed: " << (*h).index << std::endl; + //std::cout << "pushed: " << (*h).index << std::endl; int index = (*h).index; if (index < 0 || index >= cpts.size()) { std::cout << "looking for index: " << maxelement.index << std::endl; std::cout << "index out of range: " << index << std::endl; - if (!maxe) { - // Locate the face where - CDT::Locate_type lt; - int li; - CDT::Face_handle containing_face = T.locate(cpts[maxelement.index], lt, li, face_hint); - std::cout << "CDT update; point location: "; - if (lt == CDT::EDGE) { - std::cout << "on edge."; - } - else if (lt == CDT::FACE) { - std::cout << "in face."; - } - else if (lt == CDT::VERTEX) { - std::cout << "on vertex."; - } - else if (lt == CDT::OUTSIDE_CONVEX_HULL) { - std::cout << "outside convex hull."; - } - else if (lt == CDT::OUTSIDE_AFFINE_HULL) { - std::cout << "outside affine hull."; - } - std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; - } } } + else { + maxe = true; + } } heap_handle_vec().swap((*face->info().points_inside)); } + else { + pointless++; + } } + std::cout << "Pointsless faces: " << pointless << std::endl; + if (!maxe) { + std::cout << "maxelement " << maxelement.index << " not found"; + + // Locate the face where + CDT::Locate_type lt; + int li; + CDT::Face_handle containing_face = T.locate(cpts[maxelement.index], lt, li, faces[0]); + std::cout << "; point location: "; + if (lt == CDT::EDGE) { + std::cout << "on edge."; + } + else if (lt == CDT::FACE) { + std::cout << "in face."; + } + else if (lt == CDT::VERTEX) { + std::cout << "on vertex."; + } + else if (lt == CDT::OUTSIDE_CONVEX_HULL) { + std::cout << "outside convex hull."; + } + else if (lt == CDT::OUTSIDE_AFFINE_HULL) { + std::cout << "outside affine hull."; + } + std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; + std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; + + heap_handle_vec* pi = containing_face->info().points_inside; + // Delete stray handle + for (auto h = pi->begin(); h != pi->end(); h++) { + //for (auto h : pi) { + if (maxelement.index == (**h).index) { + containing_face->info().points_inside->erase(h); + break; + } + } + } + + // insert max_p in triangulation + auto face_hint = faces[0]; + + auto v = T.insert(max_p, face_hint); + face_hint = v->face(); // remove this points from the heap heap.pop(); @@ -416,7 +428,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // update the errors of affected elevation points for (auto curelement : points_to_update){ int index = (*curelement).index; - std::cout << "pulled: " << (*curelement).index << std::endl; + //std::cout << "pulled: " << (*curelement).index << std::endl; if (index > cpts.size()) { std::cout << "index out of range: " << index << std::endl; } From 01de54dfe9c803fd39e7a603b8b08f19d76db631 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Wed, 24 Apr 2019 15:33:25 +0200 Subject: [PATCH 40/53] comment all debug --- src/geomtools.cpp | 103 +++++++++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 52 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index e7216073..87114a26 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -362,63 +362,62 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { //TESTING //std::cout << "pushed: " << (*h).index << std::endl; - int index = (*h).index; - if (index < 0 || index >= cpts.size()) { - std::cout << "looking for index: " << maxelement.index << std::endl; - std::cout << "index out of range: " << index << std::endl; - } - } - else { - maxe = true; + //int index = (*h).index; + //if (index < 0 || index >= cpts.size()) { + // std::cout << "looking for index: " << maxelement.index << std::endl; + // std::cout << "index out of range: " << index << std::endl; + //} } + //else { + // maxe = true; + //} } heap_handle_vec().swap((*face->info().points_inside)); } - else { - pointless++; - } - } - std::cout << "Pointsless faces: " << pointless << std::endl; - if (!maxe) { - std::cout << "maxelement " << maxelement.index << " not found"; - - // Locate the face where - CDT::Locate_type lt; - int li; - CDT::Face_handle containing_face = T.locate(cpts[maxelement.index], lt, li, faces[0]); - std::cout << "; point location: "; - if (lt == CDT::EDGE) { - std::cout << "on edge."; - } - else if (lt == CDT::FACE) { - std::cout << "in face."; - } - else if (lt == CDT::VERTEX) { - std::cout << "on vertex."; - } - else if (lt == CDT::OUTSIDE_CONVEX_HULL) { - std::cout << "outside convex hull."; - } - else if (lt == CDT::OUTSIDE_AFFINE_HULL) { - std::cout << "outside affine hull."; - } - std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; - std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; - - heap_handle_vec* pi = containing_face->info().points_inside; - // Delete stray handle - for (auto h = pi->begin(); h != pi->end(); h++) { - //for (auto h : pi) { - if (maxelement.index == (**h).index) { - containing_face->info().points_inside->erase(h); - break; - } - } + //else { + // pointless++; + //} } + //std::cout << "Pointsless faces: " << pointless << std::endl; + //if (!maxe) { + // std::cout << "maxelement " << maxelement.index << " not found"; + + // // Locate the face where + // CDT::Locate_type lt; + // int li; + // CDT::Face_handle containing_face = T.locate(cpts[maxelement.index], lt, li, faces[0]); + // std::cout << "; point location: "; + // if (lt == CDT::EDGE) { + // std::cout << "on edge."; + // } + // else if (lt == CDT::FACE) { + // std::cout << "in face."; + // } + // else if (lt == CDT::VERTEX) { + // std::cout << "on vertex."; + // } + // else if (lt == CDT::OUTSIDE_CONVEX_HULL) { + // std::cout << "outside convex hull."; + // } + // else if (lt == CDT::OUTSIDE_AFFINE_HULL) { + // std::cout << "outside affine hull."; + // } + // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; + // std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; + + // heap_handle_vec* pi = containing_face->info().points_inside; + // // Delete stray handle + // for (auto h = pi->begin(); h != pi->end(); h++) { + // //for (auto h : pi) { + // if (maxelement.index == (**h).index) { + // containing_face->info().points_inside->erase(h); + // break; + // } + // } + //} // insert max_p in triangulation auto face_hint = faces[0]; - auto v = T.insert(max_p, face_hint); face_hint = v->face(); @@ -429,9 +428,9 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { for (auto curelement : points_to_update){ int index = (*curelement).index; //std::cout << "pulled: " << (*curelement).index << std::endl; - if (index > cpts.size()) { - std::cout << "index out of range: " << index << std::endl; - } + //if (index > cpts.size()) { + // std::cout << "index out of range: " << index << std::endl; + //} auto p = cpts[index]; //auto containing_face = T.locate(p, face_hint); CDT::Locate_type lt; From 97fa2edc28fcb6f71467a8a927cbe6365ab24012 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 12:08:54 +0200 Subject: [PATCH 41/53] change cdt error message --- src/Map3d.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Map3d.cpp b/src/Map3d.cpp index f41aed65..aef63fd9 100644 --- a/src/Map3d.cpp +++ b/src/Map3d.cpp @@ -841,7 +841,7 @@ bool Map3d::construct_CDT() { p->buildCDT(); } catch (std::exception e) { - std::cerr << std::endl << "CDT failed for " << p->get_id() << " (" << p->get_class() << ") with error: " << e.what() << std::endl; + std::cerr << std::endl << "CDT failed for object \'" << p->get_id() << "\' (class " << p->get_class() << ") with error: " << e.what() << std::endl; return false; } } From 52182f1a12823a253948966dbee0b767fd139bf9 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 12:43:03 +0200 Subject: [PATCH 42/53] change maxelement.index to maxelementindex --- src/geomtools.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 87114a26..7a4a0562 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -304,8 +304,8 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // insert points, update errors of affected triangles until threshold error is reached while (!heap.empty() && heap.top().error > threshold){ // get top element (with largest error) from heap - auto maxelement = heap.top(); - auto max_p = cpts[maxelement.index]; + auto maxelementindex = heap.top().index; + auto max_p = cpts[maxelementindex]; // get triangles that will change after inserting this max_p std::vector faces; @@ -325,7 +325,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // collect the heap_handles that need to be removed std::vector to_erase; for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); ++it) { - if (maxelement.index == (**it).index) { + if (maxelementindex == (**it).index) { to_erase.push_back(it); } } @@ -355,7 +355,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { } if (face->info().points_inside) { for (auto h : *face->info().points_inside) { - if (maxelement.index != (*h).index) { + if (maxelementindex != (*h).index) { points_to_update.push_back(h); @@ -364,7 +364,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { //int index = (*h).index; //if (index < 0 || index >= cpts.size()) { - // std::cout << "looking for index: " << maxelement.index << std::endl; + // std::cout << "looking for index: " << maxelementindex << std::endl; // std::cout << "index out of range: " << index << std::endl; //} } @@ -380,12 +380,12 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { } //std::cout << "Pointsless faces: " << pointless << std::endl; //if (!maxe) { - // std::cout << "maxelement " << maxelement.index << " not found"; + // std::cout << "maxelement " << maxelementindex << " not found"; // // Locate the face where // CDT::Locate_type lt; // int li; - // CDT::Face_handle containing_face = T.locate(cpts[maxelement.index], lt, li, faces[0]); + // CDT::Face_handle containing_face = T.locate(cpts[maxelementindex], lt, li, faces[0]); // std::cout << "; point location: "; // if (lt == CDT::EDGE) { // std::cout << "on edge."; @@ -402,14 +402,14 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // else if (lt == CDT::OUTSIDE_AFFINE_HULL) { // std::cout << "outside affine hull."; // } - // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelement.index] << std::endl; + // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelementindex] << std::endl; // std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; // heap_handle_vec* pi = containing_face->info().points_inside; // // Delete stray handle // for (auto h = pi->begin(); h != pi->end(); h++) { // //for (auto h : pi) { - // if (maxelement.index == (**h).index) { + // if (maxelementindex == (**h).index) { // containing_face->info().points_inside->erase(h); // break; // } From d58c73bb54ebfd761a6058060de5a3b9f8181ba2 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 12:46:13 +0200 Subject: [PATCH 43/53] change cpts accessor to .at() --- src/geomtools.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 7a4a0562..16738ad8 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -276,7 +276,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // compute initial point errors, build heap, store point indices in triangles { for (int i = 0; i < cpts.size(); i++) { - auto p = cpts[i]; + auto p = cpts.at(i); CDT::Locate_type lt; int li; CDT::Face_handle face = T.locate(p, lt, li); @@ -305,7 +305,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { while (!heap.empty() && heap.top().error > threshold){ // get top element (with largest error) from heap auto maxelementindex = heap.top().index; - auto max_p = cpts[maxelementindex]; + auto max_p = cpts.at(maxelementindex); // get triangles that will change after inserting this max_p std::vector faces; @@ -385,7 +385,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // // Locate the face where // CDT::Locate_type lt; // int li; - // CDT::Face_handle containing_face = T.locate(cpts[maxelementindex], lt, li, faces[0]); + // CDT::Face_handle containing_face = T.locate(cpts.at(maxelementindex), lt, li, faces[0]); // std::cout << "; point location: "; // if (lt == CDT::EDGE) { // std::cout << "on edge."; @@ -402,7 +402,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // else if (lt == CDT::OUTSIDE_AFFINE_HULL) { // std::cout << "outside affine hull."; // } - // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts[maxelementindex] << std::endl; + // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts.at(maxelementindex) << std::endl; // std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; // heap_handle_vec* pi = containing_face->info().points_inside; @@ -431,7 +431,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { //if (index > cpts.size()) { // std::cout << "index out of range: " << index << std::endl; //} - auto p = cpts[index]; + auto p = cpts.at(index); //auto containing_face = T.locate(p, face_hint); CDT::Locate_type lt; int li; From 8ba7d7a8430f2a8e2a10f0f2db1ae83d3aeca654 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 13:23:27 +0200 Subject: [PATCH 44/53] write index with fail point in greedy_insert --- src/geomtools.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 16738ad8..6b9469e8 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -453,7 +453,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { else if (lt == CDT::OUTSIDE_AFFINE_HULL) { std::cout << "outside affine hull."; } - std::cout << " Point; " << std::fixed << std::setprecision(3) << p << std::endl; + std::cout << " Point (index " << index << "); " << std::fixed << std::setprecision(3) << p << std::endl; } } } From a0ea5ba19a6fcd245516da5e3e4c5b182f265909 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 16:34:57 +0200 Subject: [PATCH 45/53] fixed error, adding points in CDT on polygon constraints are not allowed!!! --- src/geomtools.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 6b9469e8..148c3a12 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -280,7 +280,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { CDT::Locate_type lt; int li; CDT::Face_handle face = T.locate(p, lt, li); - if (lt == CDT::EDGE || lt == CDT::FACE) { + if (lt == CDT::FACE) { auto e = compute_error(p, face); auto handle = heap.push(point_error(i, e)); face->info().points_inside->push_back(handle); @@ -290,6 +290,9 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { if (lt == CDT::VERTEX) { std::cout << "on vertex."; } + else if (lt == CDT::EDGE) { + std::cout << "on edge."; + } else if (lt == CDT::OUTSIDE_CONVEX_HULL) { std::cout << "outside convex hull."; } From e4d629db4335d15f907f1721a2628e8d3ebc7286 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 16:46:34 +0200 Subject: [PATCH 46/53] greedy_insert working version --- src/geomtools.cpp | 32 +++++++++++--------------------- 1 file changed, 11 insertions(+), 21 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 148c3a12..e8c16237 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -261,17 +261,11 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { Heap heap; // Convert all elevation points to CGAL points - std::unordered_set cpts_set; std::vector cpts; cpts.reserve(pts.size()); for (auto& p : pts) { - auto r = cpts_set.insert(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); - if (!r.second) { - std::cout << "Duplicate point: " << *(r.first) << " && " << Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p)) << std::endl; - } + cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); } - cpts.assign(cpts_set.begin(), cpts_set.end()); - std::unordered_set().swap(cpts_set); // compute initial point errors, build heap, store point indices in triangles { @@ -307,7 +301,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // insert points, update errors of affected triangles until threshold error is reached while (!heap.empty() && heap.top().error > threshold){ // get top element (with largest error) from heap - auto maxelementindex = heap.top().index; + int maxelementindex = heap.top().index; auto max_p = cpts.at(maxelementindex); // get triangles that will change after inserting this max_p @@ -349,14 +343,14 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // update clear info of triangles that just changed, collect points that were inside these triangles heap_handle_vec points_to_update; - bool maxe = false; - int pointless = 0; + //bool maxe = false; + //int pointless = 0; for (auto face : faces) { if (face->info().plane) { delete face->info().plane; face->info().plane = nullptr; } - if (face->info().points_inside) { + if (face->info().points_inside && face->info().points_inside->size() > 0) { for (auto h : *face->info().points_inside) { if (maxelementindex != (*h).index) { points_to_update.push_back(h); @@ -382,6 +376,10 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { //} } //std::cout << "Pointsless faces: " << pointless << std::endl; + //if (pointless == faces.size()) { + // std::cout << "Non of the located faces contain any points, index: " << maxelementindex << std::endl; + //} + //if (!maxe) { // std::cout << "maxelement " << maxelementindex << " not found"; @@ -407,16 +405,6 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // } // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts.at(maxelementindex) << std::endl; // std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; - - // heap_handle_vec* pi = containing_face->info().points_inside; - // // Delete stray handle - // for (auto h = pi->begin(); h != pi->end(); h++) { - // //for (auto h : pi) { - // if (maxelementindex == (**h).index) { - // containing_face->info().points_inside->erase(h); - // break; - // } - // } //} // insert max_p in triangulation @@ -449,6 +437,8 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { std::cout << "CDT update; point location not in face but "; if (lt == CDT::VERTEX) { std::cout << "on vertex."; + // update error to 0.0 to disable point + heap.update(curelement, point_error(index, 0.0)); } else if (lt == CDT::OUTSIDE_CONVEX_HULL) { std::cout << "outside convex hull."; From 77ca5990894faad9898267ad660bc13a8f2af3e0 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 16:47:25 +0200 Subject: [PATCH 47/53] clean set hash code --- src/geomtools.cpp | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index e8c16237..b8bbb152 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -75,23 +75,6 @@ typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; typedef CDT::Point Point; typedef CGAL::Polygon_2 Polygon_2; -struct PointXYHash { - std::size_t operator()(Point const& p) const noexcept { - std::size_t h1 = std::hash{}(p.x()); - std::size_t h2 = std::hash{}(p.y()); - return h1 ^ (h2 << 1); - } -}; -struct PointXYEqual { - bool operator()(Point const& p1, Point const& p2) const noexcept { - return CGAL::compare_xy(p1, p2) == CGAL::Comparison_result::EQUAL; - - //auto ex = p1.x() == p2.x(); - //auto ey = p1.y() == p2.y(); - //return ex && ey; - } -}; - inline double compute_error(Point &p, CDT::Face_handle &face); void greedy_insert(CDT &T, const std::vector &pts, double threshold); From e70bcfe6a05e91b2445d1875571b7816a715387e Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Mon, 29 Apr 2019 17:33:09 +0200 Subject: [PATCH 48/53] greedy_insert; stop using cpts vector and add points to point_error struct --- src/geomtools.cpp | 59 +++++++++++++++++++++++------------------------ 1 file changed, 29 insertions(+), 30 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index b8bbb152..c5f775a9 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -39,11 +39,14 @@ #include #include +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + // fibonacci heap for greedy insertion code struct point_error { - point_error(int i, double e) : index(i), error(e){} + point_error(int i, double e, CGAL::Point_3 p) : index(i), error(e), point(p){} int index; double error; + CGAL::Point_3 point; bool operator<(point_error const & rhs) const { @@ -54,7 +57,6 @@ typedef boost::heap::fibonacci_heap Heap; typedef Heap::handle_type heap_handle; typedef std::vector heap_handle_vec; -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; typedef CGAL::Projection_traits_xy_3 Gt; typedef CGAL::Triangulation_vertex_base_with_id_2 Vb; struct FaceInfo2 @@ -243,23 +245,23 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // assumes all lidar points are inside a triangle Heap heap; - // Convert all elevation points to CGAL points - std::vector cpts; - cpts.reserve(pts.size()); - for (auto& p : pts) { - cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); - } + //// Convert all elevation points to CGAL points + //std::vector cpts; + //cpts.reserve(pts.size()); + //for (auto& p : pts) { + // cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); + //} // compute initial point errors, build heap, store point indices in triangles { - for (int i = 0; i < cpts.size(); i++) { - auto p = cpts.at(i); + for (int i = 0; i < pts.size(); i++) { + auto p = Point(bg::get<0>(pts[i]), bg::get<1>(pts[i]), bg::get<2>(pts[i])); CDT::Locate_type lt; int li; CDT::Face_handle face = T.locate(p, lt, li); if (lt == CDT::FACE) { auto e = compute_error(p, face); - auto handle = heap.push(point_error(i, e)); + auto handle = heap.push(point_error(i, e, p)); face->info().points_inside->push_back(handle); } else { @@ -284,8 +286,8 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // insert points, update errors of affected triangles until threshold error is reached while (!heap.empty() && heap.top().error > threshold){ // get top element (with largest error) from heap - int maxelementindex = heap.top().index; - auto max_p = cpts.at(maxelementindex); + point_error maxelement = heap.top(); + auto max_p = maxelement.point; // get triangles that will change after inserting this max_p std::vector faces; @@ -305,7 +307,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // collect the heap_handles that need to be removed std::vector to_erase; for (auto it = face.info().points_inside->begin(); it != face.info().points_inside->end(); ++it) { - if (maxelementindex == (**it).index) { + if (maxelement.index == (**it).index) { to_erase.push_back(it); } } @@ -335,7 +337,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { } if (face->info().points_inside && face->info().points_inside->size() > 0) { for (auto h : *face->info().points_inside) { - if (maxelementindex != (*h).index) { + if (maxelement.index != (*h).index) { points_to_update.push_back(h); @@ -343,7 +345,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { //std::cout << "pushed: " << (*h).index << std::endl; //int index = (*h).index; - //if (index < 0 || index >= cpts.size()) { + //if (index < 0 || index >= pts.size()) { // std::cout << "looking for index: " << maxelementindex << std::endl; // std::cout << "index out of range: " << index << std::endl; //} @@ -369,7 +371,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // // Locate the face where // CDT::Locate_type lt; // int li; - // CDT::Face_handle containing_face = T.locate(cpts.at(maxelementindex), lt, li, faces[0]); + // CDT::Face_handle containing_face = T.locate(max_p, lt, li, faces[0]); // std::cout << "; point location: "; // if (lt == CDT::EDGE) { // std::cout << "on edge."; @@ -386,7 +388,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // else if (lt == CDT::OUTSIDE_AFFINE_HULL) { // std::cout << "outside affine hull."; // } - // std::cout << " Point; " << std::fixed << std::setprecision(3) << cpts.at(maxelementindex) << std::endl; + // std::cout << " Point; " << std::fixed << std::setprecision(3) << max_p << std::endl; // std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; //} @@ -399,21 +401,17 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { heap.pop(); // update the errors of affected elevation points - for (auto curelement : points_to_update){ - int index = (*curelement).index; - //std::cout << "pulled: " << (*curelement).index << std::endl; - //if (index > cpts.size()) { - // std::cout << "index out of range: " << index << std::endl; - //} - auto p = cpts.at(index); - //auto containing_face = T.locate(p, face_hint); + for (heap_handle curelement : points_to_update){ + point_error element = *curelement; + auto p = element.point; CDT::Locate_type lt; int li; CDT::Face_handle containing_face = T.locate(p, lt, li, face_hint); if (lt == CDT::EDGE || lt == CDT::FACE) { const double e = compute_error(p, containing_face); - const point_error new_pe = point_error(index, e); - heap.update(curelement, new_pe); + //const point_error new_pe = point_error(element.index, e, element.point); + element.error = e; + heap.update(curelement, element); containing_face->info().points_inside->push_back(curelement); } else { @@ -421,7 +419,8 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { if (lt == CDT::VERTEX) { std::cout << "on vertex."; // update error to 0.0 to disable point - heap.update(curelement, point_error(index, 0.0)); + element.error = 0.0; + heap.update(curelement, element); } else if (lt == CDT::OUTSIDE_CONVEX_HULL) { std::cout << "outside convex hull."; @@ -429,7 +428,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { else if (lt == CDT::OUTSIDE_AFFINE_HULL) { std::cout << "outside affine hull."; } - std::cout << " Point (index " << index << "); " << std::fixed << std::setprecision(3) << p << std::endl; + std::cout << " Point (index " << element.index << "); " << std::fixed << std::setprecision(3) << p << std::endl; } } } From 7c4b8acc2a43c5a52d5ef4ac5c339dd5fc0fa061 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 30 Apr 2019 12:19:32 +0200 Subject: [PATCH 49/53] update types --- src/geomtools.cpp | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index c5f775a9..3281199c 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -245,22 +245,15 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // assumes all lidar points are inside a triangle Heap heap; - //// Convert all elevation points to CGAL points - //std::vector cpts; - //cpts.reserve(pts.size()); - //for (auto& p : pts) { - // cpts.push_back(Point(bg::get<0>(p), bg::get<1>(p), bg::get<2>(p))); - //} - // compute initial point errors, build heap, store point indices in triangles { for (int i = 0; i < pts.size(); i++) { - auto p = Point(bg::get<0>(pts[i]), bg::get<1>(pts[i]), bg::get<2>(pts[i])); + Point p = Point(bg::get<0>(pts[i]), bg::get<1>(pts[i]), bg::get<2>(pts[i])); CDT::Locate_type lt; int li; CDT::Face_handle face = T.locate(p, lt, li); if (lt == CDT::FACE) { - auto e = compute_error(p, face); + double e = compute_error(p, face); auto handle = heap.push(point_error(i, e, p)); face->info().points_inside->push_back(handle); } @@ -325,7 +318,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { continue; } - // update clear info of triangles that just changed, collect points that were inside these triangles + // clear info of triangles that just changed, collect points that were inside these triangles heap_handle_vec points_to_update; //bool maxe = false; @@ -408,9 +401,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { int li; CDT::Face_handle containing_face = T.locate(p, lt, li, face_hint); if (lt == CDT::EDGE || lt == CDT::FACE) { - const double e = compute_error(p, containing_face); - //const point_error new_pe = point_error(element.index, e, element.point); - element.error = e; + element.error = compute_error(p, containing_face); heap.update(curelement, element); containing_face->info().points_inside->push_back(curelement); } @@ -418,7 +409,7 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { std::cout << "CDT update; point location not in face but "; if (lt == CDT::VERTEX) { std::cout << "on vertex."; - // update error to 0.0 to disable point + // update error to 0.0 to disable adding point to CDT element.error = 0.0; heap.update(curelement, element); } From e4eeb92032e11096c6a90b24018325dada012c68 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 30 Apr 2019 12:20:40 +0200 Subject: [PATCH 50/53] clean debugging code --- src/geomtools.cpp | 49 ----------------------------------------------- 1 file changed, 49 deletions(-) diff --git a/src/geomtools.cpp b/src/geomtools.cpp index 3281199c..03eadf2d 100644 --- a/src/geomtools.cpp +++ b/src/geomtools.cpp @@ -321,8 +321,6 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { // clear info of triangles that just changed, collect points that were inside these triangles heap_handle_vec points_to_update; - //bool maxe = false; - //int pointless = 0; for (auto face : faces) { if (face->info().plane) { delete face->info().plane; @@ -332,58 +330,11 @@ void greedy_insert(CDT &T, const std::vector &pts, double threshold) { for (auto h : *face->info().points_inside) { if (maxelement.index != (*h).index) { points_to_update.push_back(h); - - - //TESTING - //std::cout << "pushed: " << (*h).index << std::endl; - - //int index = (*h).index; - //if (index < 0 || index >= pts.size()) { - // std::cout << "looking for index: " << maxelementindex << std::endl; - // std::cout << "index out of range: " << index << std::endl; - //} } - //else { - // maxe = true; - //} } heap_handle_vec().swap((*face->info().points_inside)); } - //else { - // pointless++; - //} } - //std::cout << "Pointsless faces: " << pointless << std::endl; - //if (pointless == faces.size()) { - // std::cout << "Non of the located faces contain any points, index: " << maxelementindex << std::endl; - //} - - //if (!maxe) { - // std::cout << "maxelement " << maxelementindex << " not found"; - - // // Locate the face where - // CDT::Locate_type lt; - // int li; - // CDT::Face_handle containing_face = T.locate(max_p, lt, li, faces[0]); - // std::cout << "; point location: "; - // if (lt == CDT::EDGE) { - // std::cout << "on edge."; - // } - // else if (lt == CDT::FACE) { - // std::cout << "in face."; - // } - // else if (lt == CDT::VERTEX) { - // std::cout << "on vertex."; - // } - // else if (lt == CDT::OUTSIDE_CONVEX_HULL) { - // std::cout << "outside convex hull."; - // } - // else if (lt == CDT::OUTSIDE_AFFINE_HULL) { - // std::cout << "outside affine hull."; - // } - // std::cout << " Point; " << std::fixed << std::setprecision(3) << max_p << std::endl; - // std::cout << " Vertex; " << std::fixed << std::setprecision(3) << containing_face->vertex(li)->point() << std::endl; - //} // insert max_p in triangulation auto face_hint = faces[0]; From 3652c596e39433ddf09dbc388abca2344e03692b Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 30 Apr 2019 12:46:25 +0200 Subject: [PATCH 51/53] up version to v1.2.1 --- src/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.cpp b/src/main.cpp index a22a249e..503d8fb1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -36,7 +36,7 @@ #include #include -std::string VERSION = "1.2"; +std::string VERSION = "1.2.1"; bool validate_yaml(const char* arg, std::set& allowedFeatures); int main(int argc, const char * argv[]); From 4c23fe1b846f9cb67aafa8c22d64458c67d1cd1a Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 7 May 2019 15:46:39 +0200 Subject: [PATCH 52/53] polyfit3d; use x0 and y0 to normalize, not x[0] and y[0] --- src/TopoFeature.cpp | 6 ++++-- thirdparty/polyfit.hpp | 40 +++++++++++++++++++--------------------- 2 files changed, 23 insertions(+), 23 deletions(-) diff --git a/src/TopoFeature.cpp b/src/TopoFeature.cpp index d8d1c464..71e51a97 100644 --- a/src/TopoFeature.cpp +++ b/src/TopoFeature.cpp @@ -1185,6 +1185,8 @@ void Boundary3D::detect_outliers(bool flatten){ // find spikes in roads (due to misclassified lidar points) and fix by averaging between previous and next vertex. std::vector idx; std::vector x, y, z, coeffs; + double x0 = ring[0].x(); + double y0 = ring[0].y(); for (int i = 0; i < ring.size(); i++) { idx.push_back(i); x.push_back(ring[i].x()); @@ -1200,7 +1202,7 @@ void Boundary3D::detect_outliers(bool flatten){ for (int i = 0; i < niter; i++) { // Fit the model std::vector correctedvalues; - polyfit3d(xtmp, ytmp, ztmp, coeffs, correctedvalues); + polyfit3d(xtmp, ytmp, ztmp, x0, y0, coeffs, correctedvalues); std::vector residuals, absResiduals; double sum = 0; @@ -1249,7 +1251,7 @@ void Boundary3D::detect_outliers(bool flatten){ // get the new values based on the coeffs of the lase equation std::vector correctedvalues; - polyval3d(x, y, coeffs, correctedvalues); + polyval3d(x, y, x0, y0, coeffs, correctedvalues); if (flatten) { for (int i = 0; i < ring.size(); i++) { _p2z[ringi][i] = correctedvalues[i]; diff --git a/thirdparty/polyfit.hpp b/thirdparty/polyfit.hpp index b3b08709..e749d519 100644 --- a/thirdparty/polyfit.hpp +++ b/thirdparty/polyfit.hpp @@ -94,39 +94,37 @@ std::vector polyvalqr(const std::vector& oCoeff, const std::vector& oX) } template -void combineXY(std::vector& oX, std::vector& oY, mathalgo::matrix& oXY) { +void combineXY(std::vector& oX, std::vector& oY, T x0, T y0, mathalgo::matrix& oXY) { size_t nCount = oX.size(); size_t nCols = 3 * 2; // 3 unknowns in 2 dimensions oXY = mathalgo::matrix(nCount, nCols); - - // normalize x and y matrix - for (size_t i = 1; i < nCount; i++) { - oX[i] = oX[i] - oX[0]; - oY[i] = oY[i] - oY[0]; - } - oX[0] = 0; - oY[0] = 0; // create the XY matrix for (size_t nRow = 0; nRow < nCount; nRow++) { + T x = oX[nRow] - x0; + T y = oY[nRow] - y0; + if (nRow == 0) { + x = 0; + y = 0; + } oXY(nRow, 0) = 1; - oXY(nRow, 1) = oX[nRow]; - oXY(nRow, 2) = oY[nRow]; - oXY(nRow, 3) = oX[nRow] * oY[nRow]; - oXY(nRow, 4) = std::pow(oX[nRow], 2); - oXY(nRow, 5) = std::pow(oY[nRow], 2); + oXY(nRow, 1) = x; + oXY(nRow, 2) = y; + oXY(nRow, 3) = x * y; + oXY(nRow, 4) = x * x; + oXY(nRow, 5) = y * y; } } // 3D plane fitting template -void polyfit3d(std::vector& oX, std::vector& oY, std::vector& oZ, std::vector& coeffs, std::vector& calculated) { +void polyfit3d(std::vector& oX, std::vector& oY, std::vector& oZ, T x0, T y0, std::vector& coeffs, std::vector& calculated) { if (oX.size() != oY.size() || oX.size() != oZ.size()) throw std::invalid_argument("X and Y or X and Z vector sizes do not match"); size_t nCount = oX.size(); mathalgo::matrix A; - combineXY(oX, oY, A); + combineXY(oX, oY, x0, y0, A); mathalgo::matrix X(nCount, 1); // copy z matrix @@ -161,16 +159,16 @@ void polyfit3d(std::vector& oX, std::vector& oY, std::vector& oZ, std:: } template -void polyval3d(std::vector& x, std::vector& y, std::vector& coeff, std::vector& calculated) { +void polyval3d(std::vector& x, std::vector& y, T x0, T y0, std::vector& coeff, std::vector& calculated) { mathalgo::matrix A; - combineXY(x, y, A); - mathalgo::matrix Y(coeff.size(), 1); + combineXY(x, y, x0, y0, A); + mathalgo::matrix YT(coeff.size(), 1); // build coeff matrix for (size_t i = 0; i < coeff.size(); i++) { - Y(i, 0) = coeff[i]; + YT(i, 0) = coeff[i]; } mathalgo::matrix AY; - A.multiply(Y, AY); + A.multiply(YT, AY); calculated = AY.data(); } #endif \ No newline at end of file From ecf1d932cd64bbf22906491b709b99cb6eb349f8 Mon Sep 17 00:00:00 2001 From: Tom Commandeur Date: Tue, 7 May 2019 15:47:05 +0200 Subject: [PATCH 53/53] version 1.2.2 --- src/main.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.cpp b/src/main.cpp index 503d8fb1..7ebffe5e 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -36,7 +36,7 @@ #include #include -std::string VERSION = "1.2.1"; +std::string VERSION = "1.2.2"; bool validate_yaml(const char* arg, std::set& allowedFeatures); int main(int argc, const char * argv[]);