diff --git a/.gitignore b/.gitignore
index 24f784a7..44e21caf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,7 +11,7 @@
 *.out.*
 *.lisp
 *.org
-*.py
+# *.py
 *.json
 /.gdb_history
 
diff --git a/Makefile b/Makefile
index f20a4557..4f681081 100644
--- a/Makefile
+++ b/Makefile
@@ -1,15 +1,16 @@
-# Copyright (c) 2020 Mark Polyakov, Karen Haining (If you edit the file, add your name here!)
-# 
+# Copyright (c) 2020 Mark Polyakov, Karen Haining, Edward Zhang
+# (If you edit the file, add your name here!)
+#
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
 # in the Software without restriction, including without limitation the rights
 # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 # copies of the Software, and to permit persons to whom the Software is
 # furnished to do so, subject to the following conditions:
-# 
+#
 # The above copyright notice and this permission notice shall be included in all
 # copies or substantial portions of the Software.
-# 
+#
 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -85,6 +86,7 @@ test: $(BIN) $(BSC) $(TEST_BIN)
 	# bash ./test/scripts/pyramid-incorrect.sh
 	bash ./test/scripts/readme-examples-test.sh
 	bash ./test/scripts/random-crap.sh
+	# bash ./test/scripts/tetra.sh
 
 $(TEST_BIN): $(TEST_OBJS)
 	$(CXX) $(LDFLAGS) -o $(TEST_BIN) $(TEST_OBJS) $(LIBS)
diff --git a/documentation/database.man b/documentation/database.man
index 3fbaa11d..fb4941e0 100644
--- a/documentation/database.man
+++ b/documentation/database.man
@@ -80,6 +80,19 @@ Sets the maximum distance \fImax\fP (degrees). Defaults to 15 if option is not s
 \fB--kvector-distance-bins\fP \fInum-bins\fP
 Sets the number of distance bins in the kvector building method to \fInum-bins\fP.  Defaults to 10000 if option is not selected, which is pretty reasonable for most cases.
 
+.SH TETRA DATABASE OPTIONS
+
+The Tetra database is needed for the Tetra star identification algorithm. Its primary payload is
+a pattern catalog , essentially a large hash table that stores 4-star patterns
+
+.TP
+\fB--tetra\fP
+Generate a Tetra database
+
+.TP
+\fB--tetra-max-angle \fImax\fP
+Sets the maximum allowable angle between stars in a Tetra pattern to \fImax\fP (degrees). Defaults to 12 if option is not selected.
+
 .SH OTHER OPTIONS
 
 .TP
diff --git a/documentation/pipeline.man b/documentation/pipeline.man
index 9ff82016..583ca78a 100644
--- a/documentation/pipeline.man
+++ b/documentation/pipeline.man
@@ -121,7 +121,7 @@ Runs the dummy centroiding algorithm (random centroid algorithm) with \fInum-sta
 
 .TP
 \fB--star-id-algo\fP \fIalgo\fP
-Runs the \fIalgo\fP star identification algorithm. Current options are "dummy", "gv", and "pyramid". Defaults to "dummy" if option is not selected.
+Runs the \fIalgo\fP star identification algorithm. Current options are "dummy", "gv", "tetra", and "pyramid". Defaults to "dummy" if option is not selected.
 
 .TP
 \fB--angular-tolerance\fP [\fItolerance\fP] Sets the estimated angular centroiding error tolerance,
