iterate_4d()

static uint32_t iterate_4d ( POINT3D *  curr, const POINT4D *  points, const uint32_t  npoints, const uint32_t  max_iter, const double  tol  )

static

Calculate denom to get the next point

Definition at line 48 of file lwgeom_median.c.

{
        uint32_t i, iter;
        double delta;
        double sum_curr = 0, sum_next = 0;
        int hit = LW_FALSE;
        double *distances = lwalloc(npoints * sizeof(double));
 
        sum_curr = calc_weighted_distances_3d(curr, points, npoints, distances);
 
        for (iter = 0; iter < max_iter; iter++)
        {
                POINT3D next = { 0, 0, 0 };
                double denom = 0;
 
                for (i = 0; i < npoints; i++)
                {
                        /* we need to use lower epsilon than in FP_IS_ZERO in the loop for calculation to converge */
                        if (distances[i] > DBL_EPSILON)
                        {
                                next.x += points[i].x / distances[i];
                                next.y += points[i].y / distances[i];
                                next.z += points[i].z / distances[i];
                                denom += 1.0 / distances[i];
                        }
                        else
                        {
                                hit = LW_TRUE;
                        }
                }
 
                if (denom < DBL_EPSILON)
                {
                        /* No movement - Final point */
                        break;
                }
 
                /* Calculate the new point */
                next.x /= denom;
                next.y /= denom;
                next.z /= denom;
 
                /* If any of the intermediate points in the calculation is found in the
                * set of input points, the standard Weiszfeld method gets stuck with a
                * divide-by-zero.
                *
                * To get ourselves out of the hole, we follow an alternate procedure to
                * get the next iteration, as described in:
                *
                * Vardi, Y. and Zhang, C. (2011) "A modified Weiszfeld algorithm for the
                * Fermat-Weber location problem."  Math. Program., Ser. A 90: 559-566.
                * DOI 10.1007/s101070100222
                *
                * Available online at the time of this writing at
                * http://www.stat.rutgers.edu/home/cunhui/papers/43.pdf
                */
                if (hit)
                {
                        double dx = 0, dy = 0, dz = 0;
                        double d_sqr;
                        hit = LW_FALSE;
 
                        for (i = 0; i < npoints; i++)
                        {
                                if (distances[i] > DBL_EPSILON)
                                {
                                        dx += (points[i].x - curr->x) / distances[i];
                                        dy += (points[i].y - curr->y) / distances[i];
                                        dz += (points[i].z - curr->z) / distances[i];
                                }
                        }
 
                        d_sqr = sqrt(dx*dx + dy*dy + dz*dz);
                        if (d_sqr > DBL_EPSILON)
                        {
                                double r_inv = FP_MAX(0, 1.0 / d_sqr);
                                next.x = (1.0 - r_inv)*next.x + r_inv*curr->x;
                                next.y = (1.0 - r_inv)*next.y + r_inv*curr->y;
                                next.z = (1.0 - r_inv)*next.z + r_inv*curr->z;
                        }
                }
 
                /* Check movement with next point */
                sum_next = calc_weighted_distances_3d(&next, points, npoints, distances);
                delta = sum_curr - sum_next;
                if (delta < tol)
                {
                        break;
                }
                else
                {
                        curr->x = next.x;
                        curr->y = next.y;
                        curr->z = next.z;
                        sum_curr = sum_next;
                }
        }
 
        lwfree(distances);
        return iter;
}

References calc_weighted_distances_3d(), FP_MAX, LW_FALSE, LW_TRUE, lwalloc(), lwfree(), POINT3D::x, POINT4D::x, POINT3D::y, POINT4D::y, POINT3D::z, and POINT4D::z.

Referenced by lwmpoint_median().

Here is the call graph for this function:

Here is the caller graph for this function: