1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
/*
# This file is part of the Astrometry.net suite.
# Licensed under a 3-clause BSD style license - see LICENSE
*/
#include <stdio.h>
#include <assert.h>
#include "cutest.h"
#include "index.h"
#include "solver.h"
#include "xylist.h"
#include "sip.h"
#include "sip-utils.h"
#include "bl.h"
#include "log.h"
#include "errors.h"
#include "sip_qfits.h"
/*
- example image: DESI Commissioning Instrument exposure 3367-CIW.
*/
void test_xscale(CuTest* ct) {
// core Astrometry solver parameters
solver_t* solver;
int imagew, imageh;
double imagecx, imagecy;
double arcmin_width_min = 5;
double arcmin_width_max = 8;
char* xyfn = "3367-W.xy";
//char* xyfn = "3367-W-110.axy";
char* indexfn = "/data1/INDEXES/5000/index-5001-31.fits";
int loglvl = LOG_MSG;
loglvl++;
log_init(loglvl);
imagew = 3072;
imageh = 2048;
imagecx = (imagew - 1.0)/2.0;
imagecy = (imageh - 1.0)/2.0;
// Here we initialize the core astrometry solver struct, telling
// it about the possible range of image scales.
solver = solver_new();
double qsf_min = 0.1;
// don't try teeny-tiny quads.
solver->quadsize_min = qsf_min * MIN(imagew, imageh);
// compute scale range in arcseconds per pixel.
// set the solver's "funits" = field (image) scale units
solver->funits_lower = 60. * arcmin_width_min / (double)imagew;
solver->funits_upper = 60. * arcmin_width_max / (double)imagew;
solver_set_keep_logodds(solver, log(1e12));
solver->logratio_toprint = log(1e6);
solver->endobj = 20;
xylist_t* xyls = xylist_open(xyfn);
starxy_t* xy = xylist_read_field(xyls, NULL);
// Feed the image source coordinates to the solver...
solver_set_field(solver, xy);
solver_set_field_bounds(solver, 0, imagew, 0, imageh);
index_t* index = index_load(indexfn, 0, NULL);
solver_add_index(solver, index);
solver->distance_from_quad_bonus = TRUE;
solver->do_tweak = TRUE;
solver->tweak_aborder = 1;
solver->tweak_abporder = 4;
solver_run(solver);
CuAssert(ct, "Should not solve on undistorted field", !solver->best_match_solves);
// "solver" will free the original "xy", so make a copy.
xy = starxy_copy(xy);
solver_set_field(solver, xy);
solver_reset_best_match(solver);
solver_reset_counters(solver);
solver->pixel_xscale = 1.1;
// *1.1?
solver_set_field_bounds(solver, 0, imagew, 0, imageh);
solver_print_to(solver, stdout);
solver_run(solver);
CuAssert(ct, "Should succeeded on scaled field", solver->best_match_solves);
double ra, dec;
double pscale;
tan_t* wcs;
logmsg("Solved using index %s with odds ratio %g\n",
solver->best_index->indexname,
solver->best_match.logodds);
// WCS is solver->best_match.wcstan
wcs = &(solver->best_match.wcstan);
// center
tan_pixelxy2radec(wcs, imagecx, imagecy, &ra, &dec);
pscale = tan_pixel_scale(wcs);
logmsg("Image center is RA,Dec = (%g,%g) degrees, size is %.2g x %.2g arcmin.\n",
ra, dec, arcsec2arcmin(pscale * imagew), arcsec2arcmin(pscale * imageh));
printf("Writing SIP solution...\n");
//solver->best_match.wcstan;
sip_write_to_file(solver->best_match.sip, "scaled-sip.wcs");
}
|
