IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25181


Ignore:
Timestamp:
Aug 24, 2009, 6:00:00 PM (17 years ago)
Author:
Paul Price
Message:

Wrote functions to merge detections from multiple skycells, and to write detections out. Compiles, not yet tested.

Location:
branches/pap_mops/ppMops/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/pap_mops/ppMops/src/ppMops.h

    r25162 r25181  
    4949    psVector *flags;                    // psphot flags
    5050    psVector *diffSkyfileId;            // Identifier for source image
     51    psVector *naxis1, *naxis2;          // Size of image
     52    psVector *mask;                     // Mask for detections
    5153} ppMopsDetections;
    5254
    5355ppMopsDetections *ppMopsDetectionsAlloc(long num);
     56
     57/// Copy a detection
     58bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index);
     59
     60/// Purge the detections list of masked detections
     61bool ppMopsDetectionsPurge(ppMopsDetections *detections);
     62
    5463
    5564/// Read detections

  • branches/pap_mops/ppMops/src/ppMopsDetections.c

    r25162 r25181  
    2828    psFree(det->flags);
    2929    psFree(det->diffSkyfileId);
     30    psFree(det->naxis1);
     31    psFree(det->naxis2);
     32    psFree(det->mask);
    3033    return;
    3134}
    … …  
    5659    det->magErr = psVectorAlloc(num, PS_TYPE_F32);
    5760    det->extended = psVectorAlloc(num, PS_TYPE_F32);
    58     det->angle = psVectorAlloc(num, PS_TYPE_F64);
    59     det->angleErr = psVectorAlloc(num, PS_TYPE_F64);
     61    det->angle = psVectorAlloc(num, PS_TYPE_F32);
     62    det->angleErr = psVectorAlloc(num, PS_TYPE_F32);
    6063    det->length = psVectorAlloc(num, PS_TYPE_F32);
    6164    det->lengthErr = psVectorAlloc(num, PS_TYPE_F32);
    6265    det->flags = psVectorAlloc(num, PS_TYPE_U32);
    6366    det->diffSkyfileId = psVectorAlloc(num, PS_TYPE_S64);
     67    det->naxis1 = psVectorAlloc(num, PS_TYPE_S32);
     68    det->naxis2 = psVectorAlloc(num, PS_TYPE_S32);
     69    det->mask = psVectorAlloc(num, PS_TYPE_U8);
    6470
    6571    return det;
    … …  
    6773
    6874
     75ppMopsDetections *ppMopsDetectionsRealloc(ppMopsDetections *det, long num)
     76{
     77    det->x = psVectorRealloc(det->x, num);
     78    det->y = psVectorRealloc(det->y, num);
     79    det->ra = psVectorRealloc(det->ra, num);
     80    det->dec = psVectorRealloc(det->dec, num);
     81    det->raErr = psVectorRealloc(det->raErr, num);
     82    det->decErr = psVectorRealloc(det->decErr, num);
     83    det->mag = psVectorRealloc(det->mag, num);
     84    det->magErr = psVectorRealloc(det->magErr, num);
     85    det->extended = psVectorRealloc(det->extended, num);
     86    det->angle = psVectorRealloc(det->angle, num);
     87    det->angleErr = psVectorRealloc(det->angleErr, num);
     88    det->length = psVectorRealloc(det->length, num);
     89    det->lengthErr = psVectorRealloc(det->lengthErr, num);
     90    det->flags = psVectorRealloc(det->flags, num);
     91    det->diffSkyfileId = psVectorRealloc(det->diffSkyfileId, num);
     92    det->naxis1 = psVectorRealloc(det->naxis1, num);
     93    det->naxis2 = psVectorRealloc(det->naxis2, num);
     94    det->mask = psVectorRealloc(det->mask, num);
     95
     96    return det;
     97}
     98
     99
     100bool ppMopsDetectionsAdd(ppMopsDetections *det, float x, float y, double ra, double dec,
     101                         double raErr, double decErr, float mag, float magErr, float extended,
     102                         float angle, float angleErr, float length, float lengthErr,
     103                         psU32 flags, psS64 diffSkyfileId, int naxis1, int naxis2)
     104{
     105    psVectorAppend(det->x, x);
     106    psVectorAppend(det->y, y);
     107    psVectorAppend(det->ra, ra);
     108    psVectorAppend(det->dec, dec);
     109    psVectorAppend(det->raErr, raErr);
     110    psVectorAppend(det->decErr, decErr);
     111    psVectorAppend(det->mag, mag);
     112    psVectorAppend(det->magErr, magErr);
     113    psVectorAppend(det->extended, extended);
     114    psVectorAppend(det->angle, angle);
     115    psVectorAppend(det->angleErr, angleErr);
     116    psVectorAppend(det->length, length);
     117    psVectorAppend(det->lengthErr, lengthErr);
     118    psVectorAppend(det->flags, flags);
     119    psVectorAppend(det->diffSkyfileId, diffSkyfileId);
     120    psVectorAppend(det->naxis1, naxis1);
     121    psVectorAppend(det->naxis2, naxis2);
     122    psVectorAppend(det->mask, 0);
     123    return true;
     124}
     125
     126
     127bool ppMopsDetectionsCopySingle(ppMopsDetections *target, const ppMopsDetections *source, long index)
     128{
     129    psVectorAppend(target->x, source->x->data.F32[index]);
     130    psVectorAppend(target->y, source->y->data.F32[index]);
     131    psVectorAppend(target->ra, source->ra->data.F64[index]);
     132    psVectorAppend(target->dec, source->dec->data.F64[index]);
     133    psVectorAppend(target->raErr, source->raErr->data.F64[index]);
     134    psVectorAppend(target->decErr, source->decErr->data.F64[index]);
     135    psVectorAppend(target->mag, source->mag->data.F32[index]);
     136    psVectorAppend(target->magErr, source->magErr->data.F32[index]);
     137    psVectorAppend(target->extended, source->extended->data.F32[index]);
     138    psVectorAppend(target->angle, source->angle->data.F32[index]);
     139    psVectorAppend(target->angleErr, source->angleErr->data.F32[index]);
     140    psVectorAppend(target->length, source->length->data.F32[index]);
     141    psVectorAppend(target->lengthErr, source->lengthErr->data.F32[index]);
     142    psVectorAppend(target->flags, source->flags->data.U32[index]);
     143    psVectorAppend(target->diffSkyfileId, source->diffSkyfileId->data.S64[index]);
     144    psVectorAppend(target->naxis1, source->naxis1->data.S32[index]);
     145    psVectorAppend(target->naxis2, source->naxis2->data.S32[index]);
     146    psVectorAppend(target->mask, 0);
     147    return true;
     148}
     149
     150
     151bool ppMopsDetectionsPurge(ppMopsDetections *det)
     152{
     153    long num = 0;
     154    for (long i = 0; i < det->num; i++) {
     155        if (!det->mask->data.U8[i]) {
     156            if (i == num) {
     157                // No need to copy
     158                num++;
     159                continue;
     160            }
     161            det->x->data.F32[num] = det->x->data.F32[i];
     162            det->y->data.F32[num] = det->y->data.F32[i];
     163            det->ra->data.F64[num] = det->ra->data.F64[i];
     164            det->dec->data.F64[num] = det->dec->data.F64[i];
     165            det->raErr->data.F64[num] = det->raErr->data.F64[i];
     166            det->decErr->data.F64[num] = det->decErr->data.F64[i];
     167            det->mag->data.F32[num] = det->mag->data.F32[i];
     168            det->magErr->data.F32[num] = det->magErr->data.F32[i];
     169            det->extended->data.F32[num] = det->extended->data.F32[i];
     170            det->angle->data.F32[num] = det->angle->data.F32[i];
     171            det->angleErr->data.F32[num] = det->angleErr->data.F32[i];
     172            det->length->data.F32[num] = det->length->data.F32[i];
     173            det->lengthErr->data.F32[num] = det->lengthErr->data.F32[i];
     174            det->flags->data.U32[num] = det->flags->data.U32[i];
     175            det->diffSkyfileId->data.S64[num] = det->diffSkyfileId->data.S64[i];
     176            det->naxis1->data.S32[num] = det->naxis1->data.S32[i];
     177            det->naxis2->data.S32[num] = det->naxis2->data.S32[i];
     178            det->mask->data.U8[num] = 0;
     179            num++;
     180        }
     181    }
     182    det->x->n = num;
     183    det->y->n = num;
     184    det->ra->n = num;
     185    det->dec->n = num;
     186    det->raErr->n = num;
     187    det->decErr->n = num;
     188    det->mag->n = num;
     189    det->magErr->n = num;
     190    det->extended->n = num;
     191    det->angle->n = num;
     192    det->angleErr->n = num;
     193    det->length->n = num;
     194    det->lengthErr->n = num;
     195    det->flags->n = num;
     196    det->diffSkyfileId->n = num;
     197    det->naxis1->n = num;
     198    det->naxis2->n = num;
     199    det->mask->n = num;
     200    det->num = num;
     201    return true;
     202}
     203

  • branches/pap_mops/ppMops/src/ppMopsMerge.c

    r25162 r25181  
    55#include <stdio.h>
    66#include <pslib.h>
     7#include <string.h>
    78
    89#include "ppMops.h"
    910
     11#define LEAF_SIZE 4                     // Size of leaf
     12#define MATCH_RADIUS 1.0/3600           // Matching radius
     13
     14// Get distance from detection to centre of image
     15static float mergeDistance(const ppMopsDetections *detections, // Detections of interest
     16                           long index                          // Index for source of interest
     17    )
     18{
     19    float dx = detections->x->data.F32[index] - detections->naxis1->data.S32[index] / 2.0;
     20    float dy = detections->y->data.F32[index] - detections->naxis2->data.S32[index] / 2.0;
     21    return PS_SQR(dx) + PS_SQR(dy);
     22}
     23
     24
    1025ppMopsDetections *ppMopsMerge(const psArray *detections)
    1126{
     27    PS_ASSERT_ARRAY_NON_NULL(detections, NULL);
    1228
    13     return NULL;
     29    ppMopsDetections *merged = psMemIncrRefCounter(detections->data[0]); // Merged list
     30    int num = 1;                                                         // Number of merged files
     31    for (int i = 1; i < detections->n; i++) {
     32        ppMopsDetections *det = detections->data[i]; // Detections of interest
     33        if (!det) {
     34            continue;
     35        }
     36        num++;
     37
     38        // XXX compare exposure properties
     39        if (strcmp(merged->raBoresight, det->raBoresight) != 0) {
     40            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure RA values differ: %s vs %s",
     41                    merged->raBoresight, det->raBoresight);
     42            return NULL;
     43        }
     44        if (strcmp(merged->decBoresight, det->decBoresight) != 0) {
     45            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure Dec values differ: %s vs %s",
     46                    merged->decBoresight, det->decBoresight);
     47            return NULL;
     48        }
     49        if (strcmp(merged->filter, det->filter) != 0) {
     50            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure filter values differ: %s vs %s",
     51                    merged->filter, det->filter);
     52            return NULL;
     53        }
     54
     55        if (merged->airmass != det->airmass) {
     56            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure airmass values differ: %f vs %f",
     57                    merged->airmass, det->airmass);
     58            return NULL;
     59        }
     60        if (merged->exptime != det->exptime) {
     61            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure exposure time values differ: %f vs %f",
     62                    merged->exptime, det->exptime);
     63            return NULL;
     64        }
     65        if (merged->posangle != det->posangle) {
     66            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure position angle values differ: %f vs %f",
     67                    merged->posangle, det->posangle);
     68            return NULL;
     69        }
     70        if (merged->alt != det->alt) {
     71            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure altitude values differ: %lf vs %lf",
     72                    merged->alt, det->alt);
     73            return NULL;
     74        }
     75        if (merged->az != det->az) {
     76            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure azimuth values differ: %lf vs %lf",
     77                    merged->az, det->az);
     78            return NULL;
     79        }
     80        if (merged->mjd != det->mjd) {
     81            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Exposure MJD values differ: %lf vs %lf",
     82                    merged->mjd, det->mjd);
     83            return NULL;
     84        }
     85
     86        merged->seeing += det->seeing;  // Taking average
     87
     88        psTree *tree = psTreePlant(2, LEAF_SIZE, PS_TREE_SPHERICAL, merged->ra, merged->dec); // kd tree
     89        if (!tree) {
     90            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate kd tree");
     91            psFree(merged);
     92            return NULL;
     93        }
     94
     95        psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of interest
     96        for (int j = 0; j < det->num; j++) {
     97            coords->data.F64[0] = det->ra->data.F64[j];
     98            coords->data.F64[1] = det->dec->data.F64[j];
     99            psVector *indices = psTreeAllWithin(tree, coords, MATCH_RADIUS); // Indices for matching sources
     100            if (!indices) {
     101                psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to search for matches");
     102                psFree(coords);
     103                psFree(tree);
     104                psFree(merged);
     105                return NULL;
     106            }
     107            if (indices->n == 0) {
     108                psFree(indices);
     109                continue;
     110            }
     111
     112            // Which one do we keep?
     113            float bestDistance = INFINITY; // Best distance to centre
     114            long bestIndex = -1;           // Index with best distance
     115            for (int k = 0; k < indices->n; k++) {
     116                long index = indices->data.S64[k]; // Index of point
     117                float distance = mergeDistance(merged, index); // Distance to centre of image
     118                if (distance < bestDistance) {
     119                    bestDistance = distance;
     120                    bestIndex = index;
     121                }
     122            }
     123
     124            float distance = mergeDistance(det, j); // Distance to centre of image
     125            if (distance < bestDistance) {
     126                // Blow away existing sources
     127                for (int k = 0; k < indices->n; k++) {
     128                    long index = indices->data.S64[k]; // Index of point
     129                    merged->mask->data.U8[index] = 0xFF;
     130                }
     131
     132                ppMopsDetectionsCopySingle(merged, det, j);
     133            }
     134            psFree(indices);
     135        }
     136
     137        psFree(tree);
     138        ppMopsDetectionsPurge(merged);
     139    }
     140
     141    merged->seeing /= num;
     142
     143    return merged;
    14144}
     145

  • branches/pap_mops/ppMops/src/ppMopsRead.c

    r25162 r25181  
    6060                             psMetadataLookupF32(NULL, header, "FWHM_MIN"));
    6161
     62        int naxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); // Number of columns
     63        int naxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); // Number of rows
     64
    6265        psFree(header);
    6366
    … …  
    8588            det->flags->data.F32[j] = psMetadataLookupU32(NULL, row, "FLAGS");
    8689            det->diffSkyfileId->data.F32[j] = diffSkyfileId;
     90            det->naxis1->data.S32[j] = naxis1;
     91            det->naxis2->data.S32[j] = naxis2;
     92            det->mask->data.U8[j] = 0;
    8793
    8894            // Calculate error in RA, Dec

  • branches/pap_mops/ppMops/src/ppMopsWrite.c

    r25162 r25181  
    88#include "ppMops.h"
    99
    10 bool ppMopsWrite(const ppMopsDetections *detections, const ppMopsArguments *args)
     10bool ppMopsWrite(const ppMopsDetections *det, const ppMopsArguments *args)
    1111{
     12    psFits *fits = psFitsOpen(args->output, "w"); // FITS file
     13    if (!fits) {
     14        psError(PS_ERR_IO, false, "Unable to open output file.");
     15        return false;
     16    }
    1217
     18    psMetadata *header = psMetadataAlloc(); // Header to write
     19    psString source = ppMopsSource(), version = ppMopsVersion();
     20    psMetadataAddStr(header, PS_LIST_TAIL, "SWSOURCE", 0, "Software source", source);
     21    psMetadataAddStr(header, PS_LIST_TAIL, "SWVERSN", 0, "Software version", version);
     22    ppMopsVersionHeader(header);
     23    psFree(source);
     24    psFree(version);
     25
     26    psMetadataAddStr(header, PS_LIST_TAIL, "EXP_NAME", 0, "Exposure name", args->exp_name);
     27    psMetadataAddS64(header, PS_LIST_TAIL, "EXP_ID", 0, "Exposure identifier", args->exp_id);
     28    psMetadataAddS64(header, PS_LIST_TAIL, "CHIP_ID", 0, "Chip stage identifier", args->chip_id);
     29    psMetadataAddS64(header, PS_LIST_TAIL, "CAM_ID", 0, "Cam stage identifier", args->cam_id);
     30    psMetadataAddS64(header, PS_LIST_TAIL, "FAKE_ID", 0, "Fake stage identifier", args->fake_id);
     31    psMetadataAddS64(header, PS_LIST_TAIL, "WARP_ID", 0, "Warp stage identifier", args->warp_id);
     32    psMetadataAddS64(header, PS_LIST_TAIL, "DIFF_ID", 0, "Diff stage identifier", args->diff_id);
     33    psMetadataAddBool(header, PS_LIST_TAIL, "DIFF_POS", 0, "Positive subtraction?", args->positive);
     34
     35    psMetadataAddF64(header, PS_LIST_TAIL, "MJD-OBS", 0, "MJD of exposure midpoint", det->mjd);
     36    psMetadataAddStr(header, PS_LIST_TAIL, "RA", 0, "Right Ascension of boresight", det->raBoresight);
     37    psMetadataAddStr(header, PS_LIST_TAIL, "DEC", 0, "Declination of boresight", det->decBoresight);
     38    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_ALT", 0, "Telescope altitude", det->alt);
     39    psMetadataAddF64(header, PS_LIST_TAIL, "TEL_AZ", 0, "Telescope azimuth", det->az);
     40    psMetadataAddF64(header, PS_LIST_TAIL, "EXPTIME", 0, "Exposure time (sec)", det->exptime);
     41    psMetadataAddF64(header, PS_LIST_TAIL, "ROTANGLE", 0, "Rotator position angle", det->posangle);
     42    psMetadataAddStr(header, PS_LIST_TAIL, "FILTER", 0, "Filter name", det->filter);
     43    psMetadataAddF32(header, PS_LIST_TAIL, "AIRMASS", 0, "Airmass of exposure", det->airmass);
     44    psMetadataAddStr(header, PS_LIST_TAIL, "OBSCODE", 0, "IAU Observatory code", OBSERVATORY_CODE);
     45    psMetadataAddF32(header, PS_LIST_TAIL, "SEEING", 0, "Mean seeing", det->seeing);
     46    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZP", 0, "Magnitude zero point", args->zp);
     47    psMetadataAddF32(header, PS_LIST_TAIL, "MAGZPERR", 0, "Error in magnitude zero point", args->zpErr);
     48    psMetadataAddF32(header, PS_LIST_TAIL, "ASTRORMS", 0, "RMS of astrometric fit", args->rmsAstrom);
     49
     50    if (det->num == 0) {
     51        // Write dummy table
     52        psMetadata *row = psMetadataAlloc(); // Output row
     53        psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", NAN);
     54        psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)", NAN);
     55        psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", NAN);
     56        psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)", NAN);
     57        psMetadataAddF64(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", NAN);
     58        psMetadataAddF64(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", NAN);
     59        psMetadataAddF64(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", NAN);
     60        psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)", NAN);
     61        psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)", NAN);
     62        psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)", NAN);
     63        psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)", NAN);
     64        psMetadataAddS32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", 0);
     65        psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile", 0);
     66        if (!psFitsWriteTableEmpty(fits, header, row, OUT_EXTNAME)) {
     67            psErrorStackPrint(stderr, "Unable to write empty table.");
     68            psFree(header);
     69            psFree(row);
     70            return false;
     71        }
     72        psFree(row);
     73    } else {
     74        psArray *table = psArrayAlloc(det->num); // Table to write
     75        for (long i = 0; i < det->num; i++) {
     76            psMetadata *row = psMetadataAlloc(); // Output row
     77            psMetadataAddF64(row, PS_LIST_TAIL, "RA", 0, "Right ascension (degrees)", det->ra->data.F64[i]);
     78            psMetadataAddF64(row, PS_LIST_TAIL, "RA_ERR", 0, "Right ascension error (degrees)",
     79                             det->raErr->data.F64[i]);
     80            psMetadataAddF64(row, PS_LIST_TAIL, "DEC", 0, "Declination (degrees)", det->dec->data.F64[i]);
     81            psMetadataAddF64(row, PS_LIST_TAIL, "DEC_ERR", 0, "Declination error (degrees)",
     82                             det->decErr->data.F64[i]);
     83            psMetadataAddF64(row, PS_LIST_TAIL, "MAG", 0, "Magnitude", det->mag->data.F64[i]);
     84            psMetadataAddF64(row, PS_LIST_TAIL, "MAG_ERR", 0, "Magnitude error", det->magErr->data.F64[i]);
     85            psMetadataAddF64(row, PS_LIST_TAIL, "STARPSF", 0, "EXT_NSIGMA", det->extended->data.F32[i]);
     86            psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE", 0, "Position angle of trail (degrees)",
     87                             det->angle->data.F32[i]);
     88            psMetadataAddF32(row, PS_LIST_TAIL, "ANGLE_ERR", 0, "Position angle error (degrees)",
     89                             det->angleErr->data.F32[i]);
     90            psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH", 0, "Length of trail (arcsec)",
     91                             det->length->data.F32[i]);
     92            psMetadataAddF32(row, PS_LIST_TAIL, "LENGTH_ERR", 0, "Length error (arcsec)",
     93                             det->lengthErr->data.F32[i]);
     94            psMetadataAddS32(row, PS_LIST_TAIL, "FLAGS", 0, "Detection bit flags", det->flags->data.S32[i]);
     95            psMetadataAddS64(row, PS_LIST_TAIL, "DIFF_SKYFILE_ID", 0, "Identifier for diff skyfile",
     96                             det->diffSkyfileId->data.S64[i]);
     97            table->data[i] = row;
     98        }
     99        if (!psFitsWriteTable(fits, header, table, OUT_EXTNAME)) {
     100            psErrorStackPrint(stderr, "Unable to write table.");
     101            psFree(header);
     102            psFree(table);
     103            return false;
     104        }
     105        psFree(table);
     106    }
     107
     108    psFree(header);
     109    psFitsClose(fits);
    13110    return true;
    14111}
Note: See TracChangeset for help on using the changeset viewer.