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-utils.h"
9 #include "starutil.h"
10 #include "codefile.h"
11 #include "starkd.h"
12 #include "errors.h"
13 #include "log.h"
14
quad_compute_star_code(const double * starxyz,double * code,int dimquads)15 void quad_compute_star_code(const double* starxyz, double* code, int dimquads) {
16 double Ax=0, Ay=0;
17 double Bx=0, By=0;
18 double ABx, ABy;
19 double scale, invscale;
20 double costheta, sintheta;
21 double midAB[3];
22 Unused anbool ok;
23 int i;
24 const double *sA, *sB;
25
26 sA = starxyz;
27 sB = starxyz + 3;
28 star_midpoint(midAB, sA, sB);
29 ok = star_coords(sA, midAB, TRUE, &Ay, &Ax);
30 assert(ok);
31 ok = star_coords(sB, midAB, TRUE, &By, &Bx);
32 assert(ok);
33 ABx = Bx - Ax;
34 ABy = By - Ay;
35 scale = (ABx * ABx) + (ABy * ABy);
36 invscale = 1.0 / scale;
37 costheta = (ABy + ABx) * invscale;
38 sintheta = (ABy - ABx) * invscale;
39
40 for (i=2; i<dimquads; i++) {
41 const double* starpos;
42 double Dx=0, Dy=0;
43 double ADx, ADy;
44 double x, y;
45 starpos = starxyz + 3*i;
46 ok = star_coords(starpos, midAB, TRUE, &Dy, &Dx);
47 assert(ok);
48 ADx = Dx - Ax;
49 ADy = Dy - Ay;
50 x = ADx * costheta + ADy * sintheta;
51 y = -ADx * sintheta + ADy * costheta;
52 code[2*(i-2)+0] = x;
53 code[2*(i-2)+1] = y;
54 }
55 }
56
quad_flip_parity(const double * code,double * flipcode,int dimcode)57 void quad_flip_parity(const double* code, double* flipcode, int dimcode) {
58 int i;
59 // swap CX <-> CY, DX <-> DY.
60 for (i=0; i<dimcode/2; i++) {
61 // use tmp in code "code" == "flipcode"
62 double tmp;
63 tmp = code[2*i+1];
64 flipcode[2*i+1] = code[2*i+0];
65 flipcode[2*i+0] = tmp;
66 }
67 }
68
quad_compute_code(const unsigned int * quad,int dimquads,startree_t * starkd,double * code)69 int quad_compute_code(const unsigned int* quad, int dimquads, startree_t* starkd,
70 double* code) {
71 int i;
72 double starxyz[3 * DQMAX];
73 for (i=0; i<dimquads; i++) {
74 if (startree_get(starkd, quad[i], starxyz + 3*i)) {
75 ERROR("Failed to get stars belonging to a quad.\n");
76 return -1;
77 }
78 }
79 quad_compute_star_code(starxyz, code, dimquads);
80 return 0;
81 }
82
quad_obeys_invariants(unsigned int * quad,double * code,int dimquads,int dimcodes)83 anbool quad_obeys_invariants(unsigned int* quad, double* code,
84 int dimquads, int dimcodes) {
85 double sum;
86 int i;
87 // check the invariant that (cx + dx + ...) / (dimquads-2) <= 1/2
88 sum = 0.0;
89 for (i=0; i<(dimquads-2); i++)
90 sum += code[2*i];
91 sum /= (dimquads-2);
92 if (sum > 0.5)
93 return FALSE;
94
95 // check the invariant that cx <= dx <= ....
96 for (i=0; i<(dimquads-3); i++)
97 if (code[2*i] > code[2*(i+1)])
98 return FALSE;
99 return TRUE;
100 }
101
quad_enforce_invariants(unsigned int * quad,double * code,int dimquads,int dimcodes)102 void quad_enforce_invariants(unsigned int* quad, double* code,
103 int dimquads, int dimcodes) {
104 double sum;
105 int i;
106
107 // here we add the invariant that (cx + dx + ...) / (dimquads-2) <= 1/2
108 sum = 0.0;
109 for (i=0; i<dimcodes/2; i++)
110 sum += code[2*i];
111 sum /= (dimcodes/2);
112 if (sum > 0.5) {
113 logdebug("Flipping code to ensure mean(x)<=0.5\n");
114 // swap the labels of A,B.
115 int tmp = quad[0];
116 quad[0] = quad[1];
117 quad[1] = tmp;
118 // rotate the code 180 degrees.
119 for (i=0; i<dimcodes; i++)
120 code[i] = 1.0 - code[i];
121 }
122
123 // here we add the invariant that cx <= dx <= ....
124 for (i=0; i<(dimquads-2); i++) {
125 int j;
126 int jsmallest;
127 double smallest;
128 double x1;
129 double dtmp;
130 int tmp;
131
132 x1 = code[2*i];
133 jsmallest = -1;
134 smallest = x1;
135 for (j=i+1; j<(dimquads-2); j++) {
136 double x2 = code[2*j];
137 if (x2 < smallest) {
138 smallest = x2;
139 jsmallest = j;
140 }
141 }
142 if (jsmallest == -1)
143 continue;
144 j = jsmallest;
145 // swap the labels.
146 tmp = quad[i+2];
147 quad[i+2] = quad[j+2];
148 quad[j+2] = tmp;
149 // swap the code values.
150 dtmp = code[2*i];
151 code[2*i] = code[2*j];
152 code[2*j] = dtmp;
153 dtmp = code[2*i+1];
154 code[2*i+1] = code[2*j+1];
155 code[2*j+1] = dtmp;
156 }
157 }
158
159