Changes backported from aptdec2

New regression, only linear as per recommendation of the Users' Guide

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)
 

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 <https://www.gnu.org/licenses/>.
+ */
+
+#include "algebra.h"
+#include <math.h>
+
+// 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;
+}

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 <https://www.gnu.org/licenses/>.
+ */
+
+#ifndef APTDEC_ALGEBRA_H
+#define APTDEC_ALGEBRA_H
+#include <stddef.h>
+
+// 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

src/image.c

@@ -23,32 +23,14 @@
 #include <stdlib.h>
 
 #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], &regr);
+			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, &regr[k]);
+			regr[k] = compute_regression(wedge);
 			for (int j = 0; j < 16; j++)
-				tele[j] = (float)rgcal((float)wedge[j], &regr[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), &regr[k]);
+			Cs = linear_calc((Cs / i), regr[k]);
 		}
 	}
 

src/libs/reg.c

@@ -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 <math.h>
-
-#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];
-}

src/libs/reg.h

@@ -1 +0,0 @@
-void polyreg(const int M, const int N, const double X[], const double Y[], double C[]);