diff --git a/Makefile b/Makefile index 6094223..b63097a 100644 --- a/Makefile +++ b/Makefile @@ -9,9 +9,9 @@ aptdec: $(OBJS) reg.o: reg.c color.o: color.c -main.o: main.c +main.o: main.c media.o: median.c -dsp.o: dsp.c +dsp.o: dsp.c filter.o: filter.c image.o: image.c pngio.o: pngio.c diff --git a/color.c b/color.c index 2a0a397..aa0a8f7 100644 --- a/color.c +++ b/color.c @@ -30,96 +30,96 @@ rgb_t applyPalette(char *palette, int val){ } rgb_t RGBcomposite(rgb_t top, float top_a, rgb_t bottom, float bottom_a){ - return (rgb_t){ - MCOMPOSITE(top.r, top_a, bottom.r, bottom_a), - MCOMPOSITE(top.g, top_a, bottom.g, bottom_a), - MCOMPOSITE(top.b, top_a, bottom.b, bottom_a) - }; + return (rgb_t){ + MCOMPOSITE(top.r, top_a, bottom.r, bottom_a), + MCOMPOSITE(top.g, top_a, bottom.g, bottom_a), + MCOMPOSITE(top.b, top_a, bottom.b, bottom_a) + }; } // The "I totally didn't just steal this from WXtoImg" palette char TempPalette[256*3] = { - "\x45\x0\x8f\x46\x0\x91\x47\x0\x92\x48\x0\x94\x49\x0\x96\x4a\x0\x98\x4b\x0\x9b\x4d\x0\x9d" - "\x4e\x0\xa0\x50\x0\xa2\x51\x0\xa5\x52\x0\xa7\x54\x0\xaa\x56\x0\xae\x57\x0\xb1" - "\x58\x0\xb4\x5a\x0\xb7\x5c\x0\xba\x5e\x0\xbd\x5f\x0\xc0\x61\x0\xc4\x64\x0\xc8" - "\x66\x0\xcb\x68\x0\xce\x69\x0\xd1\x68\x0\xd4\x65\x0\xd7\x63\x0\xda\x61\x0\xdd" - "\x5d\x0\xe1\x5b\x0\xe4\x59\x0\xe6\x56\x0\xe9\x53\x0\xeb\x50\x0\xee\x4d\x0\xf0" - "\x49\x0\xf3\x47\x0\xfc\x43\x0\xfa\x31\x0\xbf\x20\x0\x89\x20\x0\x92\x1e\x0\x95" - "\x1b\x0\x97\x19\x0\x9a\x17\x0\x9c\x15\x0\x9e\x12\x0\xa0\xf\x0\xa3\xf\x2\xa5" - "\xe\x6\xa8\xe\xa\xab\xe\xd\xad\xe\x11\xb1\xd\x15\xb4\xd\x18\xb7\xd\x1c\xba" - "\xb\x21\xbd\xa\x25\xc0\xa\x29\xc3\x9\x2d\xc6\x8\x33\xca\x7\x36\xcd\x7\x3b\xd0" - "\x7\x41\xd3\x5\x45\xd6\x4\x4b\xd9\x4\x50\xdc\x3\x55\xde\x2\x5d\xe2\x1\x61\xe5" - "\x0\x66\xe7\x0\x6c\xea\x0\x72\xec\x0\x78\xee\x0\x7d\xf0\x0\x82\xf3\x0\x8d\xfc" - "\x0\x90\xfa\x0\x71\xbf\x0\x54\x89\x0\x5c\x91\x0\x61\x94\x0\x64\x96\x0\x68\x97" - "\x0\x6d\x99\x0\x71\x9b\x0\x75\x9d\x0\x79\x9f\x0\x7e\xa0\x0\x82\xa2\x0\x87\xa4" - "\x0\x8c\xa6\x0\x92\xaa\x0\x96\xac\x0\x9c\xae\x0\xa2\xb1\x0\xa6\xb3\x0\xaa\xb5" - "\x0\xad\xb7\x0\xb1\xba\x0\xb6\xbe\x0\xba\xc0\x0\xbe\xc2\x0\xc2\xc5\x0\xc6\xc6" - "\x0\xca\xc9\x0\xcc\xca\x0\xcf\xcb\x0\xd2\xcc\x0\xd4\xcc\x0\xd6\xcc\x0\xd9\xcb" - "\x0\xdb\xcb\x0\xde\xcb\x0\xe0\xcb\x0\xe2\xcc\x0\xea\xd2\x0\xea\xcf\x0\xb9\xa4" - "\x0\x8e\x7a\x1\x94\x7c\x4\x97\x79\x7\x99\x75\x9\x9b\x71\xd\x9d\x6b\x10\x9f\x67" - "\x12\xa1\x63\x15\xa3\x5f\x17\xa5\x59\x1a\xa8\x55\x1d\xaa\x50\x20\xac\x4b\x24\xaf\x45" - "\x28\xb2\x41\x2b\xb5\x3b\x2e\xb8\x35\x31\xba\x30\x34\xbd\x2b\x39\xbf\x24\x3f\xc1\x17" - "\x49\xc5\x8\x4f\xc8\x1\x4f\xca\x0\x4e\xcd\x0\x4e\xcf\x0\x4f\xd2\x0\x54\xd5\x0" - "\x5d\xd8\x0\x68\xdb\x0\x6e\xdd\x0\x74\xdf\x0\x7a\xe2\x0\x7f\xe4\x0\x85\xe7\x0" - "\x8b\xe9\x0\x8f\xeb\x0\x9b\xf3\x0\x9e\xf2\x0\x7e\xbb\x0\x60\x8a\x0\x68\x92\x0" - "\x6d\x95\x0\x71\x96\x0\x75\x98\x0\x7b\x9a\x0\x7f\x9d\x0\x83\x9f\x0\x87\xa1\x0" - "\x8c\xa2\x0\x8f\xa5\x0\x92\xa7\x0\x96\xa9\x0\x9a\xad\x0\x9d\xb0\x0\xa1\xb2\x0" - "\xa5\xb5\x0\xa9\xb7\x0\xad\xba\x0\xb2\xbd\x0\xb6\xbf\x0\xbb\xc3\x0\xbf\xc6\x0" - "\xc3\xc8\x0\xc8\xcb\x0\xcc\xce\x0\xd0\xd1\x0\xd3\xd2\x0\xd5\xd4\x0\xd9\xd4\x0" - "\xdc\xd4\x0\xde\xd5\x0\xe1\xd5\x0\xe3\xd5\x0\xe6\xd4\x0\xe8\xd1\x0\xea\xce\x0" - "\xf2\xcf\x0\xf2\xca\x0\xbb\x99\x0\x8a\x6e\x0\x92\x72\x0\x95\x72\x0\x97\x71\x0" - "\x9a\x70\x0\x9c\x6e\x0\x9e\x6d\x0\xa0\x6b\x0\xa3\x6a\x0\xa5\x68\x0\xa8\x67\x0" - "\xab\x66\x0\xae\x65\x0\xb2\x63\x0\xb4\x61\x0\xb7\x5f\x0\xba\x5d\x0\xbd\x5c\x0" - "\xc0\x59\x0\xc3\x57\x0\xc6\x54\x0\xca\x50\x0\xcd\x4d\x0\xd0\x4a\x0\xd3\x47\x0" - "\xd6\x43\x0\xd9\x40\x0\xdc\x3d\x0\xde\x39\x0\xe2\x33\x0\xe5\x2f\x0\xe7\x2c\x0" - "\xea\x28\x0\xec\x23\x0\xef\x1f\x0\xf1\x1a\x0\xf3\x14\x0\xfb\xf\x0\xfa\xd\x0" - "\xc1\x5\x0\x8e\x0\x0\x97\x0\x0\x9b\x0\x0\x9e\x0\x0\xa1\x0\x0\xa5\x0\x0" - "\xa9\x0\x0\xad\x0\x0\xb1\x0\x0\xb6\x0\x0\xba\x0\x0\xbd\x0\x0\xc2\x0\x0" - "\xc8\x0\x0\xcc\x0\x0\xcc\x0\x0" + "\x45\x0\x8f\x46\x0\x91\x47\x0\x92\x48\x0\x94\x49\x0\x96\x4a\x0\x98\x4b\x0\x9b\x4d\x0\x9d" + "\x4e\x0\xa0\x50\x0\xa2\x51\x0\xa5\x52\x0\xa7\x54\x0\xaa\x56\x0\xae\x57\x0\xb1" + "\x58\x0\xb4\x5a\x0\xb7\x5c\x0\xba\x5e\x0\xbd\x5f\x0\xc0\x61\x0\xc4\x64\x0\xc8" + "\x66\x0\xcb\x68\x0\xce\x69\x0\xd1\x68\x0\xd4\x65\x0\xd7\x63\x0\xda\x61\x0\xdd" + "\x5d\x0\xe1\x5b\x0\xe4\x59\x0\xe6\x56\x0\xe9\x53\x0\xeb\x50\x0\xee\x4d\x0\xf0" + "\x49\x0\xf3\x47\x0\xfc\x43\x0\xfa\x31\x0\xbf\x20\x0\x89\x20\x0\x92\x1e\x0\x95" + "\x1b\x0\x97\x19\x0\x9a\x17\x0\x9c\x15\x0\x9e\x12\x0\xa0\xf\x0\xa3\xf\x2\xa5" + "\xe\x6\xa8\xe\xa\xab\xe\xd\xad\xe\x11\xb1\xd\x15\xb4\xd\x18\xb7\xd\x1c\xba" + "\xb\x21\xbd\xa\x25\xc0\xa\x29\xc3\x9\x2d\xc6\x8\x33\xca\x7\x36\xcd\x7\x3b\xd0" + "\x7\x41\xd3\x5\x45\xd6\x4\x4b\xd9\x4\x50\xdc\x3\x55\xde\x2\x5d\xe2\x1\x61\xe5" + "\x0\x66\xe7\x0\x6c\xea\x0\x72\xec\x0\x78\xee\x0\x7d\xf0\x0\x82\xf3\x0\x8d\xfc" + "\x0\x90\xfa\x0\x71\xbf\x0\x54\x89\x0\x5c\x91\x0\x61\x94\x0\x64\x96\x0\x68\x97" + "\x0\x6d\x99\x0\x71\x9b\x0\x75\x9d\x0\x79\x9f\x0\x7e\xa0\x0\x82\xa2\x0\x87\xa4" + "\x0\x8c\xa6\x0\x92\xaa\x0\x96\xac\x0\x9c\xae\x0\xa2\xb1\x0\xa6\xb3\x0\xaa\xb5" + "\x0\xad\xb7\x0\xb1\xba\x0\xb6\xbe\x0\xba\xc0\x0\xbe\xc2\x0\xc2\xc5\x0\xc6\xc6" + "\x0\xca\xc9\x0\xcc\xca\x0\xcf\xcb\x0\xd2\xcc\x0\xd4\xcc\x0\xd6\xcc\x0\xd9\xcb" + "\x0\xdb\xcb\x0\xde\xcb\x0\xe0\xcb\x0\xe2\xcc\x0\xea\xd2\x0\xea\xcf\x0\xb9\xa4" + "\x0\x8e\x7a\x1\x94\x7c\x4\x97\x79\x7\x99\x75\x9\x9b\x71\xd\x9d\x6b\x10\x9f\x67" + "\x12\xa1\x63\x15\xa3\x5f\x17\xa5\x59\x1a\xa8\x55\x1d\xaa\x50\x20\xac\x4b\x24\xaf\x45" + "\x28\xb2\x41\x2b\xb5\x3b\x2e\xb8\x35\x31\xba\x30\x34\xbd\x2b\x39\xbf\x24\x3f\xc1\x17" + "\x49\xc5\x8\x4f\xc8\x1\x4f\xca\x0\x4e\xcd\x0\x4e\xcf\x0\x4f\xd2\x0\x54\xd5\x0" + "\x5d\xd8\x0\x68\xdb\x0\x6e\xdd\x0\x74\xdf\x0\x7a\xe2\x0\x7f\xe4\x0\x85\xe7\x0" + "\x8b\xe9\x0\x8f\xeb\x0\x9b\xf3\x0\x9e\xf2\x0\x7e\xbb\x0\x60\x8a\x0\x68\x92\x0" + "\x6d\x95\x0\x71\x96\x0\x75\x98\x0\x7b\x9a\x0\x7f\x9d\x0\x83\x9f\x0\x87\xa1\x0" + "\x8c\xa2\x0\x8f\xa5\x0\x92\xa7\x0\x96\xa9\x0\x9a\xad\x0\x9d\xb0\x0\xa1\xb2\x0" + "\xa5\xb5\x0\xa9\xb7\x0\xad\xba\x0\xb2\xbd\x0\xb6\xbf\x0\xbb\xc3\x0\xbf\xc6\x0" + "\xc3\xc8\x0\xc8\xcb\x0\xcc\xce\x0\xd0\xd1\x0\xd3\xd2\x0\xd5\xd4\x0\xd9\xd4\x0" + "\xdc\xd4\x0\xde\xd5\x0\xe1\xd5\x0\xe3\xd5\x0\xe6\xd4\x0\xe8\xd1\x0\xea\xce\x0" + "\xf2\xcf\x0\xf2\xca\x0\xbb\x99\x0\x8a\x6e\x0\x92\x72\x0\x95\x72\x0\x97\x71\x0" + "\x9a\x70\x0\x9c\x6e\x0\x9e\x6d\x0\xa0\x6b\x0\xa3\x6a\x0\xa5\x68\x0\xa8\x67\x0" + "\xab\x66\x0\xae\x65\x0\xb2\x63\x0\xb4\x61\x0\xb7\x5f\x0\xba\x5d\x0\xbd\x5c\x0" + "\xc0\x59\x0\xc3\x57\x0\xc6\x54\x0\xca\x50\x0\xcd\x4d\x0\xd0\x4a\x0\xd3\x47\x0" + "\xd6\x43\x0\xd9\x40\x0\xdc\x3d\x0\xde\x39\x0\xe2\x33\x0\xe5\x2f\x0\xe7\x2c\x0" + "\xea\x28\x0\xec\x23\x0\xef\x1f\x0\xf1\x1a\x0\xf3\x14\x0\xfb\xf\x0\xfa\xd\x0" + "\xc1\x5\x0\x8e\x0\x0\x97\x0\x0\x9b\x0\x0\x9e\x0\x0\xa1\x0\x0\xa5\x0\x0" + "\xa9\x0\x0\xad\x0\x0\xb1\x0\x0\xb6\x0\x0\xba\x0\x0\xbd\x0\x0\xc2\x0\x0" + "\xc8\x0\x0\xcc\x0\x0\xcc\x0\x0" }; char PrecipPalette[256*3] = { - "\xe0\x98\x8\xec\x84\x10\xf5\x70\x1b\xfc\x5c\x29\xff\x49\x38\xff\x37\x4a" - "\xfb\x28\x5d\xf5\x1a\x71\xeb\xf\x85\xdf\x8\x99\xd0\x3\xad\xc0\x1\xbf" - "\xad\x3\xd0\x9a\x7\xdf\x86\xf\xeb\x72\x1a\xf5\x5e\x27\xfb\x4b\x37\xfe" - "\x39\x48\xff\x29\x5b\xfc\x1b\x6f\xf5\x10\x83\xec\x8\x98\xe0\x3\xab\xd2" - "\x1\xbe\xc1\x2\xcf\xaf\x7\xdd\x9c\xe\xea\x88\x19\xf4\x73\x26\xfb\x5f" - "\x35\xfe\x4c\x47\xff\x3a\x5a\xfc\x2a\x6d\xf6\x1d\x82\xed\x11\x96\xe1\x9" - "\xa9\xd3\x3\xbc\xc3\x1\xcd\xb1\x2\xdc\x9e\x6\xe9\x89\xe\xf3\x75\x18" - "\xfa\x61\x25\xfe\x4e\x34\xff\x3c\x45\xfc\x2c\x58\xf7\x1e\x6c\xee\x12\x80" - "\xe3\x9\x94\xd5\x4\xa8\xc4\x1\xba\xb2\x2\xcc\x9f\x6\xdb\x8b\xd\xe8" - "\x77\x17\xf2\x63\x23\xfa\x50\x33\xfe\x3d\x44\xff\x2d\x56\xfd\x1f\x6a\xf7" - "\x13\x7e\xef\xa\x92\xe4\x4\xa6\xd6\x1\xb9\xc6\x2\xca\xb4\x5\xda\xa1" - "\xc\xe7\x8d\x16\xf1\x79\x22\xf9\x65\x31\xfe\x51\x42\xff\x3f\x54\xfd\x2f" - "\x68\xf8\x20\x7c\xf0\x14\x90\xe5\xb\xa4\xd7\x4\xb7\xc7\x1\xc9\xb6\x1" - "\xd8\xa3\x5\xe6\x8f\xb\xf1\x7b\x15\xf9\x66\x21\xfd\x53\x30\xff\x41\x40" - "\xfd\x30\x53\xf9\x21\x66\xf1\x15\x7a\xe6\xb\x8f\xd9\x5\xa3\xc9\x1\xb6" - "\xb7\x1\xc7\xa5\x4\xd7\x91\xb\xe5\x7c\x14\xf0\x68\x20\xf8\x55\x2e\xfd" - "\x42\x3f\xff\x31\x51\xfe\x22\x64\xf9\x16\x78\xf2\xc\x8d\xe7\x5\xa1\xda" - "\x2\xb4\xca\x1\xc6\xb9\x4\xd6\xa6\xa\xe4\x93\x13\xef\x7e\x1f\xf7\x6a" - "\x2d\xfd\x56\x3d\xff\x44\x4f\xfe\x33\x63\xfa\x24\x77\xf2\x17\x8b\xe8\xd" - "\x9f\xdb\x6\xb2\xcc\x2\xc4\xbb\x1\xd4\xa8\x4\xe2\x94\x9\xee\x80\x12" - "\xf7\x6c\x1d\xfc\x58\x2b\xff\x45\x3c\xfe\x34\x4e\xfa\x25\x61\xf3\x18\x75" - "\xe9\xe\x89\xdc\x6\x9d\xcd\x2\xb1\xbc\x1\xc3\xaa\x3\xd3\x96\x9\xe1" - "\x82\x11\xed\x6e\x1c\xf6\x5a\x2a\xfc\x47\x3a\xff\x36\x4c\xfe\x26\x5f\xfb" - "\x19\x73\xf4\xf\x87\xea\x7\x9b\xde\x2\xaf\xcf\x1\xc1\xbe\x3\xd2\xac" - "\x8\xe0\x98\x10\xec\x84\x1b\xf5\x6f\x29\xfc\x5c\x39\xff\x49\x4a\xff\x37" - "\x5d\xfb\x28\x71\xf5\x1a\x85\xeb\xf\x9a\xdf\x7\xad\xd0\x3\xbf\xc0\x1" - "\xd0\xad\x3\xdf\x9a\x7\xeb\x86\xf\xf5\x71\x1a\xfb\x5d\x27\xff\x4a\x37" - "\xff\x39\x49\xfc\x29\x5c\xf5\x1b\x6f\xec\x10\x84\xe0\x8\x98\xd2\x3\xab" - "\xc1\x1\xbe\xaf\x2\xcf\x9c\x7\xde\x87\xe\xea\x73\x19\xf4\x5f\x26\xfb" - "\x4c\x36\xfe\x3a\x47\xff\x2a\x5a\xfc\x1c\x6e\xf6\x11\x82\xed\x9\x96\xe1" - "\x3\xaa\xd3\x1\xbc\xc3\x2\xcd\xb1\x6\xdc\x9d\xe\xe9\x89\x18\xf3\x75" - "\x25\xfa\x61\x34\xfe\x4e\x45\xff\x3c\x58\xfc\x2c\x6c\xf7\x1e\x80\xee\x12" - "\x94\xe2\x9\xa8\xd4\x4\xbb\xc4\x1\xcc\xb2\x2\xdb\x9f\x6\xe8\x8b\xd" - "\xf2\x77\x17\xfa\x63\x24\xfe\x4f\x33\xff\x3d\x44\xfd\x2d\x56\xf7\x1f\x6a" - "\xef\x13\x7e\xe4\xa\x92\xd6\x4\xa6\xc6\x1\xb9\xb4\x2\xca\xa1\x5\xda" - "\x8d\xc\xe7\x79\x16\xf2\x64\x22\xf9\x51\x31\xfe\x3f\x42\xff\x2e\x55\xfd" - "\x20\x68\xf8\x14\x7c\xf0\xb\x91\xe5\x4\xa4\xd7\x1\xb7\xc7\x1\xc9\xb6" - "\x5\xd9\xa3\xb\xe6\x8f\x15\xf1\x7a\x21\xf9\x66\x30\xfd\x53\x41\xff\x40" - "\x53\xfd\x30\x66\xf9\x21\x7b\xf1\x15\x8f\xe6\xb\xa3\xd8\x5\xb6\xc9\x1" - "\xc7\xb7\x1\xd7\xa4\x4\xe5\x91\xb\xf0\x7c\x14\xf8\x68\x20\xfd\x54\x2e" - "\xff\x42\x3f\xfe\x31\x51\xf9\x22\x65\xf1\x16\x79\xe7\xc\x8d\xda\x5\xa1" - "\xca\x2\xb4\xb9\x1\xc6\xa6\x4\xd6\x92\xa\xe4\x7e\x13\xef\x6a\x1f\xf7" - "\x56\x2d\xfd\x44\x3d\xff\x33\x4f\xfe\x24\x63\xfa" + "\xe0\x98\x8\xec\x84\x10\xf5\x70\x1b\xfc\x5c\x29\xff\x49\x38\xff\x37\x4a" + "\xfb\x28\x5d\xf5\x1a\x71\xeb\xf\x85\xdf\x8\x99\xd0\x3\xad\xc0\x1\xbf" + "\xad\x3\xd0\x9a\x7\xdf\x86\xf\xeb\x72\x1a\xf5\x5e\x27\xfb\x4b\x37\xfe" + "\x39\x48\xff\x29\x5b\xfc\x1b\x6f\xf5\x10\x83\xec\x8\x98\xe0\x3\xab\xd2" + "\x1\xbe\xc1\x2\xcf\xaf\x7\xdd\x9c\xe\xea\x88\x19\xf4\x73\x26\xfb\x5f" + "\x35\xfe\x4c\x47\xff\x3a\x5a\xfc\x2a\x6d\xf6\x1d\x82\xed\x11\x96\xe1\x9" + "\xa9\xd3\x3\xbc\xc3\x1\xcd\xb1\x2\xdc\x9e\x6\xe9\x89\xe\xf3\x75\x18" + "\xfa\x61\x25\xfe\x4e\x34\xff\x3c\x45\xfc\x2c\x58\xf7\x1e\x6c\xee\x12\x80" + "\xe3\x9\x94\xd5\x4\xa8\xc4\x1\xba\xb2\x2\xcc\x9f\x6\xdb\x8b\xd\xe8" + "\x77\x17\xf2\x63\x23\xfa\x50\x33\xfe\x3d\x44\xff\x2d\x56\xfd\x1f\x6a\xf7" + "\x13\x7e\xef\xa\x92\xe4\x4\xa6\xd6\x1\xb9\xc6\x2\xca\xb4\x5\xda\xa1" + "\xc\xe7\x8d\x16\xf1\x79\x22\xf9\x65\x31\xfe\x51\x42\xff\x3f\x54\xfd\x2f" + "\x68\xf8\x20\x7c\xf0\x14\x90\xe5\xb\xa4\xd7\x4\xb7\xc7\x1\xc9\xb6\x1" + "\xd8\xa3\x5\xe6\x8f\xb\xf1\x7b\x15\xf9\x66\x21\xfd\x53\x30\xff\x41\x40" + "\xfd\x30\x53\xf9\x21\x66\xf1\x15\x7a\xe6\xb\x8f\xd9\x5\xa3\xc9\x1\xb6" + "\xb7\x1\xc7\xa5\x4\xd7\x91\xb\xe5\x7c\x14\xf0\x68\x20\xf8\x55\x2e\xfd" + "\x42\x3f\xff\x31\x51\xfe\x22\x64\xf9\x16\x78\xf2\xc\x8d\xe7\x5\xa1\xda" + "\x2\xb4\xca\x1\xc6\xb9\x4\xd6\xa6\xa\xe4\x93\x13\xef\x7e\x1f\xf7\x6a" + "\x2d\xfd\x56\x3d\xff\x44\x4f\xfe\x33\x63\xfa\x24\x77\xf2\x17\x8b\xe8\xd" + "\x9f\xdb\x6\xb2\xcc\x2\xc4\xbb\x1\xd4\xa8\x4\xe2\x94\x9\xee\x80\x12" + "\xf7\x6c\x1d\xfc\x58\x2b\xff\x45\x3c\xfe\x34\x4e\xfa\x25\x61\xf3\x18\x75" + "\xe9\xe\x89\xdc\x6\x9d\xcd\x2\xb1\xbc\x1\xc3\xaa\x3\xd3\x96\x9\xe1" + "\x82\x11\xed\x6e\x1c\xf6\x5a\x2a\xfc\x47\x3a\xff\x36\x4c\xfe\x26\x5f\xfb" + "\x19\x73\xf4\xf\x87\xea\x7\x9b\xde\x2\xaf\xcf\x1\xc1\xbe\x3\xd2\xac" + "\x8\xe0\x98\x10\xec\x84\x1b\xf5\x6f\x29\xfc\x5c\x39\xff\x49\x4a\xff\x37" + "\x5d\xfb\x28\x71\xf5\x1a\x85\xeb\xf\x9a\xdf\x7\xad\xd0\x3\xbf\xc0\x1" + "\xd0\xad\x3\xdf\x9a\x7\xeb\x86\xf\xf5\x71\x1a\xfb\x5d\x27\xff\x4a\x37" + "\xff\x39\x49\xfc\x29\x5c\xf5\x1b\x6f\xec\x10\x84\xe0\x8\x98\xd2\x3\xab" + "\xc1\x1\xbe\xaf\x2\xcf\x9c\x7\xde\x87\xe\xea\x73\x19\xf4\x5f\x26\xfb" + "\x4c\x36\xfe\x3a\x47\xff\x2a\x5a\xfc\x1c\x6e\xf6\x11\x82\xed\x9\x96\xe1" + "\x3\xaa\xd3\x1\xbc\xc3\x2\xcd\xb1\x6\xdc\x9d\xe\xe9\x89\x18\xf3\x75" + "\x25\xfa\x61\x34\xfe\x4e\x45\xff\x3c\x58\xfc\x2c\x6c\xf7\x1e\x80\xee\x12" + "\x94\xe2\x9\xa8\xd4\x4\xbb\xc4\x1\xcc\xb2\x2\xdb\x9f\x6\xe8\x8b\xd" + "\xf2\x77\x17\xfa\x63\x24\xfe\x4f\x33\xff\x3d\x44\xfd\x2d\x56\xf7\x1f\x6a" + "\xef\x13\x7e\xe4\xa\x92\xd6\x4\xa6\xc6\x1\xb9\xb4\x2\xca\xa1\x5\xda" + "\x8d\xc\xe7\x79\x16\xf2\x64\x22\xf9\x51\x31\xfe\x3f\x42\xff\x2e\x55\xfd" + "\x20\x68\xf8\x14\x7c\xf0\xb\x91\xe5\x4\xa4\xd7\x1\xb7\xc7\x1\xc9\xb6" + "\x5\xd9\xa3\xb\xe6\x8f\x15\xf1\x7a\x21\xf9\x66\x30\xfd\x53\x41\xff\x40" + "\x53\xfd\x30\x66\xf9\x21\x7b\xf1\x15\x8f\xe6\xb\xa3\xd8\x5\xb6\xc9\x1" + "\xc7\xb7\x1\xd7\xa4\x4\xe5\x91\xb\xf0\x7c\x14\xf8\x68\x20\xfd\x54\x2e" + "\xff\x42\x3f\xfe\x31\x51\xf9\x22\x65\xf1\x16\x79\xe7\xc\x8d\xda\x5\xa1" + "\xca\x2\xb4\xb9\x1\xc6\xa6\x4\xd6\x92\xa\xe4\x7e\x13\xef\x6a\x1f\xf7" + "\x56\x2d\xfd\x44\x3d\xff\x33\x4f\xfe\x24\x63\xfa" }; \ No newline at end of file diff --git a/dsp.c b/dsp.c index 804d11c..d32fcae 100755 --- a/dsp.c +++ b/dsp.c @@ -74,25 +74,25 @@ static inline double Phase(double I, double Q) { if(I == 0.0 && Q == 0.0) return 0.0; - if (Q < 0) { + if (Q < 0) { s = -1; Q = -Q; - } else { + } else { s = 1; - } - - if (I >= 0) { - r = (I - Q) / (I + Q); - angle = 0.25 - 0.25 * r; - } else { - r = (I + Q) / (Q - I); - angle = 0.75 - 0.25 * r; - } - - if(s > 0){ - return angle; - }else{ - return -angle; + } + + if (I >= 0) { + r = (I - Q) / (I + Q); + angle = 0.25 - 0.25 * r; + } else { + r = (I + Q) / (Q - I); + angle = 0.75 - 0.25 * r; + } + + if(s > 0){ + return angle; + }else{ + return -angle; } } @@ -102,44 +102,44 @@ static inline double Phase(double I, double Q) { */ static double pll(double I, double Q) { // PLL coefficient - static double PhaseOsc = 0.0; - double Io, Qo; - double Ip, Qp; - double DPhi; + static double PhaseOsc = 0.0; + double Io, Qo; + double Ip, Qp; + double DPhi; // Quadrature oscillator / reference - Io = cos(PhaseOsc); - Qo = sin(PhaseOsc); + Io = cos(PhaseOsc); + Qo = sin(PhaseOsc); // Phase detector - Ip = I * Io + Q * Qo; - Qp = Q * Io - I * Qo; - DPhi = Phase(Ip, Qp); + Ip = I * Io + Q * Qo; + Qp = Q * Io - I * Qo; + DPhi = Phase(Ip, Qp); // Loop filter - PhaseOsc += 2.0 * M_PI * (K1 * DPhi + FreqOsc); - if (PhaseOsc > M_PI) + PhaseOsc += 2.0 * M_PI * (K1 * DPhi + FreqOsc); + if (PhaseOsc > M_PI) PhaseOsc -= 2.0 * M_PI; - if (PhaseOsc <= -M_PI) + if (PhaseOsc <= -M_PI) PhaseOsc += 2.0 * M_PI; - FreqOsc += K2 * DPhi; - if (FreqOsc > ((Fc + DFc) / Fe)) + FreqOsc += K2 * DPhi; + if (FreqOsc > ((Fc + DFc) / Fe)) FreqOsc = (Fc + DFc) / Fe; - if (FreqOsc < ((Fc - DFc) / Fe)) + if (FreqOsc < ((Fc - DFc) / Fe)) FreqOsc = (Fc - DFc) / Fe; - return Ip; + return Ip; } // Convert samples into pixels static int getamp(double *ampbuff, int count) { - static float inbuff[BLKIN]; - static int idxin = 0; - static int nin = 0; + static float inbuff[BLKIN]; + static int idxin = 0; + static int nin = 0; - for (int n = 0; n < count; n++) { - double I, Q; + for (int n = 0; n < count; n++) { + double I, Q; // Get some more samples when needed if (nin < IQFilterLen * 2 + 2) { @@ -147,11 +147,11 @@ static int getamp(double *ampbuff, int count) { int res; memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); idxin = 0; - + // Read some samples res = getsample(&(inbuff[nin]), BLKIN - nin); nin += res; - + // Make sure there is enough samples to continue if (nin < IQFilterLen * 2 + 2) return n; @@ -164,25 +164,25 @@ static int getamp(double *ampbuff, int count) { // Increment current sample idxin++; nin--; - } + } - return count; + return count; } // Sub-pixel offsetting + FIR compensation int getpixelv(float *pvbuff, int count) { // Amplitude buffer - static double ampbuff[BLKAMP]; - static int nam = 0; - static int idxam = 0; + static double ampbuff[BLKAMP]; + static int nam = 0; + static int idxam = 0; - double mult; + double mult; // Gaussian resampling factor - mult = (double) Fi / Fe * FreqLine; - int m = RSFilterLen / mult + 1; + mult = (double) Fi / Fe * FreqLine; + int m = RSFilterLen / mult + 1; - for (int n = 0; n < count; n++) { + for (int n = 0; n < count; n++) { int shift; if (nam < m) { @@ -203,53 +203,53 @@ int getpixelv(float *pvbuff, int count) { idxam += shift; nam -= shift; - } - - return count; + } + + return count; } // Get an entire row of pixels, aligned with sync markers // FIXME: skips noisy lines with no findable sync marker 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 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; + double corr, ecorr, lcorr; + int res; // Move the row buffer into the the image buffer - if (npv > 0) + if (npv > 0) memmove(pixelv, pixels, npv * sizeof(float)); // Get the sync line - if (npv < SyncFilterLen + 2) { + if (npv < SyncFilterLen + 2) { res = getpixelv(&(pixelv[npv]), SyncFilterLen + 2 - npv); npv += res; if (npv < SyncFilterLen + 2) return 0; - } + } - // Calculate the frequency offset - ecorr = fir(pixelv, Sync, SyncFilterLen); - corr = fir(&pixelv[1], Sync, SyncFilterLen - 1); - lcorr = fir(&pixelv[2], Sync, SyncFilterLen - 2); - FreqLine = 1.0+((ecorr-lcorr) / corr / PixelLine / 4.0); + // Calculate the frequency offset + ecorr = fir(pixelv, Sync, SyncFilterLen); + 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 if(fabs(lcorr - ecorr) < minDoppler && fabs(lcorr - ecorr) > 2){ minDoppler = fabs(lcorr - ecorr); *zenith = nrow; } - + // The point in which the pixel offset is recalculated - if (corr < 0.75 * max) { + if (corr < 0.75 * max) { synced = 0; FreqLine = 1.0; - } - max = corr; + } + max = corr; if (synced < 8) { int mshift; @@ -285,20 +285,20 @@ int getpixelrow(float *pixelv, int nrow, int *zenith) { } // Get the rest of this row - if (npv < PixelLine) { + if (npv < PixelLine) { res = getpixelv(&(pixelv[npv]), PixelLine - npv); npv += res; if (npv < PixelLine) return 0; - } + } // Move the sync lines into the output buffer with the calculated offset - if (npv == PixelLine) { + if (npv == PixelLine) { npv = 0; - } else { + } else { memmove(pixels, &(pixelv[PixelLine]), (npv - PixelLine) * sizeof(float)); npv -= PixelLine; - } + } - return 1; + return 1; } \ No newline at end of file diff --git a/filter.c b/filter.c index 9e27590..b99db14 100755 --- a/filter.c +++ b/filter.c @@ -23,45 +23,45 @@ // Finite impulse response float fir(float *buff, const float *coeff, const int len) { - double r; + double r; - r = 0.0; - for (int i = 0; i < len; i++) { - r += buff[i] * coeff[i]; - } - return r; + r = 0.0; + for (int i = 0; i < len; i++) { + r += buff[i] * coeff[i]; + } + return r; } /* IQ finite impulse response * Turn samples into a single IQ sample */ void iqfir(float *buff, const float *coeff, const int len, double *I, double *Q) { - double i = 0.0, q = 0.0; + double i = 0.0, q = 0.0; - for (int k = 0; k < len; k++) { - q += buff[2*k] * coeff[k]; - i += buff[2*k]; - } + for (int k = 0; k < len; k++) { + q += buff[2*k] * coeff[k]; + i += buff[2*k]; + } - i = buff[len-1] - (i / len); - *I = i, *Q = q; + i = buff[len-1] - (i / len); + *I = i, *Q = q; } /* Gaussian finite impulse responce compensation * https://www.recordingblogs.com/wiki/gaussian-window */ float rsfir(double *buff, const float *coeff, const int len, const double offset, const double delta) { - double out; + double out; - out = 0.0; - double n = offset; - for (int i = 0; i < (len-1)/delta-1; n += delta, i++) { - int k; - double alpha; + out = 0.0; + double n = offset; + for (int i = 0; i < (len-1)/delta-1; n += delta, i++) { + int k; + double alpha; - k = (int)floor(n); - alpha = n - k; - out += buff[i] * (coeff[k] * (1.0 - alpha) + coeff[k + 1] * alpha); - } - return out; + k = (int)floor(n); + alpha = n - k; + out += buff[i] * (coeff[k] * (1.0 - alpha) + coeff[k + 1] * alpha); + } + return out; } diff --git a/image.c b/image.c index 3f709bf..dda8d4e 100644 --- a/image.c +++ b/image.c @@ -35,22 +35,22 @@ extern void polyreg(const int m, const int n, const double x[], const double y[] // Compute regression 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 }; + // { 0.106, 0.215, 0.324, 0.433, 0.542, 0.652, 0.78, 0.87, 0.0 } + const double y[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; - polyreg(REGORDER, 9, x, y, rgpr->cf); + polyreg(REGORDER, 9, x, y, rgpr->cf); } // Convert a value to 0-255 based off the provided regression curve static double rgcal(float x, rgparam_t *rgpr) { - double y, p; + double y, p; int i; - - for (i = 0, y = 0.0, p = 1.0; i < REGORDER + 1; i++) { + + for (i = 0, y = 0.0, p = 1.0; i < REGORDER + 1; i++) { y += rgpr->cf[i] * p; p = p * x; - } - return y; + } + return y; } static double tele[16]; @@ -119,38 +119,38 @@ double teleNoise(double wedges[16]){ double noise = 0; for(int i = 0; i < 9; i++) noise += fabs(wedges[i] - pattern[i]); - + return noise; } // Get telemetry data for thermal calibration/equalization int calibrate(float **prow, int nrow, int offset, int width) { double teleline[MAX_HEIGHT] = { 0.0 }; - double wedge[16]; - rgparam_t regr[MAX_HEIGHT/FRAME_LEN + 1]; - int telestart, mtelestart = 0; + double wedge[16]; + rgparam_t regr[MAX_HEIGHT/FRAME_LEN + 1]; + int telestart, mtelestart = 0; int channel = -1; // The minimum rows required to decode a full frame - if (nrow < 192) { + if (nrow < 192) { fprintf(stderr, "Telemetry decoding error, not enough rows\n"); return 0; - } + } // Calculate average of a row of telemetry - for (int y = 0; y < nrow; y++) { + for (int y = 0; y < nrow; y++) { for (int x = 3; x < 43; x++) teleline[y] += prow[y][x + offset + width]; - + teleline[y] /= 40.0; - } + } /* Wedge 7 is white and 8 is black, this will have the largest * difference in brightness, this can be used to find the current * position within the telemetry. */ float max = 0.0f; - for (int n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) { + for (int n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) { float df; // (sum 4px below) / (sum 4px above) @@ -162,20 +162,20 @@ int calibrate(float **prow, int nrow, int offset, int width) { mtelestart = n; max = df; } - } + } - telestart = (mtelestart - 64) % FRAME_LEN; + telestart = (mtelestart - 64) % FRAME_LEN; - // Make sure that theres at least one full frame in the image - if (nrow < telestart + FRAME_LEN) { + // Make sure that theres at least one full frame in the image + if (nrow < telestart + FRAME_LEN) { fprintf(stderr, "Telemetry decoding error, not enough rows\n"); return 0; - } + } // Find the least noisy frame double minNoise = -1; int bestFrame = -1; - for (int n = telestart, k = 0; n < nrow - FRAME_LEN; n += FRAME_LEN, k++) { + for (int n = telestart, k = 0; n < nrow - FRAME_LEN; n += FRAME_LEN, k++) { int j; for (j = 0; j < 16; j++) { @@ -228,7 +228,7 @@ int calibrate(float **prow, int nrow, int offset, int width) { } Cs = rgcal(Cs / i, ®r[k]); } - } + } if(bestFrame == -1){ fprintf(stderr, "Something has gone very wrong, please file a bug report."); @@ -237,8 +237,8 @@ int calibrate(float **prow, int nrow, int offset, int width) { calibrateImage(prow, nrow, offset, width, regr[bestFrame]); - return channel + 1; - + return channel + 1; + } void distrib(options_t *opts, image_t *img, char *chid) { @@ -261,7 +261,7 @@ void distrib(options_t *opts, image_t *img, char *chid) { 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); @@ -275,7 +275,7 @@ void distrib(options_t *opts, image_t *img, char *chid) { 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); } @@ -319,21 +319,21 @@ void flipImage(image_t *img, int width, int offset){ #include "satcal.h" typedef struct { - double Nbb; - double Cs; - double Cb; - int ch; + double Nbb; + double Cs; + double Cb; + int ch; } tempparam_t; // IR channel temperature compensation static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) { - double Tbb, T[4]; - double C; + double Tbb, T[4]; + double C; - tpr->ch = ch - 4; + tpr->ch = ch - 4; - // Compute equivalent T blackbody temperature - for (int n = 0; n < 4; n++) { + // Compute equivalent T blackbody temperature + for (int n = 0; n < 4; n++) { float d0, d1, d2; C = t[9 + n] * 4.0; @@ -344,57 +344,57 @@ static void tempcomp(double t[16], int ch, int satnum, tempparam_t *tpr) { T[n] += d1 * C; C *= C; T[n] += d2 * C; - } - Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0; - Tbb = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb; + } + Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0; + Tbb = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb; - // Compute blackbody radiance temperature - C = satcal[satnum].rad[tpr->ch].vc; - tpr->Nbb = c1 * C * C * C / (expm1(c2 * C / Tbb)); + // Compute blackbody radiance temperature + C = satcal[satnum].rad[tpr->ch].vc; + tpr->Nbb = c1 * C * C * C / (expm1(c2 * C / Tbb)); - // Store blackbody count and space - tpr->Cs = Cs * 4.0; - tpr->Cb = t[14] * 4.0; + // Store blackbody count and space + tpr->Cs = Cs * 4.0; + tpr->Cb = t[14] * 4.0; } // IR channel temperature calibration static double tempcal(float Ce, int satnum, tempparam_t * rgpr) { - double Nl, Nc, Ns, Ne; - double T, vc; + double Nl, Nc, Ns, Ne; + double T, vc; - Ns = satcal[satnum].cor[rgpr->ch].Ns; - Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb); - Nc = satcal[satnum].cor[rgpr->ch].b[0] + + Ns = satcal[satnum].cor[rgpr->ch].Ns; + Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb); + Nc = satcal[satnum].cor[rgpr->ch].b[0] + satcal[satnum].cor[rgpr->ch].b[1] * Nl + satcal[satnum].cor[rgpr->ch].b[2] * Nl * Nl; - Ne = Nl + Nc; + Ne = Nl + Nc; - vc = satcal[satnum].rad[rgpr->ch].vc; - T = c2 * vc / log(c1 * vc * vc * vc / Ne + 1.0); - T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B; + vc = satcal[satnum].rad[rgpr->ch].vc; + T = c2 * vc / log(c1 * vc * vc * vc / Ne + 1.0); + T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B; // Convert to celsius T -= 273.15; - // Rescale to 0-255 for -100°C to +60°C - T = (T + 100.0) / 160.0 * 255.0; + // Rescale to 0-255 for -100°C to +60°C + T = (T + 100.0) / 160.0 * 255.0; return T; } // Temperature calibration wrapper void temperature(options_t *opts, image_t *img, int offset, int width){ - tempparam_t temp; + tempparam_t temp; - printf("Temperature... "); - fflush(stdout); + printf("Temperature... "); + fflush(stdout); - tempcomp(tele, img->chB, opts->satnum - 15, &temp); + tempcomp(tele, img->chB, opts->satnum - 15, &temp); - for (int y = 0; y < img->nrow; y++) { - for (int x = 0; x < width; x++) { + for (int y = 0; y < img->nrow; y++) { + for (int x = 0; x < width; x++) { img->prow[y][x + offset] = tempcal(img->prow[y][x + offset], opts->satnum - 15, &temp); } - } - printf("Done\n"); + } + printf("Done\n"); } \ No newline at end of file diff --git a/main.c b/main.c index 97af45c..d7537d9 100644 --- a/main.c +++ b/main.c @@ -66,7 +66,7 @@ static int processAudio(char *filename, options_t *opts); static void usage(void); int main(int argc, char **argv) { - fprintf(stderr, VERSION"\n"); + fprintf(stderr, VERSION"\n"); // Check if there are actually any input files if(argc == optind || argc == 1){ @@ -78,7 +78,7 @@ int main(int argc, char **argv) { // Parse arguments int opt; - while ((opt = getopt(argc, argv, "m:d:i:s:e:r")) != EOF) { + while ((opt = getopt(argc, argv, "m:d:i:s:e:r")) != EOF) { switch (opt) { case 'd': opts.path = optarg; @@ -105,14 +105,14 @@ int main(int argc, char **argv) { default: usage(); } - } + } // Process the files for (; optind < argc; optind++) { processAudio(argv[optind], &opts); } - exit(0); + exit(0); } static int processAudio(char *filename, options_t *opts){ @@ -124,7 +124,7 @@ static int processAudio(char *filename, options_t *opts){ char *id[7]; char *name[7]; } ch = { - { "?", "1", "2", "3A", "4", "5", "3B" }, + { "?", "1", "2", "3A", "4", "5", "3B" }, { "unknown", "visble", "near-infrared", "mid-infrared", "thermal-infrared", "thermal-infrared", "mid-infrared" } }; @@ -163,7 +163,7 @@ static int processAudio(char *filename, options_t *opts){ for (img.nrow = 0; img.nrow < MAX_HEIGHT; 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; @@ -259,40 +259,40 @@ static int processAudio(char *filename, options_t *opts){ // Distribution image if (CONTAINS(opts->type, 'd')) distrib(opts, &img, "d"); - + return 1; } static int initsnd(char *filename) { - SF_INFO infwav; - int res; + SF_INFO infwav; + int res; // Open audio file - infwav.format = 0; - audioFile = sf_open(filename, SFM_READ, &infwav); - if (audioFile == NULL) { + infwav.format = 0; + audioFile = sf_open(filename, SFM_READ, &infwav); + if (audioFile == NULL) { fprintf(stderr, "Could not open %s for reading\n", filename); return 0; - } + } - res = init_dsp(infwav.samplerate); + res = init_dsp(infwav.samplerate); printf("Input file: %s\n", filename); - if(res < 0) { + if(res < 0) { fprintf(stderr, "Input sample rate too low: %d\n", infwav.samplerate); return 0; - }else if(res > 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); + } + printf("Input sample rate: %d\n", infwav.samplerate); // TODO: accept stereo audio - if (infwav.channels != 1) { + if (infwav.channels != 1) { fprintf(stderr, "Too many channels in input file: %d\n", infwav.channels); return 0; - } + } - return 1; + return 1; } // Read samples from the wave file @@ -301,28 +301,28 @@ int getsample(float *sample, int nb) { } static void usage(void) { - fprintf(stderr, + fprintf(stderr, "Aptdec [options] audio files ...\n" "Options:\n" " -e [t|h|d|p|f|l] Effects\n" - " t: Crop telemetry\n" - " h: Histogram equalise\n" - " d: Denoise\n" - " p: Precipitation\n" - " f: Flip image\n" - " l: Linear equalise\n" + " t: Crop telemetry\n" + " h: Histogram equalise\n" + " d: Denoise\n" + " p: Precipitation\n" + " f: Flip image\n" + " l: Linear equalise\n" " -i [r|a|b|c|t|m] Output image\n" - " r: Raw\n" - " a: Channel A\n" - " b: Channel B\n" - " c: False color\n" - " t: Temperature\n" - " m: MCIR\n" - " -d Image destination directory.\n" - " -s [15-19] Satellite number\n" - " -m Map file\n" - " -r Realtime decode\n" + " r: Raw\n" + " a: Channel A\n" + " b: Channel B\n" + " c: False color\n" + " t: Temperature\n" + " m: MCIR\n" + " -d Image destination directory.\n" + " -s [15-19] Satellite number\n" + " -m Map file\n" + " -r Realtime decode\n" "\nRefer to the README for more infomation\n"); - exit(EINVAL); + exit(EINVAL); } diff --git a/median.c b/median.c index 43717e2..7e774db 100644 --- a/median.c +++ b/median.c @@ -7,24 +7,24 @@ #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]; - } + 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[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) */ @@ -51,6 +51,6 @@ float quick_select(float arr[], int n) { low = ll; if (hh >= median) high = hh - 1; - } + } } #undef ELEM_SWAP \ No newline at end of file diff --git a/pngio.c b/pngio.c index 7c00a83..bf1cd9e 100644 --- a/pngio.c +++ b/pngio.c @@ -191,8 +191,8 @@ png_text meta[] = { 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); - + sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid); + meta[1].text = desc; meta[1].text_length = sizeof(desc); @@ -208,24 +208,24 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c } // Create writer - png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if (!png_ptr) { + 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, "Could not create a PNG writer\n"); return 0; - } - png_infop info_ptr = png_create_info_struct(png_ptr); - if (!info_ptr) { + } + png_infop 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 writer\n"); return 0; - } + } if(palette == NULL && !CONTAINS(opts->effects, 'p') && !fc && opts->map[0] == '\0' && strcmp(chid, "MCIR") != 0){ greyscale = 1; // Greyscale image - png_set_IHDR(png_ptr, info_ptr, width, img->nrow, + png_set_IHDR(png_ptr, info_ptr, width, img->nrow, 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); }else{ @@ -235,17 +235,17 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); } - png_set_text(png_ptr, info_ptr, meta, 3); - png_set_pHYs(png_ptr, info_ptr, 3636, 3636, PNG_RESOLUTION_METER); + png_set_text(png_ptr, info_ptr, meta, 3); + png_set_pHYs(png_ptr, info_ptr, 3636, 3636, PNG_RESOLUTION_METER); // Init I/O - pngfile = fopen(outName, "wb"); - if (!pngfile) { + pngfile = fopen(outName, "wb"); + if (!pngfile) { fprintf(stderr, "Could not open %s for writing\n", outName); return 1; - } - png_init_io(png_ptr, pngfile); - png_write_info(png_ptr, info_ptr); + } + 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]; @@ -287,10 +287,10 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c printf("Writing %s", outName); // Build image - for (int y = 0; y < img->nrow; y++) { + for (int y = 0; y < img->nrow; y++) { png_color pix[width]; // Color png_byte mpix[width]; // Mono - + int skip = 0; for (int x = 0; x < width; x++) { if(skiptele){ @@ -320,21 +320,21 @@ int ImageOut(options_t *opts, image_t *img, int offset, int width, char *desc, c }; } } - + if(greyscale){ png_write_row(png_ptr, (png_bytep) mpix); }else{ png_write_row(png_ptr, (png_bytep) pix); } - } + } // Tidy up - png_write_end(png_ptr, info_ptr); - fclose(pngfile); - printf("\nDone\n"); - png_destroy_write_struct(&png_ptr, &info_ptr); + png_write_end(png_ptr, info_ptr); + fclose(pngfile); + printf("\nDone\n"); + png_destroy_write_struct(&png_ptr, &info_ptr); - return 1; + return 1; } // TODO: clean up everthing below this comment @@ -344,44 +344,44 @@ FILE *rt_pngfile; int initWriter(options_t *opts, image_t *img, int width, int height, char *desc, char *chid){ char outName[384]; - sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid); + sprintf(outName, "%s/%s-%s.png", opts->path, img->name, chid); meta[1].text = desc; meta[1].text_length = sizeof(desc); // Create writer - rt_png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if (!rt_png_ptr) { + rt_png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (!rt_png_ptr) { png_destroy_write_struct(&rt_png_ptr, (png_infopp) NULL); fprintf(stderr, "Could not create a PNG writer\n"); return 0; - } - rt_info_ptr = png_create_info_struct(rt_png_ptr); - if (!rt_info_ptr) { + } + rt_info_ptr = png_create_info_struct(rt_png_ptr); + if (!rt_info_ptr) { png_destroy_write_struct(&rt_png_ptr, (png_infopp) NULL); fprintf(stderr, "Could not create a PNG writer\n"); return 0; - } + } // Greyscale image - png_set_IHDR(rt_png_ptr, rt_info_ptr, width, height, + png_set_IHDR(rt_png_ptr, rt_info_ptr, width, height, 8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); - + png_set_text(rt_png_ptr, rt_info_ptr, meta, 3); // Channel = 25cm wide - png_set_pHYs(rt_png_ptr, rt_info_ptr, 3636, 3636, PNG_RESOLUTION_METER); + png_set_pHYs(rt_png_ptr, rt_info_ptr, 3636, 3636, PNG_RESOLUTION_METER); // Init I/O - rt_pngfile = fopen(outName, "wb"); - if (!rt_pngfile) { + rt_pngfile = fopen(outName, "wb"); + if (!rt_pngfile) { fprintf(stderr, "Could not open %s for writing\n", outName); return 0; - } - png_init_io(rt_png_ptr, rt_pngfile); - png_write_info(rt_png_ptr, rt_info_ptr); - + } + png_init_io(rt_png_ptr, rt_pngfile); + png_write_info(rt_png_ptr, rt_info_ptr); + return 1; } @@ -389,12 +389,12 @@ void pushRow(float *row, int width){ png_byte pix[width]; for(int i = 0; i < width; i++) pix[i] = row[i]; - + png_write_row(rt_png_ptr, (png_bytep) pix); } void closeWriter(){ png_write_end(rt_png_ptr, rt_info_ptr); - fclose(rt_pngfile); - png_destroy_write_struct(&rt_png_ptr, &rt_info_ptr); + fclose(rt_pngfile); + png_destroy_write_struct(&rt_png_ptr, &rt_info_ptr); } \ No newline at end of file diff --git a/reg.c b/reg.c index 316856d..238f74f 100644 --- a/reg.c +++ b/reg.c @@ -18,18 +18,18 @@ static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]); void polyreg(const int M, const int N, const double X[], const double Y[], double C[]) { - int R, K, J; /* Loop counters */ - double A[DMAX][DMAX]; /* A */ - double B[DMAX]; - double P[2 * DMAX + 1]; - double x, y; - double p; - - /* Zero the array */ - for (R = 0; R < M + 1; R++) + int R, K, J; /* Loop counters */ + double A[DMAX][DMAX]; /* A */ + double B[DMAX]; + double P[2 * DMAX + 1]; + double x, y; + double p; + + /* Zero the array */ + for (R = 0; R < M + 1; R++) B[R] = 0; - /* Compute the column vector */ + /* Compute the column vector */ for (K = 0; K < N; K++) { y = Y[K]; x = X[K]; @@ -38,50 +38,50 @@ void polyreg(const int M, const int N, const double X[], const double Y[], doubl B[R] += y * p; p = p * x; } - } + } - /* Zero the array */ - for (J = 1; J <= 2 * M; J++) + /* Zero the array */ + for (J = 1; J <= 2 * M; J++) P[J] = 0; - - P[0] = N; - /* Compute the sum of powers of x_(K-1) */ - for (K = 0; K < N; K++) { + P[0] = N; + + /* Compute the sum of powers of x_(K-1) */ + for (K = 0; K < N; K++) { x = X[K]; p = X[K]; for (J = 1; J <= 2 * M; J++) { P[J] += p; p = p * x; } - } + } - /* Determine the matrix entries */ - for (R = 0; R < M + 1; R++) { + /* Determine the matrix entries */ + for (R = 0; R < M + 1; R++) { for (K = 0; K < M + 1; K++) A[R][K] = P[R + K]; - } + } - /* Solve the linear system of M + 1 equations: A*C = B - for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */ - FactPiv(M + 1, A, B, C); + /* Solve the linear system of M + 1 equations: A*C = B + for the coefficient vector C = (c_1,c_2,..,c_M,c_(M+1)) */ + FactPiv(M + 1, A, B, C); } /* end main */ /*--------------------------------------------------------*/ static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) { - int K, P, C, J; /* Loop counters */ - int Row[NMAX]; /* Field with row-number */ - double X[DMAX], Y[DMAX]; - double SUM, DET = 1.0; - int T; - - /* Initialize the pointer vector */ - for (J = 0; J < NMAX; J++) + int K, P, C, J; /* Loop counters */ + int Row[NMAX]; /* Field with row-number */ + double X[DMAX], Y[DMAX]; + double SUM, DET = 1.0; + int T; + + /* Initialize the pointer vector */ + for (J = 0; J < NMAX; J++) Row[J] = J; - /* Start LU factorization */ - for (P = 0; P < N - 1; P++) { + /* Start LU factorization */ + for (P = 0; P < N - 1; P++) { /* Find pivot element */ for (K = P + 1; K < N; K++) { if (fabs(A[Row[K]][P]) > fabs(A[Row[P]][P])) { @@ -109,34 +109,34 @@ 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 */ - DET = DET * A[Row[N - 1]][N - 1]; + DET = DET * A[Row[N - 1]][N - 1]; - /* Start the forward substitution */ - for (K = 0; K < N; K++) + /* Start the forward substitution */ + for (K = 0; K < N; K++) Y[K] = B[K]; - + Y[0] = B[Row[0]]; - - for (K = 1; K < N; K++) { + + for (K = 1; K < N; K++) { SUM = 0; for (C = 0; C <= K - 1; C++) SUM += A[Row[K]][C] * Y[C]; Y[K] = B[Row[K]] - SUM; - } + } - /* Start the back substitution */ - X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1]; - for (K = N - 2; K >= 0; K--) { + /* Start the back substitution */ + X[N - 1] = Y[N - 1] / A[Row[N - 1]][N - 1]; + for (K = N - 2; K >= 0; K--) { SUM = 0; for (C = K + 1; C < N; C++) { SUM += A[Row[K]][C] * X[C]; } X[K] = (Y[K] - SUM) / A[Row[K]][K]; - } /* End of back substitution */ + } /* End of back substitution */ - /* Output */ - for (K = 0; K < N; K++) + /* Output */ + for (K = 0; K < N; K++) Cf[K] = X[K]; } diff --git a/satcal.h b/satcal.h index 6906f61..9160f5b 100644 --- a/satcal.h +++ b/satcal.h @@ -18,107 +18,107 @@ */ const struct { - float d[4][3]; - struct { - float vc, A, B; - } rad[3]; - struct { - float Ns; - float b[3]; - } cor[3]; + float d[4][3]; + struct { + float vc, A, B; + } rad[3]; + struct { + float Ns; + float b[3]; + } 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 - } + { // 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 + } }, { // 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 - } + { // 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 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 - } + { // 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 + } }, { // 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 - } + { // 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 + } }, { // 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 - } + { // 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;