Home
This package is aimed at amateur astrophotographers. Its purpose is discovery (serendipitous or targeted) of faint moving star-like objects in imaging sequences taken with a telescope and an astro camera. These could be asteroids, faint comets, potentially interstellar visitors. In most cases, the objects are too faint to show up in individual images.
With the addition of the time-binning feature (since v0.3), long image sequences, spanning days, weeks, or even months, are now possible. This opens a window into the detection/discovery of very slowly moving and faint objects such as Centaurs and Kuiper belt objects.
The requirements: a Windows or Linux computer with an NVIDIA GPU, reasonably fast, with at least 8GB of RAM. (Though you may still be able to do some work with a 4GB card).
The software is command line only (no GUI). Only FITS files are supported (RGB or monochrome).
At the core of the package is a highly optimized CUDA kernel which performs bilinear interpolation and stacking of the images very fast. For example, a sequence consisting of 134 images 20 megapixels each can be interpolated and stacked in 10 milliseconds (on H100 GPU). Additional features (spatial and time binning) can make the process much faster. This makes it practical exploring all possible moving objects paths in an image sequence.
The package consists of two principal steps/programs:
1) reconv: removing all static objects in the imaging sequence (stars, galaxies, nebulae, sky glow background etc). The majority of this is subtracted (so the data is preserved), with the rest (like residual star artifacts) masked out. At the end, the imaging sequence should look like plain noise, with some masked areas, and occasional visible moving objects - like artificial satellites and meteors.
2) asteroid_search: the most compute intensive part. Brute search for moving objects in the imaging sequence, by trying lots (thousands to millions) of different motion vectors, then stacking the whole image sequence along each vector, and looking for statistically significant brighter pixels. (This is sometimes called "synthetic tracking", though personally I don't like the term.) At the end, 4D clustering analysis is performed, to find the moving objects ("clouds" in 4D space). The 4D space consists of the x,y coordinates of the pixel in the first image, and the motion vector coordinates (Jx,Jy). The code is using GPU acceleration to make this practical. I am using CUDA, so only NVIDIA GPUs are supported.
To illustrate the software, I will be using one specific image sequence of my own: 142 images derived with SV705C camera (no cooling; with my RC8 telescope), 30 seconds exposures, taken on July 14, 2025, centered at the just-discovered interstellar comet 3I/Atlas. At that time the comet was too faint to be seen in individual images (the magnitude was ~17.5m). The images are 8 million pixels each, with FWHM=7.7 pixels.
The image series (only FITS files are supported) needs to be calibrated, registered (on stars), and optionally cropped (to get rid of edge artifacts if present). You can also do hot pixels removal if they are a problem (it is for my old Canon 6D DSLR, but not for my SV705C). Only the first image in the sequence needs to be plate solved and saved. No color calibration or background removal is needed. Finally, you should stack you sequence, to produce the "master.fit" file. I use the open source Siril package for all these steps. When stacking, I use "Average stacking with rejection" and "No normalization". In the Stacking result section, I use "Force 32b output". Do not enable "RGB equalization" - this will mess up the following steps!
This is how the master.fit file looks like for my image sequence. It's in the middle of Milky Way, there are at least 5000 stars visible here. The 3I comet is somewhere near the center, but it is too faint (17.4m) and moves too fast (~200 pixels over the image sequence) to be visible here.
This step is to subtract plus mask out (if the subtraction step leaves some artifacts) all static objects in the image sequence: stars, galaxies, static background etc. We are using the code reconv for this. It is CPU-only (doesn't need a GPU), and can take from a few seconds to ~30 seconds per image in the sequence.
It takes lots of optional arguments, but the idea is you need to figure them out only once for a given telescope + camera combination, and then simply re-use the settings.
Your goal here is to get rid of any visible signal (including stars, galaxies, and nebulae, but also remaining hot pixels and cosmic ray artifacts) from the images, leaving behind only noise and potentially some masked areas (typically no more than 3% of the total number of pixels).
Be aware that the output images are black and white, and have some Gaussian blur applied (to maximize the signal to noise ratio when detecting moving objects). I find it very convenient to analyze the output images using the free software ASTAP (it automatically rescales the brightness for each image, so it's easy to see if some images still have issues like residual star artifacts, hot pixels etc - such images will appear much darker than the rest).
For a new telescope + camera combination, you proceed in stages, until you figure out the optimal reconv parameters.
First, try the most basic operation: getting rid of static objects by subtracting a scaled and "reconvolved" version of the master.fit image from each input image. Both images are first blurred, using FWHM (an input parameter) as the blur parameter.
The algorithm works as follows. First, we derive the 2D kernel which converts the master image into the current individual image from the sequence. (The kernel is derived by doing FFT for the individual image, for the master, then dividing the former by the latter in complex space, then doing an inverse FFT.) This kernel contains all the imperfections of the individual image: guiding errors, air turbulence etc. The kernel is not axially symmetric. Next, we truncate the kernel at some fairly small radius (~2-4 FWHM units). Next, we convolve the master image with this kernel (using FFT), which produces a degraded quality version of the master image, with all the imperfections of the individual image (guiding errors etc.). We then compute the brightness scaling coefficient between the convolved master image and the individual image (only the pixels well above the noise floor are used here), and then subtract the scaled version of the former from the latter. This procedure gets rid of all the static objects (which are present both in the individual image and the master image), while leaving the moving objects data more or less untouched.
You need to measure FWHM (full width half maximum) value for the stars in your master.fit file, in pixels. You can do it in Siril: open master.fit file, go to Tools > Image Analysis > Dynamic PSF menu item, click on the detect stars icon at the bottom, then on "Average PSF parameters" icon. You will get the average FWHMx and FWHMy values, in arc seconds. You combine them as sqrt(FWHMx^2+FWHMy^2), and divide by the arc second/pixel value for your image, to get FWHM in pixels.
Run reconv with your image sequence as follows:
reconv -m master.fit -o out_ -FWHM 7.7 -R 2 -bias 0.003 -border 1 my_images_*.fit
Use your own value of FWHM. For the -bias parameter, use the value appropriate for your camera (assuming the full brightness range is 0..1). It doesn't affect the search, just the way your images are displayed in software like Siril and ASTAP (they both don't like negative brightness values, so some bias is needed). I use 0.003 for my SVbony cameras, and 0.03 for my Canon 6D camera.
The "-border 1" argument normally doesn't need adjustments. By default, it masks out the output images edges with the width = 1 in FWHM units. This is to cover some residual (and rather subtle) artifacts left from the Gaussian FFT smoothing of the images.
The most important parameter here is -R. Use -R 2 initially. The number (2 in this case; it can take fractional values) is the truncation radius (in FWHM units) of the kernel (see above). This parameter affects the quality of subtracting the static objects (like stars) from individual images, using the master image as the template. If the radius is too small or too large, you will see stellar leftover halos and/or very deep shadows around bright stars. Choose the value which minimizes the amount of artifacts around stars, especially paying attention to brighter areas (dark artifacts are much less problematic for the moving objects search).
Test multiple images in the series, with different amounts of star blur. Choose the -R value which works the best overall, most importantly - minimizing the amount of residual bright halos around stars.
For our example sequence, here is the case of using too small value of R = 1 (a zoomed in view). Star artifacts are very obvious:
Using R = 2 (the best value, star artifacts have much smaller dynamic range):
Finally, when increasing R even further (R = 4), star artifacts become bright again:
So the optimal value of R here is close to 2.
By default, the output images have the background removed, using bilinear model with 5 tiles along the shorter side of the image. If that is not enough, you can request more tiles for the background model with the switch "-bg". For example, "-bg 11" will use 11 tiles along the shorter side of the images. For the longer side of the images, more tiles will be used, keeping the shape of each tile close to a square.
Executing the command "reconv" without any arguments will print the full list of the options supported.
If your images have great quality (low turbulence, perfect guiding, identical star aberrations across the frame), this may be the only stage you will need. In most cases though you will need to proceed to the following stage(s).
If the stage 1 didn't produce an acceptable outcome (there are still some bright pixels around brighter stars), you need to mask out those artifacts.
Start by adding "-mask" argument to the reconv command. The numerical argument is the number of std units above the noise floor in the master image where any pixels brighter than that will be erased (masked out). Important point: masks are computed based on the master image, but applied to individual images. The std value is measured from the master image. The good starting point is "-mask 100", you can lower down the number if the bright artifacts around brighter stars are still present.
For our example sequence, using "-mask 100" (some limb brightening around stars is still visible):
Using "-mask 50" does a better job, but there still may be some residual limb brightening which might mess up the following motion detection step:
If a complete removal of the star leftovers requires you lowering the "-mask" argument to a very low value (say, less than 20), consider adding another argument instead, "-grow_mask". When you do that, the existing star masks are expanded, which may be enough to cover all stellar artifacts. Try "-grow_mask 1" first, then adjust the number to try to mask out all bright pixels artifacts. (E.g., try 0.5, 1, 1.25, 1.5 etc.) This approach allows one to minimize the number of masked pixels, as we are only working with the brightest star artifacts.
Here is our image sequence with "-mask 50 -grow_mask 0.25". The artifacts are pretty much gone at this point - all we see is the (blurred) Gaussian noise, plus some masked areas (1.8% of all image pixels are masked - well within the 3% range).
Here is the full-frame image - it looks very clean. Out of ~5000 stars present in the original images, the vast majority were successfully subtracted (with the data from moving objects preserved), and only the remaining ~200-300 brightest stars were masked out.
Finally, if everything else fails, you can use the "-hot_pixels" switch for reconv. As the name suggests, it can be used to deal with hot pixels (though it is always best to try to remove those during Siril images calibration.) But what it does is erasing (masking out) all the pixels a certain number of std (measured for the individual image noise) above the noise floor, so it can be used for any leftover artifacts - hot pixels, cosmic rays, bright satellites etc. It is not ideal, as some faint (but visible in individual images) moving objects may also end up being erased. Use something like "-hot_pixels 4" to effectively erase anything visible which is not noise. It is not recommended to use the values below 4, as you will start erasing noise pixels (so you will end up erasing too much valuable data).
This is the most compute intensive (and time consuming) part. It requires a GPU (only NVIDIA GPUs are supported). The whole image sequence needs to be uploaded to the GPU memory, so you need to pay attention to how much GPU memory you are using.
It is usually impractical (too slow, uses too much GPU memory) to search for moving objects directly using high quality search parameters. Instead, you should do it in two stages: draft, then fine-tune.
Here we use lower image resolution (using the -rebin argument), lower time resolution (using the -time_bin switch), and/or lower motion vector resolution (using a fairly large value for the -Step switch, say 1-2, versus the high-resolution values under 0.5) to make the search faster, and to fit the images into the GPU memory. The parameters you use should be still good enough to discover faint moving objects, may be just with a slightly less signal to noise ratio compared to a full resolution search.
In my experience, it is ok to rebin the images down to FWHM spanning ~2 pixels in the rebinned image. E.g. if FWHM=4.27 pixels in the original images, I can use "-rebin 2" (4.27/2 ~ 2.1 pixels per FWHM in the rebinned images). For my SV705C camera (FWHM=7.7 pixels), I use "-rebin 3" (7.7/3 ~ 2.6 pixels per FWHM). Rebinning is done in both dimensions, so "-rebin 3" reduces both the X and Y dimensions by a factor of 3, which accelerates the computations and reduces the GPU memory consumption by a factor of 3x3 = 9.
For the "-Step" option, a draft search value could be 1.5 or so. (For high quality search, the value 0.5 or smaller would be needed.) The computational time scales as the square of the Step argument, so going from "-Step 0.5" to "-Step 1.5" will accelerate the computations by a factor of (1.5/0.5)^2 = 9. There are no GPU memory changes when changing the Step value.
Finally, if your image sequence spans more than one night, you should definitely consider using the option "-time_bin". It takes one numerical argument: the largest amount of motion blur (in FWHM units) you can tolerate in each composite (time-binned) image. This value is measured for fastest moving objects, so the amount of time binning you achieve depends directly on the value of -RMAX (see below). If you are using fewer than 5 nights of imaging, you may want to set -fraction to a value in the 0.5...0.6 range (the default value is 0.75), to minimize the masked area in the stacked images. For two nights of imaging, use "-fraction 0.4" or so.
It is important to monitor what happens to asteroid detection when using -rebin, -time_bin and -Step options. As these parameters increase, the detection significance goes down. It is ok to loose ~5-10% of the signal-to-noise ratio when doing the draft search, but not much more.
A draft search run can look like this:
asteroid_search -draft -RMAX 0 -rebin 2 -Step 1.5 out_*fit
The option "-RMAX 0" means to search for the whole practical range of motion vectors, up to the motions spanning most of the shorter side of the image.
You can specify the range for the length of the motion vectors (in FWHM units) explicitly, using -RMIN and -RMAX switches. For example, "-RMIN 10 -RMAX 50" will test all motions when objects move between 10 and 50 FWHM units between the first and the last image. When -RMIN switch is omitted, it takes the default value of 1.
Be aware that the time spent on the draft search grows as a square of RMAX. Sometimes using the default "-RMAX 0" will result in an impractically long search, so you may have to use an explicit value of RMAX to speed things up. The good news: the program displays the estimated time for the run at the beginning (after reading the images and computing the brightness histogram), so you don't have to spend a day waiting in vain if you used a wrong (too large) RMAX value. It will also print an estimate on the number of motion vectors to be tested.
When the program finishes, it will create a few files.
list.dat: lists all pixels in stacked images above a certain threshold brightness level (p_min). Warning: the list can be very long (a million lines or more). Ideally it should not have more than 100,000 lines. If it is too long, you need to adjust the threshold (cutoff) pixel brightness p_min indirectly, using "-N_noise" option, or directly using the "-p_min" option. (See below.)
stats.dat: the result of applying the 4D clustering analysis to the above list. The file lists "clouds" (clusters of neighbouring pixels from list.dat in 4D parameter space), starting from the cloud containing the brightest pixels. Each cloud can represent a real moving object, or can be a star residual artifact, or a noise artifact.
cloud_XXX.fit: FITS files corresponding to N brightest objects/clouds found. By default N=10, you can change that with "-N_cloud" switch. E.g. "-N_cloud 25". Each file shows the stacked image corresponding to the best motion vector for the respective cloud. It also contains the special ANNOTATE keyword, which can be used to visualize the cluster location in the image inside ASTAP program. Simply open the cloud*fit image in ASTAP, and the cloud annotation will be automatically displayed. As cloud*fit files inherit all FITS keywords from the first image in the sequence (including the plate solving part), they are all plate solved. That means that in ASTAP you can use Ctrl-R key combination (or use the menu item Tools > Asteroid & comet annotations) to visualize all known asteroids and comets for this image. (If you use this ASTAP feature for the first time, you first need to download two files: MPCORB.dat and AllCometEls.txt; the instructions are right there.)
The data from stats.dat and list.dat can be visualized using the provided Octave script stats.m (see below).
If even the top of the stats.dat file is dominated by star artifact "clouds", you should go back to the "reconv" step, and apply more rigorous masking, to minimize the artifacts.
If list.dat becomes too long, you can reduce its size by using an explicit value of "-N_noise" switch. The default value is 10: "-N_noise 10". This is used as follows. Before performing the full motion vector search (which can involve a million motion vectors or even more), a preliminary short run (containing only 1000 motion vectors, randomly generated in the given interval between RMIN and RMAX) is performed. We compute a histogram of all bright pixels for these motion vectors, and choose such a cutoff brightness (p_min), that no more than 10 "noise" pixels are detected in the 1000 motion vectors. In other words, the chance of finding even one noisy pixel in any given motion vector is 10/1000 = 0.01.
If there are lots of artifacts in the image sequence, the default value for -N_noise may still be too large. You can decrease the size of list.dat file by making the argument even smaller, e.g. "-N_noise 2".
The opposite problem may occur - after running the full search, you find too few (say, <1000) or even 0 bright pixels. In this case, increase the N_noise parameter to 100 or even 1000.
The "good" number of bright pixels (and lines in list.dat) is from a few thousands to a few hundreds of thousands. More is usually better, as you can get a better sense of the noise statistics, but when the number of bright pixels becomes > 100,000, list.dat can become too large for post-processing, also the 4D cluster analysis step can become too long.
At this point, hopefully you have already found one or more clouds which either correspond to existing asteroids/comets, or are candidates for a new asteroid/comet.
Using our sequence example, I ran asteroid_search with the following arguments: "-draft -RMIN 1 -RMAX 150 -rebin 3 -Step 1.5". The computation took 5 minutes (on RTX 2080 Super), consumed 0.7 GB (out of 8 GB) of the GPU memory, processed 31,000 motion vectors, detected 21,000 bright pixels, clustered into 11,600 clusters. Visualizing the stats.dat and list.dat files (using Octave script stats.m) suggests there are two obvious moving objects (clouds) in the data, with the motion vectors R=106 and 440 pixels, with the rest of the clouds likely being noise:
The two top lines in stats.dat describe the parameters of these two clouds:
1 2 3 4 5 6 7 8 9 10 11 0 5.453109e+01 877 1531 -23.10 438.90 2.431836e-04 1261 3.909943e-01 187.84 81.67 1 4.600127e+01 1825 2602 57.75 92.40 2.431836e-04 223 7.009500e-02 114.34 40.84
Here is the explanation of the columns:
- 1: Cloud ID (brightest clouds are listed first)
- 2: The brightness of the brightest pixel in the cloud, in std units for the noise in the master image.
- 3: The x pixel coordinate of the brightest pixel in the first image of the sequence (corresponds to y coordinate in ASTAP)
- 4: The y pixel coordinate of the brightest pixel in the first image of the sequence (corresponds to x coordinate in ASTAP)
- 5: The x component of the motion vector, Jx, in pixels
- 6: The y component of the motion vector, Jy, in pixels
- 7: The p_min value used for finding the bright pixels, above the bias value (in absolute brightness units, range 0..1)
- 8: Number of pixels in the cloud
- 9: The "mass" of the cloud: sum of the brightness for all the cloud pixels.
- 10: The radial extent of the cloud in the (x,y) domain, in pixels
- 11: The radial extent of the cloud in the (Jx,Jy) domain, in pixels
This is how the stacked version of Cloud 0 looks like in ASTAP (a zoom in):
After using Ctrl-R (asteroid annotation), we can see that Cloud 0 is actually 3I/Atlas (shown here as MPEC 2026):
The Cloud 1 happens to be a known main belt asteroid 67375 (located near the bottom edge of the images). It's pretty small (3-7 km). According to the Horizons site, its visual magnitude was 17.7 mag at the time.
Now we need to measure the detected objects' parameters more accurately, using the fine tuning option of asteroid_search.
During the finetuning stage, we take one tentative object detection (near the top of stats.dat generated during the draft search), and instruct the code to search in the neighbourhood of the detected cloud (in 4D space) using much higher resolution settings - no rebinning, and smaller values of -Step (I typically use -Step 0.25 at this stage).
During the finetuning search (enabled with the -finetune option), you pretty much always want also to use another option, -crop. The -crop option results in dramatically smaller GPU memory consumption, and much faster computation. The only catch: in the cloud*fit files, only a small area along the path of the object is computed/shown, the rest is masked out. As long as the crop size is large enough, this should not be a problem.
The -finetune option takes 4 numerical arguments (4 adjacent columns in stats.dat file, columns #3,4,5,6) - x,y coordinates of the brightest pixel in the cloud, and Jx,Jy motion vector coordinates corresponding to the most significant detection. You simply paste and copy these 4 values from the draft search stats.dat file for the cloud you are interested in.
The parameters -RMIN and -RMAX have a new meaning in the -finetune mode. The first one describes the size of the square search window in the image space, while -RMAX describes the size of the square search window in the motion vectors space (both are using the units of FWHM). E.g. "-RMIN 20 -RMAX 20" means that the search window (centered on the moving object) is +-10 FWHM units relative to the object's center, while the motion vectors will be +-10 FWHM units relative to the best draft vector detection.
This is how you typically run the code in finetune mode:
asteroid_search -crop -finetune x y Jx Jy -RMIN 20 -RMAX 20 -Step 0.25 -N_noise 1000
You may have noticed the larger value of -N_noise compared to the default value of 10. This is because doing the crop search results in most of the pixels in the stacked images containing a signal, with almost no pixels containing pure noise. This affects the histogram method of estimating p_min.
An alternative to using -N_noise with a larger numerical parameter is to use the p_min value explicitly, with "-p_min value" option. You can use the p_min value computed during the draft search - the column #7 in stats.dat . For example:
asteroid_search -crop -finetune x y Jx Jy -RMIN 20 -RMAX 20 -Step 0.25 -p_min 9.150479e-05
In the finetune mode, the output files have different names: list_fine.dat, stats_fine.dat, cloud_fine_*.fit.
At the end of each fine-tune run you typically only want the top cloud information - the top line in stats_fine.dat, and the FITS file cloud_fine_000.fit . This is the high-resolution detection of the cloud requested.
Sometimes stats.dat doesn't have a single large cloud at the top, but instead has 2 or more comparable size clouds. That means that your parameters were too restrictive, not enough to capture the cloud in its entirety as a single object. Instead, you detected a few disconnected segments of the cloud. To address this issue, repeat the finetune search using more generous search parameters - larger values of RMIN and RMAX, and/or larger value of -N_noise (which lowers down p_min brightness threshold, allowing more pixels to be used for the clustering analysis).
For our example image sequence, I ran the following command for the top detection (Cloud 0; 3I/Atlas comet):
asteroid_search -crop -finetune 877 1531 -23.10 438.90 -RMIN 25 -RMAX 25 -Step 0.25 -p_min 2.431836e-04
The run was fairly short (~5 minutes), and produced the following top lines in stats_fine.dat file:
0 1.454622e+01 875 1529 -21.17 442.75 2.431836e-04 486580 1.507581e+02 155.17 144.29 1 6.815159e+00 849 1496 51.97 548.62 2.431836e-04 16311 4.535503e+00 46.28 43.56 2 6.179550e+00 818 1558 77.00 390.77 2.431836e-04 11866 3.119489e+00 53.09 49.00 ...
As you can see, the top cloud 0 is by far the largest - almost 0.5 million pixels, versus 16,000 pixels in the next one, suggesting this is our object's cloud in its entirety. The fine-tuned parameters x,y,Jx,Jy are slightly different, compared to the draft search results.
This is how Cloud 0 looks like in ASTAP:
And here is Cloud 1 (asteroid 67375), fine-tuned:
0 1.298246e+01 1823 2603 61.60 92.40 2.431836e-04 99962 3.103748e+01 91.20 74.86 1 5.918282e+00 1854 2693 80.85 205.97 2.431836e-04 718 1.900434e-01 16.33 14.97 2 5.907478e+00 1594 2298 115.50 175.18 2.431836e-04 1476 3.888894e-01 16.33 14.97
This is mostly on my to-do list. I will write a script to analyze the high significance clouds (using all pixels from list_finetune.dat), to fit a 4D Gaussian to the top of the cloud. This will produce the estimate for the most likely cluster parameters (x,y,Jx,Jy), and also the confidence intervals for each of the parameter. This is critical when trying to deduce the asteroid's orbital elements.
For now I am providing a very simple script stats.m for Octave package. It is used to visualize the list.dat and stats.dat files. In many cases it can give you a good idea if the clouds you are detecting have a high signal-to-noise ratio (so are significant detentions), or if they are likely be a noise.
Another useful tool: annotate.c program (compile it with gcc or any other C compilers). The syntax:
annotate stats.dat annotate.txt r_min
The file stats.dat is from the draft stage of asteroid_search run. The text file annotate.txt has the following format:
ANNOTATE= '2902.62;2391.41;2922.62;2411.41;-1.3;05263;⚠ obsolete;05263;' ANNOTATE= '3591.98;859.96;3611.98;879.96;-1.3;14026;⚠ obsolete;14026;' ANNOTATE= '3806.51;956.18;3826.51;976.18;-1.3;16557;⚠ obsolete;16557;' ...
This is the list of known asteroids and comets for your imaging area. To get this information, load the first image from your imaging sequence into ASTAP program. It should be already plate solved. Then press Ctrl-R to bring up the Asteroid annotation window. Inside that window, enable the option "Add as annotation to FITS header", then click on "Annotate asteroids and comets". This will look for any known asteroids or comets within the image, displays them, but also adds the ANNOTATION keywords (one per object) to the FITS header (top right area of ASTAP program). Then simply highlight all ANNOTATE lines there, and paste and copy them to a text file.
The final argument, r_min, is the maximum distance between the clouds and the known objects, in pixels. It can be 10-40 pixels or so.
The program will print all the matches between the two files: all the clouds from stats.dat which are located very close to known objects.