Index: /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Make.Common =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Make.Common (revision 24390) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Make.Common (revision 24390) @@ -0,0 +1,1 @@ +link ../Make.Common Index: /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Makefile =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Makefile (revision 24390) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/Makefile (revision 24390) @@ -0,0 +1,16 @@ +include ../Make.Common + +VERSION=0.00 +CCWARN+=$(WERROR) +CCDEFS+=$(VERSIONDEFS) +$(EXECNAME): $(OBJS) libfh.a + +include ../Make.Common + +# Dependencies by Make.Common $Revision: 2.17 $ + +$(OBJ)/sort.o: sort.c + +$(OBJ)/tiltystreak.o: tiltystreak.c fh/fh.h \ + fhreg/general.h fhreg/macros.h \ + fhreg/gpc_detector.h Index: /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/linearfit.c =================================================================== --- /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/linearfit.c (revision 24390) +++ /trunk/extsrc/gpcsw/gpcsrc/fits/tiltystreak/linearfit.c (revision 24390) @@ -0,0 +1,129 @@ +/* Routine to fit a linear function */ +/* 060304 v1.2 add wlinearfit() */ +/* 060130 v1.1 replaced static allocations with calloc's */ +/* 040208 v1.0 John Tonry (derived from fortran) */ +#include +#include +#include + +#define ABS(x) ((x)<0?(-(x)):(x)) + +static int linsolve(int n, double *y, double *a, double *x); + +int linearfit(int npt, double *y, int nx, double *x, double *param) +{ +/* Y = NPT data being fitted as a linear function of X[*][NX] + * PARAM = returned parameters of fit: + * Y[*] = PARAM[0]*X[*][0], + * + PARAM[1]*X[*][1], + * ... + * + PARAM[NX-1]*X[*][NX-1] + * + */ + int i, j, k; + double *v, *a; + v = (double *)calloc(nx, sizeof(double)); + a = (double *)calloc(nx*nx, sizeof(double)); + + if(npt < nx) { + fprintf(stderr, "linearfit: Too few points provided: %d %d\n", npt, nx); + return(-2); + } + +#if 0 /* Show us data to solve? */ + for(i=0; i rat) { + k = i; + rat = ABS(a[j+k*n]); + } + } + rowstatus[k] = 1; + row[j] = k; + if(rat == 0.0) { + fprintf(stderr, "WHOA: singular matrix, iter %d, row %d\n", j, k); + free(rowstatus); + free(row); + return(1); + } +/* Subtract away matrix below diagonal */ + for(i=0; i=0; j--) { + x[j] = y[row[j]]; + for(k=j+1; k +#include +#include + +#define MAXSTACK 256 +//#define NSTOP 15 +#define NSTOP 8 + +/* Quicksort program 1980 (John Tonry) */ +int qsort_dbl(int n, double *x) +{ + double key, kl, kr, km, temp; + int l, r, m, lstack[MAXSTACK], rstack[MAXSTACK], sp; + register int i, j, k; + int mgtl, lgtr, rgtm; + + sp = 0; + lstack[sp] = 0; + rstack[sp] = n-1; + + while(n > NSTOP && sp >= 0) { +/* Sort a subrecord off the stack */ + l = lstack[sp]; + r = rstack[sp]; + sp--; + m = (l + r) / 2; +/* Set KEY = median of X(L), X(M), X(R) */ + kl = x[l]; + km = x[m]; + kr = x[r]; + mgtl = km > kl; + rgtm = kr > km; + lgtr = kl > kr; +#ifndef STD_KEY + if(mgtl ^ rgtm) { + if(mgtl ^ lgtr) key = kr; + else key = kl; + } else { + key = km; + } +#else +/* Curiously enough, this non-median seems to work as well or better */ + if(mgtl ^ rgtm) { + key = km; + } else { + if(mgtl ^ lgtr) key = kl; + else key = kr; + } +#endif + i = l; + j = r; + while(1) { +/* Find a big record on the left */ + while(x[i] < key) i++; + +/* Find a small record on the right */ + while(x[j] > key) j--; + + if(i >= j) break; +/* Exchange records */ + temp = x[i]; + x[i] = x[j]; + x[j] = temp; + i++; + j--; + } + +/* Subfile is partitioned into two halves, left .le. right */ +/* Push the two halves on the stack */ + if(j-l+1 > NSTOP) { + sp = sp + 1; + lstack[sp] = l; + rstack[sp] = j; + } + if(r-j > NSTOP) { + sp = sp + 1; + lstack[sp] = j+1; + rstack[sp] = r; + } + if(sp >= MAXSTACK) { +// fprintf(stderr,"QSORT8: Fatal error from stack overflow\n"); +// fprintf(stderr,"Fall back on sort by insertion\n"); + break; + } + } + +/* Sorting routine that sorts the N elements of single precision */ +/* array X by straight insertion between previously sorted numbers */ + for(j=n-2; j>=0; j--) { + for(i=j+1, k=j; i NSTOP && sp >= 0) { +/* Sort a subrecord off the stack */ + l = lstack[sp]; + r = rstack[sp]; + sp--; + m = (l + r) / 2; +/* Set KEY = median of X(L), X(M), X(R) */ + kl = x[l]; + km = x[m]; + kr = x[r]; + mgtl = km > kl; + rgtm = kr > km; + lgtr = kl > kr; +#ifndef STD_KEY + if(mgtl ^ rgtm) { + if(mgtl ^ lgtr) key = kr; + else key = kl; + } else { + key = km; + } +#else +/* Curiously enough, this non-median seems to work as well or better */ + if(mgtl ^ rgtm) { + key = km; + } else { + if(mgtl ^ lgtr) key = kl; + else key = kr; + } +#endif + i = l; + j = r; + while(1) { +/* Find a big record on the left */ + while(x[i] < key) i++; + +/* Find a small record on the right */ + while(x[j] > key) j--; + + if(i >= j) break; +/* Exchange records */ + temp = x[i]; + x[i] = x[j]; + x[j] = temp; + i++; + j--; + } + +/* Subfile is partitioned into two halves, left .le. right */ +/* Push the two halves on the stack */ + if(j-l+1 > NSTOP) { + sp = sp + 1; + lstack[sp] = l; + rstack[sp] = j; + } + if(r-j > NSTOP) { + sp = sp + 1; + lstack[sp] = j+1; + rstack[sp] = r; + } + if(sp >= MAXSTACK) { +// fprintf(stderr,"QSORT8: Fatal error from stack overflow\n"); +// fprintf(stderr,"Fall back on sort by insertion\n"); + break; + } + } + +/* Sorting routine that sorts the N elements of single precision */ +/* array X by straight insertion between previously sorted numbers */ + for(j=n-2; j>=0; j--) { + for(i=j+1, k=j; i +#include +#include + +#include "math.h" +#include "fh/fh.h" +#include "fhreg/general.h" +#include "fhreg/gpc_detector.h" + +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define NINT(x) (x<0?(int)((x)-0.5):(int)((x)+0.5)) +#define ABS(a) (((a) > 0) ? (a) : -(a)) + +#define NODATA 0 /* Marker for *No Data* */ +#define NCELL 8 /* Number of cells in a row */ +#define SAT4SURE 45000 /* Blasted, we don't care about pixels above */ +#define NPOLY 3 /* N terms in of flattening polynomial */ + +//#define DEBUG /* Turn on debug output? */ + +/* Info for a cell */ +typedef struct cell_info { + HeaderUnit hu; /* Header unit for this cell */ + int cell; /* Cell number */ + int mefxt; /* Where is it in the MEF? */ + int bias; /* Bias level */ + int sky; /* Sky level */ + int rms; /* RMS in the sky */ + int nsat; /* Saturated pixel count */ + double ampl; /* Streak amplitude */ + short *data; /* Image data */ + int *mask; /* Mask bits */ +} CELL; + +/* Prototypes */ +static int cell_stats(int sx, int nx, int ny, int NX, int NY, int bzero, + CELL *cell); +static int int_median(int n, int *key); +int qsort_dbl(int n, double *x); +int qsort_int(int n, int *x); +static int streak_ratio(int sx, int nx, int ny, int NX, CELL otarow[]); +static int linear_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]); +static int bilinear_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]); +static int parabola_streak_fit(int sx, int nx, int ny, int NX, CELL otarow[]); +int linearfit(int npt, double *y, int nx, double *x, double *param); + +static int nmedian_buf=0; +static int *median_buf=NULL; +static double *Abuf=NULL, *Bbuf=NULL, *Xbuf=NULL; + +int main(int argc, const char* argv[]) +{ + int i, row, col, ext; + HeaderUnit ihu = fh_create(); + HeaderUnit ehu = NULL; + const char* ifilename; + char extname[FH_MAX_STRLEN+1], xtension[FH_MAX_STRLEN+1], cellsought[5]; + int nextend; + CELL otarow[NCELL]; + int naxis, naxis1, naxis2, naxis3; + int prescan1, ovrscan1, ovrscan2; + double bzero_d; + int bzero = 32768, firstime = 1, algo=2; + + if(argc < 2) { + fprintf(stderr, "error: syntax: tiltystreak Uncompressed_GPC1_MEF.fits\n"); + exit(EXIT_FAILURE); + } + ifilename = argv[1]; + +/* Choose different algorithms? */ + if(argc > 2 && strncmp(argv[2], "alg=", 4) == 0) { + if(strcmp(argv[2]+4, "line") == 0) { + algo = 1; + } else if(strcmp(argv[2]+4, "parabola") == 0) { + algo = 2; + } else if(strcmp(argv[2]+4, "bent") == 0) { + algo = 3; + } else { + fprintf(stderr, "error: `%s' is not a valid algorithm\n", argv[2]); + exit(EXIT_FAILURE); + } + } + + if (fh_file(ihu, ifilename, FH_FILE_RDWR) != FH_SUCCESS) { + fprintf(stderr, "error: tiltystreak could not open file `%s'\n", + ifilename); + exit(EXIT_FAILURE); + } + nextend = fh_extensions(ihu); + if (nextend < 1) { + fprintf(stderr, "error: `%s' is not a multi-extension FITS\n", + ifilename); + exit(EXIT_FAILURE); + } + +/* Loop over all rows of cells */ + for(row=0; row nextend) { + fprintf(stderr, + "error: Could not find extension `%s'\n", cellsought); + exit(EXIT_FAILURE); + } + +/* Positioned at the desired cell now */ + naxis3 = 1; + if (fh_get_PRESCAN1(ehu, &prescan1) != FH_SUCCESS || + fh_get_OVRSCAN1(ehu, &ovrscan1) != FH_SUCCESS || + fh_get_OVRSCAN2(ehu, &ovrscan2) != FH_SUCCESS || + fh_get_NAXIS(ehu, &naxis) != FH_SUCCESS || + fh_get_NAXIS1(ehu, &naxis1) != FH_SUCCESS || + fh_get_NAXIS2(ehu, &naxis2) != FH_SUCCESS) { + fprintf(stderr, "error: Cannot get NAXIS* in extension %d row %d col %d\n", ext, row, col); + exit(EXIT_FAILURE); + } +/* Cell isn't an image or is 32bpp */ + if(naxis != 2) { + otarow[col].mefxt = -1; + continue; + } + +/* Make some space now that we finally know the size */ + if(firstime) { + firstime = 0; + for(i=0; i 0) i++; + if(i == 0) continue; /* Nope, entirely wiped out */ + +/* Process this row of cells */ + +/* What's the stats in each cell, generate mask */ + for(col=0; colbias = 0; + } else { + for(j=1; jdata[sx+nx+(NX-nx-sx)/2+j*NX] + bzero; + cell->bias = int_median(ny-2, median_buf); + } + +/* Get sky stats */ + for(k=n=0; kdata[sx+i+NX*j]+bzero < SAT4SURE && + cell->data[sx+i+NX*j]+bzero != NODATA) { + median_buf[n++] = cell->data[sx+i+NX*j] + bzero; + } + } + + if(n < 20) { /* Better have hit at least 20! */ + cell->sky = cell->rms = 0.0; + return(-1); + } + +/* First pass at sky and quartile */ + cell->sky = int_median(n, median_buf); + cell->rms = 1.33*(cell->sky - median_buf[n/4]); + +/* Clip at 3 sigma */ + for(j=0; median_buf[j] < cell->sky - Q_CLIP*cell->rms && j cell->sky + Q_CLIP*cell->rms && k>n/2; k--); + cell->sky = median_buf[(j+k)/2]; + cell->rms = 1.33 * (cell->sky - median_buf[(3*j+k)/4]); + cell->sky -= cell->bias; + +/* Create mask */ + thresh = cell->bias + cell->sky + MASK_CLIP*cell->rms; + cell->nsat = 0; + for(j=0; jmask[i+NX*j] = 1; + if(cell->data[i+NX*j]+bzero > thresh) cell->mask[i+NX*j] = 0; + if(cell->data[i+NX*j]+bzero > SAT4SURE) (cell->nsat)++; + } + } + + return(0); +} + +static int int_median(int n, int *key) +{ + if(n == 0) return(0); + qsort_int(n, key); + return((key[n/2]+key[(n-1)/2]) / 2); +} + +/* Determine the relative streakiness of all cells in this row */ +/* + * Compare each cell's pixels in the difference which are greater + * than MINSTREAK with the average using pixels surviving the + * intersection of all the masks + */ +#define MINSTREAK 50 +#define EDGE 3 +static int streak_ratio(int sx, int nx, int ny, int NX, CELL otarow[]) +{ + int i, j, k, nalive=0, ncontrib; + double diff; + + for(k=0; k 0) { + nalive++; + } else { + otarow[k].ampl = 0.0; + } + } + +/* Loop over all the pixels */ + ncontrib = 0; + for(j=1+EDGE; j MINSTREAK)) { + for(k=0; k 0) { + k = int_median(n, medimg) + Abuf[j-1] * (i-nx*sx); + sy += k; + s1 += i - (sx+nx); + } + } + Abuf[j] = 0.0; + if(s1 != 0) Abuf[j] = sy / s1; + for(k=0, n=1; k 0) { + k = int_median(n, medimg); + if(i < sx+nx/3) { + s1 += k + Abuf[j-1]; + n1++; + } else { + s2 += k + Bbuf[j-1]; + n2++; + } + } + } + Abuf[j] = Bbuf[j] = 0.0; + if(n1 > 0) Abuf[j] = s1/n1; + if(n2 > 0) Bbuf[j] = s2/n2; + for(k=0, n=1; k 0) { + k = int_median(n, medimg); + sy += k; + s1 += i - (sx+nx); + } + } + if(s1 != 0) A = sy / s1; + +/* Fit a full line z = B*(i-nx-sx)+C to the left third */ + s0 = s1 = s2 = sm = B = C = 0.0; + npt = 0; + for(i=sx+EDGE; i 0) { + k = int_median(n, medimg); + s1 += k; + s0 += (i-sx-nx); + s2 += (i-sx-nx)*(i-sx-nx); + sm += (i-sx-nx)*k; + npt++; + } + } + if(npt > 0) { + B = (npt*sm-s0*s1) / (npt*s2-s0*s0); + C = (s2*s1 - s0*sm) / (npt*s2-s0*s0); + } + + if(A != B) { + i0 = C / (A-B) + sx+nx; + } else { + i0 = 0; + } + +#ifdef DEBUG + if(otarow[3].mefxt == 4) printf("%4d %4d %9.3f %9.3f %9.3f\n", + j, i0, A, B, C); +#endif + +/* Correct all the cells */ + for(k=0; k