IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 7, 2005, 4:04:22 PM (21 years ago)
Author:
Paul Price
Message:

Using autoconf; some bug fixes to the polynomial order

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/stac/src/stacTransform.c

    r5743 r5745  
    1010// Hacked the original ps_ImagePixelInterpolateBILINEAR_F32 to add variances
    1111// i.e., to square the fractions when combining.
    12 inline psF64 p_psImageErrorInterpolateBILINEAR_F32(const psImage* input, 
    13                                                    float x,
    14                                                    float y,
    15                                                    const psImage* mask,
    16                                                    unsigned int maskVal,
    17                                                    psF64 unexposedValue)
     12inline psF64 p_psImageErrorInterpolateBILINEAR_F32(const psImage* input,
     13                                                   float x,
     14                                                   float y,
     15                                                   const psImage* mask,
     16                                                   unsigned int maskVal,
     17                                                   psF64 unexposedValue)
    1818{
    19     double floorX = floor((psF64)(x) - 0.5); 
    20     double floorY = floor((psF64)(y) - 0.5); 
    21     psF64 fracX = x - 0.5 - floorX; 
    22     psF64 fracY = y - 0.5 - floorY; 
    23     int intFloorX = (int) floorX; 
    24     int intFloorY = (int) floorY; 
    25     int lastX = input->numCols - 1; 
    26     int lastY = input->numRows - 1; 
    27     psF32 V00;
    28     psF32 V01;
    29     psF32 V10;
    30     psF32 V11;
    31     bool valid00 = false; 
    32     bool valid01 = false; 
    33     bool valid10 = false; 
    34     bool valid11 = false; 
    35    
    36     if (intFloorY >= 0 && intFloorY <= lastY) { 
    37         if (intFloorX >= 0 && intFloorX <= lastX) { 
    38             V00 = input->data.F32[intFloorY][intFloorX]; 
    39             valid00 = true;
    40         } 
    41         if (intFloorX >= -1 && intFloorX < lastX) { 
    42             V10 = input->data.F32[intFloorY][intFloorX+1]; 
    43             valid10 = true;
    44         } 
    45     } 
    46     if (intFloorY >= -1 && intFloorY < lastY) { 
    47         if (intFloorX >= 0 && intFloorX <= lastX) { 
    48             V01 = input->data.F32[intFloorY+1][intFloorX]; 
    49             valid01 = true;
    50         } 
    51         if (intFloorX >= -1 && intFloorX < lastX) { 
    52             V11 = input->data.F32[intFloorY+1][intFloorX+1]; 
    53             valid11 = true;
    54         } 
    55     } 
    56    
    57     /* cover likely case of all pixels being valid more efficiently */ 
    58     if (valid00 && valid10 && valid01 && valid11) { 
    59         /* formula from the ADD */ 
     19    double floorX = floor((psF64)(x) - 0.5);
     20    double floorY = floor((psF64)(y) - 0.5);
     21    psF64 fracX = x - 0.5 - floorX;
     22    psF64 fracY = y - 0.5 - floorY;
     23    int intFloorX = (int) floorX;
     24    int intFloorY = (int) floorY;
     25    int lastX = input->numCols - 1;
     26    int lastY = input->numRows - 1;
     27    psF32 V00 = 0.0;
     28    psF32 V01 = 0.0;
     29    psF32 V10 = 0.0;
     30    psF32 V11 = 0.0;
     31    bool valid00 = false;
     32    bool valid01 = false;
     33    bool valid10 = false;
     34    bool valid11 = false;
     35
     36    if (intFloorY >= 0 && intFloorY <= lastY) {
     37        if (intFloorX >= 0 && intFloorX <= lastX) {
     38            V00 = input->data.F32[intFloorY][intFloorX];
     39            valid00 = true;
     40        }
     41        if (intFloorX >= -1 && intFloorX < lastX) {
     42            V10 = input->data.F32[intFloorY][intFloorX+1];
     43            valid10 = true;
     44        }
     45    }
     46    if (intFloorY >= -1 && intFloorY < lastY) {
     47        if (intFloorX >= 0 && intFloorX <= lastX) {
     48            V01 = input->data.F32[intFloorY+1][intFloorX];
     49            valid01 = true;
     50        }
     51        if (intFloorX >= -1 && intFloorX < lastX) {
     52            V11 = input->data.F32[intFloorY+1][intFloorX+1];
     53            valid11 = true;
     54        }
     55    }
     56
     57    /* cover likely case of all pixels being valid more efficiently */
     58    if (valid00 && valid10 && valid01 && valid11) {
     59        /* formula from the ADD */
    6060        return V00*SQUARE((1.0-fracX)*(1.0-fracY)) + V10*SQUARE(fracX*(1.0-fracY)) +
    61             V01*SQUARE(fracY*(1.0-fracX)) + V11*SQUARE(fracX*fracY);
    62     } 
    63    
    64     /* OK, at least one pixel is not valid - need to do it piecemeal */ 
    65    
    66     psF64 V0;
    67     bool valid0 = true; 
    68     if (valid00 && valid10) { 
    69         V0 = V00*SQUARE(1-fracX)+V10*SQUARE(fracX); 
    70     } else if (valid00) { 
    71         V0 = V00; 
    72     } else if (valid10) { 
    73         V0 = V10; 
    74     } else { 
    75         valid0 = false; 
    76     } 
    77    
    78     psF64 V1;
    79     bool valid1 = true; 
    80     if (valid01 && valid11) { 
    81         V1 = V01*SQUARE(1-fracX)+V11*SQUARE(fracX); 
    82     } else if (valid01) { 
    83         V1 = V01; 
    84     } else if (valid11) { 
    85         V1 = V11; 
    86     } else { 
    87         valid1 = false; 
    88     } 
    89    
    90     if (valid0 && valid1) { 
    91         return ( V0*SQUARE(1-fracY) + V1*SQUARE(fracY) ); 
    92     } else if (valid0) { 
    93         return V0; 
    94     } else if (valid1) { 
    95         return V1; 
    96     } 
    97    
    98     return unexposedValue; 
     61            V01*SQUARE(fracY*(1.0-fracX)) + V11*SQUARE(fracX*fracY);
     62    }
     63
     64    /* OK, at least one pixel is not valid - need to do it piecemeal */
     65
     66    psF64 V0 = 0.0;
     67    bool valid0 = true;
     68    if (valid00 && valid10) {
     69        V0 = V00*SQUARE(1-fracX)+V10*SQUARE(fracX);
     70    } else if (valid00) {
     71        V0 = V00;
     72    } else if (valid10) {
     73        V0 = V10;
     74    } else {
     75        valid0 = false;
     76    }
     77
     78    psF64 V1 = 0.0;
     79    bool valid1 = true;
     80    if (valid01 && valid11) {
     81        V1 = V01*SQUARE(1-fracX)+V11*SQUARE(fracX);
     82    } else if (valid01) {
     83        V1 = V01;
     84    } else if (valid11) {
     85        V1 = V11;
     86    } else {
     87        valid1 = false;
     88    }
     89
     90    if (valid0 && valid1) {
     91        return ( V0*SQUARE(1-fracY) + V1*SQUARE(fracY) );
     92    } else if (valid0) {
     93        return V0;
     94    } else if (valid1) {
     95        return V1;
     96    }
     97
     98    return unexposedValue;
    9999}
    100100
    101101
    102102
    103 bool stacTransform(psArray **outputs,   // Transformed images for output
    104                    psArray **outErrors, // Transformed error images for output
    105                    const psArray *images, // Array of images to be transformed
    106                    const psArray *maps, // Array of polynomials that do the transformation
    107                    const psArray *errors, // Array of error images to be transformed
    108                    const psArray *masks, // Masks of input images
    109                    const psImage *region, // Region of interest for transformation
    110                    const psVector *scales, // Relative scales
    111                    const psVector *offsets, // Relative offsets
    112                    int outnx, int outny // Size of output images
     103bool stacTransform(psArray **outputs,   // Transformed images for output
     104                   psArray **outErrors, // Transformed error images for output
     105                   const psArray *images, // Array of images to be transformed
     106                   const psArray *maps, // Array of polynomials that do the transformation
     107                   const psArray *errors, // Array of error images to be transformed
     108                   const psArray *masks, // Masks of input images
     109                   const psImage *region, // Region of interest for transformation
     110                   const psVector *scales, // Relative scales
     111                   const psVector *offsets, // Relative offsets
     112                   int outnx, int outny // Size of output images
    113113    )
    114114{
    115     int nImages = images->n;            // Number of images
     115    int nImages = images->n;            // Number of images
    116116
    117117    // Check input sizes
    … …  
    126126    assert(!*outputs || (*outputs)->n == nImages);
    127127    if (*outputs == NULL) {
    128         *outputs = psArrayAlloc(nImages);
    129         psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
    130         for (int i = 0; i < nImages; i++) {
    131             (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
    132             psImageInit((*outputs)->data[i], 0.0);
    133         }
     128        *outputs = psArrayAlloc(nImages);
     129        psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", outnx, outny);
     130        for (int i = 0; i < nImages; i++) {
     131            (*outputs)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
     132            psImageInit((*outputs)->data[i], 0.0);
     133        }
    134134    }
    135135
    … …  
    137137    assert(!errors || ! *outErrors || errors->n == (*outErrors)->n);
    138138    if (errors && (*outErrors == NULL)) {
    139         *outErrors = psArrayAlloc(errors->n);
    140         psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
    141         for (int i = 0; i < nImages; i++) {
    142             (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
    143         }
     139        *outErrors = psArrayAlloc(errors->n);
     140        psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", outnx, outny);
     141        for (int i = 0; i < nImages; i++) {
     142            (*outErrors)->data[i] = psImageAlloc(outnx, outny, PS_TYPE_F32);
     143        }
    144144    }
    145145
    … …  
    147147    assert(!masks || masks->n == nImages);
    148148    if (masks != NULL) {
    149         for (int i = 0; i < nImages; i++) {
    150             psImage *image = images->data[i];
    151             psImage *mask = masks->data[i];
    152             assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
    153         }
     149        for (int i = 0; i < nImages; i++) {
     150            psImage *image = images->data[i];
     151            psImage *mask = masks->data[i];
     152            assert(mask->numRows == image->numRows && mask->numCols == image->numCols);
     153        }
    154154    }
    155155
    … …  
    162162    // Iterate over the images
    163163    for (int n = 0; n < nImages; n++) {
    164         psTrace("stac.transform", 1, "Transforming image %d....\n",n);
    165 
    166         // Pull out the various stuff we're working on
    167         psImage *image = images->data[n]; // The input image
    168         psPlaneTransform *map = maps->data[n]; // The map
    169         psImage *outImage = (*outputs)->data[n]; // The output image
    170         psImage *error = NULL; // The error image
    171         psImage *outError = NULL; // The output error image
    172         if (errors) {
    173             error = errors->data[n];
    174             outError = (*outErrors)->data[n];
    175         }
    176         float offset = 0.0;             // Relative offset
    177         float scale = 1.0;              // Relative scale
    178         if (offsets) {
    179             offset = offsets->data.F32[n];
    180         }
    181         if (scales) {
    182             scale = scales->data.F32[n];
    183         }
    184 
    185         // Mask
    186         psImage *mask = NULL;
    187         if (masks != NULL) {
    188             mask = masks->data[n];
    189         }
    190 
    191         // Iterate over the output image pixels
    192         for (int y = 0; y < outny; y++) {
    193             for (int x = 0; x < outnx; x++) {
    194                 // Only transform those pixels requested
    195                 if (!region || (region && region->data.U8[y][x])) {
    196                     // Transform!
    197                     sky->x = (double)x + 0.5;
    198                     sky->y = (double)y + 0.5;
    199                     (void)psPlaneTransformApply(detector, map, sky);
    200                    
    201                     // Change PS_INTERPOLATE_BILINEAR to best available technique.
    202                     outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
    203                                                                               detector->y, mask, 1, NAN,
    204                                                                               PS_INTERPOLATE_BILINEAR);
    205                     if (error) {
    206                         // Error is actually the variance
    207                         outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
    208                                                                                                 detector->x,
    209                                                                                                 detector->y,
    210                                                                                                 mask, 1, NAN);
    211                     }
    212 
    213                     outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
    214                     outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
    215 
    216                 } // Pixels of interest
    217 
    218             }
    219         } // Iterating over output pixels
     164        psTrace("stac.transform", 1, "Transforming image %d....\n",n);
     165
     166        // Pull out the various stuff we're working on
     167        psImage *image = images->data[n]; // The input image
     168        psPlaneTransform *map = maps->data[n]; // The map
     169        psImage *outImage = (*outputs)->data[n]; // The output image
     170        psImage *error = NULL; // The error image
     171        psImage *outError = NULL; // The output error image
     172        if (errors) {
     173            error = errors->data[n];
     174            outError = (*outErrors)->data[n];
     175        }
     176        float offset = 0.0;             // Relative offset
     177        float scale = 1.0;              // Relative scale
     178        if (offsets) {
     179            offset = offsets->data.F32[n];
     180        }
     181        if (scales) {
     182            scale = scales->data.F32[n];
     183        }
     184
     185        // Mask
     186        psImage *mask = NULL;
     187        if (masks != NULL) {
     188            mask = masks->data[n];
     189        }
     190
     191        // Iterate over the output image pixels
     192        for (int y = 0; y < outny; y++) {
     193            for (int x = 0; x < outnx; x++) {
     194                // Only transform those pixels requested
     195                if (!region || (region && region->data.U8[y][x])) {
     196                    // Transform!
     197                    sky->x = (double)x + 0.5;
     198                    sky->y = (double)y + 0.5;
     199                    (void)psPlaneTransformApply(detector, map, sky);
     200
     201                    // Change PS_INTERPOLATE_BILINEAR to best available technique.
     202                    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x,
     203                                                                              detector->y, mask, 1, NAN,
     204                                                                              PS_INTERPOLATE_BILINEAR);
     205                    if (error) {
     206                        // Error is actually the variance
     207                        outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
     208                                                                                                detector->x,
     209                                                                                                detector->y,
     210                                                                                                mask, 1, NAN);
     211                    }
     212
     213                    outImage->data.F32[y][x] = (outImage->data.F32[y][x] - offset) / scale;
     214                    outImage->data.F32[y][x] = outImage->data.F32[y][x] / SQUARE(scale);
     215
     216                } // Pixels of interest
     217
     218            }
     219        } // Iterating over output pixels
    220220
    221221    } // Iterating over images
Note: See TracChangeset for help on using the changeset viewer.