Skip to content

Commit

Permalink
Add peak, flux and deviation parameters to dsp_star
Browse files Browse the repository at this point in the history
  • Loading branch information
Ilia Platone committed Nov 27, 2022
1 parent c80918c commit 01133df
Show file tree
Hide file tree
Showing 3 changed files with 136 additions and 40 deletions.
113 changes: 86 additions & 27 deletions dsp/align.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,11 +33,61 @@ static int dsp_qsort_delta_diff_desc (const void *arg1, const void *arg2)
return ((*a1).delta < (*a2).delta ? 1 : -1);
}

static int dsp_qsort_double_desc (const void *arg1, const void *arg2)
static int dsp_qsort_delta_diff_asc (const void *arg1, const void *arg2)
{
double* a = (double*)arg1;
double* b = (double*)arg2;
return ((*a) < (*b) ? 1 : -1);
delta_diff* a1 = (delta_diff*)arg1;
delta_diff* a2 = (delta_diff*)arg2;
return ((*a1).delta > (*a2).delta ? 1 : -1);
}

static int dsp_qsort_triangle_delta_desc(const void *arg1, const void *arg2)
{
dsp_triangle* a = (dsp_triangle*)arg1;
dsp_triangle* b = (dsp_triangle*)arg2;
int ret = 0;
int s = 0;
for(s = 0; s < (*a).stars_count; s++)
ret |= ((*a).sizes[s] < (*b).sizes[s]);
return (ret ? 1 : -1);
}

static int dsp_qsort_triangle_delta_asc(const void *arg1, const void *arg2)
{
dsp_triangle* a = (dsp_triangle*)arg1;
dsp_triangle* b = (dsp_triangle*)arg2;
int ret = 0;
int s = 0;
for(s = 0; s < (*a).stars_count; s++)
ret |= ((*a).sizes[s] > (*b).sizes[s]);
return (ret ? 1 : -1);
}

static int dsp_qsort_triangle_diameter_desc(const void *arg1, const void *arg2)
{
dsp_triangle* a = (dsp_triangle*)arg1;
dsp_triangle* b = (dsp_triangle*)arg2;
int ret = 0;
int s = 0;
for(s = 0; s < (*a).stars_count; s++) {
ret |= ((*a).stars[s].diameter < (*b).stars[s].diameter);
ret |= ((*a).stars[s].peak < (*b).stars[s].peak);
ret |= ((*a).stars[s].flux < (*b).stars[s].flux);
}
return (ret ? 1 : -1);
}

static int dsp_qsort_triangle_diameter_asc(const void *arg1, const void *arg2)
{
dsp_triangle* a = (dsp_triangle*)arg1;
dsp_triangle* b = (dsp_triangle*)arg2;
int ret = 0;
int s = 0;
for(s = 0; s < (*a).stars_count; s++) {
ret |= ((*a).stars[s].diameter > (*b).stars[s].diameter);
ret |= ((*a).stars[s].peak > (*b).stars[s].peak);
ret |= ((*a).stars[s].flux > (*b).stars[s].flux);
}
return (ret ? 1 : -1);
}

static int dsp_qsort_star_diameter_desc(const void *arg1, const void *arg2)
Expand All @@ -47,38 +97,50 @@ static int dsp_qsort_star_diameter_desc(const void *arg1, const void *arg2)
return ((*a).diameter < (*b).diameter ? 1 : -1);
}

static int dsp_qsort_double_asc (const void *arg1, const void *arg2)
static int dsp_qsort_star_diameter_asc(const void *arg1, const void *arg2)
{
dsp_star* a = (dsp_star*)arg1;
dsp_star* b = (dsp_star*)arg2;
return ((*a).diameter > (*b).diameter ? 1 : -1);
}

static int dsp_qsort_double_desc (const void *arg1, const void *arg2)
{
double* a = (double*)arg1;
double* b = (double*)arg2;
return ((*a) > (*b) ? 1 : -1);
return ((*a) < (*b) ? 1 : -1);
}

static int dsp_qsort_star_diameter_asc(const void *arg1, const void *arg2)
static int dsp_qsort_double_asc (const void *arg1, const void *arg2)
{
dsp_star* a = (dsp_star*)arg1;
dsp_star* b = (dsp_star*)arg2;
return ((*a).diameter > (*b).diameter ? 1 : -1);
double* a = (double*)arg1;
double* b = (double*)arg2;
return ((*a) > (*b) ? 1 : -1);
}

