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.

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 d_sqr;
89  double r_inv;
90
91  for (i = 0; i < npoints; i++)
92  {
93  if (distances[i] > 0)
94  {
95  dx += (points[i].x - curr->x) / distances[i];
96  dy += (points[i].y - curr->y) / distances[i];
97  dz += (points[i].z - curr->z) / distances[i];
98  }
99  }
100
101  d_sqr = sqrt (dx*dx + dy*dy + dz*dz);
102  r_inv = d_sqr > DBL_EPSILON ? 1.0 / d_sqr : 1.0;
103
104  next.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
105  next.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
106  next.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
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: