Compressive Sensing with Local Geometric Features
Rishi Gupta
MIT CSAILPiotr Indyk
MIT CSAILEric Price
MIT CSAILYaron Rachlin
MIT Lincoln Laboratory

We propose a framework for compressive sensing of images with local distinguishable objects, such as stars, and apply it to solve a problem in celestial navigation.
Specifically, let $x \in \R^N$ be an $N$-pixel image, consisting of a small number of local distinguishable objects plus noise. Our goal is to design an $m\times N$ measurement matrix $A$ with $m \ll N$, such that we can recover an approximation to $x$ from the measurements $Ax$.

We construct a matrix $A$ and recovery algorithm with the following properties: (i) if there are $k$ objects, the number of measurements $m$ is $O((k\log N)/(\log k))$, undercutting the best known bound of $O(k \log (N/k))$ (ii) the matrix $A$ is very sparse, which is important for hardware implementations of compressive sensing algorithms, and (iii) the recovery algorithm is empirically fast and runs in time polynomial in $k$ and $\log(N)$.

We also present a comprehensive study of the application of our algorithm to attitude determination, or finding one's orientation in space. Spacecraft typically use cameras to acquire an image of the sky, and then identify stars in the image to compute their orientation. Taking pictures is very expensive for small spacecraft, since camera sensors use a lot of power. Our algorithm optically compresses the image before it reaches the camera's array of pixels, reducing the number of sensors that are required.

§ INTRODUCTION

§.§ Compressive sensing

Traditional approaches to image acquisition first capture an entire $N$-pixel image and then process it for compression, transmission, or storage. Often, the image is captured at very high fidelity, only to be immediately compressed after digitization. In contrast, compressive sensing uses prior knowledge about a signal to obtain a compressed representation directly, by acquiring a small number of nonadaptive linear measurements of the signal in hardware [8, 11].
Formally, for an image represented by a vector $x$, we acquire the measurement vector, or sketch, $Ax$, where $A$ is an $m \times N$ matrix. The advantage of this architecture is that it uses fewer sensors, and therefore can be cheaper and use less energy than a conventional camera [9, 14, 31].

In order to reconstruct the image $x$ from the lower-dimensional sketch $Ax$, we assume that the image $x$ is $k$-sparse for some $k$ (i.e., has at most $k$ non-zero coordinates) or at least be well-approximated by a
$k$-sparse vector. Then, given $Ax$, one finds an approximation to
$x$ by performing sparse recovery.  The problem is
typically defined as follows: construct a matrix $A$ such that, for
any signal $x$, we can recover a vector $\hat{x}$ from $Ax$ satisfying

\begin{equation}
\label{e:lplq}
\norm{1}{x-\hat{x}} \le C \cdot \err^1_k(x),
\end{equation}

where $\err^1_k(x)=\min_{k\mbox{-sparse } x'}  \norm{1}{x-x'}$ and $C$ is the approximation factor.
Note that if $x$ is $k$-sparse, then $\err^1_k(x)=0$, and therefore $\hat{x}=x$.
Sparse recovery also has applications to other areas, such as data
stream computing [28, 20].

The problem of designing matrices $A$ and corresponding recovery algorithms has been a subject of extensive study over the last few years, with the goal of designing schemes that enjoy good compression rate (i.e., low values of $m$) as well as good algorithmic properties such as low encoding complexity and quick recovery times. Low encoding complexity is often achieved by using matrices that are binary (entries chosen from $\{0,1\}$ or $\{-1,1\}$), or that have low column sparsity. Column sparsity is the average number of non-zero entries per column, namely the average number of buckets into which each coordinate of the signal gets split.
It is known by now that there exist binary
matrices $A$ and associated recovery algorithms that produce
approximations $\hat{x}$ satisfying Eq. <ref> with constant approximation factor $C$ and sketch
length $m=O(k \log (N/k))$.
In particular, a
random Bernoulli matrix [8] or a random
binary matrix with column sparsity $O(\log(N/k))$ [2] has
this property with overwhelming probability.
It is also known that this sketch length is asymptotically optimal [10, 13]. See [15] for an overview of compressive sensing using matrices with low column sparsity, along with [3] for a newer algorithm that we run experiments against (Section <ref>).

§.§ Attitude determination

Spacecraft determine their attitude, or 3-axis orientation, by taking pictures of the sky ahead of them and identifying stars in the image. This function is encapsulated in a star tracker, which is essentially a digital camera connected to a processor. To acquire the initial attitude, the camera

* Takes a picture of the sky.

* Identifies a set of starlike objects in the picture, and computes the centroid of each object.

* Tries to match triangles and quadrilaterals formed by subsets of the centroids to an onboard database. A match provides approximate attitude information.

* Uses the onboard database to determine a set of stars that it expects are in the picture, along with their approximate locations. Refines the attitude information by centroiding those stars as well.

Most of the time, a star tracker has knowledge about its approximate attitude, either from the previous attitude computation or from other sensors. In that case, it goes from Step 1 directly to Step 4, in what we call tracking mode.  [22] has an overview of the process.

There are two types of sensors used in star tracker cameras, CCD (charge-coupled device), and CMOS (complimentary metal-oxide semiconductor). CCD sensors have low noise and capture a high fraction of the incoming signal, but have power usage and manufacturing costs that are super-linear in the number of pixels and that are high in absolute terms. On the other hand, CMOS sensors use little power, are cheap to manufacture, and allow random access to pixel values (important for tracking mode), but capture less of the incoming signal, and are very noisy. Most spacecraft use CCD cameras, but smaller or cheaper spacecraft use CMOS, taking a factor of 10 or higher hit in precision.

§.§ Motivation for a new algorithm

Ideally, we would like to have a camera with the precision of a high pixel CCD camera, but without the extra power and manufacturing costs that drive small spacecraft to use CMOS. The pictures taken by star trackers are naturally very sparse, in that most pixels are either empty or contain small stars not used in Steps 3 or 4. Also, the algorithms for Step 3 are very robust, and can tolerate a substantial fraction of bogus centroid information [27]. A compressive sensing solution would optically compress the incoming picture, to reduce the size of the CCD array.

However, a standard compressive sensing solution runs into several major problems. First, the $L_1$ mass of the small stars is large compared to the $L_1$ mass of the signal. In other words, $\err_k^1(x) = O(\|x\|_1)$, and so Eq. <ref> (Section <ref>) gives no guarantee at all. Second, the signal to noise ratio on star trackers is already low, and each non-zero entry of the sensing matrix $A$ will add substantial signal-independent noise to the final measurement.[18] Finally, each star is spread over multiple pixels, and makes only a small contribution to some of them. These pixels are needed to find the centroid of the star properly, but compressive sensing recovery techniques are only designed to recover the biggest pixels well.

We address many of these concerns by focusing on a compressive sensing algorithm where $A$ is very sparse. Compressive sensing algorithms with sparse measurement matrices are of general interest as well. Potential advantages in non-star tracking electronic compressive imagers include reduced interconnect complexity [25], low memory requirements for storing the measurement matrix, and gain in image acquisition speed due to reduced operations.

Unfortunately, it is known [29] that any deterministic scheme with guarantee as in Eq. <ref> requires column sparsity of $\Omega(\log(N/k))$. In the randomized case, where $A$ is a random variable, and Eq. <ref> is required to hold only with constant probability over the choice of $A$, the same paper shows that any binary matrix $A$ must have column sparsity as stated.

In this paper we overcome the above limitations by employing a two-fold approach. First, we consider a class of images where the $k$ large coefficients or $k$ local objects can be distinguished from each other. Second, we relax the recovery guarantee, by requiring that only a constant fraction of the objects are recovered correctly, and only with constant probability.

§.§ Model description
Our model for sparse images is motivated by astronomical imaging, where an image contains a small number of distinguishable objects (e.g., stars) plus some noise. We model each object as an image contained in a small $w \times w$ bounding box, for some $w = O(1)$. The image is constructed by placing $k$ objects in the image in an arbitrary fashion, subject to a minimum separation constraint. The image is then modified by adding noise. We formalize the notions of minimum separation constraint, distinguishability, and noise in the rest of this section. Some of the definitions below are illustrated in Appendix <ref>.

Let $x$ be an $N$-dimensional real vector, and assume $N=n^2$ for an integer $n$. We will treat $x$ both as a vector and as an $n \times n$ matrix, with entries $x[i,j]$ for $i, j \in [n]$.  An object $o$ is a $w \times w$ real matrix. Let $\mathcal{O}=\{o_1, \ldots, o_k\}$ be a sequence of $k$ objects, and let $\mathcal{T}=\{t_1, \ldots, t_k\}$ be a sequence of translations in $x$, i.e., elements from $[n-w]^2$. We say that $\mathcal{T}$ is valid if for any $i \neq j$ the translations $t_i$ and $t_j$ do not collide, i.e., we have $\|t_i-t_j\|_{\infty} \ge w'$ for some constant separation parameter $w' = \Omega(w)$. For $o \in \mathcal{O}$ and $t=(t_x,t_y) \in \mathcal{T}$, we define $t(o)$ to be a $w \times w$ matrix indexed by $\{t_x, \ldots, t_x+w-1\} \times \{t_y, \ldots, t_y+w-1\}$. Using somewhat sloppy notation, the ground truth image is then defined as $x=\sum_i {t_i(o_i)}$.

During our algorithm, we impose a grid $G$ on the image with cells of size $w' \times w'$. Let $x_c$ be the image (i.e., an $w'^2$-dimensional vector) corresponding to cell $c$. We then use a Does it have to be linear?projection $F$ that maps each sub-image $x_c$ into a feature vector $F(x_c)$. If $y\subset x_c$ for some cell $c$ and some set of pixels $y$, we use $F(y)$ to denote $F(x_c)$ after the entries of $x_c\setminus y$ are set to 0. If $y$ is not a cell and not contained in a cell, we leave $F(y)$ undefined.

The distinguishability property we assume is that for any two distinct $o,o'$ from the objects $\mathcal{O} \cup \{\emptyset\}$, and for any two translations $t$ and $t'$, we have $\|F(t(o))-F(t'(o'))\|_\Gamma > T$ (when it is defined) for some threshold $T>0$ and some norm $\|\cdot\|_\Gamma$. In other words, different objects need to look different under $F$. For concreteness, the features we exploit in the experimental section are the magnitude (the sum of all pixels in the cell) and centroid (the sum of all pixels in the cell, weighted by pixel coordinates), since the magnitudes of stars follow a power law, and the centroid of a star can be resolved to $.15$ times the width of a pixel in each dimension (Section <ref>). The distinguishability constraint is what ultimately allows us to undercut the usual lower bound by a factor of $\log k$.

The observed image $x'$ is equal to $x+\mu$, where $\mu$ is a noise vector. The threshold $T$ determines the total amount of noise that the algorithm tolerates. Specifically, let $\|\mu\|_F = \sum_c \|F(\mu_c)\|_\Gamma$, where $\mu_c$ is the noise corresponding to cell $c$. We assume that $\|\mu\|_F < \gamma kT$ for some small constant $\gamma>0$, and make no other assumptions about the noise.

§.§ Results and techniques

§.§.§ Theoretical result
Assume sparsity parameter $k \ge
C\log N$ for some constant $C$, and prior knowledge of $k$ and the
distinguishability parameter $T$. We construct a distribution over
random binary $m \times N$ matrices $A$, such that given $Ax'$ for
$x'$ described above, we recover, with constant probability, a set $D$
of $k$ cells, such that at least $k/2$ of the cells fully containing
an object are included in $D$[From this, a simple min or
median process can be used to recover an approximation to $x_c$ for
any $c \in D$. See Section II.A of [15] for an
explanation of the technique.].  The matrix has column sparsity
$O(\log_k N)$, and has $m=O(k \log_k N)$ rows. Note that we trade off
column sparsity and compression. If (say) $k=N^{1/2}$, then the column
sparsity is constant, and $m=O(N^{1/2})$.  The running time for the
recovery procedure is $O(k^3 \log^3_k n)$.

§.§.§ Empirical result
We implement a standard attitude determination routine, with the picture acquisition step replaced with a simplified version of the theoretical algorithm. Our algorithm performs better recovery on small numbers of measurements and is orders of magnitude faster than comparable compressive sensing methods.

§.§.§ Our techniques
Our construction of the measurement matrix resembles those of other algorithms for sparse matrices, such as
Count-Sketch [4] or Count-Min [6]: we “hash” each cell $c$ into each of $s=O(\log_k N)$ arrays of $q=O(k)$ “buckets”, and sum all the cells hashed to the same bucket. Each bucket defines one measurement of $w'^2$ pixels, which gives $m=O(k \log_k N)$.  Hashing is done using either the Reed-Solomon code or the Chinese Remainder code[Note that our use of the Chinese Remainder code does not incur any additional polylogarithmic factors.].

The recovery process is based on the following novel approach. For simplicity, assume for now that the image contains no noise, and ignore the effect of two different objects being hashed to the same bucket. In this case, all buckets containing distinct objects are distinguishable from each other. Therefore, we can group non-empty buckets into $k$ clusters of size $s$, with each cluster containing buckets with a single value. Since $q^s > N$, each cluster of buckets uniquely determines the cell in $x$ containing the object in those buckets.

In order to make this approach practical, however, we need to make it robust to errors. The errors are due to distinct objects being hashed to the same bucket, the noise vector $\mu$, and the grid cutting objects into pieces. Because of these issues, the clustering procedure aims to find clusters containing elements that are close to each other, rather than equal, and the procedure allows for some small fraction of outliers [5]. For this purpose, we use the approximation algorithm for the $k$-center problem with outliers, which correctly clusters a constant fraction of the buckets. To handle the buckets that are grouped incorrectly, we construct our hash function using a constant rate error-correcting code [17].

§ THEORETICAL RESULTS

A graphical representation of the algorithm is presented in Appendix <ref>. Our scheme works by “hashing” each cell $c$ into $s = O(\log_k N)$ different arrays of size $O(k)$.  We can think of this as a mapping $g$ from
$[N]$ to $[O(k)]^s$.  As long as each character of the mapping is
approximately pairwise independent, then (in expectation) most of the
$k$ objects will be alone in most of the array locations they map to.
Our reconstruction algorithm clusters the values in the cells,
giving us a noisy version $y'$ of the true codeword $y = g(c)$ with a
constant fraction of errors.  We then efficiently decode from
$y'$ to $c$.

The first and second sections below establish families $\mathcal{G}$ from which we will draw mappings $g$. The third uses $\mathcal{G}$ to construct a distribution over measurement matrices $A$, and the fourth presents an associated recovery algorithm. The main result is stated in Theorem <ref>.

§.§ Definitions and preliminaries

We need an efficient error correcting code that is also
approximately pairwise independent in each character. This section gives
precise definitions of our requirements, and the next section gives two codes that achieve them.

A hash family $\Hc$ of functions $h \colon A \to B$ is
pairwise-independent if, for any $x_1, x_2 \in A$ and $y_1,
y_2 \in B$ with $x_1 \neq x_2$, we have $\Pr_{h \in
\Hc}[h(x_1) = y_1 \cap h(x_2) = y_2] =
\frac{1}{\abs{B}^2}$.

For any prime $P \geq N$, the function family $\Hc_P: ax + b \pmod{P}$
for $a, b \in [P]$ is pairwise independent when viewed as a set of
functions from $[N]$ to $[P]$.

