IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 31, 2005, 5:07:12 PM (21 years ago)
Author:
Paul Price
Message:

Working very nicely now.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/stac/src/stacRead.c

    r3375 r3610  
    33#include "pslib.h"
    44#include "stac.h"
     5
     6#define BUFFER 100                      // Size of buffer for incrementally reading coordinates
    57
    68psArray *stacReadImages(stacConfig *config)
     
    4042}
    4143
     44psArray *stacReadCoords(const char *filename)
     45{
     46    FILE *file = fopen(filename, "r");
     47    if (file == NULL) {
     48        psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
     49        return NULL;
     50    }
     51
     52    psTrace("stac.read.coords", 5, "Reading coordinate file, %s\n", filename);
     53
     54    psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
     55    float x, y;                         // Coordinates to read
     56    int num = 0;                        // Number of coordinates read
     57    while (fscanf(file, "%f %f\n", &x, &y) == 2) {
     58        psPlane *coord = psPlaneAlloc();// A coordinate
     59        coord->x = x;
     60        coord->y = y;
     61        coords->data[num] = coord;
     62        num++;
     63        if (num % BUFFER) {
     64            coords = psArrayRealloc(coords, num + BUFFER);
     65        }
     66    }
     67    coords->n = num;
     68
     69    psTrace("stac.read.coords", 5, "%d coordinates read.\n", num);
     70
     71    fclose(file);
     72    return coords;
     73}
     74
     75
     76psPlaneTransform *stacReadMap(const char *filename)
     77{
     78    FILE *mapfp = fopen(filename, "r");
     79    if (mapfp == NULL) {
     80        psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
     81        return NULL;
     82    }
     83    // Read the file
     84    psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);
     85   
     86    // Format is now:
     87    // order
     88    // x coefficients
     89    // y coefficients
     90    //
     91    // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
     92    // and the coefficients are read by the following pseudo-code:
     93    //
     94    //  for (int k = 0; k < order + 1; k++)
     95    //      for (int j = 0; j < k; j++)
     96    //          int i = k - j;
     97    //              read coefficient of x^i y^j
     98    //
     99    // This produces the ordering:
     100    // 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
     101    //
     102    // This is, of course, for third order polynomials.
     103    // For lower orders, the list is truncated at the appropriate level.
     104
     105    int order = 0;                      // Polynomial order
     106    if (fscanf(mapfp, "%d", &order) != 1) {
     107        psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     108        fclose(mapfp);
     109        return NULL;
     110    }
     111   
     112    psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
     113   
     114    psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
     115    // Set coefficient masks
     116    for (int i = 0; i < order + 1; i++) {
     117        for (int j = 0; j < order + 1; j++) {
     118            if (i + j > order + 1) {
     119                map->x->mask[i][j] = 1;
     120                map->y->mask[i][j] = 1;
     121            } else {
     122                map->x->mask[i][j] = 0;
     123                map->y->mask[i][j] = 0;
     124            }
     125        }
     126    }
     127   
     128    // Read x coefficients
     129    psTrace("stac.read.map", 7, "x' = \n");
     130    for (int k = 0; k < order + 1; k++) {
     131        for (int j = 0; j < k + 1; j++) {
     132            int i = k - j;
     133            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
     134                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     135                fclose(mapfp);
     136                psFree(map);
     137                return NULL;
     138            }
     139            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
     140        }
     141    }
     142    // Read y coefficients
     143    psTrace("stac.read.maps", 7, "y' = \n");
     144    for (int k = 0; k < order + 1; k++) {
     145        for (int j = 0; j < k + 1; j++) {
     146            int i = k - j;
     147            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
     148                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     149                fclose(mapfp);
     150                psFree(map);
     151                return NULL;
     152            }
     153            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
     154        }
     155    }
     156   
     157    fclose(mapfp);
     158
     159    return map;
     160}
     161   
     162
    42163
    43164psArray *stacReadMaps(stacConfig *config)
     
    54175            exit(EXIT_FAILURE);
    55176        }
    56         // Open the file
     177        // Read the file
    57178        sprintf(mapfile,"%s.map",filenames->data[i]);
    58         FILE *mapfp = fopen(mapfile,"r");
    59         if (mapfp == NULL) {
    60             psLogMsg("stac.read.maps", PS_LOG_ERROR, "Cannot open map file, %s\n",mapfile);
    61             exit(EXIT_FAILURE);
     179        maps->data[i] = stacReadMap(mapfile);
     180        if (maps->data[i] == NULL) {
     181            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
    62182        }
    63         // Read the file
    64         psTrace("stac.read.maps",5,"Reading map file %s....\n",mapfile);
    65 
    66         // Format is now:
    67         // order
    68         // x coefficients
    69         // y coefficients
    70         //
    71         // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
    72         // and the coefficients are read by the following pseudo-code:
    73         //
    74         //      for (int k = 0; k < order + 1; k++)
    75         //          for (int j = 0; j < k; j++)
    76         //              int i = k - j;
    77         //              read coefficient of x^i y^j
    78         //
    79         // This produces the ordering:
    80         // 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
    81         //
    82         // This is, of course, for third order polynomials.
    83         // For lower orders, the list is truncated at the appropriate level.
    84 
    85         int order = 0;                  // Polynomial order
    86         if (fscanf(mapfp, "%d", &order) != 1) {
    87             psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
    88             exit(EXIT_FAILURE);
    89         }
    90 
    91         psTrace("stac.read.maps", 5, "Polynomial order: %d\n", order);
    92 
    93         psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
    94         // Set coefficient masks
    95         for (int i = 0; i < order + 1; i++) {
    96             for (int j = 0; j < order + 1; j++) {
    97                 if (i + j > order + 1) {
    98                     map->x->mask[i][j] = 1;
    99                     map->y->mask[i][j] = 1;
    100                 } else {
    101                     map->x->mask[i][j] = 0;
    102                     map->y->mask[i][j] = 0;
    103                 }
    104             }
    105         }
    106 
    107         // Read x coefficients
    108         psTrace("stac.read.maps", 7, "x' = \n");
    109         for (int k = 0; k < order + 1; k++) {
    110             for (int j = 0; j < k + 1; j++) {
    111                 int i = k - j;
    112                 if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
    113                     psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
    114                     exit(EXIT_FAILURE);
    115                 }
    116                 psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
    117             }
    118         }
    119         // Read y coefficients
    120         psTrace("stac.read.maps", 7, "y' = \n");
    121         for (int k = 0; k < order + 1; k++) {
    122             for (int j = 0; j < k + 1; j++) {
    123                 int i = k - j;
    124                 if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
    125                     psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map file %s\n", mapfile);
    126                     exit(EXIT_FAILURE);
    127                 }
    128                 psTrace("stac.read.maps", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
    129             }
    130         }
    131 
    132         fclose(mapfp);
    133 
    134         // Plug it into the array
    135         maps->data[i] = map;
    136183
    137184    }
Note: See TracChangeset for help on using the changeset viewer.