Database generation taking too long

[CelestronPDEV]

The documentation mentions that new databases can be created from the Hipparcos and Tycho catalogs. I'm trying to make a new database from the Hipparcos catalog, like this:

import sys
sys.path.append('..')
from tetra3 import Tetra3
from PIL import Image
from pathlib import Path

# Create instance with no default_database
t3 = Tetra3()

# Generate and save HIP database
t3.generate_database(max_fov=3.0, star_catalog='hip_main',save_as='hip_database',star_max_magnitude=9.0)

However, this is taking a really, really long time. I am running this on an 8-core ARM M1 Mac with 16 GB of RAM. This is a fast machine. The code above has been running for almost a full day. Here is the output:

timmyd@ubuntu:~/tetra3/examples$ python3 make_databases.py 
2023-05-23 10:59:36,283:tetra3.Tetra3-INFO: Loading catalogue hip_main with 118218 star entries.
2023-05-23 10:59:37,362:tetra3.Tetra3-INFO: Skipped 263 incomplete entries.
2023-05-23 10:59:37,369:tetra3.Tetra3-INFO: Loaded 83337 stars with magnitude below 9.0.
2023-05-23 13:32:13,258:tetra3.Tetra3-INFO: For pattern matching with at most 10 stars per FOV and no doubles: 61140.
2023-05-23 13:32:13,259:tetra3.Tetra3-INFO: For verification with at most 20 stars per FOV and no doubles: 79865.
2023-05-23 13:32:13,544:tetra3.Tetra3-INFO: Generating all possible patterns.
2023-05-23 14:18:48,094:tetra3.Tetra3-INFO: Found 42955722 patterns. Building catalogue.
2023-05-24 08:49:40,924:tetra3.Tetra3-INFO: Finished generating database.
2023-05-24 08:49:40,925:tetra3.Tetra3-INFO: Size of uncompressed star table: 1916760 Bytes.
2023-05-24 08:49:40,925:tetra3.Tetra3-INFO: Size of uncompressed pattern catalog: 687291552 Bytes.
2023-05-24 08:49:40,926:tetra3.Tetra3-INFO: Saving database to: /home/timmyd/tetra3/examples/../hip_database.npz

The database generation takes almost a full day. Am I passing an obviously bad value for the max_fov parameter, or any other input?

Eventually we would like to create a database from the Tycho catalog (and later, maybe GAIA DR3 to magnitude 12). Those catalogs contain more than 1 million stars (and roughly 3 million for GAIA to mag 12.). We need this to work for fields of view no larger than about 1 degree, seeing stars to magnitude 12. Can the database generator work with requirements like these, and generate output in a reasonable amount of time (like, less than a week of runtime?)

[gustavmpettersson]

Hi Tim

Thanks for the message, I'll try to provide some useful info. In general, this code may not be ideal for your application, as the database required would be quite large. Is this intended for some specific application?

I'd suggest you use the Tycho catalogue if you are going to magnitude 9, the initial filtering step is fast anyway.

You can reduce the size of the database (and time to generate) by reducing the number of pattern matching stars and verification stars, but the reliability will reduce. Six matching and ten verification is around the minimum you could get away with.

I'm assuming you are on the master branch, but the catalogue_generator branch I worked on last year should be much faster to generate databases, and they should be compatible with master. That branch also has other significant speed improvements to the solutions, but I never had the time to polish it up. Give it a try if you haven't, it will still take hours.

As you reduce FOV, the number of stars and patterns required in the catalogue grows fast. I have not tested at small FOVs, it should all work so long as you have enough memory for the database and accept slower solutions. But millions of stars could be difficult with this type of solver.

The python implementation has not been optimised for database size, so there could be headroom to cut down too.

[CelestronPDEV]

Some answers in line...

I'd suggest you use the Tycho catalogue if you are going to magnitude 9, the initial filtering step is fast anyway.

That makes sense. I was starting with Hipparcos because it’s a smaller catalog. But it is incomplete past mag 7.5 or so. Eventually we are going to want a database to mag 12. For that, I will have to use Tycho-2 or (better) GAIA DR3.

You can reduce the size of the database (and time to generate) by reducing the number of pattern matching stars and verification stars, but the reliability will reduce. Six matching and ten verification is around the minimum you could get away with.

OK, great! I will try that.

I'm assuming you are on the master branch, but the catalogue_generator branch I worked on last year should be much faster to generate databases, and they should be compatible with master. That branch also has other significant speed improvements to the solutions, but I never had the time to polish it up. Give it a try if you haven't, it will still take hours.

I will try that too. Thanks for mentioning.

As you reduce FOV, the number of stars and patterns required in the catalogue grows fast. I have not tested at small FOVs, it should all work so long as you have enough memory for the database and accept slower solutions. But millions of stars could be difficult with this type of solver.

OK. I guess it’s only possible to find out by trying!

Here is another possible optimization. We know the field of view of our sensor (and hence image scale) to a very high degree of accuracy. Can this be used to reject many false Pattern matches? (or stated another way, can our known image scale let us find the correct Patterns more quickly?)

The python implementation has not been optimised for database size, so there could be headroom to cut down too.