static double calc_match_score(dsp_triangle t1, dsp_triangle t2, dsp_align_info align_info)
{
int x = 0;
int d = 0;
int div = 0;
int num_baselines = t1.stars_count*(t1.stars_count-1)/2;
int stars_count = fmin(t1.stars_count, t2.stars_count);
int num_baselines = stars_count*(stars_count-1)/2;
double err = 0.0;
for(x = 0; x < fmin(t1.stars_count, t2.stars_count); x++) {
err += fabs(t2.stars[x].diameter/align_info.factor[d]-t1.stars[x].diameter)/t1.stars[x].diameter;
div++;
}
for(x = 0; x < num_baselines; x++) {
for(d = 0; d < align_info.dims; d++) {
err += fabs((t2.index-t1.index)/t1.index);
div++;
for(d = 0; d < align_info.dims; d++) {
for(x = 0; x < stars_count; x++) {
err += fabs(t2.stars[x].diameter/align_info.factor[d]-t1.stars[x].diameter)/t1.stars[x].diameter;
div++;
err += fabs(t2.stars[x].flux/align_info.factor[d]-t1.stars[x].flux)/t1.stars[x].flux;
div++;
err += fabs(t2.stars[x].peak/align_info.factor[d]-t1.stars[x].peak)/t1.stars[x].peak;
div++;
}
for(x = 0; x < num_baselines; x++) {
err += fabs(t2.sizes[x]/align_info.factor[d]-t1.sizes[x])/t1.sizes[x];
div++;
}
err += fabs(t2.ratios[x]-t1.ratios[x])/t1.ratios[x];
div++;
}
return err / div;
}
Expand Down Expand Up @@ -150,8 +212,6 @@ dsp_triangle *dsp_align_calc_triangle(dsp_star* stars, int num_stars)
triangle->index = 0.0;
qsort(stars, triangle->stars_count, sizeof(dsp_star), dsp_qsort_star_diameter_desc);
int idx = 0;
int maximum_index = 0;
double maximum = 0.0;
for(x = 0; x < triangle->stars_count; x++) {
for(y = x+1; y < triangle->stars_count; y++) {
deltadiff[idx].diff = (double*)malloc(sizeof(double)*triangle->dims);
Expand All @@ -160,10 +220,6 @@ dsp_triangle *dsp_align_calc_triangle(dsp_star* stars, int num_stars)
deltadiff[idx].delta += pow(deltadiff[idx].diff[d], 2);
}
deltadiff[idx].delta = sqrt(deltadiff[idx].delta);
if(maximum < deltadiff[idx].delta) {
maximum = deltadiff[idx].delta;
maximum_index = idx;
}
idx++;
}
}
Expand All @@ -183,6 +239,9 @@ dsp_triangle *dsp_align_calc_triangle(dsp_star* stars, int num_stars)
triangle->index /= (num_baselines+triangle->dims-2)*M_PI*2;
for(d = 0; d < triangle->stars_count; d++) {
triangle->stars[d].diameter = stars[d].diameter;
triangle->stars[d].peak = stars[d].peak;
triangle->stars[d].flux = stars[d].flux;
triangle->stars[d].theta = stars[d].theta;
triangle->stars[d].center.dims = stars[d].center.dims;
triangle->stars[d].center.location = (double*)malloc(sizeof(double)*stars[d].center.dims);
for(x = 0; x < triangle->stars[d].center.dims; x++) {
Expand All @@ -191,7 +250,7 @@ dsp_triangle *dsp_align_calc_triangle(dsp_star* stars, int num_stars)
}
for(d = 0; d < num_baselines; d++) {
triangle->sizes[d] = deltadiff[d].delta;
triangle->ratios[d] = deltadiff[d].delta / maximum;
triangle->ratios[d] = deltadiff[d].delta / deltadiff[0].delta;
free(deltadiff[d].diff);
}
free(deltadiff);
Expand All @@ -213,7 +272,7 @@ int dsp_align_get_offset(dsp_stream_p stream1, dsp_stream_p stream2, double tole
div = pow(div, 0.5);
pwarn("decimals: %lf\n", decimals);
target_score = (1.0-target_score / 100.0);
stream2->align_info.dims = 2;
stream2->align_info.dims = stream2->dims;
stream2->align_info.triangles_count = 2;
stream2->align_info.score = 1.0;
stream2->align_info.decimals = decimals;
Expand Down
6 changes: 6 additions & 0 deletions dsp/dsp.h.cmake
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,6 +240,12 @@ typedef struct dsp_star_t
dsp_point center;
/// The diameter of the star
double diameter;
/// The peak of the star
double peak;
/// The flux of the star
double flux;
/// The deviation of the star
double theta;
/// The name of the star
char name[150];
} dsp_star;
Expand Down
57 changes: 44 additions & 13 deletions dsp/stream.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,10 @@ void dsp_stream_free_buffer(dsp_stream_p stream)
free(stream->buf);
if(stream->dft.buf != NULL)
free(stream->dft.buf);
if(stream->magnitude != NULL)
dsp_stream_free_buffer(stream->magnitude);
if(stream->phase != NULL)
dsp_stream_free_buffer(stream->phase);
}

/**
Expand Down Expand Up @@ -180,6 +184,10 @@ void dsp_stream_free(dsp_stream_p stream)
free(stream->stars);
if(stream->triangles != NULL)
free(stream->triangles);
if(stream->magnitude != NULL)
dsp_stream_free(stream->magnitude);
if(stream->phase != NULL)
dsp_stream_free(stream->phase);
free(stream);
stream = NULL;
}
Expand Down Expand Up @@ -332,6 +340,9 @@ void dsp_stream_add_star(dsp_stream_p stream, dsp_star star)
stream->stars = (dsp_star*)realloc(stream->stars, sizeof(dsp_star)*(stream->stars_count+1));
strcpy(stream->stars[stream->stars_count].name, star.name);
stream->stars[stream->stars_count].diameter = star.diameter;
stream->stars[stream->stars_count].peak = star.peak;
stream->stars[stream->stars_count].flux = star.flux;
stream->stars[stream->stars_count].theta = star.theta;
stream->stars[stream->stars_count].center.dims = star.center.dims;
stream->stars[stream->stars_count].center.location = (double*)malloc(sizeof(double)*star.center.dims);
for(d = 0; d < star.center.dims; d++)
Expand Down Expand Up @@ -368,26 +379,35 @@ void dsp_stream_add_triangle(dsp_stream_p stream, dsp_triangle triangle)
{
int s;
int d;
int num_baselines = triangle.stars_count*(triangle.stars_count-1)/2;
stream->triangles = (dsp_triangle*)realloc(stream->triangles, sizeof(dsp_triangle)*(stream->triangles_count+1));
stream->triangles[stream->triangles_count].dims = triangle.dims;
stream->triangles[stream->triangles_count].index = triangle.index;
stream->triangles[stream->triangles_count].stars_count = triangle.stars_count;
stream->triangles[stream->triangles_count].theta = (double*)malloc(sizeof(double)*(stream->dims-1));
stream->triangles[stream->triangles_count].ratios = (double*)malloc(sizeof(double)*triangle.dims);
stream->triangles[stream->triangles_count].sizes = (double*)malloc(sizeof(double)*triangle.dims);
stream->triangles[stream->triangles_count].stars = (dsp_star*)malloc(sizeof(dsp_star)*triangle.dims);
stream->triangles[stream->triangles_count].ratios = (double*)malloc(sizeof(double)*num_baselines);
stream->triangles[stream->triangles_count].sizes = (double*)malloc(sizeof(double)*num_baselines);
stream->triangles[stream->triangles_count].stars = (dsp_star*)malloc(sizeof(dsp_star)*triangle.stars_count);
for (s = 0; s < triangle.dims; s++) {
if(s < stream->dims - 1) {
stream->triangles[stream->triangles_count].theta[s] = triangle.theta[s];
}
stream->triangles[stream->triangles_count].sizes[s] = triangle.sizes[s];
stream->triangles[stream->triangles_count].ratios[s] = triangle.ratios[s];
}
for(s = 0; s < triangle.stars_count; s++) {
stream->triangles[stream->triangles_count].stars[s].center.dims = triangle.stars[s].center.dims;
stream->triangles[stream->triangles_count].stars[s].diameter = triangle.stars[s].diameter;
stream->triangles[stream->triangles_count].stars[s].peak = triangle.stars[s].peak;
stream->triangles[stream->triangles_count].stars[s].flux = triangle.stars[s].flux;
stream->triangles[stream->triangles_count].stars[s].theta = triangle.stars[s].theta;
stream->triangles[stream->triangles_count].stars[s].center.location = (double*)malloc(sizeof(double)*stream->dims);
for(d = 0; d < triangle.stars[s].center.dims; d++) {
stream->triangles[stream->triangles_count].stars[s].center.location[d] = triangle.stars[s].center.location[d];
}
}
for(s = 0; s < num_baselines; s++) {
stream->triangles[stream->triangles_count].sizes[s] = triangle.sizes[s];
stream->triangles[stream->triangles_count].ratios[s] = triangle.ratios[s];
}
stream->triangles_count++;
}

Expand All @@ -398,16 +418,27 @@ void dsp_stream_add_triangle(dsp_stream_p stream, dsp_triangle triangle)
*/
void dsp_stream_del_triangle(dsp_stream_p stream, int index)
{
dsp_triangle* triangles = (dsp_triangle*)malloc(sizeof(dsp_triangle) * stream->triangles_count);
int triangles_count = stream->triangles_count;
memcpy(triangles, stream->triangles, sizeof(dsp_triangle*) * stream->triangles_count);
free(stream->triangles);
stream->triangles_count = 0;
int d;
free(stream->triangles[index].sizes);
free(stream->triangles[index].theta);
free(stream->triangles[index].ratios);
for(d = 0; d < stream->triangles[index].stars_count; d++) {
free(stream->triangles[index].stars[d].center.location);
}
free(stream->triangles[index].stars);
int i;
for(i = 0; i < triangles_count; i++) {
if(i != index) {
dsp_stream_add_triangle(stream, triangles[i]);
for(i = index; i < stream->triangles_count - 1; i++) {
stream->triangles[i] = stream->triangles[i+1];
}
stream->triangles_count--;
if(index < stream->triangles_count) {
free(stream->triangles[i].sizes);
free(stream->triangles[i].theta);
free(stream->triangles[i].ratios);
for(d = 0; d < stream->triangles[i].dims; d++) {
free(stream->triangles[i].stars[d].center.location);
}
free(stream->triangles[i].stars);
}
}

Expand Down

0 comments on commit 01133df

Please sign in to comment.