JSOC_Documentation: proj/lev1.5_hmi/apps/lev1

00001 #include <jsoc_main.h>
00002 #include <mkl.h>
00003 #include <fresize.h>
00004 #include <gapfill.h>
00005 #include "limb_fit.h"
00006 
00007 
00008 void warn(const char *fmt, ...)
00009 {
00010     va_list ap;
00011     va_start(ap, fmt);
00012     vfprintf(stderr, fmt, ap);
00013     va_end(ap);
00014 }
00015 
00016 void die(const char *fmt, ...)
00017 {
00018     va_list ap;
00019     va_start(ap, fmt);
00020     vfprintf(stderr, fmt, ap);
00021     va_end(ap);
00022     exit(1);
00023 }
00024 
00025 void ck(MKL_LONG status)
00026 {
00027     if (DftiErrorClass(status, DFTI_NO_ERROR))
00028     return;
00029     die(DftiErrorMessage(status));
00030 }
00031 
00032 //CORRECTION OF HEIGHT FORMATION (copied from ingest_dcon.c)
00033 //returns 0 if corrections were successful, 1 otherwise
00034 int heightformation(int FID, double OBSVR, float *CDELT1, float *RSUN, float *CRPIX1, float *CRPIX2, float CROTA2)
00035 {
00036   int wl=0;
00037   int status=0;
00038   float correction=0.0,correction2=0.0;
00039   
00040   wl = (FID/10)%20;  //temp is now the filter index
00041   
00042   if( (wl >= 0) && (wl < 20) )
00043     {
00044       correction  = 0.445*exp(-(wl-10.-(float)OBSVR/(0.690/6173.*3.e8/20.)-0.25)*(wl-10.-(float)OBSVR/(0.690/6173.*3.e8/20.)-0.25)/7.1);
00045       correction2 = 0.39*(-2.0*(wl-10.- (float)OBSVR/(0.690/6173.*3.e8/20.)-0.35)/6.15)*exp(-(wl-10.-(float)OBSVR/(0.690/6173.*3.e8/20.)-0.35)*(wl-10.-(float)OBSVR/(0.690/6173.*3.e8/20.)-0.35)/6.15);
00046       //printf("A CDELT1= %f CROTA2= %f OBS_VR= %f RSUN= %f correction= %f \n",*CDELT1,CROTA2,OBSVR,*RSUN,correction);
00047       *CDELT1 = *CDELT1*(*RSUN)/((*RSUN)-correction);
00048       //printf("B CDELT1= %f CROTA2= %f OBS_VR= %f \n",*CDELT1,CROTA2,OBSVR);
00049       *RSUN   = *RSUN-correction;
00050       *CRPIX1 = *CRPIX1-cos(M_PI-CROTA2*M_PI/180.)*correction2;
00051       *CRPIX2 = *CRPIX2-sin(M_PI-CROTA2*M_PI/180.)*correction2;
00052     }
00053   else status=1;
00054 
00055   return status;
00056 }
00057 
00058 char *module_name = "lev1_dcon";
00059 ModuleArgs_t module_args[] = {
00060     {ARG_STRING, "in", "", "input query"},
00061     {ARG_STRING, "out", "", "output series"},
00062     {ARG_STRING, "psf", "", "PSF query"},
00063     {ARG_INT, "iter", "25", "number of R-L iterations"},
00064     {ARG_END}
00065 };
00066 
00067 int DoIt() {
00068     const char *in = cmdparams_get_str(&cmdparams, "in", NULL);
00069     const char *out = cmdparams_get_str(&cmdparams, "out", NULL);
00070     const char *psf = cmdparams_get_str(&cmdparams, "psf", NULL);
00071     char source[100];
00072     int iter = cmdparams_get_int(&cmdparams, "iter", NULL);
00073     int status;
00074 
00075     drms_series_exists(drms_env, out, &status);
00076     if (DRMS_SUCCESS != status)
00077         die("Output series %s does not exit\n", out);
00078 
00079     // open and stage input records
00080     DRMS_RecordSet_t *rsin = drms_open_records(drms_env, in, &status);
00081     if (status || !rsin)
00082         die("Can't do drms_open_records(%s)\n", in);
00083     else
00084         warn("DRMS recordset %s opened\n", in);
00085     if (drms_stage_records(rsin, 1, 0) != DRMS_SUCCESS)
00086     {
00087         drms_close_records(rsin, DRMS_FREE_RECORD);
00088         die("Failure staging DRMS records (%s).\n", in);
00089     }
00090 
00091     int nrecs = rsin->n;
00092     if (!nrecs)
00093         die("No records matching %s in found\n", in);
00094     else
00095         warn("%d records found\n", nrecs);
00096 
00097     // read PSF
00098     DRMS_RecordSet_t *rpsf = drms_open_records(drms_env, psf, &status);
00099     if (status || !rpsf)
00100     die("Can't do drms_open_records(%s)\n", psf);
00101     DRMS_Segment_t *segpsf = drms_segment_lookup(rpsf->records[0], "psf");
00102     if (!segpsf)
00103     die("No psf segment found.\n");
00104     DRMS_Array_t *arrpsf = drms_segment_read(segpsf, DRMS_TYPE_DOUBLE, &status);
00105     if (status || !arrpsf)
00106     die("Can't read psf record %s\n", psf);
00107     else if (arrpsf->naxis != 2 || arrpsf->axis[0] != 4096 || arrpsf->axis[1] != 4096 || arrpsf->type != DRMS_TYPE_DOUBLE)
00108     die("PSF record %s does not contain 4096x4096 double precision array\n", psf);
00109     warn("PSF record %s read\n", psf);
00110 
00111     // pad PSF data array for FFT
00112     float *P = MKL_malloc(sizeof(float)*4096*4098, 64);
00113     if (!P)
00114         die("Can't allocate memory for array P\n");
00115     for (int j=0; j<4096; ++j)
00116         for (int i=0; i<4096; ++i)
00117             P[4098*j+i] = ((double *)arrpsf->data)[4096*j+i];
00118     drms_free_array(arrpsf);
00119     drms_close_records(rpsf, DRMS_FREE_RECORD);
00120 
00121     // set up MKL DFTI stuff
00122     MKL_LONG lengths[2], strdfor[3], strdbak[3];
00123     DFTI_DESCRIPTOR *p, *q;
00124     lengths[0] = lengths[1] = 4096;
00125     strdfor[0] = 0; strdfor[1] = 4098; strdfor[2] = 1;
00126     strdbak[0] = 0; strdbak[1] = 2049; strdbak[2] = 1;
00127     ck(DftiCreateDescriptor(&p, DFTI_SINGLE, DFTI_REAL, 2, lengths));
00128     ck(DftiSetValue(p, DFTI_BACKWARD_SCALE, 1.0/(4096*4096)));
00129     ck(DftiSetValue(p, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX));
00130     ck(DftiSetValue(p, DFTI_INPUT_STRIDES, strdfor));
00131     ck(DftiSetValue(p, DFTI_OUTPUT_STRIDES, strdbak));
00132     ck(DftiCommitDescriptor(p));
00133     ck(DftiCreateDescriptor(&q, DFTI_SINGLE, DFTI_REAL, 2, lengths));
00134     ck(DftiSetValue(q, DFTI_BACKWARD_SCALE, 1.0/(4096*4096)));
00135     ck(DftiSetValue(q, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX));
00136     ck(DftiSetValue(q, DFTI_INPUT_STRIDES, strdbak));
00137     ck(DftiSetValue(q, DFTI_OUTPUT_STRIDES, strdfor));
00138     ck(DftiCommitDescriptor(q));
00139 
00140     // forward transform P
00141     ck(DftiComputeForward(p, P));
00142 
00143     // prepare other work arrays
00144     // O - original image (with all MISSING values replaced by zeors)
00145     // E - current estimate
00146     // Esav - copy of E
00147     // B - Blurred estimate
00148     // R - image divided by blurred estimate (shares same location with B)
00149     // C - correction for next iteraion (shares same location with B)
00150     float *O = MKL_malloc(sizeof(float)*4096*4098, 64);
00151     float *E = MKL_malloc(sizeof(float)*4096*4098, 64);
00152     float *Esav = MKL_malloc(sizeof(float)*4096*4098, 64);
00153     float *R = MKL_malloc(sizeof(float)*4096*4098, 64);
00154     float *B = R;
00155     float *C = R;
00156     if (!O || !E || !Esav || !R)
00157         die("Can't allocate memory for array O, E, Esav, or R\n");
00158 
00159     // create all output records (will fail for large nrecs)
00160     DRMS_RecordSet_t *rsout = drms_create_records(drms_env, nrecs, out, DRMS_PERMANENT, &status);
00161     if (status || !rsout)
00162         die("Can't create %d records in output series %s", nrecs, out);
00163 
00164     //
00165     // MAIN LOOP
00166     //
00167     for (int n=0; n<nrecs; ++n) {
00168         warn("Processing record #%d...\n", n);
00169 
00170         DRMS_Record_t *recin = rsin->records[n];
00171         DRMS_Record_t *recout = rsout->records[n];
00172         drms_copykeys(recout, recin, 1, kDRMS_KeyClass_Explicit);
00173 
00174         // copy segment bad_pixel_list
00175         DRMS_Segment_t *seg1in = drms_segment_lookup(recin, "bad_pixel_list");
00176         if (!seg1in)
00177         {
00178             warn("No bad_pixel_list input segment found.\n");
00179         } else {
00180         DRMS_Array_t *arr1in = drms_segment_read(seg1in, DRMS_TYPE_INT, &status);
00181         if (status || !arr1in)
00182         die("Can't read bad_pixel_list\n");
00183         DRMS_Segment_t *seg1out = drms_segment_lookupnum(recout, 1);
00184         if (!seg1out)
00185         {
00186         die("Can't open bad_pixel_list output segment.\n");
00187         }
00188         if (DRMS_SUCCESS != drms_segment_write(seg1out, arr1in, 0))
00189         {
00190         die("Can't write bad_pixel_list output segment.\n");
00191         }
00192         drms_free_array(arr1in);
00193     }
00194 
00195         // read segment image_lev1
00196         DRMS_Segment_t *seg0in = drms_segment_lookup(recin, "image_lev1");
00197         if (!seg0in)
00198         {
00199             warn("No image_lev1 input segment found.\n");
00200         continue;
00201         }
00202         DRMS_Array_t *arr0in = drms_segment_read(seg0in, DRMS_TYPE_INT, &status);
00203         if (status || !arr0in)
00204             die("Can't read image_lev1\n");
00205 
00206         // pad image data array for FFT
00207         for (int j=0; j<4096; ++j)
00208             for (int i=0; i<4096; ++i) {
00209                 int tmp = ((int *)arr0in->data)[4096*j+i];
00210                 E[4098*j+i] = O[4098*j+i] = (tmp > 0) ? tmp : 0;
00211             }
00212 
00213         // do Richardson-Lucy iterations
00214         for (int k=0; k<iter; ++k) {
00215             memcpy(Esav, E, sizeof(float)*4096*4098);
00216             // B = conv(E,P)
00217             ck(DftiComputeForward(p, E));
00218             vcMul(2049*4096, E, P, B);
00219             ck(DftiComputeBackward(q, B));
00220             // R = O/B
00221             vsDiv(2*2049*4096, O, B, R);
00222             // C = corr(R,P) (or conv(R,P*))
00223             ck(DftiComputeForward(p, R));
00224             vcMulByConj(2049*4096, R, P, C);
00225             ck(DftiComputeBackward(q, C));
00226             // apply correction to get new estimate
00227             vsMul(2*2049*4096, Esav, C, E);
00228         }
00229 
00230         // copy final result, write segment
00231         int *imgout = malloc(sizeof(int)*4096*4096);
00232         if (!imgout)
00233             die("Can't allocate memory for output image\n");
00234         for (int j=0; j<4096; ++j)
00235             for (int i=0; i<4096; ++i) {
00236                 int tmp = ((int *)arr0in->data)[4096*j+i];
00237                 imgout[4096*j+i] = (tmp > 0) ? roundf(E[4098*j+i]) : tmp;
00238                 Esav[4096*j+i] = E[4098*j+i]; // for limb_fit()
00239             }
00240         DRMS_Segment_t *seg0out = drms_segment_lookupnum(recout, 0);
00241         if (!seg0out)
00242         {
00243             die("Can't open image_lev1 output segment.\n");
00244         }
00245         if (arr0in->data) free(arr0in->data);
00246         arr0in->data = imgout;
00247         if (DRMS_SUCCESS != drms_segment_write(seg0out, arr0in, 0))
00248         {
00249             die("Can't write image_lev1 output segment.\n");
00250         }
00251         drms_free_array(arr0in);
00252 
00253         // do limb fit and height formation correction
00254         int FID = drms_getkey_int(recin, "FID", &status);
00255         float CROTA2 = drms_getkey_float(recin, "CROTA2", &status);
00256         float CDELT1 = drms_getkey_float(recin, "CDELT1", &status);
00257         double OBSVR = drms_getkey_double(recin, "OBS_VR", &status);
00258         double tmpX0, tmpY0, tmpRSUN;
00259         status = limb_fit(recin, Esav, &tmpRSUN, &tmpX0, &tmpY0, 4096, 4096, 0);
00260         float RSUN_LF = tmpRSUN;
00261         float X0_LF = tmpX0;
00262         float Y0_LF = tmpY0;
00263         status = heightformation(FID,OBSVR,&CDELT1,&RSUN_LF,&X0_LF,&Y0_LF,-CROTA2);
00264         drms_setkey_float(recout, "CRPIX1", X0_LF+1);
00265         drms_setkey_float(recout, "CRPIX2", Y0_LF+1);
00266         drms_setkey_float(recout, "R_SUN", RSUN_LF);
00267         drms_setkey_float(recout, "CDELT1", CDELT1);
00268 
00269     // set CALVER32 bit
00270     long long calver32 = drms_getkey_longlong(recin, "CALVER32", &status);
00271     calver32 += 0x200000;
00272     drms_setkey_longlong(recout, "CALVER32", calver32);
00273 
00274     drms_setkey_double(recout, "DATE", CURRENT_SYSTEM_TIME);
00275     drms_setkey_int(recout, "ITER", iter);
00276     drms_setkey_string(recout, "PSF", psf);
00277     snprintf(source, sizeof source, "%s[:#%lld]", recin->seriesinfo->seriesname, recin->recnum);
00278     drms_setkey_string(recout, "SOURCE", source);
00279     }
00280 
00281     MKL_free(P);
00282     MKL_free(O);
00283     MKL_free(E);
00284     MKL_free(Esav);
00285     MKL_free(R);
00286     drms_close_records(rsin, DRMS_FREE_RECORD);
00287     drms_close_records(rsout, DRMS_INSERT_RECORD);
00288 
00289     return 0;
00290 }
proj/lev1.5_hmi/apps/lev1_dcon.c