OK. Last but not least … we are eventually going to need a C/C++ implementation. I was planning to translate the python code (once I can verify that it meets our requirements). But I wonder if somebody else already has a working C/C++ version. FYI, I did notice a C_Tetra set of sources files in Julian Brown’s original implementation (https://github.com/brownj4/Tetra). But it is incomplete; the code does not compile, and its catalog generator does not run successfully, either.

[gustavmpettersson]

The known FOV can be used to cull database matches very effectively, pass your known FOV as fov_estimate and your desired +/- tolerance as fov_max_error, both in degrees, to solve_from_image. This is the first check performed after all matching patterns have been unpacked from the database. Note that catalogue_generator should do this part a lot more efficiently (via bulk lookup in database instead of looping over them).

Also note you may want to consider changing the star_min_separation option for generate_database since you can probably distinguish closer doubles than the default 0.05 deg. Keep in mind that you may have a large magnitude difference.

Regarding implementation this set of code has evolved a lot since the original implementation by Julian, the core concept is the same but implementation is very different, there are not many original lines to be found.

Have fun!

[CelestronPDEV]

OK! Will do. Here is another possible optimization question. Often - but not always - we will know where our optics are looking, to a relatively high degree of accuracy (say, about 3-5 degrees). Can we use that knowledge to search a smaller subset of the Pattern database for matches? As opposed to searching an entire sky worth of Patterns?

[gustavmpettersson]

This is a pure lost-in-space algorithm, there is no benefit from prior knowledge. There are many other solvers made for updating a known state.

[CelestronPDEV]

OK, that makes sense. Can you recommend a particular solver for updating a known state that we might be able to use? If you have experience with an "update-solver" that works well, it would nicely complement tetra3 for our use-case.

[kuanjuliu]

@gustavmpettersson What are some solvers which can do such updates more quickly than Tetra3 can complete a lost-in-space search?

I found Tetra3 so much faster (orders of magnitude faster, in some cases) than, say, Astrometry’s solve-field I haven’t really looked further afield for more speedup.

[CelestronPDEV]

Here is an update. I am working on the catalogue_generator branch now. I successfully created a database using the Tycho-1 catalog of 1 million stars, to magnitude 12. (I know Tycho-1 doesn't go that faint, I just wanted to get as many stars as possible.) I did it as follows:

t3.generate_database(max_fov=1.0, star_catalog='tyc_main', pattern_stars_per_fov=6, verification_stars_per_fov=12, save_as='tyc12_database', star_max_magnitude=12.0, star_min_separation=.01)

Overnight, that generated a tyc12_database.npz:

2023-05-24 15:15:15,430:tetra3.Tetra3-INFO: Loading catalogue tyc_main with 1058332 star entries.
2023-05-24 15:15:26,074:tetra3.Tetra3-INFO: Skipped 22887 incomplete entries.
2023-05-24 15:15:26,184:tetra3.Tetra3-INFO: Loaded 1034887 stars with magnitude below 12.0.
2023-05-24 15:15:31,227:tetra3.Tetra3-INFO: Trimming database to requested star density.
2023-05-24 15:15:31,675:tetra3.Tetra3-INFO: Separating pattern stars by 0.44907311951024936deg.
2023-05-24 15:15:31,675:tetra3.Tetra3-INFO: Separating verification stars by 0.3175426480542942deg.
2023-05-24 15:15:47,712:tetra3.Tetra3-INFO: Stars for pattern matching: 311545.
2023-05-24 15:15:47,712:tetra3.Tetra3-INFO: Stars for verification: 485983.
2023-05-24 15:15:47,712:tetra3.Tetra3-INFO: Generating all possible patterns.
2023-05-24 15:33:44,072:tetra3.Tetra3-INFO: Found 38830390 patterns. Building catalogue.
2023-05-24 15:34:11,458:tetra3.Tetra3-INFO: Inserting pattern number: 1000000
...
2023-05-25 09:43:02,970:tetra3.Tetra3-INFO: Inserting pattern number: 38000000
2023-05-25 10:08:51,804:tetra3.Tetra3-INFO: Finished generating database.
2023-05-25 10:08:51,804:tetra3.Tetra3-INFO: Size of uncompressed star table: 11663592 Bytes.
2023-05-25 10:08:51,804:tetra3.Tetra3-INFO: Size of uncompressed pattern catalog: 621286240 Bytes.
2023-05-25 10:08:51,805:tetra3.Tetra3-INFO: Saving database to: tyc12_database.npz

I tried using that database to solve some of our real images on "easy" cases like M42 (the Orion Nebula) and M44 (the beehive star cluster) which contain a large number of Tycho-1 stars. In all cases, tetra3 failed to solve.

I ran each image through the get_centroids_from_image() function. I verified that with default inputs, that function is finding lots of star centroids at the correct, expected (x,y) coordinates.

I did not pass any fov_estimate to Tetra3.solve_from_image(). I do know our sensor's field of view very accurately (I cannot disclose it here due to the proprietary nature of the hardware) but it is a little bigger than a degree wide, and a little less than a degree tall.

Any suggestions for what we might be doing wrong? Gustav, could we hire you as a consultant for a couple days to get this all working? I'm quite confident it will work, now, and I think someone who knows the innards of the code will get us across the finish line faster. Please let me know your thoughts, and feel free to contact me privately at my personal email address timd@southernstars.com.

[kuanjuliu]

@CelestronTimD Candidly, I arrived at my Tetra3 settings by blasting through as many catalog builds and settings combinations as possible and looking for the “sweet spot”.

At least start by increasing match_threshold (say to 1) until you get some match — any match, no matter how poor the error — to see whether any of your existing settings even produce a result.

…

For me, I’ve always had to include fov_estimate to get good search results.

Because my images can have as few as five stars — mostly due to clouds and other obstructions like trees or even my house — I’ve found very distinct improvements by limiting catalog star magnitude to 8, rather than 9.

And in my case maybe because my star shapes aren’t great, using Hipparcos was better than Tycho. In contrast, platesolving much cleaner images off the web was better served by Tycho.

[gustavmpettersson]

I'm away from my computer now, can look more in detail tomorrow, but want to give some notes that may help immediately. I'm writing this from memory so apologies if there are errors.

Checking that you get reliable centroids is the correct first step.

Now, the database is culled to a certain star density, in you case to at least 0.45 degrees between stars. For 1x1 fov there won't be much left. It's quite possible that my suggestion of 6 pattern matching stars is too optimistic, this is something that has changed between versions, you might have to bump that back up to 10 if you want this to work reliably.

The suggestion to increase the probability threshold is very important and something I should have mentioned. With so few stars in the database the mismatch probability will never be low enough for the solution to be accepted with default threshold. I'm sorry there is no useful debugging in the solver code, but printing out all tested solutions and their probability and fov (just before the probability threshold is checked) should tell you if you get anything happening. If you there see something with the correct fov it's probably your match (and if you had more verification stars it should be accepted). Going back to 20 verification stars I think should not add too much time/size to the database, but may get your matches accepted. You can also increase the number of stars trialled for a pattern (don't recall the keyword) since many are missing in the database.

I appreciate the offer but I'm not available for paid work. I try to be available to respond to questions in my spare time, and could possibly have a quick zoom next week

[gustavmpettersson]

I pushed an update to catalogue_generator adding basic debug info for all tested matches, line 864.

The parameter for determining how many patterns are tried before giving up is pattern_checking_stars. Each pattern has 4 stars, so pattern_checking_stars choose 4 patterns are tried. By default 6 choose 4 = 15 patterns. If you have fewer patterns in the database, you may need to try more combinations of stars before "accidentally" finding one in the database, e.g. 8 (70 patterns).

[CelestronPDEV]

Gustav, I pulled your latest update to the catalogue_generator branch (that includes the new pattern_checking_stars option to solve_from_image()). I re-ran my test on a single image (of the Beehive cluster, M44, that contains lots of stars brighter than mag 10), as follows:

t3.solve_from_image(img,fov_estimate=1.257,pattern_checking_stars=10,fov_max_error=0.1,match_threshold=1.0e-2)

I increased the match_threshold_orobability from the default 1.0e-9 to 1.0e-2 (1 percent!) And also set patterh_checking_stars to 10. I believe these are the changes you suggested. Please note I am still using the same tyc12_database.npz that I generated yesterday (as described earlier in this thread.)

Tetra3 is still failing to solve:

Solution: {'RA': None, 'Dec': None, 'Roll': None, 'FOV': None, 'RMSE': None, 'Matches': None, 'Prob': None, 'T_solve': 3359.614941000018, 'T_extract': 152.92419100001098}

The tetra3 debug output does not show anything useful to me, but maybe this will provide some clues:

2023-05-26 08:37:09,632:tetra3.Tetra3-DEBUG: Tetra3 Constructor called with load_database=tyc12_database
2023-05-26 08:37:09,632:tetra3.Tetra3-DEBUG: Trying to load database
2023-05-26 08:37:09,632:tetra3.Tetra3-DEBUG: Got load database with: tyc12_database
2023-05-26 08:37:09,632:tetra3.Tetra3-DEBUG: String given, append to tetra3 directory
2023-05-26 08:37:09,632:tetra3.Tetra3-INFO: Loading database from: tyc12_database.npz
2023-05-26 08:37:09,638:tetra3.Tetra3-DEBUG: Loaded database, unpack files
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacking properties
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked pattern_mode to: edge_ratio
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked pattern_size to: 4
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked pattern_bins to: 25
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked pattern_max_error to: 0.005
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked max_fov to: 1.0
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked pattern_stars_per_fov to: 6
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked verification_stars_per_fov to: 12
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked star_max_magnitude to: 12.0
2023-05-26 08:37:11,841:tetra3.Tetra3-DEBUG: Unpacked star_min_separation to: 0.01
2023-05-26 08:37:15,640:tetra3.Tetra3-DEBUG: FAIL: Did not find a match to the stars! It took 3360 ms.

I am guessing that the problem is an inadequate database.

I will place my M44 image (and solution from astrometry.net) in a directory that you can download from.

Unfortunately, I am not able to upload the entire tyc12_database.npz (I'm currently on vacation at a remote site with very slow internet). Here is a possible workaround. Could you modify generate_database() to filter out stars more than X degrees from a specific (RA,Dec) point? That way, I could generate a much smaller pattern database of stars and quads within (say) 3 degrees of M44. This would probably allow a lot faster iteration for debugging purposes.

Let me know what you think. I will be back in civilization with fast internet on May 31st, but will check in from time to time here. Thanks again.

[gustavmpettersson]

Thanks for the note. I pushed the feature mentioned to create a partial database, will be very useful for testing when the databases get big. I created a database with (took only 30 seconds on my laptop): t3.generate_database(max_fov=1, star_catalog='tyc_main', save_as='test_m44', star_max_magnitude=12, star_min_separation=.01, range_ra=(125, 135), range_dec=(15, 25))

and ran on your image with: t3.solve_from_image(img, pattern_checking_stars=10)

Now, in the debug file we have:

2023-05-26 15:50:18,383:tetra3.Tetra3-DEBUG: Loaded database, unpack files
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacking properties
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked pattern_mode to: edge_ratio
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked pattern_size to: 4
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked pattern_bins to: 25
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked pattern_max_error to: 0.005
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked max_fov to: 1.0
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked pattern_stars_per_fov to: 10
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked verification_stars_per_fov to: 20
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked star_max_magnitude to: 12.0
2023-05-26 15:50:18,402:tetra3.Tetra3-DEBUG: Unpacked star_min_separation to: 0.01
2023-05-26 15:50:32,116:tetra3.Tetra3-DEBUG: FAIL: Did not find a match to the stars! It took 4 ms.
2023-05-26 15:51:51,894:tetra3.Tetra3-DEBUG: Possible match: Stars = 4, P_mismatch = 1.76e-03, FOV = 1.36894deg
2023-05-26 15:51:51,897:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25443deg
2023-05-26 15:51:51,904:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.05738deg
2023-05-26 15:51:51,915:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.76448deg
2023-05-26 15:51:51,916:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25443deg
2023-05-26 15:51:51,922:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25585deg
2023-05-26 15:51:51,929:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25443deg
2023-05-26 15:51:51,930:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25536deg
2023-05-26 15:51:51,938:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25443deg
2023-05-26 15:51:51,944:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25585deg
2023-05-26 15:51:51,950:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25443deg
2023-05-26 15:51:51,951:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25536deg
2023-05-26 15:51:51,954:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25585deg
2023-05-26 15:51:51,960:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25443deg
2023-05-26 15:51:51,961:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25405deg
2023-05-26 15:51:51,963:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25290deg
2023-05-26 15:51:51,966:tetra3.Tetra3-DEBUG: Possible match: Stars = 0, P_mismatch = 1.00e+00, FOV = 1.25536deg
2023-05-26 15:51:51,966:tetra3.Tetra3-DEBUG: FAIL: Did not find a match to the stars! It took 112 ms.

As you can see, we are now actually testing the correct solutions, as the FOV matches what you wrote above. Zero matches are found to verify the solution, which could be because there are not enough stars in the tycho catalogue for this small FOV.

[CelestronPDEV]

I got your latest update and I can duplicate your results. I am confused by your last sentence. There are many dozens of Tycho-1 stars in M44. Certainly many more than the 10 required to verify a pattern. I will experiment a bit more over the weekend. But I think it would be great to have a zoom call next week. Let me know generally what times are good for you. Thanks again. And have a great weekend!

[gustavmpettersson]

I think I've figured it out. The image you provided appears to be mirrored, if you add img = img.transpose(Image.FLIP_LEFT_RIGHT) it works perfectly. The solution is: {'RA': 129.95240146695699, 'Dec': 19.706620933962643, 'Roll': 314.00312457698504, 'FOV': 1.2544273167433113, 'RMSE': 1.1007408856833756, 'Matches': 20, 'Prob': 1.0157913068747216e-29, 'T_solve': 21.126100094988942, 'T_extract': 348.86879986152053}

I'll push an update to tetra3 where it will accept mirrored images too in a bit.

(For this, the catalogue was generated with 40 verification stars (separation of 0.17 degrees), to not filter out all the neighbours in M44, with: t3.generate_database(max_fov=1, verification_stars_per_fov=40, star_catalog='tyc_main', save_as='test_m44', star_max_magnitude=12, range_ra=(125, 135), range_dec=(15, 25), star_min_separation=.01)

[gustavmpettersson]

I have pushed a new update that allows mirrored images.

It also includes improvements to database generation, allowing multiscale databases (large FOV range, not relevant for you) and a new mode to simplify patterns.

The simplified patterns are such that each pattern is centred on a star with patterns created with all stars within FOV/2. The normal mode (remains the default) generates all patterns where the maximum distance between the stars is FOV. The simplification reduces the completeness of patterns but is much faster to generate, so it may be of interest for you.

Let me know how you get on with it

[CelestronPDEV]

Thanks! We are about to return home. I got your latest update and will do some experimentation on the flight home.

[CelestronPDEV]

Good news: with your latest commit, I can solve the M44 image. Cool!

When creating the database, the star_min_separation parameter is no longer recognized, so I just removed it. I assume this is deprecated?

I've begun translating your code to C++. I am able to read the .npz database with this library: https://github.com/rogersce/cnpy. That library reads the pattern_catalog and star_table arrays without any trouble, but it cannot read the props_packed table. That is a limitation of the C++ library (it cannot read numpy structured data) so will live without props_packed for now.

I understand the format of the star_table array completely. No questions there.

It looks like each row in the the pattern_catalog array is just a set of four integers defining the "quad" of stars which make up the pattern. And the integers are just row indexes into the star_table. Correct? if so, storing them as 16-bit integers may be a problem. Eventually the star_table will need to hold more than 65,536 (i.e. 2^16) entries. So I think the pattern_catalog should store the star indexes as 32-bit integers. That will avoid overflow problems. I think this will be needed when we eventually create a whole-sky database using the entire Tycho-2 (or GAIA) star catalogs to mag 12. Could you make this change, if my understanding of the pattern_catalog is correct?

Finally, I've noticed that quite a few rows in the pattern_catalog seem to be filled with all zeros. Is that expected?

[gustavmpettersson]

Great!

Correct, I noticed that star_min_separation doesn't do anything in this branch and removed it.

Correct, pattern_catalog just holds the indices in star_table of the four stars that make up the pattern. The data format for pattern_catalog is determined based on the largest index the table needs to hold, see rows 683-690 (for 99% of applications this turns out to be uint16). The shape and data type for pattern_catalog should be in the .npz header somewhere, but we could add it to props_packed too.

Yes, half of the pattern_catalog is empty to reduce collisions and therefore make lookup faster. You could experiment with this at line 682 and we could add an option to change the fill fraction if it makes a difference. All these disappear when the catalogue is compressed but they will use memory when loaded.

I experimented a bit building very large databases the other day and ran into memory limitations. I'd like to work on generating the database in chunks (of maybe 1 million patterns) saved gradually to disk, but will see when I get the time. Likewise, the option to not load the database into memory, but rather read directly from disk. Would be slower but maybe tolerable.

I assume you intend on implementing the solver only in C++, and will generate databases in Python still?

[CelestronPDEV]

Gustav, I now see the code to save the pattern_catalog as uint32 if the star table size requires it. I will modify my code to handle that (you are correct, I can get the data type from the table header.).

It does indeed feel wasteful for half the saved pattern_catalog to be empty. Perhaps this could be solved by creating a second lookup table of hashcode to index-within-pattern-catalog. That would eliminate all the empty pattern_catalog entries, and save a lot of memory, but still allow very fast pattern lookups from a hashcode?

You are also correct: for now, I only intend to port the solver to C++, and not even the star extractor (we already have one which we like in C++). The new C++ method I am writing will essentially by just "solve_from_centroids" instead of "solve_from_image". The database will still be generated in python using your code. I will probably also implement an option to use the database on disk without loading it into memory. Random file access is really fast in C++, especially with fixed-record-size file formats like we have here.

I'm sure I'll have more questions along the way :-). Thanks for being so responsive!

[CelestronPDEV]

Today's latest commit to the catalogue_generator branch looks like it's broken my ability to solve the M44 image. Here is debug output from today's code:

2023-05-31 09:00:37,495:tetra3.Tetra3-DEBUG: Tetra3 Constructor called with load_database=None
2023-05-31 09:00:37,495:tetra3.Tetra3-DEBUG: Got generate pattern catalogue with input: (1, None, 'test_m44', 'tyc_main', 10, 40, 12, 0.005, False, (125, 135), (15, 25))
2023-05-31 09:00:37,981:tetra3.Tetra3-INFO: Loading catalogue tyc_main with 1058332 star entries.
2023-05-31 09:00:44,198:tetra3.Tetra3-INFO: Skipped 22887 incomplete entries.
2023-05-31 09:00:44,298:tetra3.Tetra3-INFO: Loaded 1034887 stars with magnitude below 12.0.
2023-05-31 09:00:44,303:tetra3.Tetra3-INFO: Limited to RA range [125. 135.], keeping 36652 stars.
2023-05-31 09:00:44,303:tetra3.Tetra3-INFO: Limited to DEC range [15. 25.], keeping 1808 stars.
2023-05-31 09:00:44,312:tetra3.Tetra3-INFO: Trimming database to requested star density.
2023-05-31 09:00:44,313:tetra3.Tetra3-INFO: Generating patterns at FOV scales: [1.]
2023-05-31 09:00:44,313:tetra3.Tetra3-INFO: At FOV 1.0 separate stars by 0.18973665961010275deg.
2023-05-31 09:00:44,329:tetra3.Tetra3-INFO: Stars for creating patterns: 850.
2023-05-31 09:00:48,084:tetra3.Tetra3-INFO: Found 154300 patterns in total.
2023-05-31 09:00:48,097:tetra3.Tetra3-INFO: Total stars for verification: 1389.
2023-05-31 09:00:48,242:tetra3.Tetra3-INFO: Start building catalogue.
2023-05-31 09:00:48,318:tetra3.Tetra3-INFO: Catalog size (308600, 4) and type uint16.
2023-05-31 09:00:50,088:tetra3.Tetra3-INFO: Finished generating database.
2023-05-31 09:00:50,088:tetra3.Tetra3-INFO: Size of uncompressed star table: 33336 Bytes.
2023-05-31 09:00:50,088:tetra3.Tetra3-INFO: Size of uncompressed pattern catalog: 2468800 Bytes.
2023-05-31 09:00:50,088:tetra3.Tetra3-DEBUG: {'pattern_mode': 'edge_ratio', 'pattern_size': 4, 'pattern_bins': 50, 'pattern_max_error': 0.005, 'max_fov': 1.0, 'min_fov': 1.0, 'star_catalog': 'tyc_main', 'pattern_stars_per_fov': 10, 'verification_stars_per_fov': 40, 'star_max_magnitude': 12.0, 'simplify_pattern': False, 'range_ra': array([2.18166156, 2.35619449]), 'range_dec': array([0.26179939, 0.43633231])}
2023-05-31 09:00:50,088:tetra3.Tetra3-DEBUG: Saving generated database as: test_m44
2023-05-31 09:00:50,088:tetra3.Tetra3-DEBUG: Got save database with: test_m44
2023-05-31 09:00:50,088:tetra3.Tetra3-DEBUG: String given, append to tetra3 directory
2023-05-31 09:00:50,088:tetra3.Tetra3-INFO: Saving database to: /home/timmyd/tetra3/examples/../test_m44.npz
2023-05-31 09:00:50,090:tetra3.Tetra3-DEBUG: Packed properties into: ('edge_ratio', 4, 50, 0.005, 1., 1., 'tyc_main', 10, 40, 12., False, [2.1816616, 2.3561945], [0.2617994 , 0.43633232])
2023-05-31 09:00:50,090:tetra3.Tetra3-DEBUG: Saving as compressed numpy archive
2023-05-31 09:00:50,156:tetra3.Tetra3-DEBUG: Tetra3 Constructor called with load_database=test_m44
2023-05-31 09:00:50,156:tetra3.Tetra3-DEBUG: Trying to load database
2023-05-31 09:00:50,156:tetra3.Tetra3-DEBUG: Got load database with: test_m44
2023-05-31 09:00:50,156:tetra3.Tetra3-DEBUG: String given, append to tetra3 directory
2023-05-31 09:00:50,156:tetra3.Tetra3-INFO: Loading database from: /home/timmyd/tetra3/examples/../test_m44.npz
2023-05-31 09:00:50,157:tetra3.Tetra3-DEBUG: Loaded database, unpack files
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacking properties
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked pattern_mode to: edge_ratio
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked pattern_size to: 4
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked pattern_bins to: 50
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked pattern_max_error to: 0.005
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked max_fov to: 1.0
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked min_fov to: 1.0
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked star_catalog to: tyc_main
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked pattern_stars_per_fov to: 10
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked verification_stars_per_fov to: 40
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked star_max_magnitude to: 12.0
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked simplify_pattern to: False
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked range_ra to: [2.1816616 2.3561945]
2023-05-31 09:00:50,169:tetra3.Tetra3-DEBUG: Unpacked range_dec to: [0.2617994  0.43633232]
2023-05-31 09:00:50,176:tetra3.Tetra3-DEBUG: Got solve from image with input: (<PIL.PngImagePlugin.PngImageFile image mode=RGB size=3056x2048 at 0xFFFFAAD58BB0>, 1.257, 0.1, 20, 0.01, 0.01, {})
2023-05-31 09:00:50,483:tetra3.Tetra3-DEBUG: Extracted centroids:
[[ 392.29282  1118.23    ]
 [1033.9076   2066.8545  ]
 [ 856.3742   1814.4984  ]
 [1290.7      1832.5681  ]
 [ 419.9237   1215.9827  ]
 [ 771.974    2198.9465  ]
 [1412.3483   2364.7253  ]
 [ 943.74725  1362.287   ]
 [ 627.11816  1868.3536  ]
 [1367.1505   1297.486   ]
 [1230.4247   1849.4114  ]
 [1117.1029   2030.1577  ]
 [1882.7644    606.73175 ]
 [1781.636    2258.0884  ]
 [1241.1128   1878.0598  ]
 [ 762.1124     84.77428 ]
 [ 572.65375  2961.884   ]
 [1158.4412   2907.0107  ]
 [1484.0695   2703.4116  ]
 [ 422.1643    823.9361  ]
 [ 837.0323   1776.1843  ]
 [1030.7281    562.798   ]
 [ 567.3978   1534.2961  ]
 [ 837.3293    715.146   ]
 [ 365.6296   2311.72    ]
 [  74.191696 1715.2833  ]
 [1428.159      94.77195 ]
 [1489.8132   1465.2396  ]
 [1584.9039   2721.528   ]
 [ 168.14307  1580.7733  ]
 [ 353.69974  1201.1564  ]
 [ 193.64134   800.9888  ]
 [1796.5072   1477.8218  ]
 [ 251.46262  1857.1615  ]
 [1019.0182     97.40505 ]
 [ 345.14926   616.1117  ]
 [1274.8577    765.14197 ]
 [ 402.95642  1027.9496  ]
 [ 636.2208    729.4987  ]
 [1308.6831   1191.4276  ]]
2023-05-31 09:00:51,367:tetra3.Tetra3-DEBUG: FAIL: Did not find a match to the stars! It took 884 ms.

And here is debug output from yesterday's code, which solved the image successfully:

2023-05-30 14:51:14,255:tetra3.Tetra3-DEBUG: Tetra3 Constructor called with load_database=None
2023-05-30 14:51:14,256:tetra3.Tetra3-DEBUG: Got generate pattern catalogue with input: (1, None, 'test_m44', 'tyc_main', 10, 40, 12, 0.005, False, (125, 135), (15, 25))
2023-05-30 14:51:14,594:tetra3.Tetra3-INFO: Loading catalogue tyc_main with 1058332 star entries.
2023-05-30 14:51:20,603:tetra3.Tetra3-INFO: Skipped 22887 incomplete entries.
2023-05-30 14:51:20,696:tetra3.Tetra3-INFO: Loaded 1034887 stars with magnitude below 12.0.
2023-05-30 14:51:20,700:tetra3.Tetra3-INFO: Limited to RA range [125. 135.], keeping 36652 stars.
2023-05-30 14:51:20,701:tetra3.Tetra3-INFO: Limited to DEC range [15. 25.], keeping 1808 stars.
2023-05-30 14:51:20,710:tetra3.Tetra3-INFO: Trimming database to requested star density.
2023-05-30 14:51:20,710:tetra3.Tetra3-INFO: Generating patterns at FOV scales: [1.]
2023-05-30 14:51:20,710:tetra3.Tetra3-INFO: At FOV 1.0 separate stars by 0.18973665961010275deg.
2023-05-30 14:51:20,725:tetra3.Tetra3-INFO: Stars for creating patterns: 850.
2023-05-30 14:51:24,527:tetra3.Tetra3-INFO: Found 154300 patterns in total.
2023-05-30 14:51:24,541:tetra3.Tetra3-INFO: Total stars for verification: 1389.
2023-05-30 14:51:24,678:tetra3.Tetra3-INFO: Start building catalogue.
2023-05-30 14:51:24,753:tetra3.Tetra3-INFO: Catalog size (308600, 4) and type uint16.
2023-05-30 14:51:26,519:tetra3.Tetra3-INFO: Finished generating database.
2023-05-30 14:51:26,519:tetra3.Tetra3-INFO: Size of uncompressed star table: 33336 Bytes.
2023-05-30 14:51:26,519:tetra3.Tetra3-INFO: Size of uncompressed pattern catalog: 2468800 Bytes.
2023-05-30 14:51:26,519:tetra3.Tetra3-DEBUG: {'pattern_mode': 'edge_ratio', 'pattern_size': 4, 'pattern_bins': 50, 'pattern_max_error': 0.005, 'max_fov': 1.0, 'min_fov': 1.0, 'star_catalog': 'tyc_main', 'pattern_stars_per_fov': 10, 'verification_stars_per_fov': 40, 'star_max_magnitude': 12.0, 'simplify_pattern': False, 'range_ra': array([2.18166156, 2.35619449]), 'range_dec': array([0.26179939, 0.43633231])}
2023-05-30 14:51:26,519:tetra3.Tetra3-DEBUG: Saving generated database as: test_m44
2023-05-30 14:51:26,519:tetra3.Tetra3-DEBUG: Got save database with: test_m44
2023-05-30 14:51:26,519:tetra3.Tetra3-DEBUG: String given, append to tetra3 directory
2023-05-30 14:51:26,519:tetra3.Tetra3-INFO: Saving database to: /home/timmyd/tetra3/examples/../test_m44.npz
2023-05-30 14:51:26,519:tetra3.Tetra3-DEBUG: Packed properties into: ('edge_ratio', 4, 50, 0.005, 1., 1., 'tyc_main', 10, 40, 12., False, [2.1816616, 2.3561945], [0.2617994 , 0.43633232])
2023-05-30 14:51:26,519:tetra3.Tetra3-DEBUG: Saving as compressed numpy archive
2023-05-30 14:51:26,585:tetra3.Tetra3-DEBUG: Tetra3 Constructor called with load_database=test_m44
2023-05-30 14:51:26,585:tetra3.Tetra3-DEBUG: Trying to load database
2023-05-30 14:51:26,585:tetra3.Tetra3-DEBUG: Got load database with: test_m44
2023-05-30 14:51:26,585:tetra3.Tetra3-DEBUG: String given, append to tetra3 directory
2023-05-30 14:51:26,585:tetra3.Tetra3-INFO: Loading database from: /home/timmyd/tetra3/examples/../test_m44.npz
2023-05-30 14:51:26,585:tetra3.Tetra3-DEBUG: Loaded database, unpack files
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacking properties
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked pattern_mode to: edge_ratio
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked pattern_size to: 4
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked pattern_bins to: 50
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked pattern_max_error to: 0.005
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked max_fov to: 1.0
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked min_fov to: 1.0
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked star_catalog to: tyc_main
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked pattern_stars_per_fov to: 10
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked verification_stars_per_fov to: 40
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked star_max_magnitude to: 12.0
2023-05-30 14:51:26,597:tetra3.Tetra3-DEBUG: Unpacked simplify_pattern to: False
2023-05-30 14:51:26,598:tetra3.Tetra3-DEBUG: Unpacked range_ra to: [2.1816616 2.3561945]
2023-05-30 14:51:26,598:tetra3.Tetra3-DEBUG: Unpacked range_dec to: [0.2617994  0.43633232]
2023-05-30 14:51:26,601:tetra3.Tetra3-DEBUG: Got solve from image with input: (<PIL.PngImagePlugin.PngImageFile image mode=RGB size=3056x2048 at 0xFFFF7C833BB0>, 1.257, 0.1, 20, 0.01, 0.01, {})
2023-05-30 14:51:26,876:tetra3.Tetra3-DEBUG: Extracted centroids:
[[ 392.29282  1118.23    ]
 [1033.9076   2066.8545  ]
 [ 856.3742   1814.4984  ]
 [1290.7      1832.5681  ]
 [ 419.9237   1215.9827  ]
 [ 771.974    2198.9465  ]
 [1412.3483   2364.7253  ]
 [ 943.74725  1362.287   ]
 [ 627.11816  1868.3536  ]
 [1367.1505   1297.486   ]
 [1230.4247   1849.4114  ]
 [1117.1029   2030.1577  ]
 [1882.7644    606.73175 ]
 [1781.636    2258.0884  ]
 [1241.1128   1878.0598  ]
 [ 762.1124     84.77428 ]
 [ 572.65375  2961.884   ]
 [1158.4412   2907.0107  ]
 [1484.0695   2703.4116  ]
 [ 422.1643    823.9361  ]
 [ 837.0323   1776.1843  ]
 [1030.7281    562.798   ]
 [ 567.3978   1534.2961  ]
 [ 837.3293    715.146   ]
 [ 365.6296   2311.72    ]
 [  74.191696 1715.2833  ]
 [1428.159      94.77195 ]
 [1489.8132   1465.2396  ]
 [1584.9039   2721.528   ]
 [ 168.14307  1580.7733  ]
 [ 353.69974  1201.1564  ]
 [ 193.64134   800.9888  ]
 [1796.5072   1477.8218  ]
 [ 251.46262  1857.1615  ]
 [1019.0182     97.40505 ]
 [ 345.14926   616.1117  ]
 [1274.8577    765.14197 ]
 [ 402.95642  1027.9496  ]
 [ 636.2208    729.4987  ]
 [1308.6831   1191.4276  ]]
2023-05-30 14:51:26,888:tetra3.Tetra3-DEBUG: Possible match: Stars = 19, P_mismatch = 1.40e-30, FOV = 1.25441deg
2023-05-30 14:51:26,888:tetra3.Tetra3-DEBUG: MATCH ACCEPTED
2023-05-30 14:51:26,888:tetra3.Tetra3-DEBUG: Prob: 1.397e-30
2023-05-30 14:51:26,888:tetra3.Tetra3-DEBUG: RA:    129.95242146 deg
2023-05-30 14:51:26,888:tetra3.Tetra3-DEBUG: DEC:   19.70664516 deg
2023-05-30 14:51:26,888:tetra3.Tetra3-DEBUG: ROLL:  46.00072650 deg
2023-05-30 14:51:26,889:tetra3.Tetra3-DEBUG: FOV:   1.25441475 deg
2023-05-30 14:51:26,889:tetra3.Tetra3-DEBUG: MATCH: 19 stars
2023-05-30 14:51:26,889:tetra3.Tetra3-DEBUG: SOLVE: 12.16 ms
2023-05-30 14:51:26,889:tetra3.Tetra3-DEBUG: RESID: 1.11 asec

[gustavmpettersson]

Thanks for noticing, I've pushed a fix. I broke the FOV testing in the solver.

I think that making a separate table of lookups in the hash table without zero entries won't make the data any smaller. The hash table with uint16 mode has 8 bytes per row and 2N rows where N is number of patterns. Total size is 16N and it gets halved to 8N. If we create a separate table of 2N rows (for the hash indices) with uint32 mode (since there are millions of patterns) it has size 8N too. If you have enough entries that you need a uint32 table it could make a difference (saving 8N out of 32N bytes?) but I don't think it is that enticing. I'd rather suggest you try a table that is e.g. 1.3N long, that would bring the size down to ~10N (uint16) or ~20N (uint32) and needs no other changes. (Even better, pick next larger prime than 1.3N.)

Just a note on the mismatch probability: This is the false positive rate per test, but we can see it as if we're testing a single correct solution versus every possibility in the database. So if you want a 0.01 (1%) chance it's wrong, you need to set the testing level to 0.01/N. (If your big database is of N=100M patterns, you need 1e-10.) I'm considering adding this correction to the code, so a head's up it may change soon.

[CelestronPDEV]

I have not gotten much time to work on my C++ translation of the Tetra solver today. Hopefully, I should have more time tomorrow.

As a next step, I think I should create an "all sky" Tycho-1 database with the same star density and pattern density as we used for the M44 image ... just to get an idea of the overall database size, and the time to solve more images in different parts of the sky. Here are the settings I think I will use to create the all-sky database:

t3.generate_database(max_fov=1.0, verification_stars_per_fov=40, star_catalog='tyc_main', save_as='tyc12_allsky', star_max_magnitude=12)

Would you suggest increasing pattern_stars_per_fov from the default of 10? Should I also set the new simplify_pattern option to true?

[gustavmpettersson]

I'd suggest leaving pattern and verification stars at their default, they were not the issue before. In fact you may benefit from experimenting with making them smaller (especially pattern checking stars) depending on your tradeoff on database size and probability of success. I'd try simplify patterns too, and set the fov to close to the actual fov.

The size won't really be representative since there aren't enough stars in tycho, but it's a good next step to work out handling of very large databases.

You also may benefit from making the pattern max error smaller (makes less difference required for a pattern to get a different hash), that should speed up solving. 0.002 from the default 0.005 is something to try, it is limited by the distortion of your optics.

[CelestronPDEV]

Gustav, I am about halfway through my python-to-C++ translation of solve_from_centroids(), which is solve_from_image() that just takes a list of centroids as input (instead of the original image). So far my results agree with your results up to the inner loop here:

                valid_patterns = np.argwhere(max_edge_error < p_max_err)[:,0]
                # Go through each matching pattern and calculate further
                for index in valid_patterns:

I have not yet translated anything beyond that point.

Something which caught my eye is the method used to calculate angular distances between two vectors to two different stars in the sky. This is what you are using now, as far as I can tell:

angle = arccos ( dot_prouct ( v1, v2 ) )

(In other words, taking the inverse cosine of the dot product of the two vectors). This works fine if the two vectors point in fairly different directions (like, more than a few degrees apart.). But for really small angles, less than a degree, the inverse-cosine function becomes numerically inaccurate. A different formula, which is accurate over the entire range of zero to 180 degrees, is:

angle = 2.0 * arcsin ( distance ( v1, v2 ) / 2.0 )

in other words, twice the inverse sine of half the Cartesian distance between the vectors:

dx = v2.x - v1.x
dy = v2.y - v1.y
dz = v2.z - v1.z
distance = sqrt ( dx * dx + dy * dy + dz * dz )

It's slightly more computation, but is guaranteed to be accurate for all angles, no matter how small. I am noticing that the difference between your formula and mine is approaching 1 part in 1000 for some very small star patterns. It's enough to be worrisome.

For now, I am using your formula to stay compatible with your database. But this is something you might wish to address in the future as databases get larger (and star patterns get smaller).

[CelestronPDEV]

Just a quick note, my translation of the Tetra solver to C++ is almost complete. At least, it's complete enough that it solves the M44 image centroids now:

Loaded database with 13788 patterns and 1291 stars.
Read 136 sources from m44_centroids.txt
Solved!
RA: 129.95240784
Dec: 19.70663071
Roll: 46.01219177
FOV: 1.25435889
Matches: 10
Prob: 3.62108821e-14

I'm getting 10 matching stars using the same database.npz as the python code, which finds 12. I haven't figured out why my implementation is missing 2 stars, but it is still getting the same solution. My code needs a little cleanup, and I need to add the residuals and RMS errors. And the solution timing.

But I'm happy to share what I've got so far, if you're interested.

[gustavmpettersson]

Hi Tim,

Thanks on the note about small angle numerical error. This has not been an issue before, but I'll definitely swap out the implementation, just need to check if it impacts performance in any meaningful way.

Great that your implementation is running! Do you use your own centroiding method here? That may be the cause for slight differences. Note that the centroids passed should be ordered by brightness.

I'll also go ahead and create a solve_from_centroids method in the pyhton implementation. A simple addition to make the software more usable.

You have no obligation to with the license of tetra3, but it would be cool if you make a fork of tetra3 where you add in your solver. If it works well and runs faster maybe we'll even fold it in to the master release one day!

[CelestronPDEV]

Gustav,

I fixed the last issue in my code that was preventing me from finding those two extra stars. It wasn't the centroiding (I just copied the output of your centroid routine into a text file, and started with that.). Instead, it was the diagonal field-of-view radius calculation on line 947. That is the diameter of the FOV diagonal, not the radius. I missed where you were dividing the diameter by 2 on line 949. Now, when I do that, I get....

Loaded database with 13788 patterns and 1291 stars.
Read 136 sources from /Users/timmyd/Projects/SouthernStars/Projects/Tetra/m44_centroids.txt
Solved!
RA: 129.95240784
Dec: 19.70663261
Roll: 46.00732803
FOV: 1.25435889
Matches: 12
Prob: 1.67984294e-18
RMSE: 0.69416910 arcsec
TSolve: 2.874 ms

This exactly matches the python code output. And it seems to be about 3 times faster in C++. I consistently get solve times around 3 milliseconds, whereas the python code takes about 9 milliseconds, on the same machine (M1 MacBook Pro). I'm not a magician, it's just C++ vs python. But this is literally the only test case I've tried, so let's not celebrate quite yet.

A couple other things I've noticed.

You actually already had the second, higher-precision angle-difference formula in your code already, on line 986. This is where you're computing the residuals, and here the angles are indeed really, really small. Taking the arcsine of the magnitude of the cross product of the two vectors is almost identical to my formula. You're skipping the factor of 2 before the arcsine, and the factor of 1/2 on the inside. But for very small angles, the difference is negligible.

Finally, there are two other optimizations that might help speed up the solver. Both would require small mods to the database. In the solver, you calculate the largest angle of the catalog star patterns you are testing (lines 874-879). And later you calculate the pattern centroid vector, and sort the catalog pattern stars by distance-from-centroid (lines 915-920). This could all be done in the database generation. Save the pattern's star indices in sorted order, and save the catalog_largest_edge angle along with the pattern's indices. That angle is a floating point value in radians, but if you need to save it as an int32, you can convert it to milli-arcseconds and save the milli-arcseconds as an int32 - the precision will be good enough for all possible cases. This will avoid having to do the pattern angle calculation math (and sorting) in the solver.

Again, this stuff only really matters when you have a database of millions of stars and millions of patterns. But that's where we are heading, so I thought I would mention it.

Thanks again for sharing this amazing work with the world!

[timmyd7777]

Hi Gustav. I put the C++ port into a private GitHub repository. I need to keep it private for now. But I sent you an invitation from my personal GitHub account. Look for an invitation to collaborate on timmyd7777/Tetra3. Let me know if that invitation came through to you OK.

[KLiu-Arcadia]

@timmyd7777 I'm looking forward to your C++ port!

In the meantime: what is the speedup of the C++ vs Python, if any?

Tetra3 even six months ago was solving at >8.5 fps for 640 x 480 (ASI120MM-S @ 1280x960 binned 2x2) under actual sky conditions, making it feel real-time at least between slews.

I have yet to run the latest Tetra3 branch so I expect even greater speed (basically limited to lens/sensor speed) but if the C++ is even faster than that then I can start contemplating platesolving star trails instead of star points, making even quick manual slews real-time.

That would be mind-boggling.

[timmyd7777]

Hi KLiu. I ran a comparison of the python version vs. c++ version on a Raspberry Pi 4 with 4 GB of RAM. For both, I used the whole Tycho-1 database, and a set of 6 sample images with roughly 1-degree-wide fields of view. Resolutions are about 3K by 2K pixels.

Not every one of these images solve. The Tycho-1 database does not have enough stars. But this gives a sense of how long it takes images to fail (as well as succeed).

Keep in mind that the C++ version is not doing any star extraction, it's using the same centroids that were previously output by the python version. So you can only compare the python T_Solve times (versus the TSolve times in the C++ output).

I almost can't believe it myself, but the C++ version seems to be about 8x faster than the python version.

Here is raw output from the python version:

core@raspberrypi:~/tetra3/examples $ python3 test_origin_images.py 
2023-06-08 08:22:30,384:tetra3.Tetra3-INFO: Loading database from: ../tyc1_data.npz
Solving for image at: ../origin_images_all/Random3.png
Solution:
{'RA': None, 'Dec': None, 'Roll': None, 'FOV': None, 'RMSE': None, 'Matches': None, 'Prob': None, 'T_solve': 25156.638624000152, 'T_extract': 1794.9123629996393}
Solving for image at: ../origin_images_all/M44.png
Solution:
{'RA': 129.95240677854608, 'Dec': 19.70663064685016, 'Roll': 46.00732816566606, 'FOV': 1.2544147469552673, 'RMSE': 0.6963628574369269, 'Matches': 12, 'Prob': 1.679844087133935e-18, 'T_solve': 252.78112099931604, 'T_extract': 1810.0624680000692}
Solving for image at: ../origin_images_all/M66.png
Solution:
{'RA': None, 'Dec': None, 'Roll': None, 'FOV': None, 'RMSE': None, 'Matches': None, 'Prob': None, 'T_solve': 23355.938042999696, 'T_extract': 1720.3038810002909}
Solving for image at: ../origin_images_all/Random2.png
Solution:
{'RA': 140.15282715129635, 'Dec': -0.05241065264581809, 'Roll': 16.853378593932593, 'FOV': 1.2541711447813244, 'RMSE': 2.4811304177324223, 'Matches': 10, 'Prob': 1.1616076703095641e-14, 'T_solve': 2269.1403769995304, 'T_extract': 1762.1309719997953}
Solving for image at: ../origin_images_all/BullMono.png
Solution:
{'RA': None, 'Dec': None, 'Roll': None, 'FOV': None, 'RMSE': None, 'Matches': None, 'Prob': None, 'T_solve': 6833.120119000341, 'T_extract': 1502.9058629997962}
Solving for image at: ../origin_images_all/Random1.png
Solution:
{'RA': 154.5275568228282, 'Dec': 7.6905479807564525, 'Roll': 357.3968685429738, 'FOV': 1.2545079622928725, 'RMSE': 1.6427167104307747, 'Matches': 11, 'Prob': 7.053342420726643e-18, 'T_solve': 525.5379549998906, 'T_extract': 1797.7246720001858}
Solving for image at: ../origin_images_all/M42Mono.png
Solution:
{'RA': 84.2337572645459, 'Dec': -5.312096249165432, 'Roll': 3.906980630319052, 'FOV': 1.2550394536473033, 'RMSE': 1.1638122788864993, 'Matches': 13, 'Prob': 9.183840856943775e-21, 'T_solve': 501.87932300013927, 'T_extract': 1476.8457689997376}
Solving for image at: ../origin_images_all/Regulus.png
Solution:
{'RA': None, 'Dec': None, 'Roll': None, 'FOV': None, 'RMSE': None, 'Matches': None, 'Prob': None, 'T_solve': 21144.01852500032, 'T_extract': 1804.7082339999179}

Here is output from the C++ version. The images are processed in a different order than the python version, but results are otherwise identical:

core@raspberrypi:~/Tetra3 $ ./tetra3 ./TestData
Loaded database ./TestData/tyc1_data.npz
9310358 patterns and 634132 stars.
Read 25 sources from ./TestData/regulus_centroids.txt
Failed!
TSolve: 2969.886 ms

Read 36 sources from ./TestData/random1_centroids.txt
Solved!
RA: 154.52757263
Dec: 7.69054461
Roll: 357.39688110
FOV: 1.25441587
Matches: 11
Prob: 7.05333848e-18
RMSE: 1.62003243 arcsec
TSolve: 41.987 ms

Read 76 sources from ./TestData/m42_mono_centroids.txt
Solved!
RA: 84.23374176
Dec: -5.31209564
Roll: 3.90698147
FOV: 1.25493610
Matches: 13
Prob: 9.18383414e-21
RMSE: 1.14539719 arcsec
TSolve: 27.588 ms

Read 38 sources from ./TestData/m66_centroids.txt
Failed!
TSolve: 4071.261 ms

Read 41 sources from ./TestData/random2_centroids.txt
Solved!
RA: 140.15281677
Dec: -0.05240380
Roll: 16.85337830
FOV: 1.25430214
Matches: 10
Prob: 1.16160684e-14
RMSE: 2.44680357 arcsec
TSolve: 290.603 ms

Read 27 sources from ./TestData/random3_centroids.txt
Failed!
TSolve: 4370.096 ms

Read 18 sources from ./TestData/bull_mono_centroids.txt
Failed!
TSolve: 589.243 ms

Read 136 sources from ./TestData/m44_centroids.txt
Solved!
RA: 129.95240784
Dec: 19.70663261
Roll: 46.00732803
FOV: 1.25435889
Matches: 12
Prob: 1.67984294e-18
RMSE: 0.69416910 arcsec
TSolve: 32.149 ms

For M44, python solves in 252 ms; C++ solves in 32 ms. For Regulus, python fails 22144 ms, C++ fails in 2969 ms. Wow!

[gustavmpettersson]

I pushed a refactoring that now includes solve_from_centroids, the old solve_from_images is now little more than calling get_centroids_from_image and then solve_from_centroids on the result.

Nice to see the C++ implementation working well, and quite a bit faster is great!

I haven't looked at the optimisations you suggest yet but I'm not opposed to improving the speed. Could you take a moment to profile and analyse both database generation and the solver (maybe you already have)? I want to make sure I put effort into optimisations that actually matter.

[timmyd7777]

It's easier for me to change the C++ code at this point, than the python code. I made the first optimization on the C++ side (storing largest_edge angle for each star pattern in the database) when I read the database at program startup. Then use it later in the solver to immediately discard star patterns which are the wrong size.

That change alone speeds things up anywhere from 50% to 200%. Here is output from C++ solver, on the same Raspbery Pi 4, with the same test set of images. Compare this with yesterday's output:

core@raspberrypi:~/Tetra3 $ ./tetra3 ./TestData
Loaded database ./TestData/tyc1_data.npz
13965537 patterns and 634132 stars.
Read 25 sources from ./TestData/regulus_centroids.txt
Failed!
TSolve: 2048.167 ms

Read 36 sources from ./TestData/random1_centroids.txt
Solved!
RA: 154.52757263
Dec: 7.69054461
Roll: 357.39688110
FOV: 1.25441587
Matches: 11
Prob: 7.05333848e-18
RMSE: 1.62003243 arcsec
TSolve: 19.329 ms

Read 76 sources from ./TestData/m42_mono_centroids.txt
Solved!
RA: 84.23374176
Dec: -5.31209564
Roll: 3.90698147
FOV: 1.25493610
Matches: 13
Prob: 9.18383414e-21
RMSE: 1.14539719 arcsec
TSolve: 7.615 ms

Read 38 sources from ./TestData/m66_centroids.txt
Failed!
TSolve: 3029.438 ms

Read 41 sources from ./TestData/random2_centroids.txt
Solved!
RA: 140.15281677
Dec: -0.05240380
Roll: 16.85337830
FOV: 1.25430214
Matches: 10
Prob: 1.16160684e-14
RMSE: 2.44680357 arcsec
TSolve: 143.642 ms

Read 27 sources from ./TestData/random3_centroids.txt
Failed!
TSolve: 2963.487 ms

Read 18 sources from ./TestData/bull_mono_centroids.txt
Failed!
TSolve: 224.467 ms

Read 136 sources from ./TestData/m44_centroids.txt
Solved!
RA: 129.95240784
Dec: 19.70663261
Roll: 46.00732803
FOV: 1.25435889
Matches: 12
Prob: 1.67984294e-18
RMSE: 0.69416910 arcsec
TSolve: 12.603 ms

Regulus now fails in 2048 ms instead of 2969 ms. M44 solves in 12 ms instead of 32 ms.

This is now faster on the Raspberry Pi than the original Python code running on my 8-core M1 MacBook Pro! Holy cow.

If you want to experiment with this yourself, changes are checked into the Tetra3 C++ github repo. Look for the OPTIMIZE_DATABASE flag at the top of Tetra.hpp. Set that to 1 to turn on the database optimization or zero to turn it off (and then remember to recompile before you re-run). You can see the usages of OPTIMIZE_DATABASE to show how the changes affect the code.

[kuanjuliu]

@timmyd7777 Makes sense! Thanks for running those benchmarks.

What chance the star centroiding can also be sped up? That's the primary limiter on the RPi 4, anyway — the difference to, say, my M1 Max Macbook Pro is drastic in terms of relative runtimes between centroiding and solving.

Reason I ask is that centroiding has the potential to open a lot of interesting work, e.g. improving accuracy below pixel resolution, finding stars from star trails (think manual slews of telescopes), maybe even the tracking of satellites — but not if it takes as long as it does today.

[timmyd7777]

I guess the first step is translating the centroiding to C++. I have not even begun to look at that code yet, so I dunno.

On a multi-core machine, I would guess that making the centroiding multi-threaded should help a lot. Each thread can work on (say) 1/4 of the image at a time. On a 4-core machine, those threads can work in parallel, and - in theory - finish the whole image 4x faster.

Maybe the python version is already multi-threaded, I don't know.

[gustavmpettersson]

KJ what settings do you run for centroiding? Do you have a database of images we can use to tune the settings? Depending on what you need in terms of dimness, reliability etc there can be a lot of headroom from the default parameters.

Feel free to send stuff on the discord, Tim you can join here if you want a casual channel to discuss usage: https://discord.gg/8K9g7ePw

We should keep any bug or performance changes at GitHub though.

The centroiding method uses numpy functions for most of the heavy lifting so I'm not sure how much a 1:1 translation to C would bring

[timmyd7777]

OK. One last note here, then I will try to get onto Discord and take it from there.

I just made the second database optimization (pre-sorting star indices by distance from center of pattern, so the solver doesn't have to.) That made an even bigger difference, especially for the failure cases. Here are results (again running on the Raspberry Pi 4):

core@raspberrypi:~/Tetra3 $ ./tetra3 ./TestData
Loaded database ./TestData/tyc1_data.npz
9310358 patterns and 634132 stars.
Read 25 sources from ./TestData/regulus_centroids.txt
Failed!
TSolve: 671.260 ms

Read 36 sources from ./TestData/random1_centroids.txt
Solved!
RA: 154.52757263
Dec: 7.69054461
Roll: 357.39688110
FOV: 1.25441587
Matches: 11
Prob: 7.05333848e-18
RMSE: 1.62003243 arcsec
TSolve: 10.045 ms

Read 76 sources from ./TestData/m42_mono_centroids.txt
Solved!
RA: 84.23374176
Dec: -5.31209564
Roll: 3.90698147
FOV: 1.25493610
Matches: 13
Prob: 9.18383414e-21
RMSE: 1.14539719 arcsec
TSolve: 156.598 ms

Read 38 sources from ./TestData/m66_centroids.txt
Failed!
TSolve: 1434.101 ms

Read 41 sources from ./TestData/random2_centroids.txt
Solved!
RA: 140.15281677
Dec: -0.05240380
Roll: 16.85337830
FOV: 1.25430214
Matches: 10
Prob: 1.16160684e-14
RMSE: 2.44680357 arcsec
TSolve: 96.743 ms

Read 27 sources from ./TestData/random3_centroids.txt
Failed!
TSolve: 1342.136 ms

Read 18 sources from ./TestData/bull_mono_centroids.txt
Failed!
TSolve: 126.324 ms

Read 136 sources from ./TestData/m44_centroids.txt
Solved!
RA: 129.95240784
Dec: 19.70663261
Roll: 46.00732803
FOV: 1.25435889
Matches: 12
Prob: 1.67984294e-18
RMSE: 0.69416910 arcsec
TSolve: 7.832 ms

Regulus fails in 671 ms instead of 2048. M44 solves in 8 ms instead of 12.

Changes are checked into the C++ Tetra3 repo.

[timmyd7777]

One more big breakthrough. I was able to create a really big (224-MB) pattern database with the python3 version of Tetra3. This is based on a combined Hipparcos-2, Tycho-2, GAIA DR3 star catalog containing data for 2,875,310 stars. It is complete to magnitude 12. Of course GAIA goes much deeper than mag 12, but it's missing many of the brightest stars (so I combined it with Hipparcos and Tycho-2.)

With this big new 224 MB pattern database, both the C++ and python3 versions of Tetra3 can now solve all of my test images. As a way of saying thansk, I am providing instructions for how to recreate this work, if you would like it.

First download a ZIP file containing the combined HIP+TYC+GAIA star catalog from here. That catalog contains data for 2,875,310 stars in a CSV (comma-separated-value) format. The format is identical to other SSCore star data files, and is described here.

To build a Tetra3 pattern database from this, here's the python command:

t3 = Tetra3()
t3.generate_database(max_fov=1.0, star_catalog='HIP+TYC+GAIA.csv', save_as='hip_tyc_gaia_data', star_max_magnitude=12, simplify_pattern=True)

Oh! And you also need to add this snippet of python code to tetra3.py, so it can read the big star catalog:

        elif star_catalog == 'HIP+TYC+GAIA.csv':
            with open(catalog_file_full_pathname, 'r') as star_catalog_file:
                reader = csv.reader(star_catalog_file)
                for (i, entry) in enumerate(reader):  # star_num in range(num_entries):
                    hms = entry[1].split()
                    ra = np.deg2rad ( float ( hms[0] ) + float ( hms[1] ) / 60 + float ( hms[2] ) / 3600 ) * 15
                    dms = entry[2].split()
                    dec = np.deg2rad ( abs ( float ( dms[0] ) ) + float ( dms[1] ) / 60 + float ( dms[2] ) / 3600 )
                    if dms[0][0] == '-':
                        dec = -dec
                    if not entry[3]:
                        pmRA = 0
                    else:
                        pmRA = np.deg2rad ( float ( entry[3] ) / 3600 ) * 15
                    if not entry[4]:
                        pmDec = 0
                    else:
                        pmDec = np.deg2rad ( float ( entry[4] ) / 3600 )
                    if not entry[5]:
                        mag = star_max_magnitude + 1 # skip entries with no magnitude
                    else:
                        mag = float ( entry[5] )

                    ra += pmRA * ( current_year - 2000 )
                    dec += pmDec * ( current_year - 2000 )
                    #print ( hms, dms, ra, dec, pmRA, pmDec, mag ) 
                    star_table[i,:] = ([ra, dec, 0, 0, 0, mag])

This should go on line 544, right after the code to read the hip_main and tyc_main catalogs. I thought I would ask before pushing that change to your code.

[CelestronPDEV]

I am closing this issue. it's gotten far beyond the scope of the original question, and most importantly, the goals I was trying to achieve have been amply reached. Thanks for your help everyone!

[gustavmpettersson]

Hi guys,

Just leaving some final notes here for posterity. The points discussed here have all been pushed to master, including the considerable backlog of work that existed in the catalogue_generator branch. (The python version of) tetra3 should be about twice as fast and more flexible than the old master branch, which is still available as legacy.

Tim: Patterns are now sorted by default when generating a catalogue, this can be disabled by presort_patterns=False. There is now a key presort_patterns in Tetra3.database_properties indicating if they are sorted already (and for old databases this is automatically set to False). I also implemented saving the largest edge for faster FOV culling. On the python version these changes have negligible performance difference in my test cases (~10 deg) FOV, but it is probable that the performance gain is higher at smaller FOV since there will be (many) more patterns to reject in each hash bin.

By passing save_largest_edge=True a second array named Tetra3.pattern_largest_edge will be created, this saves the largest edge for each pattern (in the same index as the pattern_catalog) in milliradian as float16 type. I made a separate array since the datatype of pattern_catalog changes with the catalogue length, and think float16 is small and easily will meet all precision requirements for any database scale. This is different from your present solution, but handling the change should not be difficult I think. For old databases or when save_largest_edge=True is not passed, Tetra3.pattern_largest_edge returns None and the key pattern_largest_edge will be absent in the .npz archive.

I also made a test implementation using a numpy structured array to be able to combine different datatypes in one array for a cleaner implementation, but the performance was reduced since they are not able to be sliced and math-operated on in the same way as unstructured arrays.

Kind of unrelated thought, but it may be interesting to try shorter patterns (pattern_size = 3) for the small FOV application since it reduces the number of possible patterns dramatically. Of course, it also increases collisions. It's a tradeoff between database size and performance that I've never explored (and any pattern size other than 4 is not formally supported at this time).

KJ: I've improved the speed of the centroiding algorithm a little bit. Practically all compute time is now spent in numpy and scipy array operations, so there shouldn't be any performance on the table from a C/++ port of the code. In particular, when using binary_open = False the time is cut in about half, but you may have less clean output. Otherwise, a lot of speedup is possible by passing an explicit threshold level (image_th), if the application is pretty stable.

Thanks for your notes, they help making tetra3 better!

[kuanjuliu]

@gustavmpettersson Thanks! I look forward to testing the centroiding speedups — and taking out my own ham-handed attempts to use opencv to the same end.

[smroid]

Regarding faster star extraction and centroiding, I have a solution at https://github.com/smroid/StarGate that runs at better than 10ms per megapixel.

It is implemented in Rust but can be invoked from Python or other languages using gRPC. I've been using it in combination with Tetra3 with good results.