IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 15, 2013, 5:56:56 PM (13 years ago)
Author:
eugene
Message:

more work on the central pixel optimizations -- perhaps not needed (not so expensive?); add some interactive support for PCM chisq fitting; EXP and DEV are for the moment using subdivided central pixels, but this is perhaps too slow?; turn on sky fitting for the PCM model fitting

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c

    r35948 r35961  
    209209
    210210// XXX for test purposes only:
    211 # define TEST_IMAGE 1
     211# define TEST_IMAGE 0
    212212# if (TEST_IMAGE)
    213213static psImage *map = NULL;
    214214# endif
    215215
     216float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta);
     217
    216218float pmModelCP_GetFlux (pmModelCP *cp, float dx, float dy, float theta) {
    217219
    … …  
    221223# endif
    222224
    223     float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta);
    224    
     225    // float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta);
     226    // float flux = pmModelCP_GetFlux_Old (cp, dx, dy, theta);
     227    float flux = pmModelCP_GetFlux_RotSquare (cp, dx, dy, theta);
     228   
     229    // RotSquare for theta = 0.0 & Bresen give the same answer
     230    // if I count from x[0] <= ix < x[1]
     231
    225232# if (TEST_IMAGE)
    226233    psFits *fits = psFitsOpen ("map.fits", "w");
    … …  
    309316}
    310317
    311 float pmModelCP_GetFlux_RotSqaure (pmModelCP *cp, float dx, float dy, float theta) {
     318// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
     319# define COMPARE_INDEX(A,B) (y[A] < y[B])
     320# define SWAP_INDEX(TYPE,A,B) {                         \
     321        int tmp;                                        \
     322        if (A != B) {                                   \
     323            tmp = x[A];                                 \
     324            x[A] = x[B];                                \
     325            x[B] = tmp;                                 \
     326            tmp = y[A];                                 \
     327            y[A] = y[B];                                \
     328            y[B] = tmp;                                 \
     329        }                                               \
     330    }
     331
     332bool pmModelCP_SortCorners (int *x, int *y, int Npar) {
     333
     334    if (Npar < 2) return true;
     335
     336    // sort the vector set by the radius
     337    PSSORT (Npar, COMPARE_INDEX, SWAP_INDEX, NONE);
     338    return true;
     339}
     340
     341float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta) {
    312342
    313343    // the cp data is defined for the central 3x3 pixels.  we allow dx,dy to have values of
    … …  
    324354    float cs = cos(theta*PS_RAD_DEG);
    325355    float sn = sin(theta*PS_RAD_DEG);
    326 
    327356    float Nsub = 11.0;
    328     int Xsub00 = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
    329     int Ysub00 = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
    330     int Xsub01 = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
    331     int Ysub01 = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
    332     int Xsub10 = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
    333     int Ysub10 = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
    334     int Xsub11 = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
    335     int Ysub11 = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
    336 
    337     /* generic rotated square:
    338        
    339      */
    340 
    341     int Xmin, Xmax, Ymin, Ymax;
    342 
    343     Xmin = PS_MIN(Xsub00,Xsub01);
    344     Xmin = PS_MIN(Xsub10,Xmin);
    345     Xmin = PS_MIN(Xsub11,Xmin);
    346     Xmin = PS_MIN(Xmin, cp->flux->numCols - 1);
    347     Xmin = PS_MAX(Xmin, 0);
    348     Xmax = PS_MAX(Xsub00,Xsub01);
    349     Xmax = PS_MAX(Xsub10,Xmax);
    350     Xmax = PS_MAX(Xsub11,Xmax);
    351     Xmax = PS_MIN(Xmax, cp->flux->numCols - 1);
    352     Xmax = PS_MAX(Xmax, 0);
    353     Ymin = PS_MIN(Ysub00,Ysub01);
    354     Ymin = PS_MIN(Ysub10,Ymin);
    355     Ymin = PS_MIN(Ysub11,Ymin);
    356     Ymin = PS_MIN(Ymin, cp->flux->numRows - 1);
    357     Ymin = PS_MAX(Ymin, 0);
    358     Ymax = PS_MAX(Ysub00,Ysub01);
    359     Ymax = PS_MAX(Ysub10,Ymax);
    360     Ymax = PS_MAX(Ysub11,Ymax);
    361     Ymax = PS_MIN(Ymax, cp->flux->numRows - 1);
    362     Ymax = PS_MAX(Ymax, 0);
    363 
    364     // integrate pixels from Xmin,Ymin to Xmax,Ymax, only include pixels contained in the
    365     // target pixel
     357
     358    int Xsub[4], Ysub[4];
     359
     360    Xsub[0] = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
     361    Ysub[0] = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
     362    Xsub[1] = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
     363    Ysub[1] = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
     364    Xsub[2] = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub;
     365    Ysub[2] = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub;
     366    Xsub[3] = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub;
     367    Ysub[3] = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub;
     368
     369    // first, sort the corners in order of the Y coordinate:
     370    pmModelCP_SortCorners (Xsub, Ysub, 4);
    366371
    367372    float flux = 0.0;
    368     int   npix = 0;
    369     for (int i = Xmin; i < Xmax; i++) {
    370         float dX = i / Nsub - 1.5;
    371         for (int j = Ymin; j < Ymax; j++) {
    372             float dY = j / Nsub - 1.5;
    373 
    374             float Xim =  dX*cs + dY*sn;
    375             if (Xim < (dx - 0.5)) continue;
    376             if (Xim > (dx + 0.5)) continue;
    377 
    378             float Yim = -dX*sn + dY*cs;
    379             if (Yim < (dy - 0.5)) continue;
    380             if (Yim > (dy + 0.5)) continue;
    381 
    382             flux += cp->flux->data.F32[j][i];
    383             npix ++;
    384         }
    385     }
    386            
    387     float normFlux = flux / npix;
    388     return normFlux;
     373    float npix = 0.0;
     374
     375    // if (Ysub[0] == Ysub[1]), we have a simple square
     376    if (Ysub[0] == Ysub[1]) {
     377        psAssert (Ysub[2] == Ysub[3], "not square?");
     378        int Xmin = PS_MIN(Xsub[0], Xsub[1]);
     379        int Xmax = PS_MAX(Xsub[0], Xsub[1]);
     380        for (int iy = Ysub[0]; iy < Ysub[3]; iy++) {
     381            for (int ix = Xmin; ix < Xmax; ix++) {
     382                flux += cp->flux->data.F32[iy][ix];
     383                npix += 1.0;
     384# if (TEST_IMAGE)
     385                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     386                map->data.S32[iy][ix] ++;
     387# endif
     388            }
     389        }
     390        float normFlux = flux / npix;
     391        return normFlux;
     392    }
     393   
     394    // second case: Xsub[1] > Xsub[2]:
     395    if (Xsub[1] > Xsub[2]) {
     396        float dYdXp, dYdXm;
     397        // first segment, Ysub[0] to Ysub[1]:
     398        dYdXp = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]);
     399        dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
     400        for (int iy = Ysub[0]; iy < Ysub[1]; iy++) {
     401            int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0];
     402            int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0];
     403            for (int ix = Xs; ix < Xe; ix ++) {
     404                flux += cp->flux->data.F32[iy][ix];
     405                npix += 1.0;
     406# if (TEST_IMAGE)
     407                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     408                map->data.S32[iy][ix] ++;
     409# endif
     410            }
     411        }
     412        // 2nd segment, Ysub[1] to Ysub[2]:
     413        dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
     414        dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
     415        for (int iy = Ysub[1]; iy < Ysub[2]; iy++) {
     416            int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0];
     417            int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1];
     418            for (int ix = Xs; ix < Xe; ix ++) {
     419                flux += cp->flux->data.F32[iy][ix];
     420                npix += 1.0;
     421# if (TEST_IMAGE)
     422                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     423                map->data.S32[iy][ix] ++;
     424# endif
     425            }
     426        }
     427        // first segment, Ysub[0] to Ysub[1]:
     428        dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
     429        dYdXm = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]);
     430        for (int iy = Ysub[2]; iy < Ysub[3]; iy++) {
     431            int Xs = (iy - Ysub[2]) / dYdXm + Xsub[2];
     432            int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1];
     433            for (int ix = Xs; ix < Xe; ix ++) {
     434                flux += cp->flux->data.F32[iy][ix];
     435                npix += 1.0;
     436# if (TEST_IMAGE)
     437                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     438                map->data.S32[iy][ix] ++;
     439# endif
     440            }
     441        }
     442        float normFlux = flux / npix;
     443        return normFlux;
     444    }
     445
     446    // third case: Xsub[1] < Xsub[2]:
     447    if (Xsub[2] > Xsub[1]) {
     448        // first segment, Ysub[0] to Ysub[1]:
     449        float dYdXp, dYdXm;
     450        dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
     451        dYdXm = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]);
     452        for (int iy = Ysub[0]; iy < Ysub[1]; iy++) {
     453            int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0];
     454            int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0];
     455            for (int ix = Xs; ix < Xe; ix ++) {
     456                flux += cp->flux->data.F32[iy][ix];
     457                npix += 1.0;
     458# if (TEST_IMAGE)
     459                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     460                map->data.S32[iy][ix] ++;
     461# endif
     462            }
     463        }
     464        // 2nd segment, Ysub[1] to Ysub[2]:
     465        dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]);
     466        dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
     467        for (int iy = Ysub[1]; iy < Ysub[2]; iy++) {
     468            int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1];
     469            int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0];
     470            for (int ix = Xs; ix < Xe; ix ++) {
     471                flux += cp->flux->data.F32[iy][ix];
     472                npix += 1.0;
     473# if (TEST_IMAGE)
     474                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     475                map->data.S32[iy][ix] ++;
     476# endif
     477            }
     478        }
     479        // first segment, Ysub[0] to Ysub[1]:
     480        dYdXp = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]);
     481        dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]);
     482        for (int iy = Ysub[2]; iy < Ysub[3]; iy++) {
     483            int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1];
     484            int Xe = (iy - Ysub[2]) / dYdXp + Xsub[2];
     485            for (int ix = Xs; ix < Xe; ix ++) {
     486                flux += cp->flux->data.F32[iy][ix];
     487                npix += 1.0;
     488# if (TEST_IMAGE)
     489                fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]);
     490                map->data.S32[iy][ix] ++;
     491# endif
     492            }
     493        }
     494        float normFlux = flux / npix;
     495        return normFlux;
     496    }
     497    myAbort ("impossible case?");
    389498}
    390499
    … …  
    478587    }
    479588    float normFlux = flux / npix;
    480     fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux);
     589    // fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux);
    481590    return normFlux;
    482591}
    … …  
    574683    }
    575684    float normFlux = flux / npix;
    576     fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux);
     685    // fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux);
    577686    return normFlux;
    578687}
Note: See TracChangeset for help on using the changeset viewer.