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