diff --git a/reg.py b/reg.py
new file mode 100644
index 00000000..10569fbe
--- /dev/null
+++ b/reg.py
@@ -0,0 +1,65 @@
+# Script to check performance on images in given folder and prevent regression testing
+# Go over all images in some folder and check performance, log it to stdout (for now)
+# TODO: figure out some storage plan
+
+# CLI:
+# test directory name, output log name
+# (optional) attitude estimator
+# Example usage: python reg.py samples/ log.txt -att_estimator quest
+
+import subprocess
+import argparse
+
+# import os
+import glob
+import datetime
+
+parser = argparse.ArgumentParser()
+parser.add_argument("test_dir", type=str)
+parser.add_argument("log", type=str, default="log.txt")
+parser.add_argument("-att_estimator", type=str, nargs="?", default="dqm")
+args = parser.parse_args()
+
+print(f"Testing images in {args.test_dir}")
+print(f"attitude estimator: {args.att_estimator}")
+print(f"Logging outputs to {args.log}")
+
+
+def get_diff(expected, actual):
+    """Get element-wise angle difference between expected and actual (what we got)"""
+    return [expected[i] - actual[i] for i in range(len(actual))]
+
+
+output_log = open(args.log, "a+")  # append to end of log, don't overwrite
+for img_name in glob.glob(args.test_dir + "**/*.png", recursive=True):
+    print(f"===================={img_name}====================")
+    cmd = (
+        f"./lost pipeline \
+      --png {img_name} \
+      --fov 17 \
+      --centroid-algo cog \
+      --centroid-filter-brightest 8 \
+      --star-id-algo tetra \
+      --database tetra-algo-12.dat \
+      --false-stars 0 \
+      --attitude-algo {args.att_estimator} \
+      --print-attitude attitude.txt \
+      --plot-output annotated-res.png",
+    )
+    subprocess.call(cmd, shell=True)
+    # log attitude.txt
+    # Log:
+    # which image was tested
+    # Output from attitude.txt
+    dt = datetime.datetime.now().isoformat()
+    output_log.write("----------------------------\n")
+    output_log.write(f"{img_name}-{dt}-{args.att_estimator}\n")
+    output_log.write("----------------------------\n")
+    with open("attitude.txt") as att_log:
+        targets = ["attitude_ra", "attitude_de", "attitude_roll"]
+        for line in att_log:
+            line = line.split(" ")
+            if line[0] in targets:
+                output_log.write(line[1])
+
+output_log.close()
diff --git a/src/attitude-estimators.cpp b/src/attitude-estimators.cpp
index 1e80f692..34b45a05 100644
--- a/src/attitude-estimators.cpp
+++ b/src/attitude-estimators.cpp
@@ -125,7 +125,7 @@ float QuestCharPoly(float x, float a, float b, float c, float d, float s) {retur
 float QuestCharPolyPrime(float x, float a, float b, float c) {return 4*pow(x,3) - 2*(a+b)*x - c;}
 
 /**
- * Approximates roots of a real function using the Newton-Raphson algorithm 
+ * Approximates roots of a real function using the Newton-Raphson algorithm
  * @see https://www.geeksforgeeks.org/program-for-newton-raphson-method/
  */
 float QuestEigenvalueEstimator(float guess, float a, float b, float c, float d, float s) {
diff --git a/src/attitude-utils.cpp b/src/attitude-utils.cpp
index fdc8f8c1..4757e466 100644
--- a/src/attitude-utils.cpp
+++ b/src/attitude-utils.cpp
@@ -169,6 +169,16 @@ float FloatModulo(float x, float mod) {
     return result >= 0 ? result : result + mod;
 }
 
+std::ostream &operator<<(std::ostream &output, const Vec2 &vec) {
+    output << "Vec2(x: " << vec.x << ", y: " << vec.y << ")";
+    return output;
+}
+
+std::ostream &operator<<(std::ostream &output, const Vec3 &vec) {
+    output << "Vec3(x: " << vec.x << ", y: " << vec.y << ", z: " << vec.z << ")";
+    return output;
+}
+
 /// The square of the magnitude
 float Vec3::MagnitudeSq() const {
     return fma(x,x,fma(y,y, z*z));
@@ -200,12 +210,12 @@ float Vec3::operator*(const Vec3 &other) const {
     return fma(x,other.x, fma(y,other.y, z*other.z));
 }
 
-/// Dot product
+/// Vector-scalar multiplication
 Vec2 Vec2::operator*(const float &other) const {
     return { x*other, y*other };
 }
 
-/// Vector-Scalar multiplication
+/// Vector-scalar multiplication
 Vec3 Vec3::operator*(const float &other) const {
     return { x*other, y*other, z*other };
 }
@@ -225,6 +235,10 @@ Vec3 Vec3::operator-(const Vec3 &other) const {
     return { x - other.x, y - other.y, z - other.z };
 }
 
+Vec3 Vec3::operator+(const Vec3 &other) const {
+    return {x + other.x, y + other.y, z + other.z};
+}
+
 /// Usual vector cross product
 Vec3 Vec3::CrossProduct(const Vec3 &other) const {
     return {
diff --git a/src/attitude-utils.hpp b/src/attitude-utils.hpp
index 842c1220..dd7ee358 100644
--- a/src/attitude-utils.hpp
+++ b/src/attitude-utils.hpp
@@ -1,7 +1,10 @@
 #ifndef ATTITUDE_UTILS_H
 #define ATTITUDE_UTILS_H
 
+#include <algorithm>
+#include <fstream>
 #include <memory>
+#include <numeric>  // iota
 #include <vector>
 
 #include "serialize-helpers.hpp"
@@ -26,6 +29,8 @@ struct Vec2 {
     Vec2 operator*(const float &) const;
     Vec2 operator-(const Vec2 &) const;
     Vec2 operator+(const Vec2 &) const;
+
+    friend std::ostream &operator<<(std::ostream &output, const Vec2 &vec);
 };
 
 class Mat3; // define above so we can use in Vec3 class
@@ -37,6 +42,9 @@ class Vec3 {
     float y;
     float z;
 
+    Vec3(){};
+    Vec3(float x, float y, float z) : x(x), y(y), z(z){};
+
     float Magnitude() const;
     float MagnitudeSq() const;
     Vec3 Normalize() const;
@@ -45,8 +53,11 @@ class Vec3 {
     Vec3 operator*(const float &) const;
     Vec3 operator*(const Mat3 &) const;
     Vec3 operator-(const Vec3 &) const;
+    Vec3 operator+(const Vec3 &) const;
     Vec3 CrossProduct(const Vec3 &) const;
     Mat3 OuterProduct(const Vec3 &) const;
+
+    friend std::ostream &operator<<(std::ostream &output, const Vec3 &vec);
 };
 
 /// 3x3 vector with floating point components
@@ -183,8 +194,24 @@ float ArcSecToRad(float);
 /// Always returns something in [0,mod) Eg -0.8 mod 0.6 = 0.4
 float FloatModulo(float x, float mod);
 
+// Argsort function - Tetra
+// Sort first vector based on values of second vector (asc)
+template <class T, class U>
+std::vector<T> ArgsortVector(std::vector<T> arr, std::vector<U> cmp) {
+    std::vector<T> res;
+    std::vector<int> indices(arr.size());
+    std::iota(indices.begin(), indices.end(), 0);
+    std::sort(indices.begin(), indices.end(),
+              [&](int a, int b) -> bool { return cmp[a] < cmp[b]; });
+
+    for (int ind : indices) {
+        res.push_back(arr[ind]);
+    }
+    return res;
+}
+
 // TODO: quaternion and euler angle conversion, conversion between ascension/declination to rec9tu
 
-}
+}  // namespace lost
 
 #endif
diff --git a/src/camera.hpp b/src/camera.hpp
index 2c076cfd..01634329 100644
--- a/src/camera.hpp
+++ b/src/camera.hpp
@@ -41,7 +41,7 @@ class Camera {
     /// Focal length in pixels
     float FocalLength() const { return focalLength; };
     /// Horizontal field of view in radians
-    float Fov() const;
+    float Fov() const; // in radians
 
     void SetFocalLength(float focalLength) { this->focalLength = focalLength; }
 
diff --git a/src/database-options.hpp b/src/database-options.hpp
index bbad9e8f..4224531f 100644
--- a/src/database-options.hpp
+++ b/src/database-options.hpp
@@ -11,4 +11,6 @@ LOST_CLI_OPTION("kvector-max-distance"   , float      , kvectorMaxDistance    ,
 LOST_CLI_OPTION("kvector-distance-bins"  , long       , kvectorNumDistanceBins, 10000 , atol(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("swap-integer-endianness", bool       , swapIntegerEndianness , false , atobool(optarg), true)
 LOST_CLI_OPTION("swap-float-endianness"  , bool       , swapFloatEndianness   , false , atobool(optarg), true)
+LOST_CLI_OPTION("tetra"                  , bool       , tetra                 , false , atobool(optarg), true)
+LOST_CLI_OPTION("tetra-max-angle"        , float      , tetraMaxAngle         , 12    , atof(optarg)   , kNoDefaultArgument)
 LOST_CLI_OPTION("output"                 , std::string, outputPath            , "-"   , optarg         , kNoDefaultArgument)
diff --git a/src/databases.cpp b/src/databases.cpp
index 93302d61..be2bb31d 100644
--- a/src/databases.cpp
+++ b/src/databases.cpp
@@ -1,12 +1,18 @@
 #include "databases.hpp"
 
-#include <stdlib.h>
 #include <assert.h>
 #include <math.h>
+#include <stdlib.h>
 
-#include <vector>
+// TODO: check which ones are useless
 #include <algorithm>
+#include <array>
 #include <iostream>
+#include <map>
+#include <set>
+#include <unordered_map>
+#include <utility>  // pair
+#include <vector>
 
 #include "attitude-utils.hpp"
 #include "serialize-helpers.hpp"
@@ -15,6 +21,7 @@
 namespace lost {
 
 const int32_t PairDistanceKVectorDatabase::kMagicValue = 0x2536f009;
+const int32_t TetraDatabase::kMagicValue = 0x26683787;
 
 struct KVectorPair {
     int16_t index1;
@@ -213,6 +220,165 @@ void SerializePairDistanceKVector(SerializeContext *ser, const Catalog &catalog,
     }
 }
 
+std::pair<std::vector<uint16_t>, std::vector<uint16_t>> TetraPreparePattCat(const Catalog &catalog,
+                                                                            const float maxFovDeg) {
+    // Would not recommend changing these parameters!
+    // Currently optimal to generate many patterns with smaller FOV
+    // If you make the maxFOV of Tetra database larger, you need to scale the following 2 parameters
+    // up note that this will cause number of patterns to grow exponentially
+    const int pattStarsPerFOV = 10;
+    const int verificationStarsPerFOV = 20;
+
+    // To eliminate double stars, specify that star must be > 0.05 degrees apart
+    const float starMinSep = 0.05;
+
+    const float maxFOV = DegToRad(maxFovDeg);
+
+    int numEntries = catalog.size();
+    int keepForPattCount = 1;
+    std::vector<bool> keepForPatterns(numEntries, false);
+    std::vector<bool> keepForVerifying(numEntries, false);
+
+    // NOTE: pattern set will always be a subset of verification set
+    // We technically don't even need the verification set unless we plan
+    // to calculate probability of mismatch - (which we probably should)
+
+    // Definitely keep the first star
+    keepForPatterns[0] = true;
+    keepForVerifying[0] = true;
+
+    for (int i = 1; i < numEntries; i++) {
+        // Spatial vector representing new star
+        Vec3 vec = catalog[i].spatial;
+
+        bool anglesForPattOK = true;
+        int numPattStarsInFov = 0;
+
+        // We should test each new star's angular distance to all stars
+        // we've already selected to be kept for pattern construction
+        // Stop early if:
+        // a) double star: angle < min separation allowed
+        // b) Number of stars in region maxFov/2 >= pattStarsPerFOV
+        for (int j = 0; j < i; j++) {
+            if (keepForPatterns[j]) {
+                float angle = Angle(vec, catalog[j].spatial);
+                if (angle < DegToRad(starMinSep)) {
+                    anglesForPattOK = false;
+                    break;
+                }
+                // If angle between new star i and old star j is less than maxFov/2, OK
+                if (angle < maxFOV / 2) {
+                    numPattStarsInFov++;
+                    if (numPattStarsInFov >= pattStarsPerFOV) {
+                        anglesForPattOK = false;
+                        break;
+                    }
+                }
+            }
+        }
+
+        if (anglesForPattOK) {
+            keepForPatterns[i] = true;
+            keepForVerifying[i] = true;
+            keepForPattCount++;
+            continue;
+        }
+
+        bool anglesForVerifOK = true;
+        int numVerStarsInFov = 0;
+
+        // Same thing here, we should test each new star's angular distance to all
+        // stars we've already selected to be kept for verification
+        for (int j = 0; j < i; j++) {
+            if (keepForVerifying[j]) {
+                float angle = Angle(vec, catalog[j].spatial);
+                if (angle < DegToRad(starMinSep)) {
+                    anglesForVerifOK = false;
+                    break;
+                }
+                if (angle < maxFOV / 2) {
+                    numVerStarsInFov++;
+                    if (numVerStarsInFov >= verificationStarsPerFOV) {
+                        anglesForVerifOK = false;
+                        break;
+                    }
+                }
+            }
+        }
+
+        if (anglesForVerifOK) {
+            keepForVerifying[i] = true;
+        }
+    }
+
+    // List of indices, track which stars in our star catalog are OK for Tetra to use
+    std::vector<uint16_t> finalCatIndices;
+    // List of indices, indexing into finalCatIndices
+    // Track which stars in the final star table should be used for pattern construction
+    // later in Tetra's database generation step
+    std::vector<uint16_t> pattStarIndices;
+
+    // finalCat is the final version of the star table
+    for (int i = 0; i < (int)keepForVerifying.size(); i++) {
+        if (keepForVerifying[i]) {
+            finalCatIndices.push_back(i);
+        }
+    }
+
+    // Find which stars in the final star table
+    // should be used for pattern construction later in Tetra's database generation step
+    // pattStarIndices will be a double-index:
+    // our "final star table" is a list of indices into the original (unchanged) catalog
+    // pattStarIndices is a list of indices into the final star table to tell Tetra which ones
+    // to use for pattern construction
+    int cumulativeSum = -1;
+    for (int i = 0; i < (int)keepForVerifying.size(); i++) {
+        if (keepForVerifying[i]) {
+            cumulativeSum++;
+        }
+        if (keepForPatterns[i]) {
+            pattStarIndices.push_back(cumulativeSum);
+        }
+    }
+    return std::pair<std::vector<uint16_t>, std::vector<uint16_t>>{finalCatIndices, pattStarIndices};
+}
+
+TetraDatabase::TetraDatabase(DeserializeContext *des) {
+    maxAngle_ = DeserializePrimitive<float>(des);
+    pattCatSize_ = DeserializePrimitive<uint64_t>(des);
+    tetraStarCatSize_ = DeserializePrimitive<uint64_t>(des);
+    pattCats_ = DeserializeArray<uint16_t>(des, 4 * pattCatSize_);
+    starCatInds_ = DeserializeArray<uint16_t>(des, tetraStarCatSize_);
+}
+
+TetraPatt TetraDatabase::GetPattern(int index) const {
+    TetraPatt res;
+    for (int i = 0; i < 4; i++) {
+        res.push_back(pattCats_[4*index + i]);
+    }
+    return res;
+}
+
+uint16_t TetraDatabase::GetTrueCatInd(int tetraInd) const {
+    return starCatInds_[tetraInd];
+}
+
+std::vector<TetraPatt> TetraDatabase::GetPatternMatches(int index) const {
+    std::vector<TetraPatt> res;
+    for (int c = 0;; c++) {
+        int i = (index + c * c) % PattCatSize();
+
+        TetraPatt tableRow = GetPattern(i);
+
+        if (tableRow[0] == 0 && tableRow[1] == 0) {
+            break;
+        } else {
+            res.push_back(tableRow);
+        }
+    }
+    return res;
+}
+
 /// Create the database from a serialized buffer.
 PairDistanceKVectorDatabase::PairDistanceKVectorDatabase(DeserializeContext *des)
     : index(KVectorIndex(des)) {
@@ -254,42 +420,217 @@ const int16_t *PairDistanceKVectorDatabase::FindPairsExact(const Catalog &catalo
     float minQueryCos = cos(maxQueryDistance);
 
     long liberalUpperIndex;
-    long liberalLowerIndex = index.QueryLiberal(minQueryDistance, maxQueryDistance, &liberalUpperIndex);
+    long liberalLowerIndex =
+        index.QueryLiberal(minQueryDistance, maxQueryDistance, &liberalUpperIndex);
     // now we need to find the first and last index that actually matches the query
     // step the lower index forward
-    // There's no good reason to be using >= and <= for the comparison against max/min, but the tests fail otherwise (because they use angle, with its acos, instead of forward cos like us). It's an insignificant difference.
-    while (liberalLowerIndex < liberalUpperIndex
-           && catalog[pairs[liberalLowerIndex*2]].spatial * catalog[pairs[liberalLowerIndex*2+1]].spatial >= maxQueryCos
-        )
-    { liberalLowerIndex++; }
+    // There's no good reason to be using >= and <= for the comparison against max/min, but the
+    // tests fail otherwise (because they use angle, with its acos, instead of forward cos like us).
+    // It's an insignificant difference.
+    while (liberalLowerIndex < liberalUpperIndex &&
+           catalog[pairs[liberalLowerIndex * 2]].spatial *
+                   catalog[pairs[liberalLowerIndex * 2 + 1]].spatial >=
+               maxQueryCos) {
+        liberalLowerIndex++;
+    }
 
     // step the upper index backward
     while (liberalLowerIndex < liberalUpperIndex
-           // the liberalUpperIndex is past the end of the logically returned range, so we need to subtract 1
-           && catalog[pairs[(liberalUpperIndex-1)*2]].spatial * catalog[pairs[(liberalUpperIndex-1)*2+1]].spatial <= minQueryCos
-        )
-    { liberalUpperIndex--; }
+           // the liberalUpperIndex is past the end of the logically returned range, so we need to
+           // subtract 1
+           && catalog[pairs[(liberalUpperIndex - 1) * 2]].spatial *
+                      catalog[pairs[(liberalUpperIndex - 1) * 2 + 1]].spatial <=
+                  minQueryCos) {
+        liberalUpperIndex--;
+    }
 
     *end = &pairs[liberalUpperIndex * 2];
     return &pairs[liberalLowerIndex * 2];
 }
 
 /// Number of star pairs stored in the database
-long PairDistanceKVectorDatabase::NumPairs() const {
-    return index.NumValues();
-}
+long PairDistanceKVectorDatabase::NumPairs() const { return index.NumValues(); }
 
-/// Return the distances from the given star to each star it's paired with in the database (for debugging).
-std::vector<float> PairDistanceKVectorDatabase::StarDistances(int16_t star, const Catalog &catalog) const {
+/// Return the distances from the given star to each star it's paired with in the database (for
+/// debugging).
+std::vector<float> PairDistanceKVectorDatabase::StarDistances(int16_t star,
+                                                              const Catalog &catalog) const {
     std::vector<float> result;
     for (int i = 0; i < NumPairs(); i++) {
-        if (pairs[i*2] == star || pairs[i*2+1] == star) {
-            result.push_back(AngleUnit(catalog[pairs[i*2]].spatial, catalog[pairs[i*2+1]].spatial));
+        if (pairs[i * 2] == star || pairs[i * 2 + 1] == star) {
+            result.push_back(
+                AngleUnit(catalog[pairs[i * 2]].spatial, catalog[pairs[i * 2 + 1]].spatial));
         }
     }
     return result;
 }
 
+void SerializeTetraDatabase(SerializeContext *ser, const Catalog &catalog, float maxFovDeg,
+                            const std::vector<uint16_t> &pattStarIndices,
+                            const std::vector<uint16_t> &catIndices) {
+    const float maxFovRad = DegToRad(maxFovDeg);
+
+    const int pattBins = TetraConstants::numPattBins;
+    const int pattSize = TetraConstants::numPattStars;
+    const int tempBins = 4;
+
+    Catalog tetraCatalog;
+    for (int ind : catIndices) {
+        tetraCatalog.push_back(catalog[ind]);
+    }
+
+    auto spatialHash = [](const Vec3 &vec) {
+        std::hash<float> hasher;
+        return hasher(vec.x) ^ hasher(vec.y) ^ hasher(vec.z);
+    };
+    auto spatialEquals = [](const Vec3 &vec1, const Vec3 &vec2) {
+        return vec1.x == vec2.x && vec1.y == vec2.y && vec1.z == vec2.z;
+    };
+    std::unordered_map<Vec3, std::vector<uint16_t>, decltype(spatialHash), decltype(spatialEquals)>
+        tempCoarseSkyMap(8, spatialHash, spatialEquals);
+
+    for (int starID : pattStarIndices) {
+        Vec3 v{tetraCatalog[starID].spatial};
+        Vec3 hash{floor((v.x + 1) * tempBins), floor((v.y + 1) * tempBins),
+                  floor((v.z + 1) * tempBins)};
+        tempCoarseSkyMap[hash].push_back(starID);
+    }
+
+    auto tempGetNearbyStars = [&tempCoarseSkyMap, &tetraCatalog](const Vec3 &vec, float radius) {
+        std::vector<float> components{vec.x, vec.y, vec.z};
+
+        std::vector<std::vector<int>> hcSpace;
+        for (float x : components) {
+            std::vector<int> range;
+            int lo = int((x + 1 - radius) * tempBins);
+            lo = std::max(lo, 0);
+            int hi = int((x + 1 + radius) * tempBins);
+            hi = std::min(hi + 1, 2 * tempBins);
+            range.push_back(lo);
+            range.push_back(hi);
+            hcSpace.push_back(range);
+        }
+        // Hashcode space has 3 ranges, one for each of [x, y, z]
+        std::vector<uint16_t> nearbyStarIDs;  // TODO: typedef this
+        for (int a = hcSpace[0][0]; a < hcSpace[0][1]; a++) {
+            for (int b = hcSpace[1][0]; b < hcSpace[1][1]; b++) {
+                for (int c = hcSpace[2][0]; c < hcSpace[2][1]; c++) {
+                    Vec3 code{static_cast<float>(a), static_cast<float>(b), static_cast<float>(c)};
+
+                    // For each star j in partition with key=code,
+                    // see if our star and j have angle < radius. If so, they are nearby
+                    for (uint16_t starID : tempCoarseSkyMap[code]) {
+                        float dotProd = vec * tetraCatalog[starID].spatial;
+                        if (dotProd > std::cos(radius)) {
+                            nearbyStarIDs.push_back(starID);
+                        }
+                    }
+                }
+            }
+        }
+        return nearbyStarIDs;
+    };
+    using Pattern = std::array<uint16_t, 4>;
+    std::vector<Pattern> pattList;
+    Pattern patt{0, 0, 0, 0};
+
+    // Construct all possible patterns
+    for (int firstStarInd : pattStarIndices) {
+        patt[0] = firstStarInd;
+
+        Vec3 v{tetraCatalog[firstStarInd].spatial};
+        Vec3 hashCode{floor((v.x + 1) * tempBins), floor((v.y + 1) * tempBins),
+                      floor((v.z + 1) * tempBins)};
+
+        // Remove star=firstStarInd from its sky map partition
+        auto removeIt = std::find(tempCoarseSkyMap[hashCode].begin(),
+                                  tempCoarseSkyMap[hashCode].end(), firstStarInd);
+        assert(removeIt != tempCoarseSkyMap[hashCode].end());
+        tempCoarseSkyMap[hashCode].erase(removeIt);
+
+        auto nearbyStars = tempGetNearbyStars(v, maxFovRad);
+        const int n = nearbyStars.size();
+        for (int i = 0; i < n; i++) {
+            for (int j = i + 1; j < n; j++) {
+                for (int k = j + 1; k < n; k++) {
+                    patt[1] = nearbyStars[i];
+                    patt[2] = nearbyStars[j];
+                    patt[3] = nearbyStars[k];
+
+                    bool pattFits = true;
+                    for (int pair1 = 0; pair1 < pattSize; pair1++) {
+                        for (int pair2 = pair1 + 1; pair2 < pattSize; pair2++) {
+                            float dotProd = tetraCatalog[patt[pair1]].spatial *
+                                            tetraCatalog[patt[pair2]].spatial;
+                            if (dotProd <= std::cos(maxFovRad)) {
+                                pattFits = false;
+                                break;
+                            }
+                        }
+                        if (!pattFits) break;
+                    }
+
+                    if (pattFits) {
+                        pattList.push_back(patt);
+                    }
+                }
+            }
+        }
+    }
+
+    std::cerr << "Tetra found " << pattList.size() << " patterns" << std::endl;
+    // Ensure load factor just < 0.5 to prevent against cycling
+    long long pattCatalogLen = 2 * (int)pattList.size() + 1;
+    std::vector<Pattern> pattCatalog(pattCatalogLen);
+    for (Pattern patt : pattList) {
+        std::vector<double> pattEdgeLengths;
+        for (int i = 0; i < pattSize; i++) {
+            CatalogStar star1 = tetraCatalog[patt[i]];
+            for (int j = i + 1; j < pattSize; j++) {
+                // calculate distance between vectors
+                CatalogStar star2 = tetraCatalog[patt[j]];
+                double edgeLen = (star2.spatial - star1.spatial).Magnitude();
+                pattEdgeLengths.push_back(edgeLen);
+            }
+        }
+        std::sort(pattEdgeLengths.begin(), pattEdgeLengths.end());
+        double pattLargestEdge = pattEdgeLengths.back();
+        std::vector<double> pattEdgeRatios;
+        // Skip last edge since ratio is just 1
+        for (int i = 0; i < pattEdgeLengths.size() - 1; i++) {
+            pattEdgeRatios.push_back(pattEdgeLengths[i] / pattLargestEdge);
+        }
+
+        std::vector<int> key;
+        for (double edgeRatio : pattEdgeRatios) {
+            key.push_back(int(edgeRatio * pattBins));
+        }
+
+        int hashIndex = KeyToIndex(key, pattBins, pattCatalogLen);
+        long long offset = 0;
+        // Quadratic probing to find next available bucket for element with key=hashIndex
+        while (true) {
+            int index = int(hashIndex + std::pow(offset, 2)) % pattCatalogLen;
+            offset++;
+            if (pattCatalog[index][0] == 0 && pattCatalog[index][1] == 0) {
+                pattCatalog[index] = patt;
+                break;
+            }
+        }
+    }
+    SerializePrimitive<float>(ser, maxFovDeg);
+    SerializePrimitive<uint64_t>(ser, pattCatalog.size());
+    SerializePrimitive<uint64_t>(ser, catIndices.size());
+    for (Pattern patt : pattCatalog) {
+        for (int i = 0; i < pattSize; i++) {
+            SerializePrimitive<uint16_t>(ser, patt[i]);
+        }
+    }
+    for (uint16_t ind : catIndices) {
+        SerializePrimitive<uint16_t>(ser, ind);
+    }
+}
+
 /**
    MultiDatabase memory layout:
 
diff --git a/src/databases.hpp b/src/databases.hpp
index c8337c11..375d31c1 100644
--- a/src/databases.hpp
+++ b/src/databases.hpp
@@ -1,8 +1,10 @@
 #ifndef DATABASE_BUILDER_H
 #define DATABASE_BUILDER_H
 
-#include <stdlib.h>
 #include <inttypes.h>
+#include <stdlib.h>
+
+#include <utility>
 #include <vector>
 
 #include "star-utils.hpp"
@@ -14,7 +16,7 @@ const int32_t kCatalogMagicValue = 0xF9A283BC;
 
 /**
  * A data structure enabling constant-time range queries into fixed numerical data.
- * 
+ *
  * @note Not an instantiable database on its own -- used in other databases
  */
 // TODO: QueryConservative, QueryExact, QueryTrapezoidal?
@@ -74,6 +76,72 @@ class PairDistanceKVectorDatabase {
     const int16_t *pairs;
 };
 
+/*
+Pre-processing for Tetra star-id algorithm
+Return:
+(a) List of stars we can use for Tetra
+(b) Subset of (a) that we use to generate Tetra star patterns
+*/
+std::pair<std::vector<uint16_t>, std::vector<uint16_t>> TetraPreparePattCat(const Catalog &,
+                                                                            const float maxFovDeg);
+
+void SerializeTetraDatabase(SerializeContext *, const Catalog &, float maxFovDeg,
+                            const std::vector<uint16_t> &pattStarIndices,
+                            const std::vector<uint16_t> &catIndices);
+
+/// Tetra star pattern = vector of 4 star IDs
+using TetraPatt = std::vector<uint16_t>;
+
+/**
+ * A database storing Tetra star patterns
+ * TODO: implement something to prevent cycling in quadratic probe
+ * (or guarantee load factor < 0.5)
+ *
+ * Layout:
+ * | size (bytes)                     | name         | description                                                 |
+ * |----------------------------------+--------------+-------------------------------------------------------------|                                 |
+ * | sizeof float                     | maxFov       | max angle (degrees) allowed between any 2 stars             |
+ * |                                  |              | in the same pattern                                         |
+ * | sizeof(uint64_t)                 | pattCatSize  | number of rows in pattern catalog                           |
+ * | sizeof(uint64_t)                 | tetraCatSize | number of Tetra catalog indices                             |
+ * | 4*pattCatSize * sizeof(uint16_t) | pattCat      | hash table for Tetra star patternss                         |
+ * | tetraCatSize * sizeof(uint16_t)  | tetraStarCat | list of catalog indices to use for Tetra star-id algo       |
+ * |----------------------------------+--------------+-------------------------------------------------------------|
+ */
+class TetraDatabase {
+   public:
+    explicit TetraDatabase(DeserializeContext *des);
+
+    /// Max angle (in degrees) allowed between stars in the same pattern
+    float MaxAngle() const {return maxAngle_;}
+
+    /// Number of rows in pattern catalog
+    // With load factor of just under 0.5, size = numPatterns*2 + 1
+    uint64_t PattCatSize() const {return pattCatSize_;}
+
+    /// Get the 4-tuple pattern at row=index, 0-based
+    TetraPatt GetPattern(int index) const;
+
+    /// Returns a list of patterns from the pattern catalog
+    // Starts from index and does quadratic probing
+    // We assume that collisions mean the pattern stored there COULD be a match
+    std::vector<TetraPatt> GetPatternMatches(int index) const;
+
+    uint16_t GetTrueCatInd(int tetraIndex) const;
+    // TODO: should probably have a field describing number of indices for future updates to db
+
+    /// Magic value to use when storing inside a MultiDatabase
+    static const int32_t kMagicValue;
+    // static const int headerSize = sizeof(float) + sizeof(uint64_t);
+
+   private:
+    float maxAngle_;
+    uint64_t pattCatSize_;
+    uint16_t tetraStarCatSize_;
+    const uint16_t* pattCats_;
+    const uint16_t* starCatInds_;
+};
+
 // /**
 //  * @brief Stores "inner angles" between star triples
 //  * @details Unsensitive to first-order error in basic camera
@@ -114,7 +182,7 @@ class MultiDatabaseEntry {
     std::vector<unsigned char> bytes;
 };
 
-typedef std::vector<MultiDatabaseEntry> MultiDatabaseDescriptor;
+using MultiDatabaseDescriptor = std::vector<MultiDatabaseEntry>;
 
 void SerializeMultiDatabase(SerializeContext *, const MultiDatabaseDescriptor &dbs);
 
diff --git a/src/io.cpp b/src/io.cpp
index 32ef63a9..90305cc4 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -1,23 +1,23 @@
 #include "io.hpp"
 
+#include <assert.h>
 #include <cairo/cairo.h>
-#include <stdio.h>
+#include <errno.h>
 #include <inttypes.h>
+#include <limits.h>
 #include <math.h>
-#include <errno.h>
-#include <assert.h>
+#include <stdio.h>
 #include <stdlib.h>
-#include <limits.h>
 #include <unistd.h>
 
-#include <vector>
-#include <string>
+#include <algorithm>
+#include <cstring>
 #include <fstream>
 #include <iostream>
 #include <memory>
-#include <cstring>
 #include <random>
-#include <algorithm>
+#include <string>
+#include <vector>
 #include <map>
 #include <chrono>
 
@@ -101,7 +101,7 @@ const Catalog &CatalogRead() {
         std::sort(catalog.begin(), catalog.end(), [](const CatalogStar &a, const CatalogStar &b) {
             return a.spatial.x < b.spatial.x;
         });
-        for (int i = catalog.size(); i > 0; i--) {
+        for (int i = catalog.size()-1; i > 0; i--) {
             if ((catalog[i].spatial - catalog[i-1].spatial).Magnitude() < 5e-5) { // 70 stars removed at this threshold.
                 if (catalog[i].magnitude > catalog[i-1].magnitude) {
                     catalog.erase(catalog.begin() + i);
@@ -269,6 +269,28 @@ SerializeContext serFromDbValues(const DatabaseOptions &values) {
     return SerializeContext(values.swapIntegerEndianness, values.swapFloatEndianness);
 }
 
+// void BuildTetraDatabase(MultiDatabaseBuilder *builder, const Catalog &catalog, float maxAngle,
+//                         const std::vector<uint16_t> &pattStars,
+//                         const std::vector<uint16_t> &catIndices) {
+//     // long length = SerializeTetraDatabase(catalog, maxAngle, nullptr, pattStars, catIndices, false);
+//     long length = SerializeLengthTetraDatabase(catalog, maxAngle, pattStars, catIndices);
+//     unsigned char *buffer = builder->AddSubDatabase(TetraDatabase::kMagicValue, length);
+//     if (buffer == nullptr) {
+//         std::cerr << "Error: No room for Tetra database" << std::endl;
+//     }
+//     // SerializeTetraDatabase(catalog, maxAngle, buffer, pattStars, catIndices, true);
+//     SerializeTetraDatabase(catalog, maxAngle, buffer, pattStars, catIndices);
+// }
+
+// void GenerateTetraDatabases(MultiDatabaseBuilder *builder, const Catalog &catalog,
+//                             const DatabaseOptions &values, const std::vector<uint16_t> &pattStars,
+//                             const std::vector<uint16_t> &catIndices) {
+//     float maxAngle = values.tetraMaxAngle;
+//     BuildTetraDatabase(builder, catalog, maxAngle, pattStars, catIndices);
+// }
+// /// Generate and add databases to the given multidatabase builder according to the command line options in `values`
+// void GenerateDatabases(MultiDatabaseBuilder *builder, const Catalog &catalog, const DatabaseOptions &values) {
+
 MultiDatabaseDescriptor GenerateDatabases(const Catalog &catalog, const DatabaseOptions &values) {
     MultiDatabaseDescriptor dbEntries;
 
@@ -277,14 +299,36 @@ MultiDatabaseDescriptor GenerateDatabases(const Catalog &catalog, const Database
     SerializeCatalog(&catalogSer, catalog, false, true);
     dbEntries.emplace_back(kCatalogMagicValue, catalogSer.buffer);
 
+    bool dbProvided = false;
+
     if (values.kvector) {
+        dbProvided = true;
         float minDistance = DegToRad(values.kvectorMinDistance);
         float maxDistance = DegToRad(values.kvectorMaxDistance);
         long numBins = values.kvectorNumDistanceBins;
         SerializeContext ser = serFromDbValues(values);
         SerializePairDistanceKVector(&ser, catalog, minDistance, maxDistance, numBins);
         dbEntries.emplace_back(PairDistanceKVectorDatabase::kMagicValue, ser.buffer);
-    } else {
+    }
+
+    if (values.tetra) {
+        dbProvided = true;
+        float maxAngleDeg = values.tetraMaxAngle;
+        std::cerr << "Tetra max angle: " << maxAngleDeg << std::endl;
+
+        auto tetraPrepRes = TetraPreparePattCat(catalog, maxAngleDeg);
+        std::vector<uint16_t> catIndices = tetraPrepRes.first;
+        std::vector<uint16_t> pattStarsInds = tetraPrepRes.second;
+
+        std::cerr << "Tetra processed catalog has " << catIndices.size() << " stars." << std::endl;
+        std::cerr << "Number of pattern stars: " << pattStarsInds.size() << std::endl;
+
+        SerializeContext ser = serFromDbValues(values);
+        SerializeTetraDatabase(&ser, catalog, maxAngleDeg, pattStarsInds, catIndices);
+        dbEntries.emplace_back(TetraDatabase::kMagicValue, ser.buffer);
+    }
+
+    if (!dbProvided) {
         std::cerr << "No database builder selected -- no database generated." << std::endl;
         exit(1);
     }
@@ -872,7 +916,11 @@ Pipeline SetPipeline(const PipelineOptions &values) {
     } else if (values.idAlgo == "gv") {
         result.starIdAlgorithm = std::unique_ptr<StarIdAlgorithm>(new GeometricVotingStarIdAlgorithm(DegToRad(values.angularTolerance)));
     } else if (values.idAlgo == "py") {
-        result.starIdAlgorithm = std::unique_ptr<StarIdAlgorithm>(new PyramidStarIdAlgorithm(DegToRad(values.angularTolerance), values.estimatedNumFalseStars, values.maxMismatchProb, 1000));
+        result.starIdAlgorithm = std::unique_ptr<StarIdAlgorithm>(new PyramidStarIdAlgorithm(
+            DegToRad(values.angularTolerance), values.estimatedNumFalseStars,
+            values.maxMismatchProb, 1000));
+    } else if (values.idAlgo == "tetra") {
+        result.starIdAlgorithm = std::unique_ptr<StarIdAlgorithm>(new TetraStarIdAlgorithm());
     } else if (values.idAlgo != "") {
         std::cout << "Illegal id algorithm." << std::endl;
         exit(1);
@@ -903,6 +951,8 @@ PipelineOutput Pipeline::Go(const PipelineInput &input) {
     // (human centipede) until there are no more stages set.
     PipelineOutput result;
 
+    // std::cerr << *input.InputCamera() << std::endl;
+
     const Image *inputImage = input.InputImage();
     const Stars *inputStars = input.InputStars();
     const StarIdentifiers *inputStarIds = input.InputStarIds();
diff --git a/src/io.hpp b/src/io.hpp
index 0ce9dff6..c78454d3 100644
--- a/src/io.hpp
+++ b/src/io.hpp
@@ -5,27 +5,26 @@
 
 #include <cairo/cairo.h>
 
-#include <random>
-#include <vector>
-#include <map>
-#include <utility>
-#include <string>
-#include <sstream>
 #include <iostream>
+#include <map>
 #include <memory>
-
+#include <random>
+#include <sstream>
+#include <string>
+#include <utility>
+#include <vector>
 
 #ifndef CAIRO_HAS_PNG_FUNCTIONS
 #error LOST requires Cairo to be compiled with PNG support
 #endif
 
-#include "centroiders.hpp"
-#include "star-utils.hpp"
-#include "star-id.hpp"
-#include "camera.hpp"
-#include "attitude-utils.hpp"
 #include "attitude-estimators.hpp"
+#include "attitude-utils.hpp"
+#include "camera.hpp"
+#include "centroiders.hpp"
 #include "databases.hpp"
+#include "star-id.hpp"
+#include "star-utils.hpp"
 
 namespace lost {
 
@@ -296,6 +295,11 @@ SerializeContext serFromDbValues(const DatabaseOptions &values);
 /// @sa SerializeMultiDatabase
 MultiDatabaseDescriptor GenerateDatabases(const Catalog &, const DatabaseOptions &values);
 
+// // TODO:  can we avoid the split?
+// void GenerateTetraDatabases(MultiDatabaseBuilder *, const Catalog &, const DatabaseOptions &values,
+//                             const std::vector<uint16_t> &pattStars,
+//                             const std::vector<uint16_t> &catIndices);
+
 /////////////////////
 // INSPECT CATALOG //
 /////////////////////
diff --git a/src/main.cpp b/src/main.cpp
index 5d26fec4..63eb1337 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -5,21 +5,22 @@
  */
 
 #include <assert.h>
-#include <unistd.h>
 #include <getopt.h>
+#include <unistd.h>
 
-#include <string>
-#include <iostream>
-#include <fstream>
 #include <chrono>
 #include <cstring>
+#include <fstream>
+#include <iostream>
 #include <map>
+#include <string>
 
-#include "databases.hpp"
 #include "centroiders.hpp"
+#include "databases.hpp"
 #include "io.hpp"
 #include "man-database.h"
 #include "man-pipeline.h"
+#include "star-utils.hpp"
 
 namespace lost {
 
diff --git a/src/star-id.cpp b/src/star-id.cpp
index b0787313..336211e5 100644
--- a/src/star-id.cpp
+++ b/src/star-id.cpp
@@ -1,21 +1,329 @@
-#include <stdlib.h>
-#include <math.h>
+#include "star-id.hpp"
+
 #include <assert.h>
-#include <vector>
+#include <math.h>
+#include <stdlib.h>
+
 #include <algorithm>
 #include <chrono>
+#include <set>
 #include <unordered_map>
+#include <utility>  // std::pair
 
-#include "star-id.hpp"
-#include "star-id-private.hpp"
-#include "databases.hpp"
 #include "attitude-utils.hpp"
+#include "databases.hpp"
+#include "star-id-private.hpp"
+#include "star-utils.hpp"
 
 namespace lost {
 
-StarIdentifiers DummyStarIdAlgorithm::Go(
-    const unsigned char *, const Stars &stars, const Catalog &catalog, const Camera &) const {
+/**
+ * @brief Tetra only: construct a Tetra star pattern from spatial vectors
+ *
+ * This is one possible way to represent the pattern, somewhat naive but good enough in practice
+ * For each pair of spatial vectors, calculate the magnitude of their difference, the "edge length"
+ * Calculate the largest such difference, i.e. the "largest edge"
+ * Divide all edge lengths by the largest edge (except the LE itself) to get edge ratios
+ * Return the pattern = vector of edge ratios in sorted order
+ */
+static std::vector<double> TetraConstructPattern(const std::vector<Vec3> &spats);
+
+static std::vector<double> TetraConstructPattern(const std::vector<Vec3> &spats) {
+    std::vector<double> edgeLengths;
+    // C(numPattStars, 2) edge lengths calculated
+    for (int i = 0; i < (int)spats.size(); i++) {
+        for (int j = i + 1; j < (int)spats.size(); j++) {
+            Vec3 diff = spats[i] - spats[j];
+            edgeLengths.push_back(diff.Magnitude());
+        }
+    }
+
+    std::sort(edgeLengths.begin(), edgeLengths.end());
+    double largestEdge = edgeLengths.back();  // largest edge value
+
+    // size==C(numPattStars, 2) - 1
+    // We don't calculate ratio for largest edge, which would just be 1
+    std::vector<double> edgeRatios;  // size() = C(numPattStars, 2) - 1
+    for (int i = 0; i < edgeLengths.size() - 1; i++) {
+        edgeRatios.push_back(edgeLengths[i] / largestEdge);
+    }
+    return edgeRatios;
+}
+
+class TetraCentroidComboIterator {
+   public:
+    TetraCentroidComboIterator(int numPattStars, int numCentroids)
+        : firstTime_(true), pattSize_(numPattStars), numCentroids_(numCentroids) {}
+
+    /**
+     * @brief Tetra only: get indices of centroids for new star pattern, store in res
+     *
+     * Assumes that your centroids have already been sorted in decreasing order of brightness
+     * Bright stars tend to have lower centroiding error
+     */
+    bool getCentroidCombo(std::vector<int> *const res) {
+        if (numCentroids_ < pattSize_) {
+            return false;
+        }
+
+        if (firstTime_) {
+            firstTime_ = false;
+            indices_.push_back(-1);
+            for (int i = 0; i < pattSize_; i++) {
+                indices_.push_back(i);
+            }
+            indices_.push_back(numCentroids_);
+            *res = std::vector<int>(indices_.begin() + 1, indices_.end() - 1);
+            return true;
+        }
+
+        if (indices_[1] >= numCentroids_ - pattSize_) {
+            return false;
+        }
+
+        for (int i = 1; i <= pattSize_; i++) {
+            indices_[i]++;
+            if (indices_[i] < indices_[i + 1]) {
+                break;
+            }
+            indices_[i] = indices_[i - 1] + 1;
+        }
+        *res = std::vector<int>(indices_.begin() + 1, indices_.end() - 1);
+        return true;
+    }
+
+private:
+ bool firstTime_;
+ std::vector<int> indices_;
+ int pattSize_;
+ int numCentroids_;
+};
+
+StarIdentifiers TetraStarIdAlgorithm::Go(const unsigned char *database, const Stars &stars,
+                                         const Catalog &catalog, const Camera &camera) const {
+    // Result format: (centroidIndex, catalogIndex)
+    StarIdentifiers result;
+
+    MultiDatabase multiDatabase(database);
+    const unsigned char *databaseBuffer =
+        multiDatabase.SubDatabasePointer(TetraDatabase::kMagicValue);
+    if (databaseBuffer == nullptr) {
+        std::cerr << "Could not get to Tetra database" << std::endl;
+        return result;
+    }
+    DeserializeContext des(databaseBuffer);
+    TetraDatabase tetraDatabase(&des);
+
+    const long long catLength = tetraDatabase.PattCatSize();
+    const float maxFov = tetraDatabase.MaxAngle();
+
+    std::vector<int> centroidIndices;
+    for (int i = 0; i < stars.size(); i++) {
+        centroidIndices.push_back(i);
+    }
+
+    // Sort centroided stars by brightness, high to low. Larger is brighter
+    std::sort(centroidIndices.begin(), centroidIndices.end(),
+              [&stars](int a, int b) { return stars[a].magnitude > stars[b].magnitude; });
+
+    // Index of centroid indices list
+    std::vector<int> chosenCentroidIndices(TetraConstants::numPattStars);
+    TetraCentroidComboIterator tetraCentroidComboIt(TetraConstants::numPattStars, centroidIndices.size());
+
+    // TODO: In practice, maybe cap this at some number of combinations, maybe 10 or so
+    while (tetraCentroidComboIt.getCentroidCombo(&chosenCentroidIndices)) {
+        // Index in centroid indices list
+        std::vector<int> chosenIndIndices;
+        std::vector<Star> chosenStars;
+
+        for (int i = 0; i < TetraConstants::numPattStars; i++) {
+            chosenIndIndices.push_back(centroidIndices[chosenCentroidIndices[i]]);
+            chosenStars.push_back(stars[chosenIndIndices[i]]);
+        }
+
+        // Compute Vec3 spatial vectors for each star chosen
+        // to be in the Pattern
+        std::vector<Vec3> pattStarVecs;  // size==numPattStars
+        for (const Star &star : chosenStars) {
+            Vec3 spatialVec = camera.CameraToSpatial(star.position);
+            pattStarVecs.push_back(spatialVec.Normalize());  // important!
+        }
+
+        // Compute angle between each pair of stars chosen to be in the Pattern
+        // If any angle > maxFov, then we should throw away this Pattern
+        // choice, since our database will not store it
+        bool angleAcceptable = true;
+        for (int i = 0; i < (int)pattStarVecs.size(); i++) {
+            Vec3 u = pattStarVecs[i];
+            for (int j = i + 1; j < (int)pattStarVecs.size(); j++) {
+                Vec3 v = pattStarVecs[j];
+                float angle = Angle(u, v);  // in radians
+                if (RadToDeg(angle) > maxFov) {
+                    angleAcceptable = false;
+                    break;
+                }
+            }
+        }
+
+        if (!angleAcceptable) {
+            // std::cerr << "Skipping combination, some angle is too large" << std::endl;
+            // Try again, some angle is too large
+            continue;
+        }
+
+        std::vector<double> pattEdgeRatios = TetraConstructPattern(pattStarVecs);
+
+        // Binning step
+        // Account for potential error in calculated values by
+        // testing hash codes in a range
+        std::vector<std::pair<int, int>> hcSpace;
+        for (double edgeRatio : pattEdgeRatios) {
+            std::pair<int, int> range;
+            int lo =
+                int((edgeRatio - TetraConstants::pattErrorRange) * TetraConstants::numPattBins);
+            lo = std::max(lo, 0);
+            int hi =
+                int((edgeRatio + TetraConstants::pattErrorRange) * TetraConstants::numPattBins);
+            hi = std::min(hi + 1, TetraConstants::numPattBins);
+            range = std::make_pair(lo, hi);
+            hcSpace.push_back(range);
+        }
+
+        std::set<std::vector<int>> finalCodes;
+        // Go through all the actual hash codes
+        for (int a = hcSpace[0].first; a < hcSpace[0].second; a++) {
+            for (int b = hcSpace[1].first; b < hcSpace[1].second; b++) {
+                for (int c = hcSpace[2].first; c < hcSpace[2].second; c++) {
+                    for (int d = hcSpace[3].first; d < hcSpace[3].second; d++) {
+                        for (int e = hcSpace[4].first; e < hcSpace[4].second; e++) {
+                            std::vector<int> code{a, b, c, d, e};
+                            std::sort(code.begin(), code.end());
+                            finalCodes.insert(code);
+                        }
+                    }
+                }
+            }
+        }
+
+        for (std::vector<int> code : finalCodes) {
+            // Must clear results list here
+            result.clear();
+
+            int hashIndex = KeyToIndex(code, TetraConstants::numPattBins, catLength);
+            // Get a list of Pattern Catalog rows with hash code == hashIndex
+            // One of these Patterns in the database could be a match to our constructed Pattern
+            std::vector<TetraPatt> matches = tetraDatabase.GetPatternMatches(hashIndex);
+
+            if (matches.size() == 0) {
+                // No matches for this pattern
+                continue;
+            }
+
+            bool alrFoundMatch = false;
+            int numMatches = 0;
+            for (TetraPatt matchRow : matches) {
+                // Construct the pattern we found in the Pattern Catalog
+                std::vector<int> catStarInds;
+                std::vector<Vec3> catStarVecs;
+                for (int star : matchRow) {
+                    int catInd = tetraDatabase.GetTrueCatInd(star);
+                    CatalogStar catStar = catalog[catInd];
+                    Vec3 catVec = catStar.spatial;
+                    catStarInds.push_back(catInd);
+                    catStarVecs.push_back(catVec);
+                }
+
+                std::vector<double> catEdgeRatios = TetraConstructPattern(catStarVecs);
+
+                bool skipMatchRow = false;
+                for (int i = 0; i < (int)catEdgeRatios.size(); i++) {
+                    float val = catEdgeRatios[i] - pattEdgeRatios[i];
+                    val = std::abs(val);
+                    // Test that our pattern and PC pattern roughly match up
+                    // For now, just compare edge ratios
+                    if (val > TetraConstants::pattMaxError) {
+                        skipMatchRow = true;
+                    }
+                }
+
+                // Very common to skip here, database Pattern is NOT a match to ours
+                if (skipMatchRow) {
+                    continue;
+                }
+
+                // Awesome! The PC pattern does match ours
+                // Reorder star vectors in both Patterns so we can match them up by ID
+
+                numMatches++;
 
+                if (alrFoundMatch)
+                    continue;  // If we already found a match, don't add new ID'd stars
+                alrFoundMatch = true;
+
+                Vec3 pattCentroid{0, 0, 0};
+                for (Vec3 pattStarVec : pattStarVecs) {
+                    pattCentroid = pattCentroid + pattStarVec;
+                }
+                pattCentroid = pattCentroid * (1.0 / pattStarVecs.size());
+
+                std::vector<float> pattRadii;
+                for (Vec3 pattStarVec : pattStarVecs) {
+                    pattRadii.push_back((pattStarVec - pattCentroid).Magnitude());
+                }
+
+                std::vector<int> sortedCentroidIndices =
+                    ArgsortVector<int, float>(chosenIndIndices, pattRadii);
+
+                std::vector<Vec3> pattSortedVecs =
+                    ArgsortVector<Vec3, float>(pattStarVecs, pattRadii);
+
+                Vec3 catCentroid{0, 0, 0};
+                for (Vec3 catStarVec : catStarVecs) {
+                    catCentroid = catCentroid + catStarVec;
+                }
+                catCentroid = catCentroid * (1.0 / catStarVecs.size());
+
+                std::vector<float> catRadii;
+                for (Vec3 catStarVec : catStarVecs) {
+                    catRadii.push_back((catStarVec - catCentroid).Magnitude());
+                }
+
+                std::vector<int> catSortedStarIDs = ArgsortVector<int, float>(catStarInds, catRadii);
+                std::vector<Vec3> catSortedVecs = ArgsortVector<Vec3, float>(catStarVecs, catRadii);
+
+                for (int i = 0; i < TetraConstants::numPattStars; i++) {
+                    int centroidIndex = sortedCentroidIndices[i];
+
+                    int catInd = catSortedStarIDs[i];
+
+                    // std::cerr << "Centroid Index: " << centroidIndex
+                    //           << ", Result StarID: " << resultStarID << std::endl;
+
+                    result.push_back(StarIdentifier(centroidIndex, catInd));
+                }
+
+                // Note: StarIdentifier wants the catalog INDEX, not the real star ID
+                // TODO: we used to consider early return here
+            }
+
+            // If we've looked at all possible matches and found a unique one, return it
+            // TODO: in practice, this seems to have almost no effect
+            // Surprisingly almost always, if the hashed pattern has at least 1 one match in the PC,
+            // then it can find a unique match
+            if (numMatches == 1) {
+                return result;
+            }
+        }
+    }
+
+    std::cerr << "TETRA FAIL: nothing could be matched" << std::endl;
+
+    result.clear();
+    return result;
+}
+
+StarIdentifiers DummyStarIdAlgorithm::Go(const unsigned char *, const Stars &stars,
+                                         const Catalog &catalog, const Camera &) const {
     StarIdentifiers result;
 
     unsigned int randomSeed = 123456;
@@ -252,7 +560,7 @@ std::vector<int16_t> ConsumeInvolvingIterator(PairDistanceInvolvingIterator it)
 /**
  * Given the result of a pair-distance kvector query, build a hashmultimap of stars to other stars
  * that appeared with it in the query.
- * 
+ *
  * The resulting map is "symmetrical" in the sense that if a star B is in the map for star A, then
  * star A is also in the map for star B.
  */
diff --git a/src/star-id.hpp b/src/star-id.hpp
index d82727c9..d42de9eb 100644
--- a/src/star-id.hpp
+++ b/src/star-id.hpp
@@ -3,9 +3,10 @@
 
 #include <vector>
 
+#include "camera.hpp"
 #include "centroiders.hpp"
+#include "databases.hpp"
 #include "star-utils.hpp"
-#include "camera.hpp"
 
 namespace lost {
 
@@ -16,10 +17,28 @@ namespace lost {
 class StarIdAlgorithm {
 public:
     /// Actualy perform the star idenification. This is the "main" function for StarIdAlgorithm
-    virtual StarIdentifiers Go(
-        const unsigned char *database, const Stars &, const Catalog &, const Camera &) const = 0;
+    virtual StarIdentifiers Go(const unsigned char *database, const Stars &, const Catalog &,
+                               const Camera &) const = 0;
+
+    virtual ~StarIdAlgorithm(){};
+};
+
+/**
+ * Tetra star identification algorithm based on the one originally proposed by Brown & Stubis
+ * Utilizes hashing of 4-star patterns for quick lookup time, decently high accuracy
+ */
+class TetraStarIdAlgorithm : public StarIdAlgorithm {
+   public:
+    StarIdentifiers Go(const unsigned char *database, const Stars &centroids,
+                       const Catalog &catalog, const Camera &) const;
+
+   private:
+    // static const int numPattStars;
+    // // Do NOT modify the following parameters unless you know what you're doing
+    // static const int numPattBins;
+    // const float pattErrorRange = 0.001;
+    // const float pattMaxError = 0.001;
 
-    virtual ~StarIdAlgorithm() { };
 };
 
 /// A star-id algorithm that returns random results. For debugging.
diff --git a/src/star-utils.cpp b/src/star-utils.cpp
index af12dd9b..f2a0654f 100644
--- a/src/star-utils.cpp
+++ b/src/star-utils.cpp
@@ -1,14 +1,22 @@
 #include "star-utils.hpp"
 
-#include <math.h>
 #include <assert.h>
+#include <math.h>
+
 #include <algorithm>
+#include <iostream>  // TODO: remove
+#include <vector>
 #include <set>
 
 #include "serialize-helpers.hpp"
 
 namespace lost {
 
+const int TetraConstants::numPattStars = 4;
+const int TetraConstants::numPattBins = 50;
+const float TetraConstants::pattErrorRange = 0.001;
+const float TetraConstants::pattMaxError = 0.001;
+
 // brightest star first
 bool CatalogStarMagnitudeCompare(const CatalogStar &a, const CatalogStar &b) {
     return a.magnitude < b.magnitude;
@@ -17,39 +25,53 @@ bool CatalogStarMagnitudeCompare(const CatalogStar &a, const CatalogStar &b) {
 Catalog NarrowCatalog(const Catalog &catalog, int maxMagnitude, int maxStars, float minSeparation) {
     Catalog result;
     for (int i = 0; i < (int)catalog.size(); i++) {
+        // Higher magnitude implies the star is dimmer
         if (catalog[i].magnitude <= maxMagnitude) {
             result.push_back(catalog[i]);
         }
     }
+    // TODO: sort regardless?
+    std::sort(result.begin(), result.end(), CatalogStarMagnitudeCompare);
 
+    // TODO: bug? Why remove both i and j
     // remove stars that are too close to each other
     std::set<int> tooCloseIndices;
     // filter out stars that are too close together
     // easy enough to n^2 brute force, the catalog isn't that big
     for (int i = 0; i < (int)result.size(); i++) {
-        for (int j = i+1; j < (int)result.size(); j++) {
-            if (AngleUnit(result[i].spatial, result[j].spatial) < minSeparation) {
+        for (int j = i + 1; j < (int)result.size(); j++) {
+            if (AngleUnit(result[i].spatial.Normalize(), result[j].spatial.Normalize()) < minSeparation) {
                 tooCloseIndices.insert(i);
                 tooCloseIndices.insert(j);
             }
         }
     }
 
-    // Erase all the stars whose indices are in tooCloseIndices from the result.
-    // Loop backwards so indices don't get messed up as we iterate.
+    // // Erase all the stars whose indices are in tooCloseIndices from the result.
+    // // Loop backwards so indices don't get messed up as we iterate.
     for (auto it = tooCloseIndices.rbegin(); it != tooCloseIndices.rend(); it++) {
         result.erase(result.begin() + *it);
     }
 
     // and finally limit to n brightest stars
     if (maxStars < (int)result.size()) {
-        std::sort(result.begin(), result.end(), CatalogStarMagnitudeCompare);
+        // std::sort(result.begin(), result.end(), CatalogStarMagnitudeCompare);
         result.resize(maxStars);
     }
 
     return result;
 }
 
+int KeyToIndex(std::vector<int> key, int binFactor, long long maxIndex) {
+    const long long MAGIC_RAND = 2654435761;
+    long index = 0;
+    for (int i = 0; i < (int)key.size(); i++) {
+        index += key[i] * std::pow(binFactor, i);
+    }
+
+    return ((index % maxIndex) * (MAGIC_RAND % maxIndex)) % maxIndex;
+}
+
 /// Return a pointer to the star with the given name, or NULL if not found.
 Catalog::const_iterator FindNamedStar(const Catalog &catalog, int name) {
     for (auto it = catalog.cbegin(); it != catalog.cend(); ++it) {
diff --git a/src/star-utils.hpp b/src/star-utils.hpp
index 0af225bd..3566b317 100644
--- a/src/star-utils.hpp
+++ b/src/star-utils.hpp
@@ -1,6 +1,7 @@
 #ifndef STAR_UTILS_H
 #define STAR_UTILS_H
 
+#include <utility>  // pair
 #include <vector>
 
 #include "attitude-utils.hpp"
@@ -131,6 +132,18 @@ float MagToBrightness(int magnitude);
  */
 Catalog NarrowCatalog(const Catalog &, int maxMagnitude, int maxStars, float minSeparation);
 
-}
+/// Hash function for Tetra, convert from star pattern representation into index in pattern catalog
+int KeyToIndex(std::vector<int> key, int binFactor, long long maxIndex);
+
+struct TetraConstants{
+    static const int numPattStars;
+    // Do NOT modify the following parameters unless you know what you're doing
+    static const int numPattBins;
+    static const float pattErrorRange;
+    static const float pattMaxError;
+};
+
+
+}  // namespace lost
 
 #endif
diff --git a/test/scripts/pyramid-incorrect.sh b/test/scripts/pyramid-incorrect.sh
index 889c3ce9..0306b27f 100755
--- a/test/scripts/pyramid-incorrect.sh
+++ b/test/scripts/pyramid-incorrect.sh
@@ -24,7 +24,7 @@ test -e pyramid-incorrect.dat || ./lost database \
   --kvector-distance-bins 10000 \
   --output pyramid-incorrect.dat
 
-for _ in $(seq "${1:-10}"); do
+for _ in $(seq "${1:-100}"); do
   x_res=$(rand_int 100 1000)
   y_res=$(rand_int 100 1000)
   fov=$(rand_int 5 70)
diff --git a/test/scripts/tetra.sh b/test/scripts/tetra.sh
new file mode 100755
index 00000000..eb865a53
--- /dev/null
+++ b/test/scripts/tetra.sh
@@ -0,0 +1,139 @@
+#!/usr/bin/env bash
+
+# Tetra testing
+
+th="0.5"
+
+totalTime=0
+# CHANGE THIS, n = number of images to try
+n=1
+totalCorrect=0
+totalIncorrect=0
+totalNoStars=0
+
+totalAttCorrect=0
+totalAttIncorrect=0
+
+
+# compare actual to expected
+function cmpAE() {
+  ours="$1"
+  actual="$2"
+
+  # echo "ours is $ours, actual is $actual"
+
+  if (( $(echo "$ours < $actual" | bc -l) )); then
+    diff=$(echo "$actual - $ours" | bc -l )
+  elif (( $(echo "$ours > $actual" | bc -l) )); then
+    diff=$(echo "$ours - $actual" | bc -l )
+  else
+    # Because 0=true, 1=false in bash rip
+    return 0
+  fi
+
+  diff=$(echo "$diff % 360" | bc)
+
+  if (( $(echo "$diff < $th" | bc -l) || $(echo "(360 - $diff) < $th" | bc -l) )); then
+    return 0
+  else
+    return 1
+  fi
+}
+
+
+# params: low and hi
+function rand_int {
+  lo=$1
+  hi=$2
+  echo $((lo + (RANDOM % (hi - lo + 1))))
+}
+
+# This works fine
+# create database if not exists
+# test -e tetra-incorrect.dat || ./lost database \
+#   --min-mag 7 \
+#   --tetra \
+#   --tetra-max-angle 12 \
+#   --output tetra-incorrect.dat
+
+for _ in $(seq "${1:-$n}"); do
+  # ra=$(rand_int 0 359)
+  # de=$(rand_int -89 89)
+  # roll=$(rand_int 0 359)
+  # fov=$(rand_int 10 60)
+
+  fov=16
+  ra=218
+  de=12
+  roll=302
+
+  # TODO: separate image generation and actual pipeline run for accurate timing
+  # set -x
+  lost_output=$(
+    ./lost pipeline \
+      --generate 1 \
+      --generate-x-resolution 1024 \
+      --generate-y-resolution 1024 \
+      --fov "$fov" \
+      --generate-reference-brightness 100 \
+      --generate-spread-stddev 1 \
+      --generate-read-noise-stddev 0.05 \
+      --generate-ra "$ra" \
+      --generate-de "$de" \
+      --generate-roll "$roll" \
+      --database my-database-small-3.dat \
+      --centroid-mag-filter 5 \
+      --star-id-algo tetra \
+      --compare-star-ids \
+      --attitude-algo dqm \
+      --print-attitude
+  )
+
+  # --centroid-algo cog \
+  # set +x
+  num_incorrect_stars=$(grep -oP "(?<=starid_num_incorrect )\\d+" <<<"$lost_output")
+  num_correct_stars=$(grep -4oP "(?<=starid_num_correct )\\d+" <<<"$lost_output")
+
+  raOurs=$(grep -oP "(?<=attitude_ra )[-+]?[0-9]*\.?[0-9]*[eE]?[-+]?[0-9]*" <<<"$lost_output")
+  deOurs=$(grep -oP "(?<=attitude_de )[-+]?[0-9]*\.?[0-9]*[eE]?[-+]?[0-9]*" <<<"$lost_output")
+  rollOurs=$(grep -oP "(?<=attitude_roll )[-+]?[0-9]*\.?[0-9]*[eE]?[-+]?[0-9]*" <<<"$lost_output")
+
+  echo "Fov: $fov"
+  echo "Real: $ra, $de, $roll vs Calculated: $raOurs, $deOurs, $rollOurs"
+
+  # TODO: 360 and 0 should be deemed equivalent
+
+  if cmpAE "$raOurs" "$ra" && cmpAE "$deOurs" "$de" && cmpAE "$rollOurs" "$roll"; then
+    echo "TEST: [SUCCESS]: ATT All clear!"
+    totalAttCorrect=$(( $totalAttCorrect + 1 ))
+  else
+    echo "TEST: [ERROR ATT]: Incorrect attitude calculation!"
+    totalAttIncorrect=$(( $totalAttIncorrect + 1 ))
+  fi
+
+  if ((num_correct_stars > 0 && num_incorrect_stars == 0)); then
+    echo "TEST: [SUCCESS]: All clear!"
+    totalCorrect=$(( $totalCorrect + 1 ))
+  elif ((num_correct_stars == 0)); then
+    echo "TEST: [ERROR]: No correct stars identified!"
+    # echo "Values: $ra, $de, $roll"
+    totalNoStars=$(( $totalNoStars + 1 ))
+    # exit 1
+  elif ((num_incorrect_stars > 0)); then
+    echo "TEST: [ERROR]: Incorrect stars identified!"
+    # echo "Values: $ra, $de, $roll"
+    totalIncorrect=$(( $totalIncorrect + 1 ))
+    # exit 1
+  else
+    echo "TEST: [ERROR]: Unknown error!"
+    exit 1
+  fi
+done
+
+echo "Total n = $n"
+echo "Total correct: $totalCorrect"
+echo "Total incorrect: $totalIncorrect"
+echo "Total no stars: $totalNoStars"
+
+echo "Total attitude correct: $totalAttCorrect"
+echo "Total attitude incorrect: $totalAttIncorrect"