Compute outside-to-outside TLENs properly

Even if the reads overlap we need an outside-to-outside TLEN like the function claims to calculate.

src/alignment.cpp

@@ -1292,14 +1292,20 @@ pair<int32_t, int32_t> compute_template_lengths(const int64_t& pos1, const vecto
         int64_t here = pos;
         for (auto& item : cigar) {
             // Trace along the cigar
-            if (item.second == 'M') {
-                // Bases are matched. Count them in the bounds and execute the operation
-                low = min(low, here);
-                here += item.first;
-                high = max(high, here);
-            } else if (item.second == 'D') {
-                // Only other way to advance in the reference
-                here += item.first;
+            switch (item.second) {
+                case 'M':
+                case 'X':
+                case '=':
+                    // Bases are matched or mismatched. Count them in the bounds and execute the operation
+                    low = min(low, here);
+                    high = max(high, here + item.first);
+                    // Fall through.
+                case 'D':
+                case 'N':
+                    // Even for deletions/gaps, we advance on the reference.
+                    here += item.first;
+                    break;
+                // Inserts and various clips/paddings don't do anything.
             }
         }
 
@@ -1308,16 +1314,12 @@ pair<int32_t, int32_t> compute_template_lengths(const int64_t& pos1, const vecto
 
     auto bounds1 = find_bounds(pos1, cigar1);
     auto bounds2 = find_bounds(pos2, cigar2);
-
-    // Compute the separation
-    int32_t dist = 0;
-    if (bounds1.first < bounds2.second) {
-        // The reads are in order
-        dist = bounds2.second - bounds1.first;
-    } else if (bounds2.first < bounds1.second) {
-        // The reads are out of order so the other bounds apply
-        dist = bounds1.second - bounds2.first;
-    }
+
+    // We wanted the distance between the outermost points. So find them.
+    auto min_start = std::min(bounds1.first, bounds2.first);
+    auto max_end = std::max(bounds1.second, bounds2.second);
+    // And find the distance
+    int32_t dist = max_end - min_start;
 
     if (pos1 < pos2) {
         // Count read 1 as the overall "leftmost", so its value will be positive