CNG: Point cost functions are now based off of inverse square weighting

pydda/cost_functions/_cost_functions_jax.py

@@ -267,7 +267,7 @@ def calculate_smoothness_gradient(
     return y.flatten()
 
 
-def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
+def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3):
     """
     Calculates the cost function related to point observations. A mean square error cost
     function term is applied to points that are within the sphere of influence
@@ -305,14 +305,15 @@ def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
     """
     J = 0.0
     for the_point in point_list:
-        the_box = jnp.logical_and(
-            jnp.logical_and(
-                jnp.abs(x - the_point["x"]) < roi, jnp.abs(y - the_point["y"]) < roi
-            ),
-            jnp.abs(z - the_point["z"]) < roi,
+        dist = jnp.sqrt(
+            (x - the_point["x"]) ** 2
+            + (y - the_point["y"]) ** 2
+            + (z - the_point["z"]) ** 2
         )
-        the_box = jnp.where(the_box, 1.0, 0.0)
-        J += jnp.sum(((u - the_point["u"]) ** 2 + (v - the_point["v"]) ** 2) * the_box)
+        dist = jnp.maximum(dist, 1.0)
+        weight = 1 / dist**2
+        weight = weight / jnp.sum(weight)
+        J += jnp.sum(weight * ((u - the_point["u"]) ** 2 + (v - the_point["v"]) ** 2))
 
     return J * Cp
 
@@ -358,18 +359,16 @@ def calculate_point_gradient(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
     gradJ_w = jnp.zeros_like(u)
 
     for the_point in point_list:
-        the_box = jnp.where(
-            jnp.logical_and(
-                jnp.logical_and(
-                    np.abs(x - the_point["x"]) < roi, np.abs(y - the_point["y"]) < roi
-                ),
-                np.abs(z - the_point["z"]) < roi,
-            ),
-            1.0,
-            0.0,
+        dist = jnp.sqrt(
+            (x - the_point["x"]) ** 2
+            + (y - the_point["y"]) ** 2
+            + (z - the_point["z"]) ** 2
         )
-        gradJ_u += 2 * (u - the_point["u"]) * the_box
-        gradJ_v += 2 * (v - the_point["v"]) * the_box
+        dist = jnp.maximum(dist, 1.0)
+        weight = 1 / dist**2
+        weight = weight / jnp.sum(weight)
+        gradJ_u += 2 * (u - the_point["u"]) * weight
+        gradJ_v += 2 * (v - the_point["v"]) * weight
 
     gradJ = jnp.stack([gradJ_u, gradJ_v, gradJ_w], axis=0).flatten()
     return gradJ * Cp

pydda/cost_functions/_cost_functions_numpy.py

@@ -303,7 +303,7 @@ def calculate_smoothness_gradient(
     return y.flatten()
 
 
-def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
+def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3, power=2):
     """
     Calculates the cost function related to point observations. A mean square error cost
     function term is applied to points that are within the sphere of influence
@@ -339,18 +339,17 @@ def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
     for the_point in point_list:
         # Instead of worrying about whole domain, just find points in radius of influence
         # Since we know that the weight will be zero outside the sphere of influence anyways
-        the_box = np.where(
-            np.logical_and.reduce(
-                (
-                    np.abs(x - the_point["x"]) < roi,
-                    np.abs(y - the_point["y"]) < roi,
-                    np.abs(z - the_point["z"]) < roi,
-                )
-            )
-        )
-        J += np.sum(
-            ((u[the_box] - the_point["u"]) ** 2 + (v[the_box] - the_point["v"]) ** 2)
+
+        dist = np.sqrt(
+            (x - the_point["x"]) ** 2
+            + (y - the_point["y"]) ** 2
+            + (z - the_point["z"]) ** 2
         )
+        dist = np.maximum(dist, 1.0)
+        weight = 1 / dist**2
+        weight = weight / np.max(weight)
+
+        J += np.sum(weight * ((u - the_point["u"]) ** 2 + (v - the_point["v"]) ** 2))
 
     return J * Cp
 
@@ -392,17 +391,16 @@ def calculate_point_gradient(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
     gradJ_w = np.zeros_like(u)
 
     for the_point in point_list:
-        the_box = np.where(
-            np.logical_and.reduce(
-                (
-                    np.abs(x - the_point["x"]) < roi,
-                    np.abs(y - the_point["y"]) < roi,
-                    np.abs(z - the_point["z"]) < roi,
-                )
-            )
+        dist = np.sqrt(
+            (x - the_point["x"]) ** 2
+            + (y - the_point["y"]) ** 2
+            + (z - the_point["z"]) ** 2
         )
-        gradJ_u[the_box] += 2 * (u[the_box] - the_point["u"])
-        gradJ_v[the_box] += 2 * (v[the_box] - the_point["v"])
+        dist = np.maximum(dist, 1.0)
+        weight = 1 / dist**2
+        weight = weight / np.max(weight)
+        gradJ_u += 2 * weight * (u - the_point["u"])
+        gradJ_v += 2 * weight * (v - the_point["v"])
 
     gradJ = np.stack([gradJ_u, gradJ_v, gradJ_w], axis=0).flatten()
     return gradJ * Cp

pydda/cost_functions/_cost_functions_tensorflow.py

@@ -333,23 +333,18 @@ def calculate_point_cost(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
     for the_point in point_list:
         # Instead of worrying about whole domain, just find points in radius of influence
         # Since we know that the weight will be zero outside the sphere of influence anyways
-        xp = tf.ones_like(x) * the_point["x"]
-        yp = tf.ones_like(y) * the_point["y"]
-        zp = tf.ones_like(z) * the_point["z"]
         up = tf.ones_like(u) * the_point["u"]
         vp = tf.ones_like(v) * the_point["v"]
-
-        the_box = tf.where(
-            tf.math.logical_and(
-                tf.math.logical_and(
-                    tf.math.abs(x - xp) < roi, tf.math.abs(y - yp) < roi
-                ),
-                tf.math.abs(z - zp) < roi,
-            ),
-            1.0,
-            0.0,
+        dist = tf.math.sqrt(
+            (x - the_point["x"]) ** 2
+            + (y - the_point["y"]) ** 2
+            + (z - the_point["z"]) ** 2
         )
-        J.assign_add(tf.math.reduce_sum(((u - up) ** 2 + (v - vp) ** 2) * the_box))
+        dist = tf.math.maximum(dist, 1.0)
+        weight = 1 / dist**2
+        weight = weight / tf.reduce_max(weight)
+
+        J.assign_add(tf.math.reduce_sum(((u - up) ** 2 + (v - vp) ** 2) * weight))
 
     return J * Cp
 
@@ -394,24 +389,19 @@ def calculate_point_gradient(u, v, x, y, z, point_list, Cp=1e-3, roi=500.0):
     for the_point in point_list:
         # Instead of worrying about whole domain, just find points in radius of influence
         # Since we know that the weight will be zero outside the sphere of influence anyways
-        xp = tf.ones_like(x, dtype=tf.float32) * the_point["x"]
-        yp = tf.ones_like(y, dtype=tf.float32) * the_point["y"]
-        zp = tf.ones_like(z, dtype=tf.float32) * the_point["z"]
         up = tf.ones_like(u, dtype=tf.float32) * the_point["u"]
         vp = tf.ones_like(v, dtype=tf.float32) * the_point["v"]
 
-        the_box = tf.where(
-            tf.math.logical_and(
-                tf.math.logical_and(
-                    tf.math.abs(x - xp) < roi, tf.math.abs(y - yp) < roi
-                ),
-                tf.math.abs(z - zp) < roi,
-            ),
-            1.0,
-            0.0,
+        dist = tf.math.sqrt(
+            (x - the_point["x"]) ** 2
+            + (y - the_point["y"]) ** 2
+            + (z - the_point["z"]) ** 2
         )
-        gradJ_u.assign_add((2 * (u - up) * the_box))
-        gradJ_v.assign_add((2 * (v - vp) * the_box))
+        dist = tf.math.maximum(dist, 1.0)
+        weight = 1 / dist**2
+        weight = weight / tf.reduce_max(weight)
+        gradJ_u.assign_add((2 * (u - up) * weight))
+        gradJ_v.assign_add((2 * (v - vp) * weight))
     gradJ = tf.stack([gradJ_u, gradJ_v, gradJ_w], axis=0)
     gradJ = tf.reshape(gradJ, (3 * np.prod(u.shape),))
     return gradJ * Cp