PostGIS  2.5.2dev-r@@SVN_REVISION@@

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

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

49 {
50  uint32_t i, iter;
51  double delta;
52  double sum_curr = 0, sum_next = 0;
53  int hit = LW_FALSE;
54  double *distances = lwalloc(npoints * sizeof(double));
55 
56  sum_curr = calc_weighted_distances_3d(curr, points, npoints, distances);
57 
58  for (iter = 0; iter < max_iter; iter++)
59  {
60  POINT3D next = { 0, 0, 0 };
61  double denom = 0;
62 
64  for (i = 0; i < npoints; i++)
65  {
66  /* we need to use lower epsilon than in FP_IS_ZERO in the loop for calculation to converge */
67  if (distances[i] > DBL_EPSILON)
68  {
69  next.x += points[i].x / distances[i];
70  next.y += points[i].y / distances[i];
71  next.z += points[i].z / distances[i];
72  denom += 1.0 / distances[i];
73  }
74  else
75  {
76  hit = LW_TRUE;
77  }
78  }
79 
80  if (denom < DBL_EPSILON)
81  {
82  /* No movement - Final point */
83  break;
84  }
85 
86  /* Calculate the new point */
87  next.x /= denom;
88  next.y /= denom;
89  next.z /= denom;
90 
91  /* If any of the intermediate points in the calculation is found in the
92  * set of input points, the standard Weiszfeld method gets stuck with a
93  * divide-by-zero.
94  *
95  * To get ourselves out of the hole, we follow an alternate procedure to
96  * get the next iteration, as described in:
97  *
98  * Vardi, Y. and Zhang, C. (2011) "A modified Weiszfeld algorithm for the
99  * Fermat-Weber location problem." Math. Program., Ser. A 90: 559-566.
100  * DOI 10.1007/s101070100222
101  *
102  * Available online at the time of this writing at
103  * http://www.stat.rutgers.edu/home/cunhui/papers/43.pdf
104  */
105  if (hit)
106  {
107  double dx = 0, dy = 0, dz = 0;
108  double d_sqr;
109  hit = LW_FALSE;
110 
111  for (i = 0; i < npoints; i++)
112  {
113  if (distances[i] > DBL_EPSILON)
114  {
115  dx += (points[i].x - curr->x) / distances[i];
116  dy += (points[i].y - curr->y) / distances[i];
117  dz += (points[i].z - curr->z) / distances[i];
118  }
119  }
120 
121  d_sqr = sqrt(dx*dx + dy*dy + dz*dz);
122  if (d_sqr > DBL_EPSILON)
123  {
124  double r_inv = FP_MAX(0, 1.0 / d_sqr);
125  next.x = (1.0 - r_inv)*next.x + r_inv*curr->x;
126  next.y = (1.0 - r_inv)*next.y + r_inv*curr->y;
127  next.z = (1.0 - r_inv)*next.z + r_inv*curr->z;
128  }
129  }
130 
131  /* Check movement with next point */
132  sum_next = calc_weighted_distances_3d(&next, points, npoints, distances);
133  delta = sum_curr - sum_next;
134  if (delta < tol)
135  {
136  break;
137  }
138  else
139  {
140  curr->x = next.x;
141  curr->y = next.y;
142  curr->z = next.z;
143  sum_curr = sum_next;
144  }
145  }
146 
147  lwfree(distances);
148  return iter;
149 }
double x
Definition: liblwgeom.h:354
static double calc_weighted_distances_3d(const POINT3D *curr, const POINT4D *points, uint32_t npoints, double *distances)
Definition: lwgeom_median.c:33
void lwfree(void *mem)
Definition: lwutil.c:244
double y
Definition: liblwgeom.h:342
double x
Definition: liblwgeom.h:342
double z
Definition: liblwgeom.h:342
unsigned int uint32_t
Definition: uthash.h:78
#define LW_FALSE
Definition: liblwgeom.h:76
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
double z
Definition: liblwgeom.h:354
void * lwalloc(size_t size)
Definition: lwutil.c:229
double y
Definition: liblwgeom.h:354
#define FP_MAX(A, B)
Here is the call graph for this function:
Here is the caller graph for this function: