IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17047


Ignore:
Timestamp:
Mar 18, 2008, 1:55:30 PM (18 years ago)
Author:
bills
Message:

Don't need to build a mosaic if there is only one cell in the chip.
If ROI is not completely contained on the chip extract the pixels that do
overlap and setting non-existent pixels to zero. (rather than failing as
we used to do)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/pstamp/src/ppstampMakeStamp.c

    r16593 r17047  
    9898}
    9999
    100 #ifdef notdef
    101 // update the ROI to account for any change in the coordinate system in the
    102 // mosaic process
    103 // THERE isn't any change. This is something I was trying to make the megacam
    104 // transforms come out correctly.
    105 static bool updateROI(ppstampOptions *options, pmFPAfile *mosaic, pmChip *mosaicChip)
    106 {
    107     // set up the astrometry
    108     pmHDU *hdu = pmHDUFromChip(mosaicChip);
    109     PS_ASSERT_PTR_NON_NULL(hdu, 1)
    110     PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
    111 
    112 
    113     if (!pmAstromReadWCS(mosaic->fpa, mosaicChip, hdu->header, 1.0)) {
    114         psError(PS_ERR_UNKNOWN, false, "Failed to Read WCS from mosaic, cannot proceed\n");
    115         return false;
    116     }
    117     pmAstromObj *center = pmAstromObjAlloc();
    118 
    119     center->sky->r = options->centerRA;
    120     center->sky->d = options->centerDEC;
    121     center->sky->rErr = 0;
    122     center->sky->dErr = 0;
    123 
    124     skyToChip(center, mosaic->fpa, mosaicChip);
    125 
    126     int width = options->dX;
    127     int height = options->dY;
    128 
    129     fprintf(stderr, "ROI before: %s\n", psRegionToString(options->roi));
    130 
    131     options->roi.x0 = center->chip->x - width/2;
    132     options->roi.x1 = options->roi.x0 + width;
    133     options->roi.y0 = center->chip->y - height/2;
    134     options->roi.y1 = options->roi.y0 + height;
    135 
    136     fprintf(stderr, "ROI after:  %s\n", psRegionToString(options->roi));
    137 
    138 
    139     return true;
    140 }
    141 #endif
     100// Extract pixels from an image. If part of the bounds of the region are
     101// partially outside of the input those pixels are set to zero.
     102psImage *extractStamp(psImage *image, psRegion region)
     103{
     104    int width  = region.x1 - region.x0;
     105    int height = region.y1 - region.y0;
     106
     107    if (width < 0) {
     108        return NULL;
     109    }
     110    if (height < 0) {
     111        return NULL;
     112    }
     113
     114    int srcX  = region.x0;
     115    int srcY  = region.y0;
     116    int lastX = region.x1;
     117    int lastY = region.y1;
     118    if (lastX > (image->col0 + image->numCols)) {
     119        lastX = image->col0 + image->numCols;
     120    }
     121    if (lastY > (image->row0 + image->numRows)) {
     122        lastY = image->row0 + image->numRows;
     123    }
     124       
     125    int leftBlank = 0;
     126    int dstX  = 0;
     127    if (srcX < image->col0) {
     128        leftBlank = image->col0 - srcX;
     129        dstX += leftBlank;
     130        srcX = 0;
     131    }
     132
     133    int copyWidth  = lastX - srcX;
     134    int rightBlank = width - copyWidth - leftBlank;
     135    if (copyWidth <= 0) {
     136        return NULL;
     137    }
     138
     139    psImage *output = psImageAlloc(width, height, image->type.type);
     140
     141    psElemType  type = image->type.type;
     142    psS32       elementSize =  PSELEMTYPE_SIZEOF(type);
     143    int         copySize = copyWidth * elementSize;
     144
     145    for (int dstY=0 ; dstY < height; srcY++, dstY++) {
     146        psU8 *pdst = output->data.U8[dstY];
     147
     148        if ((srcY < image->row0) || (srcY >= lastY)) {
     149            // zero the row
     150            memset(pdst, width*elementSize, 0);
     151        } else {
     152            psU8 *psrc = image->data.U8[srcY]  + (srcX * elementSize);
     153
     154            if (leftBlank) {
     155                memset(pdst, leftBlank*elementSize, 0);
     156            }
     157
     158            // copy copyWidth pixels from srcX,srcY to dstX, dstY
     159            pdst += dstX*elementSize;
     160            memcpy(pdst, psrc, copySize);
     161
     162            if (rightBlank) {
     163                pdst += copyWidth * elementSize;
     164                memset(pdst, rightBlank, 0);
     165            }
     166        }
     167    }
     168
     169
     170    return output;
     171}
    142172
    143173// Build the postage stamp output file
    … …  
    169199    psMetadataAddS32(target->concepts, PS_LIST_TAIL, "CELL.YBIN", PS_META_REPLACE, "Binning in y", 1);
    170200
    171 
    172     // we extract our postage stamp from a mosaic so that
    173     // pmChipMosaic can handle the tricky bits of parity, binning, etc.
    174     pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
    175     if (mosaic == NULL) {
    176         return false;
    177     }
    178 
    179     // one cell one chip
    180     pmFPAview *mosaicView = pmFPAviewAlloc(0);
    181     mosaicView->chip = 0;
    182     mosaicView->cell = 0;
     201    pmFPAfile *srcFile;
     202    pmFPAview *srcView = pmFPAviewAlloc(0);
     203    if (psArrayLength(inChip->cells) > 1) {
     204        // we extract our postage stamp from a mosaic so that
     205        // we don't have to manage extracting from multiple cells
     206        pmFPAfile *mosaic = ppstampBuildMosaic(config, input, view);
     207        if (mosaic == NULL) {
     208            return false;
     209        }
     210        srcFile = mosaic;
     211        srcView->chip = 0;
     212    } else {
     213        srcFile = input;
     214        *srcView = *view;
     215    }
     216
     217    // At this point we know we have one cell
     218    srcView->cell = 0;
    183219
    184220    pmReadout *readout;
    185     while ((readout = pmFPAviewNextReadout (mosaicView, mosaic->fpa, 1)) != NULL) {
     221    while ((readout = pmFPAviewNextReadout (srcView, srcFile->fpa, 1)) != NULL) {
    186222        if (!readout->data_exists) {
    187             psError(PS_ERR_UNKNOWN, false, "no data in mosaic readout!\n");
     223            psError(PS_ERR_UNKNOWN, false, "no data in input readout!\n");
    188224            continue;
    189225        }
    … …  
    200236        if (ppstampMegacamWorkaround) {
    201237            // the coordinates of the mosaic are shifted 32 pixels from the chip
     238            // TODO does this always apply? For example I doubt that applies to
     239            // skycells.
    202240            extractRegion.x0 -= 32;
    203241            extractRegion.x1 -= 32;
    204242        }
    205    
    206         psImage *subsetImage = psImageSubset(readout->image, extractRegion);
    207 
    208         if (subsetImage) {
    209             // make a copy since we're going to change the image's internals
    210             outReadout->image = psImageCopy(NULL, subsetImage, subsetImage->type.type);
    211             psFree(subsetImage);
    212             if (outReadout->image) {
    213                 // The image has inherited the subset's col0 and row0. We don't want this in
    214                 // our output so override the values.
    215                 // XXX put this into a function in psImage
    216                 // (What I really needed was a function that does
    217                 //     "create a copy of this portion of an image and return it as a new image with
    218                 //      origin 0,0 ")
    219                 P_PSIMAGE_SET_COL0(outReadout->image, 0);
    220                 P_PSIMAGE_SET_ROW0(outReadout->image, 0);
    221 
    222                 outReadout->data_exists = true;
    223                 outReadout->parent->data_exists = true;
    224                 outReadout->parent->parent->data_exists = true;
    225                 status = true;
    226             } else {
    227                 psError(PS_ERR_UNKNOWN, false, "Copying of readout image failed.\n");
    228                 status = false;
    229             }
    230         } else {
    231             psError(PS_ERR_UNKNOWN, false, "Image subset failed.\n");
     243
     244        outReadout->image = extractStamp(readout->image, extractRegion);
     245        if (!outReadout->image) {
     246            psError(PS_ERR_UNKNOWN, false, "failed to allocate create postage stamp image\n");
    232247            status = false;
    233         }
     248            break;
     249        }
     250
     251        outReadout->data_exists = true;
     252        outReadout->parent->data_exists = true;
     253        outReadout->parent->parent->data_exists = true;
     254        status = true;
     255
    234256        psFree(outReadout); // drop reference
    235257    }
    236258
    237     // XXXX remove this or make writing the mosaic to a file a valid option
    238     static bool writeMosaic = false;
    239     if (writeMosaic) {
    240         mosaic->save = true;
    241         mosaicView->cell = -1;
    242         pmFPAfileIOChecks(config, mosaicView, PM_FPA_AFTER);
    243     }
    244 
    245     psFree(mosaicView);
     259    psFree(srcView);
    246260
    247261    if (status) {
    … …  
    361375    options->roip.dY = options->roi.y1 - options->roi.y0;
    362376}
    363 
    364 #ifdef notdef
    365 // I'm not using this any more. We mosaic the chip before extracting the pixels
    366 // so we don't need to deal with cells at this level.
    367 
    368 // findCell
    369 // find cell on given chip that contains the region of interest
    370 // return the status of the overlap and if the cell completely contains the ROI
    371 // set a pointer to reference the cell
    372 
    373 static ppstampOverlap findCell(pmFPAview *view, ppstampOptions *options, pmFPAfile *input,
    374                                 pmAstromObj *center, pmCell **ppCell)
    375 {
    376     pmCell *cell;
    377 
    378     view->cell = -1;
    379     while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
    380         psRegion *cellExtent = pmCellExtent(cell);
    381         if (regionContainsPoint(cellExtent, center->chip)) {
    382             if (regionContainsRegion(cellExtent, &options->roi)) {
    383                 *ppCell = cell;
    384                 psFree(cellExtent);
    385                 return PPSTAMP_ON;
    386             } else {
    387                 fprintf(stderr, "ROI must be confined to single cell. Partial overlap cell %d\n",
    388                     view->cell);
    389                 fprintf(stderr, "Partial Overlap cell %d\n", view->cell);
    390                 fprintf(stderr, "  ROI:         %s\n", psRegionToString(options->roi));
    391                 fprintf(stderr, "  Cell Extent: %s\n", psRegionToString(*cellExtent));
    392                 psFree(cellExtent);
    393                 return PPSTAMP_PARTIALLY_ON;
    394             }
    395         }
    396         psFree(cellExtent);
    397     }
    398 
    399     // This shouldn't ever happen since ROI is on the chip, it must at least partially overlap
    400     // one of the cells
    401     psError(PS_ERR_PROGRAMMING, false, "findCell couldn't find a cell containing the center\n");
    402 
    403     return PPSTAMP_OFF;
    404 }
    405 #endif
    406 
    407 
    408377
    409378// findROI
    … …  
    491460            returnval = PPSTAMP_ON;
    492461        } else {
    493             fprintf(stderr, "Partial Overlap chip %s\n", chipName);
    494             fprintf(stderr, "ROI:         %s\n", psRegionToString(options->roi));
    495             fprintf(stderr, "Chip Extent: %s\n", psRegionToString(*chipBounds));
     462            psLogMsg("ppstampMakeStamp", 2, "Partial Overlap chip %s\n", chipName);
     463            psLogMsg("ppstampMakeStamp", 2, "ROI:         %s\n", psRegionToString(options->roi));
     464            psLogMsg("ppstampMakeStamp", 2, "Chip Extent: %s\n", psRegionToString(*chipBounds));
    496465
    497466            returnval = PPSTAMP_PARTIALLY_ON;
    498 
    499467        }
    500468    }
    … …  
    550518            break;
    551519        case PPSTAMP_ON:
     520        case PPSTAMP_PARTIALLY_ON:
    552521            returnval = makeStamp(config, options, input, chip, view);
    553522            allDone = true;
    554523            foundOverlap = true;
    555524            break;
    556         case PPSTAMP_PARTIALLY_ON:
    557             psError(PS_ERR_UNKNOWN, false, "Region of interest must be confined to a single chip\n");
    558             // XXX: perhaps print a helpful message about what coordinates would be valid
    559             returnval = false;
    560             allDone = true;
    561             foundOverlap = true;
    562             break;
    563525        case PSTAMP_ERROR:
    564526            returnval = false;
Note: See TracChangeset for help on using the changeset viewer.