IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3610


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

Working very nicely now.

Location:
trunk/stac/src
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/stac/src/Makefile

    r3387 r3610  
    11SHELL = /bin/sh
    22CC = gcc
    3 CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -DCRFLUX -DTESTING
    4 PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2
     3CFLAGS += -g -std=c99 -I/home/mithrandir/price/psLib/include -D_GNU_SOURCE # -DCRFLUX -DTESTING
     4PSLIB += -L/home/mithrandir/price/psLib/lib/ -lpslib -lgsl -lgslcblas -lfftw3f -lsla -lcfitsio -lm -lxml2 -lmysqlclient
    55LDFLAGS += $(PSLIB)
    66
    77STAC = stac.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
    8         stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o
     8        stacCombine.o stacInvertMaps.o stacRejection.o stacAreaOfInterest.o stacSize.o stacScales.o
    99
    1010SHIFT = shift.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
     
    1515TEST_SHIFT = testShift.o stacConfig.o stacRead.o stacErrorImages.o stacTransform.o stacCheckMemory.o \
    1616        stacInvertMaps.o
     17
     18PHOTSCALE = photScale.o
    1719
    1820.PHONY: tags clean empty test profile optimise
     
    2628shift:          $(SHIFT)
    2729                $(CC) $(CFLAGS) -o $@ $(SHIFT) $(LDFLAGS) $(OPTFLAGS)
     30
     31photScale:      $(PHOTSCALE)
     32                $(CC) $(CFLAGS) -o $@ $(PHOTSCALE) $(LDFLAGS) $(OPTFLAGS)
    2833
    2934calcGrad:       calcGrad.o
  • trunk/stac/src/stac.c

    r3435 r3610  
    2727    // Set trace levels
    2828    psTraceSetLevel(".",0);
     29    psTraceSetLevel(".psLib", 0);
    2930    psTraceSetLevel("stac.checkMemory",10);
    3031    psTraceSetLevel("stac.config",10);
     
    3839    psTraceSetLevel("stac.area",10);
    3940    psTraceSetLevel("stac.size",10);
     41    psTraceSetLevel("stac.scales",7);
    4042
    4143    // Set logging level
     
    6567    psArray *transformed = NULL;
    6668    psArray *transformedErrors = NULL;
    67     (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, config);
     69    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, NULL, NULL,
     70                        config);
    6871
    6972    psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
     73
     74    psVector *scales = NULL;            // Relative scales between images
     75    psVector *offsets = NULL;           // Offsets between images
     76    (void)stacScales(&scales, &offsets, transformed, config);
     77    // Rescale the images
     78    (void)stacRescale(transformed, transformedErrors, NULL, scales, offsets);
    7079
    7180    // Combine with rejection
     
    122131    }
    123132
    124     // Redo transformation with the masks
     133    // Redo transformation with the masks and scales/offsets
    125134    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
    126                         combineRegion, config);
     135                        combineRegion, scales, offsets, config);
    127136
    128137    // Combine the newly-transformed CR-free images, no rejection
     
    139148    psFree(outFile);
    140149
     150    // Save shifted images
     151    if (config->saveShifts) {
     152        for (int i = 0; i < transformed->n; i++) {
     153            char shiftName[MAXCHAR];    // Filename of shifted image
     154            sprintf(shiftName, "%s.shift", config->inputs->data[i]);
     155            psFits *shiftFile = psFitsAlloc(shiftName);
     156            if (!psFitsWriteImage(shiftFile, NULL, transformed->data[i], 0, NULL)) {
     157                psErrorStackPrint(stderr, "Unable to write image: %s\n", shiftName);
     158            }
     159            psTrace("stac", 1, "Shifted image %d written to %s\n", i, shiftName);
     160            psFree(shiftFile);
     161        }
     162    }
     163
    141164    // Free everything I've used
    142165    stacConfigFree(config);
     166    psFree(scales);
     167    psFree(offsets);
    143168    psFree(combineRegion);
    144169    psFree(regions);
  • trunk/stac/src/stac.h

    r2897 r3610  
    88#define MAXCHAR 80
    99
    10 // Temporary implementation
    11 psPlaneTransform *psPlaneTransformAlloc(psS32 n1, psS32 n2);
     10// Get the current time
     11double getTime(void);
    1212
    1313/************************************************************************************************************/
     
    2121    psArray *inputs;                    // Input file names
    2222    const char *output;                 // Output file name
     23    bool saveShifts;                    // Save shifted images?
     24    const char *starFile;               // File with star coordinates
     25    const char *starMapFile;            // File with map for stars
    2326    int outnx, outny;                   // Size of output image
    2427    float xOffset, yOffset;             // Offset to get lower-left corner at 0,0
    25     float saturated;                    // Saturation level
    26     float bad;                          // Bad level
     28    psVector *saturated;                // Saturation level for each image
     29    psVector *bad;                      // Bad level for each image
    2730    float reject;                       // Rejection level
    2831    float frac;                         // Fraction of input pixel that must be masked before the pixel is
     
    3033    float grad;                         // Multiplier of the gradient
    3134    int nReject;                        // Number of rejection iterations
     35    float xMapDiff, yMapDiff;           // Difference between pure map and output image coordinates
     36    float aper;                         // Aperture size (pixels)
    3237} stacConfig;
    3338
     
    5257// Read the input files and return an array of images
    5358psArray *stacReadImages(stacConfig *config);
     59
     60// Read a file of coordinates (x,y)
     61psArray *stacReadCoords(const char *filename);
     62
     63// Read a single map file and return the transformation
     64psPlaneTransform *stacReadMap(const char *filename);
    5465
    5566// Read the map files and return an array of transformations
     
    7788                   const psArray *masks, // Masks of input images
    7889                   const psImage *region, // Region of interest for transformation
     90                   const psVector *scales, // Relative scales
     91                   const psVector *offsets, // Relative offsets
    7992                   const stacConfig *config // Configuration
    8093    );
     
    181194/************************************************************************************************************/
    182195
     196// stacScale.c
     197
     198// Calculate the background in an image
     199float stacBackground(const psImage *image, // Image for which to get the background
     200                     int sample         // Sample in increments of this value
     201    );
     202
     203
     204// Calculate the relative scales and offsets between the images
     205bool stacScales(psVector **scalesPtr,   // Scales to return
     206                psVector **offsetsPtr,  // Offsets to return
     207                const psArray *images,  // Images on which to measure the scales and offsets
     208                const stacConfig *config // Configuration
     209    );
     210
     211// Rescale images
     212bool stacRescale(psArray *images,       // Images to rescale
     213                 psArray *errImages,    // Variance images to rescale
     214                 const psImage *mask,   // Mask indicating which pixels to scale
     215                 const psVector *scales,// Scales for images
     216                 const psVector *offsets // Offsets for images
     217    );
     218
    183219#endif
  • trunk/stac/src/stacCombine.c

    r3375 r3610  
    6868    int numRows = ((psImage*)images->data[0])->numRows; // Image size
    6969    int numCols = ((psImage*)images->data[0])->numCols; // Image size
    70     float saturated = config->saturated;// Saturation limit
    71     float bad = config->bad;            // Bad pixel limit
     70    psVector *saturated = config->saturated;// Saturation limits
     71    psVector *bad = config->bad;        // Bad pixel limits
    7272    float reject = config->reject;      // Rejection (k-sigma)
    7373   
     
    134134                    pixels->data.F32[i] = pixel;
    135135                    deltas->data.F32[i] = delta;
    136                     if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
     136                    if ((pixel >= saturated->data.F32[i]) || (pixel <= bad->data.F32[i]) ||
     137                        (! isfinite(pixel)) || (! isfinite(delta))) {
    137138                        mask->data.U8[i] = (psU8)1; // Don't use!
    138139                    } else {
  • trunk/stac/src/stacConfig.c

    r3387 r3610  
    99{
    1010    fprintf (stderr, "STAC: Simultaneous Telescope Array Combination\n"
    11              "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-n NREJECT] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
     11             "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-S] [-s SAT] [-b BAD] [-p FILE MAP] [-a APER] [-k REJ] [-n NREJECT] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
    1212             "where\n"
    1313             "\t-h           Help (this info)\n"
     
    1616             "\t-r           Read noise (e) for detectors\n"
    1717             "\t-o NX NY     Size of output image (512, 512)\n"
    18              "\t-s SAT       Saturation level (65536)\n"
     18             "\t-S           Save shifted images (false)\n"
     19             "\t-s SAT       Saturation level (65535)\n"
    1920             "\t-b BAD       Bad level (0)\n"
     21             "\t-p FILE MAP  Specify file containing star coordinates, with map\n"
     22             "\t-a APER      Aperture size to use for photometry (3 pix)\n"
    2023             "\t-k REJ       Rejection level (k-sigma; 3.5)\n"
    2124             "\t-n NREJECT   Number of rejection iterations (2)\n"
     
    3942    config->inputs = NULL;
    4043    config->output = NULL;
     44    config->saveShifts = false;
     45    config->starFile = NULL;
     46    config->starMapFile = NULL;
     47    config->aper = 3.0;
    4148    config->outnx = 512;
    4249    config->outny = 512;
    43     config->saturated = 65535.0;
    44     config->bad = 0.0;
     50    config->saturated = NULL;           // Saturations; will fill this in later, when we know how many images
     51    config->bad = NULL;                 // Bad levels; will fill this in later, when we know how many images
    4552    config->reject = 3.0;
    4653    config->frac = 0.5;
     
    5461void stacConfigFree(stacConfig *config)
    5562{
    56     // Free the vector, if necessary
     63    // Free the vectors, if necessary
    5764    if (config->inputs) {
    5865        psFree(config->inputs);
     66    }
     67    if (config->saturated) {
     68        psFree(config->saturated);
     69    }
     70    if (config->bad) {
     71        psFree(config->bad);
    5972    }
    6073    // Free everything
     
    6982    stacConfig *config = stacConfigAlloc(); // Configuration values
    7083    const char *programName = argv[0];  // Program name
     84
     85    float saturated = 65535.0;          // Saturation level
     86    float bad = 0.0;                    // Bad level
    7187
    7288    /* Variables for getopt */
     
    7692
    7793    /* Parse command-line arguments using getopt */
    78     while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:n:f:G:")) != -1) {
     94    while ((opt = getopt(argc, argv, "hvSg:r:o:s:b:k:n:f:G:p:a:")) != -1) {
    7995        switch (opt) {
    8096          case 'h':
     
    87103            if (sscanf(optarg, "%f", &config->gain) != 1) {
    88104                help(programName);
     105                exit(EXIT_FAILURE);
    89106            }
    90107            break;
     
    92109            if (sscanf(optarg, "%f", &config->readnoise) != 1) {
    93110                help(programName);
     111                exit(EXIT_FAILURE);
    94112            }
    95113            break;
     
    97115            if ((sscanf(argv[optind-1], "%d", &config->outnx) != 1) ||
    98116                (sscanf(argv[optind++], "%d", &config->outny) != 1)) {
    99                 /*
    100                   Note: incrementing optind, so I can read more than
    101                   one parameter.
    102                 */
     117                // Note: incrementing optind, so I can read more than one parameter.
    103118                help(programName);
    104119                exit(EXIT_FAILURE);
     
    106121            break;
    107122          case 's':
    108             if (sscanf(optarg, "%f", &config->saturated) != 1) {
    109                 help(programName);
     123            if (sscanf(optarg, "%f", &saturated) != 1) {
     124                help(programName);
     125                exit(EXIT_FAILURE);
    110126            }
    111127            break;
    112128          case 'b':
    113             if (sscanf(optarg, "%f", &config->bad) != 1) {
    114                 help(programName);
     129            if (sscanf(optarg, "%f", &bad) != 1) {
     130                help(programName);
     131                exit(EXIT_FAILURE);
     132            }
     133            break;
     134          case 'p':
     135            if (argc < optind+1) {
     136                help(programName);
     137                exit(EXIT_FAILURE);
     138            }
     139            config->starFile = argv[optind-1];
     140            config->starMapFile = argv[optind++];
     141            // Note: incrementing optind, so I can read more than one parameter.
     142            break;
     143          case 'a':
     144            if (sscanf(optarg, "%f", &config->aper) != 1) {
     145                help(programName);
     146                exit(EXIT_FAILURE);
    115147            }
    116148            break;
     
    118150            if (sscanf(optarg, "%f", &config->reject) != 1) {
    119151                help(programName);
     152                exit(EXIT_FAILURE);
    120153            }
    121154            break;
     
    123156            if (sscanf(optarg, "%d", &config->nReject) != 1) {
    124157                help(programName);
     158                exit(EXIT_FAILURE);
    125159            }
    126160            break;
     
    128162            if (sscanf(optarg, "%f", &config->frac) != 1) {
    129163                help(programName);
     164                exit(EXIT_FAILURE);
    130165            }
    131166            break;
     
    133168            if (sscanf(optarg, "%f", &config->grad) != 1) {
    134169                help(programName);
    135             }
     170                exit(EXIT_FAILURE);
     171            }
     172            break;
     173          case 'S':
     174            config->saveShifts = true;
    136175            break;
    137176          default:
     
    156195    }
    157196
     197    // Create the saturation and bad vectors
     198    config->saturated = psVectorAlloc(argc-1, PS_TYPE_F32);
     199    config->bad = psVectorAlloc(argc-1, PS_TYPE_F32);
     200    for (int i = 0; i < argc - 1; i++) {
     201        ((psVector*)config->saturated)->data.F32[i] = saturated;
     202        ((psVector*)config->bad)->data.F32[i] = bad;
     203    }
     204
    158205    // Debugging output
    159206    psTrace("stac.config", 8, "Parsed command line for configuration\n");
  • 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    }
  • trunk/stac/src/stacSize.c

    r3469 r3610  
    9191
    9292    // Tweak the maps to account for the offset
     93    config->xMapDiff = floor(xMin);
     94    config->yMapDiff = floor(yMin);
    9395    for (int i = 0; i < nImages; i++) {
    9496        psPlaneTransform *map = maps->data[i]; // The map of interest
  • trunk/stac/src/stacTransform.c

    r3387 r3610  
    113113                   const psArray *masks, // Masks of input images
    114114                   const psImage *region, // Region of interest for transformation
     115                   const psVector *scales, // Relative scales
     116                   const psVector *offsets, // Relative offsets
    115117                   const stacConfig *config // Configuration
    116118    )
     
    123125    assert(images->n == maps->n);
    124126    assert(!errors || (images->n == errors->n));
     127    assert(!scales || scales->n == images->n);
     128    assert(!offsets || offsets->n == images->n);
     129    assert(!scales || scales->type.type == PS_TYPE_F32);
     130    assert(!offsets || offsets->type.type == PS_TYPE_F32);
    125131
    126132    // Allocate the output images if required, otherwise check the number
     
    170176        psImage *outImage = (*outputs)->data[n]; // The output image
    171177        psImage *outError = (*outErrors)->data[n]; // The output error image
     178        float offset = 0.0;             // Relative offset
     179        float scale = 1.0;              // Relative scale
     180        if (offsets) {
     181            offset = offsets->data.F32[n];
     182        }
     183        if (scales) {
     184            scale = scales->data.F32[n];
     185        }
    172186
    173187#if 0
     
    206220                                                                                            detector->y,
    207221                                                                                            mask, 1, 0.0);
     222
     223                    outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
     224                    outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
     225
    208226                } // Pixels of interest
    209227
Note: See TracChangeset for help on using the changeset viewer.