From 38d6cffc989678768926cb3b837d9c157ea8fc6f Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Fri, 6 May 2016 16:47:24 -0700 Subject: [PATCH 01/22] Skeleton version of spatrectif_000.c (v3) --- modules/source/spatrectif_000.c | 687 +++++++++++++++----------------- 1 file changed, 321 insertions(+), 366 deletions(-) diff --git a/modules/source/spatrectif_000.c b/modules/source/spatrectif_000.c index 3e8b2ec..e5fe5d3 100644 --- a/modules/source/spatrectif_000.c +++ b/modules/source/spatrectif_000.c @@ -1,366 +1,321 @@ -/* spatrectif_000.c */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -int errno; -#include "drp_structs.h" -#include "fitsio.h" - -float blame[DATA][numspec][MAXSLICE]; // Kernal for applying blame to lenslets for the first portion of the iterations -float weight[numspec][DATA]; // Weight normalization factor for distributing blame -float good[numspec][DATA]; // Weighted version of quality to decide if a pixel is valid. -float cmp[numspec][DATA]; // Weights for good - -float fblame[DATA][numspec][MAXSLICE]; // Kernal for applying blame during the final iterations -float fbasisv[DATA][numspec][MAXSLICE]; // influence matrix -float fweight[numspec][DATA]; // Weight factor for distributing final blame -float *fw; // Weight factor for distributing final blame -float *w; // Weight factor for distributing final blame -float t_image[DATA][numspec]; -float c_image[DATA][numspec]; -float *ti; -float *ci; -float dummy[DATA]; -float residual[DATA]; -float new_image[numspec][DATA]; -float adjust[MAXSLICE]; -float maxi; -int where; - -int spatrectif_000(int argc, void* argv[]) -{ - // Parameters input in IDL calling program - short int totalParmCount; - short int numiter; - float scale; - float relaxation; - float relax; - short int basesize; - short int (*hilo)[2]; - short int bottom[numspec]; - short int *effective; - float *bv; - float *bl, *fbl, *fbv; - float (*basis_vectors)[MAXSLICE][DATA]; - Module *pModule; - DataSet *pDataSet; - float (*Frames)[DATA]; - IDL_STRING *Headers; - float (*Noise)[DATA]; - unsigned char (*Quality)[DATA]; - float (*image)[DATA]; - float (*noise)[DATA]; - unsigned char (*quality)[DATA]; - // End of parameters input in IDL calling program - - double t1, t2, t3, t4, t5, t6; // Time counters used to examine execution time - short int i=0, ii=0, j=0, jj=0, sp=0, l=0; - long int in1=0, in2=0, in3=0, index[numspec]; - - // Make a temporary residual matrix - long naxes[3]; - naxes[0] = DATA; - naxes[1] = DATA; - naxes[2] = (21); - // ******************************* - - // These parameters should be set in the same order as theay are passed - // from the IDL code. This is not yet automated, and I'm not sure how to - // do it. - i = 0; - totalParmCount = *(short int * )argv[i++]; - numiter = *(short int * )argv[i++]; - relaxation = *(float * )argv[i++]; - basesize = *(short int * )argv[i++]; - hilo = (short int (*)[2] )argv[i++]; - effective = (short int * )argv[i++]; - basis_vectors = (float (*)[MAXSLICE][DATA])argv[i++]; - pModule = (Module * )argv[i++]; - pDataSet = (DataSet * )argv[i++]; - Frames = (float (*)[DATA] )argv[i++]; - Headers = (IDL_STRING * )argv[i++]; - Noise = (float (*)[DATA] )argv[i++]; - Quality = (unsigned char (*)[DATA] )argv[i++]; - image = (float (*)[DATA] )argv[i++]; - noise = (float (*)[DATA] )argv[i++]; - quality = (unsigned char (*)[DATA] )argv[i++]; - scale = *(float * )argv[i++]; // Plate scale - - /* - * Start placing items from the original rectification code here. - * This code will rectify an input data frame. - */ - // printf("Image is spatrectif is %x.\n",image); - printf( "spatrectif_000.c: Now processing RAW data...\n"); - printf("Number of iterations = %d.\n",numiter); - /* printf("Relaxation Parameter = %f.\n",relaxation); - printf("Platescale = %f.\n",scale); */ - (void)fflush(stdout); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Iterate on each lenslet. - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - memset((void *)image,0, IMAGECNT ); // initial guess of the solution - - t1=systime(); - - // Set the local pointer bv equal to the address of the lowest member of the basis_vector - bv = basis_vectors[0][0]; - for (sp = 0; sp< numspec; sp++) - { - bottom[sp]=hilo[sp][0]; - } - - printf("Calculating weights.\n"); - (void)fflush(stdout); - - for (sp=0; sp maxi ) - // { - // maxi = fbv[l]; - // where = l; - // } - // } - // bl[where]= maxi; // Alternative initial blame where only peak pixel is used. - // w[i] = bl[where]; - for (l=0; l 0.0 ) - { - for (l=0; l 0.0 ) - { - for (l=0; l 14) { - // After first set of iterations, settle down to stable solution - ci[sp]+=relax*residual[j]* fbl[jj]; - } - // if ( ii > 20) { - // After first set of iterations, settle down to stable solution - // ci[sp]+=2.0*relax*residual[j]* fbl[jj]; - //} - j++; - } - } - // For the first iterations, correction is sharpened and applied. - if ( ii < 10 ) { - for (sp=0; sp< numspec; sp++) { // Correct middle rows - if ( weight[sp][i] > 0.0 ) {// Increment guess image - ci[sp]+= ti[sp]; - if ( sp > 63 ) - if ( weight[sp-64][i] > 0.0 ) // Increment guess image - ci[sp]+= 0.4*(ti[sp]-ti[sp-64]); - if ( sp < numspec-64 ) - if ( weight[sp+64][i] > 0.0 ) // Increment guess image - ci[sp]+= 0.4*(ti[sp]-ti[sp+64]); - } - } - } - if ( (ii > 9) && (ii < 15) ) { - for (sp=0; sp< numspec; sp++) { // Correct middle rows - if ( weight[sp][i] > 0.0 ) {// Increment guess image - ci[sp]+= ti[sp]; - if ( sp > 63 ) - if ( weight[sp-64][i] > 0.0 ) // Increment guess image - ci[sp]-= 0.2*(ti[sp]-ti[sp-64]); - if ( sp < numspec-64 ) - if ( weight[sp+64][i] > 0.0 ) // Increment guess image - ci[sp]-= 0.2*(ti[sp]-ti[sp+64]); - } - } - } - - // if ( (scale > 0.099) && (ii < 2) ) { - - } // end of loop over image. - // Smooth early iterations along the spectral direction. - if ( (ii > 11) && (ii < 18) ) - { - for (sp = 0; sp 4.5 ) quality[sp][i]=9; // Implies over have of the weights have good pixels. - } // for each spectral channel sp ... - } - printf("\n"); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - t2 = systime(); - // (void)printf("Total Time = %lf\n", t2-t1 ); - - // writefitsimagefile("!/sdata1101/osiris-data/osiris2/050224/SPEC/ORP/resid.fits", FLOAT_IMG, 3, naxes, resid); - // free(resid); - - -#ifdef SAVE_INTERMEDIATE_FILES - //naxes[0] = DATA; - //naxes[1] = numspec; - //writefitsimagefile("testimage.fits", FLOAT_IMG, 2, naxes, image); -#endif - // printf("Image is spatrectif is %x.\n",image); - - return 0; -} +/* spatrectif_000.c */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +int errno; +#include "drp_structs.h" +#include "fitsio.h" + +float blame[DATA][numspec][MAXSLICE]; // Kernal for applying blame to lenslets for the first portion of the iterations +float weight[numspec][DATA]; // Weight normalization factor for distributing blame +float good[numspec][DATA]; // Weighted version of quality to decide if a pixel is valid. +float cmp[numspec][DATA]; // Weights for good +float fblame[DATA][numspec][MAXSLICE]; // Kernal for applying blame during the final iterations +float fbasisv[DATA][numspec][MAXSLICE]; // influence matrix (basis vector) +float fweight[numspec][DATA]; // Weight factor for distributing final blame +float *fw; // Weight factor for distributing final blame +float *w; // Weight factor for distributing final blame +float t_image[DATA][numspec]; // Temporary extracted spectra +float c_image[DATA][numspec]; // Cumulative extracted spectra +float *ti; // For each column accumulate delta to lenslet value +float *ci; // For each column best current estimate of lenslet value +float t_fbv[MAXSLICE]; // Temporary basis vector to manipulate if necessary +float dummy[DATA]; // Model of the flux given current lenslet estimate +float residual[DATA]; // Residual between data - model + +int where; + +int spatrectif_000(int argc, void* argv[]) +{ + // Parameters input in IDL calling program + short int totalParmCount; + short int numiter; + float scale; + float relaxation; + float relax; + short int basesize; + short int (*hilo)[2]; + short int bottom[numspec]; + short int *effective; + float *bv; + float *bl, *fbl, *fbv; + float (*basis_vectors)[MAXSLICE][DATA]; + Module *pModule; + DataSet *pDataSet; + float (*Frames)[DATA]; + IDL_STRING *Headers; + float (*Noise)[DATA]; + unsigned char (*Quality)[DATA]; + float (*image)[DATA]; + float (*noise)[DATA]; + unsigned char (*quality)[DATA]; + // End of parameters input in IDL calling program + + double t1, t2, t3, t4, t5, t6; // Time counters used to examine execution time + short int i=0, ii=0, j=0, jj=0, sp=0, l=0; + long int in1=0, in2=0, in3=0, index[numspec]; + + // Make a temporary residual matrix + long naxes[3]; + naxes[0] = DATA; + naxes[1] = DATA; + naxes[2] = (21); + // ******************************* + + // These parameters should be set in the same order as theay are passed + // from the IDL code. This is not yet automated, and I'm not sure how to + // do it. + i = 0; + totalParmCount = *(short int * )argv[i++]; + numiter = *(short int * )argv[i++]; + relaxation = *(float * )argv[i++]; + basesize = *(short int * )argv[i++]; + hilo = (short int (*)[2] )argv[i++]; + effective = (short int * )argv[i++]; + basis_vectors = (float (*)[MAXSLICE][DATA])argv[i++]; + pModule = (Module * )argv[i++]; + pDataSet = (DataSet * )argv[i++]; + Frames = (float (*)[DATA] )argv[i++]; + Headers = (IDL_STRING * )argv[i++]; + Noise = (float (*)[DATA] )argv[i++]; + Quality = (unsigned char (*)[DATA] )argv[i++]; + image = (float (*)[DATA] )argv[i++]; + noise = (float (*)[DATA] )argv[i++]; + quality = (unsigned char (*)[DATA] )argv[i++]; + scale = *(float * )argv[i++]; // Plate scale + + /* + * Start placing items from the original rectification code here. + * This code will rectify an input data frame. + */ + // printf("Image is spatrectif is %x.\n",image); + printf( "spatrectif_000.c: Now processing RAW data...\n"); + printf("Number of iterations = %d.\n",numiter); + /* printf("Relaxation Parameter = %f.\n",relaxation); + printf("Platescale = %f.\n",scale); */ + (void)fflush(stdout); + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Iterate on each lenslet. + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + memset((void *)image,0, IMAGECNT ); // initial guess of the solution + + t1=systime(); + + // Set the local pointer bv equal to the address of the lowest member of the basis_vector + bv = basis_vectors[0][0]; + for (sp = 0; sp< numspec; sp++) + { + bottom[sp]=hilo[sp][0]; + } + + printf("Calculating weights.\n"); + (void)fflush(stdout); + + // Making local copies of the basis vector + for (sp=0; sp 0.0 ) + { + for (l=0; l 0.0 ) + { + for (l=0; l 14) { + // After first set of iterations, settle down to stable solution + ti[sp]+=relax*residual[j]* fbl[jj]; + //ti[sp]+=relax*residual[j]* bl[jj]; + } + j++; // march up column to match PSF location + } + } + // Correction is applied. + for (sp=0; sp< numspec; sp++) { + if ( weight[sp][i] > 0.0 ) { + ci[sp]+= ti[sp]; + } + } + + } // end of loop over image. + } // for each iteration + + for (i=0; i 4.5 ) quality[sp][i]=9; // Implies over have of the weights have good pixels. + } // for each spectral channel sp ... + } + printf("\n"); + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + t2 = systime(); + // (void)printf("Total Time = %lf\n", t2-t1 ); + // writefitsimagefile("!/sdata1101/osiris-data/osiris2/050224/SPEC/ORP/resid.fits", FLOAT_IMG, 3, naxes, resid); + // free(resid); + + +#ifdef SAVE_INTERMEDIATE_FILES + //naxes[0] = DATA; + //naxes[1] = numspec; + //writefitsimagefile("testimage.fits", FLOAT_IMG, 2, naxes, image); +#endif + // printf("Image is spatrectif is %x.\n",image); + + return 0; +} From 049e5ed3d37bb4a593322c26cca12982a0f98317 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 12 May 2016 13:56:35 -0700 Subject: [PATCH 02/22] Changed default number of iterations for ARP_SEC from 4 to 40, in the RPBconfig.xml file. --- backbone/SupportFiles/RPBconfig.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/backbone/SupportFiles/RPBconfig.xml b/backbone/SupportFiles/RPBconfig.xml index 85ab110..19794d1 100644 --- a/backbone/SupportFiles/RPBconfig.xml +++ b/backbone/SupportFiles/RPBconfig.xml @@ -124,7 +124,7 @@ scaledskysub_COMMON___show_plots = "Yes" spatrectif_COMMON___numiter="4" - spatrectif_ARP_SPEC_numiter="4" + spatrectif_ARP_SPEC_numiter="40" spatrectif_ORP_SPEC_numiter="25" spatrectif_COMMON___relaxation="1.0" From a0b18c0db8a1fc3d6cdeef85256a6e888f77475e Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 12 May 2016 14:21:22 -0700 Subject: [PATCH 03/22] Revert "Changed default number of iterations for ARP_SEC from 4 to 40, in the RPBconfig.xml file." This reverts commit 049e5ed3d37bb4a593322c26cca12982a0f98317. Reverting this commit so I can add a new branch and make the change there. This allows me to make a new pull request separate from the current pull request for this branch (the spatrectif skeleton version). --- backbone/SupportFiles/RPBconfig.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/backbone/SupportFiles/RPBconfig.xml b/backbone/SupportFiles/RPBconfig.xml index 19794d1..85ab110 100644 --- a/backbone/SupportFiles/RPBconfig.xml +++ b/backbone/SupportFiles/RPBconfig.xml @@ -124,7 +124,7 @@ scaledskysub_COMMON___show_plots = "Yes" spatrectif_COMMON___numiter="4" - spatrectif_ARP_SPEC_numiter="40" + spatrectif_ARP_SPEC_numiter="4" spatrectif_ORP_SPEC_numiter="25" spatrectif_COMMON___relaxation="1.0" From 26cecc863fd7d1d4b0131873895f3e004c24f4f2 Mon Sep 17 00:00:00 2001 From: Jim Lyke Date: Wed, 7 Sep 2016 08:57:12 -1000 Subject: [PATCH 04/22] How I installed the DRP --- EXAMPLE.md | 1 + 1 file changed, 1 insertion(+) diff --git a/EXAMPLE.md b/EXAMPLE.md index 10bc92f..bdfb873 100644 --- a/EXAMPLE.md +++ b/EXAMPLE.md @@ -284,3 +284,4 @@ tests/test_emission_line/test_emission_line.py . ``` ## Test Successful + From 01794fe8a808a16131ce3c97508041251051f90e Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Wed, 7 Sep 2016 14:31:53 -0700 Subject: [PATCH 05/22] Test #1: commented out lines 312 - 324 thereby removing cuts that are made on rect mat pixels. --- modules/source/mkrecmatrx_000.c | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/modules/source/mkrecmatrx_000.c b/modules/source/mkrecmatrx_000.c index 7cf7d7e..055dc1b 100644 --- a/modules/source/mkrecmatrx_000.c +++ b/modules/source/mkrecmatrx_000.c @@ -309,19 +309,19 @@ int mkrecmatrx_000(int argc, void* argv[]) { if ( basis_vectors[sp][j][i] > 1.0 ) basis_vectors[sp][j][i]=0.0; pretotal += basis_vectors[sp][j][i]; - if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])*2.0 ) ) { - basis_vectors[sp][j][i]= (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])/2.0; - } - if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])*2.0 ) ) { - basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; - } - if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/2.0 ) ) { - if ( (basis_vectors[sp][j][i-1] > 0.1) ) { - if ( (basis_vectors[sp][j][i+1] > 0.1) ) { - basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; - } - } - } + //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])*2.0 ) ) { + // basis_vectors[sp][j][i]= (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])/2.0; + //} + //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])*2.0 ) ) { + // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; + //} + //if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/2.0 ) ) { + // if ( (basis_vectors[sp][j][i-1] > 0.1) ) { + // if ( (basis_vectors[sp][j][i+1] > 0.1) ) { + // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; + // } + // } + //} } // Look at the upper and lower edges of each strip. Set them to 0, if they are too large if ( fabs(basis_vectors[sp][0][i]) > (pretotal/2.0) ) basis_vectors[sp][0][i]=0.0; From 745f29e92581ce8cf328ac05cb007202b7785951 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Wed, 7 Sep 2016 14:42:09 -0700 Subject: [PATCH 06/22] Kept the first two cuts commented out (lines 312 - 317) and changed 3rd cut to do bad things! --- modules/source/mkrecmatrx_000.c | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/modules/source/mkrecmatrx_000.c b/modules/source/mkrecmatrx_000.c index 055dc1b..4ad672f 100644 --- a/modules/source/mkrecmatrx_000.c +++ b/modules/source/mkrecmatrx_000.c @@ -315,13 +315,13 @@ int mkrecmatrx_000(int argc, void* argv[]) { //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])*2.0 ) ) { // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; //} - //if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/2.0 ) ) { - // if ( (basis_vectors[sp][j][i-1] > 0.1) ) { - // if ( (basis_vectors[sp][j][i+1] > 0.1) ) { - // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; - // } - // } - //} + if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/0.2 ) ) { + if ( (basis_vectors[sp][j][i-1] > 0.1) ) { + if ( (basis_vectors[sp][j][i+1] > 0.1) ) { + basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; + } + } + } } // Look at the upper and lower edges of each strip. Set them to 0, if they are too large if ( fabs(basis_vectors[sp][0][i]) > (pretotal/2.0) ) basis_vectors[sp][0][i]=0.0; From 3e0493945d63b101fc067f58b3bc23055ff1ab94 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Wed, 7 Sep 2016 16:08:40 -0700 Subject: [PATCH 07/22] Lines 312 - 324 commented out again. --- modules/source/mkrecmatrx_000.c | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/modules/source/mkrecmatrx_000.c b/modules/source/mkrecmatrx_000.c index 4ad672f..a090498 100644 --- a/modules/source/mkrecmatrx_000.c +++ b/modules/source/mkrecmatrx_000.c @@ -309,20 +309,22 @@ int mkrecmatrx_000(int argc, void* argv[]) { if ( basis_vectors[sp][j][i] > 1.0 ) basis_vectors[sp][j][i]=0.0; pretotal += basis_vectors[sp][j][i]; + (void)printf("why"); //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])*2.0 ) ) { // basis_vectors[sp][j][i]= (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])/2.0; //} //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])*2.0 ) ) { // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; //} - if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/0.2 ) ) { - if ( (basis_vectors[sp][j][i-1] > 0.1) ) { - if ( (basis_vectors[sp][j][i+1] > 0.1) ) { - basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; - } - } - } - } + //if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/2.0 ) ) { + // if ( (basis_vectors[sp][j][i-1] > 0.1) ) { + // if ( (basis_vectors[sp][j][i+1] > 0.1) ) { + // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; + // } + // } + //} + + } // Look at the upper and lower edges of each strip. Set them to 0, if they are too large if ( fabs(basis_vectors[sp][0][i]) > (pretotal/2.0) ) basis_vectors[sp][0][i]=0.0; if ( fabs(basis_vectors[sp][(basesize-1)][i]) > (pretotal/2.0) ) basis_vectors[sp][(basesize-1)][i]=0.0; From 82f7ae39e574b32a5dfc7a73311536bf26c9cf88 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 8 Sep 2016 10:26:22 -0700 Subject: [PATCH 08/22] Adding an automated rectmat making code that is useful for quickly making and organizing matrices with different values of weightlimit, slice, and maxslice values. --- tests/test_mkrecmatrx/make_rectmat.py | 146 ++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 tests/test_mkrecmatrx/make_rectmat.py diff --git a/tests/test_mkrecmatrx/make_rectmat.py b/tests/test_mkrecmatrx/make_rectmat.py new file mode 100644 index 0000000..57612b4 --- /dev/null +++ b/tests/test_mkrecmatrx/make_rectmat.py @@ -0,0 +1,146 @@ +import os +import pyfits + +source_dir = '/Users/aboehle/osiris_dev/OsirisDRP/modules/source/' #'/Users/aboehle/osiris/drs/modules/source/' +RPB_dir = '/Users/aboehle/osiris_dev/OsirisDRP/backbone/SupportFiles/' #'/Users/aboehle/osiris/drs/backbone/SupportFiles/' + +rectmat_outdir = 'rectmat/' #rectmats_QSOdata' +reduce_dir = '/Users/aboehle/research/irlab/projects/osiris_rectmat_tests/images/160412/reduce/'#160317/reduce/' + +# make a rect mat with whatever slice values you want, update the dir_suffix + +def make_rectmat(slice,maxslice,weightlimit,dir_suffix=''): + + xml_filename = 'Kbb_50.xml' #'Hn3_100.xml' + xml_test_filename = 'testrectmat_Kbb_50_newpipeline.xml' #'testrectmat_Hn3_100_newpipeline.xml' + + #dir_suffix = 'test123' # for additional tests + + if dir_suffix: + dir_suffix = '_'+dir_suffix + + # Check that output directories exist + if not os.path.exists(reduce_dir + rectmat_outdir + '/weightlimit%1.2f_slice%i_maxslice%i%s/' % (weightlimit, slice, maxslice,dir_suffix)): + os.mkdir(reduce_dir + rectmat_outdir + '/weightlimit%1.2f_slice%i_maxslice%i%s/' % (weightlimit, slice, maxslice,dir_suffix)) + if not os.path.exists(reduce_dir + '/cubed_frames_newpipeline/weightlimit%1.2f_slice%i_maxslice%i%s/' % (weightlimit, slice, maxslice,dir_suffix)): + os.mkdir(reduce_dir + '/cubed_frames_newpipeline/weightlimit%1.2f_slice%i_maxslice%i%s/' % (weightlimit, slice, maxslice,dir_suffix)) + + setup_drp_files(slice,maxslice,weightlimit) + setup_xml_files(reduce_dir, xml_filename, xml_test_filename, rectmat_outdir, slice, maxslice, weightlimit ,dir_suffix) + recompile_ccode() + + raw_input("\n\nClose the OSIRIS DRP if it is currently running!\nThen press enter to drop the new xml files. ") + + os.chdir(reduce_dir + '/xml/') + cmds = ['run_odrp &','osirisDropDRF '+xml_filename + ' 1','osirisDropDRF '+xml_test_filename + ' 2'] + for cmd in cmds: + print cmd + os.system(cmd) + +def check_rectmats(rectmat_name = 's160412_c006___infl_Kbb_050.fits'): #'s160318_c003___infl_Hn3_100.fits'): + + slice_arr = [[14,16],[20,22],[28,30]] + weight_lim_arr = [0,0.01] + + for slice,maxslice in slice_arr: + for lmt in weight_lim_arr: + filename = reduce_dir +rectmat_outdir + '/weightlimit%1.2f_slice%i_maxslice%i/' % (lmt, slice, maxslice) + rectmat_name + print filename + mat_hdu = pyfits.open(filename) + hdr = mat_hdu[0].header + + if slice == hdr['slice']: + print 'slice test: passed' + slicefailed = False + else: + print 'slice test: FAILED' + slicefailed = True + + if lmt == hdr['wtlimit']: + print 'weightlimit test: passed' + lmtfailed = False + else: + print 'weightlimit test: FAILED' + lmtfailed = True + + if maxslice == mat_hdu[2].data.shape[1]: + print 'maxslice test: passed' + maxslice_failed = False + else: + print 'maxslice test: FAILED' + maxslice_failed = True + + if maxslice_failed or slicefailed or lmtfailed: + print 'ERROR with rectmat: '+filename + else: + print 'rectmat passed!' + + print + + #print '\nInput value/Rectmat Value:' + #print slice,hdr['slice'] + #print lmt,hdr['wtlimit'] + #print maxslice,mat_hdu[2].data.shape[1] + #print + + +def setup_drp_files(slice,maxslice,weightlimit): + + rpb_file = open(RPB_dir+'RPBconfig.xml','r') + drp_file = open(source_dir+'drp_structs.h','r') + + # Edit RPB_file lines + rpb_lines = rpb_file.readlines() + rpb_lines[74] = ' mkrecmatrx_COMMON___slice="%i"\n' % (slice) + rpb_lines[76] = ' mkrecmatrx_COMMON___weight_limit="%1.2f"\n' % (weightlimit) + rpb_file.close() + + rpb_file = open(RPB_dir+'RPBconfig.xml','w') + rpb_file.writelines(rpb_lines) + rpb_file.close() + + # Edit drp file lines + drp_lines = drp_file.readlines() + drp_lines[21] = '#define MAXSLICE %i // Maximum slice of image in pixels; original value = 16\n' % (maxslice) + drp_file.close() + + drp_file = open(source_dir+'drp_structs.h','w') + drp_file.writelines(drp_lines) + drp_file.close() + +def setup_xml_files(reduce_dir, xml_filename, xml_test_filename, rectmat_outdir, slice, maxslice, weightlimit, dir_suffix): + + xml_file = open(reduce_dir + '/xml/'+xml_filename,'r') + xml_lines = xml_file.readlines() + xml_lines[8] = ' OutputDir="'+reduce_dir + '/' +rectmat_outdir+'/weightlimit%1.2f_slice%i_maxslice%i%s/">\n' % (weightlimit, slice, maxslice, dir_suffix) + xml_lines[31] = ' OutputDir="/Users/aboehle/research/IRLab/Projects/osiris_rectmat_tests/images/160412/reduce/'+rectmat_outdir+'/weightlimit%1.2f_slice%i_maxslice%i%s/"\n' % (weightlimit, slice, maxslice, dir_suffix) + xml_file.close() + + xml_file = open(reduce_dir + '/xml/'+xml_filename,'w') # xml file to make rect mat + xml_file.writelines(xml_lines) + xml_file.close() + + xml_test_file = open(reduce_dir + '/xml/'+xml_test_filename,'r') # xml file to reduce single column scan + xml_test_lines = xml_test_file.readlines() + xml_test_lines[3] = '\n' % (weightlimit, slice, maxslice, dir_suffix) + xml_test_lines[8] = '\n' % (weightlimit, slice, maxslice, dir_suffix) + xml_test_file.close() + + xml_test_file = open(reduce_dir + '/xml/'+xml_test_filename,'w') + xml_test_file.writelines(xml_test_lines) + xml_test_file.close() + + +def recompile_ccode(): + + current_dir = os.getcwd() + + os.chdir(source_dir) + cmds = ['rm *.o','rm libosiris_drp_ext_null.so.0.0','make -f local_Makefile'] + + for cmd in cmds: + print cmd + os.system(cmd) + + os.chdir(current_dir) + From 775fcf7ef2d22fd061bc67fba2703e0a4ad40a74 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 8 Sep 2016 10:27:26 -0700 Subject: [PATCH 09/22] Adding an automated rectmat making code that is useful for quickly making and organizing matrices with different values of weightlimit, slice, and maxslice values. --- tests/test_mkrecmatrx/make_rectmat.py | 32 ++++++++++++++------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/tests/test_mkrecmatrx/make_rectmat.py b/tests/test_mkrecmatrx/make_rectmat.py index 57612b4..b3aca7d 100644 --- a/tests/test_mkrecmatrx/make_rectmat.py +++ b/tests/test_mkrecmatrx/make_rectmat.py @@ -1,22 +1,27 @@ import os import pyfits -source_dir = '/Users/aboehle/osiris_dev/OsirisDRP/modules/source/' #'/Users/aboehle/osiris/drs/modules/source/' -RPB_dir = '/Users/aboehle/osiris_dev/OsirisDRP/backbone/SupportFiles/' #'/Users/aboehle/osiris/drs/backbone/SupportFiles/' +# location of OSIRIS DRP: +drp_dir = '/Users/aboehle/osiris_dev/OsirisDRP/' -rectmat_outdir = 'rectmat/' #rectmats_QSOdata' -reduce_dir = '/Users/aboehle/research/irlab/projects/osiris_rectmat_tests/images/160412/reduce/'#160317/reduce/' +source_dir = drp_dir + '/modules/source/' +RPB_dir = drp_dir + '/backbone/SupportFiles/' -# make a rect mat with whatever slice values you want, update the dir_suffix +# location of rect mat raw scans and created rect mats +reduce_dir = '/Users/aboehle/research/irlab/projects/osiris_rectmat_tests/images/160412/reduce/' # from here, raw dir is: ../SPEC/raw/ +rectmat_outdir = 'rectmat/' # sub-directory of reduce_dir above -def make_rectmat(slice,maxslice,weightlimit,dir_suffix=''): - xml_filename = 'Kbb_50.xml' #'Hn3_100.xml' - xml_test_filename = 'testrectmat_Kbb_50_newpipeline.xml' #'testrectmat_Hn3_100_newpipeline.xml' +def make_rectmat(slice,maxslice,weightlimit,dir_suffix='', xml_filename = 'Kbb_050.xml', xml_test_filename = 'testrectmat_Kbb_50_newpipeline.xml'): + ''' + Make a rect mat with a given slice, maxslice, and weight limit value. Save the rect mat into a sub-directory of the reduce_dir (reduce_dir is hard-coded above). This sub-directory is named with the weightlimit, slice, and maxslice values followed by an optional dir_suffix. - #dir_suffix = 'test123' # for additional tests + To make rect mat: updates the given rect mat xml file (named "xml_filename") with new output dir and the appropriate drp files with new maxslice, slice, and weight limit values. Prompts the user to close the OSIRIS DRP so it can be reopened with the new drp parameter values. Recompiles the c code in drp modules/source/ directory. - if dir_suffix: + Also updates the given test filename (named "xml_test_filename") with the new data directory + ''' + + if dir_suffix and dir_suffix[0] != '_': dir_suffix = '_'+dir_suffix # Check that output directories exist @@ -37,11 +42,8 @@ def make_rectmat(slice,maxslice,weightlimit,dir_suffix=''): print cmd os.system(cmd) -def check_rectmats(rectmat_name = 's160412_c006___infl_Kbb_050.fits'): #'s160318_c003___infl_Hn3_100.fits'): - - slice_arr = [[14,16],[20,22],[28,30]] - weight_lim_arr = [0,0.01] - +def check_rectmats(rectmat_name = 's160412_c006___infl_Kbb_050.fits', slice_arr = [[14,16],[20,22],[28,30]], weight_lim_arr = [0,0.01]): #'s160318_c003___infl_Hn3_100.fits'): + for slice,maxslice in slice_arr: for lmt in weight_lim_arr: filename = reduce_dir +rectmat_outdir + '/weightlimit%1.2f_slice%i_maxslice%i/' % (lmt, slice, maxslice) + rectmat_name From 9f11afdaeb6ff4001cf2bbb2da79d4301618694e Mon Sep 17 00:00:00 2001 From: Jessica Lu Date: Thu, 8 Sep 2016 10:29:08 -0700 Subject: [PATCH 10/22] Added comments --- modules/spatrectif_000.pro | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/modules/spatrectif_000.pro b/modules/spatrectif_000.pro index d7f86a8..ede685a 100644 --- a/modules/spatrectif_000.pro +++ b/modules/spatrectif_000.pro @@ -28,12 +28,15 @@ CASE BranchID OF RETURN, ERR_BADCASE END ; CASE BadType ENDCASE - ; Now perform the common code execution for spatial rectification - ; Read in the Calibration file that will contain: - ; BASESIZE Keyword - ; hilo short int array[NUMSPEC][2] - ; effective short int array[NUMSPEC] - ; basis_vectors float array[NUMSPEC][MAXSLICE][DATAFRAMEDIM] + ; Now perform the common code execution for spatial rectification + ; Read in the Calibration file that will contain: + ; BASESIZE Keyword + ; hilo short int array[NUMSPEC][2] -- The row index of the bottom (low) and the top (high) raw for each spectral slice. + ; effective short int array[NUMSPEC] -- A flag to indicate whether this spaxel is sufficiently illuminated (0=bad, 1=good) + ; basis_vectors float array[NUMSPEC][MAXSLICE][DATAFRAMEDIM] -- The actual influence matrix. + ; Note NUMSPEC = typically 1216 spectra + ; Note MAXSLICE = 16 (default) but up to 32 (height of slice, orthogonal to dispersion direction) + ; Note DATAFRAMEDIM = 2048 (width of slice in dispersion direction) thisModuleIndex = drpModuleIndexFromCallSequence(Modules, functionName) FileName = drpXlateFileName(Modules[thisModuleIndex].CalibrationFile) pHilo = PTR_NEW(/ALLOCATE_HEAP) From 6337ad1b0f7dbd47363530e6a60ee6516e574e61 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 8 Sep 2016 10:35:00 -0700 Subject: [PATCH 11/22] Changed drp_structs and mkrectmatrx back to original form --- modules/source/drp_structs.h | 2 +- modules/source/mkrecmatrx_000.c | 27 +++++++++++++-------------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/modules/source/drp_structs.h b/modules/source/drp_structs.h index 2468a82..f9d430a 100644 --- a/modules/source/drp_structs.h +++ b/modules/source/drp_structs.h @@ -19,7 +19,7 @@ #define MAXFRAMESINDATASETS 30 #define DATA 2048 // Equivalent to BIGFRAMES' FRAMESIZE -#define MAXSLICE 16 // Maximum slice of image in pixels +#define MAXSLICE 16 // Maximum slice of image in pixels; original value = 16 #define numspec 1216 // numspec and numcolumn are same for BB and NB. #define numcolumn 19 // ... final NB frames will have 51 columns and (51*66)=3366 spectra! #define specpercol (numspec/numcolumn) // # of spectra per lenslet column diff --git a/modules/source/mkrecmatrx_000.c b/modules/source/mkrecmatrx_000.c index a090498..e1ddefe 100644 --- a/modules/source/mkrecmatrx_000.c +++ b/modules/source/mkrecmatrx_000.c @@ -309,20 +309,19 @@ int mkrecmatrx_000(int argc, void* argv[]) { if ( basis_vectors[sp][j][i] > 1.0 ) basis_vectors[sp][j][i]=0.0; pretotal += basis_vectors[sp][j][i]; - (void)printf("why"); - //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])*2.0 ) ) { - // basis_vectors[sp][j][i]= (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])/2.0; - //} - //if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])*2.0 ) ) { - // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; - //} - //if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/2.0 ) ) { - // if ( (basis_vectors[sp][j][i-1] > 0.1) ) { - // if ( (basis_vectors[sp][j][i+1] > 0.1) ) { - // basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; - // } - // } - //} + if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])*2.0 ) ) { + basis_vectors[sp][j][i]= (basis_vectors[sp][j-1][i]+basis_vectors[sp][j+1][i])/2.0; + } + if ( fabs(basis_vectors[sp][j][i]) > ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])*2.0 ) ) { + basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; + } + if ( fabs(basis_vectors[sp][j][i]) < ( (basis_vectors[sp][j][i-1]+basis_vectors[sp][j][i+1])/2.0 ) ) { + if ( (basis_vectors[sp][j][i-1] > 0.1) ) { + if ( (basis_vectors[sp][j][i+1] > 0.1) ) { + basis_vectors[sp][j][i]= (basis_vectors[sp][j][i+1]+basis_vectors[sp][j][i-1])/2.0; + } + } + } } // Look at the upper and lower edges of each strip. Set them to 0, if they are too large From 51e62b15d5c1b4fbb13d5c1b39c6afca9d672fb9 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 8 Sep 2016 10:35:44 -0700 Subject: [PATCH 12/22] Updated CFITSIOLIBDIR and IDL_INCLUDE with locations on my local machine. --- modules/source/local_Makefile | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/modules/source/local_Makefile b/modules/source/local_Makefile index 8eb925b..f340a22 100755 --- a/modules/source/local_Makefile +++ b/modules/source/local_Makefile @@ -12,8 +12,10 @@ #IDL_INCLUDE = /usr/local/idl/7.1/idl71/external/include #IDL_INCLUDE = /Applications/itt/idl70/external/include -IDL_INCLUDE = /Applications/exelis/idl/external/include -CFITSIOLIBDIR = /opt/local/lib/ +#IDL_INCLUDE = /Applications/exelis/idl/external/include +IDL_INCLUDE = /usr/local/itt/idl/external/include +#CFITSIOLIBDIR = /opt/local/lib/ +CFITSIOLIBDIR = /Applications/cfitsio/lib OUTPUTDIR = . # Include files. From 8975fac4d6b1dafa6e6f3855063c317434431e40 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 8 Sep 2016 11:03:57 -0700 Subject: [PATCH 13/22] Revert to original version of modules/source/local_Makefile --- modules/source/local_Makefile | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/modules/source/local_Makefile b/modules/source/local_Makefile index f340a22..8eb925b 100755 --- a/modules/source/local_Makefile +++ b/modules/source/local_Makefile @@ -12,10 +12,8 @@ #IDL_INCLUDE = /usr/local/idl/7.1/idl71/external/include #IDL_INCLUDE = /Applications/itt/idl70/external/include -#IDL_INCLUDE = /Applications/exelis/idl/external/include -IDL_INCLUDE = /usr/local/itt/idl/external/include -#CFITSIOLIBDIR = /opt/local/lib/ -CFITSIOLIBDIR = /Applications/cfitsio/lib +IDL_INCLUDE = /Applications/exelis/idl/external/include +CFITSIOLIBDIR = /opt/local/lib/ OUTPUTDIR = . # Include files. From 1a72706e2c3929938f62e011990491bd99f11cf0 Mon Sep 17 00:00:00 2001 From: Jessica Lu Date: Thu, 8 Sep 2016 11:48:13 -0700 Subject: [PATCH 14/22] Moved wrapper file to test directory --- .../test_misflux_singlearc}/misflux_singlearc_wrapper.pro | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {testing_scripts/flux_misassign => tests/test_misflux_singlearc}/misflux_singlearc_wrapper.pro (100%) diff --git a/testing_scripts/flux_misassign/misflux_singlearc_wrapper.pro b/tests/test_misflux_singlearc/misflux_singlearc_wrapper.pro similarity index 100% rename from testing_scripts/flux_misassign/misflux_singlearc_wrapper.pro rename to tests/test_misflux_singlearc/misflux_singlearc_wrapper.pro From 285bfda04476e7666ea1a4e5a29c97b15f41fecd Mon Sep 17 00:00:00 2001 From: Jessica Lu Date: Thu, 8 Sep 2016 11:52:57 -0700 Subject: [PATCH 15/22] Added DRF and changed file location for kelly's new test. --- .../000.Kbb_050_1col_201603.drf | 17 +++++++++++++++++ .../misflux_singlearc_wrapper.pro | 2 +- 2 files changed, 18 insertions(+), 1 deletion(-) create mode 100755 tests/test_misflux_singlearc/000.Kbb_050_1col_201603.drf diff --git a/tests/test_misflux_singlearc/000.Kbb_050_1col_201603.drf b/tests/test_misflux_singlearc/000.Kbb_050_1col_201603.drf new file mode 100755 index 0000000..6deccc5 --- /dev/null +++ b/tests/test_misflux_singlearc/000.Kbb_050_1col_201603.drf @@ -0,0 +1,17 @@ + + + + + + + + + + + + + + + + + diff --git a/tests/test_misflux_singlearc/misflux_singlearc_wrapper.pro b/tests/test_misflux_singlearc/misflux_singlearc_wrapper.pro index 18f62d2..7047a75 100644 --- a/tests/test_misflux_singlearc/misflux_singlearc_wrapper.pro +++ b/tests/test_misflux_singlearc/misflux_singlearc_wrapper.pro @@ -1,6 +1,6 @@ pro misflux_singlearc_wrapper -infile = getenv('OSIRIS_ROOT')+'/tests/test_misflux_singlearc/s160318_c017026_Kbb_050.fits' +infile = getenv('OSIRIS_ROOT')+'/tests/test_misflux_singlearc/data_misflux_single/s160318_c017026_Kbb_050.fits' refchannel = 904 From 35e38a9d366a9d3796fa3a4f0ad2434d74e31675 Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Thu, 8 Sep 2016 11:57:53 -0700 Subject: [PATCH 16/22] changed spatrectif_000.c version to the original --- modules/source/spatrectif_000.c | 687 +++++++++++++++++--------------- 1 file changed, 366 insertions(+), 321 deletions(-) diff --git a/modules/source/spatrectif_000.c b/modules/source/spatrectif_000.c index e5fe5d3..3e8b2ec 100644 --- a/modules/source/spatrectif_000.c +++ b/modules/source/spatrectif_000.c @@ -1,321 +1,366 @@ -/* spatrectif_000.c */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -int errno; -#include "drp_structs.h" -#include "fitsio.h" - -float blame[DATA][numspec][MAXSLICE]; // Kernal for applying blame to lenslets for the first portion of the iterations -float weight[numspec][DATA]; // Weight normalization factor for distributing blame -float good[numspec][DATA]; // Weighted version of quality to decide if a pixel is valid. -float cmp[numspec][DATA]; // Weights for good -float fblame[DATA][numspec][MAXSLICE]; // Kernal for applying blame during the final iterations -float fbasisv[DATA][numspec][MAXSLICE]; // influence matrix (basis vector) -float fweight[numspec][DATA]; // Weight factor for distributing final blame -float *fw; // Weight factor for distributing final blame -float *w; // Weight factor for distributing final blame -float t_image[DATA][numspec]; // Temporary extracted spectra -float c_image[DATA][numspec]; // Cumulative extracted spectra -float *ti; // For each column accumulate delta to lenslet value -float *ci; // For each column best current estimate of lenslet value -float t_fbv[MAXSLICE]; // Temporary basis vector to manipulate if necessary -float dummy[DATA]; // Model of the flux given current lenslet estimate -float residual[DATA]; // Residual between data - model - -int where; - -int spatrectif_000(int argc, void* argv[]) -{ - // Parameters input in IDL calling program - short int totalParmCount; - short int numiter; - float scale; - float relaxation; - float relax; - short int basesize; - short int (*hilo)[2]; - short int bottom[numspec]; - short int *effective; - float *bv; - float *bl, *fbl, *fbv; - float (*basis_vectors)[MAXSLICE][DATA]; - Module *pModule; - DataSet *pDataSet; - float (*Frames)[DATA]; - IDL_STRING *Headers; - float (*Noise)[DATA]; - unsigned char (*Quality)[DATA]; - float (*image)[DATA]; - float (*noise)[DATA]; - unsigned char (*quality)[DATA]; - // End of parameters input in IDL calling program - - double t1, t2, t3, t4, t5, t6; // Time counters used to examine execution time - short int i=0, ii=0, j=0, jj=0, sp=0, l=0; - long int in1=0, in2=0, in3=0, index[numspec]; - - // Make a temporary residual matrix - long naxes[3]; - naxes[0] = DATA; - naxes[1] = DATA; - naxes[2] = (21); - // ******************************* - - // These parameters should be set in the same order as theay are passed - // from the IDL code. This is not yet automated, and I'm not sure how to - // do it. - i = 0; - totalParmCount = *(short int * )argv[i++]; - numiter = *(short int * )argv[i++]; - relaxation = *(float * )argv[i++]; - basesize = *(short int * )argv[i++]; - hilo = (short int (*)[2] )argv[i++]; - effective = (short int * )argv[i++]; - basis_vectors = (float (*)[MAXSLICE][DATA])argv[i++]; - pModule = (Module * )argv[i++]; - pDataSet = (DataSet * )argv[i++]; - Frames = (float (*)[DATA] )argv[i++]; - Headers = (IDL_STRING * )argv[i++]; - Noise = (float (*)[DATA] )argv[i++]; - Quality = (unsigned char (*)[DATA] )argv[i++]; - image = (float (*)[DATA] )argv[i++]; - noise = (float (*)[DATA] )argv[i++]; - quality = (unsigned char (*)[DATA] )argv[i++]; - scale = *(float * )argv[i++]; // Plate scale - - /* - * Start placing items from the original rectification code here. - * This code will rectify an input data frame. - */ - // printf("Image is spatrectif is %x.\n",image); - printf( "spatrectif_000.c: Now processing RAW data...\n"); - printf("Number of iterations = %d.\n",numiter); - /* printf("Relaxation Parameter = %f.\n",relaxation); - printf("Platescale = %f.\n",scale); */ - (void)fflush(stdout); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Iterate on each lenslet. - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - memset((void *)image,0, IMAGECNT ); // initial guess of the solution - - t1=systime(); - - // Set the local pointer bv equal to the address of the lowest member of the basis_vector - bv = basis_vectors[0][0]; - for (sp = 0; sp< numspec; sp++) - { - bottom[sp]=hilo[sp][0]; - } - - printf("Calculating weights.\n"); - (void)fflush(stdout); - - // Making local copies of the basis vector - for (sp=0; sp 0.0 ) - { - for (l=0; l 0.0 ) - { - for (l=0; l 14) { - // After first set of iterations, settle down to stable solution - ti[sp]+=relax*residual[j]* fbl[jj]; - //ti[sp]+=relax*residual[j]* bl[jj]; - } - j++; // march up column to match PSF location - } - } - // Correction is applied. - for (sp=0; sp< numspec; sp++) { - if ( weight[sp][i] > 0.0 ) { - ci[sp]+= ti[sp]; - } - } - - } // end of loop over image. - } // for each iteration - - for (i=0; i 4.5 ) quality[sp][i]=9; // Implies over have of the weights have good pixels. - } // for each spectral channel sp ... - } - printf("\n"); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - t2 = systime(); - // (void)printf("Total Time = %lf\n", t2-t1 ); - // writefitsimagefile("!/sdata1101/osiris-data/osiris2/050224/SPEC/ORP/resid.fits", FLOAT_IMG, 3, naxes, resid); - // free(resid); - - -#ifdef SAVE_INTERMEDIATE_FILES - //naxes[0] = DATA; - //naxes[1] = numspec; - //writefitsimagefile("testimage.fits", FLOAT_IMG, 2, naxes, image); -#endif - // printf("Image is spatrectif is %x.\n",image); - - return 0; -} +/* spatrectif_000.c */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +int errno; +#include "drp_structs.h" +#include "fitsio.h" + +float blame[DATA][numspec][MAXSLICE]; // Kernal for applying blame to lenslets for the first portion of the iterations +float weight[numspec][DATA]; // Weight normalization factor for distributing blame +float good[numspec][DATA]; // Weighted version of quality to decide if a pixel is valid. +float cmp[numspec][DATA]; // Weights for good + +float fblame[DATA][numspec][MAXSLICE]; // Kernal for applying blame during the final iterations +float fbasisv[DATA][numspec][MAXSLICE]; // influence matrix +float fweight[numspec][DATA]; // Weight factor for distributing final blame +float *fw; // Weight factor for distributing final blame +float *w; // Weight factor for distributing final blame +float t_image[DATA][numspec]; +float c_image[DATA][numspec]; +float *ti; +float *ci; +float dummy[DATA]; +float residual[DATA]; +float new_image[numspec][DATA]; +float adjust[MAXSLICE]; +float maxi; +int where; + +int spatrectif_000(int argc, void* argv[]) +{ + // Parameters input in IDL calling program + short int totalParmCount; + short int numiter; + float scale; + float relaxation; + float relax; + short int basesize; + short int (*hilo)[2]; + short int bottom[numspec]; + short int *effective; + float *bv; + float *bl, *fbl, *fbv; + float (*basis_vectors)[MAXSLICE][DATA]; + Module *pModule; + DataSet *pDataSet; + float (*Frames)[DATA]; + IDL_STRING *Headers; + float (*Noise)[DATA]; + unsigned char (*Quality)[DATA]; + float (*image)[DATA]; + float (*noise)[DATA]; + unsigned char (*quality)[DATA]; + // End of parameters input in IDL calling program + + double t1, t2, t3, t4, t5, t6; // Time counters used to examine execution time + short int i=0, ii=0, j=0, jj=0, sp=0, l=0; + long int in1=0, in2=0, in3=0, index[numspec]; + + // Make a temporary residual matrix + long naxes[3]; + naxes[0] = DATA; + naxes[1] = DATA; + naxes[2] = (21); + // ******************************* + + // These parameters should be set in the same order as theay are passed + // from the IDL code. This is not yet automated, and I'm not sure how to + // do it. + i = 0; + totalParmCount = *(short int * )argv[i++]; + numiter = *(short int * )argv[i++]; + relaxation = *(float * )argv[i++]; + basesize = *(short int * )argv[i++]; + hilo = (short int (*)[2] )argv[i++]; + effective = (short int * )argv[i++]; + basis_vectors = (float (*)[MAXSLICE][DATA])argv[i++]; + pModule = (Module * )argv[i++]; + pDataSet = (DataSet * )argv[i++]; + Frames = (float (*)[DATA] )argv[i++]; + Headers = (IDL_STRING * )argv[i++]; + Noise = (float (*)[DATA] )argv[i++]; + Quality = (unsigned char (*)[DATA] )argv[i++]; + image = (float (*)[DATA] )argv[i++]; + noise = (float (*)[DATA] )argv[i++]; + quality = (unsigned char (*)[DATA] )argv[i++]; + scale = *(float * )argv[i++]; // Plate scale + + /* + * Start placing items from the original rectification code here. + * This code will rectify an input data frame. + */ + // printf("Image is spatrectif is %x.\n",image); + printf( "spatrectif_000.c: Now processing RAW data...\n"); + printf("Number of iterations = %d.\n",numiter); + /* printf("Relaxation Parameter = %f.\n",relaxation); + printf("Platescale = %f.\n",scale); */ + (void)fflush(stdout); + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Iterate on each lenslet. + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + memset((void *)image,0, IMAGECNT ); // initial guess of the solution + + t1=systime(); + + // Set the local pointer bv equal to the address of the lowest member of the basis_vector + bv = basis_vectors[0][0]; + for (sp = 0; sp< numspec; sp++) + { + bottom[sp]=hilo[sp][0]; + } + + printf("Calculating weights.\n"); + (void)fflush(stdout); + + for (sp=0; sp maxi ) + // { + // maxi = fbv[l]; + // where = l; + // } + // } + // bl[where]= maxi; // Alternative initial blame where only peak pixel is used. + // w[i] = bl[where]; + for (l=0; l 0.0 ) + { + for (l=0; l 0.0 ) + { + for (l=0; l 14) { + // After first set of iterations, settle down to stable solution + ci[sp]+=relax*residual[j]* fbl[jj]; + } + // if ( ii > 20) { + // After first set of iterations, settle down to stable solution + // ci[sp]+=2.0*relax*residual[j]* fbl[jj]; + //} + j++; + } + } + // For the first iterations, correction is sharpened and applied. + if ( ii < 10 ) { + for (sp=0; sp< numspec; sp++) { // Correct middle rows + if ( weight[sp][i] > 0.0 ) {// Increment guess image + ci[sp]+= ti[sp]; + if ( sp > 63 ) + if ( weight[sp-64][i] > 0.0 ) // Increment guess image + ci[sp]+= 0.4*(ti[sp]-ti[sp-64]); + if ( sp < numspec-64 ) + if ( weight[sp+64][i] > 0.0 ) // Increment guess image + ci[sp]+= 0.4*(ti[sp]-ti[sp+64]); + } + } + } + if ( (ii > 9) && (ii < 15) ) { + for (sp=0; sp< numspec; sp++) { // Correct middle rows + if ( weight[sp][i] > 0.0 ) {// Increment guess image + ci[sp]+= ti[sp]; + if ( sp > 63 ) + if ( weight[sp-64][i] > 0.0 ) // Increment guess image + ci[sp]-= 0.2*(ti[sp]-ti[sp-64]); + if ( sp < numspec-64 ) + if ( weight[sp+64][i] > 0.0 ) // Increment guess image + ci[sp]-= 0.2*(ti[sp]-ti[sp+64]); + } + } + } + + // if ( (scale > 0.099) && (ii < 2) ) { + + } // end of loop over image. + // Smooth early iterations along the spectral direction. + if ( (ii > 11) && (ii < 18) ) + { + for (sp = 0; sp 4.5 ) quality[sp][i]=9; // Implies over have of the weights have good pixels. + } // for each spectral channel sp ... + } + printf("\n"); + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + t2 = systime(); + // (void)printf("Total Time = %lf\n", t2-t1 ); + + // writefitsimagefile("!/sdata1101/osiris-data/osiris2/050224/SPEC/ORP/resid.fits", FLOAT_IMG, 3, naxes, resid); + // free(resid); + + +#ifdef SAVE_INTERMEDIATE_FILES + //naxes[0] = DATA; + //naxes[1] = numspec; + //writefitsimagefile("testimage.fits", FLOAT_IMG, 2, naxes, image); +#endif + // printf("Image is spatrectif is %x.\n",image); + + return 0; +} From 13065036e1042363256bcf937564a6973cb504ae Mon Sep 17 00:00:00 2001 From: Samantha Chappell Date: Thu, 8 Sep 2016 13:53:24 -0700 Subject: [PATCH 17/22] The code that computes the RMS across spaxels for a given wavelength range. --- tests/drptestbones/checkSkylines.py | 39 +++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 tests/drptestbones/checkSkylines.py diff --git a/tests/drptestbones/checkSkylines.py b/tests/drptestbones/checkSkylines.py new file mode 100644 index 0000000..318d094 --- /dev/null +++ b/tests/drptestbones/checkSkylines.py @@ -0,0 +1,39 @@ +import pyfits +import sys +import numpy as np + +def checkSkylines(filename,min_lambda,max_lambda,sumRow,outname): + + file = pyfits.open(filename) + data = file[0].data + hdr = file[0].header + sizeLambda = data.shape[2] + + min_image_lambda = hdr['CRVAL1'] + delta_lambda = hdr['CDELT1'] + + min_pixel = int(round((min_lambda - min_image_lambda) / delta_lambda)) + max_pixel = int(round((max_lambda - min_image_lambda) / delta_lambda)) + + if (min_pixel >= max_pixel): + sys.exit("Max lambda needs to be greater than min lambda") + + if (min_pixel < 0): + sys.exit("Selected minimum wavelength is off detector") + + if (max_pixel > sizeLambda): + sys.exit("Selected maximum wavelength is off detector") + + onlyLine = data[:,:,min_pixel:max_pixel] + + sumLine = np.sum(onlyLine,axis=2) + + sumRowOnly = sumLine[:,sumRow] + + totalRMS = np.std(sumLine) + rowRMS = np.std(sumRowOnly) + + + pyfits.writeto(outname,sumLine,clobber=True) + + return totalRMS, rowRMS From e98fc12a4ee9fc9c0ca7578a55429529c5314c32 Mon Sep 17 00:00:00 2001 From: Samantha Chappell Date: Thu, 8 Sep 2016 13:54:33 -0700 Subject: [PATCH 18/22] Reduces the sky frame, the resulting .fits will be compared to a _ref.fits --- tests/test_checkskylines/sky.xml | 12 ++++++++++++ 1 file changed, 12 insertions(+) create mode 100644 tests/test_checkskylines/sky.xml diff --git a/tests/test_checkskylines/sky.xml b/tests/test_checkskylines/sky.xml new file mode 100644 index 0000000..3dc84eb --- /dev/null +++ b/tests/test_checkskylines/sky.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + From 66d1795bfcfd36ddb2694a08e165b2fbdd98bd28 Mon Sep 17 00:00:00 2001 From: Samantha Chappell Date: Thu, 8 Sep 2016 13:56:28 -0700 Subject: [PATCH 19/22] Test code that calls the xml and the python code to check the reduction and the RMS for a wavelength range across all spaxels in the sky frame --- .../test_checkskylines/test_checkskylines.py | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) create mode 100644 tests/test_checkskylines/test_checkskylines.py diff --git a/tests/test_checkskylines/test_checkskylines.py b/tests/test_checkskylines/test_checkskylines.py new file mode 100644 index 0000000..4cba850 --- /dev/null +++ b/tests/test_checkskylines/test_checkskylines.py @@ -0,0 +1,26 @@ +# -*- coding: utf-8 -*- + +import pytest +import os +import glob + +from drptestbones.backbone import consume_queue_directory +from drptestbones.diff import fits_osiris_allclose +from drptestbones.fetchdata import get_test_file + +from drptestbones.checkSkylines import checkSkylines + + +def test_skyline(drf_queue): + """Test FITS sky lines""" + consume_queue_directory(drf_queue) + output_file = os.path.join(drf_queue, "s160711_a013002_Kbb_035.fits") + expected_file = os.path.join(drf_queue, "s160711_a013002_Kbb_035_ref.fits") + rms_fits = os.path.join(drf_queue,'s160711_a013002_Kbb_035_RMS.fits') + + totalRMS,lineRMS = checkSkylines(output_file,2071,2075,41,rms_fits) + fractRMS = lineRMS/totalRMS + + fits_osiris_allclose(output_file, expected_file) + assert ((fractRMS > 0.97) & (fractRMS <= 1.0)) + From 0cbda916c666921a954f4bb50b4cbc412010c112 Mon Sep 17 00:00:00 2001 From: Jim Lyke Date: Thu, 8 Sep 2016 11:30:06 -1000 Subject: [PATCH 20/22] How to install the DRP --- EXAMPLE.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/EXAMPLE.md b/EXAMPLE.md index bdfb873..38fa82e 100644 --- a/EXAMPLE.md +++ b/EXAMPLE.md @@ -283,5 +283,8 @@ tests/test_emission_line/test_emission_line.py . ========================== 1 passed in 93.84 seconds =========================== ``` ## Test Successful +<<<<<<< HEAD +======= +>>>>>>> e7e9ac5fcd93f8c47d191a24cb0360fc67d90e32 From d628cc81fd3197f72283f35db5b061f075797090 Mon Sep 17 00:00:00 2001 From: Tuan Do Date: Fri, 9 Sep 2016 13:24:45 -0700 Subject: [PATCH 21/22] slight change for bash instructions --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 7ac0abc..bf7bc6a 100644 --- a/README.md +++ b/README.md @@ -65,8 +65,8 @@ You can add these lines to your ``.bashrc`` file or other startup profile if you ``` OSIRIS_VERBOSE=0 -source scripts/osirisSetup.sh -osirisSetup /my/path/to/osiris/drp/ +source /my/path/to/osiris/drp/scripts/osirisSetup.sh +osirisSetup /my/path/to/osiris/drp ``` ### Environment Setup in CSH From 2443a5b3a687d2de500bdf83aca57d3973c3002a Mon Sep 17 00:00:00 2001 From: Anna Boehle Date: Tue, 27 Sep 2016 12:53:53 -0700 Subject: [PATCH 22/22] Skeleton version of spatrectif_000.c --- modules/source/spatrectif_000.c | 687 +++++++++++++++----------------- 1 file changed, 321 insertions(+), 366 deletions(-) diff --git a/modules/source/spatrectif_000.c b/modules/source/spatrectif_000.c index 3e8b2ec..e5fe5d3 100644 --- a/modules/source/spatrectif_000.c +++ b/modules/source/spatrectif_000.c @@ -1,366 +1,321 @@ -/* spatrectif_000.c */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -int errno; -#include "drp_structs.h" -#include "fitsio.h" - -float blame[DATA][numspec][MAXSLICE]; // Kernal for applying blame to lenslets for the first portion of the iterations -float weight[numspec][DATA]; // Weight normalization factor for distributing blame -float good[numspec][DATA]; // Weighted version of quality to decide if a pixel is valid. -float cmp[numspec][DATA]; // Weights for good - -float fblame[DATA][numspec][MAXSLICE]; // Kernal for applying blame during the final iterations -float fbasisv[DATA][numspec][MAXSLICE]; // influence matrix -float fweight[numspec][DATA]; // Weight factor for distributing final blame -float *fw; // Weight factor for distributing final blame -float *w; // Weight factor for distributing final blame -float t_image[DATA][numspec]; -float c_image[DATA][numspec]; -float *ti; -float *ci; -float dummy[DATA]; -float residual[DATA]; -float new_image[numspec][DATA]; -float adjust[MAXSLICE]; -float maxi; -int where; - -int spatrectif_000(int argc, void* argv[]) -{ - // Parameters input in IDL calling program - short int totalParmCount; - short int numiter; - float scale; - float relaxation; - float relax; - short int basesize; - short int (*hilo)[2]; - short int bottom[numspec]; - short int *effective; - float *bv; - float *bl, *fbl, *fbv; - float (*basis_vectors)[MAXSLICE][DATA]; - Module *pModule; - DataSet *pDataSet; - float (*Frames)[DATA]; - IDL_STRING *Headers; - float (*Noise)[DATA]; - unsigned char (*Quality)[DATA]; - float (*image)[DATA]; - float (*noise)[DATA]; - unsigned char (*quality)[DATA]; - // End of parameters input in IDL calling program - - double t1, t2, t3, t4, t5, t6; // Time counters used to examine execution time - short int i=0, ii=0, j=0, jj=0, sp=0, l=0; - long int in1=0, in2=0, in3=0, index[numspec]; - - // Make a temporary residual matrix - long naxes[3]; - naxes[0] = DATA; - naxes[1] = DATA; - naxes[2] = (21); - // ******************************* - - // These parameters should be set in the same order as theay are passed - // from the IDL code. This is not yet automated, and I'm not sure how to - // do it. - i = 0; - totalParmCount = *(short int * )argv[i++]; - numiter = *(short int * )argv[i++]; - relaxation = *(float * )argv[i++]; - basesize = *(short int * )argv[i++]; - hilo = (short int (*)[2] )argv[i++]; - effective = (short int * )argv[i++]; - basis_vectors = (float (*)[MAXSLICE][DATA])argv[i++]; - pModule = (Module * )argv[i++]; - pDataSet = (DataSet * )argv[i++]; - Frames = (float (*)[DATA] )argv[i++]; - Headers = (IDL_STRING * )argv[i++]; - Noise = (float (*)[DATA] )argv[i++]; - Quality = (unsigned char (*)[DATA] )argv[i++]; - image = (float (*)[DATA] )argv[i++]; - noise = (float (*)[DATA] )argv[i++]; - quality = (unsigned char (*)[DATA] )argv[i++]; - scale = *(float * )argv[i++]; // Plate scale - - /* - * Start placing items from the original rectification code here. - * This code will rectify an input data frame. - */ - // printf("Image is spatrectif is %x.\n",image); - printf( "spatrectif_000.c: Now processing RAW data...\n"); - printf("Number of iterations = %d.\n",numiter); - /* printf("Relaxation Parameter = %f.\n",relaxation); - printf("Platescale = %f.\n",scale); */ - (void)fflush(stdout); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Iterate on each lenslet. - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - memset((void *)image,0, IMAGECNT ); // initial guess of the solution - - t1=systime(); - - // Set the local pointer bv equal to the address of the lowest member of the basis_vector - bv = basis_vectors[0][0]; - for (sp = 0; sp< numspec; sp++) - { - bottom[sp]=hilo[sp][0]; - } - - printf("Calculating weights.\n"); - (void)fflush(stdout); - - for (sp=0; sp maxi ) - // { - // maxi = fbv[l]; - // where = l; - // } - // } - // bl[where]= maxi; // Alternative initial blame where only peak pixel is used. - // w[i] = bl[where]; - for (l=0; l 0.0 ) - { - for (l=0; l 0.0 ) - { - for (l=0; l 14) { - // After first set of iterations, settle down to stable solution - ci[sp]+=relax*residual[j]* fbl[jj]; - } - // if ( ii > 20) { - // After first set of iterations, settle down to stable solution - // ci[sp]+=2.0*relax*residual[j]* fbl[jj]; - //} - j++; - } - } - // For the first iterations, correction is sharpened and applied. - if ( ii < 10 ) { - for (sp=0; sp< numspec; sp++) { // Correct middle rows - if ( weight[sp][i] > 0.0 ) {// Increment guess image - ci[sp]+= ti[sp]; - if ( sp > 63 ) - if ( weight[sp-64][i] > 0.0 ) // Increment guess image - ci[sp]+= 0.4*(ti[sp]-ti[sp-64]); - if ( sp < numspec-64 ) - if ( weight[sp+64][i] > 0.0 ) // Increment guess image - ci[sp]+= 0.4*(ti[sp]-ti[sp+64]); - } - } - } - if ( (ii > 9) && (ii < 15) ) { - for (sp=0; sp< numspec; sp++) { // Correct middle rows - if ( weight[sp][i] > 0.0 ) {// Increment guess image - ci[sp]+= ti[sp]; - if ( sp > 63 ) - if ( weight[sp-64][i] > 0.0 ) // Increment guess image - ci[sp]-= 0.2*(ti[sp]-ti[sp-64]); - if ( sp < numspec-64 ) - if ( weight[sp+64][i] > 0.0 ) // Increment guess image - ci[sp]-= 0.2*(ti[sp]-ti[sp+64]); - } - } - } - - // if ( (scale > 0.099) && (ii < 2) ) { - - } // end of loop over image. - // Smooth early iterations along the spectral direction. - if ( (ii > 11) && (ii < 18) ) - { - for (sp = 0; sp 4.5 ) quality[sp][i]=9; // Implies over have of the weights have good pixels. - } // for each spectral channel sp ... - } - printf("\n"); - ////////////////////////////////////////////////////////////////////////////////////////////////////////////// - t2 = systime(); - // (void)printf("Total Time = %lf\n", t2-t1 ); - - // writefitsimagefile("!/sdata1101/osiris-data/osiris2/050224/SPEC/ORP/resid.fits", FLOAT_IMG, 3, naxes, resid); - // free(resid); - - -#ifdef SAVE_INTERMEDIATE_FILES - //naxes[0] = DATA; - //naxes[1] = numspec; - //writefitsimagefile("testimage.fits", FLOAT_IMG, 2, naxes, image); -#endif - // printf("Image is spatrectif is %x.\n",image); - - return 0; -} +/* spatrectif_000.c */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +int errno; +#include "drp_structs.h" +#include "fitsio.h" + +float blame[DATA][numspec][MAXSLICE]; // Kernal for applying blame to lenslets for the first portion of the iterations +float weight[numspec][DATA]; // Weight normalization factor for distributing blame +float good[numspec][DATA]; // Weighted version of quality to decide if a pixel is valid. +float cmp[numspec][DATA]; // Weights for good +float fblame[DATA][numspec][MAXSLICE]; // Kernal for applying blame during the final iterations +float fbasisv[DATA][numspec][MAXSLICE]; // influence matrix (basis vector) +float fweight[numspec][DATA]; // Weight factor for distributing final blame +float *fw; // Weight factor for distributing final blame +float *w; // Weight factor for distributing final blame +float t_image[DATA][numspec]; // Temporary extracted spectra +float c_image[DATA][numspec]; // Cumulative extracted spectra +float *ti; // For each column accumulate delta to lenslet value +float *ci; // For each column best current estimate of lenslet value +float t_fbv[MAXSLICE]; // Temporary basis vector to manipulate if necessary +float dummy[DATA]; // Model of the flux given current lenslet estimate +float residual[DATA]; // Residual between data - model + +int where; + +int spatrectif_000(int argc, void* argv[]) +{ + // Parameters input in IDL calling program + short int totalParmCount; + short int numiter; + float scale; + float relaxation; + float relax; + short int basesize; + short int (*hilo)[2]; + short int bottom[numspec]; + short int *effective; + float *bv; + float *bl, *fbl, *fbv; + float (*basis_vectors)[MAXSLICE][DATA]; + Module *pModule; + DataSet *pDataSet; + float (*Frames)[DATA]; + IDL_STRING *Headers; + float (*Noise)[DATA]; + unsigned char (*Quality)[DATA]; + float (*image)[DATA]; + float (*noise)[DATA]; + unsigned char (*quality)[DATA]; + // End of parameters input in IDL calling program + + double t1, t2, t3, t4, t5, t6; // Time counters used to examine execution time + short int i=0, ii=0, j=0, jj=0, sp=0, l=0; + long int in1=0, in2=0, in3=0, index[numspec]; + + // Make a temporary residual matrix + long naxes[3]; + naxes[0] = DATA; + naxes[1] = DATA; + naxes[2] = (21); + // ******************************* + + // These parameters should be set in the same order as theay are passed + // from the IDL code. This is not yet automated, and I'm not sure how to + // do it. + i = 0; + totalParmCount = *(short int * )argv[i++]; + numiter = *(short int * )argv[i++]; + relaxation = *(float * )argv[i++]; + basesize = *(short int * )argv[i++]; + hilo = (short int (*)[2] )argv[i++]; + effective = (short int * )argv[i++]; + basis_vectors = (float (*)[MAXSLICE][DATA])argv[i++]; + pModule = (Module * )argv[i++]; + pDataSet = (DataSet * )argv[i++]; + Frames = (float (*)[DATA] )argv[i++]; + Headers = (IDL_STRING * )argv[i++]; + Noise = (float (*)[DATA] )argv[i++]; + Quality = (unsigned char (*)[DATA] )argv[i++]; + image = (float (*)[DATA] )argv[i++]; + noise = (float (*)[DATA] )argv[i++]; + quality = (unsigned char (*)[DATA] )argv[i++]; + scale = *(float * )argv[i++]; // Plate scale + + /* + * Start placing items from the original rectification code here. + * This code will rectify an input data frame. + */ + // printf("Image is spatrectif is %x.\n",image); + printf( "spatrectif_000.c: Now processing RAW data...\n"); + printf("Number of iterations = %d.\n",numiter); + /* printf("Relaxation Parameter = %f.\n",relaxation); + printf("Platescale = %f.\n",scale); */ + (void)fflush(stdout); + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Iterate on each lenslet. + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + memset((void *)image,0, IMAGECNT ); // initial guess of the solution + + t1=systime(); + + // Set the local pointer bv equal to the address of the lowest member of the basis_vector + bv = basis_vectors[0][0]; + for (sp = 0; sp< numspec; sp++) + { + bottom[sp]=hilo[sp][0]; + } + + printf("Calculating weights.\n"); + (void)fflush(stdout); + + // Making local copies of the basis vector + for (sp=0; sp 0.0 ) + { + for (l=0; l 0.0 ) + { + for (l=0; l 14) { + // After first set of iterations, settle down to stable solution + ti[sp]+=relax*residual[j]* fbl[jj]; + //ti[sp]+=relax*residual[j]* bl[jj]; + } + j++; // march up column to match PSF location + } + } + // Correction is applied. + for (sp=0; sp< numspec; sp++) { + if ( weight[sp][i] > 0.0 ) { + ci[sp]+= ti[sp]; + } + } + + } // end of loop over image. + } // for each iteration + + for (i=0; i 4.5 ) quality[sp][i]=9; // Implies over have of the weights have good pixels. + } // for each spectral channel sp ... + } + printf("\n"); + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + t2 = systime(); + // (void)printf("Total Time = %lf\n", t2-t1 ); + // writefitsimagefile("!/sdata1101/osiris-data/osiris2/050224/SPEC/ORP/resid.fits", FLOAT_IMG, 3, naxes, resid); + // free(resid); + + +#ifdef SAVE_INTERMEDIATE_FILES + //naxes[0] = DATA; + //naxes[1] = numspec; + //writefitsimagefile("testimage.fits", FLOAT_IMG, 2, naxes, image); +#endif + // printf("Image is spatrectif is %x.\n",image); + + return 0; +}