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

Large changes 

Add a denoise filter
Fix temperature calibration
Tidy up some files
Improve MCIR
tags/v1.8.0
Xerbo 4 years ago
parent
87f8304f57
commit
e990602455
11 changed files with 938 additions and 726 deletions
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. +7
    -5
      Makefile
  2. +10
    -8
      README.md
  3. +159
    -0
      color.c
  4. +13
    -11
      dsp.c
  5. +1
    -1
      falsecolor.conf
  6. +42
    -61
      fcolor.c
  7. +106
    -31
      image.c
  8. +231
    -517
      main.c
  9. +56
    -0
      median.c
  10. +0
    -92
      palette.h
  11. +313
    -0
      pngio.c

+ 7
- 5
Makefile View File

@@ -1,17 +1,19 @@
CC = gcc
BIN = /usr/bin
INCLUDES = -I.
CFLAGS = -O3 -Wall $(INCLUDES)
OBJS = main.o image.o dsp.o filter.o reg.o fcolor.o
CFLAGS = -O3 -DNDEBUG -Wall $(INCLUDES)
OBJS = main.o image.o dsp.o filter.o reg.o fcolor.o pngio.o median.o color.o

aptdec: $(OBJS)
$(CC) -o $@ $(OBJS) -lm -lsndfile -lpng

main.o: main.c palette.h offsets.h messages.h
dsp.o: dsp.c filtercoeff.h filter.h
color.o: color.c
main.o: main.c 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 messages.h
image.o: image.c offsets.h messages.h offsets.h
fcolor.o: fcolor.c offsets.h
pngio.o: pngio.c offsets.h messages.h

clean:
rm -f *.o aptdec


+ 10
- 8
README.md View File

@@ -7,7 +7,7 @@ Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 20
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 audio files into .png images.
Aptdec can convert these audio files into `png` images.

For each audio file up to 6 images can be generated:

@@ -53,11 +53,11 @@ To uninstall
```
-i [r|a|b|c|t|l]
Output image type
Raw (r), Channel A (a), Channel B (b), False Color (c), Temperature (t)
Raw (r), Channel A (a), Channel B (b), False Color (c), Temperature (t) or MCIR (m)
Default: "ab"

-d <dir>
Images destination directory (optional).
Images destination directory (optional)
Default: Current directory

-s [15|16|17|18|19]
@@ -66,15 +66,15 @@ For temperature calibration
Default: "19"

-e [t|h]
Enhancements
Histogram equalise (h), Crop Telemetry (t)
Effects
Histogram equalise (h), Crop Telemetry (t) or Denoise (d)
Defaults: off

-m <file>
Map overlay generated by wxmap
Map file generated by wxmap

-c <file>
Use configuration file for false color generation.
Use configuration file for false color generation
Default: Internal defaults
```

@@ -88,11 +88,13 @@ Image names are `audiofile-x.png`, where `x` is:
- Sensor ID (1, 2, 3A, 3B, 4, 5) for channel A|B images
- `c` for false color.
- `t` for temperature calibrated images
- `m` for MCIR images

Currently there are 3 available enhancements:
Currently there are 3 available effects:

- `t` for crop telemetry, off by default, only has effects on raw images
- `h` for histogram equalise, stretch the colors in the image to black and white
- `d` for a median denoise filter

## Example



+ 159
- 0
color.c View File

