From 27d75c706e3c8f9bb99c1d92aa99f658fd9a1751 Mon Sep 17 00:00:00 2001 From: Xerbo Date: Mon, 18 Jul 2022 21:15:18 +0100 Subject: [PATCH] Changes backported from aptdec2 New regression, only linear as per recommendation of the Users' Guide --- CMakeLists.txt | 2 +- src/algebra.c | 59 ++++++++++++++++++++ src/algebra.h | 38 +++++++++++++ src/image.c | 42 +++++---------- src/libs/reg.c | 143 ------------------------------------------------- src/libs/reg.h | 1 - 6 files changed, 110 insertions(+), 175 deletions(-) create mode 100644 src/algebra.c create mode 100644 src/algebra.h delete mode 100644 src/libs/reg.c delete mode 100644 src/libs/reg.h diff --git a/CMakeLists.txt b/CMakeLists.txt index efc4443..22f3790 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,7 +11,7 @@ find_package(PNG) # libsndfile find_package(LibSndFile) -set(LIB_C_SOURCE_FILES src/color.c src/dsp.c src/filter.c src/image.c src/libs/reg.c src/libs/median.c) +set(LIB_C_SOURCE_FILES src/color.c src/dsp.c src/filter.c src/image.c src/algebra.c src/libs/median.c) set(EXE_C_SOURCE_FILES src/main.c src/pngio.c src/libs/argparse.c) set(LIB_C_HEADER_FILES src/apt.h) diff --git a/src/algebra.c b/src/algebra.c new file mode 100644 index 0000000..e9b82b6 --- /dev/null +++ b/src/algebra.c @@ -0,0 +1,59 @@ +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2019-2022 Xerbo (xerbo@protonmail.com) + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "algebra.h" +#include + +// Find the best linear equation to estimate the value of the +// dependent variable from the independent variable +linear_t linear_regression(const float *independent, const float *dependent, size_t len) { + // Calculate mean of the dependent and independent variables (this is the centoid) + float dependent_mean = 0.0f; + float independent_mean = 0.0f; + for (size_t i = 0; i < len; i++) { + dependent_mean += dependent[i] / (float)len; + independent_mean += independent[i] / (float)len; + } + + // Calculate slope + float a = 0.0f; + { + float a_numerator = 0.0f; + float a_denominator = 0.0f; + for (size_t i = 0; i < len; i++) { + a_numerator += (independent[i] - independent_mean) * (dependent[i] - dependent_mean); + a_denominator += powf(independent[i] - independent_mean, 2.0f); + } + a = a_numerator/a_denominator; + } + + // We can now solve for the y-intercept since we know the slope + // and the centoid, which the line must pass through + float b = dependent_mean - a * independent_mean; + + //printf("y(x) = %fx + %f\n", a, b); + return (linear_t){a, b}; +} + +float linear_calc(float x, linear_t line) { + return x*line.a + line.b; +} + +float quadratic_calc(float x, quadratic_t quadratic) { + return x*x*quadratic.a + x*quadratic.b + quadratic.c; +} diff --git a/src/algebra.h b/src/algebra.h new file mode 100644 index 0000000..adbf932 --- /dev/null +++ b/src/algebra.h @@ -0,0 +1,38 @@ +/* + * aptdec - A lightweight FOSS (NOAA) APT decoder + * Copyright (C) 2019-2022 Xerbo (xerbo@protonmail.com) + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef APTDEC_ALGEBRA_H +#define APTDEC_ALGEBRA_H +#include + +// A linear equation in the form of y(x) = ax + b +typedef struct { + float a, b; +} linear_t; + +// A quadratic equation in the form of y(x) = ax^2 + bx + c +typedef struct { + float a, b, c; +} quadratic_t; + +linear_t linear_regression(const float *independent, const float *dependent, size_t len); + +float linear_calc(float x, linear_t line); +float quadratic_calc(float x, quadratic_t line); + +#endif diff --git a/src/image.c b/src/image.c index 224e2c4..875265f 100644 --- a/src/image.c +++ b/src/image.c @@ -23,32 +23,14 @@ #include #include "apt.h" -#include "libs/reg.h" +#include "algebra.h" #include "image.h" -#define REGORDER 3 -typedef struct { - double cf[REGORDER + 1]; -} rgparam_t; - -// Compute regression -static void rgcomp(double x[16], rgparam_t * rgpr) { +static linear_t compute_regression(float *wedges) { // { 0.106, 0.215, 0.324, 0.433, 0.542, 0.652, 0.78, 0.87, 0.0 } - const double y[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; - - polyreg(REGORDER, 9, x, y, rgpr->cf); -} - -// Convert a value to 0-255 based off the provided regression curve -static double rgcal(float x, rgparam_t *rgpr) { - double y, p; - int i; + const float teleramp[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; - for (i = 0, y = 0.0, p = 1.0; i < REGORDER + 1; i++) { - y += rgpr->cf[i] * p; - p = p * x; - } - return y; + return linear_regression(wedges, teleramp, 9); } static double tele[16]; @@ -104,18 +86,18 @@ void apt_linearEnhance(float **prow, int nrow, int offset, int width){ } // Brightness calibrate, including telemetry -void calibrateImage(float **prow, int nrow, int offset, int width, rgparam_t regr){ +void calibrateImage(float **prow, int nrow, int offset, int width, linear_t regr){ offset -= APT_SYNC_WIDTH+APT_SPC_WIDTH; for (int y = 0; y < nrow; y++) { for (int x = 0; x < width+APT_SYNC_WIDTH+APT_SPC_WIDTH+APT_TELE_WIDTH; x++) { - float pv = (float)rgcal(prow[y][x + offset], ®r); + float pv = linear_calc(prow[y][x + offset], regr); prow[y][x + offset] = CLIP(pv, 0, 255); } } } -double teleNoise(double wedges[16]){ +double teleNoise(float *wedges){ double pattern[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; double noise = 0; for(int i = 0; i < 9; i++) @@ -127,8 +109,8 @@ double teleNoise(double wedges[16]){ // Get telemetry data for thermal calibration apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { double teleline[APT_MAX_HEIGHT] = { 0.0 }; - double wedge[16]; - rgparam_t regr[APT_MAX_HEIGHT/APT_FRAME_LEN + 1]; + float wedge[16]; + linear_t regr[APT_MAX_HEIGHT/APT_FRAME_LEN + 1]; int telestart, mtelestart = 0; int channel = -1; @@ -194,9 +176,9 @@ apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { bestFrame = k; // Compute & apply regression on the wedges - rgcomp(wedge, ®r[k]); + regr[k] = compute_regression(wedge); for (int j = 0; j < 16; j++) - tele[j] = (float)rgcal((float)wedge[j], ®r[k]); + tele[j] = linear_calc(wedge[j], regr[k]); /* Compare the channel ID wedge to the reference * wedges, the wedge with the closest match will @@ -227,7 +209,7 @@ apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { i++; } } - Cs = rgcal((float)(Cs / i), ®r[k]); + Cs = linear_calc((Cs / i), regr[k]); } } diff --git a/src/libs/reg.c b/src/libs/reg.c deleted file mode 100644 index 975993f..0000000 --- a/src/libs/reg.c +++ /dev/null @@ -1,143 +0,0 @@ -/* --------------------------------------------------------------------------- - - Polynomial regression, freely adapted from: - - NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995 - Algorithm translated to C by: Dr. Norman Fahrer - NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992 - Prentice Hall, International Editions: ISBN 0-13-625047-5 - This free software is compliments of the author. - E-mail address: in%"mathews@fullerton.edu" -*/ - -#include - -#include "reg.h" - -#define DMAX 5 /* Maximum degree of polynomial */ -#define NMAX 10 /* Maximum number of points */ -static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]); - -void polyreg(const int M, const int N, const double X[], const double Y[], double C[]) { - int R, K, J; /* Loop counters */ - double A[DMAX][DMAX]; /* A */ - double B[DMAX]; - double P[2 * DMAX + 1]; - double x, y; - double p; - - /* Zero the array */ - for (R = 0; R < M + 1; R++) - B[R] = 0; - - /* Compute the column vector */ - for (K = 0; K < N; K++) { - y = Y[K]; - x = X[K]; - p = 1.0; - for (R = 0; R < M + 1; R++) { - B[R] += y * p; - p = p * x; - } - } - - /* Zero the array */ - for (J = 1; J <= 2 * M; J++) - P[J] = 0; - - P[0] = N; - - /* Compute the sum of powers of x_(K-1) */ - for (K = 0; K < N; K++) { - x = X[K]; - p = X[K]; - for (J = 1; J <= 2 * M; J++) { - P[J] += p; - p = p * x; - } - } - - /* Determine the matrix entries */ - for (R = 0; R < M + 1; R++) { - for (K = 0; K < M + 1; K++) - A[R][K] = P[R + K]; - } - - /* Solve the linear system of M + 1 equations: A*C = B - for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */ - FactPiv(M + 1, A, B, C); -} /* end main */ - - -/*--------------------------------------------------------*/ -static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) { - int K, P, C, J; /* Loop counters */ - int Row[NMAX]; /* Field with row-number */ - double X[DMAX], Y[DMAX]; - double SUM, DET = 1.0; - int T; - - /* Initialize the pointer vector */ - for (J = 0; J < NMAX; J++) - Row[J] = J; - - /* Start LU factorization */ - for (P = 0; P < N - 1; P++) { - /* Find pivot element */ - for (K = P + 1; K < N; K++) { - if (fabs(A[Row[K]][P]) > fabs(A[Row[P]][P])) { - /* Switch the index for the p-1 th pivot row if necessary */ - T = Row[P]; - Row[P] = Row[K]; - Row[K] = T; - DET = -DET; - } - } /* End of simulated row interchange */ - if (A[Row[P]][P] == 0) { - /* The matrix is SINGULAR! */ - return; - } - - /* Multiply the diagonal elements */ - DET = DET * A[Row[P]][P]; - - /* Form multiplier */ - for (K = P + 1; K < N; K++) { - A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P]; - - /* Eliminate X_(p-1) */ - for (C = P + 1; C < N + 1; C++) { - A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C]; - } - } - } - /* End of L*U factorization routine */ - DET = DET * A[Row[N - 1]][N - 1]; - - /* Start the forward substitution */ - for (K = 0; K < N; K++) - Y[K] = B[K]; - - Y[0] = B[Row[0]]; - - for (K = 1; K < N; K++) { - SUM = 0; - for (C = 0; C <= K - 1; C++) - SUM += A[Row[K]][C] * Y[C]; - Y[K] = B[Row[K]] - SUM; - } - - /* Start the back substitution */ - X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1]; - for (K = N - 2; K >= 0; K--) { - SUM = 0; - for (C = K + 1; C < N; C++) { - SUM += A[Row[K]][C] * X[C]; - } - X[K] = (Y[K] - SUM) / A[Row[K]][K]; - } /* End of back substitution */ - - /* Output */ - for (K = 0; K < N; K++) - Cf[K] = X[K]; -} diff --git a/src/libs/reg.h b/src/libs/reg.h deleted file mode 100644 index f18371a..0000000 --- a/src/libs/reg.h +++ /dev/null @@ -1 +0,0 @@ -void polyreg(const int M, const int N, const double X[], const double Y[], double C[]);