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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
|
/*
# This file is part of the Astrometry.net suite.
# Licensed under a 3-clause BSD style license - see LICENSE
*/
#include <assert.h>
#include "quad-utils.h"
#include "starutil.h"
#include "codefile.h"
#include "starkd.h"
#include "errors.h"
#include "log.h"
void quad_compute_star_code(const double* starxyz, double* code, int dimquads) {
double Ax=0, Ay=0;
double Bx=0, By=0;
double ABx, ABy;
double scale, invscale;
double costheta, sintheta;
double midAB[3];
Unused anbool ok;
int i;
const double *sA, *sB;
sA = starxyz;
sB = starxyz + 3;
star_midpoint(midAB, sA, sB);
ok = star_coords(sA, midAB, TRUE, &Ay, &Ax);
assert(ok);
ok = star_coords(sB, midAB, TRUE, &By, &Bx);
assert(ok);
ABx = Bx - Ax;
ABy = By - Ay;
scale = (ABx * ABx) + (ABy * ABy);
invscale = 1.0 / scale;
costheta = (ABy + ABx) * invscale;
sintheta = (ABy - ABx) * invscale;
for (i=2; i<dimquads; i++) {
const double* starpos;
double Dx=0, Dy=0;
double ADx, ADy;
double x, y;
starpos = starxyz + 3*i;
ok = star_coords(starpos, midAB, TRUE, &Dy, &Dx);
assert(ok);
ADx = Dx - Ax;
ADy = Dy - Ay;
x = ADx * costheta + ADy * sintheta;
y = -ADx * sintheta + ADy * costheta;
code[2*(i-2)+0] = x;
code[2*(i-2)+1] = y;
}
}
void quad_flip_parity(const double* code, double* flipcode, int dimcode) {
int i;
// swap CX <-> CY, DX <-> DY.
for (i=0; i<dimcode/2; i++) {
// use tmp in code "code" == "flipcode"
double tmp;
tmp = code[2*i+1];
flipcode[2*i+1] = code[2*i+0];
flipcode[2*i+0] = tmp;
}
}
int quad_compute_code(const unsigned int* quad, int dimquads, startree_t* starkd,
double* code) {
int i;
double starxyz[3 * DQMAX];
for (i=0; i<dimquads; i++) {
if (startree_get(starkd, quad[i], starxyz + 3*i)) {
ERROR("Failed to get stars belonging to a quad.\n");
return -1;
}
}
quad_compute_star_code(starxyz, code, dimquads);
return 0;
}
anbool quad_obeys_invariants(unsigned int* quad, double* code,
int dimquads, int dimcodes) {
double sum;
int i;
// check the invariant that (cx + dx + ...) / (dimquads-2) <= 1/2
sum = 0.0;
for (i=0; i<(dimquads-2); i++)
sum += code[2*i];
sum /= (dimquads-2);
if (sum > 0.5)
return FALSE;
// check the invariant that cx <= dx <= ....
for (i=0; i<(dimquads-3); i++)
if (code[2*i] > code[2*(i+1)])
return FALSE;
return TRUE;
}
void quad_enforce_invariants(unsigned int* quad, double* code,
int dimquads, int dimcodes) {
double sum;
int i;
// here we add the invariant that (cx + dx + ...) / (dimquads-2) <= 1/2
sum = 0.0;
for (i=0; i<dimcodes/2; i++)
sum += code[2*i];
sum /= (dimcodes/2);
if (sum > 0.5) {
logdebug("Flipping code to ensure mean(x)<=0.5\n");
// swap the labels of A,B.
int tmp = quad[0];
quad[0] = quad[1];
quad[1] = tmp;
// rotate the code 180 degrees.
for (i=0; i<dimcodes; i++)
code[i] = 1.0 - code[i];
}
// here we add the invariant that cx <= dx <= ....
for (i=0; i<(dimquads-2); i++) {
int j;
int jsmallest;
double smallest;
double x1;
double dtmp;
int tmp;
x1 = code[2*i];
jsmallest = -1;
smallest = x1;
for (j=i+1; j<(dimquads-2); j++) {
double x2 = code[2*j];
if (x2 < smallest) {
smallest = x2;
jsmallest = j;
}
}
if (jsmallest == -1)
continue;
j = jsmallest;
// swap the labels.
tmp = quad[i+2];
quad[i+2] = quad[j+2];
quad[j+2] = tmp;
// swap the code values.
dtmp = code[2*i];
code[2*i] = code[2*j];
code[2*j] = dtmp;
dtmp = code[2*i+1];
code[2*i+1] = code[2*j+1];
code[2*j+1] = dtmp;
}
}
|
