iterate_3d()

static double iterate_3d ( POINT3D *  curr, const POINT3D *  points, uint32_t  npoints, double *  distances  )

static

Definition at line 41 of file lwgeom_median.c.

References calc_distances_3d(), distance3d_pt_pt(), FP_MAX, FP_MIN, LW_FALSE, LW_TRUE, POINT3D::x, POINT3D::y, and POINT3D::z.

Referenced by lwmpoint_median().

{
        uint32_t i;
        POINT3D next = { 0, 0, 0 };
        double delta;
        double denom = 0;
        char hit = LW_FALSE;

        calc_distances_3d(curr, points, npoints, distances);

        for (i = 0; i < npoints; i++)
        {
                if (distances[i] == 0)
                        hit = LW_TRUE;
                else
                        denom += 1.0 / distances[i];
        }

        for (i = 0; i < npoints; i++)
        {
                if (distances[i] > 0)
                {
                        next.x += (points[i].x / distances[i]) / denom;
                        next.y += (points[i].y / distances[i]) / denom;
                        next.z += (points[i].z / distances[i]) / 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;
                double dy = 0;
                double dz = 0;
                double d_sqr;
                double r_inv;

                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].z - curr->z) / distances[i];
                        }
                }

                d_sqr = sqrt (dx*dx + dy*dy + dz*dz);
                r_inv = d_sqr > DBL_EPSILON ? 1.0 / d_sqr : 1.0;

                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);

        curr->x = next.x;
        curr->y = next.y;
        curr->z = next.z;

        return delta;
}

Here is the call graph for this function:

Here is the caller graph for this function: