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.
The software is command line only (no GUI). For now, the only input images suppoted are RGB fits files (from OSC cameras). (Though it'd be a very easy fix to support monochromatic cameras as well.)
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 (with my RC8 telescope), 30 seconds exposure, 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 RGB FITS files are currently 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 calibrarion or background removal is needed. Finally, you should stack you sequence, to produce the "master.fit" file. I use Siril package for all these steps. When stacking, I use "Average stacking with rejection" and "No normalization". In the Stacking result section, I use "Output normalization" and "Force 32b output". Do not enable "RGB equilization" - 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 detectable here.
This step is to subtract (ideally) or 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 signals (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 analysie 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 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 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, 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 to your camera (when 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 (which is 2048 in internal raw units).
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 use 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, you will see bright stellar leftover halos around each star. If it is too large, you may see other artifacts becoming pronounced (like dark - negative signal - halos around 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 amount of star blur (from air turbulence and imperfect guiding). Choose the -R value which works the best overall, most importantly - minimizing teh amount of residual bright halos around the 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 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 in 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 esle 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 to erase (mask out_ all the pixels a certain number of std (measured in the individual images) above the noise floor, so it can be used for any leftoever 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 which is not noise. It is not recommended to use the values below 4, as you will start erasing noise pixels.
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.
Here we use lower image resolution (by using the -rebin argument) and lower motion vector resolution (-Step switch with a fairly large argument, 1-2, versus the high-resolution values under 0.5) to make the search faster, and to fit the images into the GPU memory. (Only -rebin option can reduce the GPU memory consumption.) 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 each pixel being ~ 2 x FWHM. 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 other camera with smaller pixels (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 the factor of 3, which accelerates the computations and reduces the GPU memory consumption by a factor of 3x3=9 times.
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 the factor of (1.5/0.5)^2 = 9. There are no GPU memory changes when changing the Step value.
It is important to monitor what happens to asteroid detection when using -rebin and -Step options. As these parameters increase, the detection significance goes down. It is ok to loose ~5-10% of signal-to-noise ratio when doing the draft search, but not much more.
A preliminary search run can look like this:
asteroid_search -RMAX 0 -rebin 2 -Step 1.5 out_*fit
The option "-RMAX 0" means to search for the whole practical range of motion vectors, starting from the shortest vector (1 in FWHM units), 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.
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 to finish the run, at the beginning of the run, 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 clustrering 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 brightest pixel in 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.)
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.
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.
Now we need to measure the object's parameters more accurately, using the fine tuning option of asteroid_search.
Durung 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).
When doing 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 around the cloud 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 an copy these 4 values from the draft search stats.dat file for the cloud you are interetsed in.
The size of the cropped area is defined by the size of the square search area, both in image space, and in the motion vectors space. 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 -N_noise is to use the p_min value explicitely, using "-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 finutune search using more generous search windows - 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).
This is mostly on my to-do list. I will write a script to analyse 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 detections), or if they are likely be a noise.
Here is an example: in this plot, there are two real asteroid detections (high peaks), with the rest being the noise. (The higher peak is 3I/Atlas interstellar comet, imaged in July 2025). The horizontal axis is the length of the motion vector, the vertical axis is the pixel brightness.
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.