@@ -1,18 +1,23 @@ | |||
CC=gcc | |||
BIN=/usr/bin | |||
INCLUDES=-I. | |||
CFLAGS= -O3 -Wall $(INCLUDES) | |||
OBJS= main.o image.o dsp.o filter.o reg.o fcolor.o | |||
atpdec: $(OBJS) | |||
aptdec: $(OBJS) | |||
$(CC) -o $@ $(OBJS) -lm -lsndfile -lpng | |||
main.o: main.c version.h temppalette.h gvipalette.h offsets.h | |||
dsp.o: dsp.c filtercoeff.h filter.h | |||
main.o: main.c version.h temppalette.h gvipalette.h offsets.h | |||
dsp.o: dsp.c filtercoeff.h filter.h | |||
filter.o: filter.c filter.h | |||
image.o: image.c satcal.h offsets.h | |||
fcolor.o : fcolor.c offsets.h | |||
image.o: image.c satcal.h offsets.h | |||
fcolor.o: fcolor.c offsets.h | |||
clean: | |||
rm -f *.o atpdec | |||
rm -f *.o aptdec | |||
install: | |||
install -m 755 aptdec $(BIN) | |||
uninstall: | |||
rm $(BIN)/aptdec |
@@ -0,0 +1,90 @@ | |||
# Aptdec | |||
Thierry Leconte F4DWV (c) 2004-2009 | |||
## Description | |||
Aptec is an FOSS program that decodes images transmitted by POES NOAA weather satellites. | |||
These satellites transmit continuously (among other things), medium resolution (1px/4km) images of the earth on 137 MHz. | |||
These transmissions could be easily received with an simple antenna and cheap SDR. | |||
Then the transimssion can easily be decoded in narrow band FM mode. | |||
Aptdec can convert these recordings into .png images. | |||
For each recording up to 6 images can be generated: | |||
1. Raw image: contains the 2 transmitted channel images + telemetry | |||
and synchronisation pulses. | |||
2. Calibrated channel A image | |||
3. Calibrated channel B image | |||
4. Temperature compensated I.R image | |||
5. False color image | |||
Input recordings must be mono with a sample rate of 11025 Hz. | |||
Aptdec uses `libsndfile` to read the input recording, so any format supported by `libsndfile` may be used (only tested with .wav files). | |||
## Compilation | |||
Aptdec is written is portable since it is written in standard C. | |||
It has only tested on Fedora and Debian, but will work on any Unix platform. | |||
Just edit the Makefile and run `make` (no configure script as of now). | |||
Aptdec uses `libsndfile`, `libpng` and `libm`. | |||
The `snd.h` and `png.h` headers must be present on your system. | |||
If they are not on standard path, edit the include path in the Makefile. | |||
## Usage | |||
To run without installing | |||
`./Aptdec [options] recordings ...` | |||
To install | |||
`sudo make install` | |||
To uninstall | |||
`sudo make uninstall` | |||
To run once installed | |||
`Aptdec [options] recordings ...` | |||
## Options | |||
``` | |||
-i [r|a|b|c|t] | |||
Toggle raw (r), channel A (a), channel B (b), false color (c), | |||
or temperature (t) output. | |||
Default: ac | |||
-d <dir> | |||
Optional images destination directory. | |||
Default: Recording directory. | |||
-s [15|16|17|18|19] | |||
Satellite number | |||
For temperature and false color generation | |||
Default: 19 | |||
-c <file> | |||
Use configuration file for false color generation. | |||
Default: Internal parameters. | |||
``` | |||
## Output | |||
Generated images are outputted in PNG, 8 bit greyscale for raw and channel A|B images, 24 bit RGB for false color. | |||
Image names are `recordingname-x.png`, where `x` is: | |||
- r for raw images | |||
- satellite instrument number (1, 2, 3A, 3B, 4, 5) for channel A|B images | |||
- c for false colors. | |||
## Example | |||
`Aptdec -d images -i ac *.wav` | |||
This will process all .wav files in the current directory, generate channel A and false color images and put them in the `images` directory. | |||
## License | |||
See LICENSE. |
@@ -1,5 +1,5 @@ | |||
/* | |||
* Atpdec | |||
* Aptec | |||
* Copyright (c) 2004 by Thierry Leconte (F4DWV) | |||
* | |||
* $Id$ | |||
@@ -23,19 +23,22 @@ | |||
#include <string.h> | |||
#include <math.h> | |||
#ifndef M_PI | |||
#define M_PI 3.14159265358979323846 /* for OS that don't know it */ | |||
#endif /* */ | |||
#define M_PI 3.14159265358979323846 /* for OS's that don't include it */ | |||
#endif | |||
#include "filter.h" | |||
#include "filtercoeff.h" | |||
extern int getsample(float *inbuff, int nb); | |||
#define BLKAMP 1024 | |||
#define BLKIN 1024 | |||
#define Fc 2400.0 | |||
#define DFc 50.0 | |||
#define PixelLine 2080 | |||
#define Fp (2*PixelLine) | |||
#define Fp (2 * PixelLine) | |||
#define RSMULT 15 | |||
#define Fi (Fp*RSMULT) | |||
#define Fi (Fp * RSMULT) | |||
static double Fe; | |||
@@ -43,72 +46,67 @@ static double offset = 0.0; | |||
static double FreqLine = 1.0; | |||
static double FreqOsc; | |||
static double K1,K2; | |||
static double K1, K2; | |||
int init_dsp(double F) | |||
{ | |||
if(F > Fi) return (1); | |||
if(F < Fp) return (-1); | |||
Fe=F; | |||
int init_dsp(double F) { | |||
if(F > Fi) return (1); | |||
if(F < Fp) return (-1); | |||
Fe = F; | |||
K1=DFc/Fe; | |||
K2=K1*K1/2.0; | |||
FreqOsc=Fc/Fe; | |||
K1 = DFc / Fe; | |||
K2 = K1 * K1 / 2.0; | |||
FreqOsc = Fc / Fe; | |||
return(0); | |||
return(0); | |||
} | |||
/* fast phase estimator */ | |||
static inline double Phase(double I,double Q) | |||
{ | |||
double angle,r; | |||
int s; | |||
if(I==0.0 && Q==0.0) | |||
return(0.0); | |||
if (Q < 0) { | |||
s=-1; | |||
Q=-Q; | |||
} else { | |||
s=1; | |||
} | |||
static inline double Phase(double I,double Q) { | |||
double angle,r; | |||
int s; | |||
if(I == 0.0 && Q == 0.0) | |||
return(0.0); | |||
if (Q < 0) { | |||
s = -1; | |||
Q = -Q; | |||
} 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); | |||
} | |||
static double pll(double I, double Q) | |||
{ | |||
/* pll coeff */ | |||
static double pll(double I, double Q) { | |||
/* pll coeff */ | |||
static double PhaseOsc = 0.0; | |||
double Io, Qo; | |||
double Ip, Qp; | |||
double DPhi; | |||
/* quadrature oscillator */ | |||
/* quadrature oscillator */ | |||
Io = cos(PhaseOsc); | |||
Qo = sin(PhaseOsc); | |||
/* phase detector */ | |||
Ip = I*Io+Q*Qo; | |||
Qp = Q*Io-I*Qo; | |||
DPhi = Phase(Ip,Qp); | |||
/* loop filter */ | |||
/* phase detector */ | |||
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; | |||
@@ -124,10 +122,7 @@ static double pll(double I, double Q) | |||
return (Ip); | |||
} | |||
static int getamp(double *ambuff, int nb) | |||
{ | |||
#define BLKIN 1024 | |||
static int getamp(double *ambuff, int nb) { | |||
static float inbuff[BLKIN]; | |||
static int idxin=0; | |||
static int nin=0; | |||
@@ -135,70 +130,66 @@ static int getamp(double *ambuff, int nb) | |||
int n; | |||
for (n = 0; n < nb; n++) { | |||
double I,Q; | |||
if (nin < IQFilterLen*2+2) { | |||
int res; | |||
memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); | |||
idxin = 0; | |||
res = getsample(&(inbuff[nin]), BLKIN - nin); | |||
nin += res; | |||
if (nin < IQFilterLen*2+2) | |||
return (n); | |||
} | |||
iqfir(&inbuff[idxin],iqfilter,IQFilterLen,&I,&Q); | |||
ambuff[n] = pll(I,Q); | |||
idxin += 1; | |||
nin -= 1; | |||
double I, Q; | |||
if (nin < IQFilterLen * 2 + 2) { | |||
int res; | |||
memmove(inbuff, &(inbuff[idxin]), nin * sizeof(float)); | |||
idxin = 0; | |||
res = getsample(&(inbuff[nin]), BLKIN - nin); | |||
nin += res; | |||
if (nin < IQFilterLen * 2 + 2) | |||
return (n); | |||
} | |||
iqfir(&inbuff[idxin], iqfilter, IQFilterLen, &I, &Q); | |||
ambuff[n] = pll(I, Q); | |||
idxin += 1; | |||
nin -= 1; | |||
} | |||
return (n); | |||
} | |||
int getpixelv(float *pvbuff, int nb) | |||
{ | |||
#define BLKAMP 1024 | |||
int getpixelv(float *pvbuff, int nb) { | |||
static double ambuff[BLKAMP]; | |||
static int nam = 0; | |||
static int idxam = 0; | |||
int n,m; | |||
int n, m; | |||
double mult; | |||
mult = (double) Fi/Fe*FreqLine; | |||
mult = (double) Fi / Fe * FreqLine; | |||
m=RSFilterLen/mult+1; | |||
m = RSFilterLen / mult + 1; | |||
for (n = 0; n < nb; n++) { | |||
int shift; | |||
if (nam < m) { | |||
int res; | |||
memmove(ambuff, &(ambuff[idxam]), nam * sizeof(double)); | |||
idxam = 0; | |||
res = getamp(&(ambuff[nam]), BLKAMP - nam); | |||
nam += res; | |||
if (nam < m) | |||
return (n); | |||
} | |||
int shift; | |||
pvbuff[n] = rsfir(&(ambuff[idxam]), rsfilter, RSFilterLen, offset, mult) * mult * 256.0; | |||
if (nam < m) { | |||
int res; | |||
memmove(ambuff, &(ambuff[idxam]), nam * sizeof(double)); | |||
idxam = 0; | |||
res = getamp(&(ambuff[nam]), BLKAMP - nam); | |||
nam += res; | |||
if (nam < m) | |||
return (n); | |||
} | |||
//printf("%g\n",pvbuff[n]); | |||
pvbuff[n] = rsfir(&(ambuff[idxam]), rsfilter, RSFilterLen, offset, mult) * mult * 256.0; | |||
shift = ((int) floor((RSMULT - offset) / mult))+1; | |||
offset = shift*mult+offset-RSMULT ; | |||
//printf("%g\n",pvbuff[n]); | |||
idxam += shift; | |||
nam -= shift; | |||
shift = ((int) floor((RSMULT - offset) / mult)) + 1; | |||
offset = shift * mult + offset - RSMULT ; | |||
idxam += shift; | |||
nam -= shift; | |||
} | |||
return (nb); | |||
} | |||
int getpixelrow(float *pixelv) | |||
{ | |||
int getpixelrow(float *pixelv) { | |||
static float pixels[PixelLine + SyncFilterLen]; | |||
static int npv = 0; | |||
static int synced = 0; | |||
@@ -208,11 +199,12 @@ int getpixelrow(float *pixelv) | |||
int res; | |||
if (npv > 0) | |||
memmove(pixelv, pixels, npv * sizeof(float)); | |||
memmove(pixelv, pixels, npv * sizeof(float)); | |||
if (npv < SyncFilterLen + 2) { | |||
res = getpixelv(&(pixelv[npv]), SyncFilterLen + 2 - npv); | |||
npv += res; | |||
if (npv < SyncFilterLen + 2) | |||
res = getpixelv(&(pixelv[npv]), SyncFilterLen + 2 - npv); | |||
npv += res; | |||
if (npv < SyncFilterLen + 2) | |||
return (0); | |||
} | |||
@@ -222,16 +214,15 @@ int getpixelrow(float *pixelv) | |||
lcorr = fir(&(pixelv[2]), Sync, SyncFilterLen); | |||
FreqLine = 1.0+((ecorr-lcorr) / corr / PixelLine / 4.0); | |||
if (corr < 0.75 * max) { | |||
synced = 0; | |||
FreqLine = 1.0; | |||
synced = 0; | |||
FreqLine = 1.0; | |||
} | |||
max = corr; | |||
if (synced < 8) { | |||
int shift, mshift; | |||
if (npv < PixelLine + SyncFilterLen) { | |||
res = | |||
getpixelv(&(pixelv[npv]), PixelLine + SyncFilterLen - npv); | |||
res = getpixelv(&(pixelv[npv]), PixelLine + SyncFilterLen - npv); | |||
npv += res; | |||
if (npv < PixelLine + SyncFilterLen) | |||
return (0); | |||
@@ -244,31 +235,30 @@ int getpixelrow(float *pixelv) | |||
corr = fir(&(pixelv[shift + 1]), Sync, SyncFilterLen); | |||
if (corr > max) { | |||
mshift = shift; | |||
max = corr; | |||
mshift = shift; | |||
max = corr; | |||
} | |||
} | |||
if (mshift != 0) { | |||
memmove(pixelv, &(pixelv[mshift]), | |||
(npv - mshift) * sizeof(float)); | |||
memmove(pixelv, &(pixelv[mshift]), (npv - mshift) * sizeof(float)); | |||
npv -= mshift; | |||
synced = 0; | |||
FreqLine = 1.0; | |||
} else | |||
synced += 1; | |||
} | |||
if (npv < PixelLine) { | |||
res = getpixelv(&(pixelv[npv]), PixelLine - npv); | |||
npv += res; | |||
if (npv < PixelLine) | |||
return (0); | |||
res = getpixelv(&(pixelv[npv]), PixelLine - npv); | |||
npv += res; | |||
if (npv < PixelLine) | |||
return (0); | |||
} | |||
if (npv == PixelLine) { | |||
npv = 0; | |||
npv = 0; | |||
} else { | |||
memmove(pixels, &(pixelv[PixelLine]), | |||
(npv - PixelLine) * sizeof(float)); | |||
npv -= PixelLine; | |||
memmove(pixels, &(pixelv[PixelLine]), (npv - PixelLine) * sizeof(float)); | |||
npv -= PixelLine; | |||
} | |||
return (1); | |||
@@ -7,4 +7,3 @@ | |||
240.0 0.6 0.4 | |||
60.0 0.6 0.2 | |||
100.0 0.1 0.5 | |||
@@ -6,15 +6,14 @@ typedef struct { | |||
float h, s, v; | |||
} hsvpix_t; | |||
static void HSVtoRGB(float *r, float *g, float *b, hsvpix_t pix) | |||
{ | |||
static void HSVtoRGB(float *r, float *g, float *b, hsvpix_t pix) { | |||
int i; | |||
double f, p, q, t, h; | |||
if (pix.s == 0) { | |||
// achromatic (grey) | |||
*r = *g = *b = pix.v; | |||
return; | |||
// achromatic (grey) | |||
*r = *g = *b = pix.v; | |||
return; | |||
} | |||
h = pix.h / 60; // sector 0 to 5 | |||
@@ -25,44 +24,44 @@ static void HSVtoRGB(float *r, float *g, float *b, hsvpix_t pix) | |||
t = pix.v * (1 - pix.s * (1 - f)); | |||
switch (i) { | |||
case 0: | |||
*r = pix.v; | |||
*g = t; | |||
*b = p; | |||
break; | |||
case 1: | |||
*r = q; | |||
*g = pix.v; | |||
*b = p; | |||
break; | |||
case 2: | |||
*r = p; | |||
*g = pix.v; | |||
*b = t; | |||
break; | |||
case 3: | |||
*r = p; | |||
*g = q; | |||
*b = pix.v; | |||
break; | |||
case 4: | |||
*r = t; | |||
*g = p; | |||
*b = pix.v; | |||
break; | |||
default: // case 5: | |||
*r = pix.v; | |||
*g = p; | |||
*b = q; | |||
break; | |||
case 0: | |||
*r = pix.v; | |||
*g = t; | |||
*b = p; | |||
break; | |||
case 1: | |||
*r = q; | |||
*g = pix.v; | |||
*b = p; | |||
break; | |||
case 2: | |||
*r = p; | |||
*g = pix.v; | |||
*b = t; | |||
break; | |||
case 3: | |||
*r = p; | |||
*g = q; | |||
*b = pix.v; | |||
break; | |||
case 4: | |||
*r = t; | |||
*g = p; | |||
*b = pix.v; | |||
break; | |||
default: // case 5: | |||
*r = pix.v; | |||
*g = p; | |||
*b = q; | |||
break; | |||
} | |||
} | |||
static struct { | |||
float Seathresold; | |||
float Landthresold; | |||
float Tthresold; | |||
float Seathreshold; | |||
float Landthreshold; | |||
float Threshold; | |||
hsvpix_t CloudTop; | |||
hsvpix_t CloudBot; | |||
hsvpix_t SeaTop; | |||
@@ -79,58 +78,47 @@ static struct { | |||
100.0, 0.0, 0.5} | |||
}; | |||
void readfconf(char *file) | |||
{ | |||
void readfconf(char *file) { | |||
FILE *fin; | |||
fin = fopen(file, "r"); | |||
if (fin == NULL) | |||
return; | |||
fscanf(fin, "%g\n", &fcinfo.Seathresold); | |||
fscanf(fin, "%g\n", &fcinfo.Landthresold); | |||
fscanf(fin, "%g\n", &fcinfo.Tthresold); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.CloudTop.h, &fcinfo.CloudTop.s, | |||
&fcinfo.CloudTop.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.CloudBot.h, &fcinfo.CloudBot.s, | |||
&fcinfo.CloudBot.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.SeaTop.h, &fcinfo.SeaTop.s, | |||
&fcinfo.SeaTop.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.SeaBot.h, &fcinfo.SeaBot.s, | |||
&fcinfo.SeaBot.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.GroundTop.h, &fcinfo.GroundTop.s, | |||
&fcinfo.GroundTop.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.GroundBot.h, &fcinfo.GroundBot.s, | |||
&fcinfo.GroundBot.v); | |||
return; | |||
fscanf(fin, "%g\n", &fcinfo.Seathreshold); | |||
fscanf(fin, "%g\n", &fcinfo.Landthreshold); | |||
fscanf(fin, "%g\n", &fcinfo.Threshold); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.CloudTop.h, &fcinfo.CloudTop.s, &fcinfo.CloudTop.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.CloudBot.h, &fcinfo.CloudBot.s, &fcinfo.CloudBot.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.SeaTop.h, &fcinfo.SeaTop.s, &fcinfo.SeaTop.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.SeaBot.h, &fcinfo.SeaBot.s, &fcinfo.SeaBot.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.GroundTop.h, &fcinfo.GroundTop.s, &fcinfo.GroundTop.v); | |||
fscanf(fin, "%g %g %g\n", &fcinfo.GroundBot.h, &fcinfo.GroundBot.s, &fcinfo.GroundBot.v); | |||
fclose(fin); | |||
}; | |||
void falsecolor(double v, double t, float *r, float *g, float *b) | |||
{ | |||
void falsecolor(double v, double t, float *r, float *g, float *b) { | |||
hsvpix_t top, bot, c; | |||
double scv, sct; | |||
if (t > fcinfo.Tthresold) { | |||
if (v < fcinfo.Seathresold) { | |||
/* sea */ | |||
top = fcinfo.SeaTop, bot = fcinfo.SeaBot; | |||
scv = v / fcinfo.Seathresold; | |||
sct = (256.0-t) / (256.0-fcinfo.Tthresold); | |||
} else { | |||
/* ground */ | |||
top = fcinfo.GroundTop, bot = fcinfo.GroundBot; | |||
scv = | |||
(v - fcinfo.Seathresold) / (fcinfo.Landthresold - | |||
fcinfo.Seathresold); | |||
sct = (256.0-t) / (256.0-fcinfo.Tthresold); | |||
} | |||
if (t > fcinfo.Threshold) { | |||
if (v < fcinfo.Seathreshold) { | |||
/* sea */ | |||
top = fcinfo.SeaTop, bot = fcinfo.SeaBot; | |||
scv = v / fcinfo.Seathreshold; | |||
sct = (256.0 - t) / (256.0 - fcinfo.Threshold); | |||
} else { | |||
/* ground */ | |||
top = fcinfo.GroundTop, bot = fcinfo.GroundBot; | |||
scv = (v - fcinfo.Seathreshold) / (fcinfo.Landthreshold - fcinfo.Seathreshold); | |||
sct = (256.0 - t) / (256.0 - fcinfo.Threshold); | |||
} | |||
} else { | |||
/* clouds */ | |||
top = fcinfo.CloudTop, bot = fcinfo.CloudBot; | |||
scv = v / 256.0; | |||
sct = (256.0-t) / 256.0; | |||
/* clouds */ | |||
top = fcinfo.CloudTop, bot = fcinfo.CloudBot; | |||
scv = v / 256.0; | |||
sct = (256.0 - t) / 256.0; | |||
} | |||
c.s = top.s + sct * (bot.s - top.s); | |||
@@ -140,11 +128,10 @@ void falsecolor(double v, double t, float *r, float *g, float *b) | |||
HSVtoRGB(r, g, b, c); | |||
}; | |||
void Ngvi(float **prow, int nrow) | |||
{ | |||
void Ngvi(float **prow, int nrow) { | |||
int n; | |||
printf("GVI ... "); | |||
printf("GVI... "); | |||
fflush(stdout); | |||
for (n = 0; n < nrow; n++) { | |||
@@ -154,11 +141,11 @@ void Ngvi(float **prow, int nrow) | |||
pixelv = prow[n]; | |||
for (i = 0; i < CH_WIDTH; i++) { | |||
float pv; | |||
double gvi; | |||
double gvi; | |||
gvi = (pixelv[i + CHA_OFFSET]-pixelv[i + CHB_OFFSET])/(pixelv[i + CHA_OFFSET]+pixelv[i + CHB_OFFSET]); | |||
gvi = (pixelv[i + CHA_OFFSET] - pixelv[i + CHB_OFFSET]) / (pixelv[i + CHA_OFFSET] + pixelv[i + CHB_OFFSET]); | |||
pv = (gvi+0.1)*340.0; | |||
pv = (gvi + 0.1) * 340.0; | |||
if (pv > 255.0) | |||
pv = 255.0; | |||
if (pv < 0.0) | |||
@@ -1,5 +1,5 @@ | |||
/* | |||
* Atpdec | |||
* Aptec | |||
* Copyright (c) 2004 by Thierry Leconte (F4DWV) | |||
* | |||
* $Id$ | |||
@@ -22,47 +22,43 @@ | |||
#include "filter.h" | |||
#include <math.h> | |||
float fir(float *buff, const float *coeff, const int len) | |||
{ | |||
float fir(float *buff, const float *coeff, const int len) { | |||
int i; | |||
double r; | |||
r = 0.0; | |||
for (i = 0; i < len; i++) { | |||
r += buff[i] * coeff[i]; | |||
r += buff[i] * coeff[i]; | |||
} | |||
return r; | |||
} | |||
void iqfir(float *buff, const float *coeff, const int len,double *I,double *Q) | |||
{ | |||
void iqfir(float *buff, const float *coeff, const int len, double *I, double *Q) { | |||
int k; | |||
double i,q; | |||
double i, q; | |||
i=q=0.0; | |||
i = q = 0.0; | |||
for (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; | |||
q += buff[2*k] * coeff[k]; | |||
i += buff[2*k]; | |||
} | |||
i= buff[len-1] - i / len; | |||
*I=i, *Q=q; | |||
} | |||
float rsfir(double *buff, const float *coeff, const int len, const double offset, | |||
const double delta) | |||
{ | |||
float rsfir(double *buff, const float *coeff, const int len, const double offset, const double delta) { | |||
int i; | |||
double n; | |||
double out; | |||
out = 0.0; | |||
for (i = 0, n = offset; i < (len-1)/delta-1; n += delta, i++) { | |||
int k; | |||
double alpha; | |||
int k; | |||
double alpha; | |||
k = (int)floor(n); | |||
alpha = n - k; | |||
out += buff[i]*(coeff[k]*(1.0-alpha)+coeff[k + 1]*alpha); | |||
k = (int)floor(n); | |||
alpha = n - k; | |||
out += buff[i] * (coeff[k] * (1.0 - alpha) + coeff[k + 1] * alpha); | |||
} | |||
return out; | |||
} | |||
@@ -1,5 +1,5 @@ | |||
/* | |||
* Atpdec | |||
* Aptec | |||
* Copyright (c) 2003 by Thierry Leconte (F4DWV) | |||
* | |||
* $Id$ | |||
@@ -20,8 +20,7 @@ | |||
* | |||
*/ | |||
float fir(float *buff,const float *coeff,const int len); | |||
void iqfir(float *buff,const float *coeff,const int len,double *I,double *Q); | |||
float rsfir(double *buff,const float *coeff,const int len ,const double offset ,const double delta); | |||
float fir(float *buff, const float *coeff, const int len); | |||
void iqfir(float *buff, const float *coeff, const int len, double *I, double *Q); | |||
float rsfir(double *buff, const float *coeff, const int len, const double offset, const double delta); | |||
@@ -1,5 +1,5 @@ | |||
/* | |||
* Atpdec | |||
* Aptec | |||
* Copyright (c) 2003 by Thierry Leconte (F4DWV) | |||
* | |||
* $Id$ | |||
@@ -21,24 +21,19 @@ | |||
*/ | |||
#define IQFilterLen 32 | |||
const float iqfilter[IQFilterLen] = | |||
{ 0.0205361, 0.0219524, 0.0235785, 0.0254648, 0.0276791, 0.0303152, | |||
const float iqfilter[IQFilterLen] = { 0.0205361, 0.0219524, 0.0235785, 0.0254648, 0.0276791, 0.0303152, | |||
0.0335063, 0.0374482, 0.0424413, 0.0489708, 0.0578745, 0.0707355, 0.0909457, 0.127324, 0.212207, 0.63662, | |||
-0.63662, -0.212207, -0.127324, -0.0909457, -0.0707355, -0.0578745, -0.0489708, -0.0424413, -0.0374482, | |||
-0.0335063, -0.0303152, -0.0276791, -0.0254648, -0.0235785, -0.0219524, -0.0205361 | |||
}; | |||
-0.0335063, -0.0303152, -0.0276791, -0.0254648, -0.0235785, -0.0219524, -0.0205361 }; | |||
#define SyncFilterLen 32 | |||
const float Sync[SyncFilterLen] = | |||
{ -14, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, | |||
18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -14 | |||
}; | |||
const float Sync[SyncFilterLen] = { -14, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, | |||
18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, 18, 18, -14, -14, -14 }; | |||
#define RSFilterLen 437 | |||
const float rsfilter[RSFilterLen] = | |||
{ -3.37279e-04, -8.80292e-06, -3.96418e-04, -1.78544e-04, -5.27511e-04, | |||
const float rsfilter[RSFilterLen] = { -3.37279e-04, -8.80292e-06, -3.96418e-04, -1.78544e-04, -5.27511e-04, | |||
-3.75376e-04, -6.95337e-04, -5.93148e-04, -8.79730e-04, -8.15327e-04, -1.05669e-03, -1.01377e-03, | |||
-1.19836e-03, -1.15443e-03, -1.26937e-03, -1.20955e-03, -1.23904e-03, -1.15302e-03, -1.08660e-03, | |||
-9.64235e-04, -8.02450e-04, -6.46202e-04, -3.95376e-04, -2.18096e-04, 1.11906e-04, 2.89567e-04, | |||
@@ -100,7 +95,6 @@ const float rsfilter[RSFilterLen] = | |||
-6.46202e-04, -8.02450e-04, -9.64235e-04, -1.08660e-03, -1.15302e-03, -1.23904e-03, -1.20955e-03, | |||
-1.26937e-03, -1.15443e-03, -1.19836e-03, -1.01377e-03, -1.05669e-03, -8.15327e-04, -8.79730e-04, | |||
-5.93148e-04, -6.95337e-04, -3.75376e-04, -5.27511e-04, -1.78544e-04, -3.96418e-04, -8.80292e-06, | |||
-3.37279e-04 | |||
}; | |||
-3.37279e-04 }; | |||
@@ -1,4 +1,4 @@ | |||
unsigned char GviPalette[256*3]= { | |||
unsigned char GviPalette[256*3] = { | |||
"\230t\17\233x\22\236{\27\241\200\33\244\203\37\247\210#\252\214'\255\220" | |||
",\260\2240\264\2305\267\2358\272\240=\274\245A\300\251E\303\255I\306\262" | |||
"M\311\266Q\314\272V\317\276Z\322\302^\325\306b\330\312g\334\317k\337\323" | |||
@@ -27,6 +27,5 @@ | |||
"\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\3g", | |||
}; | |||
"\0\33u\0\25r\0\21p\0\14l\0\7j\0\3g" | |||
}; |
@@ -1,5 +1,5 @@ | |||
/* | |||
* Atpdec | |||
* Aptec | |||
* Copyright (c) 2004 by Thierry Leconte (F4DWV) | |||
* | |||
* $Id$ | |||
@@ -29,28 +29,24 @@ | |||
#define REGORDER 3 | |||
typedef struct { | |||
double cf[REGORDER + 1]; | |||
double cf[REGORDER + 1]; | |||
} rgparam; | |||
static void rgcomp(double x[16], rgparam * rgpr) | |||
{ | |||
/*{ 0.106,0.215,0.324,0.433,0.542,0.652,0.78,0.87 ,0.0 }; */ | |||
const double y[9] = | |||
{ 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; | |||
extern void polyreg(const int m, const int n, const double x[], | |||
const double y[], double c[]); | |||
static void rgcomp(double x[16], rgparam * rgpr) { | |||
/*{ 0.106,0.215,0.324,0.433,0.542,0.652,0.78,0.87 ,0.0 }; */ | |||
const double y[9] = { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 }; | |||
extern void polyreg(const int m, const int n, const double x[], const double y[], double c[]); | |||
polyreg(REGORDER, 9, x, y, rgpr->cf); | |||
polyreg(REGORDER, 9, x, y, rgpr -> cf); | |||
} | |||
static double rgcal(float x, rgparam * rgpr) | |||
{ | |||
static double rgcal(float x, rgparam * rgpr) { | |||
double y, p; | |||
int i; | |||
for (i = 0, y = 0.0, p = 1.0; i < REGORDER + 1; i++) { | |||
y += rgpr->cf[i] * p; | |||
p = p * x; | |||
y += rgpr->cf[i] * p; | |||
p = p * x; | |||
} | |||
return (y); | |||
} | |||
@@ -60,9 +56,7 @@ static double tele[16]; | |||
static double Cs; | |||
static int nbtele; | |||
int Calibrate(float **prow, int nrow, int offset) | |||
{ | |||
int Calibrate(float **prow, int nrow, int offset) { | |||
double teleline[3000]; | |||
double wedge[16]; | |||
rgparam regr[30]; | |||
@@ -72,143 +66,142 @@ int Calibrate(float **prow, int nrow, int offset) | |||
int channel = -1; | |||
float max; | |||
printf("Calibration ... "); | |||
printf("Calibration... "); | |||
fflush(stdout); | |||
/* build telemetry values lines */ | |||
/* build telemetry values lines */ | |||
for (n = 0; n < nrow; n++) { | |||
int i; | |||
int i; | |||
teleline[n] = 0.0; | |||
for (i = 3; i < 43; i++) { | |||
teleline[n] += prow[n][i + offset + CH_WIDTH]; | |||
} | |||
teleline[n] /= 40.0; | |||
teleline[n] = 0.0; | |||
for (i = 3; i < 43; i++) { | |||
teleline[n] += prow[n][i + offset + CH_WIDTH]; | |||
} | |||
teleline[n] /= 40.0; | |||
} | |||
if (nrow < 192) { | |||
fprintf(stderr, " impossible, not enought row\n"); | |||
return (0); | |||
fprintf(stderr, " not possible, not enough rows!\n"); | |||
return (0); | |||
} | |||
/* find telemetry start in the 2nd third */ | |||
/* find telemetry start in the 2nd third */ | |||
max = 0.0; | |||
mtelestart = 0; | |||
for (n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) { | |||
float df; | |||
df = (teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + | |||
teleline[n - 1]) / (teleline[n] + teleline[n + 1] + | |||
teleline[n + 2] + teleline[n + 3]); | |||
if (df > max) { | |||
mtelestart = n; | |||
max = df; | |||
} | |||
float df; | |||
df = (teleline[n - 4] + teleline[n - 3] + teleline[n - 2] + | |||
teleline[n - 1]) / (teleline[n] + teleline[n + 1] + | |||
teleline[n + 2] + teleline[n + 3]); | |||
if (df > max) { | |||
mtelestart = n; | |||
max = df; | |||
} | |||
} | |||
mtelestart -= 64; | |||
telestart = mtelestart % 128; | |||
if (mtelestart < 0 || nrow < telestart + 128) { | |||
fprintf(stderr, " impossible, not enought row\n"); | |||
return (0); | |||
fprintf(stderr, " impossible, not enough row\n"); | |||
return (0); | |||
} | |||
/* compute wedges and regression */ | |||
/* compute wedges and regression */ | |||
for (n = telestart, k = 0; n < nrow - 128; n += 128, k++) { | |||
int j; | |||
for (j = 0; j < 16; j++) { | |||
int i; | |||
wedge[j] = 0.0; | |||
for (i = 1; i < 7; i++) | |||
wedge[j] += teleline[n + j * 8 + i]; | |||
wedge[j] /= 6; | |||
} | |||
rgcomp(wedge, &(regr[k])); | |||
int j; | |||
if (k == nrow / 256) { | |||
int i, l; | |||
for (j = 0; j < 16; j++) { | |||
int i; | |||
/* telemetry calibration */ | |||
for (j = 0; j < 16; j++) { | |||
tele[j] = rgcal(wedge[j], &(regr[k])); | |||
} | |||
/* channel ID */ | |||
for (j = 0, max = 10000.0, channel = -1; j < 6; j++) { | |||
float df; | |||
df = wedge[15] - wedge[j]; | |||
df = df * df; | |||
if (df < max) { | |||
channel = j; | |||
max = df; | |||
wedge[j] = 0.0; | |||
for (i = 1; i < 7; i++) | |||
wedge[j] += teleline[n + j * 8 + i]; | |||
wedge[j] /= 6; | |||
} | |||
} | |||
/* Cs computation */ | |||
for (Cs = 0.0, i = 0, j = n; j < n + 128; j++) { | |||
double csline; | |||
for (csline = 0.0, l = 3; l < 43; l++) | |||
csline += prow[n][l + offset -SPC_WIDTH]; | |||
csline /= 40.0; | |||
if (csline > 50.0) { | |||
Cs += csline; | |||
i++; | |||
rgcomp(wedge, &(regr[k])); | |||
if (k == nrow / 256) { | |||
int i, l; | |||
/* telemetry calibration */ | |||
for (j = 0; j < 16; j++) { | |||
tele[j] = rgcal(wedge[j], &(regr[k])); | |||
} | |||
/* channel ID */ | |||
for (j = 0, max = 10000.0, channel = -1; j < 6; j++) { | |||
float df; | |||
df = wedge[15] - wedge[j]; | |||
df = df * df; | |||
if (df < max) { | |||
channel = j; | |||
max = df; | |||
} | |||
} | |||
/* Cs computation */ | |||
for (Cs = 0.0, i = 0, j = n; j < n + 128; j++) { | |||
double csline; | |||
for (csline = 0.0, l = 3; l < 43; l++) | |||
csline += prow[n][l + offset -SPC_WIDTH]; | |||
csline /= 40.0; | |||
if (csline > 50.0) { | |||
Cs += csline; | |||
i++; | |||
} | |||
} | |||
Cs /= i; | |||
Cs = rgcal(Cs, &(regr[k])); | |||
} | |||
} | |||
Cs /= i; | |||
Cs = rgcal(Cs, &(regr[k])); | |||
} | |||
} | |||
nbtele = k; | |||
/* calibrate */ | |||
/* calibrate */ | |||
for (n = 0; n < nrow; n++) { | |||
float *pixelv; | |||
int i; | |||
pixelv = prow[n]; | |||
for (i = 0; i < CH_WIDTH; i++) { | |||
float pv; | |||
int k, kof; | |||
pv = pixelv[i + offset]; | |||
k = (n - telestart) / 128; | |||
if (k >= nbtele) | |||
k = nbtele - 1; | |||
kof = (n - telestart) % 128; | |||
if (kof < 64) { | |||
if (k < 1) { | |||
pv = rgcal(pv, &(regr[k])); | |||
} else { | |||
pv = rgcal(pv, &(regr[k])) * (64 + kof) / 128.0 + | |||
rgcal(pv, &(regr[k - 1])) * (64 - kof) / 128.0; | |||
} | |||
} else { | |||
if ((k + 1) >= nbtele) { | |||
pv = rgcal(pv, &(regr[k])); | |||
} else { | |||
pv = rgcal(pv, &(regr[k])) * (192 - kof) / 128.0 + | |||
rgcal(pv, &(regr[k + 1])) * (kof - 64) / 128.0; | |||
float *pixelv; | |||
int i; | |||
pixelv = prow[n]; | |||
for (i = 0; i < CH_WIDTH; i++) { | |||
float pv; | |||
int k, kof; | |||
pv = pixelv[i + offset]; | |||
k = (n - telestart) / 128; | |||
if (k >= nbtele) | |||
k = nbtele - 1; | |||
kof = (n - telestart) % 128; | |||
if (kof < 64) { | |||
if (k < 1) { | |||
pv = rgcal(pv, &(regr[k])); | |||
} else { | |||
pv = rgcal(pv, &(regr[k])) * (64 + kof) / 128.0 + rgcal(pv, &(regr[k - 1])) * (64 - kof) / 128.0; | |||
} | |||
} else { | |||
if ((k + 1) >= nbtele) { | |||
pv = rgcal(pv, &(regr[k])); | |||
} else { | |||
pv = rgcal(pv, &(regr[k])) * (192 - kof) / 128.0 + rgcal(pv, &(regr[k + 1])) * (kof - 64) / 128.0; | |||
} | |||
} | |||
if (pv > 255.0) | |||
pv = 255.0; | |||
if (pv < 0.0) | |||
pv = 0.0; | |||
pixelv[i + offset] = pv; | |||
} | |||
} | |||
if (pv > 255.0) | |||
pv = 255.0; | |||
if (pv < 0.0) | |||
pv = 0.0; | |||
pixelv[i + offset] = pv; | |||
} | |||
} | |||
printf("Done\n"); | |||
return (channel+1); | |||
return (channel + 1); | |||
} | |||
/* ------------------------------temperature calibration -----------------------*/ | |||
@@ -223,51 +216,47 @@ typedef struct { | |||
} tempparam; | |||
/* temperature compensation for IR channel */ | |||
static void tempcomp(double t[16], int ch, tempparam * tpr) | |||
{ | |||
static void tempcomp(double t[16], int ch, tempparam * tpr) { | |||
double Tbb, T[4]; | |||
double C; | |||
int n; | |||
tpr->ch = ch - 4; | |||
tpr -> ch = ch - 4; | |||
/* compute equivalent T black body */ | |||
for (n = 0; n < 4; n++) { | |||
float d0, d1, d2; | |||
C = t[9 + n] * 4.0; | |||
d0 = satcal[satnum].d[n][0]; | |||
d1 = satcal[satnum].d[n][1]; | |||
d2 = satcal[satnum].d[n][2]; | |||
T[n] = d0; | |||
T[n] += d1 * C; | |||
C = C * C; | |||
T[n] += d2 * C; | |||
float d0, d1, d2; | |||
C = t[9 + n] * 4.0; | |||
d0 = satcal[satnum].d[n][0]; | |||
d1 = satcal[satnum].d[n][1]; | |||
d2 = satcal[satnum].d[n][2]; | |||
T[n] = d0; | |||
T[n] += d1 * C; | |||
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 = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb; | |||
/* compute radiance Black body */ | |||
C = satcal[satnum].rad[tpr->ch].vc; | |||
tpr->Nbb = c1 * C * C * C / (exp(c2 * C / Tbb) - 1.0); | |||
/* store Count Blackbody and space */ | |||
tpr->Cs = Cs * 4.0; | |||
tpr->Cb = t[14] * 4.0; | |||
} | |||
static double tempcal(float Ce, tempparam * rgpr) | |||
{ | |||
static double tempcal(float Ce, tempparam * rgpr) { | |||
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); | |||
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; | |||
satcal[satnum].cor[rgpr->ch].b[1] * Nl + | |||
satcal[satnum].cor[rgpr->ch].b[2] * Nl * Nl; | |||
Ne = Nl + Nc; | |||
@@ -275,38 +264,37 @@ static double tempcal(float Ce, tempparam * rgpr) | |||
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; | |||
/* rescale to range 0-255 for -60 +40 °C */ | |||
/* rescale to range 0-255 for -60 +40 'C */ | |||
T = (T - 273.15 + 60.0) / 100.0 * 256.0; | |||
return (T); | |||
} | |||
void Temperature(float **prow, int nrow, int channel, int offset) | |||
{ | |||
void Temperature(float **prow, int nrow, int channel, int offset) { | |||
tempparam temp; | |||
int n; | |||
printf("Temperature ... "); | |||
printf("Temperature... "); | |||
fflush(stdout); | |||
tempcomp(tele, channel, &temp); | |||
for (n = 0; n < nrow; n++) { | |||
float *pixelv; | |||
int i; | |||
float *pixelv; | |||
int i; | |||
pixelv = prow[n]; | |||
for (i = 0; i < CH_WIDTH; i++) { | |||
float pv; | |||
pixelv = prow[n]; | |||
for (i = 0; i < CH_WIDTH; i++) { | |||
float pv; | |||
pv = tempcal(pixelv[i + offset], &temp); | |||
pv = tempcal(pixelv[i + offset], &temp); | |||
if (pv > 255.0) | |||
pv = 255.0; | |||
if (pv < 0.0) | |||
pv = 0.0; | |||
pixelv[i + offset] = pv; | |||
} | |||
if (pv > 255.0) | |||
pv = 255.0; | |||
if (pv < 0.0) | |||
pv = 0.0; | |||
pixelv[i + offset] = pv; | |||
} | |||
} | |||
printf("Done\n"); | |||
} |
@@ -1,5 +1,5 @@ | |||
/* | |||
* Atpdec | |||
* Aptdec | |||
* Copyright (c) 2004-2005 by Thierry Leconte (F4DWV) | |||
* | |||
* $Id$ | |||
@@ -45,87 +45,78 @@ extern int init_dsp(double F);; | |||
static SNDFILE *inwav; | |||
static int initsnd(char *filename) | |||
{ | |||
static int initsnd(char *filename) { | |||
SF_INFO infwav; | |||
int res; | |||
int res; | |||
/* open wav input file */ | |||
/* open wav input file */ | |||
infwav.format = 0; | |||
inwav = sf_open(filename, SFM_READ, &infwav); | |||
if (inwav == NULL) { | |||
fprintf(stderr, "could not open %s\n", filename); | |||
return (1); | |||
fprintf(stderr, "Could not open %s for reading\n", filename); | |||
return (1); | |||
} | |||
res=init_dsp(infwav.samplerate); | |||
if(res<0) { | |||
fprintf(stderr, "Sample rate too low : %d\n", infwav.samplerate); | |||
return (1); | |||
res = init_dsp(infwav.samplerate); | |||
if(res < 0) { | |||
fprintf(stderr, "Sample rate too low: %d\n", infwav.samplerate); | |||
return (1); | |||
} | |||
if(res>0) { | |||
fprintf(stderr, "Sample rate too hight : %d\n", infwav.samplerate); | |||
return (1); | |||
if(res > 0) { | |||
fprintf(stderr, "Sample rate too hight: %d\n", infwav.samplerate); | |||
return (1); | |||
} | |||
fprintf(stderr, "Sample rate : %d\n", infwav.samplerate); | |||
fprintf(stderr, "Sample rate: %d\n", infwav.samplerate); | |||
if (infwav.channels != 1) { | |||
fprintf(stderr, "Too many channels in input file : %d\n", infwav.channels); | |||
return (1); | |||
fprintf(stderr, "Too many channels in input file: %d\n", infwav.channels); | |||
return (1); | |||
} | |||
return (0); | |||
} | |||
int getsample(float *sample, int nb) | |||
{ | |||
int getsample(float *sample, int nb) { | |||
return (sf_read_float(inwav, sample, nb)); | |||
} | |||
static png_text text_ptr[] = { | |||
{PNG_TEXT_COMPRESSION_NONE, "Software", version, sizeof(version)} | |||
, | |||
{PNG_TEXT_COMPRESSION_NONE, "Channel", NULL, 0} | |||
, | |||
{PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA POES satellite Image", | |||
25} | |||
{PNG_TEXT_COMPRESSION_NONE, "Software", version, sizeof(version)}, | |||
{PNG_TEXT_COMPRESSION_NONE, "Channel", NULL, 0}, | |||
{PNG_TEXT_COMPRESSION_NONE, "Description", "NOAA POES satellite image", 25} | |||
}; | |||
static int | |||
ImageOut(char *filename, char *chid, float **prow, int nrow, | |||
int width, int offset, png_color *palette) | |||
{ | |||
static int ImageOut(char *filename, char *chid, float **prow, int nrow, int width, int offset, png_color *palette) { | |||
FILE *pngfile; | |||
png_infop info_ptr; | |||
png_structp png_ptr; | |||
int n; | |||
/* init png lib */ | |||
png_ptr = | |||
png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); | |||
/* init png lib */ | |||
png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); | |||
if (!png_ptr) { | |||
fprintf(stderr, "could not open create png_ptr\n"); | |||
return (1); | |||
fprintf(stderr, "Could not create a PNG write struct\n"); | |||
return (1); | |||
} | |||
info_ptr = png_create_info_struct(png_ptr); | |||
if (!info_ptr) { | |||
png_destroy_write_struct(&png_ptr, (png_infopp) NULL); | |||
fprintf(stderr, "could not open create info_ptr\n"); | |||
return (1); | |||
png_destroy_write_struct(&png_ptr, (png_infopp) NULL); | |||
fprintf(stderr, "Could not create a PNG info struct\n"); | |||
return (1); | |||
} | |||
if(palette==NULL) { | |||
/* grey image */ | |||
if(palette == NULL) { | |||
/* greyscale */ | |||
png_set_IHDR(png_ptr, info_ptr, width, nrow, | |||
8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, | |||
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); | |||
8, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, | |||
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); | |||
} else { | |||
/* palette color mage */ | |||
/* palette color mage */ | |||
png_set_IHDR(png_ptr, info_ptr, width, nrow, | |||
8, PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE, | |||
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); | |||
png_set_PLTE(png_ptr, info_ptr, palette, 256); | |||
8, PNG_COLOR_TYPE_PALETTE, PNG_INTERLACE_NONE, | |||
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); | |||
png_set_PLTE(png_ptr, info_ptr, palette, 256); | |||
} | |||
text_ptr[1].text = chid; | |||
@@ -133,26 +124,26 @@ ImageOut(char *filename, char *chid, float **prow, int nrow, | |||
png_set_text(png_ptr, info_ptr, text_ptr, 3); | |||
png_set_pHYs(png_ptr, info_ptr, 4000, 4000, PNG_RESOLUTION_METER); | |||
printf("Writing %s ... ", filename); | |||
printf("Writing %s... ", filename); | |||
fflush(stdout); | |||
pngfile = fopen(filename, "wb"); | |||
if (pngfile == NULL) { | |||
fprintf(stderr, "could not open %s\n", filename); | |||
return (1); | |||
fprintf(stderr, "Could not open %s for writing\n", filename); | |||
return (1); | |||
} | |||
png_init_io(png_ptr, pngfile); | |||
png_write_info(png_ptr, info_ptr); | |||
for (n = 0; n < nrow; n++) { | |||
float *pixelv; | |||
png_byte pixel[2*IMG_WIDTH]; | |||
int i; | |||
float *pixelv; | |||
png_byte pixel[2*IMG_WIDTH]; | |||
int i; | |||
pixelv = prow[n]; | |||
for (i = 0; i < width; i++) { | |||
pixel[i] = pixelv[i + offset]; | |||
} | |||
png_write_row(png_ptr, pixel); | |||
pixelv = prow[n]; | |||
for (i = 0; i < width; i++) { | |||
pixel[i] = pixelv[i + offset]; | |||
} | |||
png_write_row(png_ptr, pixel); | |||
} | |||
png_write_end(png_ptr, info_ptr); | |||
fclose(pngfile); | |||
@@ -161,72 +152,69 @@ ImageOut(char *filename, char *chid, float **prow, int nrow, | |||
return (0); | |||
} | |||
static int ImageRGBOut(char *filename, float **prow, int nrow) | |||
{ | |||
static int ImageRGBOut(char *filename, float **prow, int nrow) { | |||
FILE *pngfile; | |||
png_infop info_ptr; | |||
png_structp png_ptr; | |||
int n; | |||
extern void falsecolor(double v, double t, float *r, float *g, | |||
float *b); | |||
extern void falsecolor(double v, double t, float *r, float *g, float *b); | |||
/* init png lib */ | |||
png_ptr = | |||
png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); | |||
/* init png lib */ | |||
png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); | |||
if (!png_ptr) { | |||
fprintf(stderr, "could not open create png_ptr\n"); | |||
return (1); | |||
fprintf(stderr, "Could not create a PNG write struct\n"); | |||
return (1); | |||
} | |||
info_ptr = png_create_info_struct(png_ptr); | |||
if (!info_ptr) { | |||
png_destroy_write_struct(&png_ptr, (png_infopp) NULL); | |||
fprintf(stderr, "could not open create info_ptr\n"); | |||
return (1); | |||
png_destroy_write_struct(&png_ptr, (png_infopp) NULL); | |||
fprintf(stderr, "Could not create a PNG info struct\n"); | |||
return (1); | |||
} | |||
png_set_IHDR(png_ptr, info_ptr, CH_WIDTH , nrow , | |||
8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, | |||
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); | |||
8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, | |||
PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT); | |||
png_set_pHYs(png_ptr, info_ptr, 4000, 4000, PNG_RESOLUTION_METER); | |||
text_ptr[1].text = "False Colors"; | |||
text_ptr[1].text = "False Color"; | |||
text_ptr[1].text_length = strlen(text_ptr[1].text); | |||
png_set_text(png_ptr, info_ptr, text_ptr, 3); | |||
printf("Computing False colors & writing : %s ...", filename); | |||
printf("Computing false color & writing: %s... ", filename); | |||
fflush(stdout); | |||
pngfile = fopen(filename, "wb"); | |||
if (pngfile == NULL) { | |||
fprintf(stderr, "could not open %s\n", filename); | |||
return (1); | |||
fprintf(stderr, "Could not open %s for writing\n", filename); | |||
return (1); | |||
} | |||
png_init_io(png_ptr, pngfile); | |||
png_write_info(png_ptr, info_ptr); | |||
for (n = 0; n < nrow ; n++) { | |||
png_color pix[CH_WIDTH]; | |||
float *pixelc; | |||
int i; | |||
png_color pix[CH_WIDTH]; | |||
float *pixelc; | |||
int i; | |||
pixelc = prow[n]; | |||
pixelc = prow[n]; | |||
for (i = 0; i < CH_WIDTH - 1; i++) { | |||
float v, t; | |||
float r, g, b; | |||
for (i = 0; i < CH_WIDTH - 1; i++) { | |||
float v, t; | |||
float r, g, b; | |||
v = pixelc[i+CHA_OFFSET]; | |||
t = pixelc[i+CHB_OFFSET]; | |||
v = pixelc[i+CHA_OFFSET]; | |||
t = pixelc[i+CHB_OFFSET]; | |||
falsecolor(v, t, &r, &g, &b); | |||
falsecolor(v, t, &r, &g, &b); | |||
pix[i].red = 255.0 * r; | |||
pix[i].green = 255.0 * g; | |||
pix[i].blue = 255.0 * b; | |||
} | |||
png_write_row(png_ptr, (png_bytep) pix); | |||
pix[i].red = 255.0 * r; | |||
pix[i].green = 255.0 * g; | |||
pix[i].blue = 255.0 * b; | |||
} | |||
png_write_row(png_ptr, (png_bytep) pix); | |||
} | |||
png_write_end(png_ptr, info_ptr); | |||
fclose(pngfile); | |||
@@ -236,44 +224,41 @@ static int ImageRGBOut(char *filename, float **prow, int nrow) | |||
} | |||
static void Distrib(char *filename,float **prow,int nrow) | |||
{ | |||
unsigned int distrib[256][256]; | |||
int n; | |||
int x,y; | |||
int max=0; | |||
FILE *df; | |||
static void Distrib(char *filename,float **prow,int nrow) { | |||
unsigned int distrib[256][256]; | |||
int n; | |||
int x, y; | |||
int max = 0; | |||
FILE *df; | |||
for(y=0;y<256;y++) | |||
for(x=0;x<256;x++) | |||
distrib[y][x]=0; | |||
for(y = 0; y < 256; y++) | |||
for(x = 0; x < 256; x++) | |||
distrib[y][x] = 0; | |||
for(n=0;n<nrow;n++) { | |||
float *pixelv; | |||
int i; | |||
for(n = 0; n < nrow; n++) { | |||
float *pixelv; | |||
int i; | |||
pixelv=prow[n]; | |||
for(i=0;i<CH_WIDTH;i++) { | |||
y=(int)(pixelv[i+CHA_OFFSET]); | |||
x=(int)(pixelv[i+CHB_OFFSET]); | |||
distrib[y][x]+=1; | |||
if(distrib[y][x]> max) max=distrib[y][x]; | |||
} | |||
} | |||
df=fopen(filename,"w"); | |||
pixelv = prow[n]; | |||
for(i = 0; i < CH_WIDTH; i++) { | |||
y = (int)(pixelv[i + CHA_OFFSET]); | |||
x = (int)(pixelv[i + CHB_OFFSET]); | |||
distrib[y][x] += 1; | |||
if(distrib[y][x] > max) max=distrib[y][x]; | |||
} | |||
} | |||
df = fopen(filename,"w"); | |||
printf("Writing %s\n",filename); | |||
printf("Writing %s\n",filename); | |||
fprintf(df,"P2\n#max %d\n",max); | |||
fprintf(df,"256 256\n255\n"); | |||
for(y=0;y<256;y++) | |||
for(x=0;x<256;x++) | |||
fprintf(df,"%d\n",(int)((255.0*(double)(distrib[y][x]))/(double)max)); | |||
fclose(df); | |||
fprintf(df,"P2\n#max %d\n",max); | |||
fprintf(df,"256 256\n255\n"); | |||
for(y=0;y<256;y++) | |||
for(x=0;x<256;x++) | |||
fprintf(df, "%d\n", (int)((255.0 * (double)(distrib[y][x])) / (double)max)); | |||
fclose(df); | |||
} | |||
extern int Calibrate(float **prow, int nrow, int offset); | |||
extern void Temperature(float **prow, int nrow, int ch, int offset); | |||
extern int Ngvi(float **prow, int nrow); | |||
@@ -282,16 +267,13 @@ extern int optind, opterr; | |||
extern char *optarg; | |||
int satnum = 4; | |||
static void usage(void) | |||
{ | |||
fprintf(stderr, "atpdec [options] soundfiles ...\n"); | |||
fprintf(stderr, | |||
"options:\n-d <dir>\tDestination directory\n-i [r|a|b|c|t]\tOutput image type\n\t\t\tr: Raw\n\t\t\ta: A chan.\n\t\t\tb: B chan.\n\t\t\tc: False color\n\t\t\tt: Temperature\n-c <file>\tFalse color config file\n-s [15|16|17|18|19]\tSatellite number (for temperature and false color generation)\n"); | |||
exit(1); | |||
static void usage(void) { | |||
fprintf(stderr, "Aptdec [options] recordings ...\n"); | |||
fprintf(stderr, "Options:\n -i [r|a|b|c|t] Output image type\n r: Raw\n a: A channel\n b: B channel\n c: False color\n t: Temperature\n -d <dir> Image destination directory.\n -s [15|16|17|18|19] Satellite number\n -c <file> False color config file\n"); | |||
exit(1); | |||
} | |||
int main(int argc, char **argv) | |||
{ | |||
int main(int argc, char **argv) { | |||
char pngfilename[1024]; | |||
char name[1024]; | |||
char pngdirname[1024] = ""; | |||
@@ -299,126 +281,124 @@ int main(int argc, char **argv) | |||
float *prow[3000]; | |||
char *chid[] = { "?", "1", "2", "3A", "4", "5", "3B" }; | |||
int nrow; | |||
int chA,chB; | |||
int chA, chB; | |||
int c; | |||
printf("%s\n", version); | |||
opterr = 0; | |||
while ((c = getopt(argc, argv, "c:d:i:s:")) != EOF) { | |||
switch (c) { | |||
case 'd': | |||
strcpy(pngdirname, optarg); | |||
break; | |||
case 'c': | |||
readfconf(optarg); | |||
break; | |||
case 'i': | |||
strcpy(imgopt, optarg); | |||
break; | |||
case 's': | |||
satnum = atoi(optarg)-15; | |||
if (satnum < 0 || satnum > 4) { | |||
fprintf(stderr, "invalid satellite number : must be in [15-19]\n"); | |||
exit(1); | |||
} | |||
break; | |||
default: | |||
usage(); | |||
if(argc == 1){ | |||
usage(); | |||
} | |||
while ((c = getopt(argc, argv, "c:d:i:s:")) != EOF) { | |||
switch (c) { | |||
case 'd': | |||
strcpy(pngdirname, optarg); | |||
break; | |||
case 'c': | |||
readfconf(optarg); | |||
break; | |||
case 'i': | |||
strcpy(imgopt, optarg); | |||
break; | |||
case 's': | |||
satnum = atoi(optarg)-15; | |||
if (satnum < 0 || satnum > 4) { | |||
fprintf(stderr, "Invalid satellite number, must be in the range [15-19]\n"); | |||
exit(1); | |||
} | |||
break; | |||
default: | |||
usage(); | |||
} | |||
} | |||
for (nrow = 0; nrow < 3000; nrow++) | |||
prow[nrow] = NULL; | |||
prow[nrow] = NULL; | |||
for (; optind < argc; optind++) { | |||
chA=chB=0; | |||
strcpy(pngfilename, argv[optind]); | |||
strcpy(name, basename(pngfilename)); | |||
strtok(name, "."); | |||
if (pngdirname[0] == '\0') { | |||
strcpy(pngfilename, argv[optind]); | |||
strcpy(pngdirname, dirname(pngfilename)); | |||
} | |||
/* open snd input */ | |||
if (initsnd(argv[optind])) | |||
exit(1); | |||
/* main loop */ | |||
printf("Decoding: %s \n", argv[optind]); | |||
for (nrow = 0; nrow < 3000; nrow++) { | |||
if (prow[nrow] == NULL) | |||
prow[nrow] = (float *) malloc(sizeof(float) * 2150); | |||
if (getpixelrow(prow[nrow]) == 0) | |||
break; | |||
printf("%d\r", nrow); | |||
fflush(stdout); | |||
} | |||
printf("\nDone\n"); | |||
sf_close(inwav); | |||
for (; optind < argc; optind++) { | |||
chA = chB = 0; | |||
strcpy(pngfilename, argv[optind]); | |||
strcpy(name, basename(pngfilename)); | |||
strtok(name, "."); | |||
if (pngdirname[0] == '\0') { | |||
strcpy(pngfilename, argv[optind]); | |||
strcpy(pngdirname, dirname(pngfilename)); | |||
} | |||
/* raw image */ | |||
if (strchr(imgopt, (int) 'r') != NULL) { | |||
sprintf(pngfilename, "%s/%s-r.png", pngdirname, name); | |||
ImageOut(pngfilename, "raw", prow, nrow, IMG_WIDTH, 0,NULL); | |||
} | |||
/* open snd input */ | |||
if (initsnd(argv[optind])) | |||
exit(1); | |||
/* main loop */ | |||
printf("Decoding: %s \n", argv[optind]); | |||
for (nrow = 0; nrow < 3000; nrow++) { | |||
if (prow[nrow] == NULL) | |||
prow[nrow] = (float *) malloc(sizeof(float) * 2150); | |||
if (getpixelrow(prow[nrow]) == 0) | |||
break; | |||
printf("%d\r", nrow); | |||
fflush(stdout); | |||
} | |||
printf("\nDone\n"); | |||
sf_close(inwav); | |||
/* Channel A */ | |||
if (((strchr(imgopt, (int) 'a') != NULL) | |||
|| (strchr(imgopt, (int) 'c') != NULL) | |||
|| (strchr(imgopt, (int) 'd') != NULL))) { | |||
chA = Calibrate(prow, nrow, CHA_OFFSET); | |||
if (chA >= 0) { | |||
if (strchr(imgopt, (int) 'a') != NULL) { | |||
sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chA]); | |||
ImageOut(pngfilename, chid[chA], prow, nrow, CH_WIDTH , CHA_OFFSET,NULL); | |||
/* raw image */ | |||
if (strchr(imgopt, (int) 'r') != NULL) { | |||
sprintf(pngfilename, "%s/%s-r.png", pngdirname, name); | |||
ImageOut(pngfilename, "raw", prow, nrow, IMG_WIDTH, 0, NULL); | |||
} | |||
} | |||
} | |||
/* Channel B */ | |||
if ((strchr(imgopt, (int) 'b') != NULL) | |||
|| (strchr(imgopt, (int) 'c') != NULL) | |||
|| (strchr(imgopt, (int) 't') != NULL) | |||
|| (strchr(imgopt, (int) 'd') != NULL)) { | |||
chB = Calibrate(prow, nrow, CHB_OFFSET); | |||
if (chB >= 0) { | |||
if (strchr(imgopt, (int) 'b') != NULL) { | |||
sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chB]); | |||
ImageOut(pngfilename, chid[chB], prow, nrow, CH_WIDTH , CHB_OFFSET ,NULL); | |||
/* Channel A */ | |||
if (((strchr(imgopt, (int) 'a') != NULL) | |||
|| (strchr(imgopt, (int) 'c') != NULL) | |||
|| (strchr(imgopt, (int) 'd') != NULL))) { | |||
chA = Calibrate(prow, nrow, CHA_OFFSET); | |||
if (chA >= 0 && strchr(imgopt, (int) 'a') != NULL) { | |||
sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chA]); | |||
ImageOut(pngfilename, chid[chA], prow, nrow, CH_WIDTH , CHA_OFFSET, NULL); | |||
} | |||
} | |||
} | |||
if (chB > 3) { | |||
Temperature(prow, nrow, chB, CHB_OFFSET); | |||
if (strchr(imgopt, (int) 't') != NULL) { | |||
sprintf(pngfilename, "%s/%s-t.png", pngdirname, name); | |||
ImageOut(pngfilename, "Temperature", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)TempPalette); | |||
/* Channel B */ | |||
if ((strchr(imgopt, (int) 'b') != NULL) | |||
|| (strchr(imgopt, (int) 'c') != NULL) | |||
|| (strchr(imgopt, (int) 't') != NULL) | |||
|| (strchr(imgopt, (int) 'd') != NULL)) { | |||
chB = Calibrate(prow, nrow, CHB_OFFSET); | |||
if (chB >= 0 && strchr(imgopt, (int) 'b') != NULL) { | |||
sprintf(pngfilename, "%s/%s-%s.png", pngdirname, name, chid[chB]); | |||
ImageOut(pngfilename, chid[chB], prow, nrow, CH_WIDTH , CHB_OFFSET ,NULL); | |||
} | |||
if (chB > 3) { | |||
Temperature(prow, nrow, chB, CHB_OFFSET); | |||
if (strchr(imgopt, (int) 't') != NULL) { | |||
sprintf(pngfilename, "%s/%s-t.png", pngdirname, name); | |||
ImageOut(pngfilename, "Temperature", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)TempPalette); | |||
} | |||
} | |||
} | |||
} | |||
} | |||
/* distribution */ | |||
if (chA && chB && strchr(imgopt, (int) 'd') != NULL) { | |||
sprintf(pngfilename, "%s/%s-d.pnm", pngdirname, name); | |||
Distrib(pngfilename, prow, nrow); | |||
} | |||
/* distribution */ | |||
if (chA && chB && strchr(imgopt, (int) 'd') != NULL) { | |||
sprintf(pngfilename, "%s/%s-d.pnm", pngdirname, name); | |||
Distrib(pngfilename, prow, nrow); | |||
} | |||
/* color image */ | |||
if (chA==2 && chB==4 && strchr(imgopt, (int) 'c') != NULL) { | |||
sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); | |||
ImageRGBOut(pngfilename, prow, nrow); | |||
} | |||
/* color image */ | |||
if (chA == 2 && chB == 4 && strchr(imgopt, (int) 'c') != NULL) { | |||
sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); | |||
ImageRGBOut(pngfilename, prow, nrow); | |||
} | |||
/* GVI image */ | |||
if (chA==1 && chB==2 && strchr(imgopt, (int) 'c') != NULL) { | |||
Ngvi(prow, nrow); | |||
sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); | |||
ImageOut(pngfilename, "GVI", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)GviPalette); | |||
} | |||
/* GVI image */ | |||
if (chA == 1 && chB == 2 && strchr(imgopt, (int) 'c') != NULL) { | |||
Ngvi(prow, nrow); | |||
sprintf(pngfilename, "%s/%s-c.png", pngdirname, name); | |||
ImageOut(pngfilename, "GVI", prow, nrow, CH_WIDTH, CHB_OFFSET, (png_color*)GviPalette); | |||
} | |||
} | |||
exit(0); | |||
} |
@@ -1,7 +1,7 @@ | |||
#define SYNC_WIDTH 39 | |||
#define SPC_WIDTH 47 | |||
#define SPC_WIDTH 47 | |||
#define TELE_WIDTH 45 | |||
#define CH_WIDTH 909 | |||
#define CH_WIDTH 909 | |||
#define CH_OFFSET (SYNC_WIDTH+SPC_WIDTH+CH_WIDTH+TELE_WIDTH) | |||
#define IMG_WIDTH 2080 | |||
#define CHA_OFFSET (SYNC_WIDTH+SPC_WIDTH) | |||
@@ -1,6 +1,6 @@ | |||
/* --------------------------------------------------------------------------- | |||
Polynomial regression, freely adapted from : | |||
Polynomial regression, freely adapted from: | |||
NUMERICAL METHODS: C Programs, (c) John H. Mathews 1995 | |||
Algorithm translated to C by: Dr. Norman Fahrer | |||
@@ -10,36 +10,33 @@ | |||
E-mail address: in%"mathews@fullerton.edu" | |||
*/ | |||
#include<math.h> | |||
#include <math.h> | |||
#define DMAX 5 /* Maximum degree of polynomial */ | |||
#define NMAX 10 /* Maximum number of points */ | |||
#define DMAX 5 /* Maximum degree of polynomial */ | |||
#define NMAX 10 /* Maximum number of points */ | |||
static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]); | |||
void | |||
polyreg(const int M, const int N, const double X[], const double Y[], | |||
double C[]) | |||
{ | |||
int R, K, J; /* Loop counters */ | |||
double A[DMAX][DMAX]; /* A */ | |||
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++) | |||
B[R] = 0; | |||
B[R] = 0; | |||
/* Compute the column vector */ | |||
for (K = 0; K < N; K++) { | |||
y = Y[K]; | |||
x = X[K]; | |||
p = 1.0; | |||
for (R = 0; R < M + 1; R++) { | |||
B[R] += y * p; | |||
p = p * x; | |||
} | |||
for (K = 0; K < N; K++) { | |||
y = Y[K]; | |||
x = X[K]; | |||
p = 1.0; | |||
for (R = 0; R < M + 1; R++) { | |||
B[R] += y * p; | |||
p = p * x; | |||
} | |||
} | |||
/* Zero the array */ | |||
@@ -49,31 +46,30 @@ polyreg(const int M, const int N, const double X[], const double Y[], | |||
/* 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; | |||
} | |||
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++) { | |||
for (K = 0; K < M + 1; K++) | |||
A[R][K] = P[R + K]; | |||
for (K = 0; K < M + 1; K++) | |||
A[R][K] = P[R + K]; | |||
} | |||
/* Solve the linear system of M + 1 equations : A*C = B | |||
/* 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 */ | |||
} /* 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 */ | |||
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; | |||
@@ -81,64 +77,63 @@ static void FactPiv(int N, double A[DMAX][DMAX], double B[], double Cf[]) | |||
/* Initialize the pointer vector */ | |||
for (J = 0; J < N; J++) | |||
Row[J] = J; | |||
Row[J] = J; | |||
/* 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])) { | |||
/* Switch the index for the p-1 th pivot row if necessary */ | |||
T = Row[P]; | |||
Row[P] = Row[K]; | |||
Row[K] = T; | |||
DET = -DET; | |||
} | |||
} /* End of simulated row interchange */ | |||
if (A[Row[P]][P] == 0) { | |||
/* The matrix is SINGULAR ! */ | |||
return; | |||
} | |||
/* Multiply the diagonal elements */ | |||
DET = DET * A[Row[P]][P]; | |||
/* Form multiplier */ | |||
for (K = P + 1; K < N; K++) { | |||
A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P]; | |||
/* Eliminate X_(p-1) */ | |||
for (C = P + 1; C < N + 1; C++) { | |||
A[Row[K]][C] -= A[Row[K]][P] * A[Row[P]][C]; | |||
} | |||
} | |||
/* Find pivot element */ | |||
for (K = P + 1; K < N; K++) { | |||
if (fabs(A[Row[K]][P]) > fabs(A[Row[P]][P])) { | |||
/* Switch the index for the p-1 th pivot row if necessary */ | |||
T = Row[P]; | |||
Row[P] = Row[K]; | |||
Row[K] = T; | |||
DET = -DET; | |||
} | |||
} /* End of simulated row interchange */ | |||
if (A[Row[P]][P] == 0) { | |||
/* The matrix is SINGULAR! */ | |||
return; | |||
} | |||
/* Multiply the diagonal elements */ | |||
DET = DET * A[Row[P]][P]; | |||
/* Form multiplier */ | |||
for (K = P + 1; K < N; K++) { | |||
A[Row[K]][P] = A[Row[K]][P] / A[Row[P]][P]; | |||
/* Eliminate X_(p-1) */ | |||
for (C = P + 1; C < N + 1; C++) { | |||
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]; | |||
/* Start the forward substitution */ | |||
for (K = 0; K < N; K++) | |||
Y[K] = B[K]; | |||
Y[0] = B[Row[0]]; | |||
Y[K] = B[K]; | |||
Y[0] = B[Row[0]]; | |||
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; | |||
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--) { | |||
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 */ | |||
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 */ | |||
/* Output */ | |||
for (K = 0; K < N; K++) | |||
Cf[K] = X[K]; | |||
Cf[K] = X[K]; | |||
} |
@@ -104,7 +104,5 @@ const struct { | |||
} | |||
}; | |||
const float c1=1.1910427e-5; | |||
const float c2=1.4387752; | |||
const float c2=1.4387752; |
@@ -1,4 +1,4 @@ | |||
unsigned char TempPalette[256*3]= { | |||
unsigned char TempPalette[256*3] = { | |||
"\376\376\376\376\376\376\375\375\376\374\375\376\374\375\375\374\373\375" | |||
"\373\373\375\372\373\375\372\373\374\372\372\374\371\372\374\371\371\375" | |||
"\370\371\374\367\370\375\367\370\374\367\367\374\366\367\373\366\366\373" | |||
@@ -1,2 +1,2 @@ | |||
char version[]="Atpdec CVS version (c) 2004-2005 Thierry Leconte F4DWV"; | |||
char version[] = "Aptec CVS version (c) 2004-2005 Thierry Leconte F4DWV"; | |||
@@ -1,60 +1,57 @@ | |||
#include <stdio.h> | |||
#include <string.h> | |||
char * basename (const char *filename) | |||
{ | |||
const char *p = filename + strlen (filename); | |||
while (p != filename) | |||
{ | |||
p--; | |||
if (*p == '/' || *p == '\\') | |||
return ((char *) p + 1); | |||
} | |||
return ((char *) filename); | |||
char * basename (const char *filename) { | |||
const char *p = filename + strlen (filename); | |||
while (p != filename) { | |||
p--; | |||
if (*p == '/' || *p == '\\') | |||
return ((char *) p + 1); | |||
} | |||
return ((char *) filename); | |||
} | |||
char *dirname(const char *path) | |||
{ | |||
static char bname[1024]; | |||
register const char *endp; | |||
char *dirname(const char *path) { | |||
static char bname[1024]; | |||
register const char *endp; | |||
/* Empty or NULL string gets treated as "." */ | |||
if (path == NULL || *path == '\0') { | |||
(void)strncpy(bname, ".", sizeof bname); | |||
return(bname); | |||
} | |||
/* Empty or NULL string gets treated as "." */ | |||
if (path == NULL || *path == '\0') { | |||
(void)strncpy(bname, ".", sizeof bname); | |||
return(bname); | |||
} | |||
/* Strip trailing slashes */ | |||
endp = path + strlen(path) - 1; | |||
while (endp > path && *endp == '\\' ) | |||
endp--; | |||
/* Strip trailing slashes */ | |||
endp = path + strlen(path) - 1; | |||
while (endp > path && *endp == '\\') | |||
endp--; | |||
/* Find the start of the dir */ | |||
while (endp > path && *endp != '\\' ) | |||
endp--; | |||
while (endp > path && *endp != '\\') | |||
endp--; | |||
/* Either the dir is "/" or there are no slashes */ | |||
if (endp == path) { | |||
(void)strncpy(bname, *endp == '\\' ? "\\" : ".", sizeof bname); | |||
return(bname); | |||
(void)strncpy(bname, *endp == '\\' ? "\\": ".", sizeof bname); | |||
return(bname); | |||
} else { | |||
do { | |||
endp--; | |||
} while (endp > path && *endp == '\\'); | |||
do { | |||
endp--; | |||
} while (endp > path && *endp == '\\'); | |||
} | |||
if (endp - path + 2 > sizeof(bname)) { | |||
return(NULL); | |||
return(NULL); | |||
} | |||
strncpy(bname, path, endp - path + 2); | |||
return(bname); | |||
} | |||
int opterr = 1, /* if error message should be printed */ | |||
optind = 1, /* index into parent argv vector */ | |||
optopt, /* character checked for validity */ | |||
optreset; /* reset getopt */ | |||
char *optarg; /* argument associated with option */ | |||
int opterr = 1, /* if error message should be printed */ | |||
optind = 1, /* index into parent argv vector */ | |||
optopt, /* character checked for validity */ | |||
optreset; /* reset getopt */ | |||
char *optarg; /* argument associated with option */ | |||
#define BADCH (int)'?' | |||
#define BADARG (int)':' | |||
@@ -64,17 +61,16 @@ char *optarg; /* argument associated with option */ | |||
* getopt -- | |||
* Parse argc/argv argument vector. | |||
*/ | |||
int | |||
getopt(nargc, nargv, ostr) | |||
int getopt(nargc, nargv, ostr) | |||
int nargc; | |||
char * const *nargv; | |||
const char *ostr; | |||
{ | |||
# define __progname nargv[0] | |||
static char *place = EMSG; /* option letter processing */ | |||
char *oli; /* option letter list index */ | |||
#define __progname nargv[0] | |||
static char *place = EMSG; /* option letter processing */ | |||
char *oli; /* option letter list index */ | |||
if (optreset || !*place) { /* update scanning pointer */ | |||
if (optreset || !*place) { /* update scanning pointer */ | |||
optreset = 0; | |||
if (optind >= nargc || *(place = nargv[optind]) != '-') { | |||
place = EMSG; | |||
@@ -85,7 +81,7 @@ getopt(nargc, nargv, ostr) | |||
place = EMSG; | |||
return (-1); | |||
} | |||
} /* option letter okay? */ | |||
} /* option letter okay? */ | |||
if ((optopt = (int)*place++) == (int)':' || | |||
!(oli = strchr(ostr, optopt))) { | |||
/* | |||
@@ -100,13 +96,13 @@ getopt(nargc, nargv, ostr) | |||
(void)printf("%s: illegal option -- %c\n", __progname, optopt); | |||
return (BADCH); | |||
} | |||
if (*++oli != ':') { /* don't need argument */ | |||
if (*++oli != ':') { /* don't need argument */ | |||
optarg = NULL; | |||
if (!*place) | |||
++optind; | |||
} | |||
else { /* need an argument */ | |||
if (*place) /* no white space */ | |||
else { /* need an argument */ | |||
if (*place) /* no white space */ | |||
optarg = place; | |||
else if (nargc <= ++optind) { /* no arg */ | |||
place = EMSG; | |||
@@ -117,12 +113,11 @@ getopt(nargc, nargv, ostr) | |||
"%s: option requires an argument -- %c\n", | |||
__progname, optopt); | |||
return (BADCH); | |||
} | |||
else /* white space */ | |||
} else /* white space */ | |||
optarg = nargv[optind]; | |||
place = EMSG; | |||
++optind; | |||
} | |||
return (optopt); /* dump back option letter */ | |||
return (optopt); /* dump back option letter */ | |||
} | |||
@@ -1,12 +1,12 @@ | |||
extern char * basename (const char *filename); | |||
extern char * basename(const char *filename); | |||
extern char *dirname(const char *path); | |||
extern int opterr, /* if error message should be printed */ | |||
optind , /* index into parent argv vector */ | |||
optopt, /* character checked for validity */ | |||
optreset; /* reset getopt */ | |||
extern char *optarg; /* argument associated with option */ | |||
extern int opterr, /* if error message should be printed */ | |||
optind , /* index into parent argv vector */ | |||
optopt, /* character checked for validity */ | |||
optreset; /* reset getopt */ | |||
extern char *optarg; /* argument associated with option */ | |||
extern int getopt(int nargc, char * const *nargv, const char *ostr); |