Instrument Neutral Distributed Interface INDI  2.0.2
align.c
Go to the documentation of this file.
1 
2 /*
3 * DSP API - a digital signal processing library for astronomy usage
4 * Copyright © 2017-2021 Ilia Platone
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU Lesser General Public
8 * License as published by the Free Software Foundation; either
9 * version 3 of the License, or (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with this program; if not, write to the Free Software Foundation,
18 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20 
21 #include "dsp.h"
22 #include <time.h>
23 
24 int dsp_qsort_double_desc (const void *arg1, const void *arg2)
25 {
26  double* a1 = (double*)arg1;
27  double* a2 = (double*)arg2;
28  return ((*a1) < (*a2) ? 1 : -1);
29 }
30 
31 int dsp_qsort_double_asc (const void *arg1, const void *arg2)
32 {
33  double* a1 = (double*)arg1;
34  double* a2 = (double*)arg2;
35  return ((*a1) > (*a2) ? 1 : -1);
36 }
37 
38 int dsp_qsort_star_diameter_desc(const void *arg1, const void *arg2)
39 {
40  dsp_star* a = (dsp_star*)arg1;
41  dsp_star* b = (dsp_star*)arg2;
42  if(a->diameter < b->diameter)
43  return 1;
44  return -1;
45 }
46 
47 int dsp_qsort_star_diameter_asc(const void *arg1, const void *arg2)
48 {
49  dsp_star* a = (dsp_star*)arg1;
50  dsp_star* b = (dsp_star*)arg2;
51  if(a->diameter > b->diameter)
52  return 1;
53  return -1;
54 }
55 
56 static double calc_match_score(dsp_triangle t1, dsp_triangle t2, dsp_align_info align_info)
57 {
58  int d = 0;
59  double err = fabs((t1.index - t2.index)/t1.index);
60  double size1 = t1.sizes[0];
61  double size2 = t2.sizes[0]/align_info.factor[0];
62  err += fabs((size1 - size2)/size1);
63  for(d = 1; d < align_info.dims; d ++) {
64  double size1 = t1.sizes[d-1];
65  double size2 = t2.sizes[d-1]/align_info.factor[d-1];
66  double x1 = (t1.stars[d].center.location[d-1]-t1.stars[d-1].center.location[d-1]);
67  double y1 = (t1.stars[d].center.location[d]-t1.stars[d-1].center.location[d]);
68  double _x2 = (t2.stars[d].center.location[d-1]-align_info.offset[d-1]-align_info.center[d-1]);
69  double _y2 = (t2.stars[d].center.location[d]-align_info.offset[d]-align_info.center[d]);
70  double x2 = _x2*cos(align_info.radians[d-1])+_y2*sin(align_info.radians[d-1]);
71  double y2 = _y2*cos(align_info.radians[d-1])-_x2*sin(align_info.radians[d-1]);
72  x2 /= align_info.factor[d-1];
73  y2 /= align_info.factor[d];
74  err += fabs(((x2-x1)+(y2-y1)+(size2-size1))/size1);
75  }
76  for(d = 1 ; d < t1.dims; d++) {
77  err += fabs((t1.ratios[d] - t2.ratios[d])/t1.ratios[d]);
78  }
79  err /= 2;
80  return err / (align_info.dims + t1.dims);
81 }
82 
84 {
85  dsp_align_info align_info;
86  int dims = t1.dims - 1;
87  int d, x;
88  align_info.dims = dims;
89  align_info.center = (double*)malloc(sizeof(double)*(dims));
90  align_info.offset = (double*)malloc(sizeof(double)*(dims));
91  align_info.radians = (double*)malloc(sizeof(double)*(dims-1));
92  align_info.factor = (double*)malloc(sizeof(double)*(dims));
93  for(d = 0; d < dims; d++) {
94  align_info.factor[d] = 0;
95  align_info.offset[d] = 0;
96  align_info.center[d] = 0;
97  if(d < dims - 1)
98  align_info.radians[d] = 0;
99  }
100  for(d = 0; d < dims; d++) {
101  align_info.center[d] = t2.stars[0].center.location[d];
102  align_info.offset[d] = t2.stars[0].center.location[d] - t1.stars[0].center.location[d];
103  if(d < dims - 1) {
104  align_info.radians[d] = t1.theta[d] - t2.theta[d];
105  if(align_info.radians[d] >= M_PI*2.0)
106  align_info.radians[d] -= M_PI*2.0;
107  if(align_info.radians[d] < 0.0)
108  align_info.radians[d] += M_PI*2.0;
109  }
110  for(x = 0; x < t2.dims; x++) {
111  align_info.factor[d] += t2.sizes[x] / t1.sizes[x];
112  }
113  align_info.factor[d] /= t2.dims;
114  }
115  align_info.score = calc_match_score(t1, t2, align_info);
116  return align_info;
117 }
118 
120 {
121  int x;
122  int d;
123  dsp_triangle triangle;
124  triangle.dims = stars[0].center.dims+1;
125  double *delta = (double*)malloc(sizeof(double)*triangle.dims);
126  double **diff = (double**)malloc(sizeof(double*)*triangle.dims);
127  triangle.sizes = (double*)malloc(sizeof(double)*triangle.dims);
128  triangle.ratios = (double*)malloc(sizeof(double)*triangle.dims);
129  triangle.stars = (dsp_star*)malloc(sizeof(dsp_star)*triangle.dims);
130  triangle.theta = (double*)malloc(sizeof(double)*(triangle.dims-1));
131  triangle.index = 0.0;
132  qsort(stars, 3, sizeof(dsp_star), dsp_qsort_star_diameter_desc);
133  for(d = 0; d < triangle.dims; d++) {
134  diff[d] = (double*)malloc(sizeof(double)*stars[d].center.dims);
135  for(x = 0; x < triangle.dims-1; x++) {
136  diff[d][x] = stars[(d + 1) < triangle.dims ? d : 0].center.location[x]-stars[(d + 1) < triangle.dims ? (d + 1) : (triangle.dims - 1)].center.location[x];
137  delta[d] += pow(diff[d][x], 2);
138  }
139  delta[d] = sqrt(delta[d]);
140  }
141  for(d = 0; d < triangle.dims; d++) {
142  for(x = 0; x < stars[0].center.dims - 1; x++) {
143  if(d == 0) {
144  triangle.theta[x] = acos(diff[d][x] / delta[d]);
145  if(diff[0][x+1] < 0)
146  triangle.theta[x] = -triangle.theta[x];
147  if(triangle.theta[x] >= M_PI*2.0)
148  triangle.theta[x] -= M_PI*2.0;
149  if(triangle.theta[x] < 0)
150  triangle.theta[x] += M_PI*2.0;
151  } else {
152  double index = acos(diff[d][x] / delta[d]);
153  if(diff[d][x+1] < 0)
154  index = -index;
155  index -= triangle.theta[x];
156  triangle.index += index;
157  }
158  }
159  }
160  while(triangle.index >= M_PI)
161  triangle.index -= M_PI*2.0;
162  while(triangle.index < -M_PI)
163  triangle.index += M_PI*2.0;
164  qsort(delta, triangle.dims, sizeof(double), dsp_qsort_double_desc);
165  for(d = 0; d < triangle.dims; d++) {
166  triangle.stars[d].diameter = stars[d].diameter;
167  triangle.stars[d].center.dims = stars[d].center.dims;
168  triangle.stars[d].center.location = (double*)malloc(sizeof(double)*stars[d].center.dims);
169  triangle.sizes[d] = delta[d];
170  triangle.ratios[d] = delta[d] / delta[0];
171  for(x = 0; x < triangle.stars[d].center.dims; x++) {
172  triangle.stars[d].center.location[x] = stars[d].center.location[x];
173  }
174  free(diff[d]);
175  }
176  free(delta);
177  free(diff);
178  return triangle;
179 }
180 
181 int dsp_align_get_offset(dsp_stream_p stream1, dsp_stream_p stream2, double tolerance, double target_score)
182 {
183  dsp_align_info align_info;
184  double decimals = pow(10, tolerance);
185  double div = 0.0;
186  int d, t1, t2, x, y;
187  double phi = 0.0;
188  int dims = stream1->dims+1;
189  double min_score = 1.0;
190  double ratio = decimals*1600.0/div;
191  dsp_star *stars = (dsp_star*)malloc(sizeof(dsp_star)*dims);
192  for(d = 0; d < stream1->dims; d++) {
193  div += pow(stream2->sizes[d], 2);
194  }
195  div = pow(div, 0.5);
196  pwarn("decimals: %lf\n", decimals);
197  target_score = 1.0-target_score/100.0;
198  stream2->align_info.dims = 2;
199  stream2->align_info.score = 1.0;
200  stream2->align_info.decimals = decimals;
201  pgarb("creating triangles for reference frame...\n");
202  stream1->align_info.triangles_count = 0;
203  for(x = 0; x < stream1->stars_count; x++) {
204  for(y = x+1; y < stream1->stars_count-dims+1; y++) {
205  for(d = 0; d < dims; d++) {
206  stars[d] = stream1->stars[y+d];
207  }
209  dsp_stream_add_triangle(stream1, t);
210  }
211  }
212  pgarb("creating triangles for current frame...\n");
213  stream2->align_info.triangles_count = 0;
214  for(x = 0; x < stream2->stars_count; x++) {
215  for(y = x; y < stream2->stars_count-dims+1; y++) {
216  for(d = 0; d < dims; d++) {
217  stars[d] = stream2->stars[y+d];
218  }
220  dsp_stream_add_triangle(stream2, t);
221  }
222  }
223  free(stars);
224  for(t1 = 0; t1 < stream1->triangles_count; t1++) {
225  for(t2 = 0; t2 < stream2->triangles_count; t2++) {
226  align_info = dsp_align_fill_info(stream1->triangles[t1], stream2->triangles[t2]);
227  if(align_info.score < min_score) {
228  stream2->align_info = align_info;
229  min_score = align_info.score;
230  }
231  }
232  }
233  double radians = stream2->align_info.radians[0];
234  if(radians > M_PI)
235  radians -= M_PI*2;
236  for(d = 0; d < stream1->dims; d++) {
237  phi += pow(stream2->align_info.offset[d], 2);
238  }
239  phi = pow(phi, 0.5);
240  stream2->align_info.err = 0xf;
241  if(min_score < target_score)
242  stream2->align_info.err &= ~DSP_ALIGN_NO_MATCH;
243  if(fabs(phi * ratio * 10.0) < 1.0)
244  stream2->align_info.err &= ~DSP_ALIGN_TRANSLATED;
245  if(fabs((stream2->align_info.factor[0] - 1.0) * 100.0) * decimals < 1)
246  stream2->align_info.err &= ~DSP_ALIGN_SCALED;
247  if(fabs(radians) * decimals * 10.0 < 1)
248  stream2->align_info.err &= ~DSP_ALIGN_ROTATED;
249  return stream2->align_info.err;
250 }
#define pgarb(...)
Definition: dsp.h:155
#define DSP_ALIGN_ROTATED
The stream is rotated by the reference.
Definition: dsp.h:199
#define DSP_ALIGN_NO_MATCH
No matches were found during comparison.
Definition: dsp.h:203
#define DSP_ALIGN_SCALED
The stream is scaled by the reference.
Definition: dsp.h:195
#define DSP_ALIGN_TRANSLATED
The stream is translated by the reference.
Definition: dsp.h:191
#define pwarn(...)
Definition: dsp.h:154
int err
Errors.
Definition: dsp.h:290
dsp_point center
The center of the star.
Definition: dsp.h:240
int triangles_count
Triangles of stars or objects quantity.
Definition: dsp.h:423
int dims
Dimensions limit.
Definition: dsp.h:280
double * factor
Scaling factor.
Definition: dsp.h:278
double decimals
Decimals.
Definition: dsp.h:288
double * center
Center of rotation coordinates.
Definition: dsp.h:274
int * sizes
Sizes of each dimension.
Definition: dsp.h:373
double * location
Center of the point.
Definition: dsp.h:218
double * radians
Rotational offset.
Definition: dsp.h:276
double index
The index of the triangle.
Definition: dsp.h:253
dsp_star * stars
Stars or objects identified into the buffers - TODO.
Definition: dsp.h:417
dsp_star * stars
The stars of the triangle.
Definition: dsp.h:263
double * sizes
The sizes of the triangle.
Definition: dsp.h:259
int stars_count
Stars or objects quantity.
Definition: dsp.h:419
double diameter
The diameter of the star.
Definition: dsp.h:242
int dims
The dimensions of the triangle.
Definition: dsp.h:255
double * theta
The inclination of the triangle.
Definition: dsp.h:257
dsp_align_info align_info
Align/scale/rotation settings.
Definition: dsp.h:425
dsp_triangle * triangles
Triangles of stars or objects.
Definition: dsp.h:421
int dims
Number of dimensions of the buffers.
Definition: dsp.h:371
int triangles_count
Triangles quantity.
Definition: dsp.h:284
int dims
Dimensions limit of the point.
Definition: dsp.h:220
double score
Match score.
Definition: dsp.h:286
double * ratios
The sizes of the triangle.
Definition: dsp.h:261
double * offset
Translation offset.
Definition: dsp.h:272
DLL_EXPORT void dsp_stream_add_triangle(dsp_stream_p stream, dsp_triangle triangle)
Add a triangle to the DSP Stream passed as argument.
Definition: stream.c:367
int dsp_qsort_double_asc(const void *arg1, const void *arg2)
Callback function for qsort for double type ascending ordering.
Definition: align.c:31
dsp_align_info dsp_align_fill_info(dsp_triangle t1, dsp_triangle t2)
Fill a dsp_align_info struct by comparing two triangles.
Definition: align.c:83
int dsp_align_get_offset(dsp_stream_p stream1, dsp_stream_p stream2, double tolerance, double target_score)
Calculate offsets, rotation and scaling of two streams giving reference alignment point.
Definition: align.c:181
int dsp_qsort_star_diameter_asc(const void *arg1, const void *arg2)
Callback function for qsort for dsp_star ascending ordering by their diameters.
Definition: align.c:47
int dsp_qsort_star_diameter_desc(const void *arg1, const void *arg2)
Callback function for qsort for dsp_star descending ordering by their diameters.
Definition: align.c:38
int dsp_qsort_double_desc(const void *arg1, const void *arg2)
Callback function for qsort for double type descending ordering.
Definition: align.c:24
dsp_triangle dsp_align_calc_triangle(dsp_star *stars)
Create a dsp_triangle struct.
Definition: align.c:119
Alignment informations needed.
Definition: dsp.h:270
A star or object contained into a buffer.
Definition: dsp.h:238
Contains a set of informations and data relative to a buffer and how to use it.
Definition: dsp.h:363
A star or object contained into a buffer.
Definition: dsp.h:251