From 4d4a0c9787a27d1eba26b9299c23ae9c66e56716 Mon Sep 17 00:00:00 2001 From: Xerbo Date: Sun, 6 Nov 2022 14:42:03 +0000 Subject: [PATCH] Format codebase --- .clang-format | 4 + src/algebra.c | 13 +- src/apt.h | 50 ++-- src/calibration.c | 169 +++++------ src/calibration.h | 6 +- src/color.c | 105 +++---- src/color.h | 2 +- src/common.h | 48 +-- src/dsp.c | 358 +++++++++++----------- src/filter.c | 48 +-- src/filter.h | 6 +- src/image.c | 612 +++++++++++++++++++------------------- src/main.c | 521 ++++++++++++++++---------------- src/pngio.c | 736 +++++++++++++++++++++++----------------------- src/taps.h | 132 ++++----- src/util.c | 4 +- src/util.h | 4 +- 17 files changed, 1376 insertions(+), 1442 deletions(-) create mode 100644 .clang-format mode change 100755 => 100644 src/dsp.c mode change 100755 => 100644 src/filter.c mode change 100755 => 100644 src/filter.h diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..c62fa12 --- /dev/null +++ b/.clang-format @@ -0,0 +1,4 @@ +Language: Cpp +BasedOnStyle: Google +IndentWidth: 4 +ColumnLimit: 130 # Conservative estimate considering the size of modern displays diff --git a/src/algebra.c b/src/algebra.c index e9b82b6..bc7b9f6 100644 --- a/src/algebra.c +++ b/src/algebra.c @@ -17,6 +17,7 @@ */ #include "algebra.h" + #include // Find the best linear equation to estimate the value of the @@ -39,21 +40,17 @@ linear_t linear_regression(const float *independent, const float *dependent, siz a_numerator += (independent[i] - independent_mean) * (dependent[i] - dependent_mean); a_denominator += powf(independent[i] - independent_mean, 2.0f); } - a = a_numerator/a_denominator; + 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); + // 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 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; -} +float quadratic_calc(float x, quadratic_t quadratic) { return x * x * quadratic.a + x * quadratic.b + quadratic.c; } diff --git a/src/apt.h b/src/apt.h index 2478e9e..2e4d9ba 100644 --- a/src/apt.h +++ b/src/apt.h @@ -25,9 +25,9 @@ extern "C" { #endif -#if defined (__GNUC__) && (__GNUC__ >= 4) +#if defined(__GNUC__) && (__GNUC__ >= 4) #define APT_API __attribute__((visibility("default"))) -#elif defined (_MSC_VER) +#elif defined(_MSC_VER) #ifdef APT_API_EXPORT #define APT_API __declspec(dllexport) #elif APT_API_STATIC @@ -44,20 +44,20 @@ extern "C" { // Width in pixels of sync #define APT_SYNC_WIDTH 39 // Width in pixels of space -#define APT_SPC_WIDTH 47 +#define APT_SPC_WIDTH 47 // Width in pixels of telemetry #define APT_TELE_WIDTH 45 // Width in pixels of a single channel image -#define APT_CH_WIDTH 909 -#define APT_FRAME_LEN 128 -#define APT_CH_OFFSET (APT_SYNC_WIDTH+APT_SPC_WIDTH+APT_CH_WIDTH+APT_TELE_WIDTH) +#define APT_CH_WIDTH 909 +#define APT_FRAME_LEN 128 +#define APT_CH_OFFSET (APT_SYNC_WIDTH + APT_SPC_WIDTH + APT_CH_WIDTH + APT_TELE_WIDTH) // Width in pixels of full frame, including sync, space, images and telemetry -#define APT_IMG_WIDTH 2080 +#define APT_IMG_WIDTH 2080 // Offset in pixels to channel A -#define APT_CHA_OFFSET (APT_SYNC_WIDTH+APT_SPC_WIDTH) +#define APT_CHA_OFFSET (APT_SYNC_WIDTH + APT_SPC_WIDTH) // Offset in pixels to channel B -#define APT_CHB_OFFSET (APT_SYNC_WIDTH+APT_SPC_WIDTH+APT_CH_WIDTH+APT_TELE_WIDTH+APT_SYNC_WIDTH+APT_SPC_WIDTH) -#define APT_TOTAL_TELE (APT_SYNC_WIDTH+APT_SPC_WIDTH+APT_TELE_WIDTH+APT_SYNC_WIDTH+APT_SPC_WIDTH+APT_TELE_WIDTH) +#define APT_CHB_OFFSET (APT_SYNC_WIDTH + APT_SPC_WIDTH + APT_CH_WIDTH + APT_TELE_WIDTH + APT_SYNC_WIDTH + APT_SPC_WIDTH) +#define APT_TOTAL_TELE (APT_SYNC_WIDTH + APT_SPC_WIDTH + APT_TELE_WIDTH + APT_SYNC_WIDTH + APT_SPC_WIDTH + APT_TELE_WIDTH) // Number of rows required for apt_calibrate #define APT_CALIBRATION_ROWS 192 @@ -69,7 +69,15 @@ extern "C" { // Channel 3B: mid-infrared (3.55-3.93 um) // Channel 4: thermal-infrared (10.3-11.3 um) // Channel 5: thermal-infrared (11.5-12.5 um) -typedef enum apt_channel {APT_CHANNEL_UNKNOWN, APT_CHANNEL_1, APT_CHANNEL_2, APT_CHANNEL_3A, APT_CHANNEL_4, APT_CHANNEL_5, APT_CHANNEL_3B} apt_channel_t; +typedef enum apt_channel { + APT_CHANNEL_UNKNOWN, + APT_CHANNEL_1, + APT_CHANNEL_2, + APT_CHANNEL_3A, + APT_CHANNEL_4, + APT_CHANNEL_5, + APT_CHANNEL_3B +} apt_channel_t; // Width in elements of apt_image_t.prow arrays #define APT_PROW_WIDTH 2150 @@ -79,16 +87,16 @@ typedef enum apt_channel {APT_CHANNEL_UNKNOWN, APT_CHANNEL_1, APT_CHANNEL_2, APT typedef int (*apt_getsamples_t)(void *context, float *samples, int count); typedef struct { - float *prow[APT_MAX_HEIGHT]; // Row buffers - int nrow; // Number of rows - int zenith; // Row in image where satellite reaches peak elevation - apt_channel_t chA, chB; // ID of each channel - char name[256]; // Stripped filename - char *palette; // Filename of palette + float *prow[APT_MAX_HEIGHT]; // Row buffers + int nrow; // Number of rows + int zenith; // Row in image where satellite reaches peak elevation + apt_channel_t chA, chB; // ID of each channel + char name[256]; // Stripped filename + char *palette; // Filename of palette } apt_image_t; typedef struct { - float r, g, b; + float r, g, b; } apt_rgb_t; int APT_API apt_init(double sample_rate); @@ -96,7 +104,7 @@ int APT_API apt_getpixelrow(float *pixelv, int nrow, int *zenith, int reset, apt void APT_API apt_histogramEqualise(float **prow, int nrow, int offset, int width); void APT_API apt_linearEnhance(float **prow, int nrow, int offset, int width); -apt_channel_t APT_API apt_calibrate(float **prow, int nrow, int offset, int width) ; +apt_channel_t APT_API apt_calibrate(float **prow, int nrow, int offset, int width); void APT_API apt_denoise(float **prow, int nrow, int offset, int width); void APT_API apt_flipImage(apt_image_t *img, int width, int offset); int APT_API apt_cropNoise(apt_image_t *img); @@ -108,8 +116,8 @@ void APT_API apt_calibrate_visible(int satnum, apt_image_t *img, int offset, int apt_rgb_t APT_API apt_applyPalette(char *palette, int val); apt_rgb_t APT_API apt_RGBcomposite(apt_rgb_t top, float top_a, apt_rgb_t bottom, float bottom_a); -extern char APT_API apt_TempPalette[256*3]; -extern char APT_API apt_PrecipPalette[58*3]; +extern char APT_API apt_TempPalette[256 * 3]; +extern char APT_API apt_PrecipPalette[58 * 3]; #ifdef __cplusplus } diff --git a/src/calibration.c b/src/calibration.c index b2c7a27..e5b32fa 100644 --- a/src/calibration.c +++ b/src/calibration.c @@ -17,104 +17,85 @@ */ #include "calibration.h" + #include "util.h" -const calibration_t calibration[3] = { - { - .name = "NOAA-15", - .prt = { - { 1.36328e-06f, 0.051045f, 276.60157f }, // PRT 1 - { 1.47266e-06f, 0.050909f, 276.62531f }, // PRT 2 - { 1.47656e-06f, 0.050907f, 276.67413f }, // PRT 3 - { 1.47656e-06f, 0.050966f, 276.59258f } // PRT 4 - }, - .visible = { - { - .low = { 0.0568f, -2.1874f }, - .high = { 0.1633f, -54.9928f }, - .cutoff = 496.0f - }, { - .low = { 0.0596f, -2.4096f }, - .high = { 0.1629f, -55.2436f }, - .cutoff = 511.0f - } - }, - .rad = { - { 925.4075f, 0.337810f, 0.998719f }, // Channel 4 - { 839.8979f, 0.304558f, 0.999024f }, // Channel 5 - { 2695.9743f, 1.621256f, 0.998015f } // Channel 3B - }, - .cor = { - { -4.50f, { 0.0004524f, -0.0932f, 4.76f } }, // Channel 4 - { -3.61f, { 0.0002811f, -0.0659f, 3.83f } }, // Channel 5 - { 0.0f, { 0.0f, 0.0f , 0.0f } } // Channel 3B - } - }, { - .name = "NOAA-18", - .prt = { - { 1.657e-06f, 0.05090f, 276.601f }, // PRT 1 - { 1.482e-06f, 0.05101f, 276.683f }, // PRT 2 - { 1.313e-06f, 0.05117f, 276.565f }, // PRT 3 - { 1.484e-06f, 0.05103f, 276.615f } // PRT 4 - }, - .visible = { - { - .low = { 0.06174f, -2.434f }, - .high = { 0.1841f, -63.31f }, - .cutoff = 501.54f - }, { - .low = { 0.07514f, -2.960f }, - .high = { 0.2254f, -78.55f }, - .cutoff = 500.40f - } - }, - .rad = { - { 928.1460f, 0.436645f, 0.998607f }, // Channel 4 - { 833.2532f, 0.253179f, 0.999057f }, // Channel 5 - { 2659.7952f, 1.698704f, 0.996960f } // Channel 3B - }, - .cor = { - { -5.53f, { 0.00052337f, -0.11069f, 5.82f } }, // Channel 4 - { -2.22f, { 0.00017715f, -0.04360f, 2.67f } }, // Channel 5 - { 0.0f, { 0.0f, 0.0f, 0.0f } } // Channel 3B - } - }, { - .name = "NOAA-19", - .prt = { - { 1.405783e-06f, 0.051111f, 276.6067f }, // PRT 1 - { 1.496037e-06f, 0.051090f, 276.6119f }, // PRT 2 - { 1.496990e-06f, 0.051033f, 276.6311f }, // PRT 3 - { 1.493110e-06f, 0.051058f, 276.6268f } // PRT 4 - }, - .visible = { - { - .low = { 0.05555f, -2.159f }, - .high = { 0.1639f, -56.33f }, - .cutoff = 496.43f - }, { - .low = { 0.06614f, -2.565f }, - .high = { 0.1970f, -68.01f }, - .cutoff = 500.37f - } - }, - .rad = { - { 928.9f, 0.53959f, 0.998534f }, // Channel 4 - { 831.9f, 0.36064f, 0.998913f }, // Channel 5 - { 2670.0f, 1.67396f, 0.997364f } // Channel 3B - }, - .cor = { - { -5.49f, { 0.00054668f, -0.11187f, 5.70f } }, // Channel 4 - { -3.39f, { 0.00024985f, -0.05991f, 3.58f } }, // Channel 5 - { 0.0f, { 0.0f, 0.0f, 0.0f } } // Channel 3B - } - } -}; +const calibration_t calibration[3] = {{.name = "NOAA-15", + .prt = + { + {1.36328e-06f, 0.051045f, 276.60157f}, // PRT 1 + {1.47266e-06f, 0.050909f, 276.62531f}, // PRT 2 + {1.47656e-06f, 0.050907f, 276.67413f}, // PRT 3 + {1.47656e-06f, 0.050966f, 276.59258f} // PRT 4 + }, + .visible = {{.low = {0.0568f, -2.1874f}, .high = {0.1633f, -54.9928f}, .cutoff = 496.0f}, + {.low = {0.0596f, -2.4096f}, .high = {0.1629f, -55.2436f}, .cutoff = 511.0f}}, + .rad = + { + {925.4075f, 0.337810f, 0.998719f}, // Channel 4 + {839.8979f, 0.304558f, 0.999024f}, // Channel 5 + {2695.9743f, 1.621256f, 0.998015f} // Channel 3B + }, + .cor = + { + {-4.50f, {0.0004524f, -0.0932f, 4.76f}}, // Channel 4 + {-3.61f, {0.0002811f, -0.0659f, 3.83f}}, // Channel 5 + {0.0f, {0.0f, 0.0f, 0.0f}} // Channel 3B + }}, + {.name = "NOAA-18", + .prt = + { + {1.657e-06f, 0.05090f, 276.601f}, // PRT 1 + {1.482e-06f, 0.05101f, 276.683f}, // PRT 2 + {1.313e-06f, 0.05117f, 276.565f}, // PRT 3 + {1.484e-06f, 0.05103f, 276.615f} // PRT 4 + }, + .visible = {{.low = {0.06174f, -2.434f}, .high = {0.1841f, -63.31f}, .cutoff = 501.54f}, + {.low = {0.07514f, -2.960f}, .high = {0.2254f, -78.55f}, .cutoff = 500.40f}}, + .rad = + { + {928.1460f, 0.436645f, 0.998607f}, // Channel 4 + {833.2532f, 0.253179f, 0.999057f}, // Channel 5 + {2659.7952f, 1.698704f, 0.996960f} // Channel 3B + }, + .cor = + { + {-5.53f, {0.00052337f, -0.11069f, 5.82f}}, // Channel 4 + {-2.22f, {0.00017715f, -0.04360f, 2.67f}}, // Channel 5 + {0.0f, {0.0f, 0.0f, 0.0f}} // Channel 3B + }}, + {.name = "NOAA-19", + .prt = + { + {1.405783e-06f, 0.051111f, 276.6067f}, // PRT 1 + {1.496037e-06f, 0.051090f, 276.6119f}, // PRT 2 + {1.496990e-06f, 0.051033f, 276.6311f}, // PRT 3 + {1.493110e-06f, 0.051058f, 276.6268f} // PRT 4 + }, + .visible = {{.low = {0.05555f, -2.159f}, .high = {0.1639f, -56.33f}, .cutoff = 496.43f}, + {.low = {0.06614f, -2.565f}, .high = {0.1970f, -68.01f}, .cutoff = 500.37f}}, + .rad = + { + {928.9f, 0.53959f, 0.998534f}, // Channel 4 + {831.9f, 0.36064f, 0.998913f}, // Channel 5 + {2670.0f, 1.67396f, 0.997364f} // Channel 3B + }, + .cor = { + {-5.49f, {0.00054668f, -0.11187f, 5.70f}}, // Channel 4 + {-3.39f, {0.00024985f, -0.05991f, 3.58f}}, // Channel 5 + {0.0f, {0.0f, 0.0f, 0.0f}} // Channel 3B + }}}; calibration_t get_calibration(int satid) { switch (satid) { - case 15: return calibration[0]; - case 18: return calibration[1]; - case 19: return calibration[2]; - default: error("Invalid satid in get_calibration()"); /* following is only to shut up the compiler */ return calibration[0]; + case 15: + return calibration[0]; + case 18: + return calibration[1]; + case 19: + return calibration[2]; + default: + error("Invalid satid in get_calibration()"); /* following is only to shut up the compiler */ + return calibration[0]; } } diff --git a/src/calibration.h b/src/calibration.h index c41aa04..5580dd6 100644 --- a/src/calibration.h +++ b/src/calibration.h @@ -25,7 +25,7 @@ typedef struct { // Quadratics for calculating PRT temperature quadratic_t prt[4]; - + // Visible calibration coefficients struct { linear_t low; @@ -34,8 +34,8 @@ typedef struct { } visible[2]; // Radiance coefficients - struct { - float vc, A, B; + struct { + float vc, A, B; } rad[3]; // Non linear correction coefficients diff --git a/src/color.c b/src/color.c index 9d126b5..50401c4 100644 --- a/src/color.c +++ b/src/color.c @@ -1,4 +1,4 @@ -/* +/* * This file is part of Aptdec. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 * @@ -17,67 +17,61 @@ * */ -#include "apt.h" #include "color.h" -apt_rgb_t apt_applyPalette(char *palette, int val){ - return (apt_rgb_t){ - palette[(int)CLIP(val, 0, 255)*3 + 0], - palette[(int)CLIP(val, 0, 255)*3 + 1], - palette[(int)CLIP(val, 0, 255)*3 + 2] - }; +#include "apt.h" + +apt_rgb_t apt_applyPalette(char *palette, int val) { + return (apt_rgb_t){palette[(int)CLIP(val, 0, 255) * 3 + 0], palette[(int)CLIP(val, 0, 255) * 3 + 1], + palette[(int)CLIP(val, 0, 255) * 3 + 2]}; } -apt_rgb_t apt_RGBcomposite(apt_rgb_t top, float top_a, apt_rgb_t bottom, float bottom_a){ - return (apt_rgb_t){ - MCOMPOSITE(top.r, top_a, bottom.r, bottom_a), - MCOMPOSITE(top.g, top_a, bottom.g, bottom_a), - MCOMPOSITE(top.b, top_a, bottom.b, bottom_a) - }; +apt_rgb_t apt_RGBcomposite(apt_rgb_t top, float top_a, apt_rgb_t bottom, float bottom_a) { + return (apt_rgb_t){MCOMPOSITE(top.r, top_a, bottom.r, bottom_a), MCOMPOSITE(top.g, top_a, bottom.g, bottom_a), + MCOMPOSITE(top.b, top_a, bottom.b, bottom_a)}; } // The "I totally didn't just steal this from WXtoImg" palette -char apt_TempPalette[256*3] = { - "\x45\x0\x8f\x46\x0\x91\x47\x0\x92\x48\x0\x94\x49\x0\x96\x4a\x0\x98\x4b\x0\x9b\x4d\x0\x9d" - "\x4e\x0\xa0\x50\x0\xa2\x51\x0\xa5\x52\x0\xa7\x54\x0\xaa\x56\x0\xae\x57\x0\xb1" - "\x58\x0\xb4\x5a\x0\xb7\x5c\x0\xba\x5e\x0\xbd\x5f\x0\xc0\x61\x0\xc4\x64\x0\xc8" - "\x66\x0\xcb\x68\x0\xce\x69\x0\xd1\x68\x0\xd4\x65\x0\xd7\x63\x0\xda\x61\x0\xdd" - "\x5d\x0\xe1\x5b\x0\xe4\x59\x0\xe6\x56\x0\xe9\x53\x0\xeb\x50\x0\xee\x4d\x0\xf0" - "\x49\x0\xf3\x47\x0\xfc\x43\x0\xfa\x31\x0\xbf\x20\x0\x89\x20\x0\x92\x1e\x0\x95" - "\x1b\x0\x97\x19\x0\x9a\x17\x0\x9c\x15\x0\x9e\x12\x0\xa0\xf\x0\xa3\xf\x2\xa5" - "\xe\x6\xa8\xe\xa\xab\xe\xd\xad\xe\x11\xb1\xd\x15\xb4\xd\x18\xb7\xd\x1c\xba" - "\xb\x21\xbd\xa\x25\xc0\xa\x29\xc3\x9\x2d\xc6\x8\x33\xca\x7\x36\xcd\x7\x3b\xd0" - "\x7\x41\xd3\x5\x45\xd6\x4\x4b\xd9\x4\x50\xdc\x3\x55\xde\x2\x5d\xe2\x1\x61\xe5" - "\x0\x66\xe7\x0\x6c\xea\x0\x72\xec\x0\x78\xee\x0\x7d\xf0\x0\x82\xf3\x0\x8d\xfc" - "\x0\x90\xfa\x0\x71\xbf\x0\x54\x89\x0\x5c\x91\x0\x61\x94\x0\x64\x96\x0\x68\x97" - "\x0\x6d\x99\x0\x71\x9b\x0\x75\x9d\x0\x79\x9f\x0\x7e\xa0\x0\x82\xa2\x0\x87\xa4" - "\x0\x8c\xa6\x0\x92\xaa\x0\x96\xac\x0\x9c\xae\x0\xa2\xb1\x0\xa6\xb3\x0\xaa\xb5" - "\x0\xad\xb7\x0\xb1\xba\x0\xb6\xbe\x0\xba\xc0\x0\xbe\xc2\x0\xc2\xc5\x0\xc6\xc6" - "\x0\xca\xc9\x0\xcc\xca\x0\xcf\xcb\x0\xd2\xcc\x0\xd4\xcc\x0\xd6\xcc\x0\xd9\xcb" - "\x0\xdb\xcb\x0\xde\xcb\x0\xe0\xcb\x0\xe2\xcc\x0\xea\xd2\x0\xea\xcf\x0\xb9\xa4" - "\x0\x8e\x7a\x1\x94\x7c\x4\x97\x79\x7\x99\x75\x9\x9b\x71\xd\x9d\x6b\x10\x9f\x67" - "\x12\xa1\x63\x15\xa3\x5f\x17\xa5\x59\x1a\xa8\x55\x1d\xaa\x50\x20\xac\x4b\x24\xaf\x45" - "\x28\xb2\x41\x2b\xb5\x3b\x2e\xb8\x35\x31\xba\x30\x34\xbd\x2b\x39\xbf\x24\x3f\xc1\x17" - "\x49\xc5\x8\x4f\xc8\x1\x4f\xca\x0\x4e\xcd\x0\x4e\xcf\x0\x4f\xd2\x0\x54\xd5\x0" - "\x5d\xd8\x0\x68\xdb\x0\x6e\xdd\x0\x74\xdf\x0\x7a\xe2\x0\x7f\xe4\x0\x85\xe7\x0" - "\x8b\xe9\x0\x8f\xeb\x0\x9b\xf3\x0\x9e\xf2\x0\x7e\xbb\x0\x60\x8a\x0\x68\x92\x0" - "\x6d\x95\x0\x71\x96\x0\x75\x98\x0\x7b\x9a\x0\x7f\x9d\x0\x83\x9f\x0\x87\xa1\x0" - "\x8c\xa2\x0\x8f\xa5\x0\x92\xa7\x0\x96\xa9\x0\x9a\xad\x0\x9d\xb0\x0\xa1\xb2\x0" - "\xa5\xb5\x0\xa9\xb7\x0\xad\xba\x0\xb2\xbd\x0\xb6\xbf\x0\xbb\xc3\x0\xbf\xc6\x0" - "\xc3\xc8\x0\xc8\xcb\x0\xcc\xce\x0\xd0\xd1\x0\xd3\xd2\x0\xd5\xd4\x0\xd9\xd4\x0" - "\xdc\xd4\x0\xde\xd5\x0\xe1\xd5\x0\xe3\xd5\x0\xe6\xd4\x0\xe8\xd1\x0\xea\xce\x0" - "\xf2\xcf\x0\xf2\xca\x0\xbb\x99\x0\x8a\x6e\x0\x92\x72\x0\x95\x72\x0\x97\x71\x0" - "\x9a\x70\x0\x9c\x6e\x0\x9e\x6d\x0\xa0\x6b\x0\xa3\x6a\x0\xa5\x68\x0\xa8\x67\x0" - "\xab\x66\x0\xae\x65\x0\xb2\x63\x0\xb4\x61\x0\xb7\x5f\x0\xba\x5d\x0\xbd\x5c\x0" - "\xc0\x59\x0\xc3\x57\x0\xc6\x54\x0\xca\x50\x0\xcd\x4d\x0\xd0\x4a\x0\xd3\x47\x0" - "\xd6\x43\x0\xd9\x40\x0\xdc\x3d\x0\xde\x39\x0\xe2\x33\x0\xe5\x2f\x0\xe7\x2c\x0" - "\xea\x28\x0\xec\x23\x0\xef\x1f\x0\xf1\x1a\x0\xf3\x14\x0\xfb\xf\x0\xfa\xd\x0" - "\xc1\x5\x0\x8e\x0\x0\x97\x0\x0\x9b\x0\x0\x9e\x0\x0\xa1\x0\x0\xa5\x0\x0" - "\xa9\x0\x0\xad\x0\x0\xb1\x0\x0\xb6\x0\x0\xba\x0\x0\xbd\x0\x0\xc2\x0\x0" - "\xc8\x0\x0\xcc\x0\x0\xcc\x0\x0" -}; +char apt_TempPalette[256 * 3] = { + "\x45\x0\x8f\x46\x0\x91\x47\x0\x92\x48\x0\x94\x49\x0\x96\x4a\x0\x98\x4b\x0\x9b\x4d\x0\x9d" + "\x4e\x0\xa0\x50\x0\xa2\x51\x0\xa5\x52\x0\xa7\x54\x0\xaa\x56\x0\xae\x57\x0\xb1" + "\x58\x0\xb4\x5a\x0\xb7\x5c\x0\xba\x5e\x0\xbd\x5f\x0\xc0\x61\x0\xc4\x64\x0\xc8" + "\x66\x0\xcb\x68\x0\xce\x69\x0\xd1\x68\x0\xd4\x65\x0\xd7\x63\x0\xda\x61\x0\xdd" + "\x5d\x0\xe1\x5b\x0\xe4\x59\x0\xe6\x56\x0\xe9\x53\x0\xeb\x50\x0\xee\x4d\x0\xf0" + "\x49\x0\xf3\x47\x0\xfc\x43\x0\xfa\x31\x0\xbf\x20\x0\x89\x20\x0\x92\x1e\x0\x95" + "\x1b\x0\x97\x19\x0\x9a\x17\x0\x9c\x15\x0\x9e\x12\x0\xa0\xf\x0\xa3\xf\x2\xa5" + "\xe\x6\xa8\xe\xa\xab\xe\xd\xad\xe\x11\xb1\xd\x15\xb4\xd\x18\xb7\xd\x1c\xba" + "\xb\x21\xbd\xa\x25\xc0\xa\x29\xc3\x9\x2d\xc6\x8\x33\xca\x7\x36\xcd\x7\x3b\xd0" + "\x7\x41\xd3\x5\x45\xd6\x4\x4b\xd9\x4\x50\xdc\x3\x55\xde\x2\x5d\xe2\x1\x61\xe5" + "\x0\x66\xe7\x0\x6c\xea\x0\x72\xec\x0\x78\xee\x0\x7d\xf0\x0\x82\xf3\x0\x8d\xfc" + "\x0\x90\xfa\x0\x71\xbf\x0\x54\x89\x0\x5c\x91\x0\x61\x94\x0\x64\x96\x0\x68\x97" + "\x0\x6d\x99\x0\x71\x9b\x0\x75\x9d\x0\x79\x9f\x0\x7e\xa0\x0\x82\xa2\x0\x87\xa4" + "\x0\x8c\xa6\x0\x92\xaa\x0\x96\xac\x0\x9c\xae\x0\xa2\xb1\x0\xa6\xb3\x0\xaa\xb5" + "\x0\xad\xb7\x0\xb1\xba\x0\xb6\xbe\x0\xba\xc0\x0\xbe\xc2\x0\xc2\xc5\x0\xc6\xc6" + "\x0\xca\xc9\x0\xcc\xca\x0\xcf\xcb\x0\xd2\xcc\x0\xd4\xcc\x0\xd6\xcc\x0\xd9\xcb" + "\x0\xdb\xcb\x0\xde\xcb\x0\xe0\xcb\x0\xe2\xcc\x0\xea\xd2\x0\xea\xcf\x0\xb9\xa4" + "\x0\x8e\x7a\x1\x94\x7c\x4\x97\x79\x7\x99\x75\x9\x9b\x71\xd\x9d\x6b\x10\x9f\x67" + "\x12\xa1\x63\x15\xa3\x5f\x17\xa5\x59\x1a\xa8\x55\x1d\xaa\x50\x20\xac\x4b\x24\xaf\x45" + "\x28\xb2\x41\x2b\xb5\x3b\x2e\xb8\x35\x31\xba\x30\x34\xbd\x2b\x39\xbf\x24\x3f\xc1\x17" + "\x49\xc5\x8\x4f\xc8\x1\x4f\xca\x0\x4e\xcd\x0\x4e\xcf\x0\x4f\xd2\x0\x54\xd5\x0" + "\x5d\xd8\x0\x68\xdb\x0\x6e\xdd\x0\x74\xdf\x0\x7a\xe2\x0\x7f\xe4\x0\x85\xe7\x0" + "\x8b\xe9\x0\x8f\xeb\x0\x9b\xf3\x0\x9e\xf2\x0\x7e\xbb\x0\x60\x8a\x0\x68\x92\x0" + "\x6d\x95\x0\x71\x96\x0\x75\x98\x0\x7b\x9a\x0\x7f\x9d\x0\x83\x9f\x0\x87\xa1\x0" + "\x8c\xa2\x0\x8f\xa5\x0\x92\xa7\x0\x96\xa9\x0\x9a\xad\x0\x9d\xb0\x0\xa1\xb2\x0" + "\xa5\xb5\x0\xa9\xb7\x0\xad\xba\x0\xb2\xbd\x0\xb6\xbf\x0\xbb\xc3\x0\xbf\xc6\x0" + "\xc3\xc8\x0\xc8\xcb\x0\xcc\xce\x0\xd0\xd1\x0\xd3\xd2\x0\xd5\xd4\x0\xd9\xd4\x0" + "\xdc\xd4\x0\xde\xd5\x0\xe1\xd5\x0\xe3\xd5\x0\xe6\xd4\x0\xe8\xd1\x0\xea\xce\x0" + "\xf2\xcf\x0\xf2\xca\x0\xbb\x99\x0\x8a\x6e\x0\x92\x72\x0\x95\x72\x0\x97\x71\x0" + "\x9a\x70\x0\x9c\x6e\x0\x9e\x6d\x0\xa0\x6b\x0\xa3\x6a\x0\xa5\x68\x0\xa8\x67\x0" + "\xab\x66\x0\xae\x65\x0\xb2\x63\x0\xb4\x61\x0\xb7\x5f\x0\xba\x5d\x0\xbd\x5c\x0" + "\xc0\x59\x0\xc3\x57\x0\xc6\x54\x0\xca\x50\x0\xcd\x4d\x0\xd0\x4a\x0\xd3\x47\x0" + "\xd6\x43\x0\xd9\x40\x0\xdc\x3d\x0\xde\x39\x0\xe2\x33\x0\xe5\x2f\x0\xe7\x2c\x0" + "\xea\x28\x0\xec\x23\x0\xef\x1f\x0\xf1\x1a\x0\xf3\x14\x0\xfb\xf\x0\xfa\xd\x0" + "\xc1\x5\x0\x8e\x0\x0\x97\x0\x0\x9b\x0\x0\x9e\x0\x0\xa1\x0\x0\xa5\x0\x0" + "\xa9\x0\x0\xad\x0\x0\xb1\x0\x0\xb6\x0\x0\xba\x0\x0\xbd\x0\x0\xc2\x0\x0" + "\xc8\x0\x0\xcc\x0\x0\xcc\x0\x0"}; -char apt_PrecipPalette[58*3] = { +char apt_PrecipPalette[58 * 3] = { "\x8\x89\x41\x0\xc5\x44\x0\xd1\x2c\x0\xe3\x1c\x0\xf9\x6\x14\xff\x0\x3e\xff\x0\x5d\xff\x0" "\x80\xff\x0\xab\xff\x0\xcd\xfe\x0\xf8\xff\x0\xff\xe6\x0\xff\xb8\x0\xff\x98\x0" "\xff\x75\x0\xff\x49\x0\xfe\x26\x0\xff\x4\x0\xdf\x0\x0\xa8\x0\x0\x87\x0\x0" @@ -86,5 +80,4 @@ char apt_PrecipPalette[58*3] = { "\xc0\xc0\xc0\xce\xce\xce\xdc\xdc\xdc\xef\xef\xef\xfa\xfa\xfa\xff\xff\xff\xff\xff\xff" "\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff" "\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff" - "\xff\xff\xff" -}; + "\xff\xff\xff"}; diff --git a/src/color.h b/src/color.h index 6441fb2..1fc253e 100644 --- a/src/color.h +++ b/src/color.h @@ -1,3 +1,3 @@ #include "common.h" -#define MCOMPOSITE(m1, a1, m2, a2) (m1*a1 + m2*a2*(1-a1)) +#define MCOMPOSITE(m1, a1, m2, a2) (m1 * a1 + m2 * a2 * (1 - a1)) diff --git a/src/common.h b/src/common.h index ed253bc..7130d46 100644 --- a/src/common.h +++ b/src/common.h @@ -1,4 +1,4 @@ -/* +/* * This file is part of Aptdec. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 * @@ -22,39 +22,39 @@ // Useful macros #define CLIP(v, lo, hi) (v > hi ? hi : (v > lo ? v : lo)) -#define CONTAINS(str, char) (strchr(str, (int) char) != NULL) +#define CONTAINS(str, char) (strchr(str, (int)char) != NULL) // Typedefs #ifndef STRUCTS_DEFINED #define STRUCTS_DEFINED typedef struct { - char *type; // Output image type - char *effects; // Effects on the image - int satnum; // The satellite number - char *path; // Output directory - int realtime; // Realtime decoding - char *filename; // Output filename - char *palette; // Filename of palette - float gamma; // Gamma + char *type; // Output image type + char *effects; // Effects on the image + int satnum; // The satellite number + char *path; // Output directory + int realtime; // Realtime decoding + char *filename; // Output filename + char *palette; // Filename of palette + float gamma; // Gamma } options_t; enum imagetypes { - Raw_Image='r', - Palleted='p', - Temperature='t', - Channel_A='a', - Channel_B='b', - Distribution='d', - Visible='v' + Raw_Image = 'r', + Palleted = 'p', + Temperature = 't', + Channel_A = 'a', + Channel_B = 'b', + Distribution = 'd', + Visible = 'v' }; enum effects { - Crop_Telemetry='t', - Histogram_Equalise='h', - Denoise='d', - Precipitation_Overlay='p', - Flip_Image='f', - Linear_Equalise='l', - Crop_Noise='c' + Crop_Telemetry = 't', + Histogram_Equalise = 'h', + Denoise = 'd', + Precipitation_Overlay = 'p', + Flip_Image = 'f', + Linear_Equalise = 'l', + Crop_Noise = 'c' }; #endif diff --git a/src/dsp.c b/src/dsp.c old mode 100755 new mode 100644 index 62ff148..befffe7 --- a/src/dsp.c +++ b/src/dsp.c @@ -16,10 +16,10 @@ * along with this program. If not, see . */ -#include +#include #include +#include #include -#include #include "apt.h" #include "filter.h" @@ -34,7 +34,7 @@ #define MAX_CARRIER_OFFSET 20.0 #define RSMULT 15 -#define Fi (APT_IMG_WIDTH*2 * RSMULT) +#define Fi (APT_IMG_WIDTH * 2 * RSMULT) static float _sample_rate; @@ -47,218 +47,212 @@ static float pll_beta; // Initalise and configure PLL int apt_init(double sample_rate) { - if(sample_rate > Fi) return 1; - if(sample_rate < APT_IMG_WIDTH*2) return -1; - _sample_rate = sample_rate; + if (sample_rate > Fi) return 1; + if (sample_rate < APT_IMG_WIDTH * 2) return -1; + _sample_rate = sample_rate; - // Pll configuration - pll_alpha = 50 / _sample_rate; - pll_beta = pll_alpha * pll_alpha / 2.0; - oscillator_freq = CARRIER_FREQ/sample_rate; + // Pll configuration + pll_alpha = 50 / _sample_rate; + pll_beta = pll_alpha * pll_alpha / 2.0; + oscillator_freq = CARRIER_FREQ / sample_rate; - return 0; + return 0; } static float pll(complexf_t in) { - static float oscillator_phase = 0.0; + static float oscillator_phase = 0.0; - // Internal oscillator + // Internal oscillator #ifdef _MSC_VER - complexf_t osc = _FCbuild(cos(oscillator_phase), -sin(oscillator_phase)); + complexf_t osc = _FCbuild(cos(oscillator_phase), -sin(oscillator_phase)); in = _FCmulcc(in, osc); #else - complexf_t osc = cos(oscillator_phase) + -sin(oscillator_phase)*I; + complexf_t osc = cos(oscillator_phase) + -sin(oscillator_phase) * I; in *= osc; #endif - // Error detector - float error = cargf(in); + // Error detector + float error = cargf(in); - // Adjust frequency and phase - oscillator_freq += pll_beta*error; - oscillator_freq = clamp_half(oscillator_freq, (CARRIER_FREQ + MAX_CARRIER_OFFSET) / _sample_rate); - oscillator_phase += M_TAUf * (pll_alpha*error + oscillator_freq); - oscillator_phase = remainderf(oscillator_phase, M_TAUf); + // Adjust frequency and phase + oscillator_freq += pll_beta * error; + oscillator_freq = clamp_half(oscillator_freq, (CARRIER_FREQ + MAX_CARRIER_OFFSET) / _sample_rate); + oscillator_phase += M_TAUf * (pll_alpha * error + oscillator_freq); + oscillator_phase = remainderf(oscillator_phase, M_TAUf); - return crealf(in); + return crealf(in); } // Convert samples into pixels static int getamp(float *ampbuff, int count, apt_getsamples_t getsamples, void *context) { - static float inbuff[BLKIN]; - static int idxin = 0; - static int nin = 0; - - for (int n = 0; n < count; n++) { - float I2, Q; - - // Get some more samples when needed - if (nin < HILBERT_FILTER_SIZE * 2 + 2) { - // Number of samples read - int res; - memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); - idxin = 0; - - // Read some samples - res = getsamples(context, &(inbuff[nin]), BLKIN - nin); - nin += res; - - // Make sure there is enough samples to continue - if (nin < HILBERT_FILTER_SIZE * 2 + 2) - return n; - } - - // Process read samples into a brightness value - complexf_t sample = hilbert_transform(&inbuff[idxin], hilbert_filter, HILBERT_FILTER_SIZE); - ampbuff[n] = pll(sample); - - // Increment current sample - idxin++; - nin--; - } - - return count; + static float inbuff[BLKIN]; + static int idxin = 0; + static int nin = 0; + + for (int n = 0; n < count; n++) { + float I2, Q; + + // Get some more samples when needed + if (nin < HILBERT_FILTER_SIZE * 2 + 2) { + // Number of samples read + int res; + memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); + idxin = 0; + + // Read some samples + res = getsamples(context, &(inbuff[nin]), BLKIN - nin); + nin += res; + + // Make sure there is enough samples to continue + if (nin < HILBERT_FILTER_SIZE * 2 + 2) return n; + } + + // Process read samples into a brightness value + complexf_t sample = hilbert_transform(&inbuff[idxin], hilbert_filter, HILBERT_FILTER_SIZE); + ampbuff[n] = pll(sample); + + // Increment current sample + idxin++; + nin--; + } + + return count; } -// Sub-pixel offsetting +// Sub-pixel offsetting int getpixelv(float *pvbuff, int count, apt_getsamples_t getsamples, void *context) { - // Amplitude buffer - static float ampbuff[BLKAMP]; - static int nam = 0; - static int idxam = 0; + // Amplitude buffer + static float ampbuff[BLKAMP]; + static int nam = 0; + static int idxam = 0; - float mult; + float mult; - // Gaussian resampling factor - mult = (float) Fi / _sample_rate * FreqLine; - int m = (int)(LOW_PASS_SIZE / mult + 1); + // Gaussian resampling factor + mult = (float)Fi / _sample_rate * FreqLine; + int m = (int)(LOW_PASS_SIZE / mult + 1); - for (int n = 0; n < count; n++) { - int shift; + for (int n = 0; n < count; n++) { + int shift; - if (nam < m) { - int res; - memmove(ampbuff, &(ampbuff[idxam]), nam * sizeof(float)); - idxam = 0; - res = getamp(&(ampbuff[nam]), BLKAMP - nam, getsamples, context); - nam += res; - if (nam < m) - return n; - } + if (nam < m) { + int res; + memmove(ampbuff, &(ampbuff[idxam]), nam * sizeof(float)); + idxam = 0; + res = getamp(&(ampbuff[nam]), BLKAMP - nam, getsamples, context); + nam += res; + if (nam < m) return n; + } - pvbuff[n] = interpolating_convolve(&(ampbuff[idxam]), low_pass, LOW_PASS_SIZE, offset, mult) * mult * 256.0; + pvbuff[n] = interpolating_convolve(&(ampbuff[idxam]), low_pass, LOW_PASS_SIZE, offset, mult) * mult * 256.0; - shift = ((int) floor((RSMULT - offset) / mult)) + 1; - offset = shift * mult + offset - RSMULT; + shift = ((int)floor((RSMULT - offset) / mult)) + 1; + offset = shift * mult + offset - RSMULT; - idxam += shift; - nam -= shift; - } + idxam += shift; + nam -= shift; + } - return count; + return count; } // Get an entire row of pixels, aligned with sync markers int apt_getpixelrow(float *pixelv, int nrow, int *zenith, int reset, apt_getsamples_t getsamples, void *context) { - static float pixels[APT_IMG_WIDTH + SYNC_PATTERN_SIZE]; - static size_t npv; - static int synced = 0; - static float max = 0.0; - static float minDoppler = 1000000000, previous = 0; - - if(reset) synced = 0; - - float corr, ecorr, lcorr; - int res; - - // Move the row buffer into the the image buffer - if (npv > 0) - memmove(pixelv, pixels, npv * sizeof(float)); - - // Get the sync line - if (npv < SYNC_PATTERN_SIZE + 2) { - res = getpixelv(&(pixelv[npv]), SYNC_PATTERN_SIZE + 2 - npv, getsamples, context); - npv += res; - if (npv < SYNC_PATTERN_SIZE + 2) - return 0; - } - - // Calculate the frequency offset - ecorr = convolve(pixelv, sync_pattern, SYNC_PATTERN_SIZE); - corr = convolve(&pixelv[1], sync_pattern, SYNC_PATTERN_SIZE - 1); - lcorr = convolve(&pixelv[2], sync_pattern, SYNC_PATTERN_SIZE - 2); - FreqLine = 1.0+((ecorr-lcorr) / corr / APT_IMG_WIDTH / 4.0); - - float val = fabs(lcorr - ecorr)*0.25 + previous*0.75; - if(val < minDoppler && nrow > 10){ - minDoppler = val; - *zenith = nrow; - } - previous = fabs(lcorr - ecorr); - - // The point in which the pixel offset is recalculated - if (corr < 0.75 * max) { - synced = 0; - FreqLine = 1.0; - } - max = corr; - - if (synced < 8) { - int mshift; + static float pixels[APT_IMG_WIDTH + SYNC_PATTERN_SIZE]; + static size_t npv; + static int synced = 0; + static float max = 0.0; + static float minDoppler = 1000000000, previous = 0; + + if (reset) synced = 0; + + float corr, ecorr, lcorr; + int res; + + // Move the row buffer into the the image buffer + if (npv > 0) memmove(pixelv, pixels, npv * sizeof(float)); + + // Get the sync line + if (npv < SYNC_PATTERN_SIZE + 2) { + res = getpixelv(&(pixelv[npv]), SYNC_PATTERN_SIZE + 2 - npv, getsamples, context); + npv += res; + if (npv < SYNC_PATTERN_SIZE + 2) return 0; + } + + // Calculate the frequency offset + ecorr = convolve(pixelv, sync_pattern, SYNC_PATTERN_SIZE); + corr = convolve(&pixelv[1], sync_pattern, SYNC_PATTERN_SIZE - 1); + lcorr = convolve(&pixelv[2], sync_pattern, SYNC_PATTERN_SIZE - 2); + FreqLine = 1.0 + ((ecorr - lcorr) / corr / APT_IMG_WIDTH / 4.0); + + float val = fabs(lcorr - ecorr) * 0.25 + previous * 0.75; + if (val < minDoppler && nrow > 10) { + minDoppler = val; + *zenith = nrow; + } + previous = fabs(lcorr - ecorr); + + // The point in which the pixel offset is recalculated + if (corr < 0.75 * max) { + synced = 0; + FreqLine = 1.0; + } + max = corr; + + if (synced < 8) { + int mshift; static int lastmshift; - if (npv < APT_IMG_WIDTH + SYNC_PATTERN_SIZE) { - res = getpixelv(&(pixelv[npv]), APT_IMG_WIDTH + SYNC_PATTERN_SIZE - npv, getsamples, context); - npv += res; - if (npv < APT_IMG_WIDTH + SYNC_PATTERN_SIZE) - return 0; - } - - // Test every possible position until we get the best result - mshift = 0; - for (int shift = 0; shift < APT_IMG_WIDTH; shift++) { - float corr; - - corr = convolve(&(pixelv[shift + 1]), sync_pattern, SYNC_PATTERN_SIZE); - if (corr > max) { - mshift = shift; - max = corr; - } - } - - // Stop rows dissapearing into the void - int mshiftOrig = mshift; - if(abs(lastmshift-mshift) > 3 && nrow != 0){ - mshift = 0; - } - lastmshift = mshiftOrig; - - // If we are already as aligned as we can get, just continue - if (mshift == 0) { - synced++; - } else { - memmove(pixelv, &(pixelv[mshift]), (npv - mshift) * sizeof(float)); - npv -= mshift; - synced = 0; - FreqLine = 1.0; - } - } - - // Get the rest of this row - if (npv < APT_IMG_WIDTH) { - res = getpixelv(&(pixelv[npv]), APT_IMG_WIDTH - npv, getsamples, context); - npv += res; - if (npv < APT_IMG_WIDTH) - return 0; - } - - // Move the sync lines into the output buffer with the calculated offset - if (npv == APT_IMG_WIDTH) { - npv = 0; - } else { - memmove(pixels, &(pixelv[APT_IMG_WIDTH]), (npv - APT_IMG_WIDTH) * sizeof(float)); - npv -= APT_IMG_WIDTH; - } - - return 1; + if (npv < APT_IMG_WIDTH + SYNC_PATTERN_SIZE) { + res = getpixelv(&(pixelv[npv]), APT_IMG_WIDTH + SYNC_PATTERN_SIZE - npv, getsamples, context); + npv += res; + if (npv < APT_IMG_WIDTH + SYNC_PATTERN_SIZE) return 0; + } + + // Test every possible position until we get the best result + mshift = 0; + for (int shift = 0; shift < APT_IMG_WIDTH; shift++) { + float corr; + + corr = convolve(&(pixelv[shift + 1]), sync_pattern, SYNC_PATTERN_SIZE); + if (corr > max) { + mshift = shift; + max = corr; + } + } + + // Stop rows dissapearing into the void + int mshiftOrig = mshift; + if (abs(lastmshift - mshift) > 3 && nrow != 0) { + mshift = 0; + } + lastmshift = mshiftOrig; + + // If we are already as aligned as we can get, just continue + if (mshift == 0) { + synced++; + } else { + memmove(pixelv, &(pixelv[mshift]), (npv - mshift) * sizeof(float)); + npv -= mshift; + synced = 0; + FreqLine = 1.0; + } + } + + // Get the rest of this row + if (npv < APT_IMG_WIDTH) { + res = getpixelv(&(pixelv[npv]), APT_IMG_WIDTH - npv, getsamples, context); + npv += res; + if (npv < APT_IMG_WIDTH) return 0; + } + + // Move the sync lines into the output buffer with the calculated offset + if (npv == APT_IMG_WIDTH) { + npv = 0; + } else { + memmove(pixels, &(pixelv[APT_IMG_WIDTH]), (npv - APT_IMG_WIDTH) * sizeof(float)); + npv -= APT_IMG_WIDTH; + } + + return 1; } diff --git a/src/filter.c b/src/filter.c old mode 100755 new mode 100644 index e3adcb3..06746aa --- a/src/filter.c +++ b/src/filter.c @@ -16,45 +16,47 @@ * along with this program. If not, see . */ -#include #include "filter.h" + +#include + #include "util.h" float convolve(const float *in, const float *taps, size_t len) { - float sum = 0.0; - for (size_t i = 0; i < len; i++) { - sum += in[i] * taps[i]; - } + float sum = 0.0; + for (size_t i = 0; i < len; i++) { + sum += in[i] * taps[i]; + } - return sum; + return sum; } complexf_t hilbert_transform(const float *in, const float *taps, size_t len) { - float i = 0.0; - float q = 0.0; + float i = 0.0; + float q = 0.0; - for (size_t k = 0; k < len; k++) { - q += in[2*k] * taps[k]; - i += in[2*k]; - } + for (size_t k = 0; k < len; k++) { + q += in[2 * k] * taps[k]; + i += in[2 * k]; + } - i = in[len-1] - (i / len); + i = in[len - 1] - (i / len); #ifdef _MSC_VER - return _FCbuild(i, q); + return _FCbuild(i, q); #else - return i + q*I; + return i + q * I; #endif } float interpolating_convolve(const float *in, const float *taps, size_t len, float offset, float delta) { - float out = 0.0; - float n = offset; + float out = 0.0; + float n = offset; - for (size_t i = 0; i < (len-1)/delta-1; n += delta, i++) { - int k = (int)floor(n); - float alpha = n - k; + for (size_t i = 0; i < (len - 1) / delta - 1; n += delta, i++) { + int k = (int)floor(n); + float alpha = n - k; - out += in[i] * (taps[k] * (1.0f-alpha) + taps[k + 1] * alpha); - } - return out; + out += in[i] * (taps[k] * (1.0f - alpha) + taps[k + 1] * alpha); + } + return out; } diff --git a/src/filter.h b/src/filter.h old mode 100755 new mode 100644 index 6c290f0..124a59b --- a/src/filter.h +++ b/src/filter.h @@ -16,8 +16,8 @@ * along with this program. If not, see . */ -#include #include +#include #ifdef _MSC_VER typedef _Fcomplex complexf_t; @@ -25,6 +25,6 @@ typedef _Fcomplex complexf_t; typedef complex float complexf_t; #endif -float convolve(const float *in, const float *taps, size_t len); +float convolve(const float *in, const float *taps, size_t len); complexf_t hilbert_transform(const float *in, const float *taps, size_t len); -float interpolating_convolve(const float *in, const float *taps, size_t len, float offset, float delta); +float interpolating_convolve(const float *in, const float *taps, size_t len, float offset, float delta); diff --git a/src/image.c b/src/image.c index 3a8e4c0..8c9cbc9 100644 --- a/src/image.c +++ b/src/image.c @@ -1,4 +1,4 @@ -/* +/* * This file is part of Aptdec. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 * @@ -17,382 +17,372 @@ * */ -#include -#include +#include "image.h" + #include +#include #include +#include -#include "apt.h" #include "algebra.h" -#include "image.h" +#include "apt.h" #include "util.h" 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 float teleramp[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; + // { 0.106, 0.215, 0.324, 0.433, 0.542, 0.652, 0.78, 0.87, 0.0 } + const float teleramp[9] = {31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0}; - return linear_regression(wedges, teleramp, 9); + return linear_regression(wedges, teleramp, 9); } static float tele[16]; static float Cs; -void apt_histogramEqualise(float **prow, int nrow, int offset, int width){ - // Plot histogram - int histogram[256] = { 0 }; - for(int y = 0; y < nrow; y++) - for(int x = 0; x < width; x++) - histogram[(int)CLIP(prow[y][x+offset], 0, 255)]++; - - // Calculate cumulative frequency - long sum = 0, cf[256] = { 0 }; - for(int i = 0; i < 255; i++){ - sum += histogram[i]; - cf[i] = sum; - } - - // Apply histogram - int area = nrow * width; - for(int y = 0; y < nrow; y++){ - for(int x = 0; x < width; x++){ - int k = (int)prow[y][x+offset]; - prow[y][x+offset] = (256.0f/area) * cf[k]; - } - } +void apt_histogramEqualise(float **prow, int nrow, int offset, int width) { + // Plot histogram + int histogram[256] = {0}; + for (int y = 0; y < nrow; y++) + for (int x = 0; x < width; x++) histogram[(int)CLIP(prow[y][x + offset], 0, 255)]++; + + // Calculate cumulative frequency + long sum = 0, cf[256] = {0}; + for (int i = 0; i < 255; i++) { + sum += histogram[i]; + cf[i] = sum; + } + + // Apply histogram + int area = nrow * width; + for (int y = 0; y < nrow; y++) { + for (int x = 0; x < width; x++) { + int k = (int)prow[y][x + offset]; + prow[y][x + offset] = (256.0f / area) * cf[k]; + } + } } -void apt_linearEnhance(float **prow, int nrow, int offset, int width){ - // Plot histogram - int histogram[256] = { 0 }; - for(int y = 0; y < nrow; y++) - for(int x = 0; x < width; x++) - histogram[(int)CLIP(prow[y][x+offset], 0, 255)]++; - - // Find min/max points - int min = -1, max = -1; - for(int i = 5; i < 250; i++){ - if(histogram[i]/width/(nrow/255.0) > 0.1){ - if(min == -1) min = i; - max = i; - } - } - - // Stretch the brightness into the new range - for(int y = 0; y < nrow; y++){ - for(int x = 0; x < width; x++){ - prow[y][x+offset] = (prow[y][x+offset]-min) / (max-min) * 255.0f; - prow[y][x+offset] = CLIP(prow[y][x+offset], 0.0f, 255.0f); - } - } +void apt_linearEnhance(float **prow, int nrow, int offset, int width) { + // Plot histogram + int histogram[256] = {0}; + for (int y = 0; y < nrow; y++) + for (int x = 0; x < width; x++) histogram[(int)CLIP(prow[y][x + offset], 0, 255)]++; + + // Find min/max points + int min = -1, max = -1; + for (int i = 5; i < 250; i++) { + if (histogram[i] / width / (nrow / 255.0) > 0.1) { + if (min == -1) min = i; + max = i; + } + } + + // Stretch the brightness into the new range + for (int y = 0; y < nrow; y++) { + for (int x = 0; x < width; x++) { + prow[y][x + offset] = (prow[y][x + offset] - min) / (max - min) * 255.0f; + prow[y][x + offset] = CLIP(prow[y][x + offset], 0.0f, 255.0f); + } + } } // Brightness calibrate, including telemetry -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 = linear_calc(prow[y][x + offset], regr); - prow[y][x + offset] = CLIP(pv, 0, 255); - } - } +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 = linear_calc(prow[y][x + offset], regr); + prow[y][x + offset] = CLIP(pv, 0, 255); + } + } } -float teleNoise(float *wedges){ - float pattern[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; - float noise = 0; - for(int i = 0; i < 9; i++) - noise += fabs(wedges[i] - pattern[i]); +float teleNoise(float *wedges) { + float pattern[9] = {31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0}; + float noise = 0; + for (int i = 0; i < 9; i++) noise += fabs(wedges[i] - pattern[i]); - return noise; + return noise; } // Get telemetry data for thermal calibration apt_channel_t apt_calibrate(float **prow, int nrow, int offset, int width) { - float teleline[APT_MAX_HEIGHT] = { 0.0 }; - float wedge[16]; - linear_t regr[APT_MAX_HEIGHT/APT_FRAME_LEN + 1]; - int telestart, mtelestart = 0; - int channel = -1; - - // The minimum rows required to decode a full frame - if (nrow < APT_CALIBRATION_ROWS) { - error_noexit("Telemetry decoding error, not enough rows"); - return APT_CHANNEL_UNKNOWN; - } - - // Calculate average of a row of telemetry - for (int y = 0; y < nrow; y++) { - for (int x = 3; x < 43; x++) - teleline[y] += prow[y][x + offset + width]; - - teleline[y] /= 40.0; - } - - /* Wedge 7 is white and 8 is black, this will have the largest - * difference in brightness, this can be used to find the current - * position within the telemetry. - */ - float max = 0.0f; - for (int n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) { - float df; - - // (sum 4px below) - (sum 4px above) - df = (float)((teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + teleline[n - 1]) - - (teleline[n + 0] + teleline[n + 1] + teleline[n + 2] + teleline[n + 3])); - - // Find the maximum difference - if (df > max) { - mtelestart = n; - max = df; - } - } - - telestart = (mtelestart + 64) % APT_FRAME_LEN; - - // Make sure that theres at least one full frame in the image - if (nrow < telestart + APT_FRAME_LEN) { - error_noexit("Telemetry decoding error, not enough rows"); - return APT_CHANNEL_UNKNOWN; - } - - // Find the least noisy frame - float minNoise = -1; - int bestFrame = -1; - for (int n = telestart, k = 0; n < nrow - APT_FRAME_LEN; n += APT_FRAME_LEN, k++) { - int j; - - for (j = 0; j < 16; j++) { - int i; - - wedge[j] = 0.0; - for (i = 1; i < 7; i++) - wedge[j] += teleline[n + j * 8 + i]; - wedge[j] /= 6; - } - - float noise = teleNoise(wedge); - if(noise < minNoise || minNoise == -1){ - minNoise = noise; - bestFrame = k; - - // Compute & apply regression on the wedges - regr[k] = compute_regression(wedge); - for (int j = 0; j < 16; j++) - tele[j] = linear_calc(wedge[j], regr[k]); - - /* Compare the channel ID wedge to the reference - * wedges, the wedge with the closest match will - * be the channel ID - */ - float min = -1; - for (int j = 0; j < 6; j++) { - float df = (float)(tele[15] - tele[j]); - df *= df; - - if (df < min || min == -1) { - channel = j; - min = df; - } - } - - // Find the brightness of the minute marker, I don't really know what for - Cs = 0.0; - int i, j = n; - for (i = 0, j = n; j < n + APT_FRAME_LEN; j++) { - float csline = 0.0; - for (int l = 3; l < 43; l++) - csline += prow[n][l + offset - APT_SPC_WIDTH]; - csline /= 40.0; - - if (csline > 50.0) { - Cs += csline; - i++; - } - } - Cs = linear_calc((Cs / i), regr[k]); - } - } - - if(bestFrame == -1){ - error_noexit("Something has gone very wrong, please file a bug report"); - return APT_CHANNEL_UNKNOWN; - } - - calibrateImage(prow, nrow, offset, width, regr[bestFrame]); - - return (apt_channel_t)(channel + 1); - + float teleline[APT_MAX_HEIGHT] = {0.0}; + float wedge[16]; + linear_t regr[APT_MAX_HEIGHT / APT_FRAME_LEN + 1]; + int telestart, mtelestart = 0; + int channel = -1; + + // The minimum rows required to decode a full frame + if (nrow < APT_CALIBRATION_ROWS) { + error_noexit("Telemetry decoding error, not enough rows"); + return APT_CHANNEL_UNKNOWN; + } + + // Calculate average of a row of telemetry + for (int y = 0; y < nrow; y++) { + for (int x = 3; x < 43; x++) teleline[y] += prow[y][x + offset + width]; + + teleline[y] /= 40.0; + } + + /* Wedge 7 is white and 8 is black, this will have the largest + * difference in brightness, this can be used to find the current + * position within the telemetry. + */ + float max = 0.0f; + for (int n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) { + float df; + + // (sum 4px below) - (sum 4px above) + df = (float)((teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + teleline[n - 1]) - + (teleline[n + 0] + teleline[n + 1] + teleline[n + 2] + teleline[n + 3])); + + // Find the maximum difference + if (df > max) { + mtelestart = n; + max = df; + } + } + + telestart = (mtelestart + 64) % APT_FRAME_LEN; + + // Make sure that theres at least one full frame in the image + if (nrow < telestart + APT_FRAME_LEN) { + error_noexit("Telemetry decoding error, not enough rows"); + return APT_CHANNEL_UNKNOWN; + } + + // Find the least noisy frame + float minNoise = -1; + int bestFrame = -1; + for (int n = telestart, k = 0; n < nrow - APT_FRAME_LEN; n += APT_FRAME_LEN, k++) { + int j; + + for (j = 0; j < 16; j++) { + int i; + + wedge[j] = 0.0; + for (i = 1; i < 7; i++) wedge[j] += teleline[n + j * 8 + i]; + wedge[j] /= 6; + } + + float noise = teleNoise(wedge); + if (noise < minNoise || minNoise == -1) { + minNoise = noise; + bestFrame = k; + + // Compute & apply regression on the wedges + regr[k] = compute_regression(wedge); + for (int j = 0; j < 16; j++) tele[j] = linear_calc(wedge[j], regr[k]); + + /* Compare the channel ID wedge to the reference + * wedges, the wedge with the closest match will + * be the channel ID + */ + float min = -1; + for (int j = 0; j < 6; j++) { + float df = (float)(tele[15] - tele[j]); + df *= df; + + if (df < min || min == -1) { + channel = j; + min = df; + } + } + + // Find the brightness of the minute marker, I don't really know what for + Cs = 0.0; + int i, j = n; + for (i = 0, j = n; j < n + APT_FRAME_LEN; j++) { + float csline = 0.0; + for (int l = 3; l < 43; l++) csline += prow[n][l + offset - APT_SPC_WIDTH]; + csline /= 40.0; + + if (csline > 50.0) { + Cs += csline; + i++; + } + } + Cs = linear_calc((Cs / i), regr[k]); + } + } + + if (bestFrame == -1) { + error_noexit("Something has gone very wrong, please file a bug report"); + return APT_CHANNEL_UNKNOWN; + } + + calibrateImage(prow, nrow, offset, width, regr[bestFrame]); + + return (apt_channel_t)(channel + 1); } extern float quick_select(float arr[], int n); // Biased median denoise, pretyt ugly #define TRIG_LEVEL 40 -void apt_denoise(float **prow, int nrow, int offset, int width){ - for(int y = 2; y < nrow-2; y++){ - for(int x = offset+1; x < offset+width-1; x++){ - if(prow[y][x+1] - prow[y][x] > TRIG_LEVEL || - prow[y][x-1] - prow[y][x] > TRIG_LEVEL || - prow[y+1][x] - prow[y][x] > TRIG_LEVEL || - prow[y-1][x] - prow[y][x] > TRIG_LEVEL){ - prow[y][x] = quick_select((float[]){ - prow[y+2][x-1], prow[y+2][x], prow[y+2][x+1], - prow[y+1][x-1], prow[y+1][x], prow[y+1][x+1], - prow[y-1][x-1], prow[y-1][x], prow[y-1][x+1], - prow[y-2][x-1], prow[y-2][x], prow[y-2][x+1] - }, 12); - } - } - } +void apt_denoise(float **prow, int nrow, int offset, int width) { + for (int y = 2; y < nrow - 2; y++) { + for (int x = offset + 1; x < offset + width - 1; x++) { + if (prow[y][x + 1] - prow[y][x] > TRIG_LEVEL || prow[y][x - 1] - prow[y][x] > TRIG_LEVEL || + prow[y + 1][x] - prow[y][x] > TRIG_LEVEL || prow[y - 1][x] - prow[y][x] > TRIG_LEVEL) { + prow[y][x] = quick_select((float[]){prow[y + 2][x - 1], prow[y + 2][x], prow[y + 2][x + 1], prow[y + 1][x - 1], + prow[y + 1][x], prow[y + 1][x + 1], prow[y - 1][x - 1], prow[y - 1][x], + prow[y - 1][x + 1], prow[y - 2][x - 1], prow[y - 2][x], prow[y - 2][x + 1]}, + 12); + } + } + } } #undef TRIG_LEVEL // Flips a channel, for northbound passes -void apt_flipImage(apt_image_t *img, int width, int offset){ - for(int y = 1; y < img->nrow; y++){ - for(int x = 1; x < ceil(width / 2.0); x++){ - // Flip top-left & bottom-right - float buffer = img->prow[img->nrow - y][offset + x]; - img->prow[img->nrow - y][offset + x] = img->prow[y][offset + (width - x)]; - img->prow[y][offset + (width - x)] = buffer; - } - } +void apt_flipImage(apt_image_t *img, int width, int offset) { + for (int y = 1; y < img->nrow; y++) { + for (int x = 1; x < ceil(width / 2.0); x++) { + // Flip top-left & bottom-right + float buffer = img->prow[img->nrow - y][offset + x]; + img->prow[img->nrow - y][offset + x] = img->prow[y][offset + (width - x)]; + img->prow[y][offset + (width - x)] = buffer; + } + } } // Calculate crop to reomve noise from the start and end of an image -int apt_cropNoise(apt_image_t *img){ - #define NOISE_THRESH 180.0 - - // Average value of minute marker - float spc_rows[APT_MAX_HEIGHT] = { 0.0 }; - int startCrop = 0; int endCrop = img->nrow; - for(int y = 0; y < img->nrow; y++) { - for(int x = 0; x < APT_SPC_WIDTH; x++) { - spc_rows[y] += img->prow[y][x + (APT_CHB_OFFSET - APT_SPC_WIDTH)]; - } - spc_rows[y] /= APT_SPC_WIDTH; - - // Skip minute markings - if(spc_rows[y] < 10) { - spc_rows[y] = spc_rows[y-1]; - } - } - - // 3 row average - for(int y = 0; y < img->nrow; y++){ - spc_rows[y] = (spc_rows[y+1] + spc_rows[y+2] + spc_rows[y+3])/3; - //img.prow[y][0] = spc_rows[y]; - } - - // Find ends - for(int y = 0; y < img->nrow-1; y++) { - if(spc_rows[y] > NOISE_THRESH){ - endCrop = y; - } - } - for(int y = img->nrow; y > 0; y--) { - if(spc_rows[y] > NOISE_THRESH) { - startCrop = y; - } - } - - //printf("Crop rows: %i -> %i\n", startCrop, endCrop); - - // Remove the noisy rows at start - for(int y = 0; y < img->nrow-startCrop; y++) { - memmove(img->prow[y], img->prow[y+startCrop], sizeof(float)*APT_PROW_WIDTH); - } - - // Ignore the noisy rows at the end - img->nrow = (endCrop - startCrop); - - return startCrop; +int apt_cropNoise(apt_image_t *img) { +#define NOISE_THRESH 180.0 + + // Average value of minute marker + float spc_rows[APT_MAX_HEIGHT] = {0.0}; + int startCrop = 0; + int endCrop = img->nrow; + for (int y = 0; y < img->nrow; y++) { + for (int x = 0; x < APT_SPC_WIDTH; x++) { + spc_rows[y] += img->prow[y][x + (APT_CHB_OFFSET - APT_SPC_WIDTH)]; + } + spc_rows[y] /= APT_SPC_WIDTH; + + // Skip minute markings + if (spc_rows[y] < 10) { + spc_rows[y] = spc_rows[y - 1]; + } + } + + // 3 row average + for (int y = 0; y < img->nrow; y++) { + spc_rows[y] = (spc_rows[y + 1] + spc_rows[y + 2] + spc_rows[y + 3]) / 3; + // img.prow[y][0] = spc_rows[y]; + } + + // Find ends + for (int y = 0; y < img->nrow - 1; y++) { + if (spc_rows[y] > NOISE_THRESH) { + endCrop = y; + } + } + for (int y = img->nrow; y > 0; y--) { + if (spc_rows[y] > NOISE_THRESH) { + startCrop = y; + } + } + + // printf("Crop rows: %i -> %i\n", startCrop, endCrop); + + // Remove the noisy rows at start + for (int y = 0; y < img->nrow - startCrop; y++) { + memmove(img->prow[y], img->prow[y + startCrop], sizeof(float) * APT_PROW_WIDTH); + } + + // Ignore the noisy rows at the end + img->nrow = (endCrop - startCrop); + + return startCrop; } // --- Visible and Temperature Calibration --- // #include "calibration.h" typedef struct { - float Nbb; - float Cs; - float Cb; - int ch; + float Nbb; + float Cs; + float Cb; + int ch; } tempparam_t; // IR channel temperature compensation tempparam_t tempcomp(float t[16], int ch, int satnum) { - tempparam_t tpr; - tpr.ch = ch - 4; + tempparam_t tpr; + tpr.ch = ch - 4; - const calibration_t calibration = get_calibration(satnum); - const float Vc = calibration.rad[tpr.ch].vc; - const float A = calibration.rad[tpr.ch].A; - const float B = calibration.rad[tpr.ch].B; + const calibration_t calibration = get_calibration(satnum); + const float Vc = calibration.rad[tpr.ch].vc; + const float A = calibration.rad[tpr.ch].A; + const float B = calibration.rad[tpr.ch].B; - // Compute PRT temperature - float T[4]; - for (size_t n = 0; n < 4; n++) { - T[n] = quadratic_calc(t[n+9] * 4.0, calibration.prt[n]); - } + // Compute PRT temperature + float T[4]; + for (size_t n = 0; n < 4; n++) { + T[n] = quadratic_calc(t[n + 9] * 4.0, calibration.prt[n]); + } - float Tbb = (T[0]+T[1]+T[2]+T[3]) / 4.0; // Blackbody temperature - float Tbbstar = A + Tbb*B; // Effective blackbody temperature + float Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0; // Blackbody temperature + float Tbbstar = A + Tbb * B; // Effective blackbody temperature - tpr.Nbb = C1 * pow(Vc, 3) / (expf(C2 * Vc / Tbbstar) - 1.0f); // Blackbody radiance + tpr.Nbb = C1 * pow(Vc, 3) / (expf(C2 * Vc / Tbbstar) - 1.0f); // Blackbody radiance - tpr.Cs = 246.4 * 4.0; // FIXME - tpr.Cb = t[14] * 4.0; - return tpr; + tpr.Cs = 246.4 * 4.0; // FIXME + tpr.Cb = t[14] * 4.0; + return tpr; } // IR channel temperature calibration static float tempcal(float Ce, int satnum, tempparam_t tpr) { - const calibration_t calibration = get_calibration(satnum); - const float Ns = calibration.cor[tpr.ch].Ns; - const float Vc = calibration.rad[tpr.ch].vc; - const float A = calibration.rad[tpr.ch].A; - const float B = calibration.rad[tpr.ch].B; - - float Nl = Ns + (tpr.Nbb - Ns) * (tpr.Cs - Ce * 4.0) / (tpr.Cs - tpr.Cb); // Linear radiance estimate - float Nc = quadratic_calc(Nl, calibration.cor[tpr.ch].quadratic); // Non-linear correction - float Ne = Nl + Nc; // Corrected radiance - - float Testar = C2 * Vc / logf(C1 * powf(Vc, 3) / Ne + 1.0); // Equivlent black body temperature - float Te = (Testar - A) / B; // Temperature (kelvin) - - // Convert to celsius - Te -= 273.15; - // Rescale to 0-255 for -100°C to +60°C - return (Te + 100.0) / 160.0 * 255.0; + const calibration_t calibration = get_calibration(satnum); + const float Ns = calibration.cor[tpr.ch].Ns; + const float Vc = calibration.rad[tpr.ch].vc; + const float A = calibration.rad[tpr.ch].A; + const float B = calibration.rad[tpr.ch].B; + + float Nl = Ns + (tpr.Nbb - Ns) * (tpr.Cs - Ce * 4.0) / (tpr.Cs - tpr.Cb); // Linear radiance estimate + float Nc = quadratic_calc(Nl, calibration.cor[tpr.ch].quadratic); // Non-linear correction + float Ne = Nl + Nc; // Corrected radiance + + float Testar = C2 * Vc / logf(C1 * powf(Vc, 3) / Ne + 1.0); // Equivlent black body temperature + float Te = (Testar - A) / B; // Temperature (kelvin) + + // Convert to celsius + Te -= 273.15; + // Rescale to 0-255 for -100°C to +60°C + return (Te + 100.0) / 160.0 * 255.0; } // Temperature calibration wrapper void apt_calibrate_thermal(int satnum, apt_image_t *img, int offset, int width) { - tempparam_t temp = tempcomp(tele, img->chB, satnum); + tempparam_t temp = tempcomp(tele, img->chB, satnum); - for (int y = 0; y < img->nrow; y++) { - for (int x = 0; x < width; x++) { - img->prow[y][x + offset] = (float)tempcal(img->prow[y][x + offset], satnum, temp); - } - } + for (int y = 0; y < img->nrow; y++) { + for (int x = 0; x < width; x++) { + img->prow[y][x + offset] = (float)tempcal(img->prow[y][x + offset], satnum, temp); + } + } } float calibrate_pixel(float value, int channel, calibration_t cal) { - if (value > cal.visible[channel].cutoff) { - return linear_calc(value*4.0f, cal.visible[channel].high) * 255.0f/100.0f; - } else { - return linear_calc(value*4.0f, cal.visible[channel].low) * 255.0f/100.0f; - } + if (value > cal.visible[channel].cutoff) { + return linear_calc(value * 4.0f, cal.visible[channel].high) * 255.0f / 100.0f; + } else { + return linear_calc(value * 4.0f, cal.visible[channel].low) * 255.0f / 100.0f; + } } void apt_calibrate_visible(int satnum, apt_image_t *img, int offset, int width) { - const calibration_t calibration = get_calibration(satnum); - int channel = img->chA-1; - - for (int y = 0; y < img->nrow; y++) { - for (int x = 0; x < width; x++) { - img->prow[y][x + offset] = clamp(calibrate_pixel(img->prow[y][x + offset], channel, calibration), 255.0f, 0.0f); - } - } + const calibration_t calibration = get_calibration(satnum); + int channel = img->chA - 1; + + for (int y = 0; y < img->nrow; y++) { + for (int x = 0; x < width; x++) { + img->prow[y][x + offset] = clamp(calibrate_pixel(img->prow[y][x + offset], channel, calibration), 255.0f, 0.0f); + } + } } diff --git a/src/main.c b/src/main.c index 51d91a4..1ebc488 100644 --- a/src/main.c +++ b/src/main.c @@ -1,4 +1,4 @@ -/* +/* * This file is part of Aptdec. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 * @@ -17,26 +17,25 @@ * */ -#include #include +#include #include #ifndef _MSC_VER #include #else #include #endif +#include #include #include -#include #include -#include "argparse/argparse.h" -#include "common.h" #include "apt.h" - -#include "pngio.h" -#include "image.h" +#include "argparse/argparse.h" #include "color.h" +#include "common.h" +#include "image.h" +#include "pngio.h" #include "util.h" // Audio file @@ -51,279 +50,267 @@ static int processAudio(char *filename, options_t *opts); #ifdef _MSC_VER // Functions not supported by MSVC -static char *dirname(char *path) -{ - static char dir[MAX_PATH]; - _splitpath(path, NULL, dir, NULL, NULL); - return dir; +static char *dirname(char *path) { + static char dir[MAX_PATH]; + _splitpath(path, NULL, dir, NULL, NULL); + return dir; } -static char *basename(char *path) -{ - static char base[MAX_PATH]; - _splitpath(path, NULL, NULL, base, NULL); - return base; +static char *basename(char *path) { + static char base[MAX_PATH]; + _splitpath(path, NULL, NULL, base, NULL); + return base; } #endif int main(int argc, const char **argv) { - options_t opts = { - .type = "r", - .effects = "", - .satnum = 19, - .path = ".", - .realtime = 0, - .filename = "", - .palette = "", - .gamma = 1.0 - }; - - static const char *const usages[] = { - "aptdec [options] [[--] sources]", - "aptdec [sources]", - NULL, - }; - - struct argparse_option options[] = { - OPT_HELP(), - OPT_GROUP("Image options"), - OPT_STRING('i', "image", &opts.type, "set output image type (see the README for a list)", NULL, 0, 0), - OPT_STRING('e', "effect", &opts.effects, "add an effect (see the README for a list)", NULL, 0, 0), - OPT_FLOAT('g', "gamma", &opts.gamma, "gamma adjustment (1.0 = off)", NULL, 0, 0), - - OPT_GROUP("Satellite options"), - OPT_INTEGER('s', "satellite", &opts.satnum, "satellite ID, must be between 15 and 19", NULL, 0, 0), - - OPT_GROUP("Paths"), - OPT_STRING('p', "palette", &opts.palette, "path to a palette", NULL, 0, 0), - OPT_STRING('o', "filename", &opts.filename, "filename of the output image", NULL, 0, 0), - OPT_STRING('d', "output", &opts.path, "output directory (must exist first)", NULL, 0, 0), - - OPT_GROUP("Misc"), - OPT_BOOLEAN('r', "realtime", &opts.realtime, "decode in realtime", NULL, 0, 0), - OPT_END(), - }; - - struct argparse argparse; - argparse_init(&argparse, options, usages, 0); - argparse_describe(&argparse, "\nA lightweight FOSS NOAA APT satellite imagery decoder.", "\nSee `README.md` for a full description of command line arguments and `LICENSE` for licensing conditions."); - argc = argparse_parse(&argparse, argc, argv); - - if(argc == 0){ - argparse_usage(&argparse); - } - - // Actually decode the files - for (int i = 0; i < argc; i++) { - char *filename = strdup(argv[i]); - processAudio(filename, &opts); - } - - return 0; + options_t opts = { + .type = "r", .effects = "", .satnum = 19, .path = ".", .realtime = 0, .filename = "", .palette = "", .gamma = 1.0}; + + static const char *const usages[] = { + "aptdec [options] [[--] sources]", + "aptdec [sources]", + NULL, + }; + + struct argparse_option options[] = { + OPT_HELP(), + OPT_GROUP("Image options"), + OPT_STRING('i', "image", &opts.type, "set output image type (see the README for a list)", NULL, 0, 0), + OPT_STRING('e', "effect", &opts.effects, "add an effect (see the README for a list)", NULL, 0, 0), + OPT_FLOAT('g', "gamma", &opts.gamma, "gamma adjustment (1.0 = off)", NULL, 0, 0), + + OPT_GROUP("Satellite options"), + OPT_INTEGER('s', "satellite", &opts.satnum, "satellite ID, must be between 15 and 19", NULL, 0, 0), + + OPT_GROUP("Paths"), + OPT_STRING('p', "palette", &opts.palette, "path to a palette", NULL, 0, 0), + OPT_STRING('o', "filename", &opts.filename, "filename of the output image", NULL, 0, 0), + OPT_STRING('d', "output", &opts.path, "output directory (must exist first)", NULL, 0, 0), + + OPT_GROUP("Misc"), + OPT_BOOLEAN('r', "realtime", &opts.realtime, "decode in realtime", NULL, 0, 0), + OPT_END(), + }; + + struct argparse argparse; + argparse_init(&argparse, options, usages, 0); + argparse_describe( + &argparse, "\nA lightweight FOSS NOAA APT satellite imagery decoder.", + "\nSee `README.md` for a full description of command line arguments and `LICENSE` for licensing conditions."); + argc = argparse_parse(&argparse, argc, argv); + + if (argc == 0) { + argparse_usage(&argparse); + } + + // Actually decode the files + for (int i = 0; i < argc; i++) { + char *filename = strdup(argv[i]); + processAudio(filename, &opts); + } + + return 0; } -static int processAudio(char *filename, options_t *opts){ - // Image info struct - apt_image_t img; - - // Mapping between wedge value and channel ID - static struct { - char *id[7]; - char *name[7]; - } ch = { - { "?", "1", "2", "3A", "4", "5", "3B" }, - { "unknown", "visble", "near-infrared", "near-infrared", "thermal-infrared", "thermal-infrared", "mid-infrared" } - }; - - // Buffer for image channel - char desc[60]; - - // Parse file path - char path[256], extension[32]; - strcpy(path, filename); - strcpy(path, dirname(path)); - sscanf(basename(filename), "%255[^.].%31s", img.name, extension); - - if(opts->realtime){ - // Set output filename to current time when in realtime mode - time_t t; - time(&t); - strncpy(img.name, ctime(&t), 24); - - // Init a row writer - initWriter(opts, &img, APT_IMG_WIDTH, APT_MAX_HEIGHT, "Unprocessed realtime image", "r"); - } - - if(strcmp(extension, "png") == 0){ - // Read PNG into image buffer - printf("Reading %s\n", filename); - if(readRawImage(filename, img.prow, &img.nrow) == 0){ - exit(EPERM); - } - }else{ - // Attempt to open the audio file - if (initsnd(filename) == 0) - exit(EPERM); - - // Build image - // TODO: multithreading, would require some sort of input buffer - for (img.nrow = 0; img.nrow < APT_MAX_HEIGHT; img.nrow++) { - // Allocate memory for this row - img.prow[img.nrow] = (float *) malloc(sizeof(float) * APT_PROW_WIDTH); - - // Write into memory and break the loop when there are no more samples to read - if (apt_getpixelrow(img.prow[img.nrow], img.nrow, &img.zenith, (img.nrow == 0), getsamples, NULL) == 0) - break; - - if(opts->realtime) pushRow(img.prow[img.nrow], APT_IMG_WIDTH); - - fprintf(stderr, "Row: %d\r", img.nrow); - fflush(stderr); - } - - // Close stream - sf_close(audioFile); - } - - if(opts->realtime) closeWriter(); - - printf("Total rows: %d\n", img.nrow); - - // Calibrate - img.chA = apt_calibrate(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); - img.chB = apt_calibrate(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); - printf("Channel A: %s (%s)\n", ch.id[img.chA], ch.name[img.chA]); - printf("Channel B: %s (%s)\n", ch.id[img.chB], ch.name[img.chB]); - - // Crop noise from start and end of image - if(CONTAINS(opts->effects, Crop_Noise)){ - img.zenith -= apt_cropNoise(&img); - } - - // Denoise - if(CONTAINS(opts->effects, Denoise)){ - apt_denoise(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); - apt_denoise(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); - } - - // Flip, for northbound passes - if(CONTAINS(opts->effects, Flip_Image)){ - apt_flipImage(&img, APT_CH_WIDTH, APT_CHA_OFFSET); - apt_flipImage(&img, APT_CH_WIDTH, APT_CHB_OFFSET); - } - - // Temperature - if (CONTAINS(opts->type, Temperature) && img.chB >= 4) { - // Create another buffer as to not modify the orignal - apt_image_t tmpimg = img; - for(int i = 0; i < img.nrow; i++){ - tmpimg.prow[i] = (float *) malloc(sizeof(float) * APT_PROW_WIDTH); - memcpy(tmpimg.prow[i], img.prow[i], sizeof(float) * APT_PROW_WIDTH); - } - - // Perform temperature calibration - apt_calibrate_thermal(opts->satnum, &tmpimg, APT_CHB_OFFSET, APT_CH_WIDTH); - ImageOut(opts, &tmpimg, APT_CHB_OFFSET, APT_CH_WIDTH, "Temperature", Temperature, (char *)apt_TempPalette); - } - - // Visible - if (CONTAINS(opts->type, Visible) && img.chA <= 2) { - // Create another buffer as to not modify the orignal - apt_image_t tmpimg = img; - for(int i = 0; i < img.nrow; i++){ - tmpimg.prow[i] = (float *) malloc(sizeof(float) * APT_PROW_WIDTH); - memcpy(tmpimg.prow[i], img.prow[i], sizeof(float) * APT_PROW_WIDTH); - } - - // Perform visible calibration - apt_calibrate_visible(opts->satnum, &tmpimg, APT_CHA_OFFSET, APT_CH_WIDTH); - ImageOut(opts, &tmpimg, APT_CHA_OFFSET, APT_CH_WIDTH, "Visible", Visible, NULL); - } - - // Linear equalise - if(CONTAINS(opts->effects, Linear_Equalise)){ - apt_linearEnhance(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); - apt_linearEnhance(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); - } - - // Histogram equalise - if(CONTAINS(opts->effects, Histogram_Equalise)){ - apt_histogramEqualise(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); - apt_histogramEqualise(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); - } - - // Raw image - if (CONTAINS(opts->type, Raw_Image)) { - sprintf(desc, "%s (%s) & %s (%s)", ch.id[img.chA], ch.name[img.chA], ch.id[img.chB], ch.name[img.chB]); - ImageOut(opts, &img, 0, APT_IMG_WIDTH, desc, Raw_Image, NULL); - } - - // Palette image - if (CONTAINS(opts->type, Palleted)) { - img.palette = opts->palette; - strcpy(desc, "Palette composite"); - ImageOut(opts, &img, APT_CHA_OFFSET, APT_CH_WIDTH, desc, Palleted, NULL); - } - - // Channel A - if (CONTAINS(opts->type, Channel_A)) { - sprintf(desc, "%s (%s)", ch.id[img.chA], ch.name[img.chA]); - ImageOut(opts, &img, APT_CHA_OFFSET, APT_CH_WIDTH, desc, Channel_A, NULL); - } - - // Channel B - if (CONTAINS(opts->type, Channel_B)) { - sprintf(desc, "%s (%s)", ch.id[img.chB], ch.name[img.chB]); - ImageOut(opts, &img, APT_CHB_OFFSET, APT_CH_WIDTH, desc, Channel_B, NULL); - } - - return 1; +static int processAudio(char *filename, options_t *opts) { + // Image info struct + apt_image_t img; + + // Mapping between wedge value and channel ID + static struct { + char *id[7]; + char *name[7]; + } ch = {{"?", "1", "2", "3A", "4", "5", "3B"}, + {"unknown", "visble", "near-infrared", "near-infrared", "thermal-infrared", "thermal-infrared", "mid-infrared"}}; + + // Buffer for image channel + char desc[60]; + + // Parse file path + char path[256], extension[32]; + strcpy(path, filename); + strcpy(path, dirname(path)); + sscanf(basename(filename), "%255[^.].%31s", img.name, extension); + + if (opts->realtime) { + // Set output filename to current time when in realtime mode + time_t t; + time(&t); + strncpy(img.name, ctime(&t), 24); + + // Init a row writer + initWriter(opts, &img, APT_IMG_WIDTH, APT_MAX_HEIGHT, "Unprocessed realtime image", "r"); + } + + if (strcmp(extension, "png") == 0) { + // Read PNG into image buffer + printf("Reading %s\n", filename); + if (readRawImage(filename, img.prow, &img.nrow) == 0) { + exit(EPERM); + } + } else { + // Attempt to open the audio file + if (initsnd(filename) == 0) exit(EPERM); + + // Build image + // TODO: multithreading, would require some sort of input buffer + for (img.nrow = 0; img.nrow < APT_MAX_HEIGHT; img.nrow++) { + // Allocate memory for this row + img.prow[img.nrow] = (float *)malloc(sizeof(float) * APT_PROW_WIDTH); + + // Write into memory and break the loop when there are no more samples to read + if (apt_getpixelrow(img.prow[img.nrow], img.nrow, &img.zenith, (img.nrow == 0), getsamples, NULL) == 0) break; + + if (opts->realtime) pushRow(img.prow[img.nrow], APT_IMG_WIDTH); + + fprintf(stderr, "Row: %d\r", img.nrow); + fflush(stderr); + } + + // Close stream + sf_close(audioFile); + } + + if (opts->realtime) closeWriter(); + + printf("Total rows: %d\n", img.nrow); + + // Calibrate + img.chA = apt_calibrate(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); + img.chB = apt_calibrate(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); + printf("Channel A: %s (%s)\n", ch.id[img.chA], ch.name[img.chA]); + printf("Channel B: %s (%s)\n", ch.id[img.chB], ch.name[img.chB]); + + // Crop noise from start and end of image + if (CONTAINS(opts->effects, Crop_Noise)) { + img.zenith -= apt_cropNoise(&img); + } + + // Denoise + if (CONTAINS(opts->effects, Denoise)) { + apt_denoise(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); + apt_denoise(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); + } + + // Flip, for northbound passes + if (CONTAINS(opts->effects, Flip_Image)) { + apt_flipImage(&img, APT_CH_WIDTH, APT_CHA_OFFSET); + apt_flipImage(&img, APT_CH_WIDTH, APT_CHB_OFFSET); + } + + // Temperature + if (CONTAINS(opts->type, Temperature) && img.chB >= 4) { + // Create another buffer as to not modify the orignal + apt_image_t tmpimg = img; + for (int i = 0; i < img.nrow; i++) { + tmpimg.prow[i] = (float *)malloc(sizeof(float) * APT_PROW_WIDTH); + memcpy(tmpimg.prow[i], img.prow[i], sizeof(float) * APT_PROW_WIDTH); + } + + // Perform temperature calibration + apt_calibrate_thermal(opts->satnum, &tmpimg, APT_CHB_OFFSET, APT_CH_WIDTH); + ImageOut(opts, &tmpimg, APT_CHB_OFFSET, APT_CH_WIDTH, "Temperature", Temperature, (char *)apt_TempPalette); + } + + // Visible + if (CONTAINS(opts->type, Visible) && img.chA <= 2) { + // Create another buffer as to not modify the orignal + apt_image_t tmpimg = img; + for (int i = 0; i < img.nrow; i++) { + tmpimg.prow[i] = (float *)malloc(sizeof(float) * APT_PROW_WIDTH); + memcpy(tmpimg.prow[i], img.prow[i], sizeof(float) * APT_PROW_WIDTH); + } + + // Perform visible calibration + apt_calibrate_visible(opts->satnum, &tmpimg, APT_CHA_OFFSET, APT_CH_WIDTH); + ImageOut(opts, &tmpimg, APT_CHA_OFFSET, APT_CH_WIDTH, "Visible", Visible, NULL); + } + + // Linear equalise + if (CONTAINS(opts->effects, Linear_Equalise)) { + apt_linearEnhance(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); + apt_linearEnhance(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); + } + + // Histogram equalise + if (CONTAINS(opts->effects, Histogram_Equalise)) { + apt_histogramEqualise(img.prow, img.nrow, APT_CHA_OFFSET, APT_CH_WIDTH); + apt_histogramEqualise(img.prow, img.nrow, APT_CHB_OFFSET, APT_CH_WIDTH); + } + + // Raw image + if (CONTAINS(opts->type, Raw_Image)) { + sprintf(desc, "%s (%s) & %s (%s)", ch.id[img.chA], ch.name[img.chA], ch.id[img.chB], ch.name[img.chB]); + ImageOut(opts, &img, 0, APT_IMG_WIDTH, desc, Raw_Image, NULL); + } + + // Palette image + if (CONTAINS(opts->type, Palleted)) { + img.palette = opts->palette; + strcpy(desc, "Palette composite"); + ImageOut(opts, &img, APT_CHA_OFFSET, APT_CH_WIDTH, desc, Palleted, NULL); + } + + // Channel A + if (CONTAINS(opts->type, Channel_A)) { + sprintf(desc, "%s (%s)", ch.id[img.chA], ch.name[img.chA]); + ImageOut(opts, &img, APT_CHA_OFFSET, APT_CH_WIDTH, desc, Channel_A, NULL); + } + + // Channel B + if (CONTAINS(opts->type, Channel_B)) { + sprintf(desc, "%s (%s)", ch.id[img.chB], ch.name[img.chB]); + ImageOut(opts, &img, APT_CHB_OFFSET, APT_CH_WIDTH, desc, Channel_B, NULL); + } + + return 1; } float *samplebuf; static int initsnd(char *filename) { - SF_INFO infwav; - int res; - - // Open audio file - infwav.format = 0; - audioFile = sf_open(filename, SFM_READ, &infwav); - if (audioFile == NULL) { - error_noexit("Could not file"); - return 0; - } - - res = apt_init(infwav.samplerate); - printf("Input file: %s\n", filename); - if(res < 0) { - error_noexit("Input sample rate too low"); - return 0; - }else if(res > 0) { - error_noexit("Input sample rate too high"); - return 0; - } - printf("Input sample rate: %d\n", infwav.samplerate); - - channels = infwav.channels; - samplebuf = (float *)malloc(sizeof(float) * 32768 * channels); - - return 1; + SF_INFO infwav; + int res; + + // Open audio file + infwav.format = 0; + audioFile = sf_open(filename, SFM_READ, &infwav); + if (audioFile == NULL) { + error_noexit("Could not file"); + return 0; + } + + res = apt_init(infwav.samplerate); + printf("Input file: %s\n", filename); + if (res < 0) { + error_noexit("Input sample rate too low"); + return 0; + } else if (res > 0) { + error_noexit("Input sample rate too high"); + return 0; + } + printf("Input sample rate: %d\n", infwav.samplerate); + + channels = infwav.channels; + samplebuf = (float *)malloc(sizeof(float) * 32768 * channels); + + return 1; } // Read samples from the audio file int getsamples(void *context, float *samples, int nb) { - (void) context; - if (channels == 1){ - return (int)sf_read_float(audioFile, samples, nb); - } else if (channels == 2) { - // Stereo channels are interleaved - int samplesRead = (int)sf_read_float(audioFile, samplebuf, nb * channels); - for(int i = 0; i < nb; i++) { - samples[i] = samplebuf[i * channels]; - } - return samplesRead / channels; - } else { - printf("Only mono and stereo input files are supported\n"); - exit(1); - } + (void)context; + if (channels == 1) { + return (int)sf_read_float(audioFile, samples, nb); + } else if (channels == 2) { + // Stereo channels are interleaved + int samplesRead = (int)sf_read_float(audioFile, samplebuf, nb * channels); + for (int i = 0; i < nb; i++) { + samples[i] = samplebuf[i * channels]; + } + return samplesRead / channels; + } else { + printf("Only mono and stereo input files are supported\n"); + exit(1); + } } diff --git a/src/pngio.c b/src/pngio.c index 2e4f9bd..a061b0a 100644 --- a/src/pngio.c +++ b/src/pngio.c @@ -1,4 +1,4 @@ -/* +/* * This file is part of Aptdec. * Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2022 * @@ -17,337 +17,333 @@ * */ +#include "pngio.h" + +#include #include -#include +#include #include +#include #include -#include -#include -#include "util.h" -#include "pngio.h" +#include "util.h" int readRawImage(char *filename, float **prow, int *nrow) { - FILE *fp = fopen(filename, "rb"); + FILE *fp = fopen(filename, "rb"); printf("%s", filename); - if(!fp) { - error_noexit("Cannot open image"); - return 0; - } - - // Create reader - png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if(!png) { - fclose(fp); - return 0; - } - png_infop info = png_create_info_struct(png); - if(!info) { - fclose(fp); - return 0; - } - png_init_io(png, fp); - - // Read info from header - png_read_info(png, info); - int width = png_get_image_width(png, info); - int height = png_get_image_height(png, info); - png_byte color_type = png_get_color_type(png, info); - png_byte bit_depth = png_get_bit_depth(png, info); - - // Check the image - if(width != APT_IMG_WIDTH){ - error_noexit("Raw image must be 2080px wide"); - return 0; - }else if(bit_depth != 8){ - error_noexit("Raw image must have 8 bit color"); - return 0; - }else if(color_type != PNG_COLOR_TYPE_GRAY){ - error_noexit("Raw image must be grayscale"); - return 0; - } - - // Create row buffers - png_bytep *PNGrows = NULL; - PNGrows = (png_bytep *) malloc(sizeof(png_bytep) * height); - for(int y = 0; y < height; y++) PNGrows[y] = (png_byte *) - malloc(png_get_rowbytes(png, info)); - - // Read image - png_read_image(png, PNGrows); - - // Tidy up - fclose(fp); - png_destroy_read_struct(&png, &info, NULL); - - // Put into prow - *nrow = height; - for(int y = 0; y < height; y++) { - prow[y] = (float *) malloc(sizeof(float) * width); - - for(int x = 0; x < width; x++) - prow[y][x] = (float)PNGrows[y][x]; - } - - return 1; + if (!fp) { + error_noexit("Cannot open image"); + return 0; + } + + // Create reader + png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!png) { + fclose(fp); + return 0; + } + png_infop info = png_create_info_struct(png); + if (!info) { + fclose(fp); + return 0; + } + png_init_io(png, fp); + + // Read info from header + png_read_info(png, info); + int width = png_get_image_width(png, info); + int height = png_get_image_height(png, info); + png_byte color_type = png_get_color_type(png, info); + png_byte bit_depth = png_get_bit_depth(png, info); + + // Check the image + if (width != APT_IMG_WIDTH) { + error_noexit("Raw image must be 2080px wide"); + return 0; + } else if (bit_depth != 8) { + error_noexit("Raw image must have 8 bit color"); + return 0; + } else if (color_type != PNG_COLOR_TYPE_GRAY) { + error_noexit("Raw image must be grayscale"); + return 0; + } + + // Create row buffers + png_bytep *PNGrows = NULL; + PNGrows = (png_bytep *)malloc(sizeof(png_bytep) * height); + for (int y = 0; y < height; y++) PNGrows[y] = (png_byte *)malloc(png_get_rowbytes(png, info)); + + // Read image + png_read_image(png, PNGrows); + + // Tidy up + fclose(fp); + png_destroy_read_struct(&png, &info, NULL); + + // Put into prow + *nrow = height; + for (int y = 0; y < height; y++) { + prow[y] = (float *)malloc(sizeof(float) * width); + + for (int x = 0; x < width; x++) prow[y][x] = (float)PNGrows[y][x]; + } + + return 1; } int readPalette(char *filename, apt_rgb_t **pixels) { - FILE *fp = fopen(filename, "rb"); - if(!fp) { - error_noexit("Cannot open palette"); - return 0; - } - - // Create reader - png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if(!png) { - fclose(fp); - return 0; - } - png_infop info = png_create_info_struct(png); - if(!info) { - fclose(fp); - return 0; - } - png_init_io(png, fp); - - // Read info from header - png_read_info(png, info); - int width = png_get_image_width(png, info); - int height = png_get_image_height(png, info); - png_byte color_type = png_get_color_type(png, info); - png_byte bit_depth = png_get_bit_depth(png, info); - - // Check the image - if(width != 256 && height != 256){ - error_noexit("Palette must be 256x256"); - return 0; - }else if(bit_depth != 8){ - error_noexit("Palette must be 8 bit color"); - return 0; - }else if(color_type != PNG_COLOR_TYPE_RGB){ - error_noexit("Palette must be RGB"); - return 0; - } - - // Create row buffers - png_bytep *PNGrows = NULL; - PNGrows = (png_bytep *) malloc(sizeof(png_bytep) * height); - for(int y = 0; y < height; y++) - PNGrows[y] = (png_byte *) malloc(png_get_rowbytes(png, info)); - - // Read image - png_read_image(png, PNGrows); - - // Tidy up - fclose(fp); - png_destroy_read_struct(&png, &info, NULL); - - // Put into crow - for(int y = 0; y < height; y++) { - pixels[y] = (apt_rgb_t *) malloc(sizeof(apt_rgb_t) * width); - - for(int x = 0; x < width; x++) - pixels[y][x] = (apt_rgb_t){ - PNGrows[y][x*3], - PNGrows[y][x*3 + 1], - PNGrows[y][x*3 + 2] - }; - } - - return 1; + FILE *fp = fopen(filename, "rb"); + if (!fp) { + error_noexit("Cannot open palette"); + return 0; + } + + // Create reader + png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!png) { + fclose(fp); + return 0; + } + png_infop info = png_create_info_struct(png); + if (!info) { + fclose(fp); + return 0; + } + png_init_io(png, fp); + + // Read info from header + png_read_info(png, info); + int width = png_get_image_width(png, info); + int height = png_get_image_height(png, info); + png_byte color_type = png_get_color_type(png, info); + png_byte bit_depth = png_get_bit_depth(png, info); + + // Check the image + if (width != 256 && height != 256) { + error_noexit("Palette must be 256x256"); + return 0; + } else if (bit_depth != 8) { + error_noexit("Palette must be 8 bit color"); + return 0; + } else if (color_type != PNG_COLOR_TYPE_RGB) { + error_noexit("Palette must be RGB"); + return 0; + } + + // Create row buffers + png_bytep *PNGrows = NULL; + PNGrows = (png_bytep *)malloc(sizeof(png_bytep) * height); + for (int y = 0; y < height; y++) PNGrows[y] = (png_byte *)malloc(png_get_rowbytes(png, info)); + + // Read image + png_read_image(png, PNGrows); + + // Tidy up + fclose(fp); + png_destroy_read_struct(&png, &info, NULL); + + // Put into crow + for (int y = 0; y < height; y++) { + pixels[y] = (apt_rgb_t *)malloc(sizeof(apt_rgb_t) * width); + + for (int x = 0; x < width; x++) + pixels[y][x] = (apt_rgb_t){PNGrows[y][x * 3], PNGrows[y][x * 3 + 1], PNGrows[y][x * 3 + 2]}; + } + + return 1; } -void prow2crow(float **prow, int nrow, char *palette, apt_rgb_t **crow){ - for(int y = 0; y < nrow; y++){ - crow[y] = (apt_rgb_t *) malloc(sizeof(apt_rgb_t) * APT_IMG_WIDTH); - - for(int x = 0; x < APT_IMG_WIDTH; x++){ - if(palette == NULL) - crow[y][x].r = crow[y][x].g = crow[y][x].b = prow[y][x]; - else - crow[y][x] = apt_applyPalette(palette, prow[y][x]); - } - } +void prow2crow(float **prow, int nrow, char *palette, apt_rgb_t **crow) { + for (int y = 0; y < nrow; y++) { + crow[y] = (apt_rgb_t *)malloc(sizeof(apt_rgb_t) * APT_IMG_WIDTH); + + for (int x = 0; x < APT_IMG_WIDTH; x++) { + if (palette == NULL) + crow[y][x].r = crow[y][x].g = crow[y][x].b = prow[y][x]; + else + crow[y][x] = apt_applyPalette(palette, prow[y][x]); + } + } } -int applyUserPalette(float **prow, int nrow, char *filename, apt_rgb_t **crow){ - apt_rgb_t *pal_row[256]; - if(!readPalette(filename, pal_row)){ - error_noexit("Could not read palette"); - return 0; - } - - for(int y = 0; y < nrow; y++){ - for(int x = 0; x < APT_CH_WIDTH; x++){ - int cha = CLIP(prow[y][x + APT_CHA_OFFSET], 0, 255); - int chb = CLIP(prow[y][x + APT_CHB_OFFSET], 0, 255); - crow[y][x + APT_CHA_OFFSET] = pal_row[chb][cha]; - } - } - - return 1; +int applyUserPalette(float **prow, int nrow, char *filename, apt_rgb_t **crow) { + apt_rgb_t *pal_row[256]; + if (!readPalette(filename, pal_row)) { + error_noexit("Could not read palette"); + return 0; + } + + for (int y = 0; y < nrow; y++) { + for (int x = 0; x < APT_CH_WIDTH; x++) { + int cha = CLIP(prow[y][x + APT_CHA_OFFSET], 0, 255); + int chb = CLIP(prow[y][x + APT_CHB_OFFSET], 0, 255); + crow[y][x + APT_CHA_OFFSET] = pal_row[chb][cha]; + } + } + + return 1; } -int ImageOut(options_t *opts, apt_image_t *img, int offset, int width, char *desc, char chid, char *palette){ - char outName[512]; - if(opts->filename == NULL || opts->filename[0] == '\0'){ - sprintf(outName, "%s/%s-%c.png", opts->path, img->name, chid); - }else{ - sprintf(outName, "%s/%s", opts->path, opts->filename); - } - - png_text meta[] = { - {PNG_TEXT_COMPRESSION_NONE, "Software", VERSION, sizeof(VERSION)}, - {PNG_TEXT_COMPRESSION_NONE, "Channel", desc, sizeof(desc)}, - {PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20} - }; - - // Parse image type - int greyscale = 1; - switch (chid){ - case Palleted: - greyscale = 0; - break; - case Temperature: - greyscale = 0; - break; - case Raw_Image: break; - case Channel_A: break; - case Channel_B: break; - } - - // Parse effects - int crop_telemetry = 0; - for(unsigned long int i = 0; i < strlen(opts->effects); i++){ - switch (opts->effects[i]) { - case Crop_Telemetry: - if (width == 2080) { - width -= APT_TOTAL_TELE; - offset += APT_SYNC_WIDTH + APT_SPC_WIDTH; - crop_telemetry = 1; - } - break; - case Precipitation_Overlay: - greyscale = 0; - break; - case Flip_Image: break; - case Denoise: break; - case Histogram_Equalise: break; - case Linear_Equalise: break; - case Crop_Noise: break; - default: { - char text[100]; - sprintf(text, "Unrecognised effect, \"%c\"", opts->effects[i]); - warning(text); - break; - } - } - } - - FILE *pngfile; - - // Create writer - png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if (!png_ptr) { - png_destroy_write_struct(&png_ptr, (png_infopp) NULL); - error_noexit("Could not create a PNG writer"); - return 0; - } - png_infop info_ptr = png_create_info_struct(png_ptr); - if (!info_ptr) { - png_destroy_write_struct(&png_ptr, (png_infopp) NULL); - error_noexit("Could not create a PNG writer"); - return 0; - } - - if(greyscale){ - // Greyscale image - png_set_IHDR(png_ptr, info_ptr, width, img->nrow, - 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, - PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); - }else{ - // 8 bit RGB image - png_set_IHDR(png_ptr, info_ptr, width, img->nrow, - 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, - PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); - } - - png_set_text(png_ptr, info_ptr, meta, 3); - png_set_pHYs(png_ptr, info_ptr, 3636, 3636, PNG_RESOLUTION_METER); - - // Init I/O - pngfile = fopen(outName, "wb"); - if (!pngfile) { - error_noexit("Could not open PNG for writing"); - return 1; - } - png_init_io(png_ptr, pngfile); - png_write_info(png_ptr, info_ptr); - - // Move prow into crow, crow ~ color rows, if required - apt_rgb_t *crow[APT_MAX_HEIGHT]; - if(!greyscale){ - prow2crow(img->prow, img->nrow, palette, crow); - } - - // Apply a user provided color palette - if(chid == Palleted){ - applyUserPalette(img->prow, img->nrow, opts->palette, crow); - } - - // Precipitation overlay - if(CONTAINS(opts->effects, Precipitation_Overlay)){ - for(int y = 0; y < img->nrow; y++){ - for(int x = 0; x < APT_CH_WIDTH; x++){ - if(img->prow[y][x + APT_CHB_OFFSET] >= 198) - crow[y][x + APT_CHB_OFFSET] = crow[y][x + APT_CHA_OFFSET] = apt_applyPalette(apt_PrecipPalette, img->prow[y][x + APT_CHB_OFFSET]-198); - } - } - } - - printf("Writing %s", outName); - - // Float power macro (for gamma adjustment) - #define POWF(a, b) (b == 1.0 ? a : exp(b * log(a))) - float a = POWF(255, opts->gamma)/255; - - // Build image - for (int y = 0; y < img->nrow; y++) { - png_color pix[APT_IMG_WIDTH]; // Color - png_byte mpix[APT_IMG_WIDTH]; // Mono - - int skip = 0; - for (int x = 0; x < width; x++) { - if(crop_telemetry && x == APT_CH_WIDTH) - skip += APT_TELE_WIDTH + APT_SYNC_WIDTH + APT_SPC_WIDTH; - - if(greyscale){ - mpix[x] = POWF(img->prow[y][x + skip + offset], opts->gamma)/a; - }else{ - pix[x] = (png_color){ - POWF(crow[y][x + skip + offset].r, opts->gamma)/a, - POWF(crow[y][x + skip + offset].g, opts->gamma)/a, - POWF(crow[y][x + skip + offset].b, opts->gamma)/a - }; - } - } - - if(greyscale){ - png_write_row(png_ptr, (png_bytep) mpix); - }else{ - png_write_row(png_ptr, (png_bytep) pix); - } - } - - // Tidy up - png_write_end(png_ptr, info_ptr); - fclose(pngfile); - printf("\nDone\n"); - png_destroy_write_struct(&png_ptr, &info_ptr); - - return 1; +int ImageOut(options_t *opts, apt_image_t *img, int offset, int width, char *desc, char chid, char *palette) { + char outName[512]; + if (opts->filename == NULL || opts->filename[0] == '\0') { + sprintf(outName, "%s/%s-%c.png", opts->path, img->name, chid); + } else { + sprintf(outName, "%s/%s", opts->path, opts->filename); + } + + png_text meta[] = {{PNG_TEXT_COMPRESSION_NONE, "Software", VERSION, sizeof(VERSION)}, + {PNG_TEXT_COMPRESSION_NONE, "Channel", desc, sizeof(desc)}, + {PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20}}; + + // Parse image type + int greyscale = 1; + switch (chid) { + case Palleted: + greyscale = 0; + break; + case Temperature: + greyscale = 0; + break; + case Raw_Image: + break; + case Channel_A: + break; + case Channel_B: + break; + } + + // Parse effects + int crop_telemetry = 0; + for (unsigned long int i = 0; i < strlen(opts->effects); i++) { + switch (opts->effects[i]) { + case Crop_Telemetry: + if (width == 2080) { + width -= APT_TOTAL_TELE; + offset += APT_SYNC_WIDTH + APT_SPC_WIDTH; + crop_telemetry = 1; + } + break; + case Precipitation_Overlay: + greyscale = 0; + break; + case Flip_Image: + break; + case Denoise: + break; + case Histogram_Equalise: + break; + case Linear_Equalise: + break; + case Crop_Noise: + break; + default: { + char text[100]; + sprintf(text, "Unrecognised effect, \"%c\"", opts->effects[i]); + warning(text); + break; + } + } + } + + FILE *pngfile; + + // Create writer + png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!png_ptr) { + png_destroy_write_struct(&png_ptr, (png_infopp)NULL); + error_noexit("Could not create a PNG writer"); + return 0; + } + png_infop info_ptr = png_create_info_struct(png_ptr); + if (!info_ptr) { + png_destroy_write_struct(&png_ptr, (png_infopp)NULL); + error_noexit("Could not create a PNG writer"); + return 0; + } + + if (greyscale) { + // Greyscale image + png_set_IHDR(png_ptr, info_ptr, width, img->nrow, 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, + PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); + } else { + // 8 bit RGB image + png_set_IHDR(png_ptr, info_ptr, width, img->nrow, 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, + PNG_FILTER_TYPE_DEFAULT); + } + + png_set_text(png_ptr, info_ptr, meta, 3); + png_set_pHYs(png_ptr, info_ptr, 3636, 3636, PNG_RESOLUTION_METER); + + // Init I/O + pngfile = fopen(outName, "wb"); + if (!pngfile) { + error_noexit("Could not open PNG for writing"); + return 1; + } + png_init_io(png_ptr, pngfile); + png_write_info(png_ptr, info_ptr); + + // Move prow into crow, crow ~ color rows, if required + apt_rgb_t *crow[APT_MAX_HEIGHT]; + if (!greyscale) { + prow2crow(img->prow, img->nrow, palette, crow); + } + + // Apply a user provided color palette + if (chid == Palleted) { + applyUserPalette(img->prow, img->nrow, opts->palette, crow); + } + + // Precipitation overlay + if (CONTAINS(opts->effects, Precipitation_Overlay)) { + for (int y = 0; y < img->nrow; y++) { + for (int x = 0; x < APT_CH_WIDTH; x++) { + if (img->prow[y][x + APT_CHB_OFFSET] >= 198) + crow[y][x + APT_CHB_OFFSET] = crow[y][x + APT_CHA_OFFSET] = + apt_applyPalette(apt_PrecipPalette, img->prow[y][x + APT_CHB_OFFSET] - 198); + } + } + } + + printf("Writing %s", outName); + +// Float power macro (for gamma adjustment) +#define POWF(a, b) (b == 1.0 ? a : exp(b * log(a))) + float a = POWF(255, opts->gamma) / 255; + + // Build image + for (int y = 0; y < img->nrow; y++) { + png_color pix[APT_IMG_WIDTH]; // Color + png_byte mpix[APT_IMG_WIDTH]; // Mono + + int skip = 0; + for (int x = 0; x < width; x++) { + if (crop_telemetry && x == APT_CH_WIDTH) skip += APT_TELE_WIDTH + APT_SYNC_WIDTH + APT_SPC_WIDTH; + + if (greyscale) { + mpix[x] = POWF(img->prow[y][x + skip + offset], opts->gamma) / a; + } else { + pix[x] = (png_color){POWF(crow[y][x + skip + offset].r, opts->gamma) / a, + POWF(crow[y][x + skip + offset].g, opts->gamma) / a, + POWF(crow[y][x + skip + offset].b, opts->gamma) / a}; + } + } + + if (greyscale) { + png_write_row(png_ptr, (png_bytep)mpix); + } else { + png_write_row(png_ptr, (png_bytep)pix); + } + } + + // Tidy up + png_write_end(png_ptr, info_ptr); + fclose(pngfile); + printf("\nDone\n"); + png_destroy_write_struct(&png_ptr, &info_ptr); + + return 1; } // TODO: clean up everthing below this comment @@ -355,65 +351,61 @@ png_structp rt_png_ptr; png_infop rt_info_ptr; FILE *rt_pngfile; -int initWriter(options_t *opts, apt_image_t *img, int width, int height, char *desc, char *chid){ - char outName[384]; - sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid); - - png_text meta[] = { - {PNG_TEXT_COMPRESSION_NONE, "Software", VERSION, sizeof(VERSION)}, - {PNG_TEXT_COMPRESSION_NONE, "Channel", desc, sizeof(desc)}, - {PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20} - }; - - // Create writer - rt_png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if (!rt_png_ptr) { - png_destroy_write_struct(&rt_png_ptr, (png_infopp) NULL); - error_noexit("Could not create a PNG writer"); - return 0; - } - rt_info_ptr = png_create_info_struct(rt_png_ptr); - if (!rt_info_ptr) { - png_destroy_write_struct(&rt_png_ptr, (png_infopp) NULL); - error_noexit("Could not create a PNG writer"); - return 0; - } - - // Greyscale image - png_set_IHDR(rt_png_ptr, rt_info_ptr, width, height, - 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, - PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); - - png_set_text(rt_png_ptr, rt_info_ptr, meta, 3); - - // Channel = 25cm wide - png_set_pHYs(rt_png_ptr, rt_info_ptr, 3636, 3636, PNG_RESOLUTION_METER); - - // Init I/O - rt_pngfile = fopen(outName, "wb"); - if (!rt_pngfile) { - error_noexit("Could not open PNG for writing"); - return 0; - } - png_init_io(rt_png_ptr, rt_pngfile); - png_write_info(rt_png_ptr, rt_info_ptr); - - // Turn off compression - png_set_compression_level(rt_png_ptr, 0); - - return 1; +int initWriter(options_t *opts, apt_image_t *img, int width, int height, char *desc, char *chid) { + char outName[384]; + sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid); + + png_text meta[] = {{PNG_TEXT_COMPRESSION_NONE, "Software", VERSION, sizeof(VERSION)}, + {PNG_TEXT_COMPRESSION_NONE, "Channel", desc, sizeof(desc)}, + {PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20}}; + + // Create writer + rt_png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!rt_png_ptr) { + png_destroy_write_struct(&rt_png_ptr, (png_infopp)NULL); + error_noexit("Could not create a PNG writer"); + return 0; + } + rt_info_ptr = png_create_info_struct(rt_png_ptr); + if (!rt_info_ptr) { + png_destroy_write_struct(&rt_png_ptr, (png_infopp)NULL); + error_noexit("Could not create a PNG writer"); + return 0; + } + + // Greyscale image + png_set_IHDR(rt_png_ptr, rt_info_ptr, width, height, 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, + PNG_FILTER_TYPE_DEFAULT); + + png_set_text(rt_png_ptr, rt_info_ptr, meta, 3); + + // Channel = 25cm wide + png_set_pHYs(rt_png_ptr, rt_info_ptr, 3636, 3636, PNG_RESOLUTION_METER); + + // Init I/O + rt_pngfile = fopen(outName, "wb"); + if (!rt_pngfile) { + error_noexit("Could not open PNG for writing"); + return 0; + } + png_init_io(rt_png_ptr, rt_pngfile); + png_write_info(rt_png_ptr, rt_info_ptr); + + // Turn off compression + png_set_compression_level(rt_png_ptr, 0); + + return 1; } -void pushRow(float *row, int width){ - png_byte pix[APT_IMG_WIDTH]; - for(int i = 0; i < width; i++) - pix[i] = row[i]; +void pushRow(float *row, int width) { + png_byte pix[APT_IMG_WIDTH]; + for (int i = 0; i < width; i++) pix[i] = row[i]; - png_write_row(rt_png_ptr, (png_bytep) pix); + png_write_row(rt_png_ptr, (png_bytep)pix); } -void closeWriter(){ - png_write_end(rt_png_ptr, rt_info_ptr); - fclose(rt_pngfile); - png_destroy_write_struct(&rt_png_ptr, &rt_info_ptr); +void closeWriter() { + png_write_end(rt_png_ptr, rt_info_ptr); + fclose(rt_pngfile); + png_destroy_write_struct(&rt_png_ptr, &rt_info_ptr); } diff --git a/src/taps.h b/src/taps.h index 9d5f5cf..f33b892 100644 --- a/src/taps.h +++ b/src/taps.h @@ -16,77 +16,65 @@ * along with this program. If not, see . */ -static const float hilbert_filter[] = { 0.0205361, 0.0219524, 0.0235785, 0.0254648, 0.0276791, 0.0303152, -0.0335063, 0.0374482, 0.0424413, 0.0489708, 0.0578745, 0.0707355, 0.0909457, 0.127324, 0.212207, 0.63662, --0.63662, -0.212207, -0.127324, -0.0909457, -0.0707355, -0.0578745, -0.0489708, -0.0424413, -0.0374482, --0.0335063, -0.0303152, -0.0276791, -0.0254648, -0.0235785, -0.0219524, -0.0205361 }; -#define HILBERT_FILTER_SIZE (sizeof(hilbert_filter)/sizeof(hilbert_filter[0])) +static const float hilbert_filter[] = {0.0205361, 0.0219524, 0.0235785, 0.0254648, 0.0276791, 0.0303152, 0.0335063, + 0.0374482, 0.0424413, 0.0489708, 0.0578745, 0.0707355, 0.0909457, 0.127324, + 0.212207, 0.63662, -0.63662, -0.212207, -0.127324, -0.0909457, -0.0707355, + -0.0578745, -0.0489708, -0.0424413, -0.0374482, -0.0335063, -0.0303152, -0.0276791, + -0.0254648, -0.0235785, -0.0219524, -0.0205361}; +#define HILBERT_FILTER_SIZE (sizeof(hilbert_filter) / sizeof(hilbert_filter[0])) -static const float low_pass[] = { -3.37279e-04, -8.80292e-06, -3.96418e-04, -1.78544e-04, -5.27511e-04, --3.75376e-04, -6.95337e-04, -5.93148e-04, -8.79730e-04, -8.15327e-04, -1.05669e-03, -1.01377e-03, --1.19836e-03, -1.15443e-03, -1.26937e-03, -1.20955e-03, -1.23904e-03, -1.15302e-03, -1.08660e-03, --9.64235e-04, -8.02450e-04, -6.46202e-04, -3.95376e-04, -2.18096e-04, 1.11906e-04, 2.89567e-04, -6.67167e-04, 8.19039e-04, 1.21725e-03, 1.30556e-03, 1.69365e-03, 1.68588e-03, 2.03277e-03, 1.90159e-03, -2.18455e-03, 1.90833e-03, 2.12100e-03, 1.69052e-03, 1.77484e-03, 1.42542e-03, 1.18292e-03, 8.66979e-04, -5.54161e-04, 2.15793e-04, -1.11623e-04, -4.35173e-04, -7.27194e-04, -9.91551e-04, -1.20407e-03, --1.37032e-03, -1.46991e-03, -1.51120e-03, -1.48008e-03, -1.39047e-03, -1.23115e-03, -1.02128e-03, --7.60099e-04, -4.68008e-04, -1.46339e-04, 1.80867e-04, 5.11244e-04, 8.19243e-04, 1.09739e-03, 1.32668e-03, -1.50632e-03, 1.61522e-03, 1.66246e-03, 1.62390e-03, 1.52430e-03, 1.34273e-03, 1.10736e-03, 8.10335e-04, -4.76814e-04, 1.13622e-04, -2.64150e-04, -6.26595e-04, -9.95436e-04, -1.27846e-03, -1.54080e-03, --1.74292e-03, -1.86141e-03, -1.89318e-03, -1.83969e-03, -1.69770e-03, -1.47938e-03, -1.18696e-03, --8.37003e-04, -4.39507e-04, -1.56907e-05, 4.19904e-04, 8.43172e-04, 1.23827e-03, 1.58411e-03, -1.86382e-03, 2.06312e-03, 2.17177e-03, 2.18121e-03, 2.08906e-03, 1.89772e-03, 1.61153e-03, -1.24507e-03, 8.13976e-04, 3.29944e-04, -1.74591e-04, -6.83619e-04, -1.17826e-03, -1.61659e-03, --2.00403e-03, -2.29070e-03, -2.49179e-03, -2.56546e-03, -2.53448e-03, -2.37032e-03, -2.10060e-03, --1.72140e-03, -1.24542e-03, -7.15425e-04, -1.24964e-04, 4.83736e-04, 1.08328e-03, 1.64530e-03, -2.14503e-03, 2.55400e-03, 2.85589e-03, 3.02785e-03, 3.06271e-03, 2.95067e-03, 2.69770e-03, -2.30599e-03, 1.79763e-03, 1.18587e-03, 5.04003e-04, -2.23591e-04, -9.57591e-04, -1.66939e-03, --2.31717e-03, -2.87636e-03, -3.31209e-03, -3.60506e-03, -3.73609e-03, -3.69208e-03, -3.44913e-03, --3.06572e-03, -2.50229e-03, -1.80630e-03, -1.00532e-03, -1.22305e-04, 7.83910e-04, 1.69402e-03, -2.53826e-03, 3.30312e-03, 3.91841e-03, 4.38017e-03, 4.63546e-03, 4.68091e-03, 4.50037e-03, -4.09614e-03, 3.47811e-03, 2.67306e-03, 1.70418e-03, 6.20542e-04, -5.36994e-04, -1.70981e-03, --2.84712e-03, -3.88827e-03, -4.78659e-03, -5.48593e-03, -5.95049e-03, -6.14483e-03, -6.05118e-03, --5.65829e-03, -4.97525e-03, -4.01796e-03, -2.82224e-03, -1.43003e-03, 1.00410e-04, 1.71169e-03, -3.31983e-03, 4.87796e-03, 6.23237e-03, 7.31013e-03, 8.20642e-03, 8.67374e-03, 8.77681e-03, -8.43444e-03, 7.66794e-03, 6.46827e-03, 4.87294e-03, 2.92923e-03, 6.98913e-04, -1.72126e-03, --4.24785e-03, -6.75380e-03, -9.13309e-03, -1.12532e-02, -1.30038e-02, -1.42633e-02, -1.49338e-02, --1.49145e-02, -1.41484e-02, -1.25761e-02, -1.01870e-02, -6.97432e-03, -2.97910e-03, 1.75386e-03, -7.11899e-03, 1.30225e-02, 1.93173e-02, 2.58685e-02, 3.24965e-02, 3.90469e-02, 4.53316e-02, -5.11931e-02, 5.64604e-02, 6.09924e-02, 6.46584e-02, 6.73547e-02, 6.90049e-02, 6.97096e-02, -6.90049e-02, 6.73547e-02, 6.46584e-02, 6.09924e-02, 5.64604e-02, 5.11931e-02, 4.53316e-02, -3.90469e-02, 3.24965e-02, 2.58685e-02, 1.93173e-02, 1.30225e-02, 7.11899e-03, 1.75386e-03, --2.97910e-03, -6.97432e-03, -1.01870e-02, -1.25761e-02, -1.41484e-02, -1.49145e-02, -1.49338e-02, --1.42633e-02, -1.30038e-02, -1.12532e-02, -9.13309e-03, -6.75380e-03, -4.24785e-03, -1.72126e-03, -6.98913e-04, 2.92923e-03, 4.87294e-03, 6.46827e-03, 7.66794e-03, 8.43444e-03, 8.77681e-03, -8.67374e-03, 8.20642e-03, 7.31013e-03, 6.23237e-03, 4.87796e-03, 3.31983e-03, 1.71169e-03, -1.00410e-04, -1.43003e-03, -2.82224e-03, -4.01796e-03, -4.97525e-03, -5.65829e-03, -6.05118e-03, --6.14483e-03, -5.95049e-03, -5.48593e-03, -4.78659e-03, -3.88827e-03, -2.84712e-03, -1.70981e-03, --5.36994e-04, 6.20542e-04, 1.70418e-03, 2.67306e-03, 3.47811e-03, 4.09614e-03, 4.50037e-03, -4.68091e-03, 4.63546e-03, 4.38017e-03, 3.91841e-03, 3.30312e-03, 2.53826e-03, 1.69402e-03, -7.83910e-04, -1.22305e-04, -1.00532e-03, -1.80630e-03, -2.50229e-03, -3.06572e-03, -3.44913e-03, --3.69208e-03, -3.73609e-03, -3.60506e-03, -3.31209e-03, -2.87636e-03, -2.31717e-03, -1.66939e-03, --9.57591e-04, -2.23591e-04, 5.04003e-04, 1.18587e-03, 1.79763e-03, 2.30599e-03, 2.69770e-03, -2.95067e-03, 3.06271e-03, 3.02785e-03, 2.85589e-03, 2.55400e-03, 2.14503e-03, 1.64530e-03, -1.08328e-03, 4.83736e-04, -1.24964e-04, -7.15425e-04, -1.24542e-03, -1.72140e-03, -2.10060e-03, --2.37032e-03, -2.53448e-03, -2.56546e-03, -2.49179e-03, -2.29070e-03, -2.00403e-03, -1.61659e-03, --1.17826e-03, -6.83619e-04, -1.74591e-04, 3.29944e-04, 8.13976e-04, 1.24507e-03, 1.61153e-03, -1.89772e-03, 2.08906e-03, 2.18121e-03, 2.17177e-03, 2.06312e-03, 1.86382e-03, 1.58411e-03, -1.23827e-03, 8.43172e-04, 4.19904e-04, -1.56907e-05, -4.39507e-04, -8.37003e-04, -1.18696e-03, --1.47938e-03, -1.69770e-03, -1.83969e-03, -1.89318e-03, -1.86141e-03, -1.74292e-03, -1.54080e-03, --1.27846e-03, -9.95436e-04, -6.26595e-04, -2.64150e-04, 1.13622e-04, 4.76814e-04, 8.10335e-04, -1.10736e-03, 1.34273e-03, 1.52430e-03, 1.62390e-03, 1.66246e-03, 1.61522e-03, 1.50632e-03, -1.32668e-03, 1.09739e-03, 8.19243e-04, 5.11244e-04, 1.80867e-04, -1.46339e-04, -4.68008e-04, --7.60099e-04, -1.02128e-03, -1.23115e-03, -1.39047e-03, -1.48008e-03, -1.51120e-03, -1.46991e-03, --1.37032e-03, -1.20407e-03, -9.91551e-04, -7.27194e-04, -4.35173e-04, -1.11623e-04, 2.15793e-04, -5.54161e-04, 8.66979e-04, 1.18292e-03, 1.42542e-03, 1.77484e-03, 1.69052e-03, 2.12100e-03, -1.90833e-03, 2.18455e-03, 1.90159e-03, 2.03277e-03, 1.68588e-03, 1.69365e-03, 1.30556e-03, -1.21725e-03, 8.19039e-04, 6.67167e-04, 2.89567e-04, 1.11906e-04, -2.18096e-04, -3.95376e-04, --6.46202e-04, -8.02450e-04, -9.64235e-04, -1.08660e-03, -1.15302e-03, -1.23904e-03, -1.20955e-03, --1.26937e-03, -1.15443e-03, -1.19836e-03, -1.01377e-03, -1.05669e-03, -8.15327e-04, -8.79730e-04, --5.93148e-04, -6.95337e-04, -3.75376e-04, -5.27511e-04, -1.78544e-04, -3.96418e-04, -8.80292e-06, --3.37279e-04 }; -#define LOW_PASS_SIZE (sizeof(low_pass)/sizeof(low_pass[0])) +static const float low_pass[] = { + -3.37279e-04, -8.80292e-06, -3.96418e-04, -1.78544e-04, -5.27511e-04, -3.75376e-04, -6.95337e-04, -5.93148e-04, -8.79730e-04, + -8.15327e-04, -1.05669e-03, -1.01377e-03, -1.19836e-03, -1.15443e-03, -1.26937e-03, -1.20955e-03, -1.23904e-03, -1.15302e-03, + -1.08660e-03, -9.64235e-04, -8.02450e-04, -6.46202e-04, -3.95376e-04, -2.18096e-04, 1.11906e-04, 2.89567e-04, 6.67167e-04, + 8.19039e-04, 1.21725e-03, 1.30556e-03, 1.69365e-03, 1.68588e-03, 2.03277e-03, 1.90159e-03, 2.18455e-03, 1.90833e-03, + 2.12100e-03, 1.69052e-03, 1.77484e-03, 1.42542e-03, 1.18292e-03, 8.66979e-04, 5.54161e-04, 2.15793e-04, -1.11623e-04, + -4.35173e-04, -7.27194e-04, -9.91551e-04, -1.20407e-03, -1.37032e-03, -1.46991e-03, -1.51120e-03, -1.48008e-03, -1.39047e-03, + -1.23115e-03, -1.02128e-03, -7.60099e-04, -4.68008e-04, -1.46339e-04, 1.80867e-04, 5.11244e-04, 8.19243e-04, 1.09739e-03, + 1.32668e-03, 1.50632e-03, 1.61522e-03, 1.66246e-03, 1.62390e-03, 1.52430e-03, 1.34273e-03, 1.10736e-03, 8.10335e-04, + 4.76814e-04, 1.13622e-04, -2.64150e-04, -6.26595e-04, -9.95436e-04, -1.27846e-03, -1.54080e-03, -1.74292e-03, -1.86141e-03, + -1.89318e-03, -1.83969e-03, -1.69770e-03, -1.47938e-03, -1.18696e-03, -8.37003e-04, -4.39507e-04, -1.56907e-05, 4.19904e-04, + 8.43172e-04, 1.23827e-03, 1.58411e-03, 1.86382e-03, 2.06312e-03, 2.17177e-03, 2.18121e-03, 2.08906e-03, 1.89772e-03, + 1.61153e-03, 1.24507e-03, 8.13976e-04, 3.29944e-04, -1.74591e-04, -6.83619e-04, -1.17826e-03, -1.61659e-03, -2.00403e-03, + -2.29070e-03, -2.49179e-03, -2.56546e-03, -2.53448e-03, -2.37032e-03, -2.10060e-03, -1.72140e-03, -1.24542e-03, -7.15425e-04, + -1.24964e-04, 4.83736e-04, 1.08328e-03, 1.64530e-03, 2.14503e-03, 2.55400e-03, 2.85589e-03, 3.02785e-03, 3.06271e-03, + 2.95067e-03, 2.69770e-03, 2.30599e-03, 1.79763e-03, 1.18587e-03, 5.04003e-04, -2.23591e-04, -9.57591e-04, -1.66939e-03, + -2.31717e-03, -2.87636e-03, -3.31209e-03, -3.60506e-03, -3.73609e-03, -3.69208e-03, -3.44913e-03, -3.06572e-03, -2.50229e-03, + -1.80630e-03, -1.00532e-03, -1.22305e-04, 7.83910e-04, 1.69402e-03, 2.53826e-03, 3.30312e-03, 3.91841e-03, 4.38017e-03, + 4.63546e-03, 4.68091e-03, 4.50037e-03, 4.09614e-03, 3.47811e-03, 2.67306e-03, 1.70418e-03, 6.20542e-04, -5.36994e-04, + -1.70981e-03, -2.84712e-03, -3.88827e-03, -4.78659e-03, -5.48593e-03, -5.95049e-03, -6.14483e-03, -6.05118e-03, -5.65829e-03, + -4.97525e-03, -4.01796e-03, -2.82224e-03, -1.43003e-03, 1.00410e-04, 1.71169e-03, 3.31983e-03, 4.87796e-03, 6.23237e-03, + 7.31013e-03, 8.20642e-03, 8.67374e-03, 8.77681e-03, 8.43444e-03, 7.66794e-03, 6.46827e-03, 4.87294e-03, 2.92923e-03, + 6.98913e-04, -1.72126e-03, -4.24785e-03, -6.75380e-03, -9.13309e-03, -1.12532e-02, -1.30038e-02, -1.42633e-02, -1.49338e-02, + -1.49145e-02, -1.41484e-02, -1.25761e-02, -1.01870e-02, -6.97432e-03, -2.97910e-03, 1.75386e-03, 7.11899e-03, 1.30225e-02, + 1.93173e-02, 2.58685e-02, 3.24965e-02, 3.90469e-02, 4.53316e-02, 5.11931e-02, 5.64604e-02, 6.09924e-02, 6.46584e-02, + 6.73547e-02, 6.90049e-02, 6.97096e-02, 6.90049e-02, 6.73547e-02, 6.46584e-02, 6.09924e-02, 5.64604e-02, 5.11931e-02, + 4.53316e-02, 3.90469e-02, 3.24965e-02, 2.58685e-02, 1.93173e-02, 1.30225e-02, 7.11899e-03, 1.75386e-03, -2.97910e-03, + -6.97432e-03, -1.01870e-02, -1.25761e-02, -1.41484e-02, -1.49145e-02, -1.49338e-02, -1.42633e-02, -1.30038e-02, -1.12532e-02, + -9.13309e-03, -6.75380e-03, -4.24785e-03, -1.72126e-03, 6.98913e-04, 2.92923e-03, 4.87294e-03, 6.46827e-03, 7.66794e-03, + 8.43444e-03, 8.77681e-03, 8.67374e-03, 8.20642e-03, 7.31013e-03, 6.23237e-03, 4.87796e-03, 3.31983e-03, 1.71169e-03, + 1.00410e-04, -1.43003e-03, -2.82224e-03, -4.01796e-03, -4.97525e-03, -5.65829e-03, -6.05118e-03, -6.14483e-03, -5.95049e-03, + -5.48593e-03, -4.78659e-03, -3.88827e-03, -2.84712e-03, -1.70981e-03, -5.36994e-04, 6.20542e-04, 1.70418e-03, 2.67306e-03, + 3.47811e-03, 4.09614e-03, 4.50037e-03, 4.68091e-03, 4.63546e-03, 4.38017e-03, 3.91841e-03, 3.30312e-03, 2.53826e-03, + 1.69402e-03, 7.83910e-04, -1.22305e-04, -1.00532e-03, -1.80630e-03, -2.50229e-03, -3.06572e-03, -3.44913e-03, -3.69208e-03, + -3.73609e-03, -3.60506e-03, -3.31209e-03, -2.87636e-03, -2.31717e-03, -1.66939e-03, -9.57591e-04, -2.23591e-04, 5.04003e-04, + 1.18587e-03, 1.79763e-03, 2.30599e-03, 2.69770e-03, 2.95067e-03, 3.06271e-03, 3.02785e-03, 2.85589e-03, 2.55400e-03, + 2.14503e-03, 1.64530e-03, 1.08328e-03, 4.83736e-04, -1.24964e-04, -7.15425e-04, -1.24542e-03, -1.72140e-03, -2.10060e-03, + -2.37032e-03, -2.53448e-03, -2.56546e-03, -2.49179e-03, -2.29070e-03, -2.00403e-03, -1.61659e-03, -1.17826e-03, -6.83619e-04, + -1.74591e-04, 3.29944e-04, 8.13976e-04, 1.24507e-03, 1.61153e-03, 1.89772e-03, 2.08906e-03, 2.18121e-03, 2.17177e-03, + 2.06312e-03, 1.86382e-03, 1.58411e-03, 1.23827e-03, 8.43172e-04, 4.19904e-04, -1.56907e-05, -4.39507e-04, -8.37003e-04, + -1.18696e-03, -1.47938e-03, -1.69770e-03, -1.83969e-03, -1.89318e-03, -1.86141e-03, -1.74292e-03, -1.54080e-03, -1.27846e-03, + -9.95436e-04, -6.26595e-04, -2.64150e-04, 1.13622e-04, 4.76814e-04, 8.10335e-04, 1.10736e-03, 1.34273e-03, 1.52430e-03, + 1.62390e-03, 1.66246e-03, 1.61522e-03, 1.50632e-03, 1.32668e-03, 1.09739e-03, 8.19243e-04, 5.11244e-04, 1.80867e-04, + -1.46339e-04, -4.68008e-04, -7.60099e-04, -1.02128e-03, -1.23115e-03, -1.39047e-03, -1.48008e-03, -1.51120e-03, -1.46991e-03, + -1.37032e-03, -1.20407e-03, -9.91551e-04, -7.27194e-04, -4.35173e-04, -1.11623e-04, 2.15793e-04, 5.54161e-04, 8.66979e-04, + 1.18292e-03, 1.42542e-03, 1.77484e-03, 1.69052e-03, 2.12100e-03, 1.90833e-03, 2.18455e-03, 1.90159e-03, 2.03277e-03, + 1.68588e-03, 1.69365e-03, 1.30556e-03, 1.21725e-03, 8.19039e-04, 6.67167e-04, 2.89567e-04, 1.11906e-04, -2.18096e-04, + -3.95376e-04, -6.46202e-04, -8.02450e-04, -9.64235e-04, -1.08660e-03, -1.15302e-03, -1.23904e-03, -1.20955e-03, -1.26937e-03, + -1.15443e-03, -1.19836e-03, -1.01377e-03, -1.05669e-03, -8.15327e-04, -8.79730e-04, -5.93148e-04, -6.95337e-04, -3.75376e-04, + -5.27511e-04, -1.78544e-04, -3.96418e-04, -8.80292e-06, -3.37279e-04}; +#define LOW_PASS_SIZE (sizeof(low_pass) / sizeof(low_pass[0])) -static const float sync_pattern[] = { -14, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -14 }; -#define SYNC_PATTERN_SIZE (sizeof(sync_pattern)/sizeof(sync_pattern[0])) +static const float sync_pattern[] = {-14, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, + 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -14}; +#define SYNC_PATTERN_SIZE (sizeof(sync_pattern) / sizeof(sync_pattern[0])) diff --git a/src/util.c b/src/util.c index a808c39..612267a 100644 --- a/src/util.c +++ b/src/util.c @@ -46,6 +46,4 @@ float clamp(float x, float hi, float lo) { return x; } -float clamp_half(float x, float hi) { - return clamp(x, hi, -hi); -} +float clamp_half(float x, float hi) { return clamp(x, hi, -hi); } diff --git a/src/util.h b/src/util.h index cebcaca..7fd01ec 100644 --- a/src/util.h +++ b/src/util.h @@ -20,8 +20,8 @@ #define M_PIf 3.14159265358979323846f #define M_TAUf (M_PIf * 2.0f) -#define MIN(a,b) (((a) < (b)) ? (a) : (b)) -#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define MIN(a, b) (((a) < (b)) ? (a) : (b)) +#define MAX(a, b) (((a) > (b)) ? (a) : (b)) #ifndef UTIL_H #define UTIL_H