**
** ROUTINE: tsRunMergeOne
**
**
** Merge one Frames Pipeline Run against another. Return SH_SUCCESS
** if successful, else an error code. Only OK_RUN objects of type
** STAR, GALAXY, or UNK are matched.
**
** Assumes that nu is sufficiently small that cos(nu) ~ 1 and can thus
** be ignored. DCR corrections are applied using cosmic colors, since
** a final photometric calibration can not be guaranteed to exist for
** the run at this time.
**
**
**
**============================================================================
*/
RET_CODE
tsRunMergeOne(ooTVArray(ooHandle(arObject)) *handles
/* MODIFIED: array of object handles */,
TSCALIBOBJ *array /* MODIFIED: Array of calibrated objects */,
int nObjects /* IN: Number of objects to merge */,
double minNu /* IN: minimum nu in array (degrees) */,
double maxNu /* IN: maximum nu in array (degrees) */,
double node /* IN: ascending node of GC (degrees) */,
double incl /* IN: inclination of GC (degrees) */,
double equinox /* IN: equinox of GC coords (years) */,
ooHandle(arFramesPipeRun)& other
/* IN: Frames Pipe Run to merge against */,
double matchRadius /* IN: Match radius (arcsecs) */)
{
// Fetch the range of fields and calibrations for the other run
uint16 field0 = other->field0();
uint16 nFields = other->nFields();
uint8 camCol = other->camCol();
ooHandle(arAstromCalib) astromH = other->mergeAstrometricCalibration();
if (fabs(astromH->equinox() - equinox) > 0.000001)
{
shErrStackPush("tsRunMergeOne: mismatched equinoxes");
return SH_GENERIC_ERROR;
}
double node2 = astromH->node();
double incl2 = astromH->incl();
// All Frames Pipeline Runs are required when stuffed to have used r' as
// their reference filter. Thus, we always want the TRANS structure for
// the r' CCD.
int camRow = 1;
int ccd = 10 * camRow + camCol;
// Get number of non-overlap rows and columns on a frame.
int nRows =
other->imagingRun()->program()->ccd(camRow,other->camCol())->nRows();
// HARDWIRED: number of columns per frame.
int nCols = 2048;
// HARDWIRED: number of overlap rows per frame.
int overlapRows = 128;
// Want full size of the frame
nRows += overlapRows;
// Convert match radius to degrees
matchRadius *= at_asec2Deg;
double negMatchRadius = -matchRadius;
double matchRadiusSq = matchRadius * matchRadius;
// Range in great circle longitude covered by the array (degrees)
double minMu = array[0].mu;
double maxMu = array[nObjects-1].mu;
// Some parked scans have oddly defined great circles in which the scans
// can go through 0/360. Thus, to match-up, we'll need to force the "other"
// scan into the same range of mu as "this" scan. Thus, lets determine a
// good range of mu for "this" scan (just subtract 10 degrees from the
// minimum mu for the low mu in the range).
double muRangeLo = minMu - 10.;
double muRangeHi = muRangeLo + 360.;
// For each field ...
printf("run %d, camCol %d, rerun %d, fields %d to %d\n", other->run(),
other->camCol(), other->rerun(), field0, field0+nFields-1);
fflush(stdout);
for (int32 field = field0; field < field0 + nFields; field++)
{
// Get astrometric transformation for this field
TRANS trans = astromH->astrotoolsTrans(ccd, field);
// Get great circle coordinates of the field center (degrees).
double centerMu = 0, centerNu = 0;
atTransApply(&trans, 'r', nRows / 2., 0., nCols / 2., 0., NULL, NULL,
¢erMu, NULL, ¢erNu,NULL);
// If first field, measure the field radius (degrees).
double fieldRadius;
if (field == field0)
{
double cornerMu = 0, cornerNu = 0;
atTransApply(&trans, 'r', 0., 0., 0., 0., NULL, NULL,
&cornerMu, NULL, &cornerNu, NULL);
fieldRadius = arSphericalCoord(centerMu, centerNu).
distance(arSphericalCoord(cornerMu, cornerNu)) * at_asec2Deg;
}
// Convert center coordinates to target run great circle coordinates
// (degrees).
double ra, dec;
atGCToEq(centerMu, centerNu, &ra, &dec, node2, incl2);
atEqToGC(ra, dec, ¢erMu, ¢erNu, node, incl);
atBound(¢erMu, muRangeLo, muRangeHi);
// Skip this field if it lies outside range of mu or nu spanned by the
// target run.
if (centerMu + fieldRadius < minMu || centerMu - fieldRadius > maxMu)
continue;
if (centerNu + fieldRadius < minNu || centerNu - fieldRadius > maxNu)
continue;
printf("%4d", field);
fflush(stdout);
// Get iterator to objects for this field
ooHandle(arObjectList) objectList = other->objectList(field);
ooItr(arObject) itr;
if (objectList->objects(itr, oocUpdate) != oocSuccess)
{
shErrStackPush("tsRunMergeOne: error getting objects for field %d",
field);
return SH_GENERIC_ERROR;
}
// For each object ...
int first = 1;
while (itr.next())
{
// Only match OK_RUN objects.
if (! (itr->status() & AR_OBJECT_STATUS_OK_RUN)) continue;
// Only match stars, galaxies, and unknowns
arObjectType objType = itr->objc_type();
if (objType != AR_OBJECT_TYPE_STAR &&
objType != AR_OBJECT_TYPE_GALAXY &&
objType != AR_OBJECT_TYPE_UNK)
continue;
// Fetch calibrated great circle coordinates (degrees)
arSphericalCoord coord = itr->gcCoord(trans);
double mu = coord.longitude();
double nu = coord.latitude();
// Convert to target run great circle coordinates (degrees)
atGCToEq(mu, nu, &ra, &dec, node2, incl2);
atEqToGC(ra, dec, &mu, &nu, node, incl);
atBound(&mu, muRangeLo, muRangeHi);
// Look for match in calibrated array. We'll search based on the
// object position minus the match radius.
int index;
if (first)
{
// Do a binomial search for the first object. The binomial
// search returns an array index which is set to one below the
// lowest matching or greater than array element.
index = tsCalibObjArrayBinomialSearch(mu-matchRadius, array,
nObjects) + 1;
first = 0;
}
else
{
// Search with correlated values around the last index.
index = tsCalibObjArrayCorrelatedSearch(mu-matchRadius, array,
nObjects, index-1) + 1;
}
// Skip if all objects in the array are at a smaller mu.
if (index == nObjects) continue;
// Now check each object until we're beyond the match radius or we
// run out of the array.
int closestIndex = -1;
double closestDistance = 1000. * matchRadiusSq;
double muHi = mu + matchRadius;
for (; index < nObjects; index++)
{
// Break out if we're beyond the match radius
if (array[index].mu > muHi) break;
double nuDiff = nu - array[index].nu;
if (nuDiff > matchRadius || nuDiff < negMatchRadius) continue;
double muDiff = mu - array[index].mu;
double distance = muDiff * muDiff + nuDiff * nuDiff;
if (distance <= matchRadiusSq)
{
// This matches. Save it if its the closest matching one
// thus far.
if (distance < closestDistance)
{
closestIndex = index;
closestDistance = distance;
}
}
}
// If we had a match, link the two objects.
if (closestIndex != -1)
(*handles)[array[closestIndex].index]->addMatches(itr);
}
}
printf("\n");
fflush(stdout);
// Return.
return SH_SUCCESS;
}
// Binomial search by mu on an array of TSCALIBOBJ. Returns the index such
// that the input mu is in the range mu(index) to mu(index+1);
// the input mu is definitely greater than mu(index), but could be
// equal to mu(index+1).
int tsCalibObjArrayBinomialSearch(double mu, TSCALIBOBJ *array, int nObjects)
{
int lo = -1;
int hi = nObjects;
while (hi - lo > 1)
{
int mid = (hi + lo) / 2;
if (mu > array[mid].mu)
lo = mid;
else
hi = mid;
}
return lo;
}
// Binomial search by row on an array of TSOBJ. Returns the index such
// that the input row is in the range row(index) to row(index+1);
// the input row is definitely greater than row(index), but could be
// equal to row(index+1).
int tsObjArrayBinomialSearch(double row, TSOBJ *array, int nObjects)
{
int lo = -1;
int hi = nObjects;
while (hi - lo > 1)
{
int mid = (hi + lo) / 2;
if (row > array[mid].row)
lo = mid;
else
hi = mid;
}
return lo;
}
// Correlated search by mu on an array of TSCALIBOBJ. Returns the index such
// that the input mu is in the range mu(index) to mu(index+1);
// the input mu is definitely greater than mu(index), but could be
// equal to mu(index+1).
int tsCalibObjArrayCorrelatedSearch(double mu, TSCALIBOBJ *array, int nObjects,
int guess)
{
if (guess <= -1)
{
if (mu <= array[0].mu)
return -1;
else
guess = 0;
}
else if (guess >= nObjects)
{
if (mu > array[nObjects - 1].mu)
return nObjects - 1;
else
guess = nObjects - 1;
}
int lo, hi;
int incr = 1;
if (mu > array[guess].mu)
{
if (guess == nObjects - 1) return guess;
lo = guess;
hi = guess + incr;
while (mu > array[hi].mu)
{
lo = hi;
incr += incr;
hi += incr;
if (hi >= nObjects)
{
hi = nObjects;
break;
}
}
}
else
{
if (guess == 0) return -1;
hi = guess;
lo = guess - incr;
while (mu <= array[lo].mu)
{
hi = lo;
incr += incr;
lo -= incr;
if (lo < 0)
{
lo = -1;
break;
}
}
}
while (hi - lo > 1)
{
int mid = (hi + lo) / 2;
if (mu > array[mid].mu)
lo = mid;
else
hi = mid;
}
return lo;
}