@@ -0,0 +1,159 @@
/*
* This file is part of Aptdec.
* Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2020
*
* Aptdec is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/

#include "offsets.h"

typedef struct {
float r, g, b;
} rgb_t;

rgb_t applyPalette(char *palette, int val){
return (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]
};
}

rgb_t RGBcomposite(rgb_t top, float top_a, rgb_t bottom, float bottom_a){
rgb_t composite;
composite.r = MCOMPOSITE(top.r, top_a, bottom.r, bottom_a);
composite.g = MCOMPOSITE(top.g, top_a, bottom.g, bottom_a);
composite.b = MCOMPOSITE(top.b, top_a, bottom.b, bottom_a);
return composite;
}

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"
"o\342\330s\344\333w\350\337{\352\344\200\356\350\203\361\354\207\364\360"
"\213\367\364\220\372\370\224\376\376\230\376\375\230\376\372\227\376\366"
"\224\376\362\222\376\357\220\376\352\216\376\347\213\376\343\211\376\340"
"\207\376\334\205\376\330\202\376\325\200\376\321\177\376\315|\376\312z\376"
"\306y\376\302v\376\277t\376\273q\376\267o\376\263m\376\260k\376\254h\376"
"\250f\376\245d\376\241b\376\235_\376\232^\376\226[\376\223Y\376\217W\376"
"\213U\376\207R\376\203Q\376\200N\376}L\376zJ\376vG\376rE\376nB\376k@\376"
"g>\376d<\376`:\376\\7\376X6\376T3\376Q0\376M/\376J-\376F*\376C(\376>%\376"
";$\3767!\3763\37\3760\35\376,\32\376(\31\376$\26\376!\24\376\35\22\376\32"
"\20\376\26\15\376\22\13\376\16\11\376\12\6\376\7\4\376\0\0\373\2\0\366\3"
"\0\362\5\0\355\7\0\350\11\0\343\13\0\336\15\0\332\17\0\325\21\0\320\22\0"
"\313\24\0\307\26\0\303\27\0\275\32\0\270\33\0\264\35\0\257\37\0\253!\0\246"
"#\0\241%\0\234&\0\227)\0\223+\0\216,\0\212.\0\2050\0\2002\0}4\0w6\0s7\0n"
"9\0i;\0e>\0`?\0[A\0VC\0QE\0MG\0HH\0CK\0?M\0:N\0""5P\0""0R\0,S\0'V\0\"X\0"
"\36Z\0\31\\\0\23^\0\17_\0\12a\0\3b\0\6e\2\12i\5\17n\11\23r\14\30v\17\34{"
"\23\40\200\26$\203\31)\207\35.\215\40""1\221$6\225':\232*?\236.C\2431G\250"
"4K\2548O\260;U\265>Y\272A\\\276Eb\303Ie\307Lj\313On\320Ss\324Vw\331Y{\336"
"]\200\342`\203\347d\207\353g\213\357j\217\364n\223\371q\230\376t\222\375"
"{\207\374\205}\373\216r\371\230g\371\241\\\370\254Q\367\265F\365\300;\364"
"\3110\363\323%\362\334\32\361\347\17\357\361\3\356\371\4\353\370\3\347\365"
"\4\342\361\4\335\357\4\331\354\4\323\351\3\317\346\3\313\343\4\306\340\3"
"\301\335\3\275\332\3\270\327\3\263\324\2\257\322\3\252\316\3\245\314\3\241"
"\311\2\234\306\2\230\303\2\223\300\3\216\275\3\211\273\2\205\267\2\200\265"
"\2|\262\2w\257\2s\254\2n\251\2j\246\2e\243\2`\240\2[\235\2W\232\1S\230\2"
"M\225\1I\221\2E\217\1@\214\1;\211\1""7\206\1""1\203\1-\200\0(~\1${\0\37y"
"\0\33u\0\25r\0\21p\0\14l\0\7j\0\3"
};
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"
"\365\366\373\364\366\373\364\365\374\363\365\373\363\364\373\363\364\373"
"\362\363\372\361\363\372\361\362\372\361\361\372\360\361\372\357\361\372"
"\357\361\372\356\360\372\355\357\371\355\357\371\355\356\372\354\356\371"
"\354\355\371\354\355\371\353\355\371\353\354\370\352\353\370\351\353\370"
"\351\353\370\350\352\370\350\351\370\347\351\370\347\351\370\346\350\367"
"\346\350\367\346\350\367\345\346\367\344\346\367\344\346\367\344\345\366"
"\342\346\366\343\344\366\342\345\366\341\344\366\341\343\365\340\343\366"
"\337\342\365\337\341\366\336\341\365\337\341\365\336\340\365\335\340\365"
"\335\337\364\334\337\364\334\336\364\334\336\364\333\336\364\333\335\364"
"\332\334\364\331\334\363\331\334\364\330\333\363\330\333\363\327\332\363"
"\326\332\363\327\331\362\326\331\362\325\331\362\325\330\362\324\330\363"
"\323\327\362\323\327\362\322\326\362\322\326\362\322\325\361\321\325\361"
"\321\325\361\320\323\361\320\323\361\317\323\361\317\322\360\316\322\361"
"\316\322\360\315\321\360\314\321\360\314\320\360\313\320\360\313\317\360"
"\312\317\357\311\317\357\312\316\357\311\316\357\310\315\356\310\314\356"
"\307\314\356\307\313\357\306\313\357\306\313\356\305\312\356\305\312\355"
"\305\311\356\304\311\356\303\311\355\302\310\355\302\307\355\302\307\355"
"\301\306\355\301\306\355\300\305\354\300\305\354\277\305\354\277\304\354"
"\276\304\354\276\304\354\276\302\353\275\303\354\274\302\353\274\301\353"
"\274\301\353\273\301\353\272\300\353\272\277\353\271\277\352\271\276\352"
"\270\276\352\270\276\352\267\275\351\267\274\352\266\275\352\265\273\352"
"\265\273\351\264\273\351\265\272\351\263\272\351\263\272\351\263\271\350"
"\263\271\351\262\271\350\261\267\350\260\267\350\260\266\350\257\267\350"
"\256\270\350\255\272\350\254\273\351\252\275\351\251\275\351\247\300\351"
"\247\301\351\244\303\352\243\306\352\242\310\353\240\312\352\237\315\352"
"\236\317\352\234\322\353\233\324\353\232\330\353\230\332\354\227\336\354"
"\225\341\354\224\345\354\223\350\354\222\354\355\220\355\353\217\355\347"
"\216\355\344\214\355\337\212\356\334\211\356\330\210\357\324\207\356\317"
"\205\357\313\204\357\307\202\360\302\200\357\275\200\360\270\177\360\264"
"~\360\256|\361\251z\360\244y\361\237x\361\230v\362\223u\362\215s\362\207"
"r\362\200p\362|o\362un\362np\363ls\363kw\363i|\364g\200\363f\204\363d\210"
"\364b\215\364a\222\364`\227\365^\234\365]\241\365\\\246\365Z\254\366Y\261"
"\365W\267\366U\275\366T\303\366S\311\367Q\317\367O\325\367N\334\370L\343"
"\370J\352\370I\360\370G\367\371F\370\362E\371\354C\370\345A\371\336@\371"
"\326?\372\316=\372\307;\372\277:\372\2679\373\2577\373\2475\373\2373\373"
"\2261\373\2161\373\206.\373}-\374u,\375k*\374b(\375Y'\375O%\375E#\375;\""
"\3762\40\376(\37\376\35"
};
char PrecipPalette[256*3] = {
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0\x0"
"\x0\x0\x0\x1\x42\xc1\x3\x54\xd2\x8\x68\xe0\x10\x7c\xec\x1b\x91\xf5"
"\x29\xa4\xfc\x39\xb7\xff\x4a\xc9\xff\x5d\xd8\xfb\x71\xe6\xf5\x85\xf1\xeb"
"\x9a\xf9\xdf\xad\xfd\xd0\xbf\xff\xc0\xd0\xfd\xad\xdf\xf9\x9a\xeb\xf1\x86"
"\xf5\xe6\x71\xfb\xd9\x5d\xff\xc9\x4a\xff\xb7\x39\xfc\xa4\x29\xf5\x91\x1b"
"\xec\x7c\x10\xe0\x68\x8\xd2\x55\x3\xc1\x42\x1\xaf\x31\x2\x9c\x22\x7"
"\x87\x16\xe\x73\xc\x19\x5f\x5\x26\x4c\x2\x36\x3a\x1\x47\x2a\x4\x5a"
"\x1c\xa\x6e\x11\x13\x82\x9\x1f\x96\x3\x2d\xaa\x1\x3d\xbc\x2\x4f\xcd"
"\x6\x63\xdc\xe\x77\xe9\x18\x8b\xf3\x25\x9f\xfa\x34\xb2\xfe\x45\xc4\xff"
"\x58\xd4\xfc\x6c\xe2\xf7\x80\xee\xee\x94\xf7\xe2\xa8\xfc\xd4\xbb\xff\xc4"
"\xcc\xfe\xb2\xdb\xfa\x9f\xe8\xf3\x8b\xf2\xe9\x77\xfa\xdc\x63\xfe\xcd\x4f"
"\xff\xbc\x3d\xfd\xaa\x2d\xf7\x96\x1f\xef\x82\x13\xe4\x6e\xa\xd6\x5a\x4"
"\xc6\x47\x1\xb4\x36\x2\xa1\x26\x5\x8d\x19\xc"
};

+ 13
- 11
dsp.c View File

@@ -51,7 +51,7 @@ static double FreqLine = 1.0;
static double FreqOsc;
static double K1, K2;

// Check the input sample rate and calculate some constants
// Check the sample rate and calculate some constants
int init_dsp(double F) {
if(F > Fi) return(1);
if(F < Fp) return(-1);
@@ -72,8 +72,7 @@ static inline double Phase(double I, double Q) {
double angle, r;
int s;

if(I == 0.0 && Q == 0.0)
return(0.0);
if(I == 0.0 && Q == 0.0) return(0.0);

if (Q < 0) {
s = -1;
@@ -90,10 +89,11 @@ static inline double Phase(double I, double Q) {
angle = 0.75 - 0.25 * r;
}

if(s > 0)
if(s > 0){
return(angle);
else
}else{
return(-angle);
}
}

/* Phase locked loop
@@ -179,7 +179,7 @@ int getpixelv(float *pvbuff, int count) {

double mult;

// Sub-pixel offset
// Gaussian resampling factor
mult = (double) Fi / Fe * FreqLine;
int m = RSFilterLen / mult + 1;

@@ -205,16 +205,18 @@ int getpixelv(float *pvbuff, int count) {
idxam += shift;
nam -= shift;
}
return(count);
}

// Get an entire row of pixels, aligned with sync markers
double minDoppler = 100;
// FIXME: occasionally skips noisy lines
int getpixelrow(float *pixelv, int nrow, int *zenith) {
static float pixels[PixelLine + SyncFilterLen];
static int npv;
static int synced = 0;
static double max = 0.0;
static double minDoppler = 100;

double corr, ecorr, lcorr;
int res;
@@ -233,12 +235,12 @@ int getpixelrow(float *pixelv, int nrow, int *zenith) {

// Calculate the frequency offset
ecorr = fir(pixelv, Sync, SyncFilterLen);
corr = fir(&(pixelv[1]), Sync, SyncFilterLen - 1);
lcorr = fir(&(pixelv[2]), Sync, SyncFilterLen - 2);
corr = fir(&pixelv[1], Sync, SyncFilterLen - 1);
lcorr = fir(&pixelv[2], Sync, SyncFilterLen - 2);
FreqLine = 1.0+((ecorr-lcorr) / corr / PixelLine / 4.0);

// Find the point in which ecorr and lcorr intercept, make sure that it's not too perfect
if(fabs(lcorr - ecorr) < minDoppler && fabs(lcorr - ecorr) > 1){
// Find the point in which ecorr and lcorr intercept
if(fabs(lcorr - ecorr) < minDoppler && fabs(lcorr - ecorr) > 2){
minDoppler = fabs(lcorr - ecorr);
*zenith = nrow;
}


+ 1
- 1
falsecolor.conf View File

@@ -2,7 +2,7 @@
23 78 37
240 250 255
50
0.5
10
24
34
14


+ 42
- 61
fcolor.c View File

@@ -19,106 +19,87 @@

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

#include "offsets.h"

typedef struct {
float r, g, b, a;
} rgba;
float r, g, b;
} rgb_t;

extern rgb_t RGBcomposite(rgb_t top, float top_a, rgb_t bottom, float bottom_a);

static struct {
rgba Sea;
rgba Land;
rgba Cloud;
int Seaintensity;
float Seaoffset;
int Landthreshold;
int Landintensity;
int Landoffset;
int Cloudthreshold;
int Cloudintensity;
rgb_t Sea, Land, Cloud;
int Seaintensity, Seaoffset;
int Landthreshold, Landintensity, Landoffset;
int Cloudthreshold, Cloudintensity;
} fcinfo = {
{28, 44, 95, 1},
{23, 78, 37, 1},
{240, 250, 255, 1},
50,
0.5,
24,
34,
14,
141,
114
{28, 44, 95},
{23, 78, 37},
{240, 250, 255},
50, 10,
24, 34, 14,
141, 114
};

// Read the config file
void readfcconf(char *file) {
int readfcconf(char *file) {
FILE *fin;

fin = fopen(file, "r");

if (fin == NULL)
return;
return 0;

fscanf(fin, "%g %g %g\n", &fcinfo.Sea.r, &fcinfo.Sea.g, &fcinfo.Sea.b);
fscanf(fin, "%g %g %g\n", &fcinfo.Land.r, &fcinfo.Land.g, &fcinfo.Land.b);
fscanf(fin, "%g %g %g\n", &fcinfo.Cloud.r, &fcinfo.Cloud.g, &fcinfo.Cloud.b);
fscanf(fin, "%d\n", &fcinfo.Seaintensity);
fscanf(fin, "%g\n", &fcinfo.Seaoffset);
fscanf(fin, "%d\n", &fcinfo.Seaoffset);
fscanf(fin, "%d\n", &fcinfo.Landthreshold);
fscanf(fin, "%d\n", &fcinfo.Landintensity);
fscanf(fin, "%d\n", &fcinfo.Landoffset);
fscanf(fin, "%d\n", &fcinfo.Cloudthreshold);
fscanf(fin, "%d", &fcinfo.Cloudintensity);

fclose(fin);
};

// RGBA Composite
rgba composite(rgba top, rgba bottom){
rgba composite;
composite.r = MCOMPOSITE(top.r, top.a, bottom.r, bottom.a);
composite.g = MCOMPOSITE(top.g, top.a, bottom.g, bottom.a);
composite.b = MCOMPOSITE(top.b, top.a, bottom.b, bottom.a);
composite.a = bottom.a == 1 || top.a == 1 ? 1 : (top.a+bottom.a)/2;
return composite;
}
return 1;
};

void falsecolor(float vis, float temp, float *r, float *g, float *b){
rgba buffer;
fcinfo.Land.a = 0.0f; fcinfo.Sea.a = 0.0f;
rgb_t falsecolor(float vis, float temp, float *r, float *g, float *b){
rgb_t buffer;
float land = 0, sea, cloud;

// Calculate intensity of sea and land
fcinfo.Sea.a = CLIP(vis, 0, 20)/fcinfo.Seaintensity + fcinfo.Seaoffset;
if(vis > fcinfo.Landthreshold) fcinfo.Land.a = CLIP(vis-fcinfo.Landoffset, 0, fcinfo.Landintensity)/fcinfo.Landintensity;
// Calculate intensity of sea
sea = CLIP(vis+fcinfo.Seaoffset, 0, fcinfo.Seaintensity)/fcinfo.Seaintensity;

// Composite land on top of sea
buffer = composite(fcinfo.Land, fcinfo.Sea);
buffer.a = 1.0f;
// Land
if(vis > fcinfo.Landthreshold)
land = CLIP(vis+fcinfo.Landoffset, 0, fcinfo.Landintensity)/fcinfo.Landintensity;
// Composite land on sea
buffer = RGBcomposite(fcinfo.Land, land, fcinfo.Sea, sea);

// Composite clouds on top
fcinfo.Cloud.a = CLIP(temp-fcinfo.Cloudthreshold, 0, fcinfo.Cloudintensity)/fcinfo.Cloudintensity;
buffer = composite(fcinfo.Cloud, buffer);
cloud = CLIP(temp-fcinfo.Cloudthreshold, 0, fcinfo.Cloudintensity)/fcinfo.Cloudintensity;
buffer = RGBcomposite(fcinfo.Cloud, cloud, buffer, 1);

*r = buffer.r;
*g = buffer.g;
*b = buffer.b;
return buffer;
}

// GVI (global vegetation index) false color
void Ngvi(float **prow, int nrow) {
printf("Computing GVI false color\n");
fflush(stdout);
printf("Computing GVI false color");

for (int n = 0; n < nrow; n++) {
float *pixelv = prow[n];
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;
pv = CLIP(pv, 0, 255);
for (int i = 0; i < CH_WIDTH; i++) {
double gvi = (pixelv[i + CHA_OFFSET] - pixelv[i + CHB_OFFSET])/
(pixelv[i + CHA_OFFSET] + pixelv[i + CHB_OFFSET]);

pixelv[i + CHB_OFFSET] = pv;
gvi = (gvi + 0.1) * 340.0;
pixelv[i + CHB_OFFSET] = CLIP(gvi, 0, 255);
}
}
printf("\nDone\n");
};

+ 106
- 31
image.c View File

@@ -21,6 +21,7 @@
#include <string.h>
#include <sndfile.h>
#include <math.h>
#include <stdlib.h>

#include "offsets.h"
#include "messages.h"
@@ -28,12 +29,27 @@
#define REGORDER 3
typedef struct {
double cf[REGORDER + 1];
} rgparam;
} rgparam_t;

typedef struct {
float *prow[3000]; // Row buffers
int nrow; // Number of rows
int chA, chB; // ID of each channel
char name[256]; // Stripped filename
} image_t;

typedef struct {
char *type; // Output image type
char *effects;
int satnum; // The satellite number
char *map; // Path to a map file
char *path; // Output directory
} options_t;

extern void polyreg(const int m, const int n, const double x[], const double y[], double c[]);

// Compute regression
static void rgcomp(double x[16], rgparam * rgpr) {
static void rgcomp(double x[16], rgparam_t * rgpr) {
// { 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 };

@@ -41,7 +57,7 @@ static void rgcomp(double x[16], rgparam * rgpr) {
}

// Convert a value to 0-255 based off the provided regression curve
static double rgcal(float x, rgparam *rgpr) {
static double rgcal(float x, rgparam_t *rgpr) {
double y, p;
int i;
@@ -66,13 +82,14 @@ void histogramEqualise(float **prow, int nrow, int offset, int width){
// Find min/max points
int min = -1, max = -1;
for(int i = 5; i < 250; i++){
if(histogram[i]/width/(nrow/255.0) > 1.0){
if(min == -1) min = i;
if(histogram[i]/width/(nrow/255.0) > 0.2){
if(min == -1)
min = i;
max = i;
}
}

//printf("Min Value: %i, Max Value %i\n", min, max);
//printf("Column %i-%i: Min: %i, Max %i\n", offset, offset+width, min, max);

// Spread values to avoid overshoot
min -= 5; max += 5;
@@ -80,11 +97,11 @@ void histogramEqualise(float **prow, int nrow, int offset, int width){
// 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;
prow[y][x+offset] = CLIP((prow[y][x+offset]-min) / (max-min) * 255.0, 0, 255);
}

// Brightness calibrate, including telemetry
void calibrateBrightness(float **prow, int nrow, int offset, int width, int telestart, rgparam regr[30]){
void calibrateBrightness(float **prow, int nrow, int offset, int width, int telestart, rgparam_t regr[30]){
offset -= SYNC_WIDTH+SPC_WIDTH;

for (int n = 0; n < nrow; n++) {
@@ -94,19 +111,19 @@ void calibrateBrightness(float **prow, int nrow, int offset, int width, int tele
float pv = pixelv[i + offset];

// Blend between the calculated regression curves
/* TODO: this can actually make the image look *worse*
/* FIXME: this can actually make the image look *worse*
* if the signal has a constant input gain.
*/
int k, kof;
/*int k, kof;
k = (n - telestart) / FRAME_LEN;
if (k >= nbtele)
k = nbtele - 1;
kof = (n - telestart) % FRAME_LEN;

if (kof < 64) {
if (k < 1) {
pv = rgcal(pv, &(regr[k]));
} else {
if (k < 1) {*/
pv = rgcal(pv, &(regr[4]));
/*} else {
pv = rgcal(pv, &(regr[k])) * (64 + kof) / FRAME_LEN +
rgcal(pv, &(regr[k - 1])) * (64 - kof) / FRAME_LEN;
}
@@ -118,7 +135,7 @@ void calibrateBrightness(float **prow, int nrow, int offset, int width, int tele
rgcal(pv, &(regr[k + 1])) * (kof - 64) / FRAME_LEN;
}
}
*/
pv = CLIP(pv, 0, 255);
pixelv[i + offset] = pv;
}
@@ -129,7 +146,7 @@ void calibrateBrightness(float **prow, int nrow, int offset, int width, int tele
int calibrate(float **prow, int nrow, int offset, int width) {
double teleline[3000] = { 0.0 };
double wedge[16];
rgparam regr[30];
rgparam_t regr[30];
int n, k;
int mtelestart = 0, telestart;
int channel = -1;
@@ -204,13 +221,13 @@ int calibrate(float **prow, int nrow, int offset, int width) {
* wedges, the wedge with the closest match will
* be the channel ID
*/
float min = 10000;
float min = -1;
for (j = 0; j < 6; j++) {
float df;

df = tele[15] - tele[j];
df *= df;
if (df < min) {
if (df < min || min == -1) {
channel = j;
min = df;
}
@@ -242,7 +259,6 @@ int calibrate(float **prow, int nrow, int offset, int width) {
}

// --- Temperature Calibration --- //
extern int satnum;
#include "satcal.h"

typedef struct {
@@ -250,10 +266,10 @@ typedef struct {
double Cs;
double Cb;
int ch;
} tempparam;
} tempparam_t;

// IR channel temperature compensation
static void tempcomp(double t[16], int ch, tempparam *tpr) {
static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) {
double Tbb, T[4];
double C;

@@ -285,7 +301,7 @@ static void tempcomp(double t[16], int ch, tempparam *tpr) {
}

// IR channel temperature calibration
static double tempcal(float Ce, tempparam * rgpr) {
static double tempcal(float Ce, int satnum, tempparam_t * rgpr) {
double Nl, Nc, Ns, Ne;
double T, vc;

@@ -308,23 +324,82 @@ static double tempcal(float Ce, tempparam * rgpr) {
}

// Temperature calibration wrapper
void temperature(float **prow, int nrow, int ch, int offset){
tempparam temp;
void temperature(options_t *opts, image_t *img, int offset, int width){
tempparam_t temp;

printf("Temperature... ");
fflush(stdout);

tempcomp(tele, ch, &temp);

for (int n = 0; n < nrow; n++) {
float *pixelv = prow[n];
tempcomp(tele, img->chB, opts->satnum - 15, &temp);

for (int i = 0; i < CH_WIDTH; i++) {
float pv = tempcal(pixelv[i + offset], &temp);
for (int y = 0; y < img->nrow; y++) {
float *pixelv = img->prow[y];

pv = CLIP(pv, 0, 255);
pixelv[i + offset] = pv;
for (int x = 0; x < width; x++) {
float pv = tempcal(pixelv[x + offset], opts->satnum - 15, &temp);
pixelv[x + offset] = CLIP(pv, 0, 255);
}
}
printf("Done\n");
}

void distrib(options_t *opts, image_t *img, char *chid) {
int max = 0;

// Options
options_t options;
options.path = opts->path;

// Image options
image_t distrib;
strcpy(distrib.name, img->name);
distrib.nrow = 256;

// Assign memory
for(int i = 0; i < 256; i++)
distrib.prow[i] = (float *) malloc(sizeof(float) * 256);

for(int n = 0; n < img->nrow; n++) {
float *pixelv = img->prow[n];
for(int i = 0; i < CH_WIDTH; i++) {
int y = CLIP((int)pixelv[i + CHA_OFFSET], 0, 255);
int x = CLIP((int)pixelv[i + CHB_OFFSET], 0, 255);
distrib.prow[y][x]++;
if(distrib.prow[y][x] > max)
max = distrib.prow[y][x];
}
}

// Scale to 0-255
for(int x = 0; x < 256; x++)
for(int y = 0; y < 256; y++)
distrib.prow[y][x] = distrib.prow[y][x] / max * 255.0;
extern int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, char *chid, char *palette);
ImageOut(&options, &distrib, 0, 256, "Distribution", chid, NULL);
}

extern float quick_select(float arr[], int n);

// Recursive biased median denoise
#define TRIG_LEVEL 40
void 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

+ 231
- 517
main.c View File

@@ -24,321 +24,286 @@
#include <libgen.h>
#include <math.h>
#include <sndfile.h>
#include <png.h>
#include <errno.h>

#include "messages.h"
#include "offsets.h"
#include "palette.h"

#define FAILURE 0
#define SUCCESS 1

typedef struct {
char *type; // Output image type
char *effects;
int satnum; // The satellite number
char *map; // Path to a map file
char *path; // Output directory
} options_t;

typedef struct {
float *prow[3000]; // Row buffers
int nrow; // Number of rows
int chA, chB; // ID of each channel
char name[256]; // Stripped filename
} image_t;

// DSP
extern int getpixelrow(float *pixelv, int nrow, int *zenith);
extern int init_dsp(double F);

static SNDFILE *inwav;
// I/O
extern int readfcconf(char *file);
extern int readRawImage(char *filename, float **prow, int *nrow);
extern int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, char *chid, char *palette);

static int initsnd(char *filename) {
SF_INFO infwav;
int res;
// Image functions
extern int calibrate(float **prow, int nrow, int offset, int width);
extern void histogramEqualise(float **prow, int nrow, int offset, int width);
extern void temperature(options_t *opts, image_t *img, int offset, int width);
extern int Ngvi(float **prow, int nrow);
extern void denoise(float **prow, int nrow, int offset, int width);
extern void distrib(options_t *opts, image_t *img, char *chid);

// Open audio file
infwav.format = 0;
inwav = sf_open(filename, SFM_READ, &infwav);
if (inwav == NULL) {
fprintf(stderr, ERR_FILE_READ, filename);
return(0);
}
// Palettes
extern char GviPalette[256*3];
extern char TempPalette[256*3];

res = init_dsp(infwav.samplerate);
printf("Input file: %s\n", filename);
if(res < 0) {
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);
}
printf("Input sample rate: %d\n", infwav.samplerate);
// Row where the satellite is closest to the observer
int zenith = 0;
// Audio file
static SNDFILE *audioFile;

if (infwav.channels != 1) {
fprintf(stderr, "Too many channels in input file: %d\n", infwav.channels);
return(0);
}
static int initsnd(char *filename);
int getsample(float *sample, int nb);
static int processAudio(char *filename, options_t *opts);
static void usage(void);

return(1);
}
int main(int argc, char **argv) {
fprintf(stderr, VERSION"\n");

// Get samples from the wave file
int getsample(float *sample, int nb) {
return sf_read_float(inwav, sample, nb);
}
if(argc == 1)
usage();

static png_text text_ptr[] = {
{PNG_TEXT_COMPRESSION_NONE, "Software", VERSION},
{PNG_TEXT_COMPRESSION_NONE, "Channel", NULL, 0},
{PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20}
};

int mapOverlay(char *filename, float **crow, int nrow, int zenith, int MCIR) {
FILE *fp = fopen(filename, "rb");

// Create reader
png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
if(!png) return 0;
png_infop info = png_create_info_struct(png);
if(!info) 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 != 1040){
fprintf(stderr, "Map must be 1040px wide.\n");
return 0;
}else if(bit_depth != 16){
fprintf(stderr, "Map must be 16 bit color.\n");
return 0;
}else if(color_type != PNG_COLOR_TYPE_RGB){
fprintf(stderr, "Map must be RGB.\n");
return 0;
}else if(zenith > height/2 || nrow-zenith > height/2){
fprintf(stderr, "WARNING: Map is too short to cover entire image\n");
// Check if there are actually any input files
if(optind == argc){
fprintf(stderr, "No input files provided.\n");
usage();
}

// Create row buffers
png_bytep *mapRows = NULL;
mapRows = (png_bytep *) malloc(sizeof(png_bytep) * height);
for(int y = 0; y < height; y++) mapRows[y] = (png_byte *) malloc(png_get_rowbytes(png, info));

// Read image
png_read_image(png, mapRows);

// Tidy up
fclose(fp);
png_destroy_read_struct(&png, &info, NULL);

int mapOffset = (height/2)-zenith;
for(int y = 0; y < nrow; y++) {
for(int x = 49; x < width - 82; x++){
// Maps are 16 bit / channel
png_bytep px = &mapRows[CLIP(y + mapOffset, 0, height)][x * 6];
uint16_t r = (uint16_t)(px[0] << 8) | px[1];
uint16_t g = (uint16_t)(px[2] << 8) | px[3];
uint16_t b = (uint16_t)(px[4] << 8) | px[5];

// Pointers to the current pixel in each channel
float *cha = &crow[y][(x+36) * 3];
float *chb = &crow[y][(x+CHB_OFFSET-49) * 3];

// Fill in map
if(MCIR){
// Sea
cha[0] = 42; cha[1] = 53; cha[2] = 105;

// Land
if(g > 128){
float vegetation = (r-160) / 96.0;
cha[0] = 120; cha[1] = vegetation*60.0 + 100; cha[2] = 95;
options_t opts = { "r", "", 19, "", "." };

// Parse arguments
int opt;
while ((opt = getopt(argc, argv, "c:m:d:i:s:e:")) != EOF) {
switch (opt) {
case 'd':
opts.path = optarg;
break;
case 'c':
readfcconf(optarg);
break;
case 'm':
opts.map = optarg;
break;
case 'i':
opts.type = optarg;
break;
case 's':
opts.satnum = atoi(optarg);
if(opts.satnum < 15 || opts.satnum > 19){
fprintf(stderr, "Invalid satellite number, it must be the range 15-19\n");
exit(EPERM);
}
}

// Color -> alpha: composite
int composite = r + g + b;
// Color -> alpha: flattern color depth
float factor = (255 * 255 * 2.0) / composite;
r *= factor; g *= factor; b *= factor;
// Color -> alpha: convert black to alpha
float alpha = CLIP(abs(0 - composite) / 65535.0, 0, 1);

// Map overlay on channel A
cha[0] = MCOMPOSITE(r/257, alpha, cha[0], 1);
cha[1] = MCOMPOSITE(g/257, alpha, cha[1], 1);
cha[2] = MCOMPOSITE(b/257, alpha, cha[2], 1);

// Map overlay on channel B
if(!MCIR){
chb[0] = MCOMPOSITE(r/257, alpha, chb[0], 1);
chb[1] = MCOMPOSITE(g/257, alpha, chb[1], 1);
chb[2] = MCOMPOSITE(b/257, alpha, chb[2], 1);
}

// Cloud overlay on channel A
if(MCIR){
float cloud = (chb[0] - 110) / 120;
cha[0] = MCOMPOSITE(240, cloud, cha[0], 1);
cha[1] = MCOMPOSITE(250, cloud, cha[1], 1);
cha[2] = MCOMPOSITE(255, cloud, cha[2], 1);
}
break;
case 'e':
opts.effects = optarg;
break;
default:
usage();
}
}

// Process the files
for (; optind < argc; optind++) {
processAudio(argv[optind], &opts);
}

return 1;
exit(0);
}

// Row where to satellite reaches peak elevation
int zenith = 0;
static int processAudio(char *filename, options_t *opts){
// Image info struct
image_t img;

int applyPalette(float **crow, int nrow, int width, char *palette){
if(palette == NULL) return 0;
// 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", "mid-infrared", "thermal-infrared", "thermal-infrared", "mid-infrared" }
};

for(int y = 0; y < nrow; y++){
for(int x = 0; x < width; x++){
float *px = &crow[y][x * 3];
px[0] = palette[(int)CLIP(px[0], 0, 255)*3 + 0];
px[1] = palette[(int)CLIP(px[1], 0, 255)*3 + 1];
px[2] = palette[(int)CLIP(px[2], 0, 255)*3 + 2];
// 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), "%[^.].%s", img.name, extension);

if(strcmp(extension, "png") == 0){
// Read PNG into image buffer
printf("Reading %s", filename);
if(readRawImage(filename, img.prow, &img.nrow) == 0){
fprintf(stderr, "Skipping %s; see above.\n", img.name);
return FAILURE;
}
}

return 1;
}
}else{
// Attempt to open the audio file
if (initsnd(filename) == 0)
exit(EPERM);

// Build image
for (img.nrow = 0; img.nrow < 3000; img.nrow++) {
// Allocate memory for this row
img.prow[img.nrow] = (float *) malloc(sizeof(float) * 2150);
// Write into memory and break the loop when there are no more samples to read
if (getpixelrow(img.prow[img.nrow], img.nrow, &zenith) == 0)
break;

static int ImageOut(char *filename, char *chid, float **prow, int nrow, int width, int offset, char *palette, char *effects, char *mapFile) {
FILE *pngfile;
fprintf(stderr, "Row: %d\r", img.nrow);
fflush(stderr);
}

// Reduce the width of the image to componsate for the missing telemetry
if(CONTAINS(effects, 't')) width -= TOTAL_TELE;
// Close stream
sf_close(audioFile);
}

// 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);
fprintf(stderr, ERR_PNG_WRITE);
return(0);
}
png_infop info_ptr = png_create_info_struct(png_ptr);
if (!info_ptr) {
png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
fprintf(stderr, ERR_PNG_INFO);
return(0);
}
printf("\nTotal rows: %d\n", img.nrow);

// 8 bit RGB image
png_set_IHDR(png_ptr, info_ptr, width, nrow,
8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);

text_ptr[1].text = chid;
text_ptr[1].text_length = strlen(chid);
png_set_text(png_ptr, info_ptr, text_ptr, 3);
png_set_pHYs(png_ptr, info_ptr, 4160, 4160, PNG_RESOLUTION_METER);

// Init I/O
pngfile = fopen(filename, "wb");
if (!pngfile) {
fprintf(stderr, ERR_FILE_WRITE, filename);
return(1);
}
png_init_io(png_ptr, pngfile);
png_write_info(png_ptr, info_ptr);
// Fallback for detecting the zenith
// TODO: encode zenith in raw images
if(opts->map != NULL && opts->map[0] != '\0' && zenith == 0){
fprintf(stderr, "Guessing zenith in image, map will most likely be misaligned.\n");
zenith = img.nrow / 2;
}

// Move prow into crow, crow ~ color rows
float *crow[3000];
for(int i = 0; i < nrow; i++){
crow[i] = (float *) malloc(sizeof(float) * 2080 * 3);
// Calibrate
img.chA = calibrate(img.prow, img.nrow, CHA_OFFSET, CH_WIDTH);
img.chB = calibrate(img.prow, img.nrow, CHB_OFFSET, 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]);

for(int x = 0; x < 2080; x++)
crow[i][x*3] = crow[i][x*3 + 1] = crow[i][x*3 + 2] = prow[i][x];
// Denoise
if(CONTAINS(opts->effects, 'd')){
denoise(img.prow, img.nrow, CHA_OFFSET, CH_WIDTH);
denoise(img.prow, img.nrow, CHB_OFFSET, CH_WIDTH);
}

applyPalette(crow, nrow, 2080, palette);
// Temperature
if (CONTAINS(opts->type, 't') && img.chB >= 4) {
temperature(opts, &img, CHB_OFFSET, CH_WIDTH);
ImageOut(opts, &img, CHB_OFFSET, CH_WIDTH, "Temperature", "t", (char *)TempPalette);
}

if(mapFile != NULL && mapFile[0] != '\0'){
if(mapOverlay(mapFile, crow, nrow, zenith, strcmp(chid, "MCIR") == 0) == 0){
fprintf(stderr, "Skipping MCIR generation; see above.\n");
return 0;
// False color image
if(CONTAINS(opts->type, 'c')){
if (img.chA == 2 && img.chB >= 4) { // Normal false color
// TODO: use real MSA
// TODO: provide more than just "natural" color images
ImageOut(opts, &img, 0, CH_WIDTH, "False Color", "c", NULL);
} else if (img.chB == 2) { // GVI (global vegetation index) false color
Ngvi(img.prow, img.nrow);
ImageOut(opts, &img, CHB_OFFSET, CH_WIDTH, "GVI False Color", "c", (char *)GviPalette);
} else {
fprintf(stderr, "Skipping False Color generation; lacking required channels.\n");
}
}else if(strcmp(chid, "MCIR") == 0){
fprintf(stderr, "Skipping MCIR generation; no map provided.\n");
return 0;
}

extern void falsecolor(float vis, float temp, float *r, float *g, float *b);

printf("Writing %s", filename);

int fcimage = strcmp(chid, "False Color") == 0;

// Build RGB image
for (int n = 0; n < nrow; n++) {
png_color pix[width];

for (int i = 0; i < width; i++) {
float *px = &crow[n][offset*3 + i*3];
if(fcimage){
float r = 0, g = 0, b = 0;
falsecolor(prow[n][i + CHA_OFFSET], prow[n][i + CHB_OFFSET], &r, &g, &b);
pix[i].red = r;
pix[i].green = g;
pix[i].blue = b;
}else{
pix[i].red = px[0];
pix[i].green = px[1];
pix[i].blue = px[2];
}
}
png_write_row(png_ptr, (png_bytep) pix);
}
// MCIR
if (CONTAINS(opts->type, 'm'))
ImageOut(opts, &img, 0, IMG_WIDTH, "MCIR", "m", NULL);

// Tidy up
png_write_end(png_ptr, info_ptr);
fclose(pngfile);
printf("\nDone\n");
png_destroy_write_struct(&png_ptr, &info_ptr);
// Histogram equalise
if(CONTAINS(opts->effects, 'h')){
histogramEqualise(img.prow, img.nrow, CHA_OFFSET, CH_WIDTH);
histogramEqualise(img.prow, img.nrow, CHB_OFFSET, CH_WIDTH);
}

return(1);
}
// Raw image
if (CONTAINS(opts->type, 'r')) {
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, IMG_WIDTH, desc, "r", NULL);
}

// Outputs a image with the value distribution between channel A and B
static void distrib(char *filename, float **prow, int nrow) {
float *distrib[256];
int max = 0;

// Assign memory
for(int i = 0; i < 256; i++)
distrib[i] = (float *) malloc(sizeof(float) * 256);

for(int n = 0; n < nrow; n++) {
float *pixelv = prow[n];
for(int i = 0; i < CH_WIDTH; i++) {
int y = (int)(pixelv[i + CHA_OFFSET]);
int x = (int)(pixelv[i + CHB_OFFSET]);
distrib[y][x] += 1;
if(distrib[y][x] > max) max = distrib[y][x];
}
// Channel A
if (CONTAINS(opts->type, 'a')) {
sprintf(desc, "%s (%s)", ch.id[img.chA], ch.name[img.chA]);
ImageOut(opts, &img, CHA_OFFSET, CH_WIDTH, desc, ch.id[img.chA], NULL);
}

// Scale to 0-255
for(int x = 0; x < 256; x++)
for(int y = 0; y < 256; y++)
distrib[y][x] = distrib[y][x] / max * 255;
// Channel B
if (CONTAINS(opts->type, 'b')) {
sprintf(desc, "%s (%s)", ch.id[img.chB], ch.name[img.chB]);
ImageOut(opts, &img, CHB_OFFSET, CH_WIDTH, desc, ch.id[img.chA], NULL);
}

ImageOut(filename, "Brightness distribution", distrib, 256, 256, 0, NULL, 0, NULL);
// Distribution image
if (CONTAINS(opts->type, 'd'))
distrib(opts, &img, "d");
return SUCCESS;
}

extern int calibrate(float **prow, int nrow, int offset, int width);
extern void histogramEqualise(float **prow, int nrow, int offset, int width);
extern void temperature(float **prow, int nrow, int ch, int offset);
extern int Ngvi(float **prow, int nrow);
extern void readfcconf(char *file);
extern int optind;
extern char *optarg;
static int initsnd(char *filename) {
SF_INFO infwav;
int res;

// Default to NOAA 19
int satnum = 4;
// Open audio file
infwav.format = 0;
audioFile = sf_open(filename, SFM_READ, &infwav);
if (audioFile == NULL) {
fprintf(stderr, ERR_FILE_READ, filename);
return FAILURE;
}

res = init_dsp(infwav.samplerate);
printf("Input file: %s\n", filename);
if(res < 0) {
fprintf(stderr, "Input sample rate too low: %d\n", infwav.samplerate);
return FAILURE;
}else if(res > 0) {
fprintf(stderr, "Input sample rate too high: %d\n", infwav.samplerate);
return FAILURE;
}
printf("Input sample rate: %d\n", infwav.samplerate);

// TODO: accept stereo audio
if (infwav.channels != 1) {
fprintf(stderr, "Too many channels in input file: %d\n", infwav.channels);
return FAILURE;
}

return SUCCESS;
}

// Read samples from the wave file
int getsample(float *sample, int nb) {
return sf_read_float(audioFile, sample, nb);
}

static void usage(void) {
printf("Aptdec [options] audio files ...\n"
fprintf(stderr,
"Aptdec [options] audio files ...\n"
"Options:\n"
" -e [t|h] Enhancements\n"
" -e [t|h] Effects\n"
" t: Crop telemetry\n"
" h: Histogram equalise\n"
" -i [r|a|b|c|t] Output image type\n"
" d: Denoise\n"
" -i [r|a|b|c|t] Output image\n"
" r: Raw\n"
" a: Channel A\n"
" b: Channel B\n"
@@ -349,257 +314,6 @@ static void usage(void) {
" -s [15-19] Satellite number\n"
" -c <file> False color config file\n"
" -m <file> Map file\n");
exit(1);
}

int readRawImage(char *filename, float **prow, int *nrow) {
FILE *fp = fopen(filename, "r");

// Create reader
png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
if(!png) return 0;
png_infop info = png_create_info_struct(png);
if(!info) 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 != 2080){
fprintf(stderr, "Raw image must be 2080px wide.\n");
return 0;
}else if(bit_depth != 8){
fprintf(stderr, "Raw image must have 8 bit color.\n");
return 0;
}else if(color_type != PNG_COLOR_TYPE_GRAY){
fprintf(stderr, "Raw image must be grayscale.\n");
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 main(int argc, char **argv) {
char pngfilename[1024];
char name[128];
char pngdirname[128] = "";
char mapFile[256];
char *extension;

// Default to a raw image, with no enhancements
char imgopt[20] = "r";
char enhancements[20] = "";

// Image buffer
float *prow[3000];
int nrow;

// Mapping between telemetry wedge value and channel
static struct {
char *id[7];
char *name[7];
} ch = {
{ "?", "1", "2", "3A", "4", "5", "3B" },
{ "unknown", "visble", "near-infrared", "mid-infrared", "thermal-infrared", "thermal-infrared", "mid-infrared" }
};

// The active sensor in each channel
int chA, chB;

// Print version
printf(VERSION"\n");

// Print usage if there are no arguments
if(argc == 1)
usage();

int c;
while ((c = getopt(argc, argv, "c:m:d:i:s:e:")) != EOF) {
switch (c) {
// Output directory name
case 'd':
strcpy(pngdirname, optarg);
break;
// False color config file
case 'c':
readfcconf(optarg);
break;
// Map file
case 'm':
strcpy(mapFile, optarg);
break;
// Output image type
case 'i':
strncpy(imgopt, optarg, 20);
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, it must be the range [15-19]\n");
exit(1);
}
break;
// Enchancements
case 'e':
strncpy(enhancements, optarg, 20);
break;
default:
usage();
}
}
if(optind == argc){
printf("No input files provided.\n");
usage();
}

// Process the provided files
for (; optind < argc; optind++) {
chA = chB = 0;

// Generate output name
strcpy(pngfilename, argv[optind]);
strcpy(name, basename(pngfilename));
strtok(name, ".");
extension = strtok(NULL, ".");

if (pngdirname[0] == '\0')
strcpy(pngdirname, dirname(pngfilename));

if(strcmp(extension, "png") == 0){
printf("Reading %s", argv[optind]);
if(readRawImage(argv[optind], prow, &nrow) == 0){
fprintf(stderr, "Skipping %s; see above.\n", name);
continue;
}
}else{
// Open sound file, exit if that fails
if (initsnd(argv[optind]) == 0) exit(1);

// Main image building loop
for (nrow = 0; nrow < 3000; nrow++) {
// 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], nrow, &zenith) == 0) break;

printf("Row: %d\r", nrow);
fflush(stdout);
}

// Close sound file
sf_close(inwav);
}

if(zenith == 0 & mapFile[0] != '\0'){
fprintf(stderr, "WARNING: Guessing peak elevation in image, map will most likely not be aligned.\n");
zenith = nrow / 2;
}

printf("\nTotal rows: %d\n", nrow);

// Calibrate
chA = calibrate(prow, nrow, CHA_OFFSET, CH_WIDTH);
chB = calibrate(prow, nrow, CHB_OFFSET, CH_WIDTH);
printf("Channel A: %s (%s)\n", ch.id[chA], ch.name[chA]);
printf("Channel B: %s (%s)\n", ch.id[chB], ch.name[chB]);

// Temperature
if (CONTAINS(imgopt, 't') && chB >= 4) {
// TODO: Doesn't work with channel 4
temperature(prow, nrow, chB, CHB_OFFSET);
sprintf(pngfilename, "%s/%s-t.png", pngdirname, name);
ImageOut(pngfilename, "Temperature", prow, nrow, 909, CHB_OFFSET, (char *)TempPalette, enhancements, mapFile);
}

// 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);
ImageOut(pngfilename, "False Color", prow, nrow, CH_WIDTH, 0, NULL, enhancements, mapFile);
} else if (chB == 2) { // GVI (global vegetation index) false color
Ngvi(prow, nrow);
sprintf(pngfilename, "%s/%s-c.png", pngdirname, name);
ImageOut(pngfilename, "GVI False Color", prow, nrow, CH_WIDTH, CHB_OFFSET, (char *)GviPalette, enhancements, mapFile);
} else {
fprintf(stderr, "Skipping False Color generation; lacking required channels.\n");
}
}

// MCIR
if (CONTAINS(imgopt, 'm')) {
sprintf(pngfilename, "%s/%s-m.png", pngdirname, name);
ImageOut(pngfilename, "MCIR", prow, nrow, 909, CHA_OFFSET, NULL, enhancements, mapFile);
}

// Histogram equalise
if(CONTAINS(enhancements, 'h')){
histogramEqualise(prow, nrow, CHA_OFFSET, CH_WIDTH);
histogramEqualise(prow, nrow, CHB_OFFSET, CH_WIDTH);
}

// Raw image
if (CONTAINS(imgopt, 'r')) {
char channelstr[45];
sprintf(channelstr, "%s (%s) & %s (%s)", ch.id[chA], ch.name[chA], ch.id[chB], ch.name[chB]);

sprintf(pngfilename, "%s/%s-r.png", pngdirname, name);
ImageOut(pngfilename, channelstr, prow, nrow, IMG_WIDTH, 0, NULL, enhancements, mapFile);
}

// Channel A
if (CONTAINS(imgopt, 'a')) {
char channelstr[21];
sprintf(channelstr, "%s (%s)", ch.id[chA], ch.name[chA]);

sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, ch.id[chA]);
ImageOut(pngfilename, channelstr, prow, nrow, CH_WIDTH, CHA_OFFSET, NULL, enhancements, mapFile);
}

// Channel B
if (CONTAINS(imgopt, 'b')) {
char channelstr[21];
sprintf(channelstr, "%s (%s)", ch.id[chB], ch.name[chB]);

sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, ch.id[chB]);
ImageOut(pngfilename, channelstr, prow, nrow, CH_WIDTH , CHB_OFFSET, NULL, enhancements, mapFile);
}

// Distribution image
if (CONTAINS(imgopt, 'd')) {
sprintf(pngfilename, "%s/%s-d.png", pngdirname, name);
distrib(pngfilename, prow, nrow);
}
}

exit(0);
}
exit(EINVAL);
}

+ 56
- 0
median.c View File

@@ -0,0 +1,56 @@
/*
* This Quickselect routine is based on the algorithm described in
* "Numerical recipes in C", Second Edition,
* Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
* This code by Nicolas Devillard - 1998. Public domain.
*/

#define ELEM_SWAP(a, b) { register float t = (a); (a) = (b); (b) = t; }
float quick_select(float arr[], int n) {
int low, median, high;
int middle, ll, hh;

low = 0; high = n-1; median = (low + high) / 2;
for (;;) {
if (high <= low) /* One element only */
return arr[median] ;

if (high == low + 1) { /* Two elements only */
if (arr[low] > arr[high])
ELEM_SWAP(arr[low], arr[high]);
return arr[median];
}

/* Find median of low, middle and high items; swap into position low */
middle = (low + high) / 2;
if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]);
if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]);
if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]);

/* Swap low item (now in position middle) into position (low+1) */
ELEM_SWAP(arr[middle], arr[low+1]);

/* Nibble from each end towards middle, swapping items when stuck */
ll = low + 1;
hh = high;
for (;;) {
do ll++; while (arr[low] > arr[ll]);
do hh--; while (arr[hh] > arr[low]);

if (hh < ll)
break;

ELEM_SWAP(arr[ll], arr[hh]);
}

/* Swap middle item (in position low) back into correct position */
ELEM_SWAP(arr[low], arr[hh]);

/* Re-set active partition */
if (hh <= median)
low = ll;
if (hh >= median)
high = hh - 1;
}
}
#undef ELEM_SWAP

+ 0
- 92
palette.h View File

@@ -1,92 +0,0 @@
/*
* This file is part of Aptdec.
* Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2020
*
* Aptdec is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/

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"
"o\342\330s\344\333w\350\337{\352\344\200\356\350\203\361\354\207\364\360"
"\213\367\364\220\372\370\224\376\376\230\376\375\230\376\372\227\376\366"
"\224\376\362\222\376\357\220\376\352\216\376\347\213\376\343\211\376\340"
"\207\376\334\205\376\330\202\376\325\200\376\321\177\376\315|\376\312z\376"
"\306y\376\302v\376\277t\376\273q\376\267o\376\263m\376\260k\376\254h\376"
"\250f\376\245d\376\241b\376\235_\376\232^\376\226[\376\223Y\376\217W\376"
"\213U\376\207R\376\203Q\376\200N\376}L\376zJ\376vG\376rE\376nB\376k@\376"
"g>\376d<\376`:\376\\7\376X6\376T3\376Q0\376M/\376J-\376F*\376C(\376>%\376"
";$\3767!\3763\37\3760\35\376,\32\376(\31\376$\26\376!\24\376\35\22\376\32"
"\20\376\26\15\376\22\13\376\16\11\376\12\6\376\7\4\376\0\0\373\2\0\366\3"
"\0\362\5\0\355\7\0\350\11\0\343\13\0\336\15\0\332\17\0\325\21\0\320\22\0"
"\313\24\0\307\26\0\303\27\0\275\32\0\270\33\0\264\35\0\257\37\0\253!\0\246"
"#\0\241%\0\234&\0\227)\0\223+\0\216,\0\212.\0\2050\0\2002\0}4\0w6\0s7\0n"
"9\0i;\0e>\0`?\0[A\0VC\0QE\0MG\0HH\0CK\0?M\0:N\0""5P\0""0R\0,S\0'V\0\"X\0"
"\36Z\0\31\\\0\23^\0\17_\0\12a\0\3b\0\6e\2\12i\5\17n\11\23r\14\30v\17\34{"
"\23\40\200\26$\203\31)\207\35.\215\40""1\221$6\225':\232*?\236.C\2431G\250"
"4K\2548O\260;U\265>Y\272A\\\276Eb\303Ie\307Lj\313On\320Ss\324Vw\331Y{\336"
"]\200\342`\203\347d\207\353g\213\357j\217\364n\223\371q\230\376t\222\375"
"{\207\374\205}\373\216r\371\230g\371\241\\\370\254Q\367\265F\365\300;\364"
"\3110\363\323%\362\334\32\361\347\17\357\361\3\356\371\4\353\370\3\347\365"
"\4\342\361\4\335\357\4\331\354\4\323\351\3\317\346\3\313\343\4\306\340\3"
"\301\335\3\275\332\3\270\327\3\263\324\2\257\322\3\252\316\3\245\314\3\241"
"\311\2\234\306\2\230\303\2\223\300\3\216\275\3\211\273\2\205\267\2\200\265"
"\2|\262\2w\257\2s\254\2n\251\2j\246\2e\243\2`\240\2[\235\2W\232\1S\230\2"
"M\225\1I\221\2E\217\1@\214\1;\211\1""7\206\1""1\203\1-\200\0(~\1${\0\37y"
"\0\33u\0\25r\0\21p\0\14l\0\7j\0\3"
};
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"
"\365\366\373\364\366\373\364\365\374\363\365\373\363\364\373\363\364\373"
"\362\363\372\361\363\372\361\362\372\361\361\372\360\361\372\357\361\372"
"\357\361\372\356\360\372\355\357\371\355\357\371\355\356\372\354\356\371"
"\354\355\371\354\355\371\353\355\371\353\354\370\352\353\370\351\353\370"
"\351\353\370\350\352\370\350\351\370\347\351\370\347\351\370\346\350\367"
"\346\350\367\346\350\367\345\346\367\344\346\367\344\346\367\344\345\366"
"\342\346\366\343\344\366\342\345\366\341\344\366\341\343\365\340\343\366"
"\337\342\365\337\341\366\336\341\365\337\341\365\336\340\365\335\340\365"
"\335\337\364\334\337\364\334\336\364\334\336\364\333\336\364\333\335\364"
"\332\334\364\331\334\363\331\334\364\330\333\363\330\333\363\327\332\363"
"\326\332\363\327\331\362\326\331\362\325\331\362\325\330\362\324\330\363"
"\323\327\362\323\327\362\322\326\362\322\326\362\322\325\361\321\325\361"
"\321\325\361\320\323\361\320\323\361\317\323\361\317\322\360\316\322\361"
"\316\322\360\315\321\360\314\321\360\314\320\360\313\320\360\313\317\360"
"\312\317\357\311\317\357\312\316\357\311\316\357\310\315\356\310\314\356"
"\307\314\356\307\313\357\306\313\357\306\313\356\305\312\356\305\312\355"
"\305\311\356\304\311\356\303\311\355\302\310\355\302\307\355\302\307\355"
"\301\306\355\301\306\355\300\305\354\300\305\354\277\305\354\277\304\354"
"\276\304\354\276\304\354\276\302\353\275\303\354\274\302\353\274\301\353"
"\274\301\353\273\301\353\272\300\353\272\277\353\271\277\352\271\276\352"
"\270\276\352\270\276\352\267\275\351\267\274\352\266\275\352\265\273\352"
"\265\273\351\264\273\351\265\272\351\263\272\351\263\272\351\263\271\350"
"\263\271\351\262\271\350\261\267\350\260\267\350\260\266\350\257\267\350"
"\256\270\350\255\272\350\254\273\351\252\275\351\251\275\351\247\300\351"
"\247\301\351\244\303\352\243\306\352\242\310\353\240\312\352\237\315\352"
"\236\317\352\234\322\353\233\324\353\232\330\353\230\332\354\227\336\354"
"\225\341\354\224\345\354\223\350\354\222\354\355\220\355\353\217\355\347"
"\216\355\344\214\355\337\212\356\334\211\356\330\210\357\324\207\356\317"
"\205\357\313\204\357\307\202\360\302\200\357\275\200\360\270\177\360\264"
"~\360\256|\361\251z\360\244y\361\237x\361\230v\362\223u\362\215s\362\207"
"r\362\200p\362|o\362un\362np\363ls\363kw\363i|\364g\200\363f\204\363d\210"
"\364b\215\364a\222\364`\227\365^\234\365]\241\365\\\246\365Z\254\366Y\261"
"\365W\267\366U\275\366T\303\366S\311\367Q\317\367O\325\367N\334\370L\343"
"\370J\352\370I\360\370G\367\371F\370\362E\371\354C\370\345A\371\336@\371"
"\326?\372\316=\372\307;\372\277:\372\2679\373\2577\373\2475\373\2373\373"
"\2261\373\2161\373\206.\373}-\374u,\375k*\374b(\375Y'\375O%\375E#\375;\""
"\3762\40\376(\37\376\35"
};

+ 313
- 0
pngio.c View File

@@ -0,0 +1,313 @@
/*
* This file is part of Aptdec.
* Copyright (c) 2004-2009 Thierry Leconte (F4DWV), Xerbo (xerbo@protonmail.com) 2019-2020
*
* Aptdec is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/

#include <png.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stdint.h>

#include "offsets.h"
#include "messages.h"

typedef struct {
float r, g, b;
} rgb_t;

extern int zenith;
extern char PrecipPalette[256*3];
extern rgb_t applyPalette(char *palette, int val);
extern rgb_t RGBcomposite(rgb_t top, float top_a, rgb_t bottom, float bottom_a);

int mapOverlay(char *filename, rgb_t **crow, int nrow, int zenith, int MCIR) {
FILE *fp = fopen(filename, "rb");

// Create reader
png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
if(!png) return 0;
png_infop info = png_create_info_struct(png);
if(!info) 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 != 1040){
fprintf(stderr, "Map must be 1040px wide.\n");
return 0;
}else if(bit_depth != 16){
fprintf(stderr, "Map must be 16 bit color.\n");
return 0;
}else if(color_type != PNG_COLOR_TYPE_RGB){
fprintf(stderr, "Map must be RGB.\n");
return 0;
}else if(zenith > height/2 || nrow-zenith > height/2){
fprintf(stderr, "WARNING: Map is too short to cover entire image\n");
}

// Create row buffers
png_bytep *mapRows = NULL;
mapRows = (png_bytep *) malloc(sizeof(png_bytep) * height);
for(int y = 0; y < height; y++)
mapRows[y] = (png_byte *) malloc(png_get_rowbytes(png, info));

// Read image
png_read_image(png, mapRows);

// Tidy up
fclose(fp);
png_destroy_read_struct(&png, &info, NULL);

// Map overlay / MCIR / Precipitation
int mapOffset = (height/2)-zenith;
for(int y = 0; y < nrow; y++) {
for(int x = 49; x < width - 82; x++){
// Maps are 16 bit / channel
png_bytep px = &mapRows[CLIP(y + mapOffset, 0, height)][x * 6];
rgb_t map = {
(px[0] << 8) | px[1],
(px[2] << 8) | px[3],
(px[4] << 8) | px[5]
};

// Pixel offsets
int chb = x + CHB_OFFSET - 49;
int cha = x + 36;

// Fill in map
if(MCIR){
if(map.b < 128 && map.g > 128){
// Land
float green = CLIP((map.g-256)/32.0, 0, 1);
float blue = 1-CLIP((map.b-32)/64.0, 0, 1);
crow[y][cha] = (rgb_t){50 + blue*50, 80 + green*70, 64};
}else{
// Sea
crow[y][cha] = (rgb_t){12, 30, 85};
}
}

// Color -> alpha: composite
int composite = map.r + map.g + map.b;
// Color -> alpha: flattern and convert to 8 bits / channel
float factor = (255 * 255 * 2.0) / composite;
map.r *= factor/257.0; map.g *= factor/257.0; map.b *= factor/257.0;
// Color -> alpha: convert black to alpha
float alpha = CLIP(composite / 65535.0, 0, 1);

// Map overlay on channel A
crow[y][cha] = RGBcomposite(map, alpha, crow[y][cha], 1);

// Map overlay on channel B
if(!MCIR)
crow[y][chb] = RGBcomposite(map, alpha, crow[y][chb], 1);

// Cloud overlay on channel A
if(MCIR){
float cloud = CLIP((crow[y][chb].r - 115) / 107, 0, 1);
crow[y][cha] = RGBcomposite((rgb_t){240, 250, 255}, cloud, crow[y][cha], 1);
}

// Precipitation
if(crow[y][chb].r > 191){
float ramp = CLIP((crow[y][chb].r - 191) / 5.0, 0, 1);
crow[y][cha] = RGBcomposite(applyPalette(PrecipPalette, crow[y][chb].r), ramp, crow[y][cha], 1);
}
}
}

return 1;
}

int readRawImage(char *filename, float **prow, int *nrow) {
FILE *fp = fopen(filename, "r");

// Create reader
png_structp png = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
if(!png) return 0;
png_infop info = png_create_info_struct(png);
if(!info) 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 != IMG_WIDTH){
fprintf(stderr, "Raw image must be %ipx wide.\n", IMG_WIDTH);
return 0;
}else if(bit_depth != 8){
fprintf(stderr, "Raw image must have 8 bit color.\n");
return 0;
}else if(color_type != PNG_COLOR_TYPE_GRAY){
fprintf(stderr, "Raw image must be grayscale.\n");
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;
}

typedef struct {
float *prow[3000]; // Row buffers
int nrow; // Number of rows
int chA, chB; // ID of each channel
char name[256]; // Stripped filename
} image_t;

typedef struct {
char *type; // Output image type
char *effects;
int satnum; // The satellite number
char *map; // Path to a map file
char *path; // Output directory
} options_t;

//int ImageOut(char *filename, char *chid, float **prow, int nrow, int width, int offset, char *palette, char *effects, char *mapFile) {
int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, char *chid, char *palette){
char outName[384];
sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid);
FILE *pngfile;

png_text text_ptr[] = {
{PNG_TEXT_COMPRESSION_NONE, "Software", VERSION},
{PNG_TEXT_COMPRESSION_NONE, "Channel", desc, sizeof(desc)},
{PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA satellite image", 20}
};

// Reduce the width of the image to componsate for the missing telemetry
if(opts->effects != NULL && CONTAINS(opts->effects, 't'))
width -= TOTAL_TELE;

// 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);
fprintf(stderr, ERR_PNG_WRITE);
return(0);
}
png_infop info_ptr = png_create_info_struct(png_ptr);
if (!info_ptr) {
png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
fprintf(stderr, ERR_PNG_INFO);
return(0);
}

// 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, text_ptr, 3);

// Channel = 25cm wide
png_set_pHYs(png_ptr, info_ptr, 3636, 3636, PNG_RESOLUTION_METER);

// Init I/O
pngfile = fopen(outName, "wb");
if (!pngfile) {
fprintf(stderr, ERR_FILE_WRITE, outName);
return(1);
}
png_init_io(png_ptr, pngfile);
png_write_info(png_ptr, info_ptr);

// Move prow into crow, crow ~ color rows
rgb_t *crow[img->nrow];
for(int y = 0; y < img->nrow; y++){
crow[y] = (rgb_t *) malloc(sizeof(rgb_t) * IMG_WIDTH);

for(int x = 0; x < IMG_WIDTH; x++){
if(palette == NULL){
crow[y][x].r = crow[y][x].g = crow[y][x].b = CLIP(img->prow[y][x], 0, 255);
}else{
crow[y][x] = applyPalette(palette, img->prow[y][x]);
}
}
}

if(opts->map != NULL && opts->map[0] != '\0'){
if(mapOverlay(opts->map, crow, img->nrow, zenith, strcmp(chid, "MCIR") == 0) == 0){
fprintf(stderr, "Skipping MCIR generation; see above.\n");
return 0;
}
}else if(strcmp(chid, "MCIR") == 0){
fprintf(stderr, "Skipping MCIR generation; no map provided.\n");
return 0;
}

//extern void falsecolor(float vis, float temp, float *r, float *g, float *b);

printf("Writing %s", outName);

int fcimage = strcmp(chid, "False Color") == 0;

// Build RGB image
for (int y = 0; y < img->nrow; y++) {
png_color pix[width];

for (int x = 0; x < width; x++) {
pix[x].red = (int)crow[y][x + offset].r;
pix[x].green = (int)crow[y][x + offset].g;
pix[x].blue = (int)crow[y][x + offset].b;
}
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);
}


Write Preview
Loading…
Cancel
Save