**
** ROUTINE: tsSegmentResolve
**
**
** Resolve all objects within the specified range of great circle
** longitude (startMu <= mu < endMu) in a single
** Frames Pipeline Run, creating a new segment in the operational
** databases. The type of the segment (PRIMARY, SECONDARY, BAD), as
** specified by the argument status), determines how objects
** are resolved. For SECONDARY and BAD segments, all objects whose
** flag OK_RUN is set and which lie within
** the range of mu are considered to belong to the segment. For PRIMARY
** segments, in addition to the above, if there are abutting PRIMARY
** segments within the scanline
** (argments startMatch and endMatch), then we must
** guarantee a unique and complete set of resolved obects in the border
** region (of radius matchRadius) where the segments meet.
** This is done by removing from the segment any object with a matching
** object which belongs to an abutting PRIMARY segment, and adding to the
** segment any object outside the mu limits but within
** matchRadius which has a matching object in a Frames Pipeline
** Run containing one of the abutting PRIMARY segments but which does not
** belong to that segment. All objects belonging to a segment have
** their RESOLVED flag set, as well as the PSEGMENT flag if the segment
** is a PRIMARY segment. Additionally, if
** the object belongs to the first field in the segment, then its
** FIRST_FIELD flag is set to distinguish objects lying in
** fields shared by two segments.
**
** If the segment is BAD, nothing else is done. If it is a SECONDARY
** segment, then all RESOLVED objects (of class STAR, GALAXY, UNK, or SKY
** only) have their SECONDARY flag set.
**
** Considerably more is done for PRIMARY segments. If the object lies
** within the primary range of great circle latitude for the scanline,
** then its OK_SCANLINE flag is set. A resolution process similar to
** that performed in mu is also done here to guarantee a set of unique
** and complete detections in the border region between PRIMARY segments
** in adjacent scanlines. For objects whose OK_SCANLINE flag is set,
** if its lies within the primary range of survey latitude for the
** stripe, then its OK_STRIPE flag is set. For objects whose OK_STRIPE
** flag is set, if it lies within the primary range of survey longitude
** for the stripe then its PRIMARY flag is set. Again, a resolution
** process
** is performed to guarantee a set of unique and complete detections
** around the border region between stripes (no latitude cut is applied
** to stripes in the southern survey). If the PRIMARY flag is not
** set, then the SECONDARY flag is (the PRIMARY or SECONDARY flag is set
** for objects of class STAR, GALAXY, UNK, and SKY only).
**
** For PRIMARY segments, it is assumed that there are no overlapping
** PRIMARY segments within the same scanline already defined in the
** database.
**
** The new database segment is appended to the input array of
** segments.
**
**
** RETURNS
** SH_SUCCESS --- Successful operation.
** SH_GENERIC_ERROR --- Error occurred. No segment is created,
** however object flags may have been
** set in the database, thus a database
** rollback is required.
**
**
**
**
**============================================================================
*/
RET_CODE
tsSegmentResolve(ooHandle(arFramesPipeRun) run, /*IN: handle to run */
int startMu, /* IN: starting GC longitude (arcsec) */
int endMu, /* IN: ending GC longitude (arcsec) */
double lambdaMin, /* IN: stripe minimum lambda (deg) */
double lambdaMax, /* IN: stripe maximum lambda (deg) */
double etaMin, /* IN: stripe minimim eta (deg) */
double etaMax, /* IN: stripe maximim eta (deg) */
TSSEGMENTSTATUS status, /* IN: segment status */
const ooHandle(ooContObj)& segmentsContainer,
/* IN: clustering directive when creating database segment */
ooTVArray(ooHandle(arSegment))& segments,
/* MODIFIED: append created segment here */
ooHandle(arSegment)* startMatch,
/* IN: abutting resolved segment at head (NULL if none) */
ooHandle(arSegment)* endMatch
/* IN: abutting resolved segment at tail (NULL if none) */)
{
// Get the astrometric calibration
ooHandle(arAstromCalib) astrom = run->mergeAstrometricCalibration();
if (! astrom.isValid())
{
shErrStackPush("tsSegmentResolve: run not merged");
return SH_GENERIC_ERROR;
}
if (fabs(astrom->equinox() - at_surveyEquinox) > 0.000001)
{
shErrStackPush("tsSegmentResolve: non-survey equinox for astrometric calibration");
return SH_GENERIC_ERROR;
}
double node = astrom->node();
double incl = astrom->incl();
// Get the minimum and maximum great circle latitude (nu) bounding this
// scanline.
int camCol = run->camCol();
int scanline = (camCol - 1) * 2 +
(run->imagingRun()->strip() == AR_STRIP_NORTH ? 1 : 0);
double nuMin = (scanline - 6) * at_scanSeparation;
double nuMax = nuMin + at_scanSeparation;
// DATA: delete following comments out lines when we no longer care about
// runs 94 and 125.
// KLUDGE: fix up 94 and 125 pointings
int currun = run->run();
if (currun == 94 || currun == 125) {
nuMin += 0.015;
nuMax += 0.015;
}
if (currun == 1033 || currun == 1056) {
nuMin -= 0.024;
nuMax -= 0.024;
}
int stripe = run->imagingRun()->stripe();
double stripeNode = 95.;
double stripeIncl = at_etaMin + stripe * at_stripeSeparation +
at_surveyCenterDec;
int southern = (stripe > 45);
if (southern) stripeIncl -= 180.;
// Get the matchRadius used when this run was merged. Actually, the
// match radius to use differs on which matched pair is being examined ---
// we should use the matchRadius used when the match was actually made,
// which depends on the order the runs were matched in. In effect, we
// assume that the same match radius was used for all, and thus we'll
// use this run's merge matchRadius regardless.
ooHandle(arProduct) mergeReport = run->merge(oocRead);
if (! mergeReport.isValid())
{
shErrStackPush("tsSegmentResolve: run hasn't been merged");
return SH_GENERIC_ERROR;
}
double matchRadius = atof(mergeReport->parameter("matchRadius"));
// As a sanity test, demand that matchRadius be greater than 0 arcsecs,
// but less than (arbitrarily) 5 arcsecs.
if (matchRadius <= 0. || matchRadius > 5.)
{
shErrStackPush("tsSegmentResolve: bad value for merge parameter 'matchRadius'");
return SH_GENERIC_ERROR;
}
// HARDWIRED: number of overlap rows per frame.
int overlapRows = 128;
// Get number of rows and cols for this run. We always work with r'
// positions, since all Frames Pipeline Runs were required when stuffed to
// have used r' as their reference filter.
int camRow = 1;
int ccd = 10 * camRow + camCol;
ooHandle(arCCDProgram) program =
run->imagingRun()->program()->ccd(camRow, camCol);
int rowsPerField = program->nRows();
// HARDWIRED: number of columns per frame.
int colsPerField = 2048;
// First and last rows in a field which are considered in the primary
// range for that field.
double firstRow = overlapRows / 2.;
double lastRow = rowsPerField + overlapRows / 2.;
// Find the first and last fields
int field0 = run->field0();
int field1 = field0 + run->nFields() - 1;
// Copy endpoint info to arrays so we can loop over endpoints
int endPointField[2];
int endPointMu[2];
ooHandle(arSegment) *endPointSegment[2];
ooHandle(ooContObj) endPointCont[2];
endPointMu[0] = startMu;
endPointMu[1] = endMu;
endPointSegment[0] = startMatch;
endPointSegment[1] = endMatch;
// For each endpoint ...
for (int i = 0; i < 2; i++)
{
// Find field containing endpoint
for (int field = field0; field <= field1; field++)
{
if (endPointMu[i] <=
tsRowToMu(astrom, camCol, field, lastRow - ROW_BUFFER,
colsPerField, stripeNode, stripeIncl, southern))
break;
}
if (field > field1)
{
shErrStackPush("tsSegmentResolve: mu range beyond run");
return SH_GENERIC_ERROR;
}
if (field == field0 &&
endPointMu[i] < tsRowToMu(astrom, camCol, field,
firstRow + ROW_BUFFER, colsPerField,
stripeNode, stripeIncl, southern))
{
shErrStackPush("tsSegmentResolve: mu range beyond run");
return SH_GENERIC_ERROR;
}
endPointField[i] = field;
// Find container for matching field in abutting previously resolved
// segment.
if (status == TS_PRIMARY && endPointSegment[i] != NULL)
{
int segField;
if (i == 0)
segField = (*endPointSegment[i])->field0() +
(*endPointSegment[i])->nFields() - 1;
else
segField = (*endPointSegment[i])->field0();
endPointCont[i] =
(*endPointSegment[i])->objectList(segField)->objectCont();
}
}
// Minimize number of fields in situations where we are right on a border.
if (endPointField[0] < field1)
if (startMu >= tsRowToMu(astrom, camCol, endPointField[0]+1,
firstRow + ROW_BUFFER, colsPerField, stripeNode,
stripeIncl, southern))
endPointField[0]++;
int nFields = endPointField[1] - endPointField[0] + 1;
// For each field ...
for (int field = endPointField[0]; field <= endPointField[1]; field++)
{
printf("%4d", field);
fflush(stdout);
// Get astrometric transformation for this field (r' chip)
TRANS trans = astrom->astrotoolsTrans(ccd, field);
// Initialize iterator over all objects in the field
ooHandle(arObjectList) objectList = run->objectList(field);
ooItr(arObject) itr;
if (objectList->objects(itr, oocUpdate) != oocSuccess)
{
shErrStackPush("tsSegmentResolve: couldn't initialize iterator over objects for field %d", field);
return SH_GENERIC_ERROR;
}
// For each object ...
while (itr.next())
{
// Skip objects which have already been resolved
int currentStatus = itr->status();
if (currentStatus & AR_OBJECT_STATUS_RESOLVED) continue;
// Skip objects which are not OK_RUN
if (! (currentStatus & AR_OBJECT_STATUS_OK_RUN)) continue;
// Only objects of class STAR, GALAXY, UNK, and SKY will be classed
// as either SECONDARY or PRIMARY. Only objects of class STAR,
// GALAXY, or UNK were matched.
int objType = itr->objc_type();
int starGal = (objType == AR_OBJECT_TYPE_STAR ||
objType == AR_OBJECT_TYPE_GALAXY ||
objType == AR_OBJECT_TYPE_UNK);
int inSurvey = (starGal || objType == AR_OBJECT_TYPE_SKY);
// Get great circle coordinates
arSphericalCoord coord = itr->gcCoord(trans);
double mu = coord.longitude();
double nu = coord.latitude();
// Parked scans don't scan the stripe great circle, so they need to
// be converted.
if (fabs(node - stripeNode) > 0.00001 ||
fabs(incl - stripeIncl) > 0.00001)
{
double ra, dec;
atGCToEq(mu, nu, &ra, &dec, node, incl);
atEqToGC(ra, dec, &mu, &nu, stripeNode, stripeIncl);
}
tsMuInRange(&mu, southern);
// Does this fall within the segment's limits in mu?
int ok = (mu * at_deg2Asec >= startMu && mu * at_deg2Asec < endMu);
// If this is not a primary segment, then set the object's status
// and continue to next object.
if (status != TS_PRIMARY)
{
if (ok)
{
int newStatus = AR_OBJECT_STATUS_RESOLVED;
if (field == endPointField[0])
newStatus |= AR_OBJECT_STATUS_FIRST_FIELD;
if (status == TS_SECONDARY && inSurvey)
newStatus |= AR_OBJECT_STATUS_SECONDARY;
itr->setStatus(currentStatus | newStatus);
}
continue;
}
// This is a primary segment. If it abuts another
// primary segment, then we need to guarantee
// complete and unique detections within a region of radius
// "matchRadius" around the border between the two. Lets first
// determine whether we need to check this. Only objects of class
// STAR, GALAXY, or UNK were matched, so we need check borders only
// for those.
int checkMu[2];
checkMu[0] = 0;
checkMu[1] = 0;
if (starGal)
for (int k = 0; k < 2; k++)
if (field == endPointField[k] && endPointSegment[k] != NULL)
checkMu[k] =
(fabs(mu * at_deg2Asec - endPointMu[k]) <= matchRadius);
// Skip to next object if this one is definitely outside range in mu,
// and thus definitely does not belong to this segment.
if (! ok && ! checkMu[0] && ! checkMu[1]) continue;
// Is this object in primary nu range for this scanline?
int okScanline = (nu >= nuMin && nu < nuMax);
int checkNuMin = (starGal && scanline > 0 &&
fabs(nu - nuMin) * at_deg2Asec <= matchRadius);
int checkNuMax = (starGal && scanline < 11 &&
fabs(nu - nuMax) * at_deg2Asec <= matchRadius);
// Is it in primary eta and lambda range for this stripe? Only
// check if its possibly in the primary nu range for the scanline.
int okStripe = 0;
int okLambda = 0;
int checkEtaMin = 0;
int checkEtaMax = 0;
if (okScanline || checkNuMin || checkNuMax)
{
double lambda, eta;
atGCToSurvey(mu, nu, stripeNode, stripeIncl, &lambda, &eta);
okLambda = (lambda >= lambdaMin && lambda < lambdaMax);
if (southern)
{
// Southern stripes don't have eta cuts.
okStripe = 1;
}
else
{
// Northern stripes do.
okStripe = (eta >= etaMin && eta < etaMax);
if (starGal)
{
checkEtaMin =
(fabs(eta-etaMin)*at_deg2Asec <= matchRadius);
checkEtaMax =
(fabs(eta-etaMax)*at_deg2Asec <= matchRadius);
}
}
}
// Are we on any borders such that we need to check the matching
// objects?
if (checkMu[0] || checkMu[1] || checkNuMin || checkNuMax ||
checkEtaMin || checkEtaMax)
{
// This is so we can turn off testing once we positively ruled
// an object out of a scanline or stripe.
int testScanline = checkNuMin || checkNuMax;
int testStripe = checkEtaMin || checkEtaMax;
// Initialize iterator over matching objects
ooItr(arObject) matches;
if (itr->matches(matches, oocRead) == oocError)
{
shErrStackPush("tsSegmentResolve: couldn't initialize iterator over matching objects for field %d", field);
return SH_GENERIC_ERROR;
}
// For each matching object ...
while (matches.next())
{
// Get matching object's container, status, stripe, and
// scanline.
int matchStatus = matches->status();
ooHandle(ooContObj) contH;
matches.containedIn(contH);
ooHandle(arFramesPipeRun) matchRun;
((ooHandle(arObjectCont)&)contH)->objectList()->
framesPipelineRun(matchRun, oocRead);
int matchStripe = matchRun->imagingRun()->stripe();
int matchScanline = (matchRun->camCol() - 1) * 2 +
(matchRun->imagingRun()->strip() == AR_STRIP_NORTH ? 1:0);
// Check mu borders. No need to test other matches once a
// match to a previously resolved object in a primary
// segment is found.
if (checkMu[0])
{
if (contH == endPointCont[0])
{
if (matchStatus & AR_OBJECT_STATUS_PSEGMENT)
{
ok = 0;
break;
}
else
ok = 1;
}
}
if (checkMu[1])
{
if (contH == endPointCont[1])
{
if (matchStatus & AR_OBJECT_STATUS_PSEGMENT)
{
ok = 0;
break;
}
else
ok = 1;
}
}
// Check nu and eta borders, but only if matching object
// belongs to a primary segment. Any match to an
// OK_SCANLINE/OK_STRIPE object forces this one to be
// not OK_SCANLINE/OK_STRIPE.
if (matchStatus & AR_OBJECT_STATUS_PSEGMENT)
{
if (testScanline && (matchStripe == stripe))
{
if ((checkNuMin && (scanline-matchScanline == 1)) ||
(checkNuMax && (matchScanline-scanline == 1)))
{
if (matchStatus & AR_OBJECT_STATUS_OK_SCANLINE)
{
okScanline = 0;
testScanline = 0;
}
else
okScanline = 1;
}
}
if (testStripe &&
((checkEtaMin && (stripe - matchStripe == 1)) ||
(checkEtaMax && (matchStripe - stripe == 1))))
{
if (matchStatus & AR_OBJECT_STATUS_OK_STRIPE)
{
okStripe = 0;
testStripe = 0;
}
else
okStripe = 1;
}
}
}
}
// If outside mu range, skip to next object;
if (! ok) continue;
// Set object status.
int newStatus =
AR_OBJECT_STATUS_RESOLVED | AR_OBJECT_STATUS_PSEGMENT;
if (field == endPointField[0])
newStatus |= AR_OBJECT_STATUS_FIRST_FIELD;
if (okScanline)
{
newStatus |= AR_OBJECT_STATUS_OK_SCANLINE;
if (okStripe)
{
newStatus |= AR_OBJECT_STATUS_OK_STRIPE;
if (inSurvey)
{
if (okLambda)
newStatus |= AR_OBJECT_STATUS_PRIMARY;
else
newStatus |= AR_OBJECT_STATUS_SECONDARY;
}
}
else
{
if (inSurvey) newStatus |= AR_OBJECT_STATUS_SECONDARY;
}
}
else
{
if (inSurvey) newStatus |= AR_OBJECT_STATUS_SECONDARY;
}
itr->setStatus(currentStatus | newStatus);
}
}
printf("\n");
// Create new database segment and append it to list of segments
ooHandle(arSegment) segment = new(segmentsContainer)
arSegment(run, startMu, endMu, endPointField[0], nFields,
status == TS_PRIMARY ? AR_SEGMENT_STATUS_PRIMARY :
(status == TS_SECONDARY ? AR_SEGMENT_STATUS_SECONDARY :
AR_SEGMENT_STATUS_BAD));
if (! segment.isValid())
{
shErrStackPush("tsSegmentResolve: couldn't create database segment");
return SH_GENERIC_ERROR;
}
segments.extend(segment);
return SH_SUCCESS;
}