Update Libdsp (#1920)

Co-authored-by: Ilia Platone <[email protected]>

libs/dsp/convolution.c

@@ -31,7 +31,8 @@ void dsp_convolution_convolution(dsp_stream_p stream, dsp_stream_p matrix) {
         }
         x = dsp_stream_set_position(stream, d_pos);
         free(pos);
-        stream->magnitude->buf[x] *= sqrt(matrix->magnitude->buf[y]);
+        if(x >= 0 && x < stream->magnitude->len)
+            stream->magnitude->buf[x] *= sqrt(matrix->magnitude->buf[y]);
     }
     free(d_pos);
     dsp_fourier_idft(stream);

libs/dsp/dsp.h

@@ -24,10 +24,14 @@ extern "C" {
 #endif
 
 #ifndef DLL_EXPORT
+#ifdef _WIN32
+#define DLL_EXPORT __declspec(dllexport)
+#else
 #define DLL_EXPORT extern
 #endif
+#endif
 
-#if defined __linux__ || defined __CYGWIN__
+#ifdef __linux__
 #include <endian.h>
 #else
 #define __bswap_16(a) __builtin_bswap16(a)
@@ -70,6 +74,7 @@ typedef double dsp_t;
 typedef double complex_t[2];
 #define dsp_t_max 255
 #define dsp_t_min -dsp_t_max
+#define DSP_NAME_SIZE 128
 
 /**
 * \brief get/set the maximum number of threads allowed
@@ -240,8 +245,14 @@ 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];
+    char name[DSP_NAME_SIZE];
 } dsp_star;
 
 /**
@@ -261,6 +272,8 @@ typedef struct dsp_triangle_t
     double *ratios;
     /// The stars of the triangle
     dsp_star *stars;
+    /// The number of stars of the triangle
+    int stars_count;
 } dsp_triangle;
 
 /**
@@ -362,7 +375,7 @@ typedef void *(*dsp_func_t) (void *, ...);
 typedef struct dsp_stream_t
 {
     /// Friendly name of the stream
-    char name[128];
+    char name[DSP_NAME_SIZE];
     /// Increments by one on the copied stream
     int is_copy;
     /// The buffers length
@@ -1259,6 +1272,27 @@ DLL_EXPORT void dsp_stream_crop(dsp_stream_p stream);
 */
 DLL_EXPORT void dsp_stream_rotate(dsp_stream_p stream);
 
+/**
+* \brief Sum a stream with another
+* \param stream The stream that will be summed
+* \param in The stream that will be summed to stream
+*/
+DLL_EXPORT void dsp_stream_sum(dsp_stream_p stream, dsp_stream_p in);
+
+/**
+* \brief Multiply a stream with another
+* \param stream The stream that will be multiplied
+* \param in The stream that will be multiplied to stream
+*/
+DLL_EXPORT void dsp_stream_multiply(dsp_stream_p stream, dsp_stream_p in);
+
+/**
+* \brief Subtract a stream to another
+* \param stream The original stream
+* \param in The stream that will be subtracted to the first argument
+*/
+DLL_EXPORT void dsp_stream_subtract(dsp_stream_p stream, dsp_stream_p in);
+
 /**
 * \brief Translate a stream
 * \param stream The stream that will be translated in-place
@@ -1275,7 +1309,7 @@ DLL_EXPORT void dsp_stream_scale(dsp_stream_p stream);
 * \brief Perform scale, translate and rotate transformations in-place
 * \param stream The stream that will be transformed
 */
-DLL_EXPORT void dsp_stream_align(dsp_stream_p in);
+DLL_EXPORT void dsp_stream_align(dsp_stream_p stream);
 
 /**\}*/
 /**
@@ -1330,6 +1364,23 @@ DLL_EXPORT void dsp_modulation_frequency(dsp_stream_p stream, double samplefreq,
 */
 DLL_EXPORT void dsp_modulation_amplitude(dsp_stream_p stream, double samplefreq, double freq);
 
+/**\}*/
+/**
+ * \defgroup dsp_StreamReconstruction DSP API Stream reconstruction
+*/
+/**\{*/
+
+/**
+* \brief Try to reconstruct a stream from sub-stream matrices
+* \param stream the target DSP stream.
+* \param matrix the DSP stream with the matrices.
+*
+* The matrix stream must have all dimensions and sizes of the target plus one,
+* which is the number of matrices on to align.
+* The aligned stream should replace the target stream.
+*/
+void dsp_recons_align(dsp_stream_p stream, dsp_stream_p matrix);
+
 /**\}*/
 /**
  * \defgroup dsp_FileManagement DSP API File read/write functions
@@ -1389,6 +1440,25 @@ DLL_EXPORT void dsp_file_write_jpeg(const char* filename, int quality, dsp_strea
 */
 DLL_EXPORT void dsp_file_write_jpeg_composite(const char* filename, int components, int quality, dsp_stream_p* stream);
 
+/**
+* \brief Read a PNG file and fill a array of dsp_stream_p with its content,
+* each color channel has its own stream in this array and an additional grayscale at end will be added
+* \param filename the file name.
+* \param channels this value will be updated with the channel quantity into the picture.
+* \param stretch 1 if the buffer intensities have to be stretched
+* \return The new dsp_stream_p structure pointers array
+*/
+DLL_EXPORT dsp_stream_p* dsp_file_read_png(const char* filename, int *channels, int stretch);
+
+/**
+* \brief Write the components dsp_stream_p array into a PNG file,
+* \param filename the file name.
+* \param components the number of streams in the array to be used as components 1 or 3.
+* \param compression the compression of the output PNG 0-9.
+* \param stream the input stream array to be saved
+*/
+DLL_EXPORT void dsp_file_write_png_composite(const char* filename, int components, int compression, dsp_stream_p* stream);
+
 /**
 * \brief Convert a bayer pattern dsp_t array into a grayscale array
 * \param src the input buffer
@@ -1472,58 +1542,34 @@ DLL_EXPORT dsp_t* dsp_file_bayer_2_composite(dsp_t *src, int red, int width, int
 * \brief Fill a dsp_align_info struct by comparing two triangles
 * \param t1 the reference triangle
 * \param t2 the triangle taken for comparison
-* \return The dsp_align_info struct filled with the offsets
+* \return The dsp_align_info struct pointer filled with the offsets
 */
-DLL_EXPORT dsp_align_info dsp_align_fill_info(dsp_triangle t1, dsp_triangle t2);
+DLL_EXPORT dsp_align_info *dsp_align_fill_info(dsp_triangle t1, dsp_triangle t2);
 
 /**
 * \brief Create a dsp_triangle struct
-* \param stars the stars array meeded to build the triangle struct
-* \return A new dsp_triangle struct
+* \param stars the stars array needed to build the triangle struct
+* \param num_stars the count of the stars
+* \return A new dsp_triangle struct pointer
+*/
+DLL_EXPORT dsp_triangle *dsp_align_calc_triangle(dsp_star* stars, int num_stars);
+
+/**
+* \brief Free a dsp_triangle struct pointer
+* \param triangle pointer to an allocated dsp_triangle struct
 */
-DLL_EXPORT dsp_triangle dsp_align_calc_triangle(dsp_star* stars);
+DLL_EXPORT void dsp_align_free_triangle(dsp_triangle *triangle);
 
 /**
 * \brief Calculate offsets, rotation and scaling of two streams giving reference alignment point
 * \param ref the reference stream
 * \param to_align the stream to be aligned
 * \param tolerance number of decimals allowed
 * \param target_score the minimum matching score to reach
+* \param num_stars number of stars for each triangle
 * \return The alignment mask (bit1: translated, bit2: scaled, bit3: rotated)
 */
-DLL_EXPORT int dsp_align_get_offset(dsp_stream_p ref, dsp_stream_p to_align, double tolerance, double target_score);
-
-/**
-* \brief Callback function for qsort for double type ascending ordering
-* \param arg1 the first comparison element
-* \param arg2 the second comparison element
-* \return 1 if arg1 is greater than arg2, -1 if arg2 is greater than arg1
-*/
-DLL_EXPORT int dsp_qsort_double_asc (const void *arg1, const void *arg2);
-
-/**
-* \brief Callback function for qsort for double type descending ordering
-* \param arg1 the first comparison element
-* \param arg2 the second comparison element
-* \return 1 if arg2 is greater than arg1, -1 if arg1 is greater than arg2
-*/
-DLL_EXPORT int dsp_qsort_double_desc (const void *arg1, const void *arg2);
-
-/**
-* \brief Callback function for qsort for dsp_star ascending ordering by their diameters
-* \param arg1 the first comparison element
-* \param arg2 the second comparison element
-* \return 1 if arg1 diameter is greater than arg2, -1 if arg2 diameter is greater than arg1
-*/
-DLL_EXPORT int dsp_qsort_star_diameter_asc (const void *arg1, const void *arg2);
-
-/**
-* \brief Callback function for qsort for dsp_star descending ordering by their diameters
-* \param arg1 the first comparison element
-* \param arg2 the second comparison element
-* \return 1 if arg2 diameter is greater than arg1, -1 if arg1 diameter is greater than arg2
-*/
-DLL_EXPORT int dsp_qsort_star_diameter_desc (const void *arg1, const void *arg2);
+DLL_EXPORT int dsp_align_get_offset(dsp_stream_p ref, dsp_stream_p to_align, double tolerance, double target_score, int num_stars);
 
 /**\}*/
 /// \defgroup dsp_FitsExtensions

libs/dsp/feature.c

@@ -0,0 +1,45 @@
+/*
+*   DSP API - a digital signal processing library for astronomy usage
+*   Copyright © 2017-2021  Ilia Platone
+*
+*   This program is free software; you can redistribute it and/or
+*   modify it under the terms of the GNU Lesser General Public
+*   License as published by the Free Software Foundation; either
+*   version 3 of the License, or (at your option) any later version.
+*
+*   This program is distributed in the hope that it will be useful,
+*   but WITHOUT ANY WARRANTY; without even the implied warranty of
+*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+*   Lesser General Public License for more details.
+*
+*   You should have received a copy of the GNU Lesser General Public License
+*   along with this program; if not, write to the Free Software Foundation,
+*   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+*/
+
+#include "dsp.h"
+
+void dsp_feature_align(dsp_stream_p stream, dsp_stream_p matrix)
+{
+    int z;
+    dsp_fourier_dft(matrix, 1);
+    for(z = 0; z < matrix->sizes[stream->dims]; z++) {
+        dsp_fourier_dft(stream, 1);
+        int x, y, d;
+        dsp_t mn = dsp_stats_min(stream->buf, stream->len);
+        dsp_t mx = dsp_stats_max(stream->buf, stream->len);
+        int* d_pos = (int*)malloc(sizeof(int)*stream->dims);
+        for(y = z*stream->len; y < z*stream->len+stream->len; y++) {
+            int* pos = dsp_stream_get_position(matrix, y);
+            for(d = 0; d < stream->dims; d++) {
+                d_pos[d] = stream->sizes[d]/2+pos[d]-matrix->sizes[d]/2;
+            }
+            x = dsp_stream_set_position(stream, d_pos);
+            free(pos);
+            stream->magnitude->buf[x] *= sqrt(matrix->magnitude->buf[y]);
+        }
+        free(d_pos);
+        dsp_fourier_idft(stream);
+        dsp_buffer_stretch(stream->buf, stream->len, mn, mx);
+    }
+}

libs/dsp/fft.c

@@ -144,8 +144,9 @@ void dsp_fourier_dft(dsp_stream_p stream, int exp)
     dsp_buffer_set(stream->dft.buf, stream->len * 2, 0);
     dsp_buffer_copy(stream->buf, buf, stream->len);
     int *sizes = (int*)malloc(sizeof(int)*stream->dims);
+    dsp_buffer_copy(stream->sizes, sizes, stream->dims);
     dsp_buffer_reverse(sizes, stream->dims);
-    fftw_plan plan = fftw_plan_dft_r2c(stream->dims, stream->sizes, buf, stream->dft.pairs, FFTW_ESTIMATE_PATIENT);
+    fftw_plan plan = fftw_plan_dft_r2c(stream->dims, sizes, buf, stream->dft.pairs, FFTW_ESTIMATE_PATIENT);
     fftw_execute(plan);
     fftw_free(plan);
     free(sizes);
@@ -177,8 +178,9 @@ void dsp_fourier_idft(dsp_stream_p stream)
     dsp_buffer_set(buf, stream->len, 0);
     dsp_fourier_2complex_t(stream);
     int *sizes = (int*)malloc(sizeof(int)*stream->dims);
+    dsp_buffer_copy(stream->sizes, sizes, stream->dims);
     dsp_buffer_reverse(sizes, stream->dims);
-    fftw_plan plan = fftw_plan_dft_c2r(stream->dims, stream->sizes, stream->dft.pairs, buf, FFTW_ESTIMATE_PATIENT);
+    fftw_plan plan = fftw_plan_dft_c2r(stream->dims, sizes, stream->dft.pairs, buf, FFTW_ESTIMATE_PATIENT);
     fftw_execute(plan);
     fftw_free(plan);
     free(sizes);