IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 34442


Ignore:
Timestamp:
Sep 14, 2012, 6:10:32 PM (14 years ago)
Author:
watersc1
Message:

Still doesn't work correctly. I've simplified and corrected the problems I was having before, but I still have residuals of ~50 pixels in the warped background model relative to the science warp. I'm beginning to think that there's some issue with the interpolated background model, and that it's not quite the correct size? That's pretty much all I can think of at this point.

Location:
branches/czw_branch/20120906
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20120906/ippconfig/recipes/nightly_science.config

    r34065 r34442  
    223223  DIFFABLE  BOOL TRUE
    224224  REDUCTION STR SWEETSPOT
     225  ADDITIONAL_STACK_LABEL STR ecliptic.rp
    225226#  REDUCTION STR SSTF_4SIG
    226227#  WARP      S16 60
    … …  
    234235  DIFFABLE    BOOL TRUE
    235236  REDUCTION STR SWEETSPOT
     237  ADDITIONAL_STACK_LABEL STR ecliptic.rp
    236238#  REDUCTION STR SSTF_4SIG
    237239END
    … …  
    244246  DIFFABLE    BOOL TRUE
    245247  REDUCTION STR SWEETSPOT
     248  ADDITIONAL_STACK_LABEL STR ecliptic.rp
    246249#  REDUCTION STR SSTF_4SIG
    247250END
    … …  
    254257  DIFFABLE    BOOL TRUE
    255258  REDUCTION STR SWEETSPOT
     259  ADDITIONAL_STACK_LABEL STR ecliptic.rp
    256260#  REDUCTION STR SSTF_4SIG
    257261END

  • branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.c

    r34414 r34442  
    189189    switch (mode) {
    190190      case PS_INTERPOLATE_FLAT:
     191      case PS_INTERPOLATE_BILINEAR_SIMPLE:
    191192        // Nothing to pre-compute
    192193        break;
    … …  
    287288# endif
    288289      case PS_INTERPOLATE_BIQUADRATIC:
     290      case PS_INTERPOLATE_BILINEAR_SIMPLE:
    289291        // 2D kernel, would cost too much memory to pre-calculate
    290292        break;
     293
    291294      default:
    292295        psAbort("Unsupported interpolation mode: %x", mode);
    … …  
    937940
    938941
     942psImageInterpolateStatus interpolateJustWork(double *imageValue, double *varianceValue,
     943                                             psImageMaskType *maskValue, float x, float y,
     944                                             const psImageInterpolation *interp) {
     945  float u1,u2;
     946  int xl,xh;
     947  int yl,yh;
     948  const psImage *image = interp->image; // Image to interpolate
     949  // const psImage *mask = interp->mask; // Mask to interpolate
     950  const psImage *variance = interp->variance; // Variance to interpolate
     951
     952  xl = floor(x);
     953  if (xl < 0) { xl = 0; }
     954  if (xl >= image->numCols) {xl = image->numCols - 1; }
     955  xh = xl + 1;
     956  if (xh >= image->numCols) {xh = image->numCols - 1; }
     957  yl = floor(y);
     958  if (yl < 0) { yl = 0; }
     959  if (yl >= image->numRows) {yl = image->numRows - 1; }
     960  yh = yl + 1;
     961  if (yh >= image->numRows) {yh = image->numRows - 1; }
     962  if (varianceValue && variance) {
     963    u1 = (xh - x) * variance->data.F32[yl][xl] + (x - xl) * variance->data.F32[yl][xh];
     964    u2 = (xh - x) * variance->data.F32[yh][xl] + (x - xl) * variance->data.F32[yh][xh];
     965
     966    *varianceValue = (yh - y) * u1 + (y - yl) * u2;
     967  }
     968#if 0
     969  if (maskValue && mask) {
     970    u1 = (xh - x) * mask->data.F32[yl][xl] + (x - xl) * mask->data.F32[yl][xh];
     971    u2 = (xh - x) * mask->data.F32[yh][xl] + (x - xl) * mask->data.F32[yh][xh];
     972
     973    maskValue = (yh - y) * u1 + (y - yl) * u2;
     974  }
     975#endif
     976  if (imageValue && image) {
     977    if (image->data.F32[yl][xl] == image->data.F32[yl][xh]) {
     978      u1 = image->data.F32[yl][xh];
     979    }
     980    else {
     981      u1 = (xh - x) * image->data.F32[yl][xl] + (x - xl) * image->data.F32[yl][xh];
     982    }
     983    if (image->data.F32[yh][xl] == image->data.F32[yh][xh]) {
     984      u2 = image->data.F32[yh][xh];
     985    }
     986    else {
     987      u2 = (xh - x) * image->data.F32[yh][xl] + (x - xl) * image->data.F32[yh][xh];
     988    }
     989    if (u1 == u2) {
     990      *imageValue = u1;
     991    }
     992    else {
     993      *imageValue = (yh - y) * u1 + (y - yl) * u2;
     994    }
     995    if ((floor(x) >= image->numCols)||
     996        (floor(y) >= image->numRows)) {
     997      *imageValue = NAN;
     998    }
     999  }
     1000  if (*imageValue == 0.0) {
     1001    fprintf(stderr,"IJK: Zero!: %g %g [%d %d %d %d] %g %g (%g %g %g %g)\n",
     1002            x,y,xl,xh,yl,yh,u1,u2,
     1003            image->data.F32[yl][xl],image->data.F32[yl][xh],image->data.F32[yh][xl],image->data.F32[yh][xh]
     1004            );
     1005  }
     1006  return PS_INTERPOLATE_STATUS_GOOD;
     1007}
     1008 
    9391009
    9401010psImageInterpolateStatus psImageInterpolate(double *imageValue, double *varianceValue,
    … …  
    9691039      case PS_INTERPOLATE_BIQUADRATIC:
    9701040        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
     1041      case PS_INTERPOLATE_BILINEAR_SIMPLE:
     1042        return interpolateJustWork(imageValue, varianceValue, maskValue, x, y, interp);
    9711043      case PS_INTERPOLATE_BILINEAR:
    9721044# if (!IS_BILIN_SEPARABLE)
    973         return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
     1045        return interpolateKernel(imageValue, varianceValue, maskValue, x, y, interp);
    9741046# endif
    9751047      case PS_INTERPOLATE_GAUSS:
    … …  
    10191091    switch (mode) {
    10201092      case PS_INTERPOLATE_FLAT:
     1093      case PS_INTERPOLATE_BILINEAR_SIMPLE:
    10211094        // No smearing by design
    10221095        return 1.0;
    … …  
    10831156    if (!strcasecmp(name, "FLAT"))     return PS_INTERPOLATE_FLAT;
    10841157    if (!strcasecmp(name, "BILINEAR")) return PS_INTERPOLATE_BILINEAR;
     1158    if (!strcasecmp(name, "SIMPLEBILINEAR")) return PS_INTERPOLATE_BILINEAR_SIMPLE;
    10851159    if (!strcasecmp(name, "BIQUADRATIC")) return PS_INTERPOLATE_BIQUADRATIC;
    10861160    if (!strcasecmp(name, "GAUSS"))    return PS_INTERPOLATE_GAUSS;

  • branches/czw_branch/20120906/psLib/src/imageops/psImageInterpolate.h

    r21281 r34442  
    3333    PS_INTERPOLATE_LANCZOS3,            ///< Sinc interpolation with 6x6 pixel kernel
    3434    PS_INTERPOLATE_LANCZOS4,            ///< Sinc interpolation with 8x8 pixel kernel
     35    PS_INTERPOLATE_BILINEAR_SIMPLE,     ///< Simple manual bilinear interpolation
    3536} psImageInterpolateMode;
    3637

  • branches/czw_branch/20120906/psModules/src/detrend/pmDark.c

    r33089 r34442  
    353353
    354354    // retrieve the required parameter vectors
    355     psArray *values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
    356     psAssert(values, "values not supplied");
     355    psArray *in_values  = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES");
     356    psAssert(in_values, "values not supplied");
    357357    psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK");
    358358    psAssert(roMask, "roMask not supplied");
    359     psVector *orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
    360     psAssert(orders, "orders not supplied");
    361 
    362     psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting
     359    psVector *max_orders = psMetadataLookupPtr(&mdok, output->analysis, "DARK.ORDERS");
     360    psAssert(max_orders, "orders not supplied");
     361
     362    psArray *values_set = psArrayAlloc(max_orders->n);
     363    psArray *poly_set = psArrayAlloc(max_orders->n);
     364    psVector *logL = psVectorAlloc(max_orders->n,PS_TYPE_F64);
     365
     366    for (int i = 0; i < max_orders->n; i++) {
     367      psVector *orders = psVectorAlloc(i+1,PS_TYPE_U8);
     368      for (int j = 0; j < orders->n; j++) {
     369        orders->data.U8[j] = max_orders->data.U8[j];
     370      }
     371      poly_set->data[i] =  psPolynomialMDAlloc(orders); // Polynomial for fitting
     372     
     373      psArray *values = psArrayAlloc(in_values->n);
     374     
     375      for (int j = 0; j < values->n; j++) {
     376        psVector *these_values = psVectorAlloc(i+1,PS_TYPE_F32);
     377        psVector *input_values = in_values->data[j];
     378
     379        for (int k = 0; k < orders->n; k++) {
     380          these_values->data.F32[k] = input_values->data.F32[k];
     381        }
     382        values->data[j] = these_values;
     383      }
     384      values_set->data[i] = values;
     385    }
    363386
    364387    // retrieve the norm vector, if supplied
    … …  
    383406    }
    384407
    385     pmDarkVisualInit(values);
     408    pmDarkVisualInit(values_set->data[max_orders->n - 1]);
    386409
    387410    pmReadout *outReadout = output->readouts->data[0];
    … …  
    423446            }
    424447
    425             if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
     448            int k_best = 0;
     449            for (int k = 0; k < max_orders->n; k++) {
     450              psPolynomialMD *poly = poly_set->data[k];
     451              psArray *values = values_set->data[k];
     452             
     453              if (!psPolynomialMDClipFit(poly, pixels, NULL, mask, 0xff, values, iter, rej)) {
    426454                psErrorClear();         // Nothing we can do about it
    427455                psVectorInit(poly->coeff, NAN);
    428             }
    429 
    430             pmDarkVisualPixelFit(pixels, mask);
    431             pmDarkVisualPixelModel(poly, values);
    432 
    433             for (int k = 0; k < poly->coeff->n; k++) {
     456              }
     457
     458              pmDarkVisualPixelFit(pixels, mask);
     459              pmDarkVisualPixelModel(poly, values);
     460
     461              // Insert math here to choose optimum model.
     462              logL->data.F64[k] = 0.0;
     463              psPolynomialMD *polySig = poly_set->data[0];
     464              for (int m = 0; m < poly->deviations->n; m++) {
     465                logL->data.F64[k] += pow(poly->deviations->data.F32[m] / polySig->stdevFit,2);
     466/*              if ((xOut == 20) && (yOut == 256)) { */
     467/*                psTrace("psModules.detrend",3,"pmDarkCombine DEV: %d %d: input %d models: Norders: %d logL: %g value: %g\n", */
     468/*                        xOut,yOut,m,k,logL->data.F64[k],poly->deviations->data.F32[m]); */
     469/*              } */
     470              }
     471              if (k > 0) {
     472                if ( ( logL->data.F64[k - 1] - logL->data.F64[k] ) > 1) { // Hard coded criterion for a ~5% limit with one degree of freedom
     473                  k_best = k;
     474                }
     475              }
     476              if ((xOut <= 600) && (yOut <= 600)) {
     477                psTrace("psModules.detrend",3,"pmDarkCombine: %d %d: models: Norders: %d logL: %g BestOrders: %d\n",
     478                        xOut,yOut,k,logL->data.F64[k],k_best);
     479              }
     480            }
     481            if (k_best > 1) {
     482              k_best = 1;
     483            }
     484/*          k_best = 1; */
     485            // Select the polynomial that seems best.
     486            psPolynomialMD *poly = poly_set->data[k_best];
     487             
     488            //            for (int k = 0; k < poly->coeff->n; k++) {
     489            for (int k = 0; k < max_orders->n + 1; k++) { // There is one more coefficient than is stored here.
    434490                pmReadout *ro = output->readouts->data[k]; // Readout of interest
    435                 ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
     491                if (k < poly->coeff->n) {
     492                  ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k];
     493                }
     494                else {
     495                  ro->image->data.F32[yOut][xOut] = 0.0;
     496                }
    436497            }
    437498            counts->data.U16[yOut][xOut] = poly->numFit;

  • branches/czw_branch/20120906/psModules/src/objects/Makefile.am

    r34403 r34442  
    111111        pmSourceOutputs.h \
    112112        pmSourceIO.h \
     113        pmSourceSatstar.h \
    113114        pmSourcePlots.h \
    114115        pmSourceVisual.h \

  • branches/czw_branch/20120906/psModules/src/objects/pmFootprintCullPeaks.c

    r34198 r34442  
    9191
    9292        // max flux is above threshold for brightest peak
    93         pmPeak *maxPeak = fp->peaks->data[0];
     93      pmPeak *maxPeak = NULL;
     94      for (int i = 0; i < fp->peaks->n; i++) {
     95        pmPeak *testPeak = fp->peaks->data[i];
     96        float this_peak = useSmoothedImage ? testPeak->smoothFlux : testPeak->rawFlux;
     97       
     98        if (isfinite(this_peak)) {
     99          maxPeak = fp->peaks->data[i];
     100          break;
     101        }
     102      }
     103      psAssert(maxPeak,"maxPeak was not set in these peaks");
     104      //      = fp->peaks->data[0];
    94105        float maxFlux = useSmoothedImage ? maxPeak->smoothFlux : maxPeak->rawFlux;
    95106

  • branches/czw_branch/20120906/pswarp/src/pswarpArguments.c

    r34411 r34442  
    170170      psWarning("BACKGROUND.MODEL is not set in the %s recipe -- defaulting to FALSE.", PSWARP_RECIPE);
    171171    }
    172     int bkgXgrid = psMetadataLookupBool(&status,recipe, "BKG.XGRID"); ///< Xsize of background model
     172    int bkgXgrid = psMetadataLookupS32(&status,recipe, "BKG.XGRID"); ///< Xsize of background model
    173173    if (!status) {
    174174      bkgXgrid = 10;
    175175      psWarning("BKG.XGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgXgrid);
    176176    }
    177     int bkgYgrid = psMetadataLookupBool(&status,recipe, "BKG.YGRID"); ///< Xsize of background model
     177    int bkgYgrid = psMetadataLookupS32(&status,recipe, "BKG.YGRID"); ///< Xsize of background model
    178178    if (!status) {
    179179      bkgYgrid = 10;
    180180      psWarning("BKG.YGRID is not set in the %s recipe -- defaultint to %d.",PSWARP_RECIPE,bkgYgrid);
    181181    }
    182    
     182
    183183   
    184184    // Set recipe values in the recipe (since we've possibly altered some)

  • branches/czw_branch/20120906/pswarp/src/pswarpLoop.c

    r34411 r34442  
    159159        }
    160160    }
    161 
     161    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
     162   
     163    fprintf(stderr,"WCSFI: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n",
     164            WCSF->crval1,WCSF->crval2,
     165            WCSF->crpix1,WCSF->crpix2,
     166            WCSF->cdelt1,WCSF->cdelt2,
     167            WCSF->trans->x->coeff[1][0],
     168            WCSF->trans->x->coeff[0][1],
     169            WCSF->trans->y->coeff[1][0],
     170            WCSF->trans->y->coeff[0][1],
     171            -1.0,-1.0);
     172    psFree(WCSF);
    162173    psList *cells = psListAlloc(NULL);  // List of cells, for concepts averaging
    163174
    … …  
    197208            }
    198209        }
     210    pmAstromWCS *WCSF = pmAstromWCSfromHeader(hdu->header);
     211   
     212    fprintf(stderr,"WCSFIB: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n",
     213            WCSF->crval1,WCSF->crval2,
     214            WCSF->crpix1,WCSF->crpix2,
     215            WCSF->cdelt1,WCSF->cdelt2,
     216            WCSF->trans->x->coeff[1][0],
     217            WCSF->trans->x->coeff[0][1],
     218            WCSF->trans->y->coeff[1][0],
     219            WCSF->trans->y->coeff[0][1],
     220            -1.0,-1.0);
     221    psFree(WCSF);
    199222
    200223       
    … …  
    574597    bool bilevelAstrometry = false;
    575598    pmHDU *phu = pmFPAviewThisPHU(view, astrom->fpa);
     599
     600    pmAstromWCS *WCSF = pmAstromWCSfromHeader(phu->header);
     601   
     602    fprintf(stderr,"WCSF: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n",
     603            WCSF->crval1,WCSF->crval2,
     604            WCSF->crpix1,WCSF->crpix2,
     605            WCSF->cdelt1,WCSF->cdelt2,
     606            WCSF->trans->x->coeff[1][0],
     607            WCSF->trans->x->coeff[0][1],
     608            WCSF->trans->y->coeff[1][0],
     609            WCSF->trans->y->coeff[0][1],
     610            -1.0,-1.0);
     611/* #define CORRECT_INPUT_WCS 1 */
     612/* #if CORRECT_INPUT_WCS */
     613/*      // Correct the input WCS */
     614/*      pmAstromWCS *WCS = pmAstromWCSfromHeader(phu->header); */
     615
     616/*      double cd1f = 1.0 * 400; */
     617/*      double cd2f = 1.0 * 400; */
     618/*      fprintf(stderr,"WCS: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", */
     619/*              WCS->crval1,WCS->crval2, */
     620/*              WCS->crpix1,WCS->crpix2, */
     621/*              WCS->cdelt1,WCS->cdelt2, */
     622/*              WCS->trans->x->coeff[1][0], */
     623/*              WCS->trans->x->coeff[0][1], */
     624/*              WCS->trans->y->coeff[1][0], */
     625/*              WCS->trans->y->coeff[0][1], */
     626/*              cd1f,cd2f); */
     627
     628/*      WCS->crpix1 = WCS->crpix1 / cd1f; */
     629/*      WCS->crpix2 = WCS->crpix2 / cd2f; */
     630
     631/*      WCS->cdelt1 *= cd1f; */
     632/*      WCS->cdelt2 *= cd2f; */
     633
     634/*      // WCS->trans->x->nX/nY */
     635/*      for (int q = 0; q < WCS->trans->x->nX; q++) { */
     636/*        for (int r = 0; r < WCS->trans->x->nY; r++) { */
     637/*          WCS->trans->x->coeff[q][r] *= cd1f; */
     638/*        } */
     639/*      } */
     640/*      for (int q = 0; q < WCS->trans->y->nX; q++) { */
     641/*        for (int r = 0; r < WCS->trans->y->nY; r++) { */
     642/*          WCS->trans->y->coeff[q][r] *= cd2f; */
     643/*        } */
     644/*      } */
     645/*      fprintf(stderr,"WCS: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n", */
     646/*              WCS->crval1,WCS->crval2, */
     647/*              WCS->crpix1,WCS->crpix2, */
     648/*              WCS->cdelt1,WCS->cdelt2, */
     649/*              WCS->trans->x->coeff[1][0], */
     650/*              WCS->trans->x->coeff[0][1], */
     651/*              WCS->trans->y->coeff[1][0], */
     652/*              WCS->trans->y->coeff[0][1], */
     653/*              cd1f,cd2f); */
     654/*      pmAstromWCStoHeader (phu->header,WCS); */
     655/*      // End WCS work. */
     656/* #endif */
     657
     658
     659   
     660   
    576661    if (phu) {
    577662        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
    … …  
    610695        // read WCS data from the corresponding header
    611696        pmHDU *hdu = pmFPAviewThisHDU (view, astrom->fpa);
    612         psMetadataAddS32(config->arguments,PS_LIST_TAIL,
    613                          "INTERPOLATION.MODE",PS_META_REPLACE,"", PS_INTERPOLATE_BILINEAR); ///< Mode
    614         psMetadataAddS32(config->arguments,PS_LIST_TAIL,
    615                          "INTERPOLATION.NUM",PS_META_REPLACE,"", 64); ///< Mode
     697
     698        pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
     699
     700        double cd1f = 1.0 * 400;
     701        double cd2f = 1.0 * 400;
     702        fprintf(stderr,"WCSC: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n",
     703                WCS->crval1,WCS->crval2,
     704                WCS->crpix1,WCS->crpix2,
     705                WCS->cdelt1,WCS->cdelt2,
     706                WCS->trans->x->coeff[1][0],
     707                WCS->trans->x->coeff[0][1],
     708                WCS->trans->y->coeff[1][0],
     709                WCS->trans->y->coeff[0][1],
     710                cd1f,cd2f);
     711
     712
     713        WCS->cdelt1 *= cd1f;
     714        WCS->cdelt2 *= cd2f;
     715        WCS->crpix1 = WCS->crpix1 / cd1f + 46.0 / 400.0;
     716        WCS->crpix2 = WCS->crpix2 / cd2f + 68.0 / 400.0;
     717
     718        // WCS->trans->x->nX/nY
     719        for (int q = 0; q < WCS->trans->x->nX; q++) {
     720          for (int r = 0; r < WCS->trans->x->nY; r++) {
     721            fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->x->coeff[q][r]);
     722            WCS->trans->x->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     723            fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->x->coeff[q][r]);
     724          }
     725        }
     726        for (int q = 0; q < WCS->trans->y->nX; q++) {
     727          for (int r = 0; r < WCS->trans->y->nY; r++) {
     728            fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->y->coeff[q][r]);
     729            WCS->trans->y->coeff[q][r] *= pow(cd1f,q) * pow(cd2f,r);
     730            fprintf(stderr,"%d %d %g\n",q,r,WCS->trans->y->coeff[q][r]);
     731          }
     732        }
     733        fprintf(stderr,"WCSO: %g %g : %g %g : %g %g : %g %g %g %g : %g %g\n",
     734                WCS->crval1,WCS->crval2,
     735                WCS->crpix1,WCS->crpix2,
     736                WCS->cdelt1,WCS->cdelt2,
     737                WCS->trans->x->coeff[1][0],
     738                WCS->trans->x->coeff[0][1],
     739                WCS->trans->y->coeff[1][0],
     740                WCS->trans->y->coeff[0][1],
     741                cd1f,cd2f);
     742
     743        fprintf(stderr,"Size: %d %d %d %d\n",WCS->trans->x->nX,WCS->trans->x->nY,WCS->trans->y->nX,WCS->trans->y->nY);
     744        pmAstromWCStoHeader (hdu->header,WCS);
    616745       
    617746        if (bilevelAstrometry) {
    … …  
    624753        } else {
    625754            // we use a default FPA pixel scale of 1.0
    626             if (!pmAstromReadWCS (input->fpa, chip, hdu->header, 1.0)) {
     755            if (!pmAstromReadWCS (astrom->fpa, chip, hdu->header, 1.0)) {
    627756                psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA.");
    628757                psFree(view);
    … …  
    631760            }
    632761        }
    633 #define CORRECT_INPUT_WCS 1
    634 #if CORRECT_INPUT_WCS
    635         // Correct the input WCS
    636         pmAstromWCS *WCS = pmAstromWCSfromHeader(hdu->header);
    637 
    638         double cd1f = 1.0 * 400;
    639         double cd2f = 1.0 * 400;
    640 
    641         WCS->crpix1 = WCS->crpix1 / cd1f + 0.5;
    642         WCS->crpix2 = WCS->crpix2 / cd2f - 0.5;
    643 
    644         WCS->cdelt1 *= cd1f;
    645         WCS->cdelt2 *= cd2f;
    646 
    647         WCS->trans->x->coeff[1][0] *= cd1f;
    648         WCS->trans->x->coeff[0][1] *= cd2f;
    649         WCS->trans->y->coeff[1][0] *= cd1f;
    650         WCS->trans->y->coeff[0][1] *= cd2f;
    651 
    652 
    653         pmAstromWCStoHeader (hdu->header,WCS);
    654         // End WCS work.
    655 #endif
     762        // Modify structure here.
     763
    656764       
    657765        pmCell *cell;
    … …  
    679787                // Copy the detections from the astrometry carrier to the input, so they can be accessed by
    680788                // pswarpTransformReadout
    681                 pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry
    682                 pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry
    683                 if (detections) {
    684                     psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "Sources from input astrometry", detections);
    685                 }
    686 
     789/*                 pmReadout *astromRO = pmFPAviewThisReadout(view, astrom->fpa); // Readout for astrometry */
     790/*                 pmDetections *detections = psMetadataLookupPtr(&mdok, astromRO->analysis, "PSPHOT.DETECTIONS"); // Sources from astrometry */
     791/*                 if (detections) { */
     792/*                     psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY | PS_META_REPLACE, "Sources from input astrometry", detections); */
     793/*                 } */
     794               
    687795                for (int x = 0; x < readout->image->numCols; x++) {
    688796                  for (int y = 0; y < readout->image->numRows; y++) {
    … …  
    692800                  }
    693801                }
    694                 psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 2);
    695 /*              psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NX",PS_META_REPLACE,"",100); ///< Number of grid points in x */
    696 /*              psMetadataAddS32(config->arguments,PS_LIST_TAIL, "GRID.NY",PS_META_REPLACE,"",100); ///< Number of grid points in y */
    697 
     802               
     803                psMetadataAddS32(config->arguments,PS_LIST_TAIL, "INTERPOLATION.MODE", PS_META_REPLACE, "", 8);
     804
     805                fprintf(stderr,"Transforming Readout!\n");
    698806                pswarpTransformReadout(output, readout, config);
    699807               

  • branches/czw_branch/20120906/pswarp/src/pswarpTransformReadout.c

    r34411 r34442  
    172172        } else {
    173173            pswarpTransformTileArgs *args = job->args->data[0];
    174             // fprintf (stderr, "finished job %d,%d, Nargs: %ld\n", args->gridX, args->gridY, job->args->n);
     174            //      fprintf (stderr, "finished job %d,%d, Nargs: %ld (%d %d %d %d) (%d %d %d %d)\n", args->gridX, args->gridY, job->args->n,
     175            //               xMin,xMax,yMin,yMax,args->xMin,args->xMax,args->yMin,args->yMax);
    175176            goodPixels += args->goodPixels;
    176177            xMin = PS_MIN(args->xMin, xMin);

  • branches/czw_branch/20120906/pswarp/src/pswarpTransformSources.c

    r34411 r34442  
    4141    if (!outDetections) {
    4242        outDetections = pmDetectionsAlloc();
    43         psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY , "Warped sources", outDetections);
     43        psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY , "Warped sources", outDetections);
    4444    }
    4545    psArray *outSources = outDetections->allSources;

  • branches/czw_branch/20120906/pswarp/src/pswarpTransformTile.c

    r29331 r34442  
    8383    // Iterate over the output image pixels (parent frame)
    8484    long goodPixels = 0;                ///< Number of input pixels landing on the output image
     85    psTrace("pswarp",3,"Size: %d %d\n",xMax,yMax);
    8586    for (int y = yMin; y < yMax; y++) {
    8687        for (int x = xMin; x < xMax; x++) {
     88/*        if ((x == 931)||(x == 717)||(x == 924)) { */
     89/* /\*      if ((y == 1348)||(y == 1249)||(y == 1081)) { *\/ */
     90/* /\*        fprintf(stderr,"BEcause: %d %d [%d %d %d %d]\n",x,y,yMin,yMax,xMin,xMax); *\/ */
     91/* /\*      } *\/ */
     92/*        } */
     93            if (((x == 931)&&(y == 1348))||
     94                ((x == 717)&&(y == 1249))||
     95                ((x == 924)&&(y == 1081))) {
     96              psTrace("pswarp",3,"TT: A\n");
     97              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     98
     99            }
    87100
    88101            // Only transform those pixels requested
    … …  
    90103                continue;
    91104            }
     105            if (((x == 931)&&(y == 1348))||
     106                ((x == 717)&&(y == 1249))||
     107                ((x == 924)&&(y == 1081))) {
     108              psTrace("pswarp",3,"TT: B\n");
     109              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     110
     111            }
    92112
    93113            // pswarpMapApply converts the output coordinate (x,y) to the input coordinate.
    … …  
    95115            double xIn, yIn;            // Input pixel coordinates
    96116            pswarpMapApply(&xIn, &yIn, map, x + 0.5, y + 0.5);
     117            if (((x == 931)&&(y == 1348))||
     118                ((x == 717)&&(y == 1249))||
     119                ((x == 924)&&(y == 1081))) {
     120              psTrace("pswarp",3,"TT: Ct\n");
     121              psTrace("pswarp",3,"TT: %d %d %d %d %g %g %d %d\n",xMin,xMax,yMin,yMax,xIn,yIn,inNumCols,inNumRows);
     122
     123            }
     124
    97125            if (xIn < 0 || xIn >= inNumCols || yIn < 0 || yIn >= inNumRows) {
    98126                continue;
    99127            }
     128            if (((x == 931)&&(y == 1348))||
     129                ((x == 717)&&(y == 1249))||
     130                ((x == 924)&&(y == 1081))) {
     131              psTrace("pswarp",3,"TT: C\n");
     132              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     133
     134            }
    100135
    101136            // psImagePixelInterpolate determines the value at pixel coordinate (x,y) in child coordinates
    … …  
    109144            int xOut = x - outCol0, yOut = y - outRow0; ///< Position on output image
    110145
     146            if ((xIn >= args->interp->image->numCols - 1)||
     147                (yIn >= args->interp->image->numRows - 1)) {
     148              imageValue = NAN;
     149            }
    111150            if (outImageData) {
    112151                outImageData[yOut][xOut] = imageValue * jacobian;
    … …  
    120159
    121160            goodPixels++;
     161
     162            if ((x == 931)||(x == 717)||(x == 924)) {
     163              if ((y == 1348)||(y == 1249)||(y == 1081)) {
     164              psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld %g %g %g %g\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian,
     165                      args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp,
     166                      0.0,0.0,0.0,0.0);
     167/*                    inImage->data.F32[(int) yIn][(int) xIn], */
     168/*                    inImage->data.F32[(int) yIn+1][(int) xIn], */
     169/*                    inImage->data.F32[(int) yIn+1][(int) xIn+1], */
     170/*                    inImage->data.F32[(int) yIn][(int) xIn+1]); */
     171              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     172              }
     173            }
     174           
     175            if (((xOut == 931)&&(yOut == 1348))||
     176                ((xOut == 717)&&(yOut == 1249))||
     177                ((xOut == 924)&&(yOut == 1081))) {
     178              psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld %g %g %g %g\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian,
     179                      args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp,
     180                                      0.0,0.0,0.0,0.0);
     181/*                    inImage->data.F32[(int) yIn][(int) xIn], */
     182/*                    inImage->data.F32[(int) yIn+1][(int) xIn], */
     183/*                    inImage->data.F32[(int) yIn+1][(int) xIn+1], */
     184/*                    inImage->data.F32[(int) yIn][(int) xIn+1]); */
     185              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     186
     187            }
     188            if ((xOut ==  443)&&(yOut == 659)) {
     189              psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian,
     190                      args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp);
     191              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     192
     193            }
     194            if ((xOut ==  524)&&(yOut == 576)) {
     195              psTrace("pswarp",3,"TT: %d %d %g %g %d %d %g %g (%g) %ld\n",x,y,xIn,yIn,xOut,yOut,imageValue,jacobian,
     196                      args->input->image->data.F32[(int) yIn][(int) xIn],(long) args->interp);
     197              psTrace("pswarp",3,"TT: %d %d %d %d\n",xMin,xMax,yMin,yMax);
     198             
     199            }
    122200        }
    123201    }
Note: See TracChangeset for help on using the changeset viewer.