Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix bug in lwgeom_median and avoid division by zero #245

Closed
wants to merge 1 commit into from
Closed

Fix bug in lwgeom_median and avoid division by zero #245

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions NEWS
View file
Open in desktop
Expand Up @@ -5,6 +5,7 @@ PostGIS 2.3.8

- #4071, ST_ClusterKMeans crash on NULL/EMPTY fixed (Darafei Praliaskouski)
- #4074, Disable interrupt tests by default (Regina Obe)
- #3997, fix bug in lwgeom_median and avoid division by zero (Raúl Marín)

PostGIS 2.3.7
2018/04/06
Expand Down
19 changes: 11 additions & 8 deletions liblwgeom/lwgeom_median.c
View file
Open in desktop
Expand Up @@ -85,25 +85,28 @@ iterate_3d(POINT3D* curr, const POINT3D* points, uint32_t npoints, double* dista
double dx = 0;
double dy = 0;
double dz = 0;
double r_inv;
POINT3D alt;
double d_sqr;

for (i = 0; i < npoints; i++)
{
if (distances[i] > 0)
{
dx += (points[i].x - curr->x) / distances[i];
dy += (points[i].y - curr->y) / distances[i];
dz += (points[i].y - curr->z) / distances[i];
dz += (points[i].z - curr->z) / distances[i];
}
}

r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz );
d_sqr = sqrt (dx*dx + dy*dy + dz*dz);

alt.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
alt.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
alt.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
if (d_sqr > DBL_EPSILON)
{
double r_inv = 1.0 / d_sqr;

next = alt;
next.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
next.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
next.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
}
}

delta = distance3d_pt_pt(curr, &next);
Expand Down