IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2897


Ignore:
Timestamp:
Jan 4, 2005, 3:55:40 PM (22 years ago)
Author:
Paul Price
Message:

Updated to handle higher polynomial orders than linear, and to calculate the output size itself

Location:
trunk/stac/src
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/Makefile

    r2764 r2897  
    55LDFLAGS += $(PSLIB)
    66
    7 OBJECTS = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
    8         stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o
     7STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
     8        stacCombine.o stacInvertMaps.o stacRejection.o psPlaneTransform.o stacAreaOfInterest.o stacSize.o
    99
    10 TARGET = stac
     10FITTRANS = fitTrans.o fitTransConfig.o fitTransRead.o stacCheckMemory.o fitTransSolve.o psPlaneTransform.o
    1111
    1212.PHONY: tags clean empty test profile optimise
    … …  
    1515                $(CC) -c $(CFLAGS) $(OPTFLAGS) $<
    1616
    17 $(TARGET):      $(OBJECTS)
    18                 $(CC) $(CFLAGS) -o $@ $(OBJECTS) $(LDFLAGS) $(OPTFLAGS)
     17stac:           $(STAC)
     18                $(CC) $(CFLAGS) -o $@ $(STAC) $(LDFLAGS) $(OPTFLAGS)
    1919
    2020calcGrad:       calcGrad.o
    … …  
    2727                $(CC) $(CFLAGS) -o $@ getPixels.o $(LDFLAGS) $(OPTFLAGS)
    2828
     29fitTrans:       $(FITTRANS)
     30                $(CC) $(CFLAGS) -o $@ $(FITTRANS) $(LDFLAGS) $(OPTFLAGS)
    2931
    3032clean:

  • trunk/stac/src/stac.c

    r2783 r2897  
    3737    psTraceSetLevel("stac.time",10);
    3838    psTraceSetLevel("stac.area",10);
     39    psTraceSetLevel("stac.size",10);
    3940
    4041    // Set logging level
    … …  
    5253    psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime);
    5354
     55    // Get size
     56    stacSize(config, inputs, maps);
     57
    5458    // Invert maps
    55     psArray *inverseMaps = stacInvertMaps(maps);
     59    psArray *inverseMaps = stacInvertMaps(maps, config);
    5660
    5761    // Generate errors

  • trunk/stac/src/stac.h

    r2783 r2897  
    2222    const char *output;                 // Output file name
    2323    int outnx, outny;                   // Size of output image
     24    float xOffset, yOffset;             // Offset to get lower-left corner at 0,0
    2425    float saturated;                    // Saturation level
    2526    float bad;                          // Bad level
    … …  
    127128
    128129// Invert an array of maps
    129 psArray *stacInvertMaps(const psArray *maps // Array of maps to invert
     130psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
     131                        const stacConfig *config // Configuration
    130132    );
    131133
    … …  
    168170/************************************************************************************************************/
    169171
     172// stacSize.c
     173
     174// Calculate the size of the output image
     175bool stacSize(stacConfig *config,       // Configuration, containing the output size
     176              const psArray *images,    // Input images
     177              psArray *maps             // Transformation maps
     178    );
     179
     180
     181/************************************************************************************************************/
     182
    170183#endif

  • trunk/stac/src/stacInvertMaps.c

    r2500 r2897  
    55
    66#define MAX(x,y) ((x) > (y) ? (x) : (y))
    7 
    8 psArray *stacInvertMaps(const psArray *maps // Array of maps to invert
     7#define NUM_GRID 20
     8
     9psArray *stacInvertMaps(const psArray *maps, // Array of maps to invert
     10                        const stacConfig *config // Configuration
    911    )
    1012{
    … …  
    1416    psTrace("stac.invertMaps", 1, "Inverting maps....\n");
    1517
    16     for (int i = 0; i < nMaps; i++) {
    17         // Can't handle higher order than linear yet
    18         if (((psPlaneTransform*)maps->data[i])->x->nX != 2 ||
    19             ((psPlaneTransform*)maps->data[i])->x->nY != 2 ||
    20             ((psPlaneTransform*)maps->data[i])->y->nX != 2 ||
    21             ((psPlaneTransform*)maps->data[i])->y->nY != 2) {
    22             psError("stac.invertMaps",
    23                     "STAC cannot currently support orders other than linear.\n");
    24             psFree(inverted);
    25             return NULL;
    26         }
    27 
    28         psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
    29         psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
    30 
    31         // Now, simply do a 2x2 matrix inversion
    32 
    33         double a = oldMap->x->coeff[1][0];
    34         double b = oldMap->x->coeff[0][1];
    35         double c = oldMap->y->coeff[1][0];
    36         double d = oldMap->y->coeff[0][1];
    37         double e = oldMap->x->coeff[0][0];
    38         double f = oldMap->y->coeff[0][0];
    39 
    40         double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
    41 
    42         // Not entirely sure why this works, but it appears to do so.......................................
    43         newMap->x->coeff[1][0] = invDet * a;
    44         newMap->x->coeff[0][1] = - invDet * b;
    45         newMap->y->coeff[1][0] = - invDet * c;
    46         newMap->y->coeff[0][1] = invDet * d;
    47 
    48         newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
    49         newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
     18    // Coordinates for the transformations
     19    psPlane *inCoord = psAlloc(sizeof(psPlane));
     20    psPlane *outCoord = psAlloc(sizeof(psPlane));
     21
     22    for (int mapNum = 0; mapNum < nMaps; mapNum++) {
     23
     24        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
     25        // Check input
     26        if (oldMap->x->nX != oldMap->x->nY) {
     27            psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
     28            return NULL;
     29        }
     30        if (oldMap->y->nX != oldMap->y->nY) {
     31            psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
     32            return NULL;
     33        }
     34        if (oldMap->x->nX != oldMap->y->nX) {
     35            psError("stac.invertMaps", "Polynomial order in x and y don't match for map %d\n", mapNum);
     36            return NULL;
     37        }
     38        int order = oldMap->x->nX;      // Polynomial order
     39        psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
     40        psPlaneTransform *newMap = psPlaneTransformAlloc(order, order); // Inverted map
     41
     42        // Create fake polynomial to use in evaluation
     43        psDPolynomial2D *fakePoly = psDPolynomial2DAlloc(order, order, PS_POLYNOMIAL_ORD);
     44        for (int i = 0; i < order; i++) {
     45            for (int j = 0; j < order; j++) {
     46                fakePoly->coeff[i][j] = 1.0; // Set all coeffecients to 1
     47                fakePoly->mask[i][j] = 1; // Mask all coefficients; unmask to evaluate
     48            }
     49        }
     50
     51        // A grid of xin,yin --> xout,yout
     52        psVector *xIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     53        psVector *yIn = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     54        psVector *xOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     55        psVector *yOut = psVectorAlloc(NUM_GRID * NUM_GRID, PS_TYPE_F32);
     56
     57        // Create grid of points
     58        for (int yint = 0; yint < NUM_GRID; yint++) {
     59            inCoord->y = (float)(yint * config->outny) / (float)(NUM_GRID - 1);
     60            for (int xint = 0; xint < NUM_GRID; xint++) {
     61                inCoord->x = (float)(xint * config->outnx) / (float)(NUM_GRID - 1);
     62
     63                (void)psPlaneTransformApply(outCoord, oldMap, inCoord);
     64
     65                xOut->data.F32[yint*NUM_GRID + xint] = inCoord->x;
     66                yOut->data.F32[yint*NUM_GRID + xint] = inCoord->y;
     67                xIn->data.F32[yint*NUM_GRID + xint] = outCoord->x;
     68                yIn->data.F32[yint*NUM_GRID + xint] = outCoord->y;
     69            }
     70        }
     71
     72        // Initialise the matrix and vectors
     73        int nCoeff = order * (order + 1) / 2; // Number of polynomial coefficients
     74        psImage *matrix = psImageAlloc(nCoeff, nCoeff, PS_TYPE_F64); // Matrix for solution
     75        psVector *xVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in x
     76        psVector *yVector = psVectorAlloc(nCoeff, PS_TYPE_F64); // Vector for solution in y
     77        for (int i = 0; i < nCoeff; i++) {
     78            for (int j = 0; j < nCoeff; j++) {
     79                matrix->data.F64[i][j] = 0.0;
     80            }
     81            xVector->data.F64[i] = 0.0;
     82            yVector->data.F64[i] = 0.0;
     83        }
     84
     85        // Iterate over the grid points
     86        for (int g = 0; g < NUM_GRID*NUM_GRID; g++) {
     87
     88            // Iterate over the polynomial coefficients, accumulating the matrix and vectors
     89            for (int i = 0, ijIndex = 0; i < order; i++) {
     90                for (int j = 0; j < order - i; j++, ijIndex++) {
     91
     92                    fakePoly->mask[i][j] = 0;
     93                    double ijPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
     94                    fakePoly->mask[i][j] = 1;
     95
     96                    for (int m = 0, mnIndex = 0; m < order; m++) {
     97                        for (int n = 0; n < order - m; n++, mnIndex++) {
     98                           
     99                            fakePoly->mask[m][n] = 0;
     100                            double mnPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
     101                            fakePoly->mask[m][n] = 1;
     102
     103                            matrix->data.F64[ijIndex][mnIndex] += ijPoly * mnPoly;
     104                        }
     105                    }
     106                   
     107                    xVector->data.F64[ijIndex] += ijPoly * (double)xOut->data.F32[g];
     108                    yVector->data.F64[ijIndex] += ijPoly * (double)yOut->data.F32[g];
     109                }
     110            } // Iterating over coefficients
     111        } // Iterating over grid points
     112
     113        // Solution via LU Decomposition
     114        psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     115        psImage *luMatrix = psMatrixLUD(NULL, permutation, matrix); // LU decomposed matrix
     116        psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
     117        psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     118
     119        // Stuff coefficients into transformation
     120        for (int i = 0, ijIndex = 0; i < order; i++) {
     121            for (int j = 0; j < order - i; j++, ijIndex++) {
     122                newMap->x->coeff[i][j] = xSolution->data.F64[ijIndex];
     123                newMap->y->coeff[i][j] = ySolution->data.F64[ijIndex];
     124            }
     125        }
     126        inverted->data[mapNum] = newMap;
     127
     128#ifdef TESTING
     129        // Print x coefficients
     130        psTrace("stac.invertMaps", 7, "x' = \n");
     131        for (int i = 0; i < order; i++) {
     132            for (int j = 0; j < order - i; j++) {
     133                psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->x->coeff[i][j], i, j);
     134            }
     135        }
     136        // Print y coefficients
     137        psTrace("stac.invertMaps", 7, "y' = \n");
     138        for (int i = 0; i < order; i++) {
     139            for (int j = 0; j < order - i; j++) {
     140                psTrace("stac.invertMaps", 7, "      %f x^%d y^%d\n", newMap->y->coeff[i][j], i, j);
     141            }
     142        }
     143#endif
    50144
    51145#ifdef TESTING
    … …  
    64158#endif
    65159
    66         inverted->data[i] = newMap;     // Stuff into the array
     160        psFree(permutation);
     161        psFree(luMatrix);
     162        psFree(matrix);
     163        psFree(xVector);
     164        psFree(yVector);
     165        psFree(xSolution);
     166        psFree(ySolution);
     167        psFree(fakePoly);
     168        psFree(xIn);
     169        psFree(yIn);
     170        psFree(xOut);
     171        psFree(yOut);
    67172    }
     173
     174    psFree(inCoord);
     175    psFree(outCoord);
     176               
     177
     178#if 0
     179        // Can't handle higher order than linear yet
     180        if (((psPlaneTransform*)maps->data[i])->x->nX != 2 ||
     181            ((psPlaneTransform*)maps->data[i])->x->nY != 2 ||
     182            ((psPlaneTransform*)maps->data[i])->y->nX != 2 ||
     183            ((psPlaneTransform*)maps->data[i])->y->nY != 2) {
     184            psError("stac.invertMaps",
     185                    "STAC cannot currently support orders other than linear.\n");
     186            psFree(inverted);
     187            return NULL;
     188        }
     189
     190        psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
     191        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
     192
     193        // Now, simply do a 2x2 matrix inversion
     194
     195        double a = oldMap->x->coeff[1][0];
     196        double b = oldMap->x->coeff[0][1];
     197        double c = oldMap->y->coeff[1][0];
     198        double d = oldMap->y->coeff[0][1];
     199        double e = oldMap->x->coeff[0][0];
     200        double f = oldMap->y->coeff[0][0];
     201
     202        double invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
     203
     204        // Not entirely sure why this works, but it appears to do so.......................................
     205        newMap->x->coeff[1][0] = invDet * a;
     206        newMap->x->coeff[0][1] = - invDet * b;
     207        newMap->y->coeff[1][0] = - invDet * c;
     208        newMap->y->coeff[0][1] = invDet * d;
     209
     210        newMap->x->coeff[0][0] = - invDet * (d * e + c * f);
     211        newMap->y->coeff[0][0] = - invDet * (b * e + a * f);
     212#endif
     213
    68214
    69215    return inverted;

  • trunk/stac/src/stacRead.c

    r2764 r2897  
    33#include "pslib.h"
    44#include "stac.h"
    5 
    6 #define ORDER 2
    75
    86psArray *stacReadImages(stacConfig *config)
    … …  
    5452        // Read the file
    5553        psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
    56         psPlaneTransform *map = psPlaneTransformAlloc(ORDER,ORDER);
    5754
    58         // Format is as following:
    59         // A0  A1
    60         // B00 B01
    61         // B10 B11
    62         //
    63         // where:
    64         // x_sky = A0 + B00*x + B10*y
    65         // y_sky = A1 + B01*x + B11*y
     55        // Format is now:
     56        // order
     57        // x coefficients
     58        // y coefficients
     59        //
     60        // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
     61        // and the coefficients are read by the following pseudo-code:
     62        //
     63        //      for (int k = 0; k < order + 1; k++)
     64        //          for (int j = 0; j < k; j++)
     65        //              int i = k - j;
     66        //              read coefficient of x^i y^j
     67        //
     68        // This produces the ordering:
     69        // x^0 y^0, x^1 y^0, x^0 y^1, x^2 y^0, x^1 y^1, x^0 y^2, x^3 y^0, x^2 y^1, x^1 y^2, x^0 y^3
     70        //
     71        // This is, of course, for third order polynomials.
     72        // For lower orders, the list is truncated at the appropriate level.
    6673
    67         double a0, a1, b00, b10, b01, b11;
    68 
    69         if (fscanf(mapfp,"%lf %lf %lf %lf %lf %lf", &a0, &a1, &b00, &b01, &b10, &b11) != 6) {
    70             psLogMsg("stac.read.maps", PS_LOG_ERROR, "Couldn't read all 6 transformation values from %s\n",
    71                      mapfile);
     74        int order = 0;                  // Polynomial order
     75        if (fscanf(mapfp, "%d", &order) != 1) {
     76            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
    7277            exit(EXIT_FAILURE);
    7378        }
    74         fclose(mapfp);
    75         map->x->coeff[0][0] = a0;
    76         map->x->coeff[1][0] = b00;
    77         map->x->coeff[0][1] = b10;
    78         map->y->coeff[0][0] = a1;
    79         map->y->coeff[1][0] = b01;
    80         map->y->coeff[0][1] = b11;
    8179
    82         // Mask out the cross-terms
    83         map->x->mask[1][1] = 1;
    84         map->y->mask[1][1] = 1;
    85         // Set the cross-terms to zero, just in case...
    86         map->x->coeff[1][1] = 0.0;
    87         map->y->coeff[1][1] = 0.0;
     80        psTrace("stac.read.maps", 5, "Polynomial order: %d\n", order);
     81
     82        psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
     83        // Set coefficient masks
     84        for (int i = 0; i < order + 1; i++) {
     85            for (int j = 0; j < order + 1; j++) {
     86                if (i + j > order + 1) {
     87                    map->x->mask[i][j] = 1;
     88                    map->y->mask[i][j] = 1;
     89                } else {
     90                    map->x->mask[i][j] = 0;
     91                    map->y->mask[i][j] = 0;
     92                }
     93            }
     94        }
     95
     96        // Read x coefficients
     97        psTrace("stac.read.maps", 7, "x' = \n");
     98        for (int k = 0; k < order + 1; k++) {
     99            for (int j = 0; j < k + 1; j++) {
     100                int i = k - j;
     101                if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
     102                    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
     103                    exit(EXIT_FAILURE);
     104                }
     105                psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
     106            }
     107        }
     108        // Read y coefficients
     109        psTrace("stac.read.maps", 7, "y' = \n");
     110        for (int k = 0; k < order + 1; k++) {
     111            for (int j = 0; j < k + 1; j++) {
     112                int i = k - j;
     113                if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
     114                    psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
     115                    exit(EXIT_FAILURE);
     116                }
     117                psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
     118            }
     119        }
     120
     121        fclose(mapfp);
    88122
    89123        // Plug it into the array
    90124        maps->data[i] = map;
    91 
    92 #ifdef TESTING
    93         // Check psLib's answer is the same as I would get with old version
    94         double detx = 256.0, dety = 256.0;
    95         // My code
    96         double myx = a0 + b00*detx + b10*dety;
    97         double myy = a1 + b01*detx + b11*dety;
    98         // psLib
    99         psPlane *det = psAlloc(sizeof(psPlane));
    100         det->x = 256.0;
    101         det->y = 256.0;
    102         psPlane *pslib = psPlaneTransformApply(NULL, map, det);
    103 
    104         psTrace("stac.read.maps.test",0,"me: %lf,%lf\n",myx,myy);
    105         psTrace("stac.read.maps.test",0,"ps: %lf,%lf\n",pslib->x,pslib->y);
    106         psFree(det);
    107         psFree(pslib);
    108 #endif
    109125
    110126    }
Note: See TracChangeset for help on using the changeset viewer.