00001
00002 #include <drms.h>
00003 #include "lev0lev1.h"
00004
00005 #define XXX_BAD_GEOMETRY -1
00006 #define XXX_BAD_PIXCOUNT -2
00007
00008 #define NBINS 1048576
00009 #define MINOUT -256
00010 #define MAXOUT 1048320
00011
00012 static int hist[NBINS];
00013 static float dark2[2270912];
00014 static float flat2[2270912];
00015
00017 int do_flat_iris(LEV0LEV1 *info)
00018 {
00019 int is_dark = info->darkflag;
00020 float *flat = info->adataff;
00021 float *dark = info->adatadark;
00022 short *in = info->adata0;
00023 int *out = info->dat1.adata1A;
00024 int datavals = info->datavals;
00025 int missvals = info->missvals;
00026 int totalvals = datavals + missvals;
00027 int nx = info->nx;
00028 int ny = info->ny;
00029 int sx = info->sumx;
00030 int sy = info->sumy;
00031 int totalpix = nx * ny;
00032 int ndark, nflat;
00033 int min, max, medn;
00034 short tmp;
00035 int itmp;
00036 float ftmp,flatsum,darksum;
00037 double dtmp, s, s2, s3, s4, ss;
00038 int i, j, k, l;
00039 int idx, IDX;
00040
00041
00042 if (sx < 1 || sy < 1 || (nx*sx != 4144 && nx*sx != 2072) || ny*sy != 1096)
00043 return XXX_BAD_GEOMETRY;
00044 if (totalpix < totalvals)
00045 return XXX_BAD_PIXCOUNT;
00046
00047
00048
00049
00050 if (sx > 1 || sy > 1) {
00051 for (j=0; j<ny; ++j) {
00052 for (i=0; i<nx; ++i) {
00053 idx = j*nx + i;
00054 nflat = ndark = sx*sy;
00055 flatsum = darksum = 0;
00056 for (l=j*sy; l<(j+1)*sy; ++l) {
00057 for (k=i*sx; k<(i+1)*sx; ++k) {
00058 IDX = l*nx*sx + k;
00059 if (isnan(dark[IDX]))
00060 --ndark;
00061 else
00062 darksum += dark[IDX];
00063 if (isnan(flat[IDX]))
00064 --nflat;
00065 else
00066 flatsum += flat[IDX];
00067 }
00068 }
00069 dark2[idx] = ndark ? darksum/ndark : DRMS_MISSING_FLOAT;
00070 flat2[idx] = nflat ? flatsum/nflat : DRMS_MISSING_FLOAT;
00071 }
00072 }
00073 }
00074
00075
00076
00077
00078 memset(hist, 0, 4*NBINS);
00079 datavals = 0;
00080 for (i=0; i<totalpix; ++i) {
00081 tmp = in[i];
00082 if (DRMS_MISSING_SHORT == tmp) {
00083 out[i] = DRMS_MISSING_INT;
00084 } else {
00085 if (sx == 1 && sy == 1)
00086 ftmp = is_dark ? tmp : (tmp - dark[i]) / flat[i];
00087 else
00088 ftmp = is_dark ? tmp : (tmp - dark2[i]) / flat2[i];
00089 if (ftmp >= MINOUT && ftmp < MAXOUT) {
00090 itmp = roundf(ftmp);
00091 out[i] = itmp;
00092 ++hist[itmp-MINOUT];
00093 ++datavals;
00094 } else {
00095 out[i] = DRMS_MISSING_INT;
00096 }
00097 }
00098 }
00099 info->datavals = datavals;
00100 info->missvals = totalvals - datavals;
00101
00102
00103
00104
00105 if (datavals) {
00106 switch(info->winflip) {
00107 case 1:
00108 {
00109 int tmp[4144];
00110 for (j=0; j<ny/2; ++j) {
00111 memcpy(tmp, &out[j*nx], 4*nx);
00112 memcpy(&out[j*nx], &out[(ny-j-1)*nx], 4*nx);
00113 memcpy(&out[(ny-j-1)*nx], tmp, 4*nx);
00114 }
00115 }
00116 break;
00117 case 2:
00118 {
00119 int tmp2;
00120 for (j=0; j<ny; ++j)
00121 for (i=0; i<nx/2; ++i) {
00122 tmp2 = out[j*nx+i];
00123 out[j*nx+i] = out[j*nx+(nx-1-i)];
00124 out[j*nx+(nx-1-i)] = tmp2;
00125 }
00126 }
00127 break;
00128 case 3:
00129 {
00130 int tmp[4144], tmp2;
00131 for (j=0; j<ny/2; ++j) {
00132 memcpy(tmp, &out[j*nx], 4*nx);
00133 memcpy(&out[j*nx], &out[(ny-j-1)*nx], 4*nx);
00134 memcpy(&out[(ny-j-1)*nx], tmp, 4*nx);
00135 }
00136 for (j=0; j<ny; ++j)
00137 for (i=0; i<nx/2; ++i) {
00138 tmp2 = out[j*nx+i];
00139 out[j*nx+i] = out[j*nx+(nx-1-i)];
00140 out[j*nx+(nx-1-i)] = tmp2;
00141 }
00142 }
00143 break;
00144 }
00145 }
00146
00147
00148
00149
00150
00151 info->datamin = info->datamax = info->datamedn = DRMS_MISSING_INT;
00152 info->datamean = info->data_rms = info->dataskew = info->datakurt = DRMS_MISSING_DOUBLE;
00153
00154 min = 0;
00155 max = NBINS-1;
00156 if (datavals) {
00157 while (hist[min] == 0)
00158 ++min;
00159 info->datamin = min + MINOUT;
00160
00161 while (hist[max] == 0)
00162 --max;
00163 info->datamax = max + MINOUT;
00164
00165 medn = min;
00166 i = hist[medn];
00167 while (2*i < datavals)
00168 i += hist[++medn];
00169 info->datamedn = medn + MINOUT;
00170
00171 s = s2 = s3 = s4 = 0.0;
00172 for (i=min; i<=max; ++i) {
00173 double ii = i + MINOUT;
00174 s += (dtmp = ii*hist[i]);
00175 s2 += (dtmp *= ii);
00176 s3 += (dtmp *= ii);
00177 s4 += (dtmp *= ii);
00178 }
00179 s /= datavals;
00180 info->datamean = s;
00181 ss = s*s;
00182 s2 /= datavals;
00183 s3 /= datavals;
00184 s4 /= datavals;
00185 if (datavals > 1) {
00186 dtmp = datavals * (s2-ss) / (datavals-1);
00187 info->data_rms = sqrt(dtmp);
00188 if (dtmp > 0.0) {
00189 info->dataskew = (s3 - s * (3*s2 - 2*ss)) /
00190 (dtmp*info->data_rms);
00191 info->datakurt = (s4 - 4*s*s3 + 3*ss*(2*s2-ss)) /
00192 (dtmp*dtmp) - 3;
00193 }
00194 }
00195
00196 {
00197 int n=0, nsatpix=0, sum=0, dp[8];
00198 DRMS_Record_t *rs = info->rs1;
00199 dp[0] = 1; dp[1] = 10; dp[2] = 25; dp[3] = 75;
00200 dp[4] = 90; dp[5] = 95; dp[6] = 98; dp[7] = 99;
00201 for (i=min; i<=max; ++i) {
00202 sum += hist[i];
00203 if (sum > dp[n]*datavals/100) switch (n) {
00204 case 0:
00205 drms_setkey_float(rs, "DATAP01", 1.0 + i + MINOUT);
00206 if (sum > dp[n+1]*datavals/100) n++;
00207 case 1:
00208 drms_setkey_float(rs, "DATAP10", 1.0 + i + MINOUT);
00209 if (sum > dp[n+1]*datavals/100) n++;
00210 case 2:
00211 drms_setkey_float(rs, "DATAP25", 1.0 + i + MINOUT);
00212 if (sum > dp[n+1]*datavals/100) n++;
00213 case 3:
00214 drms_setkey_float(rs, "DATAP75", 1.0 + i + MINOUT);
00215 if (sum > dp[n+1]*datavals/100) n++;
00216 case 4:
00217 drms_setkey_float(rs, "DATAP90", 1.0 + i + MINOUT);
00218 if (sum > dp[n+1]*datavals/100) n++;
00219 case 5:
00220 drms_setkey_float(rs, "DATAP95", 1.0 + i + MINOUT);
00221 if (sum > dp[n+1]*datavals/100) n++;
00222 case 6:
00223 drms_setkey_float(rs, "DATAP98", 1.0 + i + MINOUT);
00224 if (sum > dp[n+1]*datavals/100) n++;
00225 case 7:
00226 drms_setkey_float(rs, "DATAP99", 1.0 + i + MINOUT);
00227 n++;
00228 break;
00229 }
00230 if (i>max) n++;
00231 if (i > (15000 - MINOUT)) nsatpix += hist[i];
00232 }
00233 drms_setkey_int(rs, "NSATPIX", nsatpix);
00234 }
00235 }
00236 return 0;
00237 }
00238