In many of our applications the range $B$ is the product of $s$
“symbols” $B_1 \times \dotsb \times B_s$.  For a function $f\colon A
\to B$ and $i \in [s]$, we use $f_i(x)$ to denote the $i$th coordinate
of $f$.  When $B$ is a product space, we will sometimes settle for a
weaker notion of pairwise independence.  Rather than requiring
pairwise independence for the whole range, we only require approximate
pairwise independence in each coordinate:

Let $B = B_1 \times \dotsb \times B_s$.  A hash family $\Hc$ of
functions $h \colon A \to B$ is coordinatewise
$C$-pairwise-independent if, for all $i \in [s]$, any $x_1 \neq
x_2 \in A$, and all $y_1, y_2 \in B_i$, we have $\Pr_{h \in
\Hc}[h_i(x_1) = y_1 \cap h_i(x_2) = y_2]
\leq \frac{C}{\abs{B_i}^2}$.

Let $B = B_1 \times \dotsb \times B_s$. A function $f \colon A \to
B$ is $C$-uniform if, for all $i \in [s]$ and all $y \in
B_i$, $\Pr_{x \in A}[f_i(x) = y] \leq \frac{C}{\abs{B_i}}$.

For any function $f\colon B \to D$ and family $\Hc$ of functions
$h\colon A \to B$, we use $f \circ \Hc$ to denote the family of $A \to
D$ functions $\{g(x) := f(h(x)) \mid h \in \Hc\}$.

If $\Hc$ is pairwise-independent and $f$ is $C$-uniform, then $f
\circ \Hc$ is coordinatewise $C^2$-pairwise-independent.

Let $\Hc$ be a family of functions $A \to B$ and let $f\colon B \to D =
D_1 \times \dotsb \times D_s$.  Then for any $i \in [s]$, any $x_1
\neq x_2 \in A$, and all $y_1, y_2 \in D_i$ we have:

\begin{align*}
&  \Pr_{h \in \Hc}[f_i(h(x_1)) = y_1 \cap f_i(h(x_2)) = y_2]\\
=& \sum_{z_1, z_2 \in B} \Pr_{h \in \Hc}[h(x_1) = z_1 \cap h(x_2) = z_2 \cap f_i(z_1) = y_1 \cap f_i(z_2) = y_2]\\
=& \sum_{z_1, z_2 \in B} \frac{1}{\abs{B}^2}\Pr\,[f_i(z_1) = y_1 \cap f_i(z_2) = y_2]\\
=& \Pr_{z_1, z_2 \in B} [f_i(z_1) = y_1 \cap f_i(z_2) = y_2]\\
=& \Pr_{z_1 \in B} [f_i(z_1) = y_1] \Pr_{z_2 \in B}[f_i(z_2) = y_2]\\
\leq& \frac{C^2}{\abs{B_i}^2}
\end{align*}

as desired.

We say that a function $f\colon A \to B$ for $B = B_1 \times \dotsb
\times B_s$ is an error-correcting code of distance $d$ if,
for any two distinct $x_1, x_2 \in A$, $f(x_1)$ and $f(x_2)$ differ
in at least $d$ coordinates.

We say that $f$ is efficiently decodable if we have an
algorithm $\widetilde{f}^{-1}$ running in $\log^{O(1)} |B|$ time with $\widetilde{f}^{-1}(y) =
x$ for any $x \in A$ and $y \in B$ such that $f(x)$ and $y$ differ
in fewer than $d/2$ coordinates.

Recall the hash family $\Hc_P : ax + b \pmod P$ of functions $[N] \to

If $f$ is an efficiently decodable error-correcting code with
distance $d$, then so is $f \circ h$ for every $h \in \Hc_P$ with $a
\neq P$.

Since $a \neq P$, there exists an $a^{-1}$ modulo $P$, and we can
efficiently compute it.  Hence $h$ is injective, so $f \circ h$ is
an error-correcting code of distance $d$.  Furthermore, $(f \circ
h)^{-1}(x) = a^{-1} (f^{-1}(x) - b) \pmod{P}$ is efficiently computable.

We say that a family $\Gc$ of functions $g \colon A \to B_1\times
\dotsb \times B_s$ is an $(N, s, d)_q$-independent code if $\Gc$
is coordinatewise 4-pairwise independent, $q \leq \abs{B_i} \leq 2q$
for all $i \in [s]$, $\abs{A} \geq N$, and with probability at least
$1 - 1/N$ over $g \in \Gc$ we have that $g$ is efficiently decodable
with distance $d$.

Claims <ref> and <ref> give the
following lemma:

If $f\colon [P] \to B_1 \times \dotsb \times B_s$ is $2$-uniform and
efficiently decodable with distance $d$, and $q \leq \abs{B_i} \leq 2q$ for
all $i$, then $f \circ \Hc_P$ is a $(N, s, d)_q$-independent code.

We now show that $(N, s, d)_q$-independent codes have few
collisions in expectation.

Suppose $g\colon A \to B_1 \times \dotsc \times B_s$ is drawn from a
$(N, s, d)_q$-independent code.  Let $S, S' \subset A$.  Define
the set of “colliding” symbols

