IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 20, 2004, 11:38:46 AM (22 years ago)
Author:
gusciora
Message:

Working version.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/pmSubtractSky.c

    r2755 r2760  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2004-12-18 02:27:42 $
     8 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2004-12-20 21:38:46 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2020
    2121/******************************************************************************
    22  *****************************************************************************/
    23 int p_psDetermineNumBits(unsigned int data)
    24 {
    25     int i;
    26     unsigned int tmpData = data;
    27     int numBits = 0;
    28 
    29     for (i=0;i<sizeof(unsigned int);i++) {
     22p_psDetermineNumBits(data): This routine takes an insigned int as an argument
     23and returns the number of non-zero bits.
     24 *****************************************************************************/
     25psS32 p_psDetermineNumBits(psU32 data)
     26{
     27    psS32 i;
     28    psU32 tmpData = data;
     29    psS32 numBits = 0;
     30
     31    for (i=0;i<sizeof(psU32);i++) {
    3032        if (0x0001 && tmpData) {
    3133            numBits++;
    … …  
    3739
    3840/******************************************************************************
     41getHighestPriorityStatOption(statOptions): this routine takes as in put a
     42psStats->options with multiple options set and returns one with a single
     43option set according to the precedence set in the SDRS.
    3944 *****************************************************************************/
    4045psU64 getHighestPriorityStatOption(psU64 statOptions)
    … …  
    8590    psStats *myStats = psStatsAlloc(statOptions);
    8691
    87     for (int row = 0; row < origImage->numRows ; row+=binFactor) {
    88         for (int col = 0; col < origImage->numCols ; col+=binFactor) {
    89             int count = 0;
    90             for (int binRow = 0; binRow <= binFactor ; binRow++) {
    91                 for (int binCol = 0; binCol <= binFactor ; binCol++) {
     92    for (psS32 row = 0; row < origImage->numRows ; row+=binFactor) {
     93        for (psS32 col = 0; col < origImage->numCols ; col+=binFactor) {
     94            psS32 count = 0;
     95            for (psS32 binRow = 0; binRow <= binFactor ; binRow++) {
     96                for (psS32 binCol = 0; binCol <= binFactor ; binCol++) {
    9297                    if (((row + binRow) < origImage->numRows) &&
    9398                            ((col + binCol) < origImage->numCols)) {
    … …  
    123128
    124129/******************************************************************************
     130CalculatePolyTerms(xOrder, yOrder): this routine will calculate the number of
     131coefficients (or terms) in a polynomial of order (xOrder, yOrder).
    125132 *****************************************************************************/
    126133psS32 CalculatePolyTerms(psS32 xOrder, psS32 yOrder)
    … …  
    194201 
    195202    sums[i][j] == x^i * y^j
    196  *****************************************************************************/
    197 // XXX: Use variable size arrays for polynomial sums.
     203 
     204XXX: Use a psImage for the sums data structure?
     205XXX: Check for non-NULL sums argument?
     206XXX: Check for positive x- and yOrder.
     207XXX: Use variable size arrays for polynomial sums.
     208 *****************************************************************************/
    198209#define PS_MAX_POLYNOMIAL_ORDER 20
    199210void buildSums(psF64 x,
    … …  
    221232
    222233/******************************************************************************
     234ImageFitPolynomial(myPoly, binnedImage, maskImage): this private routine takes
     235an input image along with a mask and fits a polynomial to it.  The degree of
     236the polynomial is specified by input parameter myPoly, and need not be
     237symmetrical in orders of X and Y.  The polynomial must be type
     238PS_POLYNOMIAL_ORD.  If there are not enough rows or columns in the input image
     239for the order of the polynomial, then that order is reduced.  The algorithm
     240used in this routine is based on that of the pilot project ADD, but is not
     241documented anywhere.
    223242 *****************************************************************************/
    224243psPolynomial2D *ImageFitPolynomial(psPolynomial2D *myPoly,
    … …  
    263282    psS32 aRow;
    264283    psS32 aCol;
    265     // XXX: Document this.  The myPoly->nX and ->nY terms are actual 1 larger
     284    // XXX: Document this.  The myPoly->nX and ->nY terms are actually 1 larger
    266285    // than the order of the polynomial.
    267286    psS32 **polyTerms = buildPolyTerms(myPoly->nX-1, myPoly->nY-1);
    … …  
    273292    psVector *outPerm = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
    274293
     294    //
     295    // Initialize A matrix and B vector.
     296    //
    275297    for(i=0;i<A->numRows;i++) {
    276298        for(j=0;j<A->numCols;j++) {
    … …  
    282304    }
    283305
     306
     307    //
     308    // We build the A matrix and B vector.
     309    //
    284310    for (x=0;x<binnedImage->numRows;x++) {
    285311        for (y=0;y<binnedImage->numCols;y++) {
    … …  
    315341    for (aRow=0;aRow<localPolyTerms;aRow++) {
    316342        for (aCol=0;aCol<localPolyTerms;aCol++) {
    317             psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 6,
     343            psTrace(".psModule.pmSubtractSky.ImageFitPolynomial", 8,
    318344                    "A[%d][%d] is %f\n", aRow, aCol, A->data.F64[aRow][aCol]);
    319345        }
    … …  
    325351    }
    326352
    327 
     353    //
    328354    // Solve the matrix equations for the polynomial coefficients C.
     355    //
    329356    Aout = psMatrixLUD(Aout, outPerm, A);
    330357    psVector *C = psVectorAlloc(localPolyTerms, PS_TYPE_F64);
    331358    psMatrixLUSolve(C, Aout, B, outPerm);
    332359
     360    //
    333361    // Set the appropriate coefficients in the myPoly structure.
     362    //
    334363    for (i=0;i<localPolyTerms;i++) {
    335364        myPoly->coeff[ polyTerms[i][0] ][ polyTerms[i][1] ] = C->data.F64[i];
    … …  
    338367    }
    339368
     369    //
    340370    // Free data structures that were allocated in this module.
     371    //
    341372    for (i=0;i<localPolyTerms;i++) {
    342373        psFree(polyTerms[i]);
    … …  
    349380    psFree(outPerm);
    350381
     382    //
    351383    // We restore the original size of the polynomial and set remaining
    352384    // coefficients to 0.0.
     385    //
    353386    if (oldPolyX != -1) {
    354387        myPoly->nX = oldPolyX;
    … …  
    382415                         void *fitSpec,
    383416                         psFit fit,
    384                          int binFactor,
     417                         psS32 binFactor,
    385418                         psStats *stats,
    386                          float clipSD)
    387 {
     419                         psF32 clipSD)
     420{
     421    PS_READOUT_CHECK_NULL(in, NULL);
     422    PS_READOUT_CHECK_EMPTY(in, NULL);
     423    PS_READOUT_CHECK_TYPE(in, PS_TYPE_F32, NULL);
    388424    psTrace(".psModule.pmSubtractSky", 4,
    389425            "---- pmSubtractSky() begin ----\n");
    … …  
    396432    psImage *origImage = in->image;
    397433    psImage *binnedImage = NULL;
    398     psPolynomial2D *myPoly;
     434    psPolynomial2D *myPoly = NULL;
    399435    psImage *binnedMaskImage = NULL;
    400 
     436    psU32 oldStatOptions = 0;
     437
     438    //
     439    // Determine which statistic to use when binning pixels, if any.
     440    //
    401441    psStatsOptions statOptions = stats->options;
    402     if (1 > p_psDetermineNumBits(statOptions)) {
     442    if (1 < p_psDetermineNumBits(statOptions)) {
    403443        //XXX  psWarning(PS_ERR_UNKNOWN,true, "Multiple statistical options have been requested.\n");
    404444        statOptions = getHighestPriorityStatOption(statOptions);
    405445        // XXX: Don't modify input parameter.
     446        oldStatOptions = stats->options;
    406447        stats->options = statOptions;
    407         if (statOptions <= 0) {
    408             psLogMsg(__func__, PS_LOG_WARN,
    409                      "WARNING: pmSubtractSky(): unallowable stats->options was requested\n");
    410             return(in);
    411         }
     448    }
     449    if (0 >= p_psDetermineNumBits(statOptions)) {
     450        psLogMsg(__func__, PS_LOG_WARN,
     451                 "WARNING: pmSubtractSky(): no stats->options was requested\n");
     452        return(in);
    412453    }
    413454
    … …  
    430471    //
    431472    // Bin the input image according to input parameters.
    432     //
    433     if ((binFactor <= 1) ||
    434             (stats == NULL) ||
    435             (0 == p_psDetermineNumBits(statOptions))) {
     473    // Create a new binned image mask.
     474    //
     475    if ((binFactor <= 1) || (stats == NULL) || (0 == p_psDetermineNumBits(statOptions))) {
     476        // No binning is required here.  Simply create a copy of the image
     477        // and a mask.
    436478        binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
     479        if (in->mask != NULL) {
     480            binnedMaskImage = psImageCopy(binnedMaskImage, in->mask, PS_TYPE_U8);
     481        } else {
     482            binnedMaskImage = psImageAlloc(binnedImage->numCols,
     483                                           binnedImage->numRows,
     484                                           PS_TYPE_U8);
     485            PS_IMAGE_SET_U8(binnedMaskImage, 0);
     486        }
    437487    } else {
    438         // Add the original image mask in here.
    439         binnedImage = psImageCopy(binnedImage, origImage, PS_TYPE_F32);
    440         binnedImage = psImageRebin(NULL, binnedImage, NULL, 0, binFactor, stats);
     488        binnedImage = psImageRebin(NULL, origImage, in->mask, 0, binFactor, stats);
     489        binnedMaskImage = psImageAlloc(binnedImage->numCols,
     490                                       binnedImage->numRows,
     491                                       PS_TYPE_U8);
     492        PS_IMAGE_SET_U8(binnedMaskImage, 0);
    441493    }
    442494    psTrace(".psModule.pmSubtractSky", 4,
    … …  
    456508        myStats =  psImageStats(myStats, binnedImage, NULL, 0);
    457509        p_psGetStatValue(myStats, &binnedMean);
    458         psTrace(".psModule.pmSubtractSky", 8,
     510        psTrace(".psModule.pmSubtractSky", 6,
    459511                "binned Mean is %f\n", binnedMean);
    460512
    … …  
    463515        p_psGetStatValue(myStats, &binnedStdev);
    464516        psFree(myStats);
    465         psTrace(".psModule.pmSubtractSky", 8,
     517        psTrace(".psModule.pmSubtractSky", 6,
    466518                "binned StDev is %f\n", binnedStdev);
    467519
    468520        // Clip all pixels which are more than clipSD sigmas from the mean.
    469         // XXX: Is this correct?  We simply set the mask.
    470         // XXX: we must unset this later since we modify the image mask.
    471         // XXX: Determine which pixels, mask or image, should be clipped.
    472 
    473         binnedMaskImage = psImageAlloc(binnedImage->numCols,
    474                                        binnedImage->numRows,
    475                                        PS_TYPE_U8);
    476 
    477         psTrace(".psModule.pmSubtractSky", 8,
     521        psTrace(".psModule.pmSubtractSky", 6,
    478522                "clipSD is %f\n", clipSD);
    479523
    480         for (int row = 0; row < binnedImage->numRows ; row++) {
    481             for (int col = 0; col < binnedImage->numCols ; col++) {
     524        for (psS32 row = 0; row < binnedImage->numRows ; row++) {
     525            for (psS32 col = 0; col < binnedImage->numCols ; col++) {
    482526                if (fabs(binnedImage->data.F32[row][col] - binnedMean) >
    483527                        (clipSD * binnedStdev)) {
    … …  
    499543
    500544        if (myPoly != NULL) {
    501             // XXX:Add ordinary polynomials to psImageEvalPolynomial()
    502             for (int row = 0; row < binnedImage->numRows ; row++) {
    503                 for (int col = 0; col < binnedImage->numCols ; col++) {
    504                     binnedImage->data.F32[row][col] =
    505                         psPolynomial2DEval(myPoly, (psF32) row, (psF32) col);
    506                     psTrace(".psModule.pmSubtractSky", 8,
    507                             "binned Image[%d][%d] is %f\n",
    508                             row, col, binnedImage->data.F32[row][col]);
    509                 }
    510             }
     545            // Set the pixels in the binned image to that of the polynomial.
     546            binnedImage = psImageEvalPolynomial(binnedImage, myPoly);
    511547        } else {
    512548            psLogMsg(__func__, PS_LOG_WARN,
    … …  
    516552                psFree(binnedImage);
    517553            }
     554            if (oldStatOptions != 0) {
     555                stats->options = statOptions;
     556            }
    518557            return(in);
    519558        }
    … …  
    527566    if (binFactor <= 1) {
    528567        // The binned image is the same size as the original image.
    529         for (int row = 0; row < origImage->numRows ; row++) {
    530             for (int col = 0; col < origImage->numCols ; col++) {
     568        for (psS32 row = 0; row < origImage->numRows ; row++) {
     569            for (psS32 col = 0; col < origImage->numCols ; col++) {
    531570                origImage->data.F32[row][col]-= binnedImage->data.F32[row][col];
    532571            }
    533572        }
    534573    } else {
    535         for (int row = 0; row < origImage->numRows ; row++) {
    536             for (int col = 0; col < origImage->numCols ; col++) {
     574        for (psS32 row = 0; row < origImage->numRows ; row++) {
     575            for (psS32 col = 0; col < origImage->numCols ; col++) {
    537576                // We calculate the F32 value of the pixel coordinates in the
    538577                // binned image and then use a pixel interpolation routine to
    … …  
    548587
    549588                psF32 binPixel = (psF32) psImagePixelInterpolate(
    550                                      binnedImage, binRowF64, binColF64,
     589                                     binnedImage, binColF64, binRowF64,
    551590                                     NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
    552591                origImage->data.F32[row][col]-= binPixel;
    553592
    554593                psTrace(".psModule.pmSubtractSky", 8,
    555                         "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f (%f)\n",
    556                         row, col, binRowF64-0.5, binColF64-0.5, binPixel,
    557                         binnedImage->data.F32[(psS32)binRowF64][(psS32)binColF64]);
    558             }
    559         }
    560 
    561     }
    562 
     594                        "image[%d][%d] <--> binnedImage[%.2f][%.2f]: %f\n",
     595                        row, col, binRowF64-0.5, binColF64-0.5, binPixel);
     596            }
     597        }
     598
     599    }
    563600    psFree(binnedMaskImage);
    564     if (!((binFactor <= 1) || (stats == NULL))) {
    565         psFree(binnedImage);
     601    psFree(binnedImage);
     602    if (oldStatOptions != 0) {
     603        stats->options = statOptions;
    566604    }
    567605
Note: See TracChangeset for help on using the changeset viewer.