IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5745


Ignore:
Timestamp:
Dec 7, 2005, 4:04:22 PM (20 years ago)
Author:
Paul Price
Message:

Using autoconf; some bug fixes to the polynomial order

Location:
trunk/stac
Files:
4 added
1 deleted
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/stac.c

    r5743 r5745  
    77int main(int argc, char **argv)
    88{
    9     double startTime = getTime();
     9    double startTime;
     10    startTime = getTime();
    1011
    1112#if 0
    … …  
    155156        for (int i = 0; i < transformed->n; i++) {
    156157            char shiftName[MAXCHAR];    // Filename of shifted image
    157             sprintf(shiftName, "%s.shift", config->inputs->data[i]);
     158            sprintf(shiftName, "%s.shift", (const char*)config->inputs->data[i]);
    158159            psFits *shiftFile = psFitsOpen(shiftName, "w");
    159160            image = psImageCopy(NULL, transformed->data[i], PS_TYPE_F32);

  • trunk/stac/src/stacCheckMemory.c

    r5743 r5745  
    44
    55
    6 #define LEAKS "leaks.dat"               // File to which to write leaks data
     6#define LEAKS "leaks.dat"               // File to which to write leaks data
    77
    88
    … …  
    1010{
    1111    psLogMsg("stac.memoryPrint", PS_LOG_INFO,
    12              "Memory block %d:\n"
    13              "\tFile %s, line %d, size %d\n"
    14              "\tPosts: %x %x %x\n",
    15              ptr->id, ptr->file, ptr->lineno, ptr->userMemorySize, ptr->startblock, ptr->endblock,
    16              *(void**)((int8_t *)(ptr + 1) + ptr->userMemorySize));
     12             "Memory block %d:\n"
     13             "\tFile %s, line %d, size %d\n"
     14             "\tPosts: %x %x %x\n",
     15             ptr->id, ptr->file, ptr->lineno, ptr->userMemorySize, ptr->startblock, ptr->endblock,
     16             *(void**)((int8_t *)(ptr + 1) + ptr->userMemorySize));
    1717    return 0;
    1818}
    … …  
    2020
    2121void stacMemoryProblem(const psMemBlock* ptr, ///< the pointer to the problematic memory block.
    22                        const char *file, ///< the file in which the problem originated
    23                        psS32 lineno     ///< the line number in which the problem originated
     22                       const char *file, ///< the file in which the problem originated
     23                       psS32 lineno     ///< the line number in which the problem originated
    2424    )
    2525{
    2626    psLogMsg("stac.checkMemory.corruption", PS_LOG_WARN,
    27              "Memory corruption detected in memBlock %d\n"
    28              "\tFile %s, line %d, size %d\n"
    29              "\tPosts: %x %x %x\n",
    30              ptr->id, file, lineno, ptr->userMemorySize, ptr->startblock, ptr->endblock,
    31              (ptr + 1 + ptr->userMemorySize));
     27             "Memory corruption detected in memBlock %d\n"
     28             "\tFile %s, line %d, size %d\n"
     29             "\tPosts: %x %x %x\n",
     30             ptr->id, file, lineno, ptr->userMemorySize, ptr->startblock, ptr->endblock,
     31             (ptr + 1 + ptr->userMemorySize));
    3232}
    3333
    … …  
    3636void stacCheckMemory(void)
    3737{
    38     psMemBlock **leaks = NULL;          // List of leaks
    39     FILE *leakFile;                     // File to write leaks to
     38    psMemBlock **leaks = NULL;          // List of leaks
     39    FILE *leakFile;                     // File to write leaks to
    4040
    4141    psTrace("stac.checkMemory", 1, "Checking for memory problems....\n");
    … …  
    4444
    4545    if ((leakFile = fopen(LEAKS, "w")) == NULL) {
    46         fprintf(stderr, "Unable to open leaks file, %s\n", LEAKS);
    47         return;
     46        fprintf(stderr, "Unable to open leaks file, %s\n", LEAKS);
     47        return;
    4848    }
    4949
    … …  
    5151    psTrace("stac.checkMemory", 1, "%d leaks found.\n", nLeaks);
    5252    for (int i = 0; i < nLeaks; i++) {
    53         psLogMsg("stac.checkMemory.leaks", PS_LOG_WARN,
    54                 "Memory leak detection: memBlock %d\n"
    55                 "\tFile %s, line %d, size %d\n",
    56                 leaks[i]->id, leaks[i]->file, leaks[i]->lineno, leaks[i]->userMemorySize);
     53        psLogMsg("stac.checkMemory.leaks", PS_LOG_WARN,
     54                "Memory leak detection: memBlock %d\n"
     55                "\tFile %s, line %d, size %d\n",
     56                leaks[i]->id, leaks[i]->file, leaks[i]->lineno, leaks[i]->userMemorySize);
    5757    }
    5858
    59     int nCorrupted = psMemCheckCorruption(false); // Number of corrupted
     59    int nCorrupted;                     // Number of corrupted memory blocks
     60    nCorrupted = psMemCheckCorruption(false);
    6061    psTrace("stac.checkMemory", 1, "%d memory blocks corrupted.\n", nCorrupted);
    6162

  • trunk/stac/src/stacCombine.c

    r5743 r5745  
    88#define ABS(x) ((x) >= 0 ? (x) : -(x))
    99
    10 float stacCombineMean(psVector *values, // Values for which to take the mean
    11                       psVector *errors, // Errors in the values
    12                       psVector *masks   // Masks for the values, 1 = don't use, 0 = use
     10float stacCombineMean(psVector *values, // Values for which to take the mean
     11                      psVector *errors, // Errors in the values
     12                      psVector *masks   // Masks for the values, 1 = don't use, 0 = use
    1313    )
    1414{
    … …  
    2626    int num = values->n;
    2727    for (int i = 0; i < num; i++) {
    28         if (! masks->data.U8[i]) {
    29             // "error" here is the variance
    30             sum += values->data.F32[i] / errors->data.F32[i];
    31             weights += 1.0 / errors->data.F32[i];
    32         }
     28        if (! masks->data.U8[i]) {
     29            // "error" here is the variance
     30            sum += values->data.F32[i] / errors->data.F32[i];
     31            weights += 1.0 / errors->data.F32[i];
     32        }
    3333    }
    3434    if (weights > 0.0) {
    35         return (float)(sum/weights);
     35        return (float)(sum/weights);
    3636    } else {
    37         return NAN;
     37        return NAN;
    3838    }
    3939#endif
    … …  
    4343
    4444float stacCombineMedian(psVector *values, // Values for which to take the median
    45                         psVector *errors, // Errors in the values
    46                         psVector *masks // Masks for the values, 0 = don't use, 1 = use
     45                        psVector *errors, // Errors in the values
     46                        psVector *masks // Masks for the values, 0 = don't use, 1 = use
    4747    )
    4848{
    … …  
    5656
    5757
    58 bool stacCombine(psImage **combined,    // The combined image for output
    59                  psArray **rejected,    // Array of rejection masks
    60                  psArray *images,       // Array of transformed images
    61                  psArray *errors,       // Array of transformed error images
    62                  int nReject,           // Number of rejection iterations
    63                  psImage *region,       // Region to combine
    64                  psVector *saturated,   // Saturation limits for each image
    65                  psVector *bad,         // Bad pixel limits for each image
    66                  float reject           // Rejection (k-sigma)
     58bool stacCombine(psImage **combined,    // The combined image for output
     59                 psArray **rejected,    // Array of rejection masks
     60                 psArray *images,       // Array of transformed images
     61                 psArray *errors,       // Array of transformed error images
     62                 int nReject,           // Number of rejection iterations
     63                 psImage *region,       // Region to combine
     64                 psVector *saturated,   // Saturation limits for each image
     65                 psVector *bad,         // Bad pixel limits for each image
     66                 float reject           // Rejection (k-sigma)
    6767    )
    6868{
    … …  
    7676    assert(bad->n == images->n);
    7777
    78     int nImages = images->n;            // Number of images
    79     int numRows = ((psImage*)images->data[0])->numRows; // Image size
     78    int nImages = images->n;            // Number of images
     79    int numRows = ((psImage*)images->data[0])->numRows; // Image size
    8080    int numCols = ((psImage*)images->data[0])->numCols; // Image size
    81    
     81
    8282    // Check dimensions for consistency
    8383    for (int i = 0; i < nImages; i++) {
    84         psImage *image = (psImage *)images->data[i]; // The image
    85         psImage *error = (psImage *)errors->data[i]; // The error image
    86 
    87         assert(image->numCols == numCols && image->numRows == numRows);
    88         assert(error->numCols == numCols && error->numRows == numRows);
     84        psImage *image = (psImage *)images->data[i]; // The image
     85        psImage *error = (psImage *)errors->data[i]; // The error image
     86
     87        assert(image->numCols == numCols && image->numRows == numRows);
     88        assert(error->numCols == numCols && error->numRows == numRows);
    8989    }
    9090
    9191    // Check combined image
    92     assert(!*combined || ((*combined)->numRows == numRows) && ((*combined)->numCols == numCols));
     92    assert(!*combined || ((*combined)->numRows == numRows && (*combined)->numCols == numCols));
    9393    if (*combined == NULL) {
    94         *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
     94        *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
    9595    }
    9696
    9797    // Check area of interest
    98     assert(!region || (region->numRows == numRows) && (region->numCols == numCols));
     98    assert(!region || (region->numRows == numRows && region->numCols == numCols));
    9999
    100100    psTrace("stac.combine", 1, "Combining images....\n");
    … …  
    106106    // Set up rejection masks
    107107    if (nReject > 0) {
    108         if (*rejected == NULL) {
    109             // Allocate the rejection masks, if required
    110             *rejected = psArrayAlloc(nImages);
    111         } else {
    112             assert((*rejected)->n != nImages);
    113         }
    114 
    115         // Create and initialise rejection masks
    116         for (int i = 0; i < nImages; i++) {
    117             (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
    118             for (int r = 0; r < numRows; r++) {
    119                 for (int c = 0; c < numCols; c++) {
    120                     ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
    121                 }
    122             }
    123         }
     108        if (*rejected == NULL) {
     109            // Allocate the rejection masks, if required
     110            *rejected = psArrayAlloc(nImages);
     111        } else {
     112            assert((*rejected)->n != nImages);
     113        }
     114
     115        // Create and initialise rejection masks
     116        for (int i = 0; i < nImages; i++) {
     117            (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
     118            for (int r = 0; r < numRows; r++) {
     119                for (int c = 0; c < numCols; c++) {
     120                    ((psImage*)((*rejected)->data[i]))->data.U8[r][c] = 0;
     121                }
     122            }
     123        }
    124124    }
    125125
    … …  
    127127    // chi^2 image
    128128    psImage *chi2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    129     static int iteration = 0;           // Number of times function has been called
    130 #endif
    131    
     129    static int iteration = 0;           // Number of times function has been called
     130#endif
     131
    132132    for (int y = 0; y < numRows; y++) {
    133         for (int x = 0; x < numCols; x++) {
    134 
    135             // Only combine those pixels requested
    136             if (!region || (region && region->data.U8[y][x])) {
    137            
    138                 // Export pixels into the vector and get stats
    139                 for (int i = 0; i < nImages; i++) {
    140                     float pixel = ((psImage*)images->data[i])->data.F32[y][x];
    141                     float delta = ((psImage*)errors->data[i])->data.F32[y][x]; // This is the variance!
    142                     pixels->data.F32[i] = pixel;
    143                     deltas->data.F32[i] = delta;
    144                     if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
    145                         (! isfinite(pixel)) || (! isfinite(delta))) {
    146                         mask->data.U8[i] = (psU8)1; // Don't use!
    147                     } else {
    148                         mask->data.U8[i] = (psU8)0; // Use.
    149                     }
    150                 }
    151                
    152                 float average = stacCombineMean(pixels, deltas, mask); // Combined value
    153 
    154                 // We set the value BEFORE the rejection iteration because not all pixels that we reject
    155                 // here will be rejected in the final cut.
    156                 (*combined)->data.F32[y][x] = average;
     133        for (int x = 0; x < numCols; x++) {
     134
     135            // Only combine those pixels requested
     136            if (!region || (region && region->data.U8[y][x])) {
     137
     138                // Export pixels into the vector and get stats
     139                for (int i = 0; i < nImages; i++) {
     140                    float pixel = ((psImage*)images->data[i])->data.F32[y][x];
     141                    float delta = ((psImage*)errors->data[i])->data.F32[y][x]; // This is the variance!
     142                    pixels->data.F32[i] = pixel;
     143                    deltas->data.F32[i] = delta;
     144                    if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
     145                        (! isfinite(pixel)) || (! isfinite(delta))) {
     146                        mask->data.U8[i] = (psU8)1; // Don't use!
     147                    } else {
     148                        mask->data.U8[i] = (psU8)0; // Use.
     149                    }
     150                }
     151
     152                float average = stacCombineMean(pixels, deltas, mask); // Combined value
     153
     154                // We set the value BEFORE the rejection iteration because not all pixels that we reject
     155                // here will be rejected in the final cut.
     156                (*combined)->data.F32[y][x] = average;
    157157
    158158#ifdef TESTING
    159                 // Calculate chi^2
    160                 chi2->data.F32[y][x] = 0.0;
    161                 int numGoodPix = 0;
    162                 for (int i = 0; i < nImages; i++) {
    163                     if (! mask->data.U8[i]) {
    164                         chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
    165                         numGoodPix++;
    166                     }
    167                 }
    168                 chi2->data.F32[y][x] /= (float)numGoodPix;
    169 #endif
    170                
    171                 // Rejection iterations
    172                 bool keepGoing = true;  // Keep going with rejection?
    173                 for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
    174                     float max = 0.0;    // Maximum deviation
    175                     int maxIndex = -1;  // Index of the maximum deviation
    176                     for (int i = 0; i < nImages; i++) {
    177                         if (!mask->data.U8[i] &&
    178                             ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
    179                             max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
    180                             maxIndex = i;
    181                         }
    182                     }
    183                     // Reject the pixel with the maximum deviation
    184                     if (max > reject) {
    185                         mask->data.U8[maxIndex] = 1;
    186                         ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
    187                         // Re-do combination following rejection
    188                         average = stacCombineMean(pixels, deltas, mask);
    189                     } else {
    190                         keepGoing = false;
    191                     }
    192                 }
    193            
    194             } // Pixels of interest
    195 
    196         }
     159                // Calculate chi^2
     160                chi2->data.F32[y][x] = 0.0;
     161                int numGoodPix = 0;
     162                for (int i = 0; i < nImages; i++) {
     163                    if (! mask->data.U8[i]) {
     164                        chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
     165                        numGoodPix++;
     166                    }
     167                }
     168                chi2->data.F32[y][x] /= (float)numGoodPix;
     169#endif
     170
     171                // Rejection iterations
     172                bool keepGoing = true;  // Keep going with rejection?
     173                for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
     174                    float max = 0.0;    // Maximum deviation
     175                    int maxIndex = -1;  // Index of the maximum deviation
     176                    for (int i = 0; i < nImages; i++) {
     177                        if (!mask->data.U8[i] &&
     178                            ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
     179                            max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
     180                            maxIndex = i;
     181                        }
     182                    }
     183                    // Reject the pixel with the maximum deviation
     184                    if (max > reject) {
     185                        mask->data.U8[maxIndex] = 1;
     186                        ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
     187                        // Re-do combination following rejection
     188                        average = stacCombineMean(pixels, deltas, mask);
     189                    } else {
     190                        keepGoing = false;
     191                    }
     192                }
     193
     194            } // Pixels of interest
     195
     196        }
    197197    } // Iterating over output pixels
    198198
    … …  
    200200    // Write chi^2 image
    201201    iteration++;
    202     char chiName[MAXCHAR];              // Filename of chi^2 image
     202    char chiName[MAXCHAR];              // Filename of chi^2 image
    203203    sprintf(chiName,"chi2_%d.fits",iteration);
    204204    psFits *chiFile = psFitsAlloc(chiName);
    205205    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0)) {
    206         psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
     206        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
    207207    }
    208208    psTrace("stac.combine", 1, "Chi^2 image written to %s\n", chiName);

  • trunk/stac/src/stacInvertMaps.c

    r5743 r5745  
    99
    1010psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
    11                         int outnx, int outny // Size of output image
     11                        int outnx, int outny // Size of output image
    1212    )
    1313{
    14     int nMaps = maps->n;                // Number of maps
     14    int nMaps = maps->n;                // Number of maps
    1515    psArray *inverted = psArrayAlloc(nMaps); // Array of inverted maps for output
    1616
    … …  
    2323    for (int mapNum = 0; mapNum < nMaps; mapNum++) {
    2424
    25         psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
    26         // Check input
    27         assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
    28                oldMap->x->nX == oldMap->y->nX);
    29         int order = oldMap->x->nX + 1;  // Polynomial order
    30         psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
    31         psPlaneTransform *newMap = psPlaneTransformAlloc(order + 1, order + 1); // Inverted map
    32 
    33         // Create fake polynomial to use in evaluation
    34         psPolynomial2D *fakePoly = psPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
    35         for (int i = 0; i < order; i++) {
    36             for (int j = 0; j < order; j++) {
    37                 fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
    38                 fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
    39             }
    40         }
    41 
    42         // A grid of xin,yin --> xout,yout
    43         psVector *xIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
    44         psVector *yIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
    45         psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
    46         psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
    47 
    48         // Create grid of points
    49         for (int yint = 0; yint < NUM_GRID; yint++) {
    50             inCoord->y = (float)(yint * outny) / (float)(NUM_GRID - 1);
    51             for (int xint = 0; xint < NUM_GRID; xint++) {
    52                 inCoord->x = (float)(xint * outnx) / (float)(NUM_GRID - 1);
    53 
    54                 (void)psPlaneTransformApply(outCoord, oldMap, inCoord);
    55 
    56                 xOut->data.F32[yint*NUM_GRID + xint] = inCoord->x;
    57                 yOut->data.F32[yint*NUM_GRID + xint] = inCoord->y;
    58                 xIn->data.F32[yint*NUM_GRID + xint] = outCoord->x;
    59                 yIn->data.F32[yint*NUM_GRID + xint] = outCoord->y;
    60             }
    61         }
    62 
    63         // Initialise the matrix and vectors
    64         int nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
    65         psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
    66         psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
    67         psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
    68         for (int i = 0; i < nCoeff; i++) {
    69             for (int j = 0; j < nCoeff; j++) {
    70                 matrix->data.F64[i][j] = 0.0;
    71             }
    72             xVector->data.F64[i] = 0.0;
    73             yVector->data.F64[i] = 0.0;
    74         }
    75 
    76         // Iterate over the grid points
    77         for (int g = 0; g < NUM_GRID*NUM_GRID; g++) {
    78 
    79             // Iterate over the polynomial coefficients, accumulating the matrix and vectors
    80             for (int i = 0, ijIndex = 0; i < order; i++) {
    81                 for (int j = 0; j < order - i; j++, ijIndex++) {
    82 
    83                     fakePoly->mask[i][j] = 0;
    84                     double ijPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
    85                     fakePoly->mask[i][j] = 1;
    86 
    87                     for (int m = 0, mnIndex = 0; m < order; m++) {
    88                         for (int n = 0; n < order - m; n++, mnIndex++) {
    89                            
    90                             fakePoly->mask[m][n] = 0;
    91                             double mnPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
    92                             fakePoly->mask[m][n] = 1;
    93 
    94                             matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
    95                         }
    96                     }
    97                    
    98                     xVector->data.F64[ijIndex] += ijPoly * (double)xOut->data.F32[g];
    99                     yVector->data.F64[ijIndex] += ijPoly * (double)yOut->data.F32[g];
    100                 }
    101             } // Iterating over coefficients
    102         } // Iterating over grid points
    103 
    104         // Solution via LU Decomposition
    105         psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
    106         psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
    107         psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
    108         psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
    109 
    110         // Stuff coefficients into transformation
    111         for (int i = 0, ijIndex = 0; i < order; i++) {
    112             for (int j = 0; j < order - i; j++, ijIndex++) {
    113                 newMap->x->coeff[i][j] = xSolution->data.F64[ijIndex];
    114                 newMap->y->coeff[i][j] = ySolution->data.F64[ijIndex];
    115             }
    116         }
    117         inverted->data[mapNum] = newMap;
     25        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
     26        // Check input
     27        assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
     28               oldMap->x->nX == oldMap->y->nX);
     29        int order = oldMap->x->nX;      // Polynomial order
     30        psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
     31        psPlaneTransform *newMap = psPlaneTransformAlloc(order, order); // Inverted map
     32
     33        // Create fake polynomial to use in evaluation
     34        psPolynomial2D *fakePoly = psPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     35        for (int i = 0; i <= order; i++) {
     36            for (int j = 0; j <= order; j++) {
     37                fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
     38                fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
     39            }
     40        }
     41
     42        // A grid of xin,yin --> xout,yout
     43        psVector *xIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     44        psVector *yIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     45        psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     46        psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     47
     48        // Create grid of points
     49        for (int yint = 0; yint < NUM_GRID; yint++) {
     50            inCoord->y = (float)(yint * outny) / (float)(NUM_GRID - 1);
     51            for (int xint = 0; xint < NUM_GRID; xint++) {
     52                inCoord->x = (float)(xint * outnx) / (float)(NUM_GRID - 1);
     53
     54                (void)psPlaneTransformApply(outCoord, oldMap, inCoord);
     55
     56                xOut->data.F32[yint*NUM_GRID + xint] = inCoord->x;
     57                yOut->data.F32[yint*NUM_GRID + xint] = inCoord->y;
     58                xIn->data.F32[yint*NUM_GRID + xint] = outCoord->x;
     59                yIn->data.F32[yint*NUM_GRID + xint] = outCoord->y;
     60            }
     61        }
     62
     63        // Initialise the matrix and vectors
     64        int nCoeff = (order + 1) * (order + 2) / 2; // Number of polynomial coefficients
     65        psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
     66        psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
     67        psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     68        for (int i = 0; i < nCoeff; i++) {
     69            for (int j = 0; j < nCoeff; j++) {
     70                matrix->data.F64[i][j] = 0.0;
     71            }
     72            xVector->data.F64[i] = 0.0;
     73            yVector->data.F64[i] = 0.0;
     74        }
     75
     76        // Iterate over the grid points
     77        for (int g = 0; g < NUM_GRID*NUM_GRID; g++) {
     78
     79            // Iterate over the polynomial coefficients, accumulating the matrix and vectors
     80            for (int i = 0, ijIndex = 0; i <= order; i++) {
     81                for (int j = 0; j <= order - i; j++, ijIndex++) {
     82
     83                    fakePoly->mask[i][j] = 0;
     84                    double ijPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
     85                    fakePoly->mask[i][j] = 1;
     86
     87                    for (int m = 0, mnIndex = 0; m <= order; m++) {
     88                        for (int n = 0; n <= order - m; n++, mnIndex++) {
     89
     90                            fakePoly->mask[m][n] = 0;
     91                            double mnPoly = psPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
     92                            fakePoly->mask[m][n] = 1;
     93
     94                            matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
     95                        }
     96                    }
     97
     98                    xVector->data.F64[ijIndex] += ijPoly * (double)xOut->data.F32[g];
     99                    yVector->data.F64[ijIndex] += ijPoly * (double)yOut->data.F32[g];
     100                }
     101            } // Iterating over coefficients
     102        } // Iterating over grid points
     103
     104        // Solution via LU Decomposition
     105        psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     106        psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
     107        psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
     108        psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     109
     110        // Stuff coefficients into transformation
     111        for (int i = 0, ijIndex = 0; i <= order; i++) {
     112            for (int j = 0; j <= order - i; j++, ijIndex++) {
     113                newMap->x->coeff[i][j] = xSolution->data.F64[ijIndex];
     114                newMap->y->coeff[i][j] = ySolution->data.F64[ijIndex];
     115            }
     116        }
     117        inverted->data[mapNum] = newMap;
    118118
    119119#ifdef TESTING
    120         // Print x coefficients
    121         psTrace("stac.invertMaps", 7, "x' = \n");
    122         for (int i = 0; i < order; i++) {
    123             for (int j = 0; j < order - i; j++) {
    124                 psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->x->coeff[i][j], i, j);
    125             }
    126         }
    127         // Print y coefficients
    128         psTrace("stac.invertMaps", 7, "y' = \n");
    129         for (int i = 0; i < order; i++) {
    130             for (int j = 0; j < order - i; j++) {
    131                 psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->y->coeff[i][j], i, j);
    132             }
    133         }
     120        // Print x coefficients
     121        psTrace("stac.invertMaps", 7, "x' = \n");
     122        for (int i = 0; i <= order; i++) {
     123            for (int j = 0; j <= order - i; j++) {
     124                psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->x->coeff[i][j], i, j);
     125            }
     126        }
     127        // Print y coefficients
     128        psTrace("stac.invertMaps", 7, "y' = \n");
     129        for (int i = 0; i <= order; i++) {
     130            for (int j = 0; j <= order - i; j++) {
     131                psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->y->coeff[i][j], i, j);
     132            }
     133        }
    134134#endif
    135135
    136136#ifdef TESTING
    137         // Go forward then backward
    138         double x = 123.4;
    139         double y = 432.1;
    140         psPlane *oldCoords = psAlloc(sizeof(psPlane));
    141         oldCoords->x = x;
    142         oldCoords->y = y;
    143         psPlane *newCoords = psPlaneTransformApply(NULL, oldMap, oldCoords);
    144         psTrace("stac.invertMaps.test", 5, "%f,%f --> %f,%f\n", x, y, newCoords->x, newCoords->y);
    145         (void)psPlaneTransformApply(oldCoords, newMap, newCoords);
    146         psTrace("stac.invertMaps.test", 5, "--------> %f,%f\n", oldCoords->x, oldCoords->y);
    147         psFree(newCoords);
    148         psFree(oldCoords);
     137        // Go forward then backward
     138        double x = 123.4;
     139        double y = 432.1;
     140        psPlane *oldCoords = psAlloc(sizeof(psPlane));
     141        oldCoords->x = x;
     142        oldCoords->y = y;
     143        psPlane *newCoords = psPlaneTransformApply(NULL, oldMap, oldCoords);
     144        psTrace("stac.invertMaps.test", 5, "%f,%f --> %f,%f\n", x, y, newCoords->x, newCoords->y);
     145        (void)psPlaneTransformApply(oldCoords, newMap, newCoords);
     146        psTrace("stac.invertMaps.test", 5, "--------> %f,%f\n", oldCoords->x, oldCoords->y);
     147        psFree(newCoords);
     148        psFree(oldCoords);
    149149#endif
    150150
    151         psFree(permutation);
    152         psFree(luMatrix);
    153         psFree(matrix);
    154         psFree(xVector);
    155         psFree(yVector);
    156         psFree(xSolution);
    157         psFree(ySolution);
    158         psFree(fakePoly);
    159         psFree(xIn);
    160         psFree(yIn);
    161         psFree(xOut);
    162         psFree(yOut);
     151        psFree(permutation);
     152        psFree(luMatrix);
     153        psFree(matrix);
     154        psFree(xVector);
     155        psFree(yVector);
     156        psFree(xSolution);
     157        psFree(ySolution);
     158        psFree(fakePoly);
     159        psFree(xIn);
     160        psFree(yIn);
     161        psFree(xOut);
     162        psFree(yOut);
    163163    }
    164164
    165165    psFree(inCoord);
    166166    psFree(outCoord);
    167                
     167
    168168
    169169#if 0
    170         // Can't handle higher order than linear yet
    171         assert(((psPlaneTransform*)maps->data[i])->x->nX == 2 &&
    172                ((psPlaneTransform*)maps->data[i])->x->nY == 2 &&
    173                ((psPlaneTransform*)maps->data[i])->y->nX == 2 &&
    174                ((psPlaneTransform*)maps->data[i])->y->nY == 2);
    175         psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
    176         psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
    177 
    178         // Now, simply do a 2x2 matrix inversion
    179 
    180         double a = oldMap->x->coeff[1][0];
    181         double b = oldMap->x->coeff[0][1];
    182         double c = oldMap->y->coeff[1][0];
    183         double d = oldMap->y->coeff[0][1];
    184         double e = oldMap->x->coeff[0][0];
    185         double f = oldMap->y->coeff[0][0];
    186 
    187         double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
    188 
    189         // Not entirely sure why this works, but it appears to do so.......................................
    190         newMap->x->coeff[1][0] = invDet * a;
    191         newMap->x->coeff[0][1] = - invDet * b;
    192         newMap->y->coeff[1][0] = - invDet * c;
    193         newMap->y->coeff[0][1] = invDet * d;
    194 
    195         newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
    196         newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
     170        // Can't handle higher order than linear yet
     171        assert(((psPlaneTransform*)maps->data[i])->x->nX == 2 &&
     172               ((psPlaneTransform*)maps->data[i])->x->nY == 2 &&
     173               ((psPlaneTransform*)maps->data[i])->y->nX == 2 &&
     174               ((psPlaneTransform*)maps->data[i])->y->nY == 2);
     175        psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
     176        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
     177
     178        // Now, simply do a 2x2 matrix inversion
     179
     180        double a = oldMap->x->coeff[1][0];
     181        double b = oldMap->x->coeff[0][1];
     182        double c = oldMap->y->coeff[1][0];
     183        double d = oldMap->y->coeff[0][1];
     184        double e = oldMap->x->coeff[0][0];
     185        double f = oldMap->y->coeff[0][0];
     186
     187        double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
     188
     189        // Not entirely sure why this works, but it appears to do so.......................................
     190        newMap->x->coeff[1][0] = invDet * a;
     191        newMap->x->coeff[0][1] = - invDet * b;
     192        newMap->y->coeff[1][0] = - invDet * c;
     193        newMap->y->coeff[0][1] = invDet * d;
     194
     195        newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
     196        newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
    197197#endif
    198198

  • trunk/stac/src/stacRead.c

    r5743 r5745  
    124124    psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
    125125
    126     psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
     126    psPlaneTransform *map = psPlaneTransformAlloc(order, order); // The transformation
    127127    // Set coefficient masks
    128     for (int i = 0; i < order + 1; i++) {
    129         for (int j = 0; j < order + 1; j++) {
     128    for (int i = 0; i <= order; i++) {
     129        for (int j = 0; j <= order; j++) {
    130130            if (i + j > order + 1) {
    131131                map->x->mask[i][j] = 1;
    … …  
    140140    // Read x coefficients
    141141    psTrace("stac.read.map", 7, "x' = \n");
    142     for (int k = 0; k < order + 1; k++) {
    143         for (int j = 0; j < k + 1; j++) {
     142    for (int k = 0; k <= order; k++) {
     143        for (int j = 0; j <= k; j++) {
    144144            int i = k - j;
    145145            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
    … …  
    154154    // Read y coefficients
    155155    psTrace("stac.read.maps", 7, "y' = \n");
    156     for (int k = 0; k < order + 1; k++) {
    157         for (int j = 0; j < k + 1; j++) {
     156    for (int k = 0; k <= order; k++) {
     157        for (int j = 0; j <= k; j++) {
    158158            int i = k - j;
    159159            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
    … …  
    188188        }
    189189        // Read the file
    190         sprintf(mapfile,"%s.map",filenames->data[i]);
     190        sprintf(mapfile, "%s.map", (const char*)filenames->data[i]);
    191191        maps->data[i] = stacReadMap(mapfile);
    192192        if (maps->data[i] == NULL) {

  • trunk/stac/src/stacScales.c

    r5743 r5745  
    66#define SQUARE(x) ((x)*(x))
    77
    8 #define SAMPLE 10                       // Subsample rate for images
     8#define SAMPLE 10                       // Subsample rate for images
    99
    1010float stacBackground(const psImage *image, // Image for which to get the background
    11                      int sample         // Sample in increments of this value
     11                     int sample         // Sample in increments of this value
    1212    )
    1313{
    … …  
    1515
    1616#if 1
    17     int size = image->numCols * image->numRows; // Number of pixels in image
    18     int numSamples = size / sample;     // Number of samples in image
     17    int size = image->numCols * image->numRows; // Number of pixels in image
     18    int numSamples = size / sample;     // Number of samples in image
    1919    psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
    20     psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
    21 
    22     int offset = 0;                     // Offset from start of the row
    23     int index = 0;                      // Sample number
     20    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8); // Mask for sample
     21
     22    int offset = 0;                     // Offset from start of the row
     23    int index = 0;                      // Sample number
    2424    for (int row = 0; row < image->numRows; row++) {
    25         // I'd cast this as a "for", but this makes it a bit easier to understand.
    26         int col = offset;
    27         while (col < image->numCols) {
    28             values->data.F32[index] = image->data.F32[row][col];
    29             if (isnan(values->data.F32[index])) {
    30                 mask->data.U8[index] = 1;
    31             } else {
    32                 mask->data.U8[index] = 0;
    33             }
    34             col += sample;
    35             index++;
    36         }
    37         offset = col - image->numCols;
     25        // I'd cast this as a "for", but this makes it a bit easier to understand.
     26        int col = offset;
     27        while (col < image->numCols) {
     28            values->data.F32[index] = image->data.F32[row][col];
     29            if (isnan(values->data.F32[index])) {
     30                mask->data.U8[index] = 1;
     31            } else {
     32                mask->data.U8[index] = 0;
     33            }
     34            col += sample;
     35            index++;
     36        }
     37        offset = col - image->numCols;
    3838    }
    3939
    … …  
    5858
    5959
    60 bool stacScales(psVector **scalesPtr,   // Scales to return
    61                 psVector **offsetsPtr,  // Offsets to return
    62                 const psArray *images,  // Images on which to measure the scales and offsets
    63                 const char *starFile,   // File containing coordinates to photometer
    64                 const char *starMapFile, // Map for coodinates to the common output frame
    65                 float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
    66                 float aper              // Aperture to use for photometry (radius)
     60bool stacScales(psVector **scalesPtr,   // Scales to return
     61                psVector **offsetsPtr,  // Offsets to return
     62                const psArray *images,  // Images on which to measure the scales and offsets
     63                const char *starFile,   // File containing coordinates to photometer
     64                const char *starMapFile, // Map for coodinates to the common output frame
     65                float xMapDiff, float yMapDiff, // Difference from the map to apply (due to size difference)
     66                float aper              // Aperture to use for photometry (radius)
    6767    )
    6868{
    … …  
    7171    assert(images);
    7272    for (int i = 0; i < images->n; i++) {
    73         psImage *image = images->data[i];
    74         assert(image->type.type == PS_TYPE_F32);
     73        psImage *image = images->data[i];
     74        assert(image->type.type == PS_TYPE_F32);
    7575    }
    7676
    … …  
    7878    psVector *offsets = NULL; // Offsets between images (ADU)
    7979    if (*scalesPtr) {
    80         scales = *scalesPtr;
    81         assert(scales);
    82         assert(scales->n == images->n);
    83         assert(scales->type.type == PS_TYPE_F32);
     80        scales = *scalesPtr;
     81        assert(scales);
     82        assert(scales->n == images->n);
     83        assert(scales->type.type == PS_TYPE_F32);
    8484    } else {
    85         scales = psVectorAlloc(images->n, PS_TYPE_F32);
    86         *scalesPtr = scales;
     85        scales = psVectorAlloc(images->n, PS_TYPE_F32);
     86        *scalesPtr = scales;
    8787    }
    8888    if (*offsetsPtr) {
    89         offsets = *offsetsPtr;
    90         assert(offsets);
    91         assert(offsets->n == images->n);
    92         assert(offsets->type.type == PS_TYPE_F32);
     89        offsets = *offsetsPtr;
     90        assert(offsets);
     91        assert(offsets->n == images->n);
     92        assert(offsets->type.type == PS_TYPE_F32);
    9393    } else {
    94         offsets = psVectorAlloc(images->n, PS_TYPE_F32);
    95         *offsetsPtr = offsets;
     94        offsets = psVectorAlloc(images->n, PS_TYPE_F32);
     95        *offsetsPtr = offsets;
    9696    }
    9797
    … …  
    100100    double startTime = getTime();
    101101    for (int i = 0; i < images->n; i++) {
    102         offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE);
    103         psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]);
    104         double time = getTime();
    105         psTrace("stac.scales", 10, "Took %f sec\n", time - startTime);
    106         startTime = time;
     102        offsets->data.F32[i] = stacBackground(images->data[i], SAMPLE);
     103        psTrace("stac.scales", 5, "Background in image %d is %f\n", i, offsets->data.F32[i]);
     104        double time = getTime();
     105        psTrace("stac.scales", 10, "Took %f sec\n", time - startTime);
     106        startTime = time;
    107107    }
    108108
    109109    // Now the scales
    110110    if (starFile == NULL || starMapFile == NULL) {
    111         psLogMsg("stac.scales", PS_LOG_INFO,
    112                 "No coordinates available to set scales --- assuming all are identical.\n");
    113         for (int i = 0; i < images->n; i++) {
    114             scales->data.F32[i] = 1.0;
    115             psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
    116         }
     111        psLogMsg("stac.scales", PS_LOG_INFO,
     112                "No coordinates available to set scales --- assuming all are identical.\n");
     113        for (int i = 0; i < images->n; i++) {
     114            scales->data.F32[i] = 1.0;
     115            psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
     116        }
    117117    } else {
    118         // Read star coordinates and map
    119         psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
    120         psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
    121 
    122         // Transform the stellar positions to match the transformed reference frame
    123         psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
    124         // Fix up difference between map and output frame
    125         starMap->x->coeff[0][0] -= xMapDiff;
    126         starMap->y->coeff[0][0] -= yMapDiff;
    127         for (int i = 0; i < starCoords->n; i++) {
    128             starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
    129         }
    130         psFree(starMap);
    131 
    132         psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
    133         psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
    134 
    135         // Set scales relative to the first image
    136         scales->data.F32[0] = 1.0;
    137 
    138         // Iterate over images
    139         for (int i = 0; i < scales->n; i++) {
    140             psImage *image = images->data[i]; // The image we're working with
    141 
    142             float background = offsets->data.F32[i]; // Background in image
    143 
    144             // Do photometry on transformed image i
    145             psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i);
    146             psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
    147             psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
    148             for (int j = 0; j < starCoordsTransformed->n; j++) {
    149                 psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
    150                
    151                 if (coords->x < aper || coords->y < aper ||
    152                     coords->x + aper > image->numCols - 1 ||
    153                     coords->y + aper > image->numRows) {
    154                     mask->data.U8[j] = 1;
    155                 } else {
    156                     // Sum flux within the aperture
    157                     float sum = 0.0;
    158                     int numPix = 0;
    159                     float aper2 = SQUARE(aper);
    160                     for (int y = (int)floorf(coords->y - aper);
    161                          y <= (int)ceilf(coords->y + aper); y++) {
    162                         for (int x = (int)floorf(coords->x - aper);
    163                              x <= (int)ceilf(coords->x + aper); x++) {
    164                             if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
    165                                 aper2) {
    166                                 sum += image->data.F32[y][x];
    167                                 numPix++;
    168                             }
    169                         }
    170                     }
    171                     // Subtract background, renormalise to account for circular aperture
    172                     if (numPix > 0 && sum > 0) {
    173                         sum -= offsets->data.F32[i] * (float)numPix;
    174                         photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
    175                         if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
    176                             mask->data.U8[j] = 1;
    177                             psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
    178                         } else {
    179                             mask->data.U8[j] = 0;
    180                         }
    181                     } else {
    182                         mask->data.U8[j] = 0;
    183                     }
    184                 }
    185             }
    186             stars->data[i] = photometry;
    187             masks->data[i] = mask;
    188         }
    189         psFree(starCoords);
    190         psFree(starCoordsTransformed);
    191 
    192         // Get the scales
    193         psVector *ref = stars->data[0]; // The reference photometry
    194         psVector *refMask = masks->data[0]; // The reference mask
     118        // Read star coordinates and map
     119        psArray *starCoords = stacReadCoords(starFile); // Array of star coordinates
     120        psPlaneTransform *starMap = stacReadMap(starMapFile); // Transformation for star coordinates
     121
     122        // Transform the stellar positions to match the transformed reference frame
     123        psArray *starCoordsTransformed = psArrayAlloc(starCoords->n); // Transformed positions
     124        // Fix up difference between map and output frame
     125        starMap->x->coeff[0][0] -= xMapDiff;
     126        starMap->y->coeff[0][0] -= yMapDiff;
     127        for (int i = 0; i < starCoords->n; i++) {
     128            starCoordsTransformed->data[i] = psPlaneTransformApply(NULL, starMap, starCoords->data[i]);
     129        }
     130        psFree(starMap);
     131
     132        psArray *stars = psArrayAlloc(images->n); // Array of stellar photometry vectors
     133        psArray *masks = psArrayAlloc(images->n); // Array of masks for stars
     134
     135        // Set scales relative to the first image
     136        scales->data.F32[0] = 1.0;
     137
     138        // Iterate over images
     139        for (int i = 0; i < scales->n; i++) {
     140            psImage *image = images->data[i]; // The image we're working with
     141
     142            // Do photometry on transformed image i
     143            psTrace("stac.scales", 3, "Doing photometry on image %d....\n", i);
     144            psVector *photometry = psVectorAlloc(starCoords->n, PS_TYPE_F32); // Photometry of the stars
     145            psVector *mask = psVectorAlloc(starCoords->n, PS_TYPE_U8); // Mask for the photometry
     146            for (int j = 0; j < starCoordsTransformed->n; j++) {
     147                psPlane *coords = starCoordsTransformed->data[j]; // The coordinates of the star
     148
     149                if (coords->x < aper || coords->y < aper ||
     150                    coords->x + aper > image->numCols - 1 ||
     151                    coords->y + aper > image->numRows) {
     152                    mask->data.U8[j] = 1;
     153                } else {
     154                    // Sum flux within the aperture
     155                    float sum = 0.0;
     156                    int numPix = 0;
     157                    float aper2 = SQUARE(aper);
     158                    for (int y = (int)floorf(coords->y - aper);
     159                         y <= (int)ceilf(coords->y + aper); y++) {
     160                        for (int x = (int)floorf(coords->x - aper);
     161                             x <= (int)ceilf(coords->x + aper); x++) {
     162                            if (SQUARE((float)x + 0.5 - coords->x) + SQUARE((float)y + 0.5 - coords->y) <=
     163                                aper2) {
     164                                sum += image->data.F32[y][x];
     165                                numPix++;
     166                            }
     167                        }
     168                    }
     169                    // Subtract background, renormalise to account for circular aperture
     170                    if (numPix > 0 && sum > 0) {
     171                        sum -= offsets->data.F32[i] * (float)numPix;
     172                        photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
     173                        if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
     174                            mask->data.U8[j] = 1;
     175                            psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
     176                        } else {
     177                            mask->data.U8[j] = 0;
     178                        }
     179                    } else {
     180                        mask->data.U8[j] = 0;
     181                    }
     182                }
     183            }
     184            stars->data[i] = photometry;
     185            masks->data[i] = mask;
     186        }
     187        psFree(starCoords);
     188        psFree(starCoordsTransformed);
     189
     190        // Get the scales
     191        psVector *ref = stars->data[0]; // The reference photometry
     192        psVector *refMask = masks->data[0]; // The reference mask
    195193#if 0
    196         psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
    197         stats->clipSigma = 2.5;
    198         stats->clipIter = 3;
     194        psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
     195        stats->clipSigma = 2.5;
     196        stats->clipIter = 3;
    199197#else
    200         psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
    201 #endif
    202         for (int i = 1; i < scales->n; i++) {
    203             psVector *input = stars->data[i];   // The comparison photometry
    204             psVector *inputMask = masks->data[i]; // The comparison mask
    205 
    206             psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
    207             psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
    208             (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
     198        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
     199#endif
     200        for (int i = 1; i < scales->n; i++) {
     201            psVector *input = stars->data[i];   // The comparison photometry
     202            psVector *inputMask = masks->data[i]; // The comparison mask
     203
     204            psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
     205            psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
     206            (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
    209207
    210208#if 0
    211             psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
    212             for (int j = 0; j < compare->n; j++) {
    213                 if (compareMask->data.U8[j] & 1) {
    214                     psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
    215                             ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
    216                 } else {
    217                     psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
    218                             ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
    219                 }
    220             }
    221 #endif
    222 
    223             psVectorStats(stats, compare, NULL, compareMask, 1);
     209            psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
     210            for (int j = 0; j < compare->n; j++) {
     211                if (compareMask->data.U8[j] & 1) {
     212                    psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
     213                            ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
     214                } else {
     215                    psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
     216                            ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
     217                }
     218            }
     219#endif
     220
     221            psVectorStats(stats, compare, NULL, compareMask, 1);
    224222
    225223#if 0
    226             scales->data.F32[i] = stats->clippedMean;
     224            scales->data.F32[i] = stats->clippedMean;
    227225#else
    228             scales->data.F32[i] = stats->sampleMedian;
    229 #endif
    230             psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
    231             psFree(compare);
    232             psFree(compareMask);
    233         }
    234         psFree(stats);
    235 
    236         psFree(stars);
    237         psFree(masks);
     226            scales->data.F32[i] = stats->sampleMedian;
     227#endif
     228            psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
     229            psFree(compare);
     230            psFree(compareMask);
     231        }
     232        psFree(stats);
     233
     234        psFree(stars);
     235        psFree(masks);
    238236    }
    239237
    … …  
    242240
    243241
    244 bool stacRescale(psArray *images,       // Images to rescale
    245                  psArray *errImages,    // Variance images to rescale
    246                  const psImage *mask,   // Mask for pixels to scale
    247                 const psVector *scales,// Scales for images
    248                 const psVector *offsets // Offsets for images
     242bool stacRescale(psArray *images,       // Images to rescale
     243                 psArray *errImages,    // Variance images to rescale
     244                 const psImage *mask,   // Mask for pixels to scale
     245                const psVector *scales,// Scales for images
     246                const psVector *offsets // Offsets for images
    249247    )
    250248{
    … …  
    258256    assert(offsets->type.type == PS_TYPE_F32);
    259257    for (int i = 0; i < images->n; i++) {
    260         psImage *image = images->data[i]; // Image of interest
    261         assert(image->type.type == PS_TYPE_F32);
    262         assert(scales->data.F32[i] != 0);
    263         if (mask) {
    264             assert(mask->type.type == PS_TYPE_U8);
    265             assert(image->numCols == mask->numCols && image->numRows == mask->numRows);
    266         }
    267         if (errImages) {
    268             psImage *errImage = errImages->data[i];
    269             assert(errImage->type.type == PS_TYPE_F32);
    270             assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows);
    271         }
    272     }
    273 
    274     for (int i = 0; i < images->n; i++) {
    275         psImage *image = images->data[i]; // Image to rescale
    276         psImage *errImage = NULL;       // Variance image to rescale
    277         if (errImages) {
    278             errImage = errImages->data[i];
    279         }
    280         float scale = scales->data.F32[i]; // Scale to use
    281         float offset = offsets->data.F32[i]; // Offset to use
    282         for (int y = 0; y < image->numRows; y++) {
    283             for (int x = 0; x < image->numCols; x++) {
    284                 if (!mask || mask->data.F32[y][x]) {
    285                     image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale;
    286                     if (errImage) {
    287                         errImage->data.F32[y][x] /= SQUARE(scale);
    288                     }
    289                 }
    290             }
    291         }
     258        psImage *image = images->data[i]; // Image of interest
     259        assert(image->type.type == PS_TYPE_F32);
     260        assert(scales->data.F32[i] != 0);
     261        if (mask) {
     262            assert(mask->type.type == PS_TYPE_U8);
     263            assert(image->numCols == mask->numCols && image->numRows == mask->numRows);
     264        }
     265        if (errImages) {
     266            psImage *errImage = errImages->data[i];
     267            assert(errImage->type.type == PS_TYPE_F32);
     268            assert(errImage->numCols == image->numCols && errImage->numRows == image->numRows);
     269        }
     270    }
     271
     272    for (int i = 0; i < images->n; i++) {
     273        psImage *image = images->data[i]; // Image to rescale
     274        psImage *errImage = NULL;       // Variance image to rescale
     275        if (errImages) {
     276            errImage = errImages->data[i];
     277        }
     278        float scale = scales->data.F32[i]; // Scale to use
     279        float offset = offsets->data.F32[i]; // Offset to use
     280        for (int y = 0; y < image->numRows; y++) {
     281            for (int x = 0; x < image->numCols; x++) {
     282                if (!mask || mask->data.F32[y][x]) {
     283                    image->data.F32[y][x] = (image->data.F32[y][x] - offset) / scale;
     284                    if (errImage) {
     285                        errImage->data.F32[y][x] /= SQUARE(scale);
     286                    }
     287                }
     288            }
     289        }
    292290    }
    293291

  • trunk/stac/src/stacTransform.c

    r5743 r5745  
    1010// Hacked the original ps_ImagePixelInterpolateBILINEAR_F32 to add variances
    1111// i.e., to square the fractions when combining.
    12 inline psF64 p_psImageErrorInterpolateBILINEAR_F32(const psImage* input, 
    13                                                    float x,
    14                                                    float y,
    15                                                    const psImage* mask,
    16                                                    unsigned int maskVal,
    17                                                    psF64 unexposedValue)
     12inline psF64 p_psImageErrorInterpolateBILINEAR_F32(const psImage* input,
     13                                                   float x,
     14                                                   float y,
     15                                                   const psImage* mask,
     16                                                   unsigned int maskVal,
     17                                                   psF64 unexposedValue)
    1818{
    19     double floorX = floor((psF64)(x) - 0.5); 
    20     double floorY = floor((psF64)(y) - 0.5); 
    21     psF64 fracX = x - 0.5 - floorX; 
    22     psF64 fracY = y - 0.5 - floorY; 
    23     int intFloorX = (int) floorX; 
    24     int intFloorY = (int) floorY; 
    25     int lastX = input->numCols - 1; 
    26     int lastY = input->numRows - 1; 
    27     psF32 V00;
    28     psF32 V01;
    29     psF32 V10;
    30     psF32 V11;
    31     bool valid00 = false; 
    32     bool valid01 = false; 
    33     bool valid10 = false; 
    34     bool valid11 = false; 
    35    
    36     if (intFloorY >= 0 && intFloorY <= lastY) { 
    37         if (intFloorX >= 0 && intFloorX <= lastX) { 
    38             V00 = input->data.F32[intFloorY][intFloorX]; 
    39             valid00 = true;
    40         } 
    41         if (intFloorX >= -1 && intFloorX < lastX) { 
    42             V10 = input->data.F32[intFloorY][intFloorX+1]; 
    43             valid10 = true;
    44         } 
    45     } 
    46     if (intFloorY >= -1 && intFloorY < lastY) { 
    47         if (intFloorX >= 0 && intFloorX <= lastX) { 
    48             V01 = input->data.F32[intFloorY+1][intFloorX]; 
    49             valid01 = true;
    50         } 
    51         if (intFloorX >= -1 && intFloorX < lastX) { 
    52             V11 = input->data.F32[intFloorY+1][intFloorX+1]; 
    53             valid11 = true;
    54         } 
    55     } 
    56    
    57     /* cover likely case of all pixels being valid more efficiently */ 
    58     if (valid00 && valid10 && valid01 && valid11) { 
    59         /* formula from the ADD */ 
     19    double floorX = floor((psF64)(x) - 0.5);
     20    double floorY = floor((psF64)(y) - 0.5);
     21    psF64 fracX = x - 0.5 - floorX;
     22    psF64 fracY = y - 0.5 - floorY;
     23    int intFloorX = (int) floorX;
     24    int intFloorY = (int) floorY;
     25    int lastX = input->numCols - 1;
     26    int lastY = input->numRows - 1;
     27    psF32 V00 = 0.0;
     28    psF32 V01 = 0.0;
     29    psF32 V10 = 0.0;
     30    psF32 V11 = 0.0;
     31    bool valid00 = false;
     32    bool valid01 = false;
     33    bool valid10 = false;
     34    bool valid11 = false;
     35
     36    if (intFloorY >= 0 && intFloorY <= lastY) {
     37        if (intFloorX >= 0 && intFloorX <= lastX) {
     38            V00 = input->data.F32[intFloorY][intFloorX];
     39            valid00 = true;
     40        }
     41        if (intFloorX >= -1 && intFloorX < lastX) {
     42            V10 = input->data.F32[intFloorY][intFloorX+1];
     43            valid10 = true;
     44        }
     45    }
     46    if (intFloorY >= -1 && intFloorY < lastY) {
     47        if (intFloorX >= 0 && intFloorX <= lastX) {
     48            V01 = input->data.F32[intFloorY+1][intFloorX];
     49            valid01 = true;
     50        }
     51        if (intFloorX >= -1 && intFloorX < lastX) {
     52            V11 = input->data.F32[intFloorY+1][intFloorX+1];
     53            valid11 = true;
     54        }
     55    }
     56
     57    /* cover likely case of all pixels being valid more efficiently */
     58    if (valid00 && valid10 && valid01 && valid11) {
     59        /* formula from the ADD */
    6060        return V00*SQUARE((1.0-fracX)*(1.0-fracY)) + V10*SQUARE(fracX*(1.0-fracY)) +
    61             V01*SQUARE(fracY*(1.0-fracX)) + V11*SQUARE(fracX*fracY);
    62     } 
    63    
    64     /* OK, at least one pixel is not valid - need to do it piecemeal */ 
    65    
    66     psF64 V0;
    67     bool valid0 = true; 
    68     if (valid00 && valid10) { 
    69         V0 = V00*SQUARE(1-fracX)+V10*SQUARE(fracX); 
    70     } else if (valid00) { 
    71         V0 = V00; 
    72     } else if (valid10) { 
    73         V0 = V10; 
    74     } else { 
    75         valid0 = false; 
    76     } 
    77    
    78     psF64 V1;
    79     bool valid1 = true; 
    80     if (valid01 && valid11) { 
    81         V1 = V01*SQUARE(1-fracX)+V11*SQUARE(fracX); 
    82     } else if (valid01) { 
    83         V1 = V01; 
    84     } else if (valid11) { 
    85         V1 = V11; 
    86     } else { 
    87         valid1 = false; 
    88     } 
    89    
    90     if (valid0 && valid1) { 
    91         return ( V0*SQUARE(1-fracY) + V1*SQUARE(fracY) ); 
    92     } else if (valid0) { 
    93         return V0; 
    94     } else if (valid1) { 
    95         return V1; 
    96     } 
    97    
    98     return unexposedValue; 
     61            V01*SQUARE(fracY*(1.0-fracX)) + V11*SQUARE(fracX*fracY);
     62    }
     63
     64    /* OK, at least one pixel is not valid - need to do it piecemeal */
     65
     66    psF64 V0 = 0.0;
     67    bool valid0 = true;
     68    if (valid00 && valid10) {
     69        V0 = V00*SQUARE(1-fracX)+V10*SQUARE(fracX);
     70    } else if (valid00) {
     71        V0 = V00;
     72    } else if (valid10) {
     73        V0 = V10;
     74    } else {
     75        valid0 = false;
     76    }
     77
     78    psF64 V1 = 0.0;
     79    bool valid1 = true;
     80    if (valid01 && valid11) {
     81        V1 = V01*SQUARE(1-fracX)+V11*SQUARE(fracX);
     82    } else if (valid01) {
     83        V1 = V01;
     84    } else if (valid11) {
     85        V1 = V11;
     86    } else {
     87        valid1 = false;
     88    }
     89
     90    if (valid0 && valid1) {
     91        return ( V0*SQUARE(1-fracY) + V1*SQUARE(fracY) );
     92    } else if (valid0) {
     93        return V0;
     94    } else if (valid1) {
     95        return V1;
     96    }
     97
     98    return unexposedValue;
    9999}
    100100
    101101
    102102
    103 bool stacTransform(psArray **outputs,   // Transformed images for output
    104                    psArray **outErrors, // Transformed error images for output
    105                    const psArray *images, // Array of images to be transformed
    106                    const psArray *maps, // Array of polynomials that do the transformation
    107                    const psArray *errors, // Array of error images to be transformed
    108                    const psArray *masks, // Masks of input images
    109                    const psImage *region, // Region of interest for transformation
    110                    const psVector *scales, // Relative scales
    111                    const psVector *offsets, // Relative offsets
    112                    int outnx, int outny // Size of output images
     103bool stacTransform(psArray **outputs,   // Transformed images for output
     104                   psArray **outErrors, // Transformed error images for output
     105                   const psArray *images, // Array of images to be transformed
     106                   const psArray *maps, // Array of polynomials that do the transformation
     107                   const psArray *errors, // Array of error images to be transformed
     108                   const psArray *masks, // Masks of input images
     109                   const psImage *region, // Region of interest for transformation
     110                   const psVector *scales, // Relative scales
     111                   const psVector *offsets, // Relative offsets
     112                   int outnx, int outny // Size of output images
    113113    )
    114114{
    115     int nImages = images->n;            // Number of images
     115    int nImages = images->n;            // Number of images
    116116
    117117    // Check input sizes
    … …  
    126126    assert(!*outputs || (*outputs)->n == nImages);
    127127    if (*outputs == NULL) {
    128         *outputs = psArrayAlloc(nImages);
    129         psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
    130         for (int i = 0; i < nImages; i++) {
    131             (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
    132             psImageInit((*outputs)->data[i], 0.0);
    133         }
     128        *outputs = psArrayAlloc(nImages);
     129        psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
     130        for (int i = 0; i < nImages; i++) {
     131            (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
     132            psImageInit((*outputs)->data[i], 0.0);
     133        }
    134134    }
    135135
    … …  
    137137    assert(!errors || ! *outErrors || errors->n == (*outErrors)->n);
    138138    if (errors && (*outErrors == NULL)) {
    139         *outErrors = psArrayAlloc(errors->n);
    140         psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
    141         for (int i = 0; i < nImages; i++) {
    142             (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
    143         }
     139        *outErrors = psArrayAlloc(errors->n);
     140        psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
     141        for (int i = 0; i < nImages; i++) {
     142            (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
     143        }
    144144    }
    145145
    … …  
    147147    assert(!masks || masks->n == nImages);
    148148    if (masks != NULL) {
    149         for (int i = 0; i < nImages; i++) {
    150             psImage *image = images->data[i];
    151             psImage *mask = masks->data[i];
    152             assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
    153         }
     149        for (int i = 0; i < nImages; i++) {
     150            psImage *image = images->data[i];
     151            psImage *mask = masks->data[i];
     152            assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
     153        }
    154154    }
    155155
    … …  
    162162    // Iterate over the images
    163163    for (int n = 0; n < nImages; n++) {
    164         psTrace("stac.transform", 1, "Transforming image %d....\n",n);
    165 
    166         // Pull out the various stuff we're working on
    167         psImage *image = images->data[n]; // The input image
    168         psPlaneTransform *map = maps->data[n]; // The map
    169         psImage *outImage = (*outputs)->data[n]; // The output image
    170         psImage *error = NULL; // The error image
    171         psImage *outError = NULL; // The output error image
    172         if (errors) {
    173             error = errors->data[n];
    174             outError = (*outErrors)->data[n];
    175         }
    176         float offset = 0.0;             // Relative offset
    177         float scale = 1.0;              // Relative scale
    178         if (offsets) {
    179             offset = offsets->data.F32[n];
    180         }
    181         if (scales) {
    182             scale = scales->data.F32[n];
    183         }
    184 
    185         // Mask
    186         psImage *mask = NULL;
    187         if (masks != NULL) {
    188             mask = masks->data[n];
    189         }
    190 
    191         // Iterate over the output image pixels
    192         for (int y = 0; y < outny; y++) {
    193             for (int x = 0; x < outnx; x++) {
    194                 // Only transform those pixels requested
    195                 if (!region || (region && region->data.U8[y][x])) {
    196                     // Transform!
    197                     sky->x = (double)x + 0.5;
    198                     sky->y = (double)y + 0.5;
    199                     (void)psPlaneTransformApply(detector, map, sky);
    200                    
    201                     // Change PS_INTERPOLATE_BILINEAR to best available technique.
    202                     outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
    203                                                                               detector->y, mask, 1, NAN,
    204                                                                               PS_INTERPOLATE_BILINEAR);
    205                     if (error) {
    206                         // Error is actually the variance
    207                         outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
    208                                                                                                 detector->x,
    209                                                                                                 detector->y,
    210                                                                                                 mask, 1, NAN);
    211                     }
    212 
    213                     outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
    214                     outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
    215 
    216                 } // Pixels of interest
    217 
    218             }
    219         } // Iterating over output pixels
     164        psTrace("stac.transform", 1, "Transforming image %d....\n",n);
     165
     166        // Pull out the various stuff we're working on
     167        psImage *image = images->data[n]; // The input image
     168        psPlaneTransform *map = maps->data[n]; // The map
     169        psImage *outImage = (*outputs)->data[n]; // The output image
     170        psImage *error = NULL; // The error image
     171        psImage *outError = NULL; // The output error image
     172        if (errors) {
     173            error = errors->data[n];
     174            outError = (*outErrors)->data[n];
     175        }
     176        float offset = 0.0;             // Relative offset
     177        float scale = 1.0;              // Relative scale
     178        if (offsets) {
     179            offset = offsets->data.F32[n];
     180        }
     181        if (scales) {
     182            scale = scales->data.F32[n];
     183        }
     184
     185        // Mask
     186        psImage *mask = NULL;
     187        if (masks != NULL) {
     188            mask = masks->data[n];
     189        }
     190
     191        // Iterate over the output image pixels
     192        for (int y = 0; y < outny; y++) {
     193            for (int x = 0; x < outnx; x++) {
     194                // Only transform those pixels requested
     195                if (!region || (region && region->data.U8[y][x])) {
     196                    // Transform!
     197                    sky->x = (double)x + 0.5;
     198                    sky->y = (double)y + 0.5;
     199                    (void)psPlaneTransformApply(detector, map, sky);
     200
     201                    // Change PS_INTERPOLATE_BILINEAR to best available technique.
     202                    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
     203                                                                              detector->y, mask, 1, NAN,
     204                                                                              PS_INTERPOLATE_BILINEAR);
     205                    if (error) {
     206                        // Error is actually the variance
     207                        outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
     208                                                                                                detector->x,
     209                                                                                                detector->y,
     210                                                                                                mask, 1, NAN);
     211                    }
     212
     213                    outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
     214                    outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
     215
     216                } // Pixels of interest
     217
     218            }
     219        } // Iterating over output pixels
    220220
    221221    } // Iterating over images

  • trunk/stac/src/stacWrite.c

    r5743 r5745  
    44#include "stac.h"
    55
    6 bool stacWriteMap(const char *mapName,  // Filename to write to
    7                   psPlaneTransform *map // Map to write
     6bool stacWriteMap(const char *mapName,  // Filename to write to
     7                  psPlaneTransform *map // Map to write
    88    )
    99{
    … …  
    1212    FILE *mapFile = fopen(mapName, "w");
    1313    if (!mapFile) {
    14         fprintf(stderr, "Unable to open map file: %s\n", mapName);
    15         return false;
     14        fprintf(stderr, "Unable to open map file: %s\n", mapName);
     15        return false;
    1616    }
    1717
    18     psPolynomial2D *xMap = map->x;      // x transform
    19     psPolynomial2D *yMap = map->y;      // y transform
     18    psPolynomial2D *xMap = map->x;      // x transform
     19    psPolynomial2D *yMap = map->y;      // y transform
    2020
    2121    // A crucial limitation of the current system --- the order of each polynomial must be the same
    2222    assert(xMap->nX == xMap->nY && yMap->nX == yMap->nY && xMap->nX == yMap->nX);
    23     int order = xMap->nX;       // The polynomial order
     23    int order = xMap->nX;       // The polynomial order
    2424    fprintf(mapFile, "%d\n", order);
    25    
     25
    2626    // x coefficients
    27     for (int k = 0; k < order + 1; k++) {
    28         for (int j = 0; j < k + 1; j++) {
    29             int i = k - j;
    30             if (xMap->mask[i][j]) {
    31                 fprintf(mapFile, "0.0 ");
    32             } else {
    33                 fprintf(mapFile, "%g ", xMap->coeff[i][j]);
    34             }
    35         }
     27    for (int k = 0; k <= order; k++) {
     28        for (int j = 0; j <= k; j++) {
     29            int i = k - j;
     30            if (xMap->mask[i][j]) {
     31                fprintf(mapFile, "0.0 ");
     32            } else {
     33                fprintf(mapFile, "%g ", xMap->coeff[i][j]);
     34            }
     35        }
    3636    }
    3737    fprintf(mapFile, "\n");
    3838
    3939    // y coefficients
    40     for (int k = 0; k < order + 1; k++) {
    41         for (int j = 0; j < k + 1; j++) {
    42             int i = k - j;
    43             if (yMap->mask[i][j]) {
    44                 fprintf(mapFile, "0.0 ");
    45             } else {
    46                 fprintf(mapFile, "%g ", yMap->coeff[i][j]);
    47             }
    48         }
     40    for (int k = 0; k <= order; k++) {
     41        for (int j = 0; j <= k; j++) {
     42            int i = k - j;
     43            if (yMap->mask[i][j]) {
     44                fprintf(mapFile, "0.0 ");
     45            } else {
     46                fprintf(mapFile, "%g ", yMap->coeff[i][j]);
     47            }
     48        }
    4949    }
    5050    fprintf(mapFile, "\n");
    51    
     51
    5252    fclose(mapFile);
    5353
    … …  
    5757
    5858bool stacWriteMaps(const psArray *names, // Filenames of the input images (will add ".map")
    59                    const psArray *maps  // Maps to write
     59                   const psArray *maps  // Maps to write
    6060    )
    6161{
    … …  
    6565
    6666    for (int i = 0; i < names->n; i++) {
    67         char mapName[MAXCHAR];          // Filename of error image
    68         sprintf(mapName, "%s.map", names->data[i]);
    69         if (!stacWriteMap(mapName, maps->data[i])) {
    70             return false;
    71         }
     67        char mapName[MAXCHAR];          // Filename of error image
     68        sprintf(mapName, "%s.map", (const char*)names->data[i]);
     69        if (!stacWriteMap(mapName, maps->data[i])) {
     70            return false;
     71        }
    7272    }
    7373

  • trunk/stac/src/sum.c

    r5743 r5745  
    88int main(int argc, char *argv[])
    99{
    10     const char *programName = argv[0];  // Program name
    1110    const char *outputName = argv[1];   // Output file name
    1211    psArray *inputNames = psArrayAlloc(argc-2);
Note: See TracChangeset for help on using the changeset viewer.