\[
X = \{(a, i) \mid a\in S, i \in [s], \exists a' \in S' \text{
s.t. } g_i(a) = g_i(a'), a \neq a'\}.
\]

With probability at least $7/8$, $\abs{X} \leq 32\abs{S}\abs{S'}s/q$.

We observe that

\begin{align*}
\E[\abs{X}] &= \sum_{i \in [s]} \sum_{a \in S} \Pr[(a,i) \in X]\\
&\leq  \sum_{i \in [s]} \sum_{a \in S} \sum_{\substack{a' \in S'\\a' \neq a}} \Pr[g_i(a) = g_i(a')]\\
&=  \sum_{i \in [s]} \sum_{a \in S} \sum_{\substack{a' \in S'\\a' \neq a}} \sum_{z \in B_i} \Pr[g_i(a) = z \cap g_i(a') = z]\\
&\leq  \sum_{i \in [s]} \sum_{a \in S} \sum_{\substack{a' \in S'\\a' \neq a}} \sum_{z \in B_i} \frac{4}{\abs{B_i}^2}\\
&\leq s\abs{S}\abs{S'} 4/ q.
\end{align*}

Hence, by Markov's inequality, $\abs{X} \leq
32\abs{S}\abs{S'}s/q$ with probability at least $7/8$.

§.§ Two code constructions

We explicitly give two $(N, s, s-r)_q$-independent codes.  Both are
achievable for any parameters with $2N < q^r$ and $s < q / \log q$
(and the first code allows any $s < q$).  We let $P$ be a prime in
$\{\frac{1}{2}q^r, \dotsc, q^r\}$.

§.§.§ Reed-Solomon code

Let $q\ge s$. The Reed-Solomon code $f_{RS} \colon [q^r] \to [q]^s$ is defined for
$f(x)$ by (i) interpreting $x$ as an element of $\mathbb{F}_q^r$, (ii) defining
$\chi_x \in \mathbb{F}_q[\xi]$ to be the degree $r-1$ polynomial with
coefficients corresponding to $x$, and (iii) outputting $f(x) =
(\chi_x(1), \dotsc, \chi_x(s))$.  It is well known to have distance
$s-r$ and to be efficiently decodable [21].

Let $f: [P] \to [q]^s$ be the restriction of $f_{RS}$ to $[P]$. Then
$f$ is $2$-uniform, so $\Gc_{RS} = f \circ \Hc_P$ is a $(N, s,
s-r)_q$-independent code.

Basic facts about polynomials give that $f_{RS}$ is $1$-uniform.
Since $P \geq q^r/2$, $f$ is $2$-uniform.
Lemma <ref> then gives the result.

§.§.§ Chinese remainder theorem (CRT) code

Let $p_1, \dotsc, p_s \in [q, 2q]$ be distinct primes; note that the asymptotic distribution of prime numbers
implies $q/\log q =\Omega(s)$.  Hence for any $x \in [N]$, any $r$ of
the residues mod $p_1, \dotsc, p_s$ uniquely identify $x$.  The CRT
code $f_{CRT}\colon [P] \to [p_1] \times \dotsc \times [p_s]$ is
defined by taking the residues modulo each prime.  It has distance
$s-r$ and is efficiently decodable [16].

The CRT code $f_{CRT}$ is $2$-uniform.  Hence $\Gc_{CRT} = f_{CRT}
\circ \Hc_P$ is a $(N, s, s-r)_q$-independent code.

Let $i \in [s]$.  The projection of $f_{CRT}(x)$ onto its $i$th
coordinate is $x \bmod p_i$.  Hence over the domain $[P]$, the ratio
between the likelihood of the most common and the least common
values in the range is $\frac{\ceil{P / p_i}}{\floor{P / p_i}} \leq
2$.  Thus $f_{CRT}$ is $2$-uniform, and
Lemma <ref> gives the result.

§.§ The measurement matrix

In this section we present the measurement matrix $A$. A graphical
representation of the measurement process is presented on the first
page of Appendix <ref>. Let $\mathcal{O}=\{o_1, \dotsc,
o_k\}$ be a sequence of $k$ features, and let $\mathcal{T}=\{t_1,
\dotsc, t_k\}$ be a sequence of (non-colliding) translations in
$x$. Let $\mu$ be the noise vector, and let $x'$ be the noisy
image. Finally, let $\alpha,\beta,\delta, \eta>0$ be (small)
constants whose values will be determined in the course of the

At the beginning, we impose a square grid $G$ with  $w' \times w'$ cells on the image $x'$, such that $w'=w/\alpha$. The grid is shifted by a vector $v$ chosen uniformly at random from $[w']^2$.

Let $S'$ be the set of cells that intersect or contain some object $t_i(o_i)$, and $S \subset S'$ be the set of cells that fully contain some object $t_i(o_i)$.
Observe that a fixed object is fully contained in some cell with probability $(1-w/w')^2 > 1-2\alpha$, since each axis of the grid intersects the object with probability $w/w'$. This implies that the expected number of cells in $S'-S$ is at most
$2 \alpha k$, and by Markov's inequality $|S'-S| \le 16\alpha k$ with probability $7/8$.
From now on, we will assume the latter event holds. Let $k'=|S'|$. We choose $\alpha>0$ such that $k' \le 2k$.

Our measurement matrix $A$ is defined by the following linear
mapping. Let $G$ denote the set of cells.  Let $g \colon G \to B = B_1
\times \dotsb \times B_s$ be drawn from a $(N, s,
4(3\delta+\beta)s)_q$-independent code (such as either $\Gc_{RS}$ or
$\Gc_{CRT}$).  Moreover, we require that $k/q \le \eta$; such a code
is achievable per Section <ref> with
$s=\Theta(\log_k N)$ as long as $k > C \log N$ for some constant $C$
(such that both $q^{(1-4(3\delta + \beta))s} > k^{s/2} > 2N$ and $s <
\log N / \log k \leq q / \log q$).  For each $i =1, \dotsc, s$, we define a
$\abs{B_i}$-dimensional vector $z^i$ whose entries are elements in
$\R^{w'^2}$, such that for any $j$

\[  z^i_j = \sum_{g_i(c)=j} x'_c. \]

That is, we “hash” all cells into $|B_i| \ge q$ buckets, and sum all cells
hashed to the same bucket. The measurement vector $z=Ax'$ is now equal
to a concatenation of vectors $z^1, \dotsc, z^s$.  Note that the
dimension of $z$ is equal to $m= w'^2 \sum \abs{B_i} = O(qs) = O(k
\log_k N)$.

§.§ Recovery algorithm

A graphical representation of the recovery process is presented on the second page of Appendix <ref>.
The recovery algorithm starts by identifying the buckets that likely contain the cells from $S$, and labels them consistently (i.e., two buckets containing cells from $S$ should receive the same label), allowing for a small fraction of errors. We then use the labels to identify the cells.

The algorithm runs as follows. For a set $X \subset [s] \times [2q]$ of
pairs of indices, let $F(X)$ denote $\{F(z^i_j) : (i,j) \in X \}$.

* Identify $R = \{(i,j) : \|F(z^i_j)\|_\Gamma \ge T/2\}$ (that is, $R$ contains the “heavy cells” of the measurement vector $z$).

* Partition $R$ into sets $R',R^1,\dots,R^k$ such that $|R'| \le \delta sk$, and such that for each $1 \le l \le k$ the diameter of $F(R^l)$ is at most $T/2$.

* For each label $l=1, \dotsc, k$, create a vector $u^l \in B$ such that for each $i=1, \dotsc, s$, $u^l_i=j$ if $(i,j) \in R^l$ (if there are many such $j$, ties are broken arbitrarily), or $u_i^l=\erasure$ (an arbitrary erasure symbol) if no such $j$ exists.

* For each label $l=1, \dotsc, k$ apply the decoding algorithm[Technically, we replace each $\erasure$ in $u^l$ with an arbitrary $j$ before running the decoding algorithm, since the decoding algorithms don't know about $\erasure$.] for $g$ to $u^l$, obtaining a (possibly invalid) decoded cell $d^l$.

We analyze the algorithm by keeping track of the errors at each step.

Step 1  For any cell $c \in S$ and $i =1, \dotsc, s$, we say that $i$ preserves $c$ if $\|F(z^i_{g_i(c)}) - F(x_{c}) \|_\Gamma \le T/24$ and $g_i(c') \ne g_i(c)$ for all other $c' \in S$. That is, there is no collision from the hashing process, and the total amount of distortion due to the noise $\mu$ is small. Let $P = \{(i,g_i(c)) : i \mbox{ preserves } c\}$. Note that $P\subset R$. We show that $P$ is large and that most of $R$ is in $P$.

With probability at least $7/8$,
\[ |P| \ge (1-\beta)sk. \]

Consider any pair $(c,i) \in S \times \{1, \dotsc, s\}$, and let $j=g_i(c)$. If $i$ does not preserve $c$, it must be because either (i) there is another cell $c' \in S'$, $c' \neq c$ such that $g_i(c')=j$, or because (ii) the total noise affecting $z^i_j$, equal to $F(\mu^i_j)\le \sum_{g_i(c)=j} F(\mu_c)$, has norm at least $T/24$.

By Lemma <ref> with probability at least $7/8$
the number of pairs affected by (i) is at most $32k s |S'|/q$.  The
event (ii) is determined by the noise vector $\mu$. However, for each
$i$, there are at most $\frac{\sum_c \norm{\Gamma}{F(\mu_c)}}{T/24} = \frac{\|\mu\|_F }{T/24} < 24 \gamma k$
additional cells $c \in S$ that are not preserved under $i$ due to
this reason, where the latter inequality follows from the assumption that $\norm{F}{\mu} < \gamma k T$ (Section <ref>).

Altogether, the total number of pairs $(c,i)$ such that $c$ is not preserved by $i$ is at most

\[ 32sk\abs{S'}/q + 24 \gamma sk \leq [32 \eta (1+16\alpha) +24\gamma]  s k =\beta sk \]

for some small constant $\beta$, as desired.

With probability at least $3/4$,
\[ |R \setminus P| \le \delta sk. \]

Any element $(i,j)$ of $R\setminus P$ (“heavy but not preserved”) must belong to one of the following three categories:

* $j=g_i(c)$ for $c \in S$ such that $c$ is not preserved by $i$. By the previous lemma, there are at most $\beta s k$ such pairs $(c,i)$ with probability at least $7/8$.

* $j=g_i(c)$ for some cell $c \in S'\setminus S$. There are at most $16\alpha s k $ such pairs $(c,i)$, with probability at least $7/8$.

* The vector $F(\mu^i_j)=\sum_{g_i(c)=j} F(\mu_c)$ has norm at least $T/2$. There are at most $2 \gamma s k $ such pairs $(i,j)$.

This implies that with probability at least $3/4$ the total number of pairs $(i,j) \in R\setminus P$ is at most

\[ (\beta + 16\alpha  + 2 \gamma) sk = \delta sk \]

for some small constant $\delta$, as desired.

Step 2  Observe that the elements of $F(P)$ can be
clustered into $k$ clusters of diameter $T/12$. Thus, by the previous
lemma, there is a $k$-clustering of all but $\delta sk$ elements of
$F(R)$ such that the diameter of each cluster is at most $T/12$. We
now apply a $6$-approximation algorithm for this problem, finding a
$k$-clustering of $F(R)$ such that the diameter of each cluster is at
most $T/2$. Such an approximation algorithm follows immediately from
the $3$-approximation algorithm for $k$-center with outliers
in [5], which gives a $k$-clustering with radius at most $T/4$
and hence diameter at most $T/2$.

Step 3  Consider cells $c,c' \in S$ such that  $c$ is preserved by $i$ and $c'$ is preserved by $i'$.
If $F(z^i_{g_i(c)})$ and  $F(z^{i'}_{g_{i'}(c')})$ belong to the same cluster, then it must be the case that $c=c'$, since otherwise the distance between them would be at least $T - 2T/24>T/2$. In other words, for each $l$, if $u^l \subset P \cap R^l$ contains at least one element of $P$, then all the elements of $u^l$ are “derived” from the same cell.

With probability at least 3/4, $u^1,\dotsc, u^k$ contain a total of at most $2\delta sk$ errors and $(\delta+\beta)sk$ erasures ($i,l$ such that $u_i^l = \erasure$).

Let $R'' = R \setminus R' = R^1\cup\cdots\cup R^k$. Let $P' = P \cap R'$, and $P'' = P \cap R''$.

Note that $|P'| \le |R'| \le \delta sk$.
Each error in $u^1, \dotsc, u^k$ corresponds to a unique element of $R'' \setminus P''$, and we have

\[ |R'' \setminus P''| \le |R'' \setminus P| + |P \setminus P''| \le |R \setminus P| + |P'| \le \delta sk + \delta sk = 2\delta sk. \]

Additionally, $\{u^1, \dotsc, u^k\}$ contains at least $P''$ elements total, and so the number of erasures is at most $sk - |P''| = sk - |P| + |P'| \le \beta sk + \delta sk$, where we use $|P| \ge (1-\beta) sk$ from Lemma <ref>.

Step 4  We can replace erasures by errors, and conclude that $u^1, \dotsc, u^k$ have a total of at most $(3\delta+\beta)sk$ errors. It follows that at least $k/2$ of them have at most $2(3\delta+\beta)s$ errors, and therefore can be decoded. Therefore, the set $D=\{d^1, \dotsc, d^k\}$ contains at least $k/2$ elements of $S$.

The running time of the recovery algorithm is dominated by Step 2, where we approximate $k$-median with outliers via the method in [5]. This takes $O((ks)^3) = O(k^3 \log^3_k n)$ time.

Assume $k \ge C\log N$ for some constant $C$, a signal $x$ with $k$
objects, and a noise vector $\mu$, all subject to the constraints
delineated in the Model description of
Section <ref>.  There is a distribution over random
binary $m \times N$ matrices $A$, $m=O(k \log_k N)$, and an
associated recovery algorithm with the following property.  Suppose
that the algorithm is given $Ax'$ for $x'=x+\mu$. Then the algorithm
recovers (with probability at least $3/4$) a set $D$ of $k$ cells,
such that at least $k/2$ of the cells fully containing an object are
included in $D$.  Moreover, the algorithm runs in $O(k^3 \log^3_k
n)$ time and the matrix has column sparsity $O(\log_k

§ APPLICATIONS TO ATTITUDE DETERMINATION

Star trackers determine their attitude, or 3-axis orientation, by taking pictures of the sky and identifying stars in the image. We provide a detailed review of the current technology in the first section below. We then present a compressive sensing algorithm for attitude determination, along with a discussion about hardware implementation. Finally, we present results from a software simulation of the algorithm.

§.§ Current star trackers

A star tracker is essentially a digital camera, called a star camera, connected to a microprocessor. We describe various characteristics of the camera hardware and star identification algorithms.

§.§.§ Numbers

We first provide some numbers from [22] and [36] to give a sense of scale. As of 2001, a typical CCD star tracker consumes 5-15W of power[22]. A small spacecraft uses 200W of power, and a minimal one uses less than 100WciteRL, so this can be a substantial amount.

A high-end star tracker can resolve approximately the same set of stars that an unaided human can on a moonless night away from all light pollution. The number of stars in a star tracker's database varies from 58 to many thousands. The camera's field of view can vary from $2\times 2$ degrees to $30\times 30$ degrees, or anywhere from .01% to 4% of the sky. For comparison, the full moon is about .5 degrees across, and an adult fist held at arm's length is about 10 degrees across. A CCD camera can have up to a million pixels, and the accuracy of the final attitude is usually around .001 degrees (1 standard deviation), compared to .01 degrees for the next best sensors. The attitude is updated anywhere from 0.5 to 10 times a second.

§.§.§ CCD and CMOS

CCD (charge-coupled device), and CMOS/APS (complimentary metal-oxide semiconductor/active pixel sensor) are the two different types of sensors used in star cameras. We abuse notation and use CCD/CMOS to refer to the sensor, an $n\times n$ array of sensors, the architecture of this array, and the star cameras they are a part of. Both CCD and CMOS turn incoming photons into electrons using the photoelectric effect, read the electrons as a single voltage or charge, and then digitize the charge with an analog to digital converter (ADC).

The CCD has a single ADC located in the corner of the pixel array. A CCD array is read as follows: each pixel repeatedly transfers its charge into a neighboring pixel, so that the charge from any given pixel eventually travels a taxicab path to the ADC. Charges from different pixels are never combined, so there are a total of $\Theta(n^3)$ charge transfers. Since the ADC only digitizes one pixel at a time, it also takes $\Theta(n^2)$ time to read the whole array. In addition, each charge transfer leaves a fraction $\epsilon$ of the electrons behind, where $1-\epsilon$ equals the charge transfer efficiency. The electrons in the farthest pixels undergo $\Theta(n)$ charge transfers, and in practice it is costly to achieve $\epsilon < 10^{-5}$, which puts a bound on the maximum size of a CCD array [19, 12]. Even if future technology were to allow a better charge transfer efficiency, it is worth noting that each charge transfer uses a large constant amount of power, and that the total number of charge transfers is super-linear in the number of pixels.

On the other hand, CMOS devices have an ADC built into every pixel. This solves all of the problems noted above, and adds another important feature: random access reading. In other words, we can choose to read and digitize only a subset of the pixels, and in practice, that is done, saving power and subsequent digital processing costs [22]. However, the ADCs take up valuable real estate, and reduce the percentage of the chip that is available to collect photons. CMOS devices also generate substantially more noise than CCDs, further reducing the signal to noise ratio [23].

In practice, many consumer products such as cell phone cameras use CMOS, while scientific instruments use CCDs. Nevertheless, star trackers on small or low-power-budget satellites are starting to use CMOS, forgoing factors of 10 and higher in precision. We give the specification of a CCD tracker and a CMOS tracker in current (2011) production to illustrate the difference, as well as some of the numbers in highlighted in Section <ref>. The CT-602 Star Tracker has a CCD camera and is made by Ball Aerospace & Technologies. It uses 8-9W of power, weighs 5.5kg, has 6000 stars in its database, an $8\times 8$ degree field of view, $512\times 512$ pixels, an attitude accuracy of .0008 degrees, and updates 10 times a second [1]. Comtech AeroAstro's Miniature Star Tracker has a CMOS camera, uses $<2$W of power, weighs .4-.8 kg, has 1200 stars in its database, a $24\times 30$ degree field of view, and $1024\times 1280$ pixels, but has an accuracy of only .03 degrees and updates only 2 times a second [7].

§.§ Star tracker operation

We now study the process of attitude determination in more detail. As noted in Section <ref>, we often have an approximate attitude and go directly from Step 1 to Step 4 below.

0. Celestial intuitions
Let the apparent mass (hereafter, mass) of a star be the number of photons from the star that reach our camera. The masses of the stars in the night sky follow a power-law distribution, with exponent $-1.17$ [22].

The stars are essentially infinitely far away, and can be treated as point sources of light. In particular, if a star camera were perfectly focused, the light from any given star would land on exactly 1 pixel, and we would not be able to resolve the centroid of any star at higher than pixel resolution. Star camera lenses are therefore intentionally out of focus, so that the light from a star lands on multiple pixels, which we can then average to get a centroid with sub-pixel resolution. The blurring can be mimicked by convolving the image with a Gaussian of radius .5 pixels[Technically, we want to use an Airy function, not a Gaussian. But for a blur of radius .5 pixels a Gaussian is a good approximation.] [22], after which most of a star's mass lands on 4-6 pixels. Note that this number is independent of the star's mass, the field of view, or the total number of pixels.

Additionally, all the stars that are visible to star trackers (or humans) are a part of the Milky Way galaxy, which is shaped like a disk. This means that the stars are not distributed uniformly over the sky. The density of stars varies by a factor of 4 [24]; Figure <ref> in Section <ref> gives some intuition. Most star trackers with a small field of view are not able to resolve the area perpendicular to the galaxy.

1. Image acquisition
The expected number of photons (i.e. mass) captured by a pixel is proportional to the area of the lens and the exposure time. Photons hit a pixel, and eject an electron with some constant probability via the photoelectric effect. The pixel is then read and digitized (Section <ref>).

Almost all the noise is introduced at this step. The signal-dependent portion of the noise, shot noise, is due to random photon arrival times. This is best modelled by a Poisson distribution, which can be approximated as a Gaussian with a standard deviation equal to the square root of the expectation [19]. For example, if a pixel expects to see 900 photons, it will see 900 photons with a standard deviation of $\sqrt{900} = 30$ photons.

There is also a large per-pixel component of the noise, which depends on the hardware, temperature, and other factors. See [19] for a summary.

2. Centroiding
We locate a set of stars $S$ in the image, by finding pixels with a value several standard deviations above the mean. We then centroid each star, by either taking the weighted mean of a few neighboring pixels, or by doing a more complicated Gaussian fitting. For a bright star, this can give a centroid resolution of .1 pixels in each dimension, even under moderate noise.

There are several errors that can occur here. If the threshold for being a star is low, we may get stars that are composed only of noise. If our catalog is incomplete (as it almost certainly will be), we may get stars in $S$ that aren't in the catalog. Finally, radiation or dust can occasionally cause a pixel to have an arbitrary value, in which case it might get classified as a star.

3. Star identification
In this step, we match the stars in $S$ to an onboard database. If $S$ has zero or one stars, there is nothing to do, and we give up.
Some star trackers may be able to make an identification when $|S|=2$, but we assume that $|S|\ge 3$. We explain the algorithm from one common star identification algorithm [27], since almost all of them follow the same outline [32].

Preprocessing  Depending on the size of the lens, the capabilities of the processor, and the expected signal to noise ratio, anywhere from a few dozen to tens of thousands of stars are selected from a catalog. The star tracker firmware is loaded with a data structure $D$ that stores all pairs of stars that can appear in the same picture, along with their distances.

In space  We check subsets of three stars from $S$ to see if they can form a valid triangle using edges from $D$. Once we find one, (and if $|S|\ge 4$), we try to find a fourth star from $S$ such that all $\binom{4}{2}$ distances are consistent with $D$, in which case we declare a match. We then use the 3-4 matched stars to determine an approximate attitude.

Even with tens of thousands of stars, it is extremely unlikely that we find a match that is consistent with $D$ that ends up being wrong. In other words, most quadrilaterals consistent with $D$ are far away from all other such quadrilaterals, and most quadrilaterals formed using erroneous stars match nothing.
However, it is possible for there to be no noise in Steps 1-2 and still not have a match, due to catalog errors and omissions.

4. Precise attitude computation
We use the approximate attitude information and the onboard database to obtain the rough location of other stars in the picture. The attitude is then refined using the centroids of those stars as well.

Note that with a CMOS imager, if we start with an approximate attitude, we can actually skip most of Step 1 as well, and only read and digitize the handful of pixels under the stars we want to centroid.

§.§ Specialization of the general algorithm

We specialize the theoretical algorithm presented in Section <ref> to obtain a new algorithm for attitude determination, which we call Attitude Determination Under Analog Folding (ADUAF). The local geometric objects turn into stars, and the features we use are star centroid and mass. We use $f_{CRT}$ as our underlying error correcting code, rather than $f_{RS}$ (Section <ref>).

The biggest change from the theoretical algorithm to the practical algorithm is that we assume the stars are randomly rather than adversarially distributed. This means we no longer need to compose our error correcting code with the hash family $\Hc_p$. Additionally, we no longer need to shift the grid $G$ by the random vector $v$, as we do in Section <ref>. In fact, the notion of a grid is no longer needed at all, and we allow the decoded cell $d^l$ (Section <ref>) to be any $w' \times w'$ patch of the original image.

Given the physical cost of splitting the signal, the column sparsity $s$ is set to the smallest value our algorithm can tolerate, which is 2. This has the additional effect of turning the clustering process from Step 2 of Section <ref> into a much simpler bipartite matching process.

Finally, rather than just recovering the $w' \times w'$ patch $d^l$, we recover the centroid of $d^l$, or the centroid of the star on $d^l$. We then run the recovered centroids through the star identification process of Section <ref>.
The full algorithm is presented below.

§.§.§ Measurements

If $p_1$ and $p_2$ are the primes used to construct the CRT code in Section <ref>, we find primes $p'_i$ such that $p'_i \approx \sqrt{p_i}$. Note that we don't literally need $p'_1$ and $p'_2$ to be prime, as long as they are relatively prime, since the Chinese Remainder Theorem from Step 4 of the recovery algorithm applies for relatively prime numbers as well. We will use the word prime to refer to $p'_1$ and $p'_2$ even when they are just relatively prime.

Thinking of the noisy incoming signal $x'$ as an $n$-by-$n$ image, we define $z^i$ to be a $p'_i$-by-$p'_i$ image with

\[ z^i[j_1,j_2] = \sum_{\substack{c_1 \equiv j_1 \pmod{p'_i} \\ c_2 \equiv j_2 \pmod{p'_i}}} x'[c_1,c_2]. \]

We say we fold $x'$ to obtain $z^i$. Figure <ref> provides some intuition.

Taking measurements with column sparsity 2. Picture idea from [35].

We actually think of each $z^i$ as a torus, and identify the top/bottom and left/right edges together, allowing recovery of stars that are on the edge of the folded picture. Since the stars were assumed to be randomly distributed in $x'$, they are randomly distributed within each $z^i$ as well. Also, one could define the “measurement vector” $z$ to be a 1-dimensional representation of the pair $(z^1,z^2)$, and construct a measurement matrix $A$ accordingly. However, it is more useful to think of each $z^i$ as a 2-dimensional image.

§.§.§ Recovery algorithm

We follow the presentation of the corresponding paragraph in Section <ref>. For concreteness, we use a few specific numbers from our software implementation in Section <ref>.

* In each $z^i$, we identify ten $3\times 3$ ($w\times w$) cells with high mass such that no two cells collide in more than four pixels. We allow the $3\times 3$ cells to wrap around the edges of $z^i$, as in a torus.

* In place of the $k$-centering algorithm, we greedily choose up to eight pairs of cells $(c^1,c^2)$ from $(z^1,z^2)$ such that the feature vectors $F(c^1)$ and $F(c^2)$ are close. In our case, $F(c)$ is a triple of values: the mass of $c$, along with two coordinates for the centroid of $c$.

* No longer needed. Each pair from Step 2 corresponds to one of the $u^l$.

* In place of the error correcting code, we simply apply the Chinese Remainder Theorem in each dimension to recover a $3\times 3$ region of the original image for each pair from Step 2.

We compute centroids of the recovered $3\times 3$ regions, by averaging the centroids (weighted by mass) of the corresponding $3\times 3$ regions of $z^1$ and $z^2$. This is now $S$ in Step 3 of Section <ref>.

Step 1 above takes time linear in the number of measurements, or time $p_1+p_2\approx n$. Steps 2-4 take time at most quadratic in the number of cells chosen in Step 1 (ten, in our example).

§.§.§ Potential hardware

We have identified several existing hardware architectures that implement the folding mechanism above. In our case, we would use $s$ copies ($s=2$ above) of the folding hardware, where $s$ is the column sparsity. We call each of the stacked squares in Figure <ref> a piece of the image. We denote the number of pieces by $\tau$ ($\tau=9$ in both $z^1$ and $z^2$ in Figure <ref>). After being folded, each piece of the image lands on the focal plane, which has an array of $z^i$ pixels.

The most plausible architecture uses mirrors to directly reflect the pieces of the image onto the focal plane. The single-pixel camera [9] implements a generalized version of this for $z^i=1$ using $n$ moving mirrors. We would use $\tau$ rigid mirrors for each $z^i$, which is a substantial simplification of what they implemented.

Another possibility is to use $\tau$ different lenses, each of which focuses one piece of the image onto the focal plane [35]. The challenge with this design will be making all the lenses face the same direction.

The remaining two architectures are designed for unrelated earth-based tasks, and are more wasteful of the signal than is appropriate for a star tracker. Nevertheless, they show the extent to which researchers are already thinking about hardware designs for folding.
Both of the designs use beamsplitters and mirrors to combine the signal, and lose a factor of $\tau$ from the signal in the process. Figure <ref> depicts an element from each design. To extend <ref>(a) to $\tau$ pieces one would use $\log{\tau}$ different such image combiners, and to extend <ref>(b) to $\tau$ pieces one would use $\sqrt{\tau}-1$ stacked beamsplitters plus a mirror to fold in one direction, and $\sqrt{\tau}-1$ stacked beamsplitters plus a mirror to fold in the other.

[From [35].]
[From [34].]
Two examples of how to combine pieces of an image using a mirror and beamsplitter. Dashed lines indicate lost signal.

Finally, it is possible we could split the signal and fold it after it hits the focal plane but before it gets digitized, which would not save on sensor costs, but would reduce the number of ADC computations, and reduce the load on the onboard processors. CCD images query an entire row of the image before digitizing [23], so folding in at least the first dimension could be relatively easy. CMOS imagers are already built directly onto integrated circuits, so additional circuitry for folding is likely to be cheap [30]. Several (non-satellite) CMOS imagers have already been built that use dense compressive sensing matrices to reduce the number of measurements that need to be digitized [30, 26].

§.§ Software implementation

We implement the algorithm presented in Section <ref>, and run experiments on simulated images. Source code is available at <http://web.mit.edu/rishig/papers/local-geo/>.

For simplicity, we project the stars onto a rectangular interval, rather than operating on a sphere (Figure <ref>). We index the rectangle by right ascension (longitude) $\alpha \in [-\pi,\pi]$ and declination (latitude) $\delta \in [-\pi/2,\pi/2]$; for example, $\delta = 0$ is the equator, and $\delta = \pi/2$ is the north pole. So that the approximation makes sense, we ignore the portion of the sky where $|\delta| > \pi/2 - \pi/8$ (the dashed blue lines in Figure <ref>). This also has the effect of removing the portion of the sky that has the fewest stars, which some star trackers don't operate on anyway. We assume without loss of generality that the camera is axis-aligned.

Mercator projection of the night sky. The dense omega-shaped region is the central disk of the galaxy. We test our algorithm on the area between the dashed blue lines.

We fix $n = 800$ ($N = 640000$) for all experiments. We expose the camera to a .08 radian by .08 radian (4.6 by 4.6 degree) patch of the sky, which means that a typical image will have 3 to 10 bright stars (10%-ile to 90%-ile) and 50 to 150 stars that act as background noise. Above, we defined the mass (apparent mass) of a star to be the number of photons from the star that hit our camera. In our pictures, if the median total star mass is scaled to 1, the 10%-ile mass is .6, and the 90%-ile mass is 2.25. These numbers were determined empirically from the Smithsonian Astrophysical Observatory (SAO) Star Catalog [33]. Recall the mass of the $j^{th}$ brightest star in the sky is $\Theta(j^{-1.17})$ [22].

We choose the stars for the preprocessed database $D$ (Step 3 of Section <ref>) as follows. We extract a subset SAO$'$ from the full SAO Catalog, by taking the 10 most massive stars in every ball of radius .08 radians. SAO$'$ has 17100 stars, compared to 259000 stars in the SAO catalog.

To generate test images, we randomly select axis-aligned patches of the sky from the area between the blue lines of Figure <ref>, and use the SAO star catalog to simulate star positions and mass. We convolve the stars with a Gaussian of radius $.5$ pixels, which causes the vast majority of a star's mass to fall within a $3\times 3$ box. We then apply Poisson noise to account for photon arrival times, and fold the image using two relatively prime numbers. We add Gaussian noise (amount varying by experiment) to each pixel of both folded images to account for all signal-independent sources of noise. To keep the recovery process simple, we do not account for the Poisson noise introduced by splitting the signal; in other words, we apply Poisson noise once to each star, whereas in a hardware implementation it is likely that Poisson noise will be applied independently to each folded copy of a star. =-1

Figure <ref> has a few pictures to help build intuition about star density and what stars look like.

It is generally possible to see where the biggest stars are located, though some fraction of them get occluded by small stars or noise.

[Underlying noiseless signal $x$ (zoom in).]

[$p'= 29$, No Noise.]
[$p' = 29$, Noise = 150]
log(mass) in sample images from a representative part of the sky. The legend on the right applies to all three images.  We cut off pixels with values below 0 before taking the log in (c).

We run our algorithm, ADUAF, as well as Sequential Sparse Matching Pursuit (SSMP) on 159 images of the sky. SSMP includes a linear (in $N$) time sparse binary compressive sensing recovery algorithm, with recovery quality on par with other known methods [3]. We use the same folded images, or equivalently, the same measurement matrix, for both algorithms.

ADUAF recovers a list of centroids, and we find and centroid stars in the image recovered from SSMP. We run the star identification algorithm from Step 3 of Section <ref> on the outputs of both of them, and declare success if they identify the stars correctly. We report our results as a function of the standard deviation of the added Gaussian noise (Figure <ref>).

Experimental results. Each point on the figure is computed using the same 159 underlying images.

The first observation we make is that ADUAF works very well down to an almost minimal number of measurements. The product $p'_1p'_2$ has to be greater than 800, and the minimal set of primes is 26 and 31. As the number of measurements increases, SSMP catches up and surpasses ADUAF, but we note that running SSMP (implemented in C) takes 2.4 seconds per trial on a 2.3 GHz laptop, while ADUAF (implemented in Octave/Matlab) takes .03 seconds per trial. Computation power on a satellite is substantially lower than that of a low end laptop, and given that the entire acquisition has to happen in .1 to 2 seconds, it seems unlikely that any algorithm linear or near linear in $N$ is going to be practical. Finally, we note that the plot lines for both SSMP and ADUAF could be improved by a more sophisticated implementation of the star identification algorithm.

§ ACKNOWLEDGEMENTS

The authors would like to thank Tye Brady and Ben Lane from Draper Laboratory for numerous conversations and for help with the data. We would also like to thank the anonymous reviewers for their thorough reviews and for helping clarify the presentation.

This research has been supported in part by a David and
Lucille Packard Fellowship, MADALGO (Center for Massive
Data Algorithmics, funded by the Danish National Research
Association) and NSF grant CCF-0728645. R. Gupta has been supported in part by a Draper Laboratory Fellowship. E. Price has been supported in part by an NSF Graduate Research Fellowship.

[1]
Ball Aerospace & Technologies Corp.
CT-602 star tracker.
Linked from <http://www.ballaerospace.com/page.jsp?page=104>.

[2]
R. Berinde, A.C. Gilbert, P. Indyk, H. Karloff, and M.J. Strauss.
Combining geometry and combinatorics: A unified approach to sparse
signal recovery.
In 46th Annual Allerton Conf. on Communication, Control, and
Computing, pages 798–805. IEEE, 2008.

[3]
R. Berinde and P. Indyk.
Sequential sparse matching pursuit.
In 47th Annual Allerton Conf. on Communication, Control, and
Computing, pages 36–43. IEEE, 2009.

[4]
M. Charikar, K. Chen, and M. Farach-Colton.
Finding frequent items in data streams.
Theoretical Comp. Sci., 312(1):3–15, 2004.

[5]
M. Charikar, S. Khuller, D.M. Mount, and G. Narasimhan.
Algorithms for facility location problems with outliers.
In Proc. 12th Annual Symp. on Discrete Algorithms, pages
642–651. SIAM, 2001.

[6]
G. Cormode and S. Muthukrishnan.
An improved data stream summary: the count-min sketch and its
J. Algorithms, 55(1):58–75, 2005.

[7]
Comtech AeroAstro.
Miniature star tracker data sheet.
Linked from

[8]
E. J. Candès, J. Romberg, and T. Tao.
Stable signal recovery from incomplete and inaccurate measurements.
Comm. on Pure and Applied Math., 59(8):1208–1223, 2006.

[9]
M.F. Duarte, M.A. Davenport, D. Takhar, J.N. Laska, T. Sun, K.F. Kelly, and
R.G. Baraniuk.
Single-pixel imaging via compressive sampling.
Signal Processing Magazine, IEEE, 25(2):83–91, 2008.

[10]
K. Do Ba, P. Indyk, E. Price, and D.P. Woodruff.
Lower bounds for sparse recovery.
In Proc. 21st Annual Symp. on Discrete Algorithms, pages
1190–1197. SIAM, 2010.

[11]
D. L. Donoho.
Compressed Sensing.
IEEE Trans. Info. Theory, 52(4):1289–1306, 2006.

[12]
E. R. Fossum.
Active Pixel Sensors: Are CCD's dinosaurs?
SPIE, 1900:2–14, 1993.

[13]
S. Foucart, A. Pajor, H. Rauhut, and T. Ullrich.
The gelfand widths of lp-balls for $0< p <= 1$.
J. Complexity, 26(6):629–640, 2010.

[14]
R. Fergus, A. Torralba, and W. T. Freeman.
Random lens imaging.
MIT CSAIL Technical Report 2006-058, 2006.

[15]
A. Gilbert and P. Indyk.
Sparse recovery using sparse matrices.
Proc. IEEE, 98(6):937–947, 2010.

[16]
O. Goldreich, D. Ron, and M. Sudan.
Chinese remaindering with errors.
IEEE Trans. Information Theory, 46(4):1330–1338, 2000.

[17]
V. Guruswami.
Introduction to coding theory.
Course notes, Lecture 1, 2010.
Available at
<http://www.cs.cmu.edu/ venkatg/teaching/codingtheory/notes/notes1.pdf>.

[18]
G. C. Holst and T. S. Lomheim.
CMOS/CCD Sensors and Camera Systems.
JCD Publishing and SPIE Press, 2007.

[19]
G. C. Holst.
CCD Arrays, Cameras, and Displays, pages 79–82, 91–93,
JCD Publishing and SPIE Optical Engineering Press, second edition,

[20]
P. Indyk.
Sketching, streaming and sublinear-space algorithms.
Graduate course notes, 2007.
Available at <http://stellar.mit.edu/S/course/6/fa07/6.895/>.

[21]
J. Justesen.
On the complexity of decoding Reed-Solomon codes (Corresp.).
IEEE Trans. Information Theory, 22(2):237–238, 1976.

[22]
C. C. Liebe.
Accuracy performance of star trackers — a tutorial.
IEEE Trans. Aerospace and Electronic Systems, 38(2):587–599,

[23]
Dave Litwiller.
CMOS vs. CCD: Facts and fiction.
Photonics Spectra, January 2001.

[24]
S. Lee, G. G. Ortiz, and J. W. Alexander.
Star tracker-based acquisition, tracking, and pointing technology for
deep-space optical communications.
Interplanetary Network Progress Report, (42-161), 2005.

[25]
James D. Meindl.
Beyond moore's law: The interconnect era.
Computing in Science and Engineering, 5:20–24, 2003.

[26]
V Majidzadeh, L Jacques, A Schmid, P Vandergheynst, and Y Leblebici.
A (256x256) Pixel 76.7mW CMOS Imager/Compressor Based
on Real-Time In-Pixel Compressive Sensing.
In IEEE Int. Symp. Circuits and Systems, 2010.

[27]
D. Mortari.
Search-less algorithm for star pattern recognition.
J. Astronautical Sciences, 45(2):179–194, 1997.

[28]
S. Muthukrishnan.
Data streams: Algorithms and applications.
Foundations and Trends in Theoretical Comp. Sci., 2005.

[29]
M. Nachin.
Lower bounds on the column sparsity of sparse recovery matrices.
MIT Undergraduate Thesis, 2010.

[30]
R. Robucci, J.D. Gray, L.K. Chiu, J. Romberg, and P. Hasler.
Compressive sensing on a CMOS separable-transform image sensor.
Proc. IEEE, 98(6):1089–1101, 2010.

[31]
J. Romberg.
Compressive sampling by random convolution.
SIIMS, 2009.

[32]
B. B. Spratling IV and D. Mortari.
A survey on star identification algorithms.
Algorithms, 2:93–107, 2009.

[33]
Smithsonian astrophysical observatory star catalog.
Available at

[34]
V. Treeaporn, A. Ashok, and M. A. Neifeld.
Increased field of view through optical multiplexing.
Optics Express, 18(21), 2010.

[35]
S. Uttam, A. Goodman, M. Neifeld, C. Kim, R. John, J. Kim, and D. Brady.
Optically multiplexed imaging with superposition space tracking.
Optics Express, 17(3), 2009.

[36]
J. R. Wertz and W. J. Larson.
Space mission analysis and design, pages 315–316, 321–322,
Space technology series. Microcosm press and Kluwer academic
publishers, third edition, 1999.