1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <assert.h>
7
8 #include "quad-builder.h"
9 #include "quad-utils.h"
10 #include "mathutil.h"
11 #include "errors.h"
12 #include "log.h"
13
14 struct quad {
15 unsigned int star[DQMAX];
16 };
17 typedef struct quad quad;
18
check_scale(quadbuilder_t * qb,pquad_t * pq)19 static void check_scale(quadbuilder_t* qb, pquad_t* pq) {
20 double *sA, *sB;
21 double s2;
22 double Bx=0, By=0;
23 double invscale;
24 double ABx, ABy;
25 Unused anbool ok;
26 if (!(qb->check_scale_low || qb->check_scale_high)) {
27 logverb("Not checking scale\n");
28 // ??
29 pq->scale_ok = TRUE;
30 return;
31 }
32 sA = qb->starxyz + pq->iA * 3;
33 sB = qb->starxyz + pq->iB * 3;
34 // s2: squared AB dist
35 s2 = distsq(sA, sB, 3);
36 pq->scale_ok = TRUE;
37 if (qb->check_scale_low && s2 < qb->quadd2_low)
38 pq->scale_ok = FALSE;
39 if (pq->scale_ok && qb->check_scale_high && s2 > qb->quadd2_high)
40 pq->scale_ok = FALSE;
41 if (!pq->scale_ok) {
42 qb->nbadscale++;
43 return;
44 }
45 star_midpoint(pq->midAB, sA, sB);
46 pq->scale_ok = TRUE;
47 pq->staridA = qb->starinds[pq->iA];
48 pq->staridB = qb->starinds[pq->iB];
49 ok = star_coords(sA, pq->midAB, TRUE, &pq->Ay, &pq->Ax);
50 assert(ok);
51 ok = star_coords(sB, pq->midAB, TRUE, &By, &Bx);
52 assert(ok);
53 ABx = Bx - pq->Ax;
54 ABy = By - pq->Ay;
55 invscale = 1.0 / (ABx*ABx + ABy*ABy);
56 pq->costheta = (ABy + ABx) * invscale;
57 pq->sintheta = (ABy - ABx) * invscale;
58 //nabok++;
59 }
60
61 static int
check_inbox(pquad_t * pq,int * inds,int ninds,double * stars)62 check_inbox(pquad_t* pq, int* inds, int ninds, double* stars) {
63 int i, ind;
64 double* starpos;
65 double Dx=0, Dy=0;
66 double ADx, ADy;
67 double x, y;
68 int destind = 0;
69 anbool ok;
70 for (i=0; i<ninds; i++) {
71 double r;
72 ind = inds[i];
73 starpos = stars + ind*3;
74 logverb("Star position: [%.4f, %.4f, %.4f]\n",
75 starpos[0], starpos[1], starpos[2]);
76 logverb("MidAB: [%.4f, %.4f, %.4f]\n",
77 pq->midAB[0], pq->midAB[1], pq->midAB[2]);
78
79 ok = star_coords(starpos, pq->midAB, TRUE, &Dy, &Dx);
80 if (!ok) {
81 logverb("star coords not ok\n");
82 continue;
83 }
84 ADx = Dx - pq->Ax;
85 ADy = Dy - pq->Ay;
86 x = ADx * pq->costheta + ADy * pq->sintheta;
87 y = -ADx * pq->sintheta + ADy * pq->costheta;
88 // make sure it's in the circle centered at (0.5, 0.5)...
89 // (x-1/2)^2 + (y-1/2)^2 <= r^2
90 // x^2-x+1/4 + y^2-y+1/4 <= (1/sqrt(2))^2
91 // x^2-x + y^2-y + 1/2 <= 1/2
92 // x^2-x + y^2-y <= 0
93 r = (x*x - x) + (y*y - y);
94 if (r > 0.0) {
95 logverb("star not in circle\n");
96 continue;
97 }
98 inds[destind] = ind;
99 destind++;
100 }
101 return destind;
102 }
103
104 /**
105 inbox, ninbox: the stars we have to work with.
106 starinds: the star identifiers (indexed by the contents of 'inbox')
107 - ie, starinds[inbox[0]] is an externally-recognized star identifier.
108 q: where we record the star identifiers
109 starnum: which star we're adding: eg, A=0, B=1, C=2, ... dimquads-1.
110 beginning: the first index in "inbox" to assign to star 'starnum'.
111 */
add_interior_stars(quadbuilder_t * qb,int ninbox,int * inbox,quad * q,int starnum,int dimquads,int beginning)112 static void add_interior_stars(quadbuilder_t* qb,
113 int ninbox, int* inbox, quad* q,
114 int starnum, int dimquads, int beginning) {
115 int i;
116 for (i=beginning; i<ninbox; i++) {
117 int iC = inbox[i];
118 q->star[starnum] = qb->starinds[iC];
119 // Did we just add the last star?
120 if (starnum == dimquads-1) {
121 if (qb->check_full_quad &&
122 !qb->check_full_quad(qb, q->star, dimquads, qb->check_full_quad_token))
123 continue;
124
125 qb->add_quad(qb, q->star, qb->add_quad_token);
126 } else {
127 if (qb->check_partial_quad &&
128 !qb->check_partial_quad(qb, q->star, starnum+1, qb->check_partial_quad_token))
129 continue;
130 // Recurse.
131 add_interior_stars(qb, ninbox, inbox, q, starnum+1, dimquads, i+1);
132 }
133 if (qb->stop_creating)
134 return;
135 }
136 }
137
quadbuilder_init()138 quadbuilder_t* quadbuilder_init() {
139 quadbuilder_t* qb = calloc(1, sizeof(quadbuilder_t));
140 return qb;
141 }
142
quadbuilder_free(quadbuilder_t * qb)143 void quadbuilder_free(quadbuilder_t* qb) {
144 free(qb->inbox);
145 free(qb->pquads);
146 free(qb);
147 }
148
quadbuilder_create(quadbuilder_t * qb)149 int quadbuilder_create(quadbuilder_t* qb) {
150 int iA=0, iB, iC, iD, newpoint;
151 int ninbox;
152 int i, j;
153 int iAalloc;
154 quad q;
155 pquad_t* qb_pquads;
156
157 // ensure the arrays are large enough...
158 if (qb->Nstars > qb->Ncq) {
159 // (free and malloc rather than realloc because we don't care about
160 // the previous contents)
161 free(qb->inbox);
162 free(qb->pquads);
163 qb->Ncq = qb->Nstars;
164 qb->inbox = calloc(qb->Nstars, sizeof(int));
165 qb->pquads = calloc(qb->Nstars * qb->Nstars, sizeof(pquad_t));
166 if (!qb->inbox || !qb->pquads) {
167 ERROR("quad-builder: failed to malloc qb->inbox or qb->pquads. Nstars=%i.\n", qb->Nstars);
168 return -1;
169 }
170 }
171
172 qb_pquads = qb->pquads;
173
174 /*
175 Each time through the "for" loop below, we consider a new
176 star ("newpoint"). First, we try building all quads that
177 have the new star on the diagonal (star B). Then, we try
178 building all quads that have the star not on the diagonal
179 (star D).
180
181 Note that we keep the invariants iA < iB and iC < iD.
182 */
183 memset(&q, 0, sizeof(quad));
184 for (newpoint=0; newpoint<qb->Nstars; newpoint++) {
185 pquad_t* pq;
186 logverb("Adding new star %i\n", newpoint);
187 // quads with the new star on the diagonal:
188 iB = newpoint;
189 for (iA = 0; iA < newpoint; iA++) {
190 pq = qb_pquads + iA*qb->Nstars + iB;
191 pq->inbox = NULL;
192 pq->ninbox = 0;
193 pq->iA = iA;
194 pq->iB = iB;
195
196 check_scale(qb, pq);
197 if (!pq->scale_ok) {
198 logverb("Dropping pair %i, %i based on scale\n", newpoint, iA);
199 continue;
200 }
201
202 q.star[0] = pq->staridA;
203 q.star[1] = pq->staridB;
204
205 pq->check_ok = TRUE;
206 if (qb->check_AB_stars)
207 pq->check_ok = qb->check_AB_stars(qb, pq, qb->check_AB_stars_token);
208 if (!pq->check_ok) {
209 logverb("Failed check for AB stars\n");
210 continue;
211 }
212
213 // list the possible internal stars...
214 ninbox = 0;
215 for (iC = 0; iC < newpoint; iC++) {
216 if ((iC == iA) || (iC == iB))
217 continue;
218 qb->inbox[ninbox] = iC;
219 ninbox++;
220 }
221
222 logverb("Number of possible internal stars for pair %i, %i: %i\n", newpoint, iA, ninbox);
223
224 // check which ones are inside the box...
225 ninbox = check_inbox(pq, qb->inbox, ninbox, qb->starxyz);
226 logverb("Number of stars in the box: %i\n", ninbox);
227
228 //if (!ninbox)
229 //continue;
230 if (ninbox && qb->check_internal_stars)
231 ninbox = qb->check_internal_stars(qb, q.star[0], q.star[1], qb->inbox, ninbox, qb->check_internal_stars_token);
232 //if (!ninbox)
233 //continue;
234
235 logverb("Number of stars in the box after checking: %i\n", ninbox);
236
237 add_interior_stars(qb, ninbox, qb->inbox, &q, 2, qb->dimquads, 0);
238 if (qb->stop_creating)
239 goto theend;
240
241 pq->inbox = malloc(qb->Nstars * sizeof(int));
242 if (!pq->inbox) {
243 ERROR("hpquads: failed to malloc pq->inbox.\n");
244 exit(-1);
245 }
246 pq->ninbox = ninbox;
247 memcpy(pq->inbox, qb->inbox, ninbox * sizeof(int));
248 debug("iA=%i, iB=%i: saved %i 'inbox' entries.\n", iA, iB, ninbox);
249 }
250 iAalloc = iA;
251
252 // quads with the new star not on the diagonal:
253 iD = newpoint;
254 for (iA = 0; iA < newpoint; iA++) {
255 for (iB = iA + 1; iB < newpoint; iB++) {
256 pq = qb_pquads + iA*qb->Nstars + iB;
257 if (!(pq->scale_ok && pq->check_ok))
258 continue;
259 // check if this new star is in the box.
260 qb->inbox[0] = iD;
261 ninbox = check_inbox(pq, qb->inbox, 1, qb->starxyz);
262 if (!ninbox)
263 continue;
264 if (qb->check_internal_stars)
265 ninbox = qb->check_internal_stars(qb, q.star[0], q.star[1], qb->inbox, ninbox, qb->check_internal_stars_token);
266 if (!ninbox)
267 continue;
268
269 pq->inbox[pq->ninbox] = iD;
270 pq->ninbox++;
271
272 q.star[0] = pq->staridA;
273 q.star[1] = pq->staridB;
274
275 add_interior_stars(qb, pq->ninbox, pq->inbox, &q, 2, qb->dimquads, 0);
276 if (qb->stop_creating) {
277 iA = iAalloc;
278 goto theend;
279 }
280 }
281 }
282 }
283 theend:
284 for (i=0; i<imin(qb->Nstars, newpoint+1); i++) {
285 int lim = (i == newpoint) ? iA : i;
286 for (j=0; j<lim; j++) {
287 pquad_t* pq = qb_pquads + j*qb->Nstars + i;
288 free(pq->inbox);
289 pq->inbox = NULL;
290 }
291 }
292 return 0;
293 }
294
295