No solution with generated database

[astropatel]

Hello Gustav,

I have been trying to make this work, and have used your suggestions previously in generating a database that should find a solution, but I am at my wit's end. I was hoping you could point me in the right direction. I don't see your email on your profile so this was the only forum on which to contact you. (I've attached raw data in the dropbox link below.)

I have generated a fake 40x40 star field (full_catalog...fits) along with an associated datafile with star locations in pixel and ra,dec space and magnitudes. Using this data, I generated a database using the snippets of code:

But when I load the test fits file (onlythis_cutout_RAst59.49_DECst28.27.fits) data and pass it to solve_from_image, all the solutions lead to NaNs. I don't quite understand why a solution isn't being generated. Any assistance would be greatly appreciated.

Sincerely, Rahul

https://www.dropbox.com/scl/fo/qsemiisiwz8n02ixdxkw1/h?rlkey=v0trny3iv5m3c0lgj4ocohdsa&dl=0

[gustavmpettersson]

Hi Rahul,

You need to provide far more detail about what you are trying to do and provide all the code for how you are running tetra3. I understand you are making several modifications for an application yet to be explained.

What is this database? What modifications have you done to tetra3? Where are your fits files from? (They are fake? So how did you generate them?) How are you building the database? How are you getting the centroids? How are you calling tetra3?

As this does not appear to be an issue with tetra3 itself I will close this issue, please don't open any more. I can help in this thread as time permits.

[astropatel]

Apologies, Gustav for the lack of clarity. I'll start over.

I am trying to test out tetra3 using a simulated database of stars. I am randomly generating stars in a 40x40 deg FOV, using astropy to transform pixel locations to RA and DEC assuming some pixel scale. The database (.dat file in the dropbox folder) has the locations for all the point sources that I generated. I have also saved the associated FITS file showing all the stars that I generated (in the dropbox folder). This was all done using the astropy fits module. The centroids are generated using the following script (with astropy and photutils). They are just gaussian sources, with randomly generated pixel locations:

import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
import pandas as pd
from astropy import units as u
from astropy.io import fits
from photutils.datasets import make_gaussian_sources_image, make_random_gaussians_table
from astropy.nddata import Cutout2D

def stars(image, number, mag_min, mag_max):
    """
    Add some stars to the image.
    """

    # Most of the code below is a direct copy/paste from
    # https://photutils.readthedocs.io/en/stable/_modules/photutils/datasets/make.html#make_100gaussians_image

    flux_min = 10**mag_min
    flux_max = 10**mag_max
    flux_range = [flux_min, flux_max]

    y_max, x_max = image.shape
    xmean_range = [0.01 * x_max, 0.99 * x_max]
    ymean_range = [0.01 * y_max, 0.99 * y_max]
    xstddev_range = [4, 4]
    ystddev_range = [4, 4]
    params = dict([('flux', flux_range),
                  ('x_mean', xmean_range),
                  ('y_mean', ymean_range),
                  ('x_stddev', xstddev_range),
                  ('y_stddev', ystddev_range),
                  ('theta', [0, 2*np.pi])])

    sources = make_random_gaussians_table(number, params, seed=12345)

    star_im = make_gaussian_sources_image(image.shape, sources)

    return sources, star_im

The database I created was done using the following piece of code

https://user-images.githubusercontent.com/3648312/258617732-b608c691-61bf-429d-9560-ab9ef152bab6.png

with the attached version of Tetra3 I changed to include the use of a new database in generate_database

I then created a test image by taking a 10x10 deg cutout of the original 40x40 deg image to see if the new database would be able to identify the central pixel location in the test image using solve_from_image.

To run it, I"m using the following script. This works just fine using the test data and default database you provide in the repo, but provides NaNs for the database I created and the test image I use on it.::

import sys
sys.path.append('..')
from tetra3 import get_centroids_from_image
from tetra3 import Tetra3
from PIL import Image
from astropy.io import fits
from pathlib import Path

database = 'DP'

stuff = {'default_database':
             {'path':Path('../test_data/'),'ext':'*.tiff','database':'default_database'},
         'DP':
             {'path':Path('../../generated_star_catalog_files'), 'ext':'onlythis*cutout*.fits',
                  'database':Path('../hashtables/random_DP_catalog_default.npz')}
         }

t3 = Tetra3(stuff[database]['database'])

path = stuff[database]['path']
kw_args = {'return_images':True,'return_moments':True}

for impath in path.glob(stuff[database]['ext']):
    print('Solving for image at: ' + str(impath))
    if database == 'default_database':
        with Image.open(str(impath)) as img:
            # solved = t3.solve_from_image(img)
            solved = t3.solve_from_image(img)# Adding e.g. fov_estimate=11.4, fov_max_error=.1 improves performance
            centroids = get_centroids_from_image(img, **kw_args)
    else:
        hdu = fits.open(impath)[0].data.T
        solved = t3.solve_from_image(hdu)
        centroids = get_centroids_from_image(hdu, **kw_args)
    print('Solution: ' + str(solved))

I hope this information is sufficient. Please let me know if this information/explanation is incomplete.

tetra3_rahul.zip

[gustavmpettersson]

Thanks Rahul, I now understand what you are trying to do.

Your database generation script shows 10deg FOV but the debug print shows 40, I assume you have tried both?

Can you share the details of how you generate the star catalogue (ra/dec) from your randomised star positions?

Ideally, if you can put together the whole process from generating random image to running tetra3 in a script that I can run on my end, that would possibly be the fastest way forward.

[astropatel]

Hi Gustav,

Yes, I've tried both 10 and 40deg FOV for the max_fov parameter.

Below I've pasted the scripts you'll need to run it on your end.

0. Imports

   
   import numpy as np
   from pathlib import Path
   import matplotlib.pyplot as plt
   import pandas as pd
   from astropy.wcs import WCS
   from astropy import units as u
   from astropy.io import fits
   from astropy.visualization import AsinhStretch, LogStretch
   from astropy.visualization.mpl_normalize import ImageNormalize
   from convenience_functions import show_image
   from photutils.datasets import make_gaussian_sources_image, make_random_gaussians_table
   from astropy.nddata import Cutout2D

2. Generate random stars 

cwd = Path.cwd() fits_pixel_file = Path(cwd,'generated_star_catalog_files','full_catalog_base_pixelimage{}x{}deg.fits'.format(D_RA, D_DEC)) star_cat_file = Path(cwd,'generated_star_catalog_files','star_catalog_random_DP.dat')

RA_0, RA_F = 30, 70 DEC_0, DEC_F = 1,41 MAG_0, MAG_F = 3, 8 N_STARS = 250

D_RA = np.abs(RA_0 - RA_F) D_DEC = np.abs(DEC_0 - DEC_F)

for 40 deg x 40 deg FOV, use 20,000 pixels on a side

GRID_LENGTH = int(2e3) # N PIXELS PER SIDE

WCS STUFF

CRPIX = (int(GRID_LENGTH / 2), int(GRID_LENGTH / 2)) CRVAL = (RA_0 + D_RA / 2., DEC_0 + D_DEC / 2.)

cdelt_x = D_RA / GRID_LENGTH cdelt_y = D_DEC / GRID_LENGTH CDELT = (cdelt_x, cdelt_y) CTYPE = ["RA-TAN", "DEC-TAN"] RADESYS = 'FK5'

base_pix_full_image = np.zeros((GRID_LENGTH,GRID_LENGTH)) star_data, star_cat_img = stars(base_pix_full_image, N_STARS, MAG_0, MAG_F)

WCS

wcs_header = add_wcs(CRPIX, CRVAL, CDELT, CTYPE, RADESYS) star_catalog_data = star_data.to_pandas()

ra_cat, dec_cat = WCS(wcs_header).all_pix2world(star_catalog_data.x_mean, star_catalog_data.y_mean, 0) star_catalog_data['RA'] = ra_cat star_catalog_data['DEC'] = dec_cat star_catalog_data['MAG'] = np.log10(star_catalog_data['flux']) star_catalog_data['ID'] = np.uint16(star_catalog_data.index)

hdr = fits.Header() hdu_cat = fits.PrimaryHDU(header=wcs_header,data=star_cat_img)

WRITE FITS PIXEL FILE

hdu_cat.writeto(fits_pixel_file, overwrite=True) print('WRITING OUT NEW {}'.format(star_cat_file)) star_catalog_data.to_csv(star_cat_file,index=False)

2.5 Load Data

print('LOADING PRE-GENERATED "{}"'.format(fits_pixel_file)) hdu_cat = fits.open(fits_pixel_file) wcs_header = WCS(hdu_cat[0].header) star_cat_img = hdu_cat[0].data

print('LOADING PRE-GENERATED STAR DATA: {}'.format(star_cat_file)) star_catalog_data = pd.read_csv(star_cat_file) ra_cat = star_catalog_data['RA'] dec_cat = star_catalog_data['DEC'] mag_cat = star_catalog_data['MAG'] x_mean_cat = star_catalog_data['x_mean'] y_mean_cat = star_catalog_data['y_mean']

3. Create test image

roi_window_size = 10 ra_0_rand = np.random.uniform(RA_0, RA_F) dec_0_rand = np.random.uniform(DEC_0, DEC_F)

if ra_0_rand + roi_window_size > RA_F: ra_0_rand -= roi_window_size ra_f_rand = ra_0_rand + roi_window_size

if dec_0_rand + roi_window_size > DEC_F: dec_0_rand -= roi_window_size dec_f_rand = dec_0_rand + roi_window_size

Define the RA and Dec ranges of the region of interest (ROI)

ra_start, ra_end = ra_0_rand, ra_f_rand # RA range in degrees dec_start, dec_end = dec_0_rand, dec_f_rand # Dec range in degrees print("\n Random Box: RA:{:.2f}-{:.2f}, DEC:{:.2f}-{:.2f}".format(ra_start, ra_end, dec_start, dec_end))

Convert the RA and Dec ranges to pixel coordinates

x_start, y_start = wcs_header.all_world2pix(ra_start, dec_start, 0) x_end, y_end = wcs_header.all_world2pix(ra_end, dec_end, 0) xc = x_start + (x_end - x_start) 0.5 yc = y_start + (y_end - y_start) 0.5

Round the pixel coordinates to integers

x_start, x_end = int(np.round(x_start)), int(np.round(x_end)) y_start, y_end = int(np.round(y_start)), int(np.round(y_end)) roi_window_size_pix = x_end - x_start

ix_dat = np.where( (x_mean_cat>x_start) & (x_mean_cat<x_end) & (y_mean_cat>y_start) & (y_mean_cat<y_end))

CROP AND WRITE OUT CROPPED IMAGE FITS FILE WITH UPDATED WCS

hdu_cat_copy = hdu_cat[0].copy() hdu = hdu_cat_copy.data hdr = hdu_cat_copy.header

hdu_crop = Cutout2D(star_cat_img,position=(yc,xc),size=roi_window_size_pix,wcs=WCS(hdr))

wcs_cropped = hdu_crop.wcs

hdr.update(wcs_cropped.to_header()) hdr['COMMENT'] = "= Cropped fits file" fits_file_crop = Path(cwd,'generated_star_catalog_files', 'cutout_RAst{}_DECst{}.fits'.format(round(ra_start,2),round(dec_start,2))) fits.writeto(fits_file_crop, hdu_crop.data,hdr) print('Saving cropped fits file {}'.format(fits_file_crop))

WRITE OUT CROPPED STAR DATA

data_crop = star_catalog_data.loc[ix_dat].reset_index(drop=True) data_fits_file_crop = Path(cwd,'generated_star_catalog_files', 'cutout_RAst{}_DECst{}.csv'.format(round(ra_start,2),round(dec_start,2))) data_crop.to_csv(data_fits_file_crop, index=False)

4. Generate database/hashtable

import tetra3 t3 = tetra3.Tetra3(load_database=None) t3.generate_database(max_fov=10, save_as='hashtables/random_DP_catalog_default.npz', star_catalog='DP', star_max_magnitude=11 #MAG_F+3, pattern_max_error=0.005)

5. Use `solve_from_image`

import sys sys.path.append('..') from tetra3 import get_centroids_from_image from tetra3 import Tetra3 from PIL import Image from astropy.io import fits from pathlib import Path

database = 'DP'

stuff = {'default_database': {'path':Path('../test_data/'),'ext':'.tiff','database':'default_database'}, 'DP': {'path':Path('../../generated_star_catalog_files'), 'ext':'onlythiscutout*.fits', 'database':Path('../hashtables/random_DP_catalog_default.npz')} }

t3 = Tetra3(stuff[database]['database'])

path = stuff[database]['path'] kw_args = {'return_images':True,'return_moments':True}

for impath in path.glob(stuff[database]['ext']): print('Solving for image at: ' + str(impath)) if database == 'default_database': with Image.open(str(impath)) as img:

solved = t3.solve_from_image(img)

        solved = t3.solve_from_image(img)# Adding e.g. fov_estimate=11.4, fov_max_error=.1 improves performance
        centroids = get_centroids_from_image(img, **kw_args)
else:
    hdu = fits.open(impath)[0].data.T
    solved = t3.solve_from_image(hdu)
    centroids = get_centroids_from_image(hdu, **kw_args)
print('Solution: ' + str(solved))

[gustavmpettersson]

I've looked briefly but I'm not familiar with the wcs/fits ecosystem enough to give any immediate advice. It seems that tetra3 never tries the correct solution, as the number of matches is 3-4 at best. If you had the correct solution, then because this synthetic data I'd expect all centroids to be matched.

My inkling is that the projections from sky coordinates to camera image plane are not correct, and therefore the patterns are not as tetra3 assumes them. I would start with narrow FOV (say 5-20 degrees) centered at dec=0 to minimise the impact of these distortions, then you can set pattern_max_error = .02 or thereabouts when generating the database to relax the pattern accuracy required. I'd also just centre crop for the test images. Further, if you set max_fov to the size of your generated data and min_fov to a smaller value, you can try many crops with the same database.

[astropatel]

I tried those options but it still didn't work. I went ahead and did a line by line debugging to see what was going on. One issue I was able to address was that get_centroids_from_image was not producing any centroids, even though the diagnostic images showed clear centroids. What I realized was that max_area was smaller than the actual area calculated in calc_stats. I changed that, and that worked (set it to 300 pixels). But that still doesnt' produce any results.

But now, I found out that valid_patterns in solve_from_centroids is null (valid_patterns = np.argwhere(max_edge_error < p_max_err)[:, 0]) . I am trying to set p_max_err to 0.9 to in the database generation but that's taking forever to generate.

[gustavmpettersson]

Thanks for the note on get_centroids..., not showing the final result in any of the return images is a bug. They used to be in 'label_statistics', I think, that was last changed years ago. I've opened a separate issue to fix this soon.

pattern_max_error = 0.9 is way too large. All patterns will end up with the same hash index (0) and the algorithm grinds to a halt (as you've seen). The default is 0.005, you can try maybe up to 0.05.

[astropatel]

Thanks, Gustav.

At this point, I'm at my wit's end. I've verified the centroids extracted are correct from the test images, and confirmed the positions of the extracted centroids and as far as I can tell this should work. I really don't know what's keeping this from working.