Detecting small low emission radiating sources

The article addresses the possibility of robust detection of geometrically small, low emission sources on a significantly stronger background. This problem is important for homeland security. A technique of detecting such sources using Compton type cameras is developed, which is shown on numerical examples to have high sensitivity and specificity and also allows to assign confidence probabilities of the detection. 2D case is considered in detail.


Introduction
One of the missions of the Department of Homeland Security is to prevent smuggling of weapon-grade nuclear materials (e.g., [25]). It is expected that such materials, unlike those needed for a "dirty bomb", will have low emission rates and will be well shielded, so that very few gamma photons or neutrons would escape, and even less would be detected non-scattered (ballistic). An additional hurdle for the detection of illicit nuclear substances is the strong natural radiation background. If the particles radiated by the source were physically distinct (e.g., in terms of their energies) from the ones prevalent in the background, then the detection would be easier. We thus assume that the particles from the source are identical to those coming from the background in terms of their nature, energies, etc., so the discrimination using such parameters is impossible. This leads to the situation when the signal to noise ratio (SNR), i.e. the count of the detected non-scattered particles emitted by the source versus the total number of detected particles, can be as low as 10 −3 , if not lower. It is rather clear that if only the locations of the hits are detected, with no information about the directions of the incoming particles available, there is no chance to detect the presence of such a weak source. Thus, one needs to employ detectors that can provide some directional information. For instance, the γ-cameras most commonly used in SPECT (Single Photon Emission Computed Tomography) medical imaging [4,27], are mechanically collimated, and thus they "count" the particles coming from a narrow cone of directions and discard the rest. However, mechanical collimation dramatically reduces the particle count, and thus, in the case of an extremely low emission source, can essentially eliminate the useful signal (this difficulty, although in much less extreme form, arises also in SPECT). Typically, only one out of thousands emitted photons is detected by a collimated camera. When dealing with low emission sources, one wants to capture as many particles as possible, and thus mechanical collimation is unsuitable for the problem at hand.
Another option is to use the so called Compton γ-cameras (see Section 5), which do not discard any incoming particles by collimation. The price one pays for this is that the directional information provided by such a camera is more limited. Namely, one obtains just a hollow (i.e., surface rather than solid) cone of possible directions rather than a single direction. Still, Compton cameras are good candidates for the applications we have in mind. While Compton cameras are available for detecting γ-photons, neutron detectors that provide a similar "incoming cone" information are being currently developed as well [25,40]. The methods we will discuss do not depend upon a specific kind of particles; we will thus call all detectors that can provide the cone information Compton type detectors.
Even if some directional information is available, the task of detecting extremely low emission coming from a geometrically large source in the presence of a significant radiation background would be very hard, if not impossible. So, another important condition that we impose, besides availability of directional information, is the small size of the possible source, which is usually a safe assumption in the applications we are interested in.
Simulations of radiation from a small high-enriched uranium source placed in a cargo container suggest that only about 0.1% of the signal received by detectors might be due to the ballistic (non-scattered) particles emitted from the source. The remaining 99.9% of detected particles come either from extraneous (natural) sources, or from scattered source particles from the source [5].
It is natural to try to use the rather standard SPECT techniques in the present situation, albeit the chances of success (and even the applicability of the tomographic models) are questionable. This circle of issues is discussed analytically and tested numerically in Section 2, where the conclusion is made that backprojection technique might be the best bet here, while the usual filtration parts in SPECT reconstructions do not do any good. This conclusion gets its foundation in the probabilistic discussions of Section 3 and then in numerical examples of Section 4. However, to check the viability of the approach, all these considerations dealt with the simplest, collimated detectors. Section 5 is devoted to the case of Compton type cameras. Here tomographic and backprojection techniques are considered. The analysis of the 2D Compton camera case is provided in Section 6. The corresponding numerical tests are conducted and discussed in Section 7. The overall conclusion is that the suggested backprojection Compton type cameras techniques allow a robust detection of presence of geometrically small low emission sources with extremely low SNR, with confidence levels attached.
Like in the case of SPECT [4], one can try to use statistical methods. E.g., the recent paper [46] describes a Bayesian approach to the same problem.
The paper ends with the sections devoted to remarks and conclusions, acknowledgments, and bibliography.

Source detection using tomographic techniques
We will briefly present here some relevant mathematical formulas and conclusions from the SPECT version of emission tomography, which is the closest to the problem of our interest. One can find more details in [4,21,22,27,28].
Let f (x) be the intensity distribution of the sources of particles of certain type (γ-photons, neutrons, etc.) inside an object (e.g., cargo container, or a truck). Let also µ(x) be the attenuation coefficient (usually not known the applications we have in mind). Then, in the approximation of sufficiently high emission rate, the radiation transport equation implies that the particle count per unit of time at a detector collimated in the direction L is Here L x is the segment of the line L between the emission point x and the detector and dy denotes the standard linear measure on L. The operator T µ is said to be the attenuated Radon (or X-ray) transform (with attenuation µ(x)) of the function f (x).

Figure 1: The set-up of the SPECT emission tomography
We will fix an origin O and assign the normal coordinates (ω, s) to any line L. Here ω = (cos θ, sin θ) is a unit vector orthogonal to L and s is the (signed) distance from the origin to L (see Fig. 1). The equation of this line is x · ω = s and the line itself will be sometimes denoted as L (ω,s) . The value of (1) provides the mathematical expectation of the counts at a detector per unit of time, which is not a good approximation of the measured data, if the emission rate is low. Let us still accept for the time being this integral-geometric formulation and see where it leads us.
The attenuated backprojection operator with attenuation ν(x, ω) > 0 (which may or may not be related to µ(x)) acting on data g(ω, s) is In particular, if ν = 1 (which we will assume in our experiments), this is the standard backprojection operator [20,27,28].
Let H be the Hilbert transform acting in variable s on functions g(ω, s): Here v.p. indicates that the integral is considered in the principal value meaning.
The following result holds: (see details in [22]) Let f (x) be a function (our distribution of radiating sources) and µ be sufficiently smooth 1 , then the operator applied to f preserves all singularities of f and does not create new ones. In particular, if f has a jump across an interface, the location of this jump is preserved in the image is not a reconstruction of f from its measured projections T µ f (which would be impossible anyway with the low signal strength and the attenuation µ being unknown), but it gives a correct reconstruction of all singularities (e.g., interfaces) of f .
Since we assume that our source is very small geometrically, it is rather singular, and thus one might expect that any guess for the unknown attenuation (e.g., that there is none) should reconstruct the location of the source correctly.
The operator is exactly what is applied to the data in our experiments described below. We thus will not worry about the presence of an unknown attenuation. The operator (4) is often called a filtered backprojection (FBP), where H d ds and T # 1 are the filtration and backprojection parts, correspondingly. Remark 2.
• Since we are looking for a rather singular object, it is reasonable to recall that there is a way in tomography to emphasize singularities of the function to be reconstructed. This is what the so called local tomography does (e.g., [10,22]). For instance, the singularities will look brighter if one replaces the filtration part H d ds in (4) with a stronger filter d 2 ds 2 : This option will also be explored.
• In the detection problem we are discussing, different types of particles (e.g., γphotons of different energies) will be present, and each of them could have different distribution f j of sources, as well as different attenuation distribution µ j in the cargo. Thus, what we get if we lump all these particle counts together (which seems to be a strange thing to do), is the sum 1 The detection procedures we will eventually come to, will not require smoothness of attenuation.
Can this complicate things? We claim that it should not. Indeed, all these functions have the same boundary of the source as their interface. One can derive now from the above theorem the following Corollary 3. Applying the same reconstruction operator T # ν H d ds to the whole sum, one recovers the correct source interface.
It is still not clear whether such a procedure will work in our problem, since there are at least two potentially serious obstacles: • The (attenuated) X-ray transform model does not apply when significant statistical noise is present, in particular for a very weak source.
• This model does not take into account presence of a very large (from standard tomographic point of view, enormous) background noise.
Still, we will try to apply this approach and see where it can lead us. We show below some sample results of numerical experimentation with the 2D and 3D X-ray transforms, as well as 3D Radon transform.

2D X-ray reconstructions
Although the problem is three-dimensional, for the first trial we consider the 2D problem, which is more straightforward to handle numerically.
In the first reconstruction we have modeled random γ-photon emission from a shielded ball of radius 4cm of HEU 235 for realistic values of the emission and attenuation rates. The ball was placed inside in the center 2 of a much (about 100 times) larger "cargo". A uniform and isotropic random background was also added. The model included 64 detectors, each of which could detect 64 directional sectors (about 2.8 degrees each). In this reconstruction we assumed the SNR to be 1%, i.e. 99% of detected hits were generated by background, or scattered source particles. Then, the described above tomographic X-ray transform inversion (see (3) with ν = 1) was applied to the resulting matrix of counts. The pictures below show the results of the reconstructions in 44 × 44 pixels (both the gray scale density and surface plots are presented). The result, shown in Fig. 2, shows clear detection of the source.
We then try to lower the SNR towards our benchmark value 0.1%, and the success seems to evaporate. E.g., Fig. 3 shows the failed attempt of reconstruction with the SNR being 0.4% and with only about a hundred ballistic particles detected from the source. Now let us consider a fixed signal-to-noise ratio (0.1%) and varying total number of detected particles. We will first describe the set up for the numerical results presented in rest of this subsection.
Four direction sensitive detector arrays, each consisting of 100 detectors (bins), were placed along the sides of the square [−1, 1] 2 . Random background radiation was created by simulating particles propagating along straight lines in such a way that these lines were uniformly randomly distributed in the square. More precisely, to each particle we assigned a direction of propagation by choosing the normal coordinates (θ, s) of a line such that θ and s were uniformly distributed in (−π/2, π/2) and [− √ 2, √ 2], correspondingly. Lines which pass through the imaged square necessarily intersect two  detector arrays, and one of the two intersection points was randomly chosen as the site of detection of the particle. Finally, for each such particle we recorded the bin on the detector array and the exact incoming direction α (see Fig 4). We also simulated an isotropic point source which emitted a number of particles roughly equal to 0.1% of the background. In order to use the filtered backprojection (FBP) inversion formula (4) we first computed the normal coordinates of the (straight-line) trajectories of the particles: where (x i , y i ) is the center of the bin on the detector array. After that, the coordinates θ and s were put into bins of size 0.02. A two dimensional data array was constructed in which every element equaled the number of particles with trajectories in the corresponding θ, s bin. Different tomographic inversions were applied to such data. The criterion for successful detection was whether the highest pick of the reconstructed image occurred in a small neighbourhood of the source location. Figure 4: A particle propagates along a random line with normal coordinates (θ, s) and collides with a detector array at the point (x, y). We record the bin with center (x i , y i ) and the direction α = θ + π/2.
In our first experiment, we simulate about 300, 000 background particles and additional 300 particles coming from a source located at (0.401, −0.133). In most cases, filtered backprojection reconstruction (4) would fail to detect the source: only 7 out of 20 experiments were successful. However, increasing the total number of detected particles, while keeping the SNR at 0.1%, improved the results significantly. Setting the background to 400, 000 and then to 500, 000 particles resulted in 13 and 18 detections out of 20 experiments, correspondingly.
As we have already mentioned, the local tomography (which uses a stronger filter, see (5)) emphasizes the singularities of an image, and thus one could expect it to detect geometrically small sources better. This, however, turns out not to be the case. In our experiments local tomography always showed inferior results to filtered backprojection, often failing to detect a source when an FBP detection from the same data was successful. A stronger high-pass filter just increases the contribution from the noise, so the reconstructed images, in fact, worsened. This suggests that maybe one should try to eliminate the filtration step altogether and relay upon the backprojection alone.
Indeed, our results using backprojection with no filtration were very encouraging. Backprojection of the same data used for FBP and local reconstructions resulted in 16, 18 and 20 out of 20 successful detections for 300, 000, 400, 000 and 500, 000 detected background particles, correspondingly. The results of these experiments are summarized in Table 2.1.   , only peaks deviating more than 3.2 standard deviations from the mean (right). A number of peaks in the image are higher than the peak at the source. Bottom row: backprojection reconstruction (left), only peaks deviating more than 3.5 standard deviations from the mean (right). The correct location of the source is found.

3D X-ray and Radon reconstructions
While in 2D there is a single transform that can be called interchangeably Radon transform or X-ray transform, these notions part in higher dimensions. In 3D, there is the X-ray transform, which integrates a function f (x) over all straight lines in R 3 , as well as the Radon transform, which integrates functions over all 2D planes. Both are of interest in tomography, the first arising in X-ray CT, SPECT, and PET scanners and the second in MRI imaging. In the situation we consider in this text, X-ray transform represents, the same way as it does in 2D, the case of collimated detectors, while the Radon transform will have some relation to the Compton camera case considered in Section 5. Indeed, in both 3D Radon and Compton cases, the data represents surface integrals (over planes and cones respectively).

3D X-ray reconstructions
The same issues that we have discussed in 2D case remain: the X-ray transform model is valid only for strong sources (or long observation time) and it does not address the background noise. However, in 3D there is another feature that one has to take into account. Namely, while in 2D the set of lines is also 2-dimensional (the lines can be parametrized by their two normal coordinates), in 3D the space of lines is 4dimensional. Thus, the X-ray data in 3D is redundant. In tomography, this redundancy is usually dealt with by selecting a 3-dimensional sub-set of data (the rest might not even been collected). For instance, one can use only rays parallel to a fixed plane and do the reconstruction slice-by-slice, or one can pick only lines intersecting a given curve (often a helix). In the detection situation we are facing, this is not an option, since any attempt to reduce dimension of the data will essentially remove all signal and leave only noise in its wake. Thus, we need to use all, overdetermined data. This requires corresponding analytic formulas (which are easy to obtain and thus will not be introduced) and much heavier calculations. Without further ado, we present an example (Fig. 6) of a successful FBP reconstruction, where 10 3 ballistic particles from the source and 1.35 × 10 6 background ones were detected.
Further experimenting with examples leads to the same conclusion as in 2D: backprojection is superior to FBP in this situation, while local reconstruction leads to deterioration.

3D Radon reconstructions
In this section, we present some Radon transform reconstructions in 3D. The reader should recall that the 3D Radon transform integrates functions over affine 2D planes, rather than lines. This is a preliminary test of the principle before doing 3D reconstructions for Compton cameras, where surface integrals (over cones) are also involved.
In the figures 7 -9 below, there were 10 6 background particles detected, while the number of ballistic particles from the source was varying, thus changing the SNR and the level of detectability. Density and surface plots of 2D sections through the source locations are shown. One sees that with the SNR of 0.5% the source can be detected, while at 0.1% and lower the detection deteriorates.
Let us see now whether backprojection alone does a better job. Figures 10 and 11 show the backprojection reconstruction from the same data as in Fig. 8 and 9. One  clearly sees that the backprojection alone detects the source when the FBP method fails to do so.
Again, the examples above confirm the general feasibility of the approach and superiority of the backprojection.

Discussion of the tomographic reconstructions
The experiments described in the previous sections beg several questions:  1. Why do tomographic methods still work to some degree, although the assumptions that lead to the X-ray transform model are not satisfied (i.e., the source is too weak) and a strong random background is present? Indeed, the use of standard tomography assumes that the data represents line integrals of the source distribution function. This is the same as to say that we measure the expectation of the number of hits, per unit time, from particles moving along a line. This assumption is reasonable when a substantial number of source particles is detected, i.e. when the source is strong, or the observation time is long. However, it fails in the case of low emission sources that cannot be observed for a really long time. Hence, the X-ray/Radon transform type models and techniques appear to be inappropriate.
2. Why does strengthening the filter (i.e., local tomography), which should increase detectability of singularities, makes the reconstruction worse?
3. Why do tomographic techniques fail to reliably detect sources at the levels of SNR and total number of particles for which the probabilistic consideration of Figure 10: BP reconstruction from 3D Radon data. 1000 ballistic source particles, SNR 0.1%. Top: Cross section through source. Bottom: Isosurfaces at 50% (blue), 60% (green) and 70% (red) of peak value.
the next section suggests that sources should be detectable?
As we might have guessed already, the answer to the first question is that in fact, tomographic inversions do not work. Or, to put it differently, only a part of the algorithm (backprojection) works, while its other part (filtration) only hurts.
The answer to the second question is not hard to guess. Indeed, high-pass filters required in X-ray/Radon transform inversions amplify noise. This is a very significant factor, as the noise constitutes 99.9% of our data. Using local tomography can only aggravate this difficulty, and we have seen that it indeed does.
The answer to the third question is probably the same as the previous one: highpass filters might be the culprits.
These considerations suggest to try to eliminate high-pass filtering completely, and thus to use backprojection alone (not a good idea in the standard imaging, where it leads to blurring [20,27]). We will see in the next sections that this is indeed the solution. Simultaneously, we will eliminate another drawback of the standard tomographic techniques, which is the difficulty of quantifying the confidence level of detection. Figure 11: BP reconstruction from 3D Radon data. 700 ballistic source particles, SNR 0.07%. Top: Cross section through source. Bottom: Isosurfaces at 55% (blue), 70% (green) and 80% (red) of peak value.

A probabilistic discussion
Now we investigate why it might be possible to detect geometrically small low emission sources by backprojection. Let us suppose that the detectors are collimated and thus can determine information about the incident direction of particles. Then backprojection essentially "sends" the detected particles back along their incoming straight lines. In particular, these backpropagated trajectories of all ballistic source particles will pass through the source region. The background particles or the source particles that have scattered before reaching the detector, will be backprojected along straight lines, which are different from their original zig-zag trajectories, about which we have no knowledge. Therefore, in the considerations that follow, we assume that the backpropagated trajectories of these particles are random straight lines. We also assume for now that the random distribution of trajectories is uniform.
In a nutshell, our simple argument, which will be made more precise below, is that if the number of lines passing through a tiny region exceeds significantly the mean of the background, then most probably there is a source at this location (see Fig. 12).
Let us try to quantify this and consider N random particle trajectories in a ball B R of radius R. What is the probability of at least n particles out of the total N contained Figure 12: If the local density of lines is significantly higher in the gray region, this is likely due to a source being present at this location. The large circle indicates the region of interest (ROI) rather than any physical boundary.
in B R to pass through a small ball B r with radius r? Let, as before, L (ω,s) denote the line in R 2 defined by its unit normal vector ω ∈ S 1 and a signed distance s from the origin, i.e.
Thus, the lines on the plane are mapped 1-to-1 to the points on the half-cylinder , and these points are uniformly distributed on the half-cylinder. All lines intersecting B R correspond to points on S 1 + × (−R, R), and those intersecting B r are bijectively mapped to points on S 1 + × (−r, r). Thus, the probability of a random line in B R passing through B r is The probability that n out of the total N lines cross B r is given by the binomial distribution B(N, p), and is equal to N n p n (1 − p) N −n . Since we are in the situation when p is fixed by the dimension of the source and N is large, the Central Limit Theorem applies. It follows (since the values of N p we will see will range in thousands) that the binomial distribution is approximated well by the Normal distribution N (µ, σ 2 ) with the mean µ = N p and standard deviation σ = N p(1 − p). If it happens that the number of lines n crossing B r exceeds significantly (in comparison with σ) the mean µ, the probability of this occurring just due to random reasons is very small. Therefore, one can be almost certain that this clustering of trajectories is the result of a radioactive source located at B r . Our numerical experiments of source detection by backprojection agree with these expectations.
As an example, in Fig. 13 we present a typical histogram of the number of lines crossing a pixel in the square D = [−1, 1] 2 with a uniform grid of 100 2 pixels. As described in section 2.1, we generated random background and an isotropic point source. The particles were detected by four arrays of collimated detectors placed on the sides of the square. The exact incoming directions of particles were recorded. For the result shown below, we backprojected each particle using a line-drawing algorithm. Namely, each particle was back propagated along a straight line into D by adding a value of 1 to all pixels in D intersected by the line. After all particles were backprojected, the value at each pixel of D is equal to the number of lines intersecting it. On figure 13 we show a histogram of the pixel values after bakprojecting particles coming from a random background of about 10 6 particles and 10 3 source particles. On the far right of the histogram one sees the outlier contribution from the source, which deviates about 10.3 standard deviations from the mean 3 . This result complies with discussion above: the appearance of such a highly unlikely outlier indicates presence of a source. Now we will make these statements more precise. Suppose that the area we are imaging is divided into N pix pixels. Suppose further that n lines intersect a pixel B r with radius r. We will write n = n s + n b , where n s and n b are the number of ballistic source and background trajectories, respectively, which cross B r . We will describe a method to determine whether there is a source located at this pixel.
Let us choose a threshold value k t , such that the probability of a normal variable to reach more than k t standard deviations above the mean is very small. E.g., k t = 5 or higher suffices. If an abnormally high number n of lines passing through a pixel B r is detected, i.e. if n > n t := µ + k t σ, we will claim that the pixel contains a source. Otherwise, such claim will not be made. The probability of at least n t lines crossing B r due to random reasons is approximately r := 0.5 erfc(k t / √ 2). We also have to take into account the total number N pix of pixels and the possibility of such clustering of lines occurring inside at least one of them. The probability of at least n t lines crossing at least one of the pixels due to random reasons (and thus causing a "false alarm") can be approximated by Here fp rate stands for false positive rate, which is the rate of false detections [11]. The above estimate assumed independence across different pixels of the events "the pixel is crossed by n lines". It is easy to see that this is in fact not the case: If a pixel is crossed by a large number of lines, there is a good chance that the surrounding pixels would also have a many lines intersecting. However, we claim that treating them as independent does not have a significant impact on the resulting estimates. To confirm this, we show some statistics from numerical simulations of random backgrounds at the end of this section. The true negative rate (or specificity), is given by tn rate = 1 − fp rate and represents the rate of correctly classified negative outcomes. The true negative rate is also the confidence probability that we have found a true source. We will write tn rate = confidence : In particular, if the threshold value k t is set to 5 and there are 100 2 pixels, the true negative rate is .997, which gives also this number as the confidence that the alleged sources are indeed true. This confidence probability depends only on the pre-set threshold value k t of the method and the number of pixels N pix . However, if in some pixel the number of lines deviates from the mean even further, this increases the confidence probability that there is indeed a source there. For instance, if n ≥ µ + 7σ, the confidence reaches, for all practical purposes, 100%. Let us discuss now the probability of missing a source. We assume that an apriori ballpark estimate of possible value of n s (the number of ballistic particles coming from the source) is available. We can measure n s in the σ units: Here k s is a positive constant and σ is the standard deviation of the distribution of the number of background trajectories passing through a pixel 4 . Now suppose that there is indeed a source present in a pixel B r . How can we miss it? Since our test identifies a source whenever n = n s + n b > µ + k t σ, the necessary condition for missing it is n b < µ − (k s − k t )σ. The probability of this happening (i.e., the false negative rate (fn rate)) is given by For instance, if k s ≥ k t + 3, the probability of missing this source is at most 0.0013. The sensitivity, or true positive rate (tp rate), is given by tp rate = 1 − fn rate.
Note that the sensitivity of the method depends on the threshold value k t as well as on our estimate of detected source particles.
From the above discussion we see that in order to find sources reliably, it suffices to detect n s > (k t + 3)σ ballistic particles from the source.
We claim that for any given SNR s > 0 and for any k t > 0 it is possible to detect the source, provided that the total number of detected particles, N , is sufficiently large.
Indeed, if a signal to noise ratio s is known, and if p 1, the number of detected source particles will be n s ≈ sN and σ ≈ N p(1 − p). Then, we require the following inequality to be satisfied: Since the left hand side of (9) grows linearly with N and the right hand side grows only as √ N , the inequality will be satisfied for a sufficiently large N . We come up with the following rule of thumb: One can expect to detect reliably the source with high sensitivity and specificity, if the total particle count N at the detectors is on the order of or higher, where s is the SNR (ballistic particles from the source vs the total hits) and p = r/R, where r and R are the linear dimensions of the source and the whole area correspondingly. Let us play with some practically possible values of the parameters. A nuclear source with shielding could have a radius of the order of several centimeters, while a vehicle or a cargo container has size in the order of meters. Thus, we will assume that R = 1, r = 0.01 ⇒ p = 10 −2 .
Substituting this value of p and s ≈ 10 −3 into 10, one expects that when N reaches close to 640, 000 and directional information at detectors is available, the sources should be reliably detectable. This squares well with the results of our previous computations with the backprojection methods. The conclusion from this probabilistic discussion is that the simple backprojection should reliably detect sources with low SNR, as long as the total number N of detected particles is high enough (see (10)). E.g., for the SNR of about 0.1% the values of N around 640, 000 should be sufficient. Besides, the confidence probability of detection (or absence of a source) can be computed.
The use of backprojection should also smooth out unstructured noise. In the next section we present numeric examples of such backprojection detections.
It is worthwhile to notice that, as long as the threshold parameter k t (equivalently, the desired confidence probability) is fixed, the whole detection procedure can be done automatically, without necessity of a human eye assessment of the reconstruction.
Comparison to simulation data. Recall that equation (6) and most of the following discussion assumed (incorrectly) independence for different pixels of the events "a pixel is crossed by n lines". To confirm that the confidences we estimated under this assumption are not vastly optimistic, we simulate a large number of random background fields and compare the statistics of these samples to our theoretical estimates. Each sample consists of 10 6 uniformly distributed random lines on a grid of 100×100 pixels. For this setting the mean and standard deviation of the number of lines crossing a pixel are approximately µ = 10, 000 and σ = 99.5. For different values of the threshold parameter k t ranging from 4 to 5, we record the true negative rate r t of the sample population, which is the ratio of samples in which each pixel had less or equal than n t = µ + k t σ lines intersecting. We compare this to the estimated confidences given by (7). The results from 50,000 samples are shown in table 3. It can be seen that for all values of k t the true negative rate from the sample population is slightly lower but very close to the estimated confidences. We conclude from this that the simplifying independence assumption made in this section is reasonable and the estimated probabilities agree well the true probabilities.  Table 2: Comparison between estimated confidences and statistics from 50,000 negative samples

Source Detection Using Backprojection
We illustrate the conclusions of the previous Section on a couple of typical examples of 2D X-ray data. The data was generated in the way described in Section 2.1. The incoming trajectories of detected particles were backprojected using a line drawing algorithm, as was explained in Section 3. Recall that we recorded approximate locations of the site at which a detection occurs, since in reality detector arrays consist of a number of bins (detectors).
In the first example, we simulated 639, 954 background particles and 640 particles from a source located at (−0.43, −0.11). All particles were backprojected and the source was detected with 100% confidence (see Fig.14).
For our second example, we used only three detector arrays along three sides of a square. Such a "detector gate" is a more realistic configuration than detectors surrounding the whole object. Figure 15 shows a typical backprojection reconstruction. The drop of the values at one side of the square is due to detectors missing there. Thus, effectively, there is a non-uniform distribution of the background, which contradicts to our assumption of uniformity. However, this non-uniformity is smooth (slowly varying), so we can assume that locally the distribution of background trajectories is uniform. In order to apply our probabilistic arguments described in Section 3 we estimate the local mean and standard deviation at every point. For the reconstruction shown on Fig.15, this was done by analyzing the data on a 7 × 7 sub-grid surrounding each pixel. Then, the local means were subtracted and the resulting values plotted (see the bottom left picture in Fig. 15) in the units that are equal to the local standard deviation σ. Finally, the values below the threshold k t = 4 were erased (bottom right in Fig. 15). The result clearly shows the detected source. The confidence probability was calculated as where k source is the height of the highest peak. The calculated confidence level for this example was 99.99%. Notice that if we used the threshold value k t = 4 instead of the actual peak's height, we would have significantly underestimated the confidence.

Compton Cameras
Compton cameras, also called Compton scatter cameras or electronically collimated cameras, have received a lot of attention since they were first proposed in [41] and a prototype was developed in [36,37]. Compton cameras were first suggested for use in nuclear medicine (SPECT), but they also have applications in astrophysics [14,38,44], monitoring nuclear power plants [33,34], and other areas. Compton cameras offer several advantages over the conventional mechanically collimated (Anger) γ-cameras.
The most significant of these is the dramatic increase of sensitivity -Compton cameras are reported to count orders of magnitude more photons than Anger cameras [6,24]. Another advantage is the flexibility of geometrical design [33,34], which even allows for hand-held devices [17,23]. A Compton camera consists of two detectors, which, in the standard setup, are planar and are placed one behind the other (see Fig. 16). A photon incident on the camera undergoes Compton scattering in the first detector, which records the location x 1 and the energy of interaction E 1 . After the scattering, the photon is absorbed in the second detector, where the position x 2 and energy of absorption E 2 are measured. From the knowledge of E 1 and E 2 , the scattering angle ψ can be determined by the formula (e.g. [7,41]) Here m is the mass of an electron and c is the speed of light. The direction into which a photon scatters is given by −β := From the knowledge of β and the scattering angle ψ we conclude that the photon originated from the surface of the cone with central axis β, vertex x 1 and opening angle 2ψ. Therefore, although the exact Figure 16: Schematic representation of a Compton camera.
incoming direction of the detected particle is not available, one knows a cone of such possible directions. The goal of Compton camera imaging is to recover the distribution of the radiation sources from this data. The description above addresses the Compton γ-cameras only. However, neutron detectors are currently being developed [25,40] that (employing a different physics principle) provide the similar cone information. We thus will not make distinction between these and call them all "Compton type cameras." We will now address briefly the mathematics behind the Compton type measurements.

The Cone Transform
Suppose that we want to image the distribution f (y) of radioactivity sources inside a certain object in R 3 . From now on, we will assume that f (y) is supported on one side of the Compton camera and that suppf does not intersect the camera. Under the assumption of a sufficiently long observation time (and thus a large number of particles detected), the Compton camera provides us with the projections [1,7,43], which we denote by Cf (x, β, ψ), i.e. the integrals of f (y) over cones parameterized by a vertex x lying on the detector, central axis vector β from the unit sphere S 2 , and half-angle ψ ∈ [0, π] (see Fig.17): Here K(ψ) is the Klein-Nishina distribution of scattering angles 5 . The function Cf (x, β, ψ) provides the expectation of the number of hits, per unit time, from the particles moving along the cone. In particular, when the particle count is low, these values are not well determined by the data. This issue was addressed in the preceding Sections. An important observation is that the problem of inverting the cone transform Cf (x, β, ψ) is even more overdetermined than the 3D X-ray transform. Indeed, f (y) is a function of three variables, while Cf (x, β, ψ) depends on five parameters (here x is restricted to a Figure 17: A Compton cone with apex x, central axis β and half-angle ψ. 2D detector surface). Most authors consider restricted versions of the cone transform, reducing the number of parameters from five to three [1,7,29] . The paper [39] contains a closed form inversion formula, which uses four of the parameters 6 . We, however, have to use the full 5-dimensional dataset, since, as it was mentioned before, any attempt to restrict the data would wipe out the signal.

Reconstruction Techniques in Compton Camera Imaging
We provide here a brief survey of some known reconstruction methods in Compton cameras imaging. Many of them reduce cone projections to Radon projections, i.e. integrals of f over planes.
Algebraic Reconstruction Techniques, as elsewhere in tomography, are widely used in Compton cameras imaging. This is partially due to the fact that analytic inversion formulas were not available until the second half of 1990's. Maximum likelihood methods were developed in [3,8,15]. A maximum likelihood method was also used in imaging with the COMPTEL telescope [38]. In [35] a matrix inversion technique was utilized to solve the problem for a specific detector geometry. Fast backprojection algorithms were developed in [32,45].
A spherical harmonics expansion solution was proposed in [1]. For every point x on the detector array, and a fixed in advance half-angle ψ the authors relate the conical projection Cf (x, β, ψ) to the Radon projection (integral) of f over the plane passing through x and perpendicular to β. A fast algorithm for computing the spherical harmonic series was also developed in this paper. Several other authors provided spherical harmonics solutions, however their model for Compton data differs from (13), e.g. [19,31,42].
Probably the first closed form analytic inversion formula was obtained in [7]. The authors considered only cones perpendicular to the detector plane, i.e. cones with central axis β = z. Thus, the Compton data depends only on the vertex x := (x, y) and the opening angle ψ. Let us denote g(x, y, ψ) := C f (x, z, ψ) and let F (u, v, ·), G(u, v, ·) be the two-dimensional (x, y)-Fourier transforms of f (x, y, z) and g(x, y, ψ) respectively. Then the following relation holds: where ξ = z √ u 2 + v 2 and t = tan ψ. Here H 0 is the zero-order Hankel transform and J 0 (γ) = (2π) −1 2π 0 e iγ cos φ dφ is the zero-order Bessel function. Later, in [30], more properties of the same restricted cone transform were given. In [26] the authors use yet another mathematical formulation for the forward problem, which accounts for the efficiency of the detector at different incident angles of particles. They extended the approach of [7] to show that the radioactivity distribution can be reconstructed from any family of cones that have central axes β such that the angle between β and the z-axis is fixed. Then averaging of the solutions for different values of this angle was employed, thus making use of all available data and reducing noise.
The following useful relation between the cone and Radon transforms was found in [39]: Cf (x, β, ψ)h(cos ψ) dψ (15) Here H denotes the Hilbert transform and R is the Radon transform in three dimensions. Note that four out of five variables are used and that the HRf (β, x · β) is constant on the line L c := {x ∈ detector plane|x · β = c}, for any fixed constant c. The last observation shows that it is possible to reconstruct f if the vertices of cones x are restricted to a curve. This also allows for averaging of the solution on the lines L c . In Section 6.1 we will consider a two-dimensional analog of this formula.

Compton Cameras in Two Dimensions
In the two dimensional setting we assume that the detectors x lie on a line, which we call the detector line. A cone in two dimensions is defined by a point x that serves as its vertex, a central axis unit vector β, and a half-angle ψ. Such cones simply consist of two rays with a common vertex. More precisely, let β = (cos β, sin β) ∈ S 1 , β ∈ [0, π] be the central axis of a cone with opening half-angle ψ ∈ [0, π]. The two half-lines that constitute the cone pass through the detector point x and are given by where α 1 = (cos(β − ψ), sin(β − ψ)) and α 2 = (cos(β + ψ), sin(β + ψ)), see Fig. 18.
For convenience, we will sometimes express Cf and Df as functions of the angle β instead of the unit vector β. Thus, we can write (16) as As in 3D, the problem of inverting the two-dimensional transform is overdetermined: Cf (x, β, ψ) depends on three parameters, while f (y) only on two.
The authors are aware of few papers devoted to Compton cameras in two dimensions. In [2,13] subsets of the cone projections of f (x) are represented as line integrals of the function f added with its mirrored shear transformation. Namely, let us introduce the variables k 1 = cot(β + ψ), k 2 = cot(β − ψ). The cone transform can be written as Cf where the x− axis coincides with the detector line and the z− axis is perpendicular to the detectors. Let us fix the sum k 1 + k 2 = K and define the function k 1 z, z). Thus, for each fixed value of K, the cone projections are represented as line integrals of a certain function, so all that remains is to invert the Radon transform. The reconstructed image will contain the function f above the detector line and a mirrored shear transformation of f below the detectors. Obviously, by fixing the parameter K, not all available data is used. The case when K = 0 corresponds to using only those cones with central axes orthogonal to the detector line. Since we need to use all overdetermined data, this approach is inconvenient. Besides, as our previous consideration of tomographic methods shows, in our specific situation, backprojection methods, rather than attempting the "full" reconstruction should be used.

Method A
The following relation, proven in [18], makes use of all parameters and is an analog of the three dimensional result presented in [39].
Cf (x, β, ψ)h(cos ψ) dψ (20) Here H is the Hilbert transform in the first (scalar) variable and the integrals are understood in the principal value sense.
Here B 2 is the unit disk in the plane and H s 0 is the standard notation for Sobolev spaces (e.g., [9]).

Corollary 5. Theorem 4 provides an inversion formula for the cone transform.
Indeed, computing the integral in the right hand side of (20), one recovers HRf , where R is the 2D Radon transform, and H is the Hilbert transform with respect to the linear variable. Then, the filtered backprojection formula for inversion of the two-dimensional Radon transform, implies that after differentiating the right hand side of (20) with respect to the linear variable and backprojecting, one recovers the function f .

Remark 6.
In the left hand side of (20), one sees Rf (s, β), where s = x · β. Thus, if one infinite detector array is used and β is perpendicular to the array, we would only know Rf (s, β) for a single value of s. If the detector array is finite, as is the case in all implementations, we would miss a much bigger chunk of data. One possible way to obtain the complete Radon data in (s, β) ∈ [−1, 1] × S 1 is, for example, to use three finite size detector arrays placed along the sides of a square containing the object.

Method B
Another, simpler way of obtaining Radon data from cone data is presented below. Let us denote by u(x, α) the integral of f along the ray starting at the point x in direction of the unit vector α (i.e. the fanbeam projection of f ): Let us fix a detector position x and a direction α 0 and try to recover the fanbeam data u(x, α 0 ). For any ψ ∈ [−π, π] we can write Cf (x, β, |ψ|) = u(x, β − |ψ|) + u(x, β + |ψ|). Thus, the following relation holds u(x, α 0 ) = Cf (x, α 0 + ψ, |ψ|) − u(x, α 0 + 2ψ) Figure 19: In order to determine the integral of f along the line l 0 , add the integral over the cone consisting of l 0 and l 1 to the integral over the cone determined by l 0 and l 2 . Then subtract the integral over the cone determined by l 1 and l 2 .
It is easily seen now (Fig. 19) that for any two half-angles ψ 1 , ψ 2 ∈ [−π, π] the fanbeam data u(x, α 0 ) can be obtained from cone data as follows: The angles ψ 1 and ψ 2 were arbitrarily chosen, so one could average on ψ 1 and ψ 2 in order to use all available data. This averaging also helps reducing the effects of the background noise.
If averaging over half-angles ψ 1 and ψ 2 is used, this method is computationally expensive. We propose a third reconstruction method, which is fast and makes use of all available data.

Method C
Let us recall that the source intensity function f (y) is compactly supported and suppf lies on one side and away from the detector array. Let us fix a detector location a. Then the function u(a, α) (see (22)) is supported in (0, π), so it can be extended periodically in α from [0, 2π] to (−∞, ∞). Let us definẽ We can write, In the above equation, the periodicity of u(a, α) with respect to α was used. Then, u(a, α) = u(a, α) Here, as before, R # denotes the backprojection operator The second term in the right-hand side of (26) can be written as where δ is the Dirac's delta-function. Therefore, we come to the following conclusion: Theorem 7. When the inverse fan-beam transform is applied toũ, the result equals the sum of the function f (y) and a distribution supported on the detector array only.
Hence, functions f supported away from the detector arrays are recovered correctly.
In essence, this method approximates an integral of a function on a given line by averaging all its cone integrals which have one side lying on the line.

2D Compton data backprojection examples
In order to generate Compton camera data, we first generated X-ray data in the manner described in Section 2.1. Recall that the X-ray data for each particle consists of a bin on a detector array and exact incoming direction α. Then, for every particle we chose a scattering angle ψ drawn from a uniform distribution in (−α, π − α). The central axis of the scattering cone was computed as the sum β = α + ψ. The information provided by the Compton camera is three dimensional and consists of the bin on the detector array, the central axis of the cone β and the absolute value of the scattering angle ψ.
For the reconstructions from Compton camera data presented in this section, we used Method C, described in Section 6.1.3, to convert the cone data into X-ray data. Method C was chosen as the least computationally expensive and also using all cone data. Note that for the particular type of data we are handling, the central axis and scattering angle of cones are continuous random variables. Thus, for a finite number detected particles, one cannot average on cones with a common side, because such cones are a rare occurrence. Therefore, the method is equivalent to simply treating a particle coming from a cone (x, β, ψ) as two separate particles detected at x and having incoming directions β − ψ and β + ψ. After converting cone data into X-ray data in this way, we employed the usual reconstruction procedures for X-ray data.
It is also important to make clear that in our experiments we assume that exact cone axes and scattering angles are known. The effects of uncertainty of these quantities are to be investigated in the future.
First, we used filtered backprojection to reconstruct the source distribution from data provided by four Compton detector arrays placed on the sides of the square [−1, 1]. The signal to noise ratio was fixed at 0.1% and we varied the bulk number of detected particles. As for the case of X-ray data, we defined a detection successful provided that highest peak of the reconstructed image occurred at the location of the source. The number of successful FBP detections increased slowly as the bulk number of detected particles went from 600, 000 to 1, 000, 000. The results are reported in Table 7. An example of a successful detection using Method C combined with FBP is shown in Fig. 20 Table 3: Number of successful detections from cone data of a source located at (0.311, −0.433), with a fixed SNR 0.1% and varying number of detected particles, using different reconstruction methods. For each level of background particles 20 sets of random data were generated. In each case we applied Method C combined with FBP and backprojection followed by local grids of size 9 × 9. In all but one experiment, if an FBP detection was successful, so was the corresponding backprojection detection. The one example in which FBP detection occurred but backprojection failed was at the level of 700, 000 background particles. Next, we use Method C for converting cone data into X-ray data, followed by backprojection based on the line drawing algorithm described above. Thus, for every detected particle, two particles are backprojected along the sides of the cone. Therefore, the value of the resulting image at every pixel equals the number of cones intersecting it. On Fig. 21(top) we show a typical backprojection reconstruction from four Compton detector arrays from exactly the same data used for the previous example (Fig. 20). The background is not uniformly distributed, so we employ the local grids approach introduced in Section 4. The size of the local grids was chosen to be 9 × 9 pixels. On Fig. 21(bottom) a plot of the number of local standard deviations away from the local mean is shown. The value at the source location deviates 4.93 standard deviations from the mean, and the estimated confidence is 99.59%.
In Table 7 we compared FBP reconstructions with backprojection reconstructions. For all backprojection results we used local grids of size 9 × 9 pixels on which the local means and standard deviations were estimated. In order to compare backprojection to FBP, we deemed a detection successful, if the number of local standard deviations was highest at the location of the source. Our results show that backprojection outperforms Figure 21: Unprocessed backprojection reconstruction from Compton data (top). The number of background particles is 900020 versus 900 source particles. The bottom left shows deviations from the (local) mean, measured in the units of (local) standard deviations. The size of local grids is 9 × 9 pixels. Only the peaks that exceed the local mean more than 4.3 (local) standard deviations are shown on bottom right. Source located at (0.311, −0.433) is detected with confidence 99.59%.
FBP, as in the case of X-ray data. The results are summarized in Table 7.
As we have already seen, detection based on backprojection allows for estimation of confidence. Similarly to detection by X-ray backprojection, one could set a threshold for detection at 4.1 or 4.3 local standard deviations above the local mean. According to our estimate for the confidence (7), such thresholds would result in confidences at least 81% or 91.8% respectively. In Table 7 the number of successful detections for these two values of the threshold parameter are given. We used the same simulation data as for Table 7. threshold \ bkgd particles 600, 000 700, 000 800, 000 900, 000 1, 000, 000 4.  Table 4: Number of successful detections from cone data with a fixed SNR 0.1% and varying number of detected particles. Backprojection was used and local means and standard deviations were estimated on local grids of size 9 × 9. Two different threshold values for detection were set, resulting in confidences at least 81% and 91.8% respectively.
In comparison with X-ray data, a larger number of background particles is needed for detection with Compton cameras. This reflects the somewhat lower quality of Compton data in comparison with the data from collimated detectors. A couple of examples of reconstructions from three Compton cameras, which form a "gate" of detectors, are provided below.
The source in Fig. 22 is close to the side of the square where there is no detector, which makes it harder to observe. The number of ballistic source particles was 1050 versus 999, 925 background particles. The local mean and standard deviation were estimated on 9 × 9 local grids. The maximum number of local standard deviations was 5.56 and occurred at the location of the source. The estimated confidence of detection is 99.98%. Figure 22: Backprojection reconstruction from Compton data (top). Detected source particles -1050, background particles -999, 925, source location -(−0.503, 0.11). Number of (local) standard deviations above the (local) mean at each pixel is shown on bottom left. The size of local grids is Peaks that deviate more than 4.3 (local) standard deviations from the mean are shown on bottom right. The source is detected with confidence 99.98%.
Finally, we present an example of a successful detection of several sources (Fig. 23). Three sources were placed inside the imaged area, and three Compton-type detector arrays were used to detect particles. The total number of detected ballistic source particles was 3, 229, each of the sources having roughly the same contribution. Additional 999, 325 random background particles were detected. We have set up the threshold to k t = 4.3σ, which gives confidence of detection at least 91.97%. Note that this confidence probability is estimated based on the threshold, and not on the actual number of deviations at each peak, as in the previous examples. Thus, individual confidence probabilities for the three sources discovered are even higher. The local grids we used for estimating the local means and standard deviations had dimensions 9 × 9.

Final remarks and conclusions
Here are the main conclusions that we have reached: 1. Direction sensitive detectors are needed to be able to uncover existence of a low emission source in the presence of dominating background noise. Collimated cameras, although providing the most valuable direction information, are not applicable, since collimation would decimate the already low signal. Compton type cameras can be used instead.
2. Although standard integral-geometric models of emission tomography (X-ray and Radon transforms and their attenuated versions) are not applicable to the situation of a very weak source, they work to some extent in the detection problem. However, they fail way before reaching the low values of SNR needed in applications.
3. In fact, only the backprojection part of the tomographic reconstruction does the detection work, and the filtering part worsens the results. 4. The (unknown) attenuation might lower the SNR, but otherwise is irrelevant for the validity of the backprojection method.
5. The conclusions above hold both for the collimated and Compton type cameras.
6. The backprojection procedure can be justified by a simple probabilistic consideration, which also provides confidence probabilities of detection.

7.
A simple algorithm is provided that recovers the source from Compton data at SNRs on the order 0.1% and beyond. It also provides confidence probabilities. In theoretical considerations and numerical tests, the algorithm demonstrates high sensitivity and specificity. The procedure can be made completely automatic, without the need of a human eye assessment of the image.
8. The algorithm is based upon the assumption of an uniform random background. However, the non-uniformity that arises when the object is only partially surrounded by detector arrays, can be handled successfully by using local grids.
9. As it is done in SPECT, one can try to use Bayesian methods to approach the same detection problem. This was done in 2D case in [46].
Certainly, there remain quite a few questions that need to be resolved. Some of them will be discussed in the next publication. In particular, 3D Compton backprojection algorithms will be developed and studied analytically and numerically; more complex cargo structures will be modeled and tested, which will involve non-uniform attenuation and scattering; local grid version of the algorithm will be tested on more complex nonuniform random backgrounds than the ones arising due to the partial view; dependence of the results on the precision of Compton type cameras will also be addressed. We assumed that the exact incoming direction of particles is detected. In practice this is usually not the case, as only a probability distribution of the directions may be known. The effect of directional uncertainty is important from practical standpoint and we plan to consider it in the future.