IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37050 for trunk/psLib


Ignore:
Timestamp:
Jul 17, 2014, 10:22:22 AM (12 years ago)
Author:
eugene
Message:

minor changes from EAM testing of pswarp bug; nothing substantive changed

Location:
trunk/psLib/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/fits/psFitsImage.c

    r36537 r37050  
    345345        if (newScaleZero) {
    346346            // Choose an appropriate BSCALE and BZERO
    347           if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) {
    348                 // We can't have the write dying for this reason --- try to save it somehow!
    349                 psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
    350                 psErrorClear();
    351                 return psMemIncrRefCounter((psImage*)image);
    352             }
     347            if (!psFitsScaleDetermine(bscale, bzero, boffset, bsoften, blank, image, mask, maskVal, fits)) {
     348                // We can't have the write dying for this reason --- try to save it somehow!
     349                psWarning("Unable to determine BSCALE and BZERO for image --- refusing to quantise.");
     350                psErrorClear();
     351                return psMemIncrRefCounter((psImage*)image);
     352            }
    353353        } else {
    354             // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
    355             int status = 0;                 // Status of cfitsio
    356             if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
     354            // Don't want to muck with the current BSCALE and BZERO.  Get the current values and use those.
     355            int status = 0;                 // Status of cfitsio
     356            if (fits_read_key_dbl(fits->fd, "BSCALE", bscale, NULL, &status) && status != KEY_NO_EXIST) {
    357357                psFitsError(status, true, "Unable to read header.");
    358358                return NULL;
    … …  
    371371        }
    372372
    373         return psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng);
     373        // XXX this test triggers the pswarp memory corruption problem.  you need to include "psThread.h" above to activate this
     374        // for (int nx = 0; nx < 3; nx++) {
     375        //     psList *pending = psThreadGetPendingQueue();
     376        //     if (pending && pending->head) {
     377        //      fprintf (stderr, "thread jobs are still pending...\n");
     378        //     }
     379        //     // XXX can we mess up the memory?
     380        //     psImage *temp1 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32);
     381        //     fprintf (stderr, "corruption check set 1\n");
     382        //     psMemCheckCorruption (stderr, true);
     383        //     psImageInit (temp1, 2.0);
     384        //     psFree (temp1);
     385        //     fprintf (stderr, "corruption check set 2\n");
     386        //     psMemCheckCorruption (stderr, true);
     387        // }
     388
     389        // fprintf (stderr, "corruption check psFitsInsertImageWithMask\n");
     390        // psMemCheckCorruption (stderr, true);
     391
     392        psImage *out = psFitsScaleForDisk(image, fits, *bscale, *bzero, *boffset, *bsoften, rng);
     393        return out;
    374394    }
    375395
    … …  
    638658    long blank = 0;                     // Blank (undefined) value for image
    639659    psFitsFloat floatType;              // Custom floating-point convention type
    640     psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, image,
    641                                                    mask, maskVal, NULL, true); // Image to write out
     660
     661    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType,
     662                                                   fits, image, mask, maskVal, NULL, true); // Image to write out
     663
    642664    if (!diskImage) {
    643665        psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
    … …  
    666688        goto INSERT_DONE;
    667689    }
     690
    668691    if (cfitsioBzero != 0.0) {
    669692        psAssert(bzero == 0.0 && bscale == 1.0,
    … …  
    840863        }
    841864    }
     865
    842866    if (psFitsError(status, true, "Could not write image to file.")) {
    843867        success = false;
    … …  
    856880 INSERT_DONE:
    857881    psFree(diskImage);
     882
    858883    if (compress) {
    859884        // Restore compression state
    … …  
    906931    psImage *diskImage = imageToDiskRepresentation(&bscale, &bzero, &boffset, &bsoften, &blank, &floatType, fits, input,
    907932                                                   mask, maskVal, NULL, false); // Image to write out
     933   
    908934    if (!diskImage) {
    909935        psError(PS_ERR_UNKNOWN, false, "Unable to convert image to desired disk format.");
    … …  
    932958        goto UPDATE_DONE;
    933959    }
     960
    934961    if (cfitsioBzero != 0.0) {
    935962        psAssert(bzero == 0.0 && bscale == 1.0,
    … …  
    9791006        goto UPDATE_DONE;
    9801007    }
    981 
     1008   
    9821009    // We only want cfitsio to do the scale and zero if the type conversion requires it (e.g., input type is
    9831010    // an unsigned integer type).  In all other cases, we have already converted the image to use the

  • trunk/psLib/src/fits/psFitsScale.c

    r36760 r37050  
    929929
    930930
    931 psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften,
    932                             psRandom *rng)
     931psImage *psFitsScaleForDisk(const psImage *image, const psFits *fits, double bscale, double bzero, double boffset, double bsoften, psRandom *rng)
    933932{
    934933    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
    … …  
    971970    psImage *out = psImageAlloc(numCols, numRows, outType); // Output image
    972971
    973     if (!psMemIncrRefCounter(rng) && options->fuzz) {
    974         // Don't blab about which seed we're going to get --- it's not necessary for this purpose
    975         rng = psRandomAlloc(PS_RANDOM_TAUS);
    976     }
    977 
    978 #define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE) \
    979     case PS_TYPE_##OUTTYPE: { \
    980         double scale = 1.0 / bscale; \
    981         double zero = bzero; \
     972    // fprintf (stderr, "corruption check psFitsScaleForDisk 1\n");
     973    // psMemCheckCorruption (stderr, true);
     974
     975    if (options->fuzz && !psMemIncrRefCounter (rng)) {
     976      // Don't blab about which seed we're going to get --- it's not necessary for this purpose
     977      rng = psRandomAlloc(PS_RANDOM_TAUS);
     978    }
     979
     980#define SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, OUTTYPE)                  \
     981    case PS_TYPE_##OUTTYPE: {                                           \
     982        double scale = 1.0 / bscale;                                    \
     983        double zero = bzero;                                            \
    982984        /* Note: BITPIX=8 treated differently, since it uses unsigned values; the rest use signed */ \
    983         double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1); \
     985        double min = bitpix == 8 ? 0 : -pow(2.0, options->bitpix - 1);  \
    984986        double max = bitpix == 8 ? 255 : (pow(2.0, options->bitpix - 1) - 1.0); \
    985         for (int y = 0; y < numRows; y++) { \
    986             for (int x = 0; x < numCols; x++) { \
    987               ps##INTYPE value; \
    988               if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)||       \
    989                 (options->scaling == PS_FITS_SCALE_LOG_MANUAL)||        \
    990                   (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \
    991                   (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
    992                   (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
    993                 if (isfinite((IN)->data.INTYPE[y][x])) {                \
    994                   value = log10( (IN)->data.INTYPE[y][x] - boffset );   \
     987        for (int y = 0; y < numRows; y++) {                             \
     988            for (int x = 0; x < numCols; x++) {                         \
     989                ps##INTYPE value;                                       \
     990                if ((options->scaling == PS_FITS_SCALE_LOG_RANGE)||     \
     991                    (options->scaling == PS_FITS_SCALE_LOG_MANUAL)||    \
     992                    (options->scaling == PS_FITS_SCALE_LOG_STDEV_POSITIVE)|| \
     993                    (options->scaling == PS_FITS_SCALE_LOG_STDEV_NEGATIVE)|| \
     994                    (options->scaling == PS_FITS_SCALE_LOG_STDEV_BOTH)) { \
     995                    if (isfinite((IN)->data.INTYPE[y][x])) {            \
     996                        value = log10( (IN)->data.INTYPE[y][x] - boffset ); \
     997                    }                                                   \
     998                    else {                                              \
     999                        value = (IN)->data.INTYPE[y][x];                \
     1000                    }                                                   \
     1001                }                                                       \
     1002                else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)|| \
     1003                         (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)|| \
     1004                         (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \
     1005                         (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \
     1006                         (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) { \
     1007                    if (isfinite((IN)->data.INTYPE[y][x])) {            \
     1008                        value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \
     1009                    }                                                   \
     1010                    else {                                              \
     1011                        value = (IN)->data.INTYPE[y][x];                \
     1012                    }                                                   \
    9951013                }                                                       \
    9961014                else {                                                  \
    997                   value = (IN)->data.INTYPE[y][x];                      \
     1015                    value = (IN)->data.INTYPE[y][x];                    \
    9981016                }                                                       \
    999               }                                                         \
    1000               else if ((options->scaling == PS_FITS_SCALE_ASINH_RANGE)||        \
    1001                 (options->scaling == PS_FITS_SCALE_ASINH_MANUAL)||        \
    1002                   (options->scaling == PS_FITS_SCALE_ASINH_STDEV_POSITIVE)|| \
    1003                   (options->scaling == PS_FITS_SCALE_ASINH_STDEV_NEGATIVE)|| \
    1004                   (options->scaling == PS_FITS_SCALE_ASINH_STDEV_BOTH)) {       \
    1005                 if (isfinite((IN)->data.INTYPE[y][x])) {                \
    1006                   value = 1.0857362 * (asinh( ((IN)->data.INTYPE[y][x] - boffset) / (2.0 * bsoften))); \
     1017                if (!isfinite(value)) {                                 \
     1018                    /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \
     1019                    /* something, and "min" would produce holes in the cores of bright stars. */ \
     1020                    (OUT)->data.OUTTYPE[y][x] = max;                    \
     1021                } else {                                                \
     1022                    value = (value - zero) * scale;                     \
     1023                    if (options->fuzz && (value - (int)value != 0.0)) { \
     1024                        /* Add random factor [0.0,1.0): adds a variance of 1/12, */ \
     1025                        /* but preserves the expectation value given the floor() */ \
     1026                        value += psRandomUniform(rng) ;                 \
     1027                    }                                                   \
     1028                    /* Check for underflow and overflow; set either to max */ \
     1029                    (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : floor(value)); \
    10071030                }                                                       \
    1008                 else {                                                  \
    1009                   value = (IN)->data.INTYPE[y][x];                      \
    1010                 }                                                       \
    1011               }                                                         \
    1012                 else { \
    1013                   value = (IN)->data.INTYPE[y][x];      \
    1014                 }                                           \
    1015                 if (!isfinite(value)) { \
    1016                     /* This choice of "max" for non-finite pixels is mainly cosmetic --- it has to be */ \
    1017                     /* something, and "min" would produce holes in the cores of bright stars. */ \
    1018                     (OUT)->data.OUTTYPE[y][x] = max; \
    1019                 } else { \
    1020                     value = (value - zero) * scale; \
    1021                     if (options->fuzz && (value - (int)value != 0.0)) { \
    1022                        /* Add random factor [-0.5,0.5): adds a variance of 1/12, */ \
    1023                        /* but preserves the expectation value */ \
    1024                         value += psRandomUniform(rng) ; \
    1025                     } \
    1026                     /* Check for underflow and overflow; set either to max */ \
    1027                     (OUT)->data.OUTTYPE[y][x] = (value < min || value > max ? max : floor(value)); \
    1028                 } \
    1029             } \
    1030         } \
    1031         break; \
    1032     }
    1033 
    1034 #define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT) \
    1035     case PS_TYPE_##INTYPE: { \
    1036         switch (outType) { \
    1037           SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);;   \
    1038           SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);;  \
    1039           SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);;  \
    1040           SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);;  \
    1041           default: \
    1042             psAbort("Should be unreachable."); \
    1043         } \
    1044         break; \
    1045     }
    1046 
     1031            }                                                           \
     1032        }                                                               \
     1033        break;                                                          \
     1034    }
     1035   
     1036#define SCALE_WRITE_IN_CASE(IN, INTYPE, OUT)                    \
     1037    case PS_TYPE_##INTYPE: {                                    \
     1038        switch (outType) {                                      \
     1039            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, U8);;         \
     1040            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S16);;        \
     1041            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S32);;        \
     1042            SCALE_WRITE_OUT_CASE(IN, INTYPE, OUT, S64);;        \
     1043          default:                                              \
     1044            psAbort("Should be unreachable.");                  \
     1045        }                                                       \
     1046        break;                                                  \
     1047    }
     1048   
    10471049    switch (image->type.type) {
    1048         SCALE_WRITE_IN_CASE(image, F32, out);
    1049         SCALE_WRITE_IN_CASE(image, F64, out);
     1050        SCALE_WRITE_IN_CASE(image, F32, out);
     1051        SCALE_WRITE_IN_CASE(image, F64, out);
    10501052      default:
    10511053        psAbort("Should be unreachable.");
    … …  
    10531055
    10541056    psFree(rng);
    1055 
    10561057    return out;
    10571058}

  • trunk/psLib/src/imageops/psImageBackground.c

    r31152 r37050  
    124124        }
    125125    } else if (nSubset >= nGoodPixels) {
    126             psTrace ("psLib.imageops", 4, "case 2: nSubset >= nGoodPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
     126        psTrace ("psLib.imageops", 4, "case 2: nSubset >= nGoodPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels);
    127127        // in this case, we have to select from all masked pixels just to get the desired
    128128        // sample size
    … …  
    157157            double frnd = psRandomUniform(rng);
    158158            int pixel = pixelVector->n * frnd;
    159                
     159            // psAssert (pixel < pixelVector->nalloc, "oops, we went too far (1)");
     160            // psAssert (pixel >= 0, "oops, we went too far (2)");
     161            // psAssert (pixelVector->n - 1 < pixelVector->nalloc, "oops, we went too far (3)");
     162            // psAssert (pixelVector->n - 1 >= 0, "oops, we went too far (4)");
     163
    160164            psVectorAppend(values, pixelVector->data.F32[pixel]);
    161165            pixelVector->data.F32[pixel] = pixelVector->data.F32[pixelVector->n - 1];

  • trunk/psLib/src/sys/psMemory.c

    r32746 r37050  
    13151315    checkingForCorruption = true;
    13161316
     1317    // int nPrint = 0;
    13171318    psS32 nbad = 0;               // number of bad blocks
    13181319    for (psMemBlock *memBlock = (psMemBlock *) lastMemBlockAllocated; memBlock != NULL; memBlock = memBlock->nextBlock) {
     1320        // if (nPrint < 20) {
     1321        //     psMemBlockPrint (stderr, memBlock);
     1322        //     nPrint ++;
     1323        // }
    13191324        if (isBadMemBlock(output, memBlock, __FILE__, __LINE__, __func__)) {
    13201325            nbad++;

  • trunk/psLib/src/sys/psThread.c

    r34153 r37050  
    152152    psThreadJob *job = psListGetAndRemove((psList *) done, PS_LIST_HEAD);
    153153    return job;
     154}
     155
     156psList *psThreadGetPendingQueue () {
     157  return (psList *) pending;
    154158}
    155159

  • trunk/psLib/src/sys/psThread.h

    r33089 r37050  
    9090psThreadJob *psThreadJobGetDone(void);
    9191
     92// utility function to access the pending queue
     93psList *psThreadGetPendingQueue ();
     94
    9295/// Allocate a thread task
    9396psThreadTask *psThreadTaskAlloc(const char *type, // Type of task
Note: See TracChangeset for help on using the changeset viewer.