IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3375


Ignore:
Timestamp:
Mar 4, 2005, 11:37:01 AM (21 years ago)
Author:
Paul Price
Message:

Updated to psLib rel5alpha1

Location:
trunk/stac/src
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • trunk/stac/src/Makefile

    r2897 r3375  
    11SHELL = /bin/sh
    22CC = gcc
    3 CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib3/psLib/include -DCRFLUX -DTESTING
    4 PSLIB += -L/home/mithrandir/price/psLib3/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm
     3CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -DCRFLUX -DTESTING
     4PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2
    55LDFLAGS += $(PSLIB)
    66
    77STAC = 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
     8        stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o
    99
    10 FITTRANS = fitTrans.o fitTransConfig.o fitTransRead.o stacCheckMemory.o fitTransSolve.o psPlaneTransform.o
     10SHIFT = shift.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
     11        stacCombine.o stacInvertMaps.o stacSize.o
     12
     13FITTRANS = fitTrans.o fitTransConfig.o fitTransRead.o stacCheckMemory.o fitTransSolve.o
    1114
    1215.PHONY: tags clean empty test profile optimise
     
    1720stac:           $(STAC)
    1821                $(CC) $(CFLAGS) -o $@ $(STAC) $(LDFLAGS) $(OPTFLAGS)
     22
     23shift:          $(SHIFT)
     24                $(CC) $(CFLAGS) -o $@ $(SHIFT) $(LDFLAGS) $(OPTFLAGS)
    1925
    2026calcGrad:       calcGrad.o
  • trunk/stac/src/psPlaneTransform.c

    r2500 r3375  
    1 #include <stdio.h>
    2 #include "pslib.h"
    3 #include "stac.h"
    4 
    5 #define PS_INT_CHECK_NON_NEGATIVE(NAME, RVAL) \
    6 if (NAME < 0) { \
    7     psError("pslib.badValue", \
    8             "Error: %s is less than 0.", #NAME); \
    9     return(RVAL); \
    10 }
    11 
    12 static void planeTransformFree(psPlaneTransform *pt)
    13 {
    14     psFree(pt->x);
    15     psFree(pt->y);
    16 }
    17 
    18 psPlaneTransform* psPlaneTransformAlloc(psS32 n1, psS32 n2)
    19 {
    20     PS_INT_CHECK_NON_NEGATIVE(n1, NULL);
    21     PS_INT_CHECK_NON_NEGATIVE(n2, NULL);
    22 
    23     psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
    24     pt->x = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
    25     pt->y = psDPolynomial2DAlloc(n1, n2, PS_POLYNOMIAL_ORD);
    26 
    27     p_psMemSetDeallocator(pt, (psFreeFcn) planeTransformFree);
    28     return(pt);
    29 }
    30 
  • trunk/stac/src/stac.c

    r2897 r3375  
    7777
    7878#ifdef TESTING
    79     char prefile[MAXCHAR];              // Filename of precombined image
    80     sprintf(prefile,"%s.pre",config->output);
    81     psImageWriteSection(combined,0,0,0,NULL,0,prefile);
     79    char preName[MAXCHAR];              // Filename of precombined image
     80    sprintf(preName,"%s.pre",config->output);
     81    psFits *preFile = psFitsAlloc(preName);
     82    if (!psFitsWriteImage(preFile, NULL, combined, 0, NULL)) {
     83        psErrorStackPrint(stderr, "Unable to write image: %s\n", preName);
     84    }
     85    psTrace("stac", 1, "Pre-combined image written to %s\n", preName);
     86    psFree(preFile);
    8287#endif
    8388
     
    8994                                              ((psImage*)(inputs->data[i]))->numRows);
    9095#ifdef TESTING
    91         char regionFile[MAXCHAR];       // Filename of region image
    92         sprintf(regionFile,"%s.region",config->inputs->data[i]);
    93         psImageWriteSection(regions->data[i], 0, 0, 0, NULL, 0, regionFile);
     96        char regionName[MAXCHAR];       // Filename of region image
     97        sprintf(regionName,"%s.region",config->inputs->data[i]);
     98        psFits *regionFile = psFitsAlloc(regionName);
     99        if (!psFitsWriteImage(regionFile, NULL, regions->data[i], 0, NULL)) {
     100            psErrorStackPrint(stderr, "Unable to write image: %s\n", regionName);
     101        }
     102        psTrace("stac", 1, "Region image written to %s\n", regionName);
     103        psFree(regionFile);
    94104#endif
    95105    }
     
    122132
    123133    // Write output image
    124     psImageWriteSection(combined,0,0,0,NULL,0,config->output);
     134    psFits *outFile = psFitsAlloc(config->output);
     135    if (!psFitsWriteImage(outFile, NULL, combined, 0, NULL)) {
     136        psErrorStackPrint(stderr, "Unable to write image: %s\n", config->output);
     137    }
     138    psTrace("stac", 1, "Combined image written to %s\n", config->output);
     139    psFree(outFile);
    125140
    126141    // Free everything I've used
  • trunk/stac/src/stacAreaOfInterest.c

    r2783 r3375  
    11#include <stdio.h>
     2#include <assert.h>
    23#include "pslib.h"
    34#include "stac.h"
     
    1516    psDPolynomial2D *yPoly = transform->y;
    1617
    17     if (xPoly->type != PS_POLYNOMIAL_ORD || yPoly->type != PS_POLYNOMIAL_ORD) {
    18         psError("stac.area.deriv", "Transformation polynomials aren't ordinary polynomials.\n");
    19         return false;
    20     }
     18    assert(xPoly->type == PS_POLYNOMIAL_ORD);
     19    assert(yPoly->type == PS_POLYNOMIAL_ORD);
    2120
    2221    // Initialise derivatives
     
    6463{
    6564    // Input checks
    66     if (mask == NULL) {
    67         psError("stac.area", "Input mask is NULL for some weird reason.\n");
    68         return NULL;
    69     }
    70     if (map == NULL) {
    71         psError("stac.area", "Input map is NULL for some weird reason.\n");
    72         return NULL;
    73     }
    74     if (mask->type.type != PS_TYPE_U8) {
    75         psError("stac.area", "Input mask must be of type U8.\n");
    76         return NULL;
    77     }
     65    assert(mask);
     66    assert(map);
     67    assert(mask->type.type == PS_TYPE_U8);
    7868   
    7969    // Number of pixels in x and y
  • trunk/stac/src/stacCheckMemory.c

    r2500 r3375  
    4444
    4545    if ((leakFile = fopen(LEAKS, "w")) == NULL) {
    46         psError("stac.checkMemory", "Unable to open leaks file, %s\n",LEAKS);
     46        fprintf(stderr, "Unable to open leaks file, %s\n", LEAKS);
    4747        return;
    4848    }
    4949
    50     int nLeaks = psMemCheckLeaks(0, &leaks, leakFile); // Number of leaks
     50    int nLeaks = psMemCheckLeaks(0, &leaks, leakFile, false); // Number of leaks
    5151    psTrace("stac.checkMemory", 1, "%d leaks found.\n", nLeaks);
    5252    for (int i = 0; i < nLeaks; i++) {
  • trunk/stac/src/stacCombine.c

    r2783 r3375  
    11#include <stdio.h>
     2#include <assert.h>
    23#include <math.h>
    34#include "pslib.h"
     
    1516    // Would like to use psVectorStats, but it doesn't have errors built in yet
    1617    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    17     (void)psVectorStats(stats, values, masks, 1);
     18    (void)psVectorStats(stats, values, NULL, masks, 1);
    1819    float mean = stats->sampleMean;
    1920    psFree(stats);
     
    4748{
    4849    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    49     (void)psVectorStats(stats, values, masks, 1);
     50    (void)psVectorStats(stats, values, NULL, masks, 1);
    5051    float median = stats->sampleMedian;
    5152    psFree(stats);
     
    7677        psImage *error = (psImage *)errors->data[i]; // The error image
    7778
    78         if ((image->numCols != numCols) || (image->numRows != numRows)) {
    79             psError("stac.combine",
    80                     "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
    81                     i, numCols, numRows, image->numCols, image->numRows);
    82             return false;
    83         }
    84         if ((error->numCols != numCols) || (error->numRows != numRows)) {
    85             psError("stac.combine",
    86                     "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
    87                     i, numCols, numRows, error->numCols, error->numRows);
    88             return false;
    89         }
     79        assert(image->numCols == numCols && image->numRows == numRows);
     80        assert(error->numCols == numCols && error->numRows == numRows);
    9081    }
    9182
    9283    // Check combined image
     84    assert(!*combined || ((*combined)->numRows == numRows) && ((*combined)->numCols == numCols));
    9385    if (*combined == NULL) {
    9486        *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
    95     } else if (((*combined)->numRows != numRows) || ((*combined)->numCols != numCols)) {
    96         psError("stac.combine", "Size of combined image (%dx%d) does not match inputs (%dx%d)\n",
    97                 (*combined)->numCols, (*combined)->numRows, numCols, numRows);
    98         return false;
    9987    }
    10088
    10189    // Check area of interest
    102     if (region && ((region->numRows != numRows) || (region->numCols != numCols))) {
    103         psError("stac.combine", "Size of area of interest (%dx%d) does not match inputs (%dx%d)\n",
    104                 region->numCols, region->numRows, numCols, numRows);
    105         return false;
    106     }
     90    assert(!region || (region->numRows == numRows) && (region->numCols == numCols));
    10791
    10892    psTrace("stac.combine", 1, "Combining images....\n");
     
    117101            // Allocate the rejection masks, if required
    118102            *rejected = psArrayAlloc(nImages);
    119         } else if ((*rejected)->n != nImages) {
    120             psError("stac.combine", "Number of rejection masks (%d) does not match number of input"
    121                     "images (%d).\n", (*rejected)->n, nImages);
    122             return NULL;
     103        } else {
     104            assert((*rejected)->n != nImages);
    123105        }
    124106
     
    208190    if (nReject > 0) {
    209191        for (int i = 0; i < nImages; i++) {
    210             char rejfile[MAXCHAR];      // Filename of rejection image
    211             sprintf(rejfile,"%s.shiftrej",config->inputs->data[i]);
    212             psImageWriteSection((psImage*)((*rejected)->data[i]),0,0,0,NULL,0,rejfile);
     192            char rejName[MAXCHAR];      // Filename of rejection image
     193            sprintf(rejName,"%s.shiftrej",config->inputs->data[i]);
     194
     195            psFits *rejFile = psFitsAlloc(rejName);
     196            if (!psFitsWriteImage(rejFile, NULL, (*rejected)->data[i], 0, NULL)) {
     197                psErrorStackPrint(stderr, "Unable to write image: %s\n", rejName);
     198            }
     199            psTrace("stac", 1, "Rejection image written to %s\n", rejName);
     200            psFree(rejFile);
    213201        }
    214202    }
    215203    // Write chi^2 image
    216204    iteration++;
    217     char chifile[MAXCHAR];              // Filename of chi^2 image
    218     sprintf(chifile,"chi2_%d.fits",iteration);
    219     psImageWriteSection(chi2,0,0,0,NULL,0,chifile);
     205    char chiName[MAXCHAR];              // Filename of chi^2 image
     206    sprintf(chiName,"chi2_%d.fits",iteration);
     207    psFits *chiFile = psFitsAlloc(chiName);
     208    if (!psFitsWriteImage(chiFile, NULL, chi2 , 0, NULL)) {
     209        psErrorStackPrint(stderr, "Unable to write image: %s\n", chiName);
     210    }
     211    psTrace("stac", 1, "Chi^2 image written to %s\n", chiName);
     212    psFree(chiFile);
    220213    psFree(chi2);
    221214#endif
  • trunk/stac/src/stacErrorImages.c

    r2500 r3375  
    3535#ifdef TESTING
    3636        // Write error image out to check
    37         char errfile[MAXCHAR];          // Filename of error image
    38         sprintf(errfile,"%s.err",config->inputs->data[i]);
    39         psImageWriteSection(error,0,0,0,NULL,0,errfile);
     37        char errName[MAXCHAR];          // Filename of error image
     38        sprintf(errName,"%s.err",config->inputs->data[i]);
     39        psFits *errorFile = psFitsAlloc(errName);
     40        if (!psFitsWriteImage(errorFile, NULL, error, 0, NULL)) {
     41            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
     42        }
     43        psTrace("stac", 1, "Error image written to %s\n", errName);
     44        psFree(errorFile);
    4045#endif
    4146
  • trunk/stac/src/stacInvertMaps.c

    r2897 r3375  
    11#include <stdio.h>
     2#include <assert.h>
    23#include "pslib.h"
    34#include "stac.h"
     
    2425        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[mapNum]; // Uninverted map
    2526        // 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         }
     27        assert(oldMap->x->nX == oldMap->x->nY && oldMap->y->nX == oldMap->y->nY &&
     28               oldMap->x->nX == oldMap->y->nX);
    3829        int order = oldMap->x->nX;      // Polynomial order
    3930        psTrace("stac.invertMaps", 4, "Generating order %d polynomial inverse transformation.\n", order);
     
    9182
    9283                    fakePoly->mask[i][j] = 0;
    93                     double ijPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
     84                    double ijPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
    9485                    fakePoly->mask[i][j] = 1;
    9586
     
    9889                           
    9990                            fakePoly->mask[m][n] = 0;
    100                             double mnPoly = psDPolynomial2DEval(xIn->data.F32[g], yIn->data.F32[g], fakePoly);
     91                            double mnPoly = psDPolynomial2DEval(fakePoly, xIn->data.F32[g], yIn->data.F32[g]);
    10192                            fakePoly->mask[m][n] = 1;
    10293
     
    113104        // Solution via LU Decomposition
    114105        psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
    115         psImage *luMatrix = psMatrixLUD(NULL, permutation, matrix); // LU decomposed matrix
     106        psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
    116107        psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
    117108        psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
     
    178169#if 0
    179170        // 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 
     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);
    190175        psPlaneTransform *newMap = psPlaneTransformAlloc(2, 2); // Inverted map
    191176        psPlaneTransform *oldMap = (psPlaneTransform*)maps->data[i]; // Uninverted map
  • trunk/stac/src/stacRead.c

    r2897 r3375  
    88    psArray *filenames = config->inputs;// The file names
    99    int nFiles = filenames->n;          // The number of input files
    10     psArray *files = psArrayAlloc(nFiles); // The input files, to be returned
     10    psArray *images = psArrayAlloc(nFiles); // The input files, to be returned
     11    psRegion *imageRegion = psRegionAlloc(0, 0, 0, 0); // Region of image to read
    1112
    12     psTrace("stac.read.images",1,"Reading input images....\n");
     13    psTrace("stac.read.images", 1, "Reading input images....\n");
    1314    for (int i = 0; i < nFiles; i++) {
     15        psTrace("stac.read.images", 2, "Reading input image %s....\n",filenames->data[i]);
     16        psFits *imageFile = psFitsAlloc(filenames->data[i]);
    1417        // We only read PHUs --- not mucking around with extensions for now
    15         psTrace("stac.read.images",2,"Reading input image %s....\n",filenames->data[i]);
    16         files->data[i] = psImageReadSection(NULL, 0, 0, 0, 0, 0, NULL, 0, filenames->data[i]);
    17         if (files->data[i] == NULL) {
     18        psImage *image = psFitsReadImage(NULL, imageFile, *imageRegion, 0);
     19        if (image == NULL) {
    1820            psErrorStackPrint(stderr,"Fatal error: Unable to read %s\n",filenames->data[i]);
    1921            exit(EXIT_FAILURE);
    2022        }
    21         // Big assumption: input images are 32-bit.
    22         if (((psImage*)(files->data[i]))->type.type != PS_TYPE_F32) {
    23             psError("stac.read.images", "Image %d is not 32-bit\n",i);
     23        psTrace("stac.read.images", 4, "Image %s is %dx%d\n", filenames->data[i], image->numCols,
     24                image->numRows);
     25        // Convert to 32-bit floating point, in necessary
     26        if (image->type.type != PS_TYPE_F32) {
     27            psTrace("stac.read.images", 3, "Converting %s to floating point in memory....",
     28                    filenames->data[i]);
     29            images->data[i] = psImageCopy(NULL, image, PS_TYPE_F32);
     30            psFree(image);
     31        } else {
     32            images->data[i] = image;
    2433        }
     34        psFree(imageFile);
    2535    }
    2636    psTrace("stac.read.images",1,"%d input images read.\n",nFiles);
     37    psFree(imageRegion);
    2738
    28     return files;
     39    return images;
    2940}
    3041
  • trunk/stac/src/stacRejection.c

    r2765 r3375  
    11#include <stdio.h>
     2#include <assert.h>
    23#include "pslib.h"
    34#include "stac.h"
     
    3435    // Get the median
    3536    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    36     (void)psVectorStats(stats, pixels, mask, 1);
     37    (void)psVectorStats(stats, pixels, NULL, mask, 1);
    3738    float median = stats->sampleMedian;
    3839    psFree(stats);
     
    5354
    5455    // Check inputs
    55     if (inputs->n != rejected->n) {
    56         psError("stac.rejection", "Number of input images (%d) and rejection masks (%d) does not match.\n",
    57                 inputs->n, rejected->n);
    58         return NULL;
    59     }
    60     if (inputs->n != regions->n) {
    61         psError("stac.rejection", "Number of input images (%d) and region masks (%d) does not match.\n",
    62                 inputs->n, regions->n);
    63         return NULL;
    64     }
    65     if (inputs->n != maps->n) {
    66         psError("stac.rejection", "Number of input images (%d) and transformation maps (%d) does "
    67                 "not match.\n", inputs->n, maps->n);
    68         return NULL;
    69     }
    70     if (inputs->n != inverseMaps->n) {
    71         psError("stac.rejection", "Number of input images (%d) and inverse transformation maps (%d) does "
    72                 "not match.\n", inputs->n, inverseMaps->n);
    73         return NULL;
    74     }
     56    assert(inputs->n == rejected->n);
     57    assert(inputs->n == regions->n);
     58    assert(inputs->n == maps->n);
     59    assert(inputs->n == inverseMaps->n);
     60
    7561    for (int i = 0; i < nImages; i++) {
    76         if ((((psImage*)(inputs->data[i]))->numRows != ((psImage*)(regions->data[i]))->numRows) ||
    77             (((psImage*)(inputs->data[i]))->numCols != ((psImage*)(regions->data[i]))->numCols)) {
    78             psError("stac.rejection",
    79                     "Sizes of input image (%dx%d) and region mask (%dx%d) for input %d do not match.\n",
    80                     ((psImage*)(inputs->data[i]))->numCols, ((psImage*)(inputs->data[i]))->numRows,
    81                     ((psImage*)(regions->data[i]))->numCols, ((psImage*)(regions->data[i]))->numRows, i);
    82             return NULL;
    83         }
     62        psImage *input = inputs->data[i];
     63        psImage *region = regions->data[i];
     64        assert(input->numRows == region->numRows && input->numCols == region->numCols);
    8465    }
    8566
     
    11495        sprintf(crfile,"%s.cr",config->inputs->data[i]);
    11596        if ((crs = fopen(crfile, "w")) == NULL) {
    116             psError("stac.rejection","Unable to open file for detailed output, %s\n");
     97            fprintf(stderr, "Unable to open file for detailed output, %s\n");
    11798            return NULL;
    11899        }
     
    198179#ifdef TESTING
    199180        // Write error image out to check
    200         char maskfile[MAXCHAR];         // Filename of mask image
    201         char rejmapfile[MAXCHAR];       // Filename of rejection image
    202         char gradfile[MAXCHAR];         // Filename of gradient image
    203         sprintf(maskfile,"%s.mask",config->inputs->data[i]);
    204         sprintf(rejmapfile,"%s.rejmap",config->inputs->data[i]);
    205         sprintf(gradfile,"%s.grad",config->inputs->data[i]);
    206         psImageWriteSection(mask,0,0,0,NULL,0,maskfile);
    207         psImageWriteSection(rejmap,0,0,0,NULL,0,rejmapfile);
    208         psImageWriteSection(grad,0,0,0,NULL,0,gradfile);
     181        char maskName[MAXCHAR];         // Filename of mask image
     182        char rejmapName[MAXCHAR];       // Filename of rejection image
     183        char gradName[MAXCHAR];         // Filename of gradient image
     184        sprintf(maskName,"%s.mask",config->inputs->data[i]);
     185        sprintf(rejmapName,"%s.rejmap",config->inputs->data[i]);
     186        sprintf(gradName,"%s.grad",config->inputs->data[i]);
     187
     188        psFits *maskFile = psFitsAlloc(maskName);
     189        psFits *rejmapFile = psFitsAlloc(rejmapName);
     190        psFits *gradFile = psFitsAlloc(gradName);
     191        if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
     192            psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
     193        }
     194        psTrace("stac", 1, "Mask image written to %s\n", maskName);
     195        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
     196            psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
     197        }
     198        psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
     199        if (!psFitsWriteImage(gradFile, NULL, grad, 0, NULL)) {
     200            psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
     201        }
     202        psTrace("stac", 1, "Gradient image written to %s\n", gradName);
     203        psFree(maskFile);
     204        psFree(rejmapFile);
     205        psFree(gradFile);
    209206        psFree(rejmap);
    210207        psFree(grad);
  • trunk/stac/src/stacSize.c

    r2897 r3375  
    11#include <stdio.h>
     2#include <assert.h>
    23#include "pslib.h"
    34#include "stac.h"
     
    910{
    1011    int nImages = images->n;            // Number of input images
    11     if (nImages != maps->n) {
    12         psError("stac.size", "Number of images (%d) and maps (%d) do not match.\n", nImages, maps->n);
    13         return false;
    14     }
     12    assert(maps->n == nImages);
    1513
    1614    psPlane *inCoord = psAlloc(sizeof(psPlane)); // Input coordinates
  • trunk/stac/src/stacTransform.c

    r2783 r3375  
    11#include <stdio.h>
     2#include <assert.h>
    23#include "pslib.h"
    34#include "stac.h"
     
    120121
    121122    // Check input sizes
    122     if (images->n != maps->n) {
    123         psError("stac.transform", "Number of maps (%d) does not match number of images (%d).\n",
    124                 maps->n, images->n);
    125         return false;
    126     }
    127     if (errors && (images->n != errors->n)) {
    128         psError("stac.transform", "Number of error images (%d) does not match number of images (%d).\n",
    129                 errors->n, images->n);
    130         return false;
    131     }
     123    assert(images->n == maps->n);
     124    assert(!errors || (images->n == errors->n));
    132125
    133126    // Allocate the output images if required, otherwise check the number
     127    assert(!*outputs || (*outputs)->n == nImages);
    134128    if (*outputs == NULL) {
    135129        *outputs = psArrayAlloc(nImages);
     
    138132            (*outputs)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
    139133        }
    140     } else if ((*outputs)->n != nImages) {
    141         psError("stac.transform", "Number of output images (%d) does not match number of input images "
    142                 "(%d).\n", (*outputs)->n, nImages);
    143         return false;
    144134    }
    145135
    146136    // Allocate the output error images, if required, otherwise check the number
     137    assert(!errors || ! *outErrors || errors->n == (*outErrors)->n);
    147138    if (errors && (*outErrors == NULL)) {
    148139        *outErrors = psArrayAlloc(errors->n);
     
    151142            (*outErrors)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
    152143        }
    153     } else if (errors->n != (*outErrors)->n) {
    154         psError("stac.transform", "Number of error output images (%d) does not match number of input"
    155                 "images (%d).\n", (*outErrors)->n, errors->n);
    156         return false;
    157144    }
    158145
    159146    // Check the masks, if specified
     147    assert(!masks || masks->n == nImages);
    160148    if (masks != NULL) {
    161         if (masks->n != nImages) {
    162             psError("stac.transform", "Number of masks (%d) does not match number of input images (%d).\n",
    163                     masks->n, nImages);
    164             return false;
    165         }
    166149        for (int i = 0; i < nImages; i++) {
    167150            psImage *image = images->data[i];
    168151            psImage *mask = masks->data[i];
    169             if ((mask->numRows != image->numRows) || (mask->numCols != image->numCols)) {
    170                 psError ("stac.transform",
    171                          "Size of input mask (%dx%d) does not match that of input image (%dx%d).\n",
    172                          mask->numCols, mask->numRows, image->numCols, image->numRows);
    173                 return false;
    174             }
     152            assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
    175153        }
    176154    }
     
    236214#ifdef TESTING
    237215        // Write error image out to check
    238         char shiftfile[MAXCHAR];        // Filename of shift image
    239         char errfile[MAXCHAR];          // Filename of error image
    240         sprintf(shiftfile,"%s.shift.%d",config->inputs->data[n],numTransforms);
    241         sprintf(errfile,"%s.shifterr.%d",config->inputs->data[n],numTransforms);
     216        char shiftName[MAXCHAR];        // Filename of shift image
     217        char errName[MAXCHAR];          // Filename of error image
     218        sprintf(shiftName,"%s.shift.%d",config->inputs->data[n],numTransforms);
     219        sprintf(errName,"%s.shifterr.%d",config->inputs->data[n],numTransforms);
    242220        psTrace("stac.transform.test", 6,
    243221                "Output files have size: %dx%d\n",outImage->numCols,outImage->numRows);
    244         psImageWriteSection(outImage,0,0,0,NULL,0,shiftfile);
    245         psImageWriteSection(outError,0,0,0,NULL,0,errfile);
     222
     223        psFits *shiftFile = psFitsAlloc(shiftName);
     224        psFits *errFile = psFitsAlloc(errName);
     225        if (!psFitsWriteImage(shiftFile, NULL, outImage, 0, NULL)) {
     226            psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
     227        }
     228        psTrace("stac", 1, "Shifted image written to %s\n", shiftName);
     229        if (!psFitsWriteImage(errFile, NULL, outError, 0, NULL)) {
     230            psErrorStackPrint(stderr, "Unable to write image: %s\n", errName);
     231        }
     232        psTrace("stac", 1, "Shifted error image written to %s\n", errName);
     233        psFree(shiftFile);
     234        psFree(errFile);
    246235#endif
    247236
Note: See TracChangeset for help on using the changeset viewer.