hms
/
aptdec
mirror of https://github.com/Xerbo/aptdec
1
0
Fork 0
Code Issues 0 Releases 7 Wiki Activity
Browse Source

Format codebase

tags/v1.8.0
Xerbo 1 year ago
parent
3e8297a912
commit
4d4a0c9787
No known key found for this signature in database GPG Key ID: 34103F6D8F11CEB0
17 changed files with 1376 additions and 1442 deletions
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +4
    -0
      .clang-format
  2. +5
    -8
      src/algebra.c
  3. +29
    -21
      src/apt.h
  4. +75
    -94
      src/calibration.c
  5. +3
    -3
      src/calibration.h
  6. +49
    -56
      src/color.c
  7. +1
    -1
      src/color.h
  8. +24
    -24
      src/common.h
  9. +176
    -182
      src/dsp.c
  10. +25
    -23
      src/filter.c
  11. +3
    -3
      src/filter.h
  12. +301
    -311
      src/image.c
  13. +254
    -267
      src/main.c
  14. +364
    -372
      src/pngio.c
  15. +60
    -72
      src/taps.h
  16. +1
    -3
      src/util.c
  17. +2
    -2
      src/util.h

+ 4
- 0
.clang-format View File

@@ -0,0 +1,4 @@
Language: Cpp
BasedOnStyle: Google
IndentWidth: 4
ColumnLimit: 130 # Conservative estimate considering the size of modern displays

+ 5
- 8
src/algebra.c View File

@@ -17,6 +17,7 @@
*/

#include "algebra.h"

#include <math.h>

// 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; }

+ 29
- 21
src/apt.h View File

@@ -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
}


+ 75
- 94
src/calibration.c View File

@@ -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];
}
}

+ 3
- 3
src/calibration.h View File

@@ -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


+ 49
- 56
src/color.c View File

@@ -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"};

+ 1
- 1
src/color.h View File

@@ -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))

+ 24
- 24
src/common.h View File

@@ -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

+ 176
- 182
src/dsp.c View File

@@ -16,10 +16,10 @@
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*/

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#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;
}

+ 25
- 23
src/filter.c View File

@@ -16,45 +16,47 @@
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*/

#include <math.h>
#include "filter.h"

#include <math.h>

#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;
}

+ 3
- 3
src/filter.h View File

@@ -16,8 +16,8 @@
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*/
#include <stddef.h>
#include <complex.h>
#include <stddef.h>
#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);

+ 301
- 311
src/image.c View File

@@ -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 <stdio.h>
#include <string.h>
#include "image.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#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);
}
}
}

+ 254
- 267
src/main.c View File

@@ -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 <stdlib.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#ifndef _MSC_VER
#include <libgen.h>
#else
#include <windows.h>
#endif
#include <errno.h>
#include <math.h>
#include <sndfile.h>
#include <errno.h>
#include <time.h>
#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);
}
}

+ 364
- 372
src/pngio.c View File

@@ -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 <math.h>
#include <png.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <math.h>
#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);
}

+ 60
- 72
src/taps.h View File

@@ -16,77 +16,65 @@
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*/

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

+ 1
- 3
src/util.c View File

@@ -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); }

+ 2
- 2
src/util.h View File

@@ -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


Write Preview
Loading…
Cancel
Save