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[]);