IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 7, 2005, 2:57:14 PM (20 years ago)
Author:
Paul Price
Message:

Importing from PAP cvs tree again. This will now be the working tree.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/stacRead.c

    r3681 r5743  
    11#include <stdio.h>
    22#include <string.h>
     3#include <assert.h>
    34#include "pslib.h"
    45#include "stac.h"
    56
    6 #define BUFFER 100                      // Size of buffer for incrementally reading coordinates
    7 
    8 psArray *stacReadImages(psArray *filenames // The file names of the images
     7#define BUFFER 100                      // Size of buffer for incrementally reading coordinates
     8
     9psArray *stacReadImages(psArray **headers, // The image headers, to be returned
     10    psArray *filenames // The file names of the images
    911    )
    1012{
    11     int nFiles = filenames->n;          // The number of input files
     13    int nFiles = filenames->n;          // The number of input files
    1214    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
    13     psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
     15    assert(!headers || ! *headers || (*headers)->n == nFiles);
     16    if (headers && ! *headers) {
     17        *headers = psArrayAlloc(nFiles);
     18    }
     19
     20    psRegion imageRegion = {0, 0, 0, 0}; // Region of image to read
    1421
    1522    psTrace("stac.read.images", 1, "Reading input images....\n");
    1623    for (int i = 0; i < nFiles; i++) {
    17         psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
    18         psFits *imageFile = psFitsAlloc(filenames->data[i]);
    19         // We only read PHUs --- not mucking around with extensions for now
    20         psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
    21         if (image == NULL) {
    22             psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
    23             exit(EXIT_FAILURE);
    24         }
    25         psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
    26                 image->numRows);
    27         // Convert to 32-bit floating point, in necessary
    28         if (image->type.type != PS_TYPE_F32) {
    29             psTrace("stac.read.images", 3, "Converting %s to floating point in memory....",
    30                     filenames->data[i]);
    31             images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
    32             psFree(image);
    33         } else {
    34             int numNaN = psImageClipNaN(image, FLT_MAX);
    35             if (numNaN) {
    36                 psTrace("stac.read.images", 5, "Clipped %d NaN pixels.\n", numNaN);
    37             }
    38             images->data[i] = image;
    39         }
    40         psFree(imageFile);
     24        psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
     25        psFits *imageFile = psFitsOpen(filenames->data[i], "r");
     26        // We only read PHUs --- not mucking around with extensions for now
     27        if (headers) {
     28            (*headers)->data[i] = psFitsReadHeader(NULL, imageFile);
     29        }
     30        psImage *image = psFitsReadImage(NULL, imageFile, imageRegion, 0);
     31        if (image == NULL) {
     32            psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
     33            exit(EXIT_FAILURE);
     34        }
     35        psFitsClose(imageFile);
     36        psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
     37                image->numRows);
     38        // Convert to 32-bit floating point, in necessary
     39        if (image->type.type != PS_TYPE_F32) {
     40            psTrace("stac.read.images", 3, "Converting %s to floating point in memory....\n",
     41                    filenames->data[i]);
     42            images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
     43            psFree(image);
     44        }
     45        int numNaN = psImageClipNaN(image, 0.0);
     46        if (numNaN) {
     47            psTrace("stac.read.images", 5, "Clipped %d NaN pixels to zero.\n", numNaN);
     48        }
     49        images->data[i] = image;
    4150    }
    4251    psTrace("stac.read.images",1,"%d input images read.\n",nFiles);
    43     psFree(imageRegion);
    4452
    4553    return images;
    … …  
    5058    FILE *file = fopen(filename, "r");
    5159    if (file == NULL) {
    52         psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
    53         return NULL;
     60        psLogMsg("stac.read.coords", PS_LOG_ERROR, "Cannot open coordinate file, %s\n", filename);
     61        return NULL;
    5462    }
    5563
    … …  
    5765
    5866    psArray *coords = psArrayAlloc(BUFFER); // The array of coordinates to be returned
    59     float x, y;                         // Coordinates to read
    60     int num = 0;                        // Number of coordinates read
     67    float x, y;                         // Coordinates to read
     68    int num = 0;                        // Number of coordinates read
    6169    while (fscanf(file, "%f %f\n", &x, &y) == 2) {
    62         psPlane *coord = psPlaneAlloc();// A coordinate
    63         coord->x = x;
    64         coord->y = y;
    65         coords->data[num] = coord;
    66         num++;
    67         if (num % BUFFER) {
    68             coords = psArrayRealloc(coords, num + BUFFER);
    69         }
     70        psPlane *coord = psPlaneAlloc();// A coordinate
     71        coord->x = x;
     72        coord->y = y;
     73        coords->data[num] = coord;
     74        num++;
     75        if (num % BUFFER) {
     76            coords = psArrayRealloc(coords, num + BUFFER);
     77        }
    7078    }
    7179    coords->n = num;
    … …  
    8290    FILE *mapfp = fopen(filename, "r");
    8391    if (mapfp == NULL) {
    84         psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
    85         return NULL;
     92        psLogMsg("stac.read.map", PS_LOG_ERROR, "Cannot open map file, %s\n", filename);
     93        return NULL;
    8694    }
    8795    // Read the file
    8896    psTrace("stac.read.map", 5, "Reading map file %s....\n", filename);
    89    
     97
    9098    // Format is now:
    9199    // order
    92100    // x coefficients
    93101    // y coefficients
    94     // 
     102    //
    95103    // where the order is 1 for linear, 2 for quadratic, 3 for cubic.
    96104    // and the coefficients are read by the following pseudo-code:
    97     // 
    98     //  for (int k = 0; k < order + 1; k++)
    99     //      for (int j = 0; j < k; j++)
    100     //          int i = k - j;
     105    //
     106    //  for (int k = 0; k < order + 1; k++)
     107    //      for (int j = 0; j < k; j++)
     108    //          int i = k - j;
    101109    //              read coefficient of x^i y^j
    102     // 
     110    //
    103111    // This produces the ordering:
    104112    // 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
    105     // 
     113    //
    106114    // This is, of course, for third order polynomials.
    107115    // For lower orders, the list is truncated at the appropriate level.
    108116
    109     int order = 0;                      // Polynomial order
     117    int order = 0;                      // Polynomial order
    110118    if (fscanf(mapfp, "%d", &order) != 1) {
    111         psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
    112         fclose(mapfp);
    113         return NULL;
    114     }
    115    
     119        psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     120        fclose(mapfp);
     121        return NULL;
     122    }
     123
    116124    psTrace("stac.read.map", 5, "Polynomial order: %d\n", order);
    117    
     125
    118126    psPlaneTransform *map = psPlaneTransformAlloc(order + 1, order + 1); // The transformation
    119127    // Set coefficient masks
    120128    for (int i = 0; i < order + 1; i++) {
    121         for (int j = 0; j < order + 1; j++) {
    122             if (i + j > order + 1) {
    123                 map->x->mask[i][j] = 1;
    124                 map->y->mask[i][j] = 1;
    125             } else {
    126                 map->x->mask[i][j] = 0;
    127                 map->y->mask[i][j] = 0;
    128             }
    129         }
    130     }
    131    
     129        for (int j = 0; j < order + 1; j++) {
     130            if (i + j > order + 1) {
     131                map->x->mask[i][j] = 1;
     132                map->y->mask[i][j] = 1;
     133            } else {
     134                map->x->mask[i][j] = 0;
     135                map->y->mask[i][j] = 0;
     136            }
     137        }
     138    }
     139
    132140    // Read x coefficients
    133141    psTrace("stac.read.map", 7, "x' = \n");
    134142    for (int k = 0; k < order + 1; k++) {
    135         for (int j = 0; j < k + 1; j++) {
    136             int i = k - j;
    137             if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
    138                 psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
    139                 fclose(mapfp);
    140                 psFree(map);
    141                 return NULL;
    142             }
    143             psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
    144         }
     143        for (int j = 0; j < k + 1; j++) {
     144            int i = k - j;
     145            if (fscanf(mapfp, "%lf", &map->x->coeff[i][j]) != 1) {
     146                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     147                fclose(mapfp);
     148                psFree(map);
     149                return NULL;
     150            }
     151            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->x->coeff[i][j], i, j);
     152        }
    145153    }
    146154    // Read y coefficients
    147155    psTrace("stac.read.maps", 7, "y' = \n");
    148156    for (int k = 0; k < order + 1; k++) {
    149         for (int j = 0; j < k + 1; j++) {
    150             int i = k - j;
    151             if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
    152                 psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
    153                 fclose(mapfp);
    154                 psFree(map);
    155                 return NULL;
    156             }
    157             psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
    158         }
    159     }
    160    
     157        for (int j = 0; j < k + 1; j++) {
     158            int i = k - j;
     159            if (fscanf(mapfp, "%lf", &map->y->coeff[i][j]) != 1) {
     160                psLogMsg("stac.read.map", PS_LOG_ERROR, "Unable to read map file %s\n", filename);
     161                fclose(mapfp);
     162                psFree(map);
     163                return NULL;
     164            }
     165            psTrace("stac.read.map", 7, "      %f x^%d y^%d\n", map->y->coeff[i][j], i, j);
     166        }
     167    }
     168
    161169    fclose(mapfp);
    162170
    163171    return map;
    164172}
    165    
     173
    166174
    167175
    … …  
    169177    )
    170178{
    171     int nFiles = filenames->n;          // The number of input files
     179    int nFiles = filenames->n;          // The number of input files
    172180    psArray *maps = psArrayAlloc(nFiles); // The maps, to be returned
    173     char mapfile[MAXCHAR];              // Filename of map
    174    
     181    char mapfile[MAXCHAR];              // Filename of map
     182
    175183    psTrace("stac.read.maps",1,"Reading maps....\n");
    176184    for (int i = 0; i < nFiles; i++) {
    177         if (strlen(filenames->data[i]) > MAXCHAR - 4) {
    178             psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
    179             exit(EXIT_FAILURE);
    180         }
    181         // Read the file
    182         sprintf(mapfile,"%s.map",filenames->data[i]);
    183         maps->data[i] = stacReadMap(mapfile);
    184         if (maps->data[i] == NULL) {
    185             psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
    186         }
     185        if (strlen(filenames->data[i]) > MAXCHAR - 4) {
     186            psLogMsg("stac.read.maps",PS_LOG_ERROR,"Filename %s is too long.\n",filenames->data[i]);
     187            exit(EXIT_FAILURE);
     188        }
     189        // Read the file
     190        sprintf(mapfile,"%s.map",filenames->data[i]);
     191        maps->data[i] = stacReadMap(mapfile);
     192        if (maps->data[i] == NULL) {
     193            psLogMsg("stac.read.maps", PS_LOG_ERROR, "Unable to read map: %s\n", mapfile);
     194        }
    187195
    188196    }
Note: See TracChangeset for help on using the changeset viewer.