Add logic to pass from baseline UVW to station UVW

include/UvwUtils.h

+// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+
+// This code was part of DP3 (https://git.astron.nl/RD/DP3.git)
+
+#include <Baseline.h>
+#include <algorithm>
+#include <numeric>
+#include <vector>
+#include <xtensor/xtensor.hpp>
+
+std::vector<int> NSetupSplitUvw(unsigned int n_antennas,
+                                const std::vector<int> &ant1,
+                                const std::vector<int> &ant2) {
+  // Get the indices of the baselines needed to split the baseline UVWs into
+  // station UVWs. They are in such an order that the UVW of a station is known
+  // before used in another baseline to derive the UVW of the other station.
+  // It can handle cases where baselines occur in disjoint station groups
+  // like 0-1, 0-2, 1-2 and 3-4, 4-5, 5-6.
+  // Note that the first station of a group gets UVW=0. All other station UVWs
+  // are relative to it using the baseline UVWs.
+  // Also note that nr of groups can be derived from the size of the returned
+  // vector (because it contains no entry for the first antenna in a group).
+  std::vector<int> uvw_bl;
+  uvw_bl.reserve(n_antennas);
+
+  std::vector<bool> known(n_antennas, false);
+  unsigned int nset = 0;
+  // Loop until all stations are set.
+  while (nset < n_antennas) {
+    // Disjoint groups might exist, so keep a vector containing related antennae
+    // which are members of the same group.
+    std::vector<unsigned int> members(1);
+    // Set first unset station as the reference station (which gets UVW=0).
+    for (unsigned int i = 0; i < n_antennas; ++i) {
+      if (!known[i]) {
+        members[0] = i;
+        known[i] = true;
+        ++nset;
+        break;
+      }
+    }
+    // Loop through all members in the group.
+    // Note that new members can be appended in this loop.
+    for (unsigned int j = 0; j < members.size(); ++j) {
+      int refst = members[j];
+      // Find all stations having a baseline with the reference station.
+      for (unsigned int i = 0; i < ant1.size(); ++i) {
+        int a1 = ant1[i];
+        int a2 = ant2[i];
+        // Only take baselines into account for which one station is known,
+        // so the other can be derived from it.
+        // The unknown station becomes member of the group.
+        if (known[a1] != known[a2]) {
+          if (a1 == refst) {
+            uvw_bl.push_back(i);
+            members.push_back(a2);
+            known[a2] = true;
+            ++nset;
+          } else if (a2 == refst) {
+            uvw_bl.push_back(-(i + 1));
+            members.push_back(a1);
+            known[a1] = true;
+            ++nset;
+          }
+        }
+      }
+    }
+  }
+  return uvw_bl;
+}
+
+std::vector<int>
+NSetupSplitUvw(unsigned int nant, const std::vector<int> &antennas1,
+               const std::vector<int> &antennas2,
+               const std::vector<std::array<double, 3>> &antenna_positions) {
+  // sort baselines by length
+  std::vector<double> baseline_length_squared(antennas1.size());
+  for (size_t i = 0; i < baseline_length_squared.size(); ++i) {
+    int ant1 = antennas1[i];
+    int ant2 = antennas2[i];
+    double u = antenna_positions[ant2][0] - antenna_positions[ant1][0];
+    double v = antenna_positions[ant2][1] - antenna_positions[ant1][1];
+    double w = antenna_positions[ant2][2] - antenna_positions[ant1][2];
+    baseline_length_squared[i] = u * u + v * v + w * w;
+  }
+
+  auto compare_by_length = [&baseline_length_squared](size_t bl1, size_t bl2) {
+    return baseline_length_squared[bl1] > baseline_length_squared[bl2];
+  };
+
+  std::vector<size_t> bl_idx_sorted(antennas1.size());
+  std::iota(bl_idx_sorted.begin(), bl_idx_sorted.end(), 0);
+  std::sort(bl_idx_sorted.begin(), bl_idx_sorted.end(), compare_by_length);
+
+  // Initialize data structure to keep track of which baselines are selected,
+  // to which group an antenna belongs and which antennas are in a group
+  std::vector<size_t> bl_selection;
+  std::vector<int> ant_to_group_id_map(nant, -1);
+  std::vector<std::vector<int>> groups;
+
+  for (size_t bl : bl_idx_sorted) {
+    int ant1 = antennas1[bl];
+    int ant2 = antennas2[bl];
+    if (ant_to_group_id_map[ant1] == ant_to_group_id_map[ant2]) {
+      if (ant_to_group_id_map[ant1] == -1) {
+        // both antennas are unkown
+        // add both to a new group
+        ant_to_group_id_map[ant1] = groups.size();
+        ant_to_group_id_map[ant2] = groups.size();
+        groups.push_back({ant1, ant2});
+      } else {
+        // both antennas are known, and in the same group
+        // this baseline adds no new information
+        continue;
+      }
+    } else {
+      if (ant_to_group_id_map[ant1] == -1) {
+        // ant1 is unknown, add it to the group of ant2
+        ant_to_group_id_map[ant1] = ant_to_group_id_map[ant2];
+        groups[ant_to_group_id_map[ant1]].push_back(ant1);
+      } else if (ant_to_group_id_map[ant2] == -1) {
+        // ant2 is unknown, add it to the group of ant1
+        ant_to_group_id_map[ant2] = ant_to_group_id_map[ant1];
+        groups[ant_to_group_id_map[ant2]].push_back(ant2);
+      } else {
+        // both antennas are known, but in different groups
+        // Add the group of ant2 to the group of ant1
+        for (int ant : groups[ant_to_group_id_map[ant2]]) {
+          ant_to_group_id_map[ant] = ant_to_group_id_map[ant1];
+        }
+        groups[ant_to_group_id_map[ant1]].insert(
+            groups[ant_to_group_id_map[ant1]].end(),
+            groups[ant_to_group_id_map[ant2]].begin(),
+            groups[ant_to_group_id_map[ant2]].end());
+        groups[ant_to_group_id_map[ant2]].clear();
+      }
+    }
+    // Add baseline to selection
+    bl_selection.push_back(bl);
+    // Exit early if selected baselines already form a single spanning tree
+    if (bl_selection.size() == size_t(nant - 1))
+      break;
+  }
+
+  // Select the entries from vectors antenna1 and antenna2 that correspons to
+  // selected baselines
+  std::vector<int> antennas1_bl_selection;
+  std::transform(bl_selection.begin(), bl_selection.end(),
+                 std::back_inserter(antennas1_bl_selection),
+                 [&antennas1](size_t bl) -> int { return antennas1[bl]; });
+  std::vector<int> antennas2_bl_selection;
+  std::transform(bl_selection.begin(), bl_selection.end(),
+                 std::back_inserter(antennas2_bl_selection),
+                 [&antennas2](size_t bl) -> int { return antennas2[bl]; });
+
+  // Compute the indices for the selected baselines
+  std::vector<int> bl_idx =
+      NSetupSplitUvw(nant, antennas1_bl_selection, antennas2_bl_selection);
+
+  // Map the indices in the selection to indices in the original baselines
+  std::transform(bl_idx.begin(), bl_idx.end(), bl_idx.begin(),
+                 [&bl_selection](int bl) -> int {
+                   return bl >= 0 ? bl_selection[bl]
+                                  : -bl_selection[-bl - 1] - 1;
+                 });
+  return bl_idx;
+}
+
+void NSplitUVW(const std::vector<int> &baseline_indices,
+               const std::vector<dp3::base::Baseline> &baselines,
+               const xt::xtensor<double, 2> &uvw_bl,
+               xt::xtensor<double, 2> &uvw_ant) {
+  uvw_ant.fill(0.0);
+  for (unsigned int i = 0; i < baseline_indices.size(); ++i) {
+    int index = baseline_indices[i];
+    if (index < 0) {
+      // Ant2 is known.
+      index = -index - 1;
+      const size_t first_baseline = baselines[index].first;
+      const size_t second_baseline = baselines[index].second;
+      for (int j = 0; j < 3; ++j) {
+        uvw_ant(first_baseline, j) =
+            uvw_ant(second_baseline, j) - uvw_bl(index, j);
+      }
+    } else {
+      // Ant1 is known.
+      const size_t first_baseline = baselines[index].first;
+      const size_t second_baseline = baselines[index].second;
+      for (int j = 0; j < 3; ++j) {
+        uvw_ant(second_baseline, j) =
+            uvw_ant(first_baseline, j) + uvw_bl(index, j);
+      }
+    }
+  }
+}

tests/test_uvwutils.cpp

+
+// Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+/// @author Sebastiaan van der Tol
+
+#define BOOST_TEST_MODULE UVU_UTILS_TEST
+#include <UvwUtils.h>
+#include <algorithm>
+#include <array>
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_SUITE(uvwutils)
+
+BOOST_AUTO_TEST_CASE(nsetupsplituvw) {
+  //
+  // Simple test for nsetupSplitUVW with antenna position info
+  //
+  // The antenna layout is as follows
+  //
+  // 1  3  5---6
+  // |  |  |  /
+  // |  |  | /
+  // |  |  |/
+  // 0  2--4
+  //
+  // One of the baselines in the 4-5-6 cycle is redundant
+  // Because baseline 5-6 is the shortest in the cycle, it should not be
+  // selected All other baselines are needed to make the spanning trees
+  //
+
+  std::vector<std::array<double, 3>> antenna_pos = {
+      {0, 0, 0}, {0, 4, 0}, {2, 0, 0}, {2, 4, 0},
+      {4, 0, 0}, {4, 4, 0}, {7, 4, 0},
+  };
+  std::vector<int> antenna1 = {0, 3, 2, 4, 5, 4};
+  std::vector<int> antenna2 = {1, 2, 4, 5, 6, 6};
+
+  size_t nant = antenna_pos.size();
+
+  std::vector<int> bl_indices =
+      NSetupSplitUvw(nant, antenna1, antenna2, antenna_pos);
+
+  // Baseline nr. 4 (antennas 5-6) should be excluded
+  std::vector<int> good_bl_indices = {0, 1, 2, 3, 5};
+
+  BOOST_CHECK_EQUAL(bl_indices.size(), good_bl_indices.size());
+
+  // All good_indices must be present in the bl_indices found by nsetupSplitUVW
+  for (int good_bl_idx : good_bl_indices) {
+    // Reversed baselines are represented by negative indices
+    // These are a matches too
+    auto indices_match = [good_bl_idx](int idx) -> bool {
+      return ((good_bl_idx == idx) || (good_bl_idx == -(idx + 1)));
+    };
+    BOOST_CHECK(std::find_if(bl_indices.begin(), bl_indices.end(),
+                             indices_match) != bl_indices.end());
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()