Added ZScale for scaling zmax

cbassa · Dec 16, 2024 · 65136df · 65136df
1 parent 6023636
commit 65136df
Show file tree

Hide file tree

Showing 7 changed files with 235 additions and 32 deletions.
diff --git a/Makefile.linux b/Makefile.linux
@@ -19,20 +19,20 @@ all:
 rffit: rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o rfsites.o rftles.o
 	gfortran -o rffit rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o rfsites.o rftles.o $(LFLAGS)
 
-rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o
-	gfortran -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o $(LFLAGS)
+rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o
+	gfortran -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o $(LFLAGS)
 
-rfedit: rfedit.o rfio.o rftime.o
-	$(CC) -o rfedit rfedit.o rfio.o rftime.o -lm
+rfedit: zscale.o rfedit.o rfio.o rftime.o
+	$(CC) -o rfedit rfedit.o zscale.o rfio.o rftime.o -lm
 
-rffind: rffind.o rfio.o rftime.o
-	$(CC) -o rffind rffind.o rfio.o rftime.o -lm
+rffind: rffind.o rfio.o rftime.o zscale.o
+	$(CC) -o rffind rffind.o rfio.o rftime.o zscale.o -lm
 
-rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o
-	$(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm
+rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o zscale.o
+	$(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o zscale.o -lm
 
-rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o
-	gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o $(LFLAGS)
+rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o
+	gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o $(LFLAGS)
 
 rffft: rffft.o rffft_internal.o rftime.o
 	$(CC) -o rffft rffft.o rffft_internal.o rftime.o -lfftw3f -lm -lsox

diff --git a/Makefile.osx b/Makefile.osx
@@ -33,20 +33,20 @@ all:
 rffit: rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o rfsites.o rftles.o
 	$(CC) -o rffit rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o rfsites.o rftles.o $(LFLAGS)
 
-rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o
-	$(CC) -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o $(LFLAGS)
+rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o
+	$(CC) -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o $(LFLAGS)
 
-rfedit: rfedit.o rfio.o rftime.o
-	$(CC) -o rfedit rfedit.o rfio.o rftime.o -lm
+rfedit: rfedit.o rfio.o rftime.o zscale.o
+	$(CC) -o rfedit rfedit.o rfio.o rftime.o zscale.o -lm
 
-rffind: rffind.o rfio.o rftime.o
-	$(CC) -o rffind rffind.o rfio.o rftime.o -lm
+rffind: rffind.o rfio.o rftime.o zscale.o
+	$(CC) -o rffind rffind.o rfio.o rftime.o zscale.o -lm
 
-rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o
-	$(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm
+rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o zscale.o
+	$(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o zscale.o -lm
 
-rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o
-	$(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o $(LFLAGS)
+rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o zscale.o
+	$(CC) -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o zscale.o $(LFLAGS)
 
 rffft: rffft.o rffft_internal.o rftime.o
 	$(CC) -o rffft rffft.o rffft_internal.o rftime.o -lfftw3f -lm -lsox $(LFLAGS)

diff --git a/makefile b/makefile
@@ -19,23 +19,23 @@ all:
 rffit: rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o rfsites.o rftles.o
 	gfortran -o rffit rffit.o sgdp4.o satutl.o deep.o ferror.o dsmin.o simplex.o versafit.o rfsites.o rftles.o $(LFLAGS)
 
-rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o
-	gfortran -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o $(LFLAGS)
+rfpng: rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o
+	gfortran -o rfpng rfpng.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o $(LFLAGS)
 
-rfdop: rfdop.o rftrace.o rfio.o rftime.o sgdp4.o satutl.o deep.o ferror.o rftles.o
-	$(CC) -o rfdop rfdop.o rftrace.o rfio.o rftime.o sgdp4.o satutl.o deep.o ferror.o rftles.o -lm
+rfdop: rfdop.o rftrace.o rfio.o rftime.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o
+	$(CC) -o rfdop rfdop.o rftrace.o rfio.o rftime.o sgdp4.o satutl.o deep.o ferror.o rftles.o zscale.o -lm
 
-rfedit: rfedit.o rfio.o rftime.o
-	$(CC) -o rfedit rfedit.o rfio.o rftime.o -lm
+rfedit: rfedit.o rfio.o rftime.o zscale.o
+	$(CC) -o rfedit rfedit.o rfio.o rftime.o zscale.o -lm
 
-rffind: rffind.o rfio.o rftime.o
-	$(CC) -o rffind rffind.o rfio.o rftime.o -lm
+rffind: rffind.o rfio.o rftime.o zscale.o
+	$(CC) -o rffind rffind.o rfio.o rftime.o zscale.o -lm
 
-rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o
-	$(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o -lm
+rftrack: rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o zscale.o
+	$(CC) -o rftrack rftrack.o rfio.o rftime.o rftrace.o sgdp4.o satutl.o deep.o ferror.o zscale.o -lm
 
-rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o
-	gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o $(LFLAGS)
+rfplot: rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o zscale.o
+	gfortran -o rfplot rfplot.o rftime.o rfio.o rftrace.o sgdp4.o satutl.o deep.o ferror.o versafit.o dsmin.o simplex.o rftles.o zscale.o $(LFLAGS)
 
 rffft: rffft.o rffft_internal.o rftime.o
 	$(CC) -o rffft rffft.o rffft_internal.o rftime.o -lfftw3f -lm -lsox

diff --git a/rfio.c b/rfio.c
@@ -4,6 +4,7 @@
 #include <math.h>
 #include "rftime.h"
 #include "rfio.h"
+#include "zscale.h"
 
 struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin,double foff)
 {
@@ -192,6 +193,11 @@ struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,dou
     }
   }
 
+  double z1,z2;
+  zscale(&s, s.nsub, 0.25,&z1, &z2);
+  printf("z1 = %f, z2 = %f\n", z1, z2);
+  s.zmin = z1;
+  s.zmax = z2;
   // Free 
   free(z);
   free(cz);

diff --git a/rfio.h b/rfio.h
@@ -1,3 +1,5 @@
+#ifndef RFIO_H
+#define RFIO_H
 struct spectrogram {
   int nsub,nchan;
   double *mjd;
@@ -10,3 +12,4 @@ struct spectrogram {
 struct spectrogram read_spectrogram(char *prefix,int isub,int nsub,double f0,double df0,int nbin,double foff);
 void write_spectrogram(struct spectrogram s,char *prefix);
 void free_spectrogram(struct spectrogram s);
+#endif
diff --git a/zscale.c b/zscale.c
@@ -0,0 +1,186 @@
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include "rftime.h"
+#include "rfio.h"
+
+#define MAX_REJECT 0.5
+#define MIN_NPIXELS 5
+#define GOOD_PIXEL 0
+#define BAD_PIXEL 1
+#define KREJ 2.5
+#define MAX_ITERATIONS 5
+
+void zsc_sample(struct spectrogram *image, int maxpix, double *samples, int *nsamples) {
+    int nc = image->nchan;
+    int nl = image->nsub;
+    double stride = fmax(1.0, sqrt((nc - 1) * (nl - 1) / (double)maxpix));
+    int istride = (int)stride;
+
+    int count = 0;
+    for (int i = 0; i < nc; i += istride) {
+        for (int j = 0; j < nl; j += istride) {
+            if (count >= maxpix) {
+                *nsamples = count;
+                return;
+            }
+            samples[count++] = image->z[i * nl + j];
+        }
+    }
+    *nsamples = count;
+}
+
+void zsc_compute_sigma(double *flat, int *badpix, int npix, int *ngoodpix, double *mean, double *sigma) {
+    double sumz = 0.0;
+    double sumsq = 0.0;
+    *ngoodpix = 0;
+
+    for (int i = 0; i < npix; i++) {
+        if (badpix[i] == GOOD_PIXEL) {
+            sumz += flat[i];
+            sumsq += flat[i] * flat[i];
+            (*ngoodpix)++;
+        }
+    }
+
+    if (*ngoodpix == 0) {
+        *mean = NAN;
+        *sigma = NAN;
+    } else if (*ngoodpix == 1) {
+        *mean = sumz;
+        *sigma = NAN;
+    } else {
+        *mean = sumz / *ngoodpix;
+        double temp = sumsq / (*ngoodpix - 1) - (sumz * sumz) / (*ngoodpix * (*ngoodpix - 1));
+        *sigma = temp < 0.0 ? 0.0 : sqrt(temp);
+    }
+}
+
+void zsc_fit_line(double *samples, int npix, double krej, int ngrow, int maxiter,
+                  int *ngoodpix_out, double *zstart, double *zslope) {
+    double xscale = 2.0 / (npix - 1);
+    double *xnorm = (double *)malloc(npix * sizeof(double));
+    int *badpix = (int *)calloc(npix, sizeof(int));
+
+    for (int i = 0; i < npix; i++) {
+        xnorm[i] = i * xscale - 1.0;
+    }
+
+    int ngoodpix = npix;
+    int minpix = fmax(MIN_NPIXELS, (int)(npix * MAX_REJECT));
+    int last_ngoodpix = npix + 1;
+
+    double intercept = 0.0;
+    double slope = 0.0;
+
+    for (int niter = 0; niter < maxiter; niter++) {
+        if (ngoodpix >= last_ngoodpix || ngoodpix < minpix) {
+            break;
+        }
+
+        double sumx = 0.0, sumxx = 0.0, sumxy = 0.0, sumy = 0.0;
+        int count = 0;
+        for (int i = 0; i < npix; i++) {
+            if (badpix[i] == GOOD_PIXEL) {
+                sumx += xnorm[i];
+                sumxx += xnorm[i] * xnorm[i];
+                sumxy += xnorm[i] * samples[i];
+                sumy += samples[i];
+                count++;
+            }
+        }
+
+        double delta = count * sumxx - sumx * sumx;
+        intercept = (sumxx * sumy - sumx * sumxy) / delta;
+        slope = (count * sumxy - sumx * sumy) / delta;
+
+        double *fitted = (double *)malloc(npix * sizeof(double));
+        double *flat = (double *)malloc(npix * sizeof(double));
+        for (int i = 0; i < npix; i++) {
+            fitted[i] = xnorm[i] * slope + intercept;
+            flat[i] = samples[i] - fitted[i];
+        }
+
+        double mean, sigma;
+        zsc_compute_sigma(flat, badpix, npix, &ngoodpix, &mean, &sigma);
+
+        double threshold = sigma * krej;
+        double lcut = -threshold;
+        double hcut = threshold;
+
+        for (int i = 0; i < npix; i++) {
+            if (flat[i] < lcut || flat[i] > hcut) {
+                badpix[i] = BAD_PIXEL;
+            }
+        }
+
+        for (int i = 0; i < npix; i++) {
+            if (badpix[i] == BAD_PIXEL) {
+                for (int j = -ngrow; j <= ngrow; j++) {
+                    int idx = i + j;
+                    if (idx >= 0 && idx < npix) {
+                        badpix[idx] = BAD_PIXEL;
+                    }
+                }
+            }
+        }
+
+        free(fitted);
+        free(flat);
+
+        last_ngoodpix = ngoodpix;
+        ngoodpix = 0;
+        for (int i = 0; i < npix; i++) {
+            if (badpix[i] == GOOD_PIXEL) {
+                ngoodpix++;
+            }
+        }
+    }
+
+    *zstart = intercept - slope;
+    *zslope = slope * xscale;
+
+    free(xnorm);
+    free(badpix);
+}
+
+int compare_doubles(const void *a, const void *b) {
+    double diff = (*(double *)a - *(double *)b);
+    return (diff > 0) - (diff < 0);
+}
+
+void zscale(struct spectrogram *image, int nsamples, double contrast, double *z1, double *z2) {
+    double *samples = (double *)malloc(nsamples * sizeof(double));
+    int npix;
+
+    zsc_sample(image, nsamples, samples, &npix);
+
+    qsort(samples, npix, sizeof(double), compare_doubles);
+
+    double zmin = samples[0];
+    double zmax = samples[npix - 1];
+    double median;
+
+    if (npix % 2 == 1) {
+        median = samples[npix / 2];
+    } else {
+        median = 0.5 * (samples[npix / 2 - 1] + samples[npix / 2]);
+    }
+
+    int ngoodpix;
+    double zstart, zslope;
+    zsc_fit_line(samples, npix, KREJ, fmax(1, npix * 0.01), MAX_ITERATIONS,
+                 &ngoodpix, &zstart, &zslope);
+
+    if (ngoodpix < fmax(MIN_NPIXELS, npix * MAX_REJECT)) {
+        *z1 = zmin;
+        *z2 = zmax;
+    } else {
+        if (contrast > 0) zslope /= contrast;
+        *z1 = fmax(zmin, median - (npix / 2 - 1) * zslope);
+        *z2 = fmin(zmax, median + (npix - npix / 2) * zslope);
+    }
+
+    free(samples);
+}
diff --git a/zscale.h b/zscale.h
@@ -0,0 +1,8 @@
+#ifndef ZSCALE_H
+#define ZSCALE_H
+
+#include "rfio.h"
+
+void zscale(struct spectrogram *image, int nsamples, double contrast, double *z1, double *z2);
+
+#endif