inverse distance gauss weighted nearest neighbour, work in progess

src/interpolation_method_nnn.c

@@ -758,6 +758,65 @@ generate_interpolation_plan_nnn(struct intermediate_method_data * data) {
 
    switch (reduction_type) {
 
+      case (NNN_GAUSS_WEIGHTED):
+
+         for (unsigned i = 0; i < num_interpolated_points; ++i) {
+
+        	unsigned j;
+
+            double scale_d = 1.0; // shall be user defined in xml
+        	double distance_sum = 0.0;
+
+        	for (j = 0; j < num_src_per_tgt; ++j) {
+               distances[j] = acos(tgt_results[i].results[j].cos_angle);
+               if (distances[j] <= yac_angle_tol) break;
+               distances[j] *= distances[j];
+               distance_sum += 1.0 / distances[j];
+            }
+
+            if (j == num_src_per_tgt) {
+
+            	// if no source point directly overlaps with the target point
+
+            	// a) compute average distance between source points
+
+            	double sum_d = 0.0;
+
+            	for (int j = 0; j < num_src_per_tgt; ++j) {
+                   for (int k = j + 1; k < num_src_per_tgt; ++k) {
+                      double d =
+                         get_vector_angle(tgt_results[i].results[j].coordinate_xyz,
+                                          tgt_results[i].results[k].coordinate_xyz);
+                         sum_d += d;
+                   }
+                }
+
+            	sum_d = sum_d * sum_d / (double)num_src_per_tgt;;
+
+            	// b) calculate weights from Gauss weighted inverse distances
+
+            	distance_sum = 1.0 / distance_sum;
+
+                double sum_weights = 0.0;
+
+                for (unsigned k = 0; k < num_src_per_tgt; ++k) {
+                   weights[i*num_src_per_tgt+k] = exp ( -0.5 * distances[k] / scale_d * sum_d );
+                   sum_weights += weights[i*num_src_per_tgt+k];
+                }
+                for (unsigned k = 0; k < num_src_per_tgt; ++k)
+                   weights[i*num_src_per_tgt+k] = weights[i*num_src_per_tgt+k] / sum_weights;
+
+            } else {
+
+              // special case: source and target at same location
+              for ( unsigned k = 0; k < num_src_per_tgt; ++k )
+                 weights[i*num_src_per_tgt+k] = 0.0;
+              weights[i*num_src_per_tgt+j] = 1.0;
+            }
+
+         }
+      break;
+
       case (NNN_DISTANCE_WEIGHTED): 
          for (unsigned i = 0; i < num_interpolated_points; ++i) {
 

src/interpolation_method_nnn.h

@@ -49,6 +49,7 @@
 enum yac_nnn_reduction_type {
   NNN_AVERAGE,           //!< average of n source points
   NNN_DISTANCE_WEIGHTED, //!< distance weighted average of n source points
+  NNN_GAUSS_WEIGHTED,    //!< distance with Gauss weights of n source points
   NNN_RBF,               //!< radial basis functions
 };