PostGIS  2.3.8dev-r@@SVN_REVISION@@

◆ 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().

42 {
43  uint32_t i;
44  POINT3D next = { 0, 0, 0 };
45  double delta;
46  double denom = 0;
47  char hit = LW_FALSE;
48 
49  calc_distances_3d(curr, points, npoints, distances);
50 
51  for (i = 0; i < npoints; i++)
52  {
53  if (distances[i] == 0)
54  hit = LW_TRUE;
55  else
56  denom += 1.0 / distances[i];
57  }
58 
59  for (i = 0; i < npoints; i++)
60  {
61  if (distances[i] > 0)
62  {
63  next.x += (points[i].x / distances[i]) / denom;
64  next.y += (points[i].y / distances[i]) / denom;
65  next.z += (points[i].z / distances[i]) / denom;
66  }
67  }
68 
69  /* If any of the intermediate points in the calculation is found in the
70  * set of input points, the standard Weiszfeld method gets stuck with a
71  * divide-by-zero.
72  *
73  * To get ourselves out of the hole, we follow an alternate procedure to
74  * get the next iteration, as described in:
75  *
76  * Vardi, Y. and Zhang, C. (2011) "A modified Weiszfeld algorithm for the
77  * Fermat-Weber location problem." Math. Program., Ser. A 90: 559-566.
78  * DOI 10.1007/s101070100222
79  *
80  * Available online at the time of this writing at
81  * http://www.stat.rutgers.edu/home/cunhui/papers/43.pdf
82  */
83  if (hit)
84  {
85  double dx = 0;
86  double dy = 0;
87  double dz = 0;
88  double r_inv;
89  POINT3D alt;
90  for (i = 0; i < npoints; i++)
91  {
92  if (distances[i] > 0)
93  {
94  dx += (points[i].x - curr->x) / distances[i];
95  dy += (points[i].y - curr->y) / distances[i];
96  dz += (points[i].y - curr->z) / distances[i];
97  }
98  }
99 
100  r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz );
101 
102  alt.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
103  alt.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
104  alt.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
105 
106  next = alt;
107  }
108 
109  delta = distance3d_pt_pt(curr, &next);
110 
111  curr->x = next.x;
112  curr->y = next.y;
113  curr->z = next.z;
114 
115  return delta;
116 }
double y
Definition: liblwgeom.h:339
double x
Definition: liblwgeom.h:339
#define FP_MIN(A, B)
double z
Definition: liblwgeom.h:339
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:795
static void calc_distances_3d(const POINT3D *curr, const POINT3D *points, uint32_t npoints, double *distances)
Definition: lwgeom_median.c:31
#define FP_MAX(A, B)
Here is the call graph for this function:
Here is the caller graph for this function: