diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a75882b --- /dev/null +++ b/.gitignore @@ -0,0 +1,53 @@ +# Created by https://www.gitignore.io/api/c +# Edit at https://www.gitignore.io/?templates=c + +# Object files +*.o +*.ko +*.obj +*.elf + +# Linker output +*.ilk +*.map +*.exp + +# Precompiled Headers +*.gch +*.pch + +# Libraries +*.lib +*.a +*.la +*.lo + +# Shared objects (inc. Windows DLLs) +*.dll +*.so +*.so.* +*.dylib + +# Executables +*.exe +*.out +*.app +*.i*86 +*.x86_64 +*.hex + +# Debug files +*.dSYM/ +*.su +*.idb +*.pdb + +# Program specifics +*.png +*.wav +!textlogo.png +*.pnm +aptdec + +# VSCode +.vscode \ No newline at end of file diff --git a/Makefile b/Makefile index b521492..c6df7be 100644 --- a/Makefile +++ b/Makefile @@ -1,16 +1,16 @@ -CC=gcc -BIN=/usr/bin -INCLUDES=-I. -CFLAGS= -O3 -Wall $(INCLUDES) -OBJS= main.o image.o dsp.o filter.o reg.o fcolor.o +CC = gcc +BIN = /usr/bin +INCLUDES = -I. +CFLAGS = -O3 -Wall $(INCLUDES) +OBJS = main.o image.o dsp.o filter.o reg.o fcolor.o aptdec: $(OBJS) $(CC) -o $@ $(OBJS) -lm -lsndfile -lpng -main.o: main.c version.h temppalette.h gvipalette.h offsets.h +main.o: main.c temppalette.h gvipalette.h offsets.h messages.h dsp.o: dsp.c filtercoeff.h filter.h filter.o: filter.c filter.h -image.o: image.c satcal.h offsets.h +image.o: image.c satcal.h offsets.h messages.h fcolor.o: fcolor.c offsets.h clean: diff --git a/README b/README deleted file mode 100644 index 1b6b6d2..0000000 --- a/README +++ /dev/null @@ -1,79 +0,0 @@ - -ATPDEC README (Thierry Leconte F4DWV (c) 2004-2009) - -DESCRIPTION - -Atpdec is an open source program that decodes images transmitted by POES -NOAA weather satellite series. -These satellites transmit continuously, among other things, medium -resolution images of the earth on 137Mhz. -These transmissions could be easily received with an inexpensive antenna -and dedicated receiver. -Output from such a receiver, is an audio signal that could be recorded -into a soundfile with any soundcard. - -Atpdec will convert these sounfiles into .png images. - -For each soundfile up to 6 images could be generated : - - 1. Raw image : contains the 2 transmitted channel images + telemetry -and synchro pulses. - 2. Calibrated channel A image - 3. Calibrated channel B image - 4. Temperature compensed I.R image - 5. False color image - -Input soundfiles must be mono signal sampled at 11025 Hz. -Atpdec use libsndfile to read soundfile, so any sound file format supported by libsndfile -could be read.(Only tested with .wav file). - -COMPILATION - -Atpdec is written in plain standart C and must be very portable. -It was only tested on Linux Fedora , but must work on any Unix platform. -Just adapt the Makefile and type make (sorry no configure). - -atpdec use libsndfile, libpng and libm. -snd.h and png.h header must be present on your system. -If they are not on standard path, edit the include path in the Makefile. - -USAGE - -atpdec [options] soundfiles ... - -OPTIONS - - -i [r|a|b|c|t] - Toggle raw (r) , channel A (a) , channel B (b) , false color (c) , - or temperature (t) output. - Default : "ac" - - -d directory - Optional images destination directory. - Default : soundfile directory. - - -s n - Satellite number 15 to 19 - Used for Temperature compensation. - Default : NOAA-19 - - -c conf_file - Use configuration file for false color generation. - Default : Internal parameters. - -OUTPUT - -Generated image are in png format, 8bits greyscale for raw and channel A|B images, -24bits RVB for false color. - -Image names are soundfilename-x.png, where x is : - -r for raw images - -satellite instrument number (1,2,3A,3B,4,5) for channel A|B images - -c for false colors. - -EXAMPLE - -atpdec -d image -i ac *.wav - -Will process all .wav files in the current directory, generate only channel A and false color images and put them in the image directory. - diff --git a/README.md b/README.md index 977e21a..ed449e5 100644 --- a/README.md +++ b/README.md @@ -1,33 +1,31 @@ -# Aptdec +![Aptdec logo](textlogo.png) Thierry Leconte F4DWV (c) 2004-2009 ## Description -Aptec is an FOSS program that decodes images transmitted by POES NOAA weather satellites. -These satellites transmit continuously (among other things), medium resolution (1px/4km) images of the earth on 137 MHz. -These transmissions could be easily received with an simple antenna and cheap SDR. -Then the transimssion can easily be decoded in narrow band FM mode. +Aptdec is an FOSS program that decodes images transmitted by NOAA weather satellites. These satellites transmit continuously (among other things), medium resolution (1px/4km) images of the earth on 137 MHz. +These transmissions could be easily received with an simple antenna and cheap SDR. Then the transmission can easily be decoded in narrow band FM mode. -Aptdec can convert these recordings into .png images. +Aptdec can convert these audio files into .png images. -For each recording up to 6 images can be generated: +For each audio file up to 6 images can be generated: -1. Raw image: contains the 2 transmitted channel images + telemetry -and synchronisation pulses. +1. Raw image: contains the 2 transmitted channel images + telemetry and synchronization pulses. 2. Calibrated channel A image 3. Calibrated channel B image -4. Temperature compensated I.R image +4. Temperature compensated IR image 5. False color image +6. Layered image, boosts cloud visibility -Input recordings must be mono with a sample rate of 11025 Hz. -Aptdec uses `libsndfile` to read the input recording, so any format supported by `libsndfile` may be used (only tested with .wav files). +The input audio file must be mono with a sample rate in the range of 4160-62400 Hz, lower samples rates will process faster. +Aptdec uses `libsndfile` to read the input audio, so any format supported by `libsndfile` may be used (however it has only tested with `.wav` files). ## Compilation -Aptdec is written is portable since it is written in standard C. -It has only tested on Fedora and Debian, but will work on any Unix platform. -Just edit the Makefile and run `make` (no configure script as of now). +Aptdec is portable since it is written in standard C. +It has successfully compiled and ran on Debian with both `gcc` and `clang` and will most likely work on any Unix platform. +Just edit the Makefile and run `make` (no configure script as of right now). Aptdec uses `libsndfile`, `libpng` and `libm`. The `snd.h` and `png.h` headers must be present on your system. @@ -35,55 +33,75 @@ If they are not on standard path, edit the include path in the Makefile. ## Usage +To compile +`make` + To run without installing -`./Aptdec [options] recordings ...` +`./aptdec [options] audio files...` To install `sudo make install` +To run once installed +`aptdec [options] audio files...` + To uninstall `sudo make uninstall` -To run once installed -`Aptdec [options] recordings ...` - ## Options ``` --i [r|a|b|c|t] -Toggle raw (r), channel A (a), channel B (b), false color (c), -or temperature (t) output. -Default: ac +-i [r|a|b|c|t|l] +Output image type +Raw (r), Channel A (a), Channel B (b), False Color (c), Temperature (t), Layered (l) +Default: ab -d -Optional images destination directory. -Default: Recording directory. +Images destination directory (optional). +Default: Current directory -s [15|16|17|18|19] Satellite number For temperature and false color generation Default: 19 +-e [c|t] +Enhancements +Contrast (c) or Crop Telemetry (t) +Defaults: ct + -c Use configuration file for false color generation. -Default: Internal parameters. +Default: Internal parameters ``` ## Output Generated images are outputted in PNG, 8 bit greyscale for raw and channel A|B images, 24 bit RGB for false color. -Image names are `recordingname-x.png`, where `x` is: +Image names are `audiofile-x.png`, where `x` is: + + - `r` for raw images + - Sensor ID (1, 2, 3A, 3B, 4, 5) for channel A|B images + - `c` for false color. + - `t` for temperature calibrated images + - `l` for layered images + +Currently there are 2 available enchancements: - - r for raw images - - satellite instrument number (1, 2, 3A, 3B, 4, 5) for channel A|B images - - c for false colors. + - `c` for contrast enhancements, on by default + - `t` for crop telemetry, on by default, only has effects on raw images ## Example -`Aptdec -d images -i ac *.wav` +`aptdec -d images -i ab *.wav` + +This will process all `.wav` files in the current directory, generate contrast enhanced channel A and B images and put them in the `images` directory. + +## Further reading -This will process all .wav files in the current directory, generate channel A and false color images and put them in the `images` directory. +[https://noaasis.noaa.gov/NOAASIS/pubs/Users_Guide-Building_Receive_Stations_March_2009.pdf](https://noaasis.noaa.gov/NOAASIS/pubs/Users_Guide-Building_Receive_Stations_March_2009.pdf) +[https://web.archive.org/web/20141220021557/https://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/tables.htm](https://web.archive.org/web/20141220021557/https://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/tables.htm) ## License diff --git a/dsp.c b/dsp.c index 670b3bf..c9ec059 100755 --- a/dsp.c +++ b/dsp.c @@ -1,5 +1,5 @@ /* - * Aptec + * Aptdec * Copyright (c) 2004 by Thierry Leconte (F4DWV) * * $Id$ @@ -20,17 +20,25 @@ * */ #include +#include #include #include #ifndef M_PI -#define M_PI 3.14159265358979323846 /* for OS's that don't include it */ +#define M_PI 3.1415926535 // For OS's that don't include it #endif #include "filter.h" #include "filtercoeff.h" +/* WARNING + * The comments in this file are at best educated guesses. + * I am not a DSP god and can not figure out what a complicated + * math function does just by looking at it. + */ + extern int getsample(float *inbuff, int nb); #define BLKAMP 1024 +// Block size #define BLKIN 1024 #define Fc 2400.0 @@ -48,10 +56,10 @@ static double FreqLine = 1.0; static double FreqOsc; static double K1, K2; - +// Check if the input sample rate is correct int init_dsp(double F) { - if(F > Fi) return (1); - if(F < Fp) return (-1); + if(F > Fi) return(1); + if(F < Fp) return(-1); Fe = F; K1 = DFc / Fe; @@ -61,9 +69,10 @@ int init_dsp(double F) { return(0); } -/* fast phase estimator */ -static inline double Phase(double I,double Q) { - double angle,r; +// Fast phase estimator +// Calculates the phase angle of a signal from a IQ sample +static inline double Phase(double I, double Q) { + double angle, r; int s; if(I == 0.0 && Q == 0.0) @@ -90,68 +99,79 @@ static inline double Phase(double I,double Q) { return(-angle); } +// Phase locked loop +// Used to get value from an IQ sample with noise reduction static double pll(double I, double Q) { - /* pll coeff */ + // PLL coefficient static double PhaseOsc = 0.0; double Io, Qo; double Ip, Qp; double DPhi; - /* quadrature oscillator */ + // Quadrature oscillator / reference Io = cos(PhaseOsc); Qo = sin(PhaseOsc); - /* phase detector */ + // Phase detector Ip = I * Io + Q * Qo; Qp = Q * Io - I * Qo; DPhi = Phase(Ip, Qp); - /* loop filter */ + // Loop filter PhaseOsc += 2.0 * M_PI * (K1 * DPhi + FreqOsc); if (PhaseOsc > M_PI) - PhaseOsc -= 2.0 * M_PI; + PhaseOsc -= 2.0 * M_PI; if (PhaseOsc <= -M_PI) - PhaseOsc += 2.0 * M_PI; + PhaseOsc += 2.0 * M_PI; FreqOsc += K2 * DPhi; if (FreqOsc > ((Fc + DFc) / Fe)) - FreqOsc = (Fc + DFc) / Fe; + FreqOsc = (Fc + DFc) / Fe; if (FreqOsc < ((Fc - DFc) / Fe)) - FreqOsc = (Fc - DFc) / Fe; + FreqOsc = (Fc - DFc) / Fe; - return (Ip); + return(Ip); } +// Convert audio samples into a pixel static int getamp(double *ambuff, int nb) { static float inbuff[BLKIN]; - static int idxin=0; - static int nin=0; - - int n; + static int idxin = 0; + static int nin = 0; + int n; for (n = 0; n < nb; n++) { double I, Q; + // If the amount of samples is small enough to be processed if (nin < IQFilterLen * 2 + 2) { + // Number of samples read int res; memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); idxin = 0; - res = getsample(&(inbuff[nin]), BLKIN - nin); + + // Read some samples + res = getsample(&(inbuff[nin]), BLKIN - nin); nin += res; + // If we haven't read any more samples, return how far we got if (nin < IQFilterLen * 2 + 2) - return (n); + return(n); } + // Process read samples into a brightness value iqfir(&inbuff[idxin], iqfilter, IQFilterLen, &I, &Q); ambuff[n] = pll(I, Q); idxin += 1; nin -= 1; } - return (n); + + return(n); } +// Get an entire row of pixels, without alignment int getpixelv(float *pvbuff, int nb) { + // Amplitude buffer static double ambuff[BLKAMP]; static int nam = 0; static int idxam = 0; @@ -160,7 +180,6 @@ int getpixelv(float *pvbuff, int nb) { double mult; mult = (double) Fi / Fe * FreqLine; - m = RSFilterLen / mult + 1; for (n = 0; n < nb; n++) { @@ -173,23 +192,24 @@ int getpixelv(float *pvbuff, int nb) { res = getamp(&(ambuff[nam]), BLKAMP - nam); nam += res; if (nam < m) - return (n); + return(n); } + // Denoise pvbuff[n] = rsfir(&(ambuff[idxam]), rsfilter, RSFilterLen, offset, mult) * mult * 256.0; - //printf("%g\n",pvbuff[n]); - shift = ((int) floor((RSMULT - offset) / mult)) + 1; - offset = shift * mult + offset - RSMULT ; + offset = shift * mult + offset - RSMULT; idxam += shift; nam -= shift; } - return (nb); + return(nb); } +// Align this line based off of the synchronisation markers int getpixelrow(float *pixelv) { + // Create an array for this row static float pixels[PixelLine + SyncFilterLen]; static int npv = 0; static int synced = 0; @@ -205,61 +225,70 @@ int getpixelrow(float *pixelv) { res = getpixelv(&(pixelv[npv]), SyncFilterLen + 2 - npv); npv += res; if (npv < SyncFilterLen + 2) - return (0); + return(0); } - /* test sync */ + // Test current synchronisation ecorr = fir(pixelv, Sync, SyncFilterLen); - corr = fir(&(pixelv[1]), Sync, SyncFilterLen); - lcorr = fir(&(pixelv[2]), Sync, SyncFilterLen); + corr = fir(&(pixelv[1]), Sync, SyncFilterLen - 1); + lcorr = fir(&(pixelv[2]), Sync, SyncFilterLen - 2); + // Calculate the per pixel offset FreqLine = 1.0+((ecorr-lcorr) / corr / PixelLine / 4.0); + // Maximum acceptable offset if (corr < 0.75 * max) { synced = 0; FreqLine = 1.0; } max = corr; + if (synced < 8) { - int shift, mshift; + int mshift; - if (npv < PixelLine + SyncFilterLen) { - res = getpixelv(&(pixelv[npv]), PixelLine + SyncFilterLen - npv); - npv += res; - if (npv < PixelLine + SyncFilterLen) - return (0); - } + if (npv < PixelLine + SyncFilterLen) { + res = getpixelv(&(pixelv[npv]), PixelLine + SyncFilterLen - npv); + npv += res; + if (npv < PixelLine + SyncFilterLen) + return(0); + } - /* lookup sync start */ - mshift = 0; - for (shift = 1; shift < PixelLine; shift++) { - double corr; + // Shift this line until we see the best results + mshift = 0; + for (int shift = 1; shift < PixelLine; shift++) { + double corr; - corr = fir(&(pixelv[shift + 1]), Sync, SyncFilterLen); - if (corr > max) { - mshift = shift; - max = corr; - } + corr = fir(&(pixelv[shift + 1]), Sync, SyncFilterLen); + if (corr > max) { + mshift = shift; + max = corr; + } + } + + // Shift memory, shifting this row + if (mshift != 0) { + memmove(pixelv, &(pixelv[mshift]), (npv - mshift) * sizeof(float)); + npv -= mshift; + synced = 0; + FreqLine = 1.0; + } else + synced += 1; } - if (mshift != 0) { - memmove(pixelv, &(pixelv[mshift]), (npv - mshift) * sizeof(float)); - npv -= mshift; - synced = 0; - FreqLine = 1.0; - } else - synced += 1; - } + // If there are not enough pixels try to grab some more if (npv < PixelLine) { res = getpixelv(&(pixelv[npv]), PixelLine - npv); npv += res; + // If we fail this then exit with 0 (which breaks the loop at main.c:338) if (npv < PixelLine) - return (0); + return(0); } + + // If we're finished reset npv to 0 if (npv == PixelLine) { npv = 0; - } else { + } else { // Move the pixel build buffer to the output buffer memmove(pixels, &(pixelv[PixelLine]), (npv - PixelLine) * sizeof(float)); npv -= PixelLine; } - return (1); + return(1); } diff --git a/fcolor.c b/fcolor.c index 10aafbb..e3608c8 100644 --- a/fcolor.c +++ b/fcolor.c @@ -104,18 +104,18 @@ void falsecolor(double v, double t, float *r, float *g, float *b) { if (t > fcinfo.Threshold) { if (v < fcinfo.Seathreshold) { - /* sea */ + // Sea top = fcinfo.SeaTop, bot = fcinfo.SeaBot; scv = v / fcinfo.Seathreshold; sct = (256.0 - t) / (256.0 - fcinfo.Threshold); } else { - /* ground */ + // Ground top = fcinfo.GroundTop, bot = fcinfo.GroundBot; scv = (v - fcinfo.Seathreshold) / (fcinfo.Landthreshold - fcinfo.Seathreshold); sct = (256.0 - t) / (256.0 - fcinfo.Threshold); } } else { - /* clouds */ + // Clouds top = fcinfo.CloudTop, bot = fcinfo.CloudBot; scv = v / 256.0; sct = (256.0 - t) / 256.0; @@ -129,27 +129,21 @@ void falsecolor(double v, double t, float *r, float *g, float *b) { }; void Ngvi(float **prow, int nrow) { - int n; - printf("GVI... "); fflush(stdout); - for (n = 0; n < nrow; n++) { + for (int n = 0; n < nrow; n++) { float *pixelv; - int i; pixelv = prow[n]; - for (i = 0; i < CH_WIDTH; i++) { + for (int i = 0; i < CH_WIDTH; i++) { float pv; double gvi; gvi = (pixelv[i + CHA_OFFSET] - pixelv[i + CHB_OFFSET]) / (pixelv[i + CHA_OFFSET] + pixelv[i + CHB_OFFSET]); pv = (gvi + 0.1) * 340.0; - if (pv > 255.0) - pv = 255.0; - if (pv < 0.0) - pv = 0.0; + pv = CLIP(pv, 0, 255); pixelv[i + CHB_OFFSET] = pv; } diff --git a/filter.c b/filter.c index b821533..b655cda 100755 --- a/filter.c +++ b/filter.c @@ -1,5 +1,5 @@ /* - * Aptec + * Aptdec * Copyright (c) 2004 by Thierry Leconte (F4DWV) * * $Id$ @@ -22,37 +22,39 @@ #include "filter.h" #include +// Sum of a matrix multiplication of 2 arrays float fir(float *buff, const float *coeff, const int len) { - int i; double r; r = 0.0; - for (i = 0; i < len; i++) { + for (int i = 0; i < len; i++) { r += buff[i] * coeff[i]; } - return r; + return(r); } +// Create an IQ sample from a sample buffer void iqfir(float *buff, const float *coeff, const int len, double *I, double *Q) { - int k; double i, q; i = q = 0.0; - for (k = 0; k < len; k++) { + for (int k = 0; k < len; k++) { q += buff[2*k] * coeff[k]; + // Average out the I samples, which gives us the DC offset i += buff[2*k]; } - i= buff[len-1] - i / len; - *I=i, *Q=q; + // Grab the peak value of the wave and subtract the DC offset + i = buff[len-1] - (i / len); + *I = i, *Q = q; } +// Denoise, I don't know how it works, but it does float rsfir(double *buff, const float *coeff, const int len, const double offset, const double delta) { - int i; - double n; double out; out = 0.0; - for (i = 0, n = offset; i < (len-1)/delta-1; n += delta, i++) { + double n = offset; + for (int i = 0; i < (len-1)/delta-1; n += delta, i++) { int k; double alpha; @@ -60,6 +62,5 @@ float rsfir(double *buff, const float *coeff, const int len, const double offset alpha = n - k; out += buff[i] * (coeff[k] * (1.0 - alpha) + coeff[k + 1] * alpha); } - return out; + return(out); } - diff --git a/filter.h b/filter.h index 4a3ecb9..8d94d5d 100755 --- a/filter.h +++ b/filter.h @@ -1,5 +1,5 @@ /* - * Aptec + * Aptdec * Copyright (c) 2003 by Thierry Leconte (F4DWV) * * $Id$ diff --git a/filtercoeff.h b/filtercoeff.h index dfed057..d76a573 100755 --- a/filtercoeff.h +++ b/filtercoeff.h @@ -1,5 +1,5 @@ /* - * Aptec + * Aptdec * Copyright (c) 2003 by Thierry Leconte (F4DWV) * * $Id$ @@ -26,8 +26,9 @@ const float iqfilter[IQFilterLen] = { 0.0205361, 0.0219524, 0.0235785, 0.0254648 -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 }; - +// Length of Sync #define SyncFilterLen 32 +// Pattern of a NOAA sync line const float Sync[SyncFilterLen] = { -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 }; @@ -95,6 +96,4 @@ const float rsfilter[RSFilterLen] = { -3.37279e-04, -8.80292e-06, -3.96418e-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 }; - - +-3.37279e-04 }; \ No newline at end of file diff --git a/gvipalette.h b/gvipalette.h index 8efaed1..a9ab41f 100644 --- a/gvipalette.h +++ b/gvipalette.h @@ -1,4 +1,4 @@ - unsigned char GviPalette[256*3] = { +unsigned char GviPalette[256*3] = { "\230t\17\233x\22\236{\27\241\200\33\244\203\37\247\210#\252\214'\255\220" ",\260\2240\264\2305\267\2358\272\240=\274\245A\300\251E\303\255I\306\262" "M\311\266Q\314\272V\317\276Z\322\302^\325\306b\330\312g\334\317k\337\323" diff --git a/image.c b/image.c index 34c5323..744ef76 100644 --- a/image.c +++ b/image.c @@ -1,5 +1,5 @@ /* - * Aptec + * Aptdec * Copyright (c) 2004 by Thierry Leconte (F4DWV) * * $Id$ @@ -26,6 +26,7 @@ #include #include "offsets.h" +#include "messages.h" #define REGORDER 3 typedef struct { @@ -33,7 +34,7 @@ typedef struct { } rgparam; static void rgcomp(double x[16], rgparam * rgpr) { - /*{ 0.106,0.215,0.324,0.433,0.542,0.652,0.78,0.87 ,0.0 }; */ + //const double y[9] = { 0.106, 0.215, 0.324, 0.433, 0.542, 0.652, 0.78, 0.87, 0.0 }; const double y[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; extern void polyreg(const int m, const int n, const double x[], const double y[], double c[]); @@ -48,15 +49,53 @@ static double rgcal(float x, rgparam * rgpr) { y += rgpr->cf[i] * p; p = p * x; } - return (y); + return(y); } - static double tele[16]; static double Cs; static int nbtele; -int Calibrate(float **prow, int nrow, int offset) { +// Contrast enchance +void equalise(float **prow, int nrow, int offset, int telestart, rgparam regr[30]){ + for (int n = 0; n < nrow; n++) { + float *pixelv; + int i; + + pixelv = prow[n]; + for (i = 0; i < CH_WIDTH; i++) { + float pv; + int k, kof; + + pv = pixelv[i + offset]; + + k = (n - telestart) / 128; + if (k >= nbtele) + k = nbtele - 1; + kof = (n - telestart) % 128; + + if (kof < 64) { + if (k < 1) { + pv = rgcal(pv, &(regr[k])); + } else { + pv = rgcal(pv, &(regr[k])) * (64 + kof) / 128.0 + rgcal(pv, &(regr[k - 1])) * (64 - kof) / 128.0; + } + } else { + if ((k + 1) >= nbtele) { + pv = rgcal(pv, &(regr[k])); + } else { + pv = rgcal(pv, &(regr[k])) * (192 - kof) / 128.0 + rgcal(pv, &(regr[k + 1])) * (kof - 64) / 128.0; + } + } + + pv = CLIP(pv, 0, 255); + pixelv[i + offset] = pv; + } + } +} + +// Get telemetry data for thermal calibration +int calibrate(float **prow, int nrow, int offset, int contrastBoost) { double teleline[3000]; double wedge[16]; rgparam regr[30]; @@ -66,58 +105,77 @@ int Calibrate(float **prow, int nrow, int offset) { int channel = -1; float max; - printf("Calibration... "); - fflush(stdout); - - /* build telemetry values lines */ + // Extract telemetry data, for a single pixel row for (n = 0; n < nrow; n++) { int i; teleline[n] = 0.0; + // Average the 40 center pixels from the telemetry block for (i = 3; i < 43; i++) { teleline[n] += prow[n][i + offset + CH_WIDTH]; } + // Compute the average teleline[n] /= 40.0; } + // A good minimum amount of pixels to find the telemetry start if (nrow < 192) { - fprintf(stderr, " not possible, not enough rows!\n"); - return (0); + fprintf(stderr, ERR_TELE_ROW); + return(0); } - /* find telemetry start in the 2nd third */ + // Find the biggest contrast in the telemetry max = 0.0; mtelestart = 0; + + // Only check the center of the image, where the signal is most likely strongest for (n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) { float df; + // Calculate the contrast, in other words + // (sum 4px below) / (sum 4px above) df = (teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + teleline[n - 1]) / (teleline[n] + teleline[n + 1] + teleline[n + 2] + teleline[n + 3]); + // Find the biggest contrast if (df > max) { mtelestart = n; max = df; } } - mtelestart -= 64; - telestart = mtelestart % 128; + // Calculate relative offset + telestart = (mtelestart - 64) % 128; + // If we cannot find the start of the telemetry or if there is not enough of it if (mtelestart < 0 || nrow < telestart + 128) { - fprintf(stderr, " impossible, not enough row\n"); - return (0); + fprintf(stderr, ERR_TELE_ROW); + return(0); } - /* compute wedges and regression */ + /* Compute wedges and regression + * + * This loop loops every 128 pixels after the relative telemetry start. + * Gets the values of where the wedges should be and then feeds into a + * regression algorithm which calculates the amount of noise on the + * telemetry. + * + * It then finds the part of the telemetry data with the least noise and + * turns it into digital values. + */ for (n = telestart, k = 0; n < nrow - 128; n += 128, k++) { int j; + // Loop through the 16 wedges 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]; + // Center 6 pixels + for (i = 1; i < 7; i++){ + wedge[j] += teleline[(j * 8) + n + i]; + } + // Average wedge[j] /= 6; } @@ -126,12 +184,12 @@ int Calibrate(float **prow, int nrow, int offset) { if (k == nrow / 256) { int i, l; - /* telemetry calibration */ + // Telemetry calibration for (j = 0; j < 16; j++) { tele[j] = rgcal(wedge[j], &(regr[k])); } - /* channel ID */ + // Channel ID for (j = 0, max = 10000.0, channel = -1; j < 6; j++) { float df; @@ -162,49 +220,13 @@ int Calibrate(float **prow, int nrow, int offset) { } nbtele = k; - /* calibrate */ - for (n = 0; n < nrow; n++) { - float *pixelv; - int i; - - pixelv = prow[n]; - for (i = 0; i < CH_WIDTH; i++) { - float pv; - int k, kof; - - pv = pixelv[i + offset]; - - k = (n - telestart) / 128; - if (k >= nbtele) - k = nbtele - 1; - kof = (n - telestart) % 128; - - if (kof < 64) { - if (k < 1) { - pv = rgcal(pv, &(regr[k])); - } else { - pv = rgcal(pv, &(regr[k])) * (64 + kof) / 128.0 + rgcal(pv, &(regr[k - 1])) * (64 - kof) / 128.0; - } - } else { - if ((k + 1) >= nbtele) { - pv = rgcal(pv, &(regr[k])); - } else { - pv = rgcal(pv, &(regr[k])) * (192 - kof) / 128.0 + rgcal(pv, &(regr[k + 1])) * (kof - 64) / 128.0; - } - } + // Image contrast + if(contrastBoost) equalise(prow, nrow, offset, telestart, regr); - if (pv > 255.0) - pv = 255.0; - if (pv < 0.0) - pv = 0.0; - pixelv[i + offset] = pv; - } - } - printf("Done\n"); - return (channel + 1); + return(channel + 1); } -/* ------------------------------temperature calibration -----------------------*/ +// --- Temperature Calibration --- // extern int satnum; #include "satcal.h" @@ -241,7 +263,7 @@ static void tempcomp(double t[16], int ch, tempparam * tpr) { /* compute radiance Black body */ C = satcal[satnum].rad[tpr->ch].vc; - tpr->Nbb = c1 * C * C * C / (exp(c2 * C / Tbb) - 1.0); + tpr->Nbb = c1 * C * C * C / (expm1(c2 * C / Tbb)); /* store Count Blackbody and space */ tpr->Cs = Cs * 4.0; @@ -261,13 +283,13 @@ static double tempcal(float Ce, tempparam * rgpr) { Ne = Nl + Nc; vc = satcal[satnum].rad[rgpr->ch].vc; - T = c2 * vc / log(c1 * vc * vc * vc / Ne + 1.0); + T = c2 * vc / log1p(c1 * vc * vc * vc / Ne); T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B; /* rescale to range 0-255 for -60 +40 'C */ T = (T - 273.15 + 60.0) / 100.0 * 256.0; - return (T); + return(T); } void Temperature(float **prow, int nrow, int channel, int offset) { @@ -289,10 +311,7 @@ void Temperature(float **prow, int nrow, int channel, int offset) { pv = tempcal(pixelv[i + offset], &temp); - if (pv > 255.0) - pv = 255.0; - if (pv < 0.0) - pv = 0.0; + pv = CLIP(pv, 0, 255); pixelv[i + offset] = pv; } } diff --git a/main.c b/main.c index f390f1f..169c6fe 100644 --- a/main.c +++ b/main.c @@ -25,16 +25,12 @@ #include #include -#ifdef WIN32 -#include "w32util.h" -#else #include -#endif #include #include -#include "version.h" +#include "messages.h" #include "offsets.h" #include "temppalette.h" @@ -49,70 +45,74 @@ static int initsnd(char *filename) { SF_INFO infwav; int res; - /* open wav input file */ + // Open audio file infwav.format = 0; inwav = sf_open(filename, SFM_READ, &infwav); if (inwav == NULL) { - fprintf(stderr, "Could not open %s for reading\n", filename); - return (1); + fprintf(stderr, ERR_FILE_READ, filename); + return(0); } res = init_dsp(infwav.samplerate); + printf("Input file: %s\n", filename); if(res < 0) { - fprintf(stderr, "Sample rate too low: %d\n", infwav.samplerate); - return (1); + fprintf(stderr, "Input sample rate too low: %d\n", infwav.samplerate); + return(0); + }else if(res > 0) { + fprintf(stderr, "Input sample rate too high: %d\n", infwav.samplerate); + return(0); } - if(res > 0) { - fprintf(stderr, "Sample rate too hight: %d\n", infwav.samplerate); - return (1); - } - fprintf(stderr, "Sample rate: %d\n", infwav.samplerate); + printf("Input sample rate: %d\n", infwav.samplerate); if (infwav.channels != 1) { fprintf(stderr, "Too many channels in input file: %d\n", infwav.channels); - return (1); + return(0); } - return (0); + return(1); } +// Get a sample from the wave file int getsample(float *sample, int nb) { - return (sf_read_float(inwav, sample, nb)); + return(sf_read_float(inwav, sample, nb)); } static png_text text_ptr[] = { - {PNG_TEXT_COMPRESSION_NONE, "Software", version, sizeof(version)}, + {PNG_TEXT_COMPRESSION_NONE, "Software", VERSION}, {PNG_TEXT_COMPRESSION_NONE, "Channel", NULL, 0}, - {PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA POES satellite image", 25} + {PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20} }; -static int ImageOut(char *filename, char *chid, float **prow, int nrow, int width, int offset, png_color *palette) { +static int ImageOut(char *filename, char *chid, float **prow, int nrow, int width, int offset, png_color *palette, int croptele, int layered) { FILE *pngfile; png_infop info_ptr; png_structp png_ptr; - int n; - /* init png lib */ + // Reduce the width of the image to componsate for the missing telemetry + if(croptele) width -= TOTAL_TELE; + + // Initalise the PNG writer png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); if (!png_ptr) { - fprintf(stderr, "Could not create a PNG write struct\n"); - return (1); + fprintf(stderr, ERR_PNG_WRITE); + return(1); } + // Metadata info_ptr = png_create_info_struct(png_ptr); if (!info_ptr) { png_destroy_write_struct(&png_ptr, (png_infopp) NULL); - fprintf(stderr, "Could not create a PNG info struct\n"); - return (1); + fprintf(stderr, ERR_PNG_INFO); + return(1); } if(palette == NULL) { - /* greyscale */ + // Greyscale image png_set_IHDR(png_ptr, info_ptr, width, nrow, 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); } else { - /* palette color mage */ + // Palleted color image png_set_IHDR(png_ptr, info_ptr, width, nrow, 8, PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); @@ -124,57 +124,77 @@ static int ImageOut(char *filename, char *chid, float **prow, int nrow, int widt png_set_text(png_ptr, info_ptr, text_ptr, 3); png_set_pHYs(png_ptr, info_ptr, 4000, 4000, PNG_RESOLUTION_METER); - printf("Writing %s... ", filename); + printf("Writing %s ", filename); fflush(stdout); pngfile = fopen(filename, "wb"); if (pngfile == NULL) { - fprintf(stderr, "Could not open %s for writing\n", filename); - return (1); + fprintf(stderr, ERR_FILE_WRITE, filename); + return(1); } png_init_io(png_ptr, pngfile); png_write_info(png_ptr, info_ptr); - for (n = 0; n < nrow; n++) { + for (int n = 0; n < nrow; n++) { float *pixelv; png_byte pixel[2*IMG_WIDTH]; - int i; pixelv = prow[n]; - for (i = 0; i < width; i++) { - pixel[i] = pixelv[i + offset]; + int f = 0; + for (int i = 0; i < width; i++) { + // Skip parts of the image that are telemtry + if(croptele){ + switch (i) { + case 0: + f += SYNC_WIDTH + SPC_WIDTH; + break; + case CH_WIDTH: + f += TELE_WIDTH + SYNC_WIDTH + SPC_WIDTH; + break; + case CH_WIDTH*2: + f += TELE_WIDTH; + } + } + + if(layered){ + // Layered image, basically overlay highlights in channel B over channel A + float cloud = CLIP(pixelv[i+CHB_OFFSET]-141, 0, 255)/114*255; + pixel[i] = CLIP(pixelv[i+CHA_OFFSET] + cloud, 0, 255); + }else{ + pixel[i] = pixelv[i + offset + f]; + } } png_write_row(png_ptr, pixel); + } png_write_end(png_ptr, info_ptr); fclose(pngfile); - printf("Done\n"); + printf("\nDone\n"); png_destroy_write_struct(&png_ptr, &info_ptr); - return (0); + return(0); } static int ImageRGBOut(char *filename, float **prow, int nrow) { FILE *pngfile; png_infop info_ptr; png_structp png_ptr; - int n; extern void falsecolor(double v, double t, float *r, float *g, float *b); - /* init png lib */ + // Initalise the PNG writer png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); if (!png_ptr) { - fprintf(stderr, "Could not create a PNG write struct\n"); - return (1); + fprintf(stderr, ERR_PNG_WRITE); + return(1); } info_ptr = png_create_info_struct(png_ptr); if (!info_ptr) { png_destroy_write_struct(&png_ptr, (png_infopp) NULL); - fprintf(stderr, "Could not create a PNG info struct\n"); - return (1); + fprintf(stderr, ERR_PNG_WRITE); + return(1); } - png_set_IHDR(png_ptr, info_ptr, CH_WIDTH , nrow , + png_set_IHDR(png_ptr, info_ptr, CH_WIDTH, nrow, 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); @@ -188,20 +208,19 @@ static int ImageRGBOut(char *filename, float **prow, int nrow) { fflush(stdout); pngfile = fopen(filename, "wb"); if (pngfile == NULL) { - fprintf(stderr, "Could not open %s for writing\n", filename); - return (1); + fprintf(stderr, ERR_FILE_WRITE, filename); + return(1); } png_init_io(png_ptr, pngfile); png_write_info(png_ptr, info_ptr); - for (n = 0; n < nrow ; n++) { + for (int n = 0; n < nrow ; n++) { png_color pix[CH_WIDTH]; - float *pixelc; - int i; + float *pixelc; pixelc = prow[n]; - for (i = 0; i < CH_WIDTH - 1; i++) { + for (int i = 0; i < CH_WIDTH - 1; i++) { float v, t; float r, g, b; @@ -220,11 +239,11 @@ static int ImageRGBOut(char *filename, float **prow, int nrow) { fclose(pngfile); printf("Done\n"); png_destroy_write_struct(&png_ptr, &info_ptr); - return (0); + return(0); } - -static void Distrib(char *filename,float **prow,int nrow) { +// Distribution of values between Channel A and Channel B +static void Distrib(char *filename, float **prow, int nrow) { unsigned int distrib[256][256]; int n; int x, y; @@ -253,152 +272,229 @@ static void Distrib(char *filename,float **prow,int nrow) { fprintf(df,"P2\n#max %d\n",max); fprintf(df,"256 256\n255\n"); - for(y=0;y<256;y++) - for(x=0;x<256;x++) + for(y = 0; y < 256; y++) + for(x = 0; x < 256; x++) fprintf(df, "%d\n", (int)((255.0 * (double)(distrib[y][x])) / (double)max)); fclose(df); } -extern int Calibrate(float **prow, int nrow, int offset); +extern int calibrate(float **prow, int nrow, int offset, int contrastBoost); extern void Temperature(float **prow, int nrow, int ch, int offset); extern int Ngvi(float **prow, int nrow); extern void readfconf(char *file); -extern int optind, opterr; +extern int optind; extern char *optarg; + +// Default to NOAA 19 int satnum = 4; static void usage(void) { - fprintf(stderr, "Aptdec [options] recordings ...\n"); - fprintf(stderr, "Options:\n -i [r|a|b|c|t] Output image type\n r: Raw\n a: A channel\n b: B channel\n c: False color\n t: Temperature\n -d Image destination directory.\n -s [15|16|17|18|19] Satellite number\n -c False color config file\n"); + printf("Aptdec [options] audio files ...\n" + "Options:\n" + " -e [c|t] Enhancements\n" + " c: Contrast enhance\n" + " t: Crop telemetry\n" + " -i [r|a|b|c|t] Output image type\n" + " r: Raw\n" + " a: A channel\n" + " b: B channel\n" + " c: False color\n" + " t: Temperature\n" + " l: Layered\n" + " -d Image destination directory.\n" + " -s [15-19] Satellite number\n" + " -c False color config file\n"); exit(1); } int main(int argc, char **argv) { char pngfilename[1024]; - char name[1024]; - char pngdirname[1024] = ""; - char imgopt[20] = "ac"; + char name[500]; + char pngdirname[500] = ""; + + // Images to create, default to a channel A and channel B image with contrast enhancement and cropped telemetry + char imgopt[20] = "ab"; + char enchancements[20] = "ct"; + + // Image buffer float *prow[3000]; - char *chid[] = { "?", "1", "2", "3A", "4", "5", "3B" }; - int nrow; + int nrow; + + // Mapping between wedge value and channel ID + char *chid[] = { "?", "1", "2", "3A", "4", "5", "3B" }; + char *chname[] = { "unknown", "visble", "near-infrared", "near-infrared", "thermal-infrared", "thermal-infrared", "thermal-infrared" }; + + // The active sensor in each channel int chA, chB; - int c; - printf("%s\n", version); + // Print version + printf("%s\n", VERSION); - opterr = 0; - if(argc == 1){ + // Print usage if there are no arguments + if(argc == 1) usage(); - } - while ((c = getopt(argc, argv, "c:d:i:s:")) != EOF) { + + int c; + while ((c = getopt(argc, argv, "c:d:i:s:e:")) != EOF) { switch (c) { + // Output directory name case 'd': strcpy(pngdirname, optarg); break; + // False color config file case 'c': readfconf(optarg); break; + // Output image type case 'i': strcpy(imgopt, optarg); break; + // Satellite number (for calibration) case 's': satnum = atoi(optarg)-15; + // Check if it's within the valid range if (satnum < 0 || satnum > 4) { - fprintf(stderr, "Invalid satellite number, must be in the range [15-19]\n"); + fprintf(stderr, "Invalid satellite number, it must be the range [15-19]\n"); exit(1); } break; + // Enchancements + case 'e': + strcpy(enchancements, optarg); + break; default: usage(); } } - for (nrow = 0; nrow < 3000; nrow++) - prow[nrow] = NULL; + if(optind == argc){ + printf("No audio files provided.\n"); + usage(); + } + + // Pull this in from filtercoeff.h + extern float Sync[32]; + // Loop through the provided files for (; optind < argc; optind++) { chA = chB = 0; + // Generate output name strcpy(pngfilename, argv[optind]); strcpy(name, basename(pngfilename)); strtok(name, "."); - if (pngdirname[0] == '\0') { - strcpy(pngfilename, argv[optind]); + + if (pngdirname[0] == '\0') strcpy(pngdirname, dirname(pngfilename)); - } - /* open snd input */ - if (initsnd(argv[optind])) + // Open sound file, exit if that fails + if (initsnd(argv[optind]) == 0) exit(1); - /* main loop */ - printf("Decoding: %s \n", argv[optind]); + // Main image building loop for (nrow = 0; nrow < 3000; nrow++) { - if (prow[nrow] == NULL) + // Allocate 2150 floats worth of memory for every line of the image prow[nrow] = (float *) malloc(sizeof(float) * 2150); + + // Read into prow and break the loop once we reach the end of the image if (getpixelrow(prow[nrow]) == 0) break; + + // Signal level + double signal = 0; + for (int i = 0; i < 32; i++) signal += prow[nrow][i] - Sync[i]; - printf("%d\r", nrow); + // Signal level bar + char bar[30]; + for (int i = 0; i < 30; i++) { + bar[i] = ' '; + if(i < signal/2000*30) bar[i] = '#'; + } + + printf("Row: %d, Signal Strength: %s\r", nrow, bar); fflush(stdout); } - printf("\nDone\n"); + printf("\nTotal rows: %d\n", nrow); + sf_close(inwav); - /* raw image */ - if (strchr(imgopt, (int) 'r') != NULL) { + // Layered images require contrast enhancements + int contrastBoost = CONTAINS(enchancements, 'c') || CONTAINS(imgopt, 'l'); + + chA = calibrate(prow, nrow, CHA_OFFSET, 0); + chB = calibrate(prow, nrow, CHB_OFFSET, 0); + printf("Channel A: %s (%s)\n", chid[chA], chname[chA]); + printf("Channel B: %s (%s)\n", chid[chB], chname[chB]); + + // Temperature + if (CONTAINS(imgopt, 't') && chB >= 4) { + Temperature(prow, nrow, chB, CHB_OFFSET); + sprintf(pngfilename, "%s/%s-t.png", pngdirname, name); + ImageOut(pngfilename, "Temperature", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)TempPalette, 0, 0); + } + + // We have to run the contrast enhance here because the temperature function requires real data + // Yes, this is bodgy, yes I should replace this with something more elegant + chA = calibrate(prow, nrow, CHA_OFFSET, contrastBoost); + chB = calibrate(prow, nrow, CHB_OFFSET, contrastBoost); + + // Layered + if (CONTAINS(imgopt, 'l')){ + sprintf(pngfilename, "%s/%s-l.png", pngdirname, name); + ImageOut(pngfilename, "Layered", prow, nrow, CH_WIDTH, 0, NULL, 0, 1); + } + + // Raw image + if (CONTAINS(imgopt, 'r')) { + int croptele = CONTAINS(enchancements, 't'); + char channelstr[45]; + sprintf(channelstr, "%s (%s) & %s (%s)", chid[chA], chname[chA], chid[chB], chname[chB]); + sprintf(pngfilename, "%s/%s-r.png", pngdirname, name); - ImageOut(pngfilename, "raw", prow, nrow, IMG_WIDTH, 0, NULL); + ImageOut(pngfilename, channelstr, prow, nrow, IMG_WIDTH, 0, NULL, croptele, 0); } - /* Channel A */ - if (((strchr(imgopt, (int) 'a') != NULL) - || (strchr(imgopt, (int) 'c') != NULL) - || (strchr(imgopt, (int) 'd') != NULL))) { - chA = Calibrate(prow, nrow, CHA_OFFSET); - if (chA >= 0 && strchr(imgopt, (int) 'a') != NULL) { - sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chA]); - ImageOut(pngfilename, chid[chA], prow, nrow, CH_WIDTH , CHA_OFFSET, NULL); - } + // Channel A + if (CONTAINS(imgopt, 'a')) { + char channelstr[21]; + sprintf(channelstr, "%s (%s)", chid[chA], chname[chA]); + + sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chA]); + ImageOut(pngfilename, channelstr, prow, nrow, CH_WIDTH, CHA_OFFSET, NULL, 0, 0); } - /* Channel B */ - if ((strchr(imgopt, (int) 'b') != NULL) - || (strchr(imgopt, (int) 'c') != NULL) - || (strchr(imgopt, (int) 't') != NULL) - || (strchr(imgopt, (int) 'd') != NULL)) { - chB = Calibrate(prow, nrow, CHB_OFFSET); - if (chB >= 0 && strchr(imgopt, (int) 'b') != NULL) { - sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chB]); - ImageOut(pngfilename, chid[chB], prow, nrow, CH_WIDTH , CHB_OFFSET ,NULL); - } - if (chB > 3) { - Temperature(prow, nrow, chB, CHB_OFFSET); - if (strchr(imgopt, (int) 't') != NULL) { - sprintf(pngfilename, "%s/%s-t.png", pngdirname, name); - ImageOut(pngfilename, "Temperature", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)TempPalette); - } - } + // Channel B + if (CONTAINS(imgopt, 'b')) { + char channelstr[21]; + sprintf(channelstr, "%s (%s)", chid[chB], chname[chB]); + + sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chB]); + ImageOut(pngfilename, channelstr, prow, nrow, CH_WIDTH , CHB_OFFSET, NULL, 0, 0); } - /* distribution */ - if (chA && chB && strchr(imgopt, (int) 'd') != NULL) { + // Distribution map + if (CONTAINS(imgopt, 'd')) { sprintf(pngfilename, "%s/%s-d.pnm", pngdirname, name); Distrib(pngfilename, prow, nrow); } - /* color image */ - if (chA == 2 && chB == 4 && strchr(imgopt, (int) 'c') != NULL) { - sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); - ImageRGBOut(pngfilename, prow, nrow); + // False color image + if(CONTAINS(imgopt, 'c')){ + if (chA == 2 && chB == 4) { // Normal false color + sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); + ImageRGBOut(pngfilename, prow, nrow); + } else if (chA == 1 && chB == 2) { // GVI false color + Ngvi(prow, nrow); + sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); + ImageOut(pngfilename, "GVI False Color", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)GviPalette, 0, 0); + } else { + fprintf(stderr, "Lacking channels required for false color computation.\n"); + } } - /* GVI image */ - if (chA == 1 && chB == 2 && strchr(imgopt, (int) 'c') != NULL) { - Ngvi(prow, nrow); - sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); - ImageOut(pngfilename, "GVI", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)GviPalette); - } - } + // Free the allocated memory + for(int i = 0; i < 3000; i++) free(prow[i]); + } + exit(0); } diff --git a/messages.h b/messages.h new file mode 100644 index 0000000..702ade2 --- /dev/null +++ b/messages.h @@ -0,0 +1,6 @@ +#define ERR_FILE_WRITE "Could not open %s for writing\n" +#define ERR_FILE_READ "Could not open %s for reading\n" +#define ERR_PNG_WRITE "Could not create a PNG write struct\n" +#define ERR_PNG_INFO "Could not create a PNG info struct\n" +#define ERR_TELE_ROW "Telemetry decoding error, not enough rows.\n" +#define VERSION "Aptdec CVS version (c) 2004-2005 Thierry Leconte F4DWV" \ No newline at end of file diff --git a/offsets.h b/offsets.h index f034cb9..d8980bb 100644 --- a/offsets.h +++ b/offsets.h @@ -6,3 +6,7 @@ #define IMG_WIDTH 2080 #define CHA_OFFSET (SYNC_WIDTH+SPC_WIDTH) #define CHB_OFFSET (SYNC_WIDTH+SPC_WIDTH+CH_WIDTH+TELE_WIDTH+SYNC_WIDTH+SPC_WIDTH) +#define TOTAL_TELE (SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH+SYNC_WIDTH+SPC_WIDTH+TELE_WIDTH) + +#define CLIP(val, bottom, top) (val > top ? top : (val > bottom ? val : bottom)) +#define CONTAINS(str, char) (strchr(str, (int) char) != NULL) \ No newline at end of file diff --git a/reg.c b/reg.c index 508d96c..316856d 100644 --- a/reg.c +++ b/reg.c @@ -12,10 +12,11 @@ #include -#define DMAX 5 /* Maximum degree of polynomial */ -#define NMAX 10 /* Maximum number of points */ +#define DMAX 5 /* Maximum degree of polynomial */ +#define NMAX 10 /* Maximum number of points */ static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]); + void polyreg(const int M, const int N, const double X[], const double Y[], double C[]) { int R, K, J; /* Loop counters */ double A[DMAX][DMAX]; /* A */ @@ -41,7 +42,8 @@ void polyreg(const int M, const int N, const double X[], const double Y[], doubl /* Zero the array */ for (J = 1; J <= 2 * M; J++) - P[J] = 0; + P[J] = 0; + P[0] = N; /* Compute the sum of powers of x_(K-1) */ @@ -74,9 +76,8 @@ static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) { double SUM, DET = 1.0; int T; - /* Initialize the pointer vector */ - for (J = 0; J < N; J++) + for (J = 0; J < NMAX; J++) Row[J] = J; /* Start LU factorization */ @@ -108,13 +109,15 @@ static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) { A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C]; } } - } /* End of L*U factorization routine */ + } + /* End of L*U factorization routine */ DET = DET * A[Row[N - 1]][N - 1]; /* Start the forward substitution */ for (K = 0; K < N; K++) Y[K] = B[K]; - Y[0] = B[Row[0]]; + + Y[0] = B[Row[0]]; for (K = 1; K < N; K++) { SUM = 0; diff --git a/satcal.h b/satcal.h index 643d173..ebe897c 100644 --- a/satcal.h +++ b/satcal.h @@ -6,103 +6,101 @@ const struct { struct { float Ns; float b[3]; - } cor[3]; -} satcal[] = -{ -/* calibration coeff from NOAA KLM POES satellite user guide */ -/* http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/tables.htm */ - -{/* NOAA-15 */ -{ /* PRT coeff d0,d1,d2 */ -{276.60157 , 0.051045 , 1.36328E-06}, -{276.62531 , 0.050909 , 1.47266E-06}, -{276.67413 , 0.050907 , 1.47656E-06}, -{276.59258 , 0.050966 , 1.47656E-06} -}, -{ /* channel radiance coeff vc,A,B*/ -{925.4075 , 0.337810 , 0.998719}, /* channel 4 */ -{839.8979 , 0.304558 , 0.999024}, /* channel 5 */ -{2695.9743 , 1.621256 , 0.998015} /* channel 3B */ -}, -{ /* nonlinear radiance correction Ns,b0,b1,b2 */ -{-4.50 , {4.76 , -0.0932 , 0.0004524}}, /* channel 4 */ -{-3.61 , {3.83 , -0.0659 , 0.0002811}}, /* channel 5*/ -{0.0,{0.0,0.0,0.0}} /* channel 3B*/ -} -}, -{/* NOAA-16 */ -{ /* PRT coeff d0,d1,d2 */ -{276.355 , 5.562E-02 ,-1.590E-05}, -{276.142 , 5.605E-02 ,-1.707E-05}, -{275.996 , 5.486E-02 ,-1.223E-05}, -{276.132 , 5.494E-02 ,-1.344E-05} -}, -{ /* channel radiance coeff vc,A,B*/ -{917.2289 ,0.332380 ,0.998522}, /* channel 4 */ -{838.1255 ,0.674623 ,0.998363}, /* channel 5 */ -{2700.1148 ,1.592459,0.998147} /* channel 3B */ -}, -{ /* nonlinear radiance correction Ns,b0,b1,b2 */ -{-2.467, {2.96,-0.05411,0.00024532}}, /* channel 4 */ -{-2.009,{2.25,-0.03665,0.00014854}}, /* channel 5*/ -{0.0,{0.0,0.0,0.0}} /* channel 3B*/ -} -}, -{/* NOAA 17 */ -{ /* PRT coeff d0,d1,d2 */ -{276.628 , 0.05098 , 1.371e-06}, -{276.538 , 0.05098 , 1.371e-06}, -{276.761 , 0.05097 , 1.369e-06}, -{276.660 , 0.05100 , 1.348e-06} -}, -{ /* channel radiance coeff vc,A,B*/ -{926.2947 , 0.271683 , 0.998794}, /* channel 4 */ -{839.8246 , 0.309180 , 0.999012}, /* channel 5 */ -{2669.3554 , 1.702380 , 0.997378} /* channel 3B */ -}, -{ /* nonlinear radiance correction Ns,b0,b1,b2 */ -{-8.55 , {8.22 , -0.15795 , 0.00075579}}, /* channel 4 */ -{-3.97 , {4.31 , -0.07318 , 0.00030976}}, /* channel 5 */ -{0.0,{0.0,0.0,0.0}} /* channel 3B*/ -} -}, -{/* NOAA 18 */ -{ /* PRT coeff d0,d1,d2 */ -{276.601 , 0.05090 , 1.657e-06}, -{276.683 , 0.05101 , 1.482e-06}, -{276.565 , 0.05117 , 1.313e-06}, -{276.615 , 0.05103 , 1.484e-06} -}, -{ /* channel radiance coeff vc,A,B*/ -{928.1460 , 0.436645 , 0.998607}, /* channel 4 */ -{833.2532 , 0.253179 , 0.999057}, /* channel 5 */ -{2659.7952 , 1.698704 , 0.996960} /* channel 3B */ + } cor[3]; +/* Calibration corficcients taken from the NOAA KLM satellite user guide + * https://web.archive.org/web/20141220021557/https://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/klm/tables.htm + */ +} satcal[] = { +{ // NOAA 15 + { // PRT coefficient d0, d1, d2 + {276.60157, 0.051045, 1.36328E-06}, + {276.62531, 0.050909, 1.47266E-06}, + {276.67413, 0.050907, 1.47656E-06}, + {276.59258, 0.050966, 1.47656E-06} + }, + { // Channel radiance coefficient vc, A, B + {925.4075, 0.337810, 0.998719}, // Channel 4 + {839.8979, 0.304558, 0.999024}, // Channel 5 + {2695.9743, 1.621256, 0.998015} // Channel 3B + }, + { // Nonlinear radiance correction Ns, b0, b1, b2 + {-4.50, {4.76, -0.0932, 0.0004524}}, // Channel 4 + {-3.61, {3.83, -0.0659, 0.0002811}}, // Channel 5 + {0.0, {0.0, 0.0, 0.0}} // Channel 3B + } }, -{ /* nonlinear radiance correction Ns,b0,b1,b2 */ -{-5.53 , {5.82 , -0.11069 , 0.00052337}}, /* channel 4 */ -{-2.22 , {2.67 , -0.04360 , 0.00017715}}, /* channel 5 */ -{0.0,{0.0,0.0,0.0}} /* channel 3B*/ -} +{ // NOAA 16 + { // PRT coeff d0, d1, d2 + {276.355, 5.562E-02, -1.590E-05}, + {276.142, 5.605E-02, -1.707E-05}, + {275.996, 5.486E-02, -1.223E-05}, + {276.132, 5.494E-02, -1.344E-05} + }, + { // Channel radiance coefficient vc, A, B + {917.2289, 0.332380, 0.998522}, // Channel 4 + {838.1255, 0.674623, 0.998363}, // Channel 5 + {2700.1148, 1.592459, 0.998147} // Channel 3B + }, + { // Nonlinear radiance correction Ns, b0, b1, b2 + {-2.467, {2.96, -0.05411, 0.00024532}}, // Channel 4 + {-2.009, {2.25, -0.03665, 0.00014854}}, // Channel 5 + {0.0, {0.0, 0.0, 0.0}} // Channel 3B + } }, -{/* NOAA 19 */ -{ /* PRT coeff d0,d1,d2 */ -{276.6067 , 0.051111 , 1.405783E-06}, -{276.6119 , 0.051090 , 1.496037E-06}, -{276.6311 , 0.051033 , 1.496990E-06}, -{276.6268 , 0.051058 , 1.493110E-06} +{ // NOAA 17 + { // PRT coefficient d0, d1, d2 + {276.628, 0.05098, 1.371e-06}, + {276.538, 0.05098, 1.371e-06}, + {276.761, 0.05097, 1.369e-06}, + {276.660, 0.05100, 1.348e-06} + }, + { // Channel radiance coefficient vc, A, B + {926.2947, 0.271683, 0.998794}, // Channel 4 + {839.8246, 0.309180, 0.999012}, // Channel 5 + {2669.3554, 1.702380, 0.997378} // Channel 3B + }, + { // Nonlinear radiance correction Ns, b0, b1, b2 + {-8.55, {8.22, -0.15795, 0.00075579}}, // Channel 4 + {-3.97, {4.31, -0.07318, 0.00030976}}, // Channel 5 + {0.0, {0.0, 0.0, 0.0}} // Channel 3B + } }, -{ /* channel radiance coeff vc,A,B*/ -{928.9 , 0.53959 , 0.998534}, /* channel 4 */ -{831.9 , 0.36064 , 0.998913}, /* channel 5 */ -{2670.0 , 1.67396 , 0.997364} /* channel 3B */ +{ // NOAA 18 + { // PRT coefficient d0, d1, d2 + {276.601, 0.05090, 1.657e-06}, + {276.683, 0.05101, 1.482e-06}, + {276.565, 0.05117, 1.313e-06}, + {276.615, 0.05103, 1.484e-06} + }, + { // Channel radiance coefficient vc, A, B + {928.1460, 0.436645, 0.998607}, // Channel 4 + {833.2532, 0.253179, 0.999057}, // Channel 5 + {2659.7952, 1.698704, 0.996960} // Channel 3B + }, + { // Nonlinear radiance correction Ns, b0, b1, b2 + {-5.53, {5.82, -0.11069, 0.00052337}}, // Channel 4 + {-2.22, {2.67, -0.04360, 0.00017715}}, // Channel 5 + {0.0, {0.0, 0.0, 0.0}} // Channel 3B + } }, -{ /* nonlinear radiance correction Ns,b0,b1,b2 */ -{-5.49 , {5.70 -0.11187 , 0.00054668}}, /* channel 4 */ -{-3.39 , {3.58 -0.05991 , 0.00024985}}, /* channel 5 */ -{0.0,{0.0,0.0,0.0}} /* channel 3B*/ -} -} -}; +{ // NOAA 19 + { // PRT coefficient d0, d1, d2 + {276.6067, 0.051111, 1.405783E-06}, + {276.6119, 0.051090, 1.496037E-06}, + {276.6311, 0.051033, 1.496990E-06}, + {276.6268, 0.051058, 1.493110E-06} + }, + { // Channel radiance coefficient vc, A, B + {928.9, 0.53959, 0.998534}, // Channel 4 + {831.9, 0.36064, 0.998913}, // Channel 5 + {2670.0, 1.67396, 0.997364} // Channel 3B + }, + { // Nonlinear radiance correction Ns, b0, b1, b2 + {-5.49, {5.70 -0.11187, 0.00054668}}, // Channel 4 + {-3.39, {3.58 -0.05991, 0.00024985}}, // Channel 5 + {0.0, {0.0, 0.0, 0.0}} // Channel 3B + } +}}; -const float c1=1.1910427e-5; -const float c2=1.4387752; \ No newline at end of file +const float c1 = 1.1910427e-5; +const float c2 = 1.4387752; \ No newline at end of file diff --git a/temppalette.h b/temppalette.h index 954a683..6e6d716 100644 --- a/temppalette.h +++ b/temppalette.h @@ -1,4 +1,4 @@ -unsigned char TempPalette[256*3] = { +unsigned char TempPalette[256*3] = { "\376\376\376\376\376\376\375\375\376\374\375\376\374\375\375\374\373\375" "\373\373\375\372\373\375\372\373\374\372\372\374\371\372\374\371\371\375" "\370\371\374\367\370\375\367\370\374\367\367\374\366\367\373\366\366\373" @@ -40,4 +40,3 @@ unsigned char TempPalette[256*3] = { "\2261\373\2161\373\206.\373}-\374u,\375k*\374b(\375Y'\375O%\375E#\375;\"" "\3762\40\376(\37\376\35\35", }; - diff --git a/textlogo.png b/textlogo.png new file mode 100644 index 0000000..95207f9 Binary files /dev/null and b/textlogo.png differ diff --git a/version.h b/version.h deleted file mode 100644 index 835ebdb..0000000 --- a/version.h +++ /dev/null @@ -1,2 +0,0 @@ -char version[] = "Aptec CVS version (c) 2004-2005 Thierry Leconte F4DWV"; - diff --git a/w32util.c b/w32util.c deleted file mode 100755 index 7c23419..0000000 --- a/w32util.c +++ /dev/null @@ -1,123 +0,0 @@ -#include -#include - -char * basename (const char *filename) { - const char *p = filename + strlen (filename); - while (p != filename) { - p--; - if (*p == '/' || *p == '\\') - return ((char *) p + 1); - } - return ((char *) filename); -} - -char *dirname(const char *path) { - static char bname[1024]; - register const char *endp; - - /* Empty or NULL string gets treated as "." */ - if (path == NULL || *path == '\0') { - (void)strncpy(bname, ".", sizeof bname); - return(bname); - } - - /* Strip trailing slashes */ - endp = path + strlen(path) - 1; - while (endp > path && *endp == '\\') - endp--; - - /* Find the start of the dir */ - while (endp > path && *endp != '\\') - endp--; - - /* Either the dir is "/" or there are no slashes */ - if (endp == path) { - (void)strncpy(bname, *endp == '\\' ? "\\": ".", sizeof bname); - return(bname); - } else { - do { - endp--; - } while (endp > path && *endp == '\\'); - } - - if (endp - path + 2 > sizeof(bname)) { - return(NULL); - } - strncpy(bname, path, endp - path + 2); - return(bname); -} - -int opterr = 1, /* if error message should be printed */ - optind = 1, /* index into parent argv vector */ - optopt, /* character checked for validity */ - optreset; /* reset getopt */ -char *optarg; /* argument associated with option */ - -#define BADCH (int)'?' -#define BADARG (int)':' -#define EMSG "" - -/* - * getopt -- - * Parse argc/argv argument vector. - */ -int getopt(nargc, nargv, ostr) - int nargc; - char * const *nargv; - const char *ostr; -{ - #define __progname nargv[0] - static char *place = EMSG; /* option letter processing */ - char *oli; /* option letter list index */ - - if (optreset || !*place) { /* update scanning pointer */ - optreset = 0; - if (optind >= nargc || *(place = nargv[optind]) != '-') { - place = EMSG; - return (-1); - } - if (place[1] && *++place == '-') { /* found "--" */ - ++optind; - place = EMSG; - return (-1); - } - } /* option letter okay? */ - if ((optopt = (int)*place++) == (int)':' || - !(oli = strchr(ostr, optopt))) { - /* - * if the user didn't specify '-' as an option, - * assume it means -1. - */ - if (optopt == (int)'-') - return (-1); - if (!*place) - ++optind; - if (opterr && *ostr != ':') - (void)printf("%s: illegal option -- %c\n", __progname, optopt); - return (BADCH); - } - if (*++oli != ':') { /* don't need argument */ - optarg = NULL; - if (!*place) - ++optind; - } - else { /* need an argument */ - if (*place) /* no white space */ - optarg = place; - else if (nargc <= ++optind) { /* no arg */ - place = EMSG; - if (*ostr == ':') - return (BADARG); - if (opterr) - (void)printf( - "%s: option requires an argument -- %c\n", - __progname, optopt); - return (BADCH); - } else /* white space */ - optarg = nargv[optind]; - place = EMSG; - ++optind; - } - return (optopt); /* dump back option letter */ -} - diff --git a/w32util.h b/w32util.h deleted file mode 100755 index c6cd387..0000000 --- a/w32util.h +++ /dev/null @@ -1,12 +0,0 @@ -extern char * basename(const char *filename); - -extern char *dirname(const char *path); - -extern int opterr, /* if error message should be printed */ - optind , /* index into parent argv vector */ - optopt, /* character checked for validity */ - optreset; /* reset getopt */ -extern char *optarg; /* argument associated with option */ - - -extern int getopt(int nargc, char * const *nargv, const char *ostr);