Commits

Renaud Miel authored 6886c290470
CAS-10682 Code cleanup and comments SDGrid::getXY
No tags

casa5/code/synthesis/MeasurementComponents/SDGrid.cc

Modified
1418 1418 return i;
1419 1419 }
1420 1420 }
1421 1421 }
1422 1422 // No match!
1423 1423 return -1;
1424 1424 }
1425 1425
1426 1426 Bool SDGrid::getXYPos(const VisBuffer& vb, Int row) {
1427 1427
1428 - Bool dointerp;
1429 - Bool nullPointingTable = false;
1430 - // const MSPointingColumns& act_mspc =
1431 - // interpolationRefFrame == InterpolationRefFrame::POINTING ?
1432 - // vb.msColumns().pointing() : convertedPointing_p ;
1433 - const MSPointingColumns& act_mspc = vb.msColumns().pointing();
1434 - nullPointingTable = (act_mspc.nrow() < 1);
1435 - Int pointIndex = -1;
1436 - if (!nullPointingTable) {
1437 - ///if(vb.newMS()) vb.newMS does not work well using msid
1438 - if (vb.msId() != msId_p) {
1439 - lastIndex_p = 0;
1440 - if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1441 - lastIndexPerAnt_p.resize(vb.numberAnt());
1442 - }
1443 - lastIndexPerAnt_p = 0;
1444 - msId_p = vb.msId();
1445 - lastAntID_p = -1;
1428 + // Check POINTING table.
1429 + // If the calling code is iterating over millions of rows,
1430 + // we'll do that check millions of times ...
1431 + const MSPointingColumns& act_mspc = vb.msColumns().pointing();
1432 + const auto nPointings = act_mspc.nrow();
1433 + Bool havePointings = (nPointings >= 1);
1434 +
1435 + // We'll need to call these many times, so let's call them once for good
1436 + const auto rowTime = vb.time()(row);
1437 + const auto rowTimeInterval = vb.timeInterval()(row);
1438 + const auto rowAntenna1 = vb.antenna1()(row);
1439 +
1440 + // 1. Try to find the index of a pointing recorded:
1441 + // - for the antenna of specified row,
1442 + // - at a time close enough to the time at which
1443 + // data of specified row was taken using that antenna
1444 + Int pointingIndex = -1;
1445 + if (havePointings) {
1446 + // if (vb.newMS()) vb.newMS does not work well using msid
1447 + // Note about above comment:
1448 + // - vb.newMS probably works well
1449 + // - but if the calling code is iterating over the rows of a subchunk
1450 + // vb.newMS returns true for all rows belonging to the first subchunk
1451 + // of the first chunk of a new MS.
1452 +
1453 + // ???
1454 + // What if vb changed since we were last called ?
1455 + // What if the calling code calls put and get, with different VisBuffers ?
1456 + if (vb.msId() != msId_p) {
1457 + lastIndex_p = 0; // no longer used
1458 + if (lastIndexPerAnt_p.nelements() < (size_t)vb.numberAnt()) {
1459 + lastIndexPerAnt_p.resize(vb.numberAnt());
1460 + }
1461 + lastIndexPerAnt_p = 0;
1462 + msId_p = vb.msId();
1463 + lastAntID_p = -1;
1464 + }
1465 +
1466 + // Try to locate a pointing verifying for specified row:
1467 + // | POINTING.TIME - MAIN.TIME | <= 0.5*(MAIN.INTERVAL + tolerance)
1468 +
1469 + // Try first using a tiny tolerance
1470 + constexpr Double useTinyTolerance = -1.0;
1471 + pointingIndex = getIndex(act_mspc, rowTime, useTinyTolerance , rowAntenna1);
1472 +
1473 + const Bool foundPointing = (pointingIndex >= 0);
1474 + if (not foundPointing) {
1475 + // Try again using tolerance = MAIN.INTERVAL
1476 + pointingIndex = getIndex(act_mspc, rowTime, rowTimeInterval, rowAntenna1);
1477 + }
1478 +
1479 + // Making the implicit type conversion explicit.
1480 + // Conversion is safe because it occurs only when pointingIndex >= 0.
1481 + const Bool foundValidPointing = (foundPointing and (static_cast<rownr_t>(pointingIndex) < nPointings));
1482 +
1483 + if (not foundValidPointing) {
1484 + ostringstream o;
1485 + o << "Failed to find pointing information for time "
1486 + << MVTime(rowTime/86400.0) << " : Omitting this point";
1487 +
1488 + logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
1489 +
1490 + return false;
1491 + }
1446 1492 }
1447 - pointIndex = getIndex(act_mspc, vb.time()(row), -1.0, vb.antenna1()(row));
1448 - //Try again to locate a pointing within the integration
1449 - if (pointIndex < 0)
1450 - pointIndex = getIndex(act_mspc, vb.time()(row), vb.timeInterval()(row), vb.antenna1()(row));
1451 - }
1452 - if (!nullPointingTable && ((pointIndex < 0) || (pointIndex >= Int(act_mspc.time().nrow())))) {
1453 - ostringstream o;
1454 - o << "Failed to find pointing information for time " <<
1455 - MVTime(vb.time()(row)/86400.0) << ": Omitting this point";
1456 - logIO_p << LogIO::DEBUGGING << String(o) << LogIO::POST;
1457 - // logIO_p << String(o) << LogIO::POST;
1458 - return false;
1459 - }
1460 1493
1461 - dointerp = false;
1462 - if (!nullPointingTable && (vb.timeInterval()(row) < act_mspc.interval()(pointIndex))) {
1463 - dointerp = true;
1464 - if (!isSplineInterpolationReady) {
1465 - interpolator = new SDPosInterpolator(vb, pointingDirCol_p);
1466 - isSplineInterpolationReady = true;
1467 - } else {
1468 - if (!interpolator->inTimeRange(vb.time()(row), vb.antenna1()(row))) {
1469 - // setup spline interpolator for the current dataset (CAS-11261, 2018/5/22 WK)
1470 - delete interpolator;
1471 - interpolator = 0;
1472 - interpolator = new SDPosInterpolator(vb, pointingDirCol_p);
1473 - }
1494 + // 2. At this stage we have a valid pointingIndex.
1495 + // Decide now if we need to interpolate antenna's pointing direction
1496 + // at data-taking time:
1497 + // we'll do so when data is sampled faster than pointings are recorded
1498 + const auto pointingInterval = act_mspc.interval()(pointingIndex);
1499 + const auto needInterpolation = (rowTimeInterval < pointingInterval);
1500 +
1501 + // 3. Create interpolator if needed
1502 + Bool dointerp = false;
1503 + if (havePointings && needInterpolation) {
1504 + dointerp = true;
1505 + // Known points are the directions of the specified
1506 + // POINTING table column,
1507 + // relative to the reference frame of the POINTING table
1508 + if (!isSplineInterpolationReady) {
1509 + interpolator = new SDPosInterpolator(vb, pointingDirCol_p);
1510 + isSplineInterpolationReady = true;
1511 + } else {
1512 + if (not interpolator->inTimeRange(rowTime, rowAntenna1)) {
1513 + // setup spline interpolator for the current dataset (CAS-11261, 2018/5/22 WK)
1514 + delete interpolator;
1515 + interpolator = 0;
1516 + interpolator = new SDPosInterpolator(vb, pointingDirCol_p);
1517 + // Missing isSplineInterpolationReady = true; ?
1518 + }
1519 + }
1474 1520 }
1475 - }
1476 1521
1477 - MEpoch epoch(Quantity(vb.time()(row), "s"));
1478 - if (!pointingToImage) {
1479 - // Set the frame
1480 - MPosition pos;
1481 - lastAntID_p = vb.antenna1()(row);
1482 - pos = vb.msColumns().antenna().positionMeas()(lastAntID_p);
1483 - mFrame_p = MeasFrame(epoch, pos);
1484 - if (!nullPointingTable) {
1485 - if (dointerp) {
1486 - worldPosMeas = directionMeas(act_mspc, pointIndex, vb.time()(row));
1487 - } else {
1488 - worldPosMeas = directionMeas(act_mspc, pointIndex);
1489 - }
1522 + // 4. If it does not already exists, create the machine to convert pointings directions
1523 + // and update the frame holding the measurements for this row
1524 + const MEpoch rowEpoch(Quantity(rowTime, "s"));
1525 + if (not pointingToImage) {
1526 + // Set the frame
1527 + const MPosition rowAntenna1Position =
1528 + vb.msColumns().antenna().positionMeas()(rowAntenna1);
1529 +
1530 + mFrame_p = MeasFrame(rowEpoch, rowAntenna1Position);
1531 +
1532 + // Remember antenna id for next call,
1533 + // which may be done using a different VisBuffer ...
1534 + lastAntID_p = rowAntenna1;
1535 +
1536 + // Not clear why we compute directions at this stage
1537 + if (havePointings) {
1538 + worldPosMeas = dointerp ? directionMeas(act_mspc, pointingIndex, rowTime)
1539 + : directionMeas(act_mspc, pointingIndex);
1540 + } else {
1541 + // Without pointings, this sets the direction to the phase center
1542 + worldPosMeas = vb.direction1()(row);
1543 + }
1544 +
1545 + // Make a machine to convert from the worldPosMeas to the output
1546 + // Direction Measure type for the relevant frame
1547 + MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1548 + pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1549 + if (not pointingToImage) {
1550 + logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
1551 + }
1490 1552 } else {
1491 - worldPosMeas = vb.direction1()(row);
1492 - }
1553 + // Reset the frame
1554 + // Always reset epoch
1555 + mFrame_p.resetEpoch(rowEpoch);
1556 +
1557 + // Reset antenna position only if antenna changed since we were last called
1558 + if (lastAntID_p != rowAntenna1) {
1559 + // Debug messages
1560 + if (lastAntID_p == -1) {
1561 + // antenna ID is unset
1562 + logIO_p << LogIO::DEBUGGING
1563 + << "updating antenna position for conversion: new MS ID " << msId_p
1564 + << ", antenna ID " << rowAntenna1 << LogIO::POST;
1565 + } else {
1566 + logIO_p << LogIO::DEBUGGING
1567 + << "updating antenna position for conversion: MS ID " << msId_p
1568 + << ", last antenna ID " << lastAntID_p
1569 + << ", new antenna ID " << rowAntenna1 << LogIO::POST;
1570 + }
1571 + const MPosition rowAntenna1Position =
1572 + vb.msColumns().antenna().positionMeas()(rowAntenna1);
1493 1573
1494 - //worldPosMeas=directionMeas(act_mspc, pointIndex);
1495 - // Make a machine to convert from the worldPosMeas to the output
1496 - // Direction Measure type for the relevant frame
1497 - MDirection::Ref outRef(directionCoord.directionType(), mFrame_p);
1498 - pointingToImage = new MDirection::Convert(worldPosMeas, outRef);
1499 - if (!pointingToImage) {
1500 - logIO_p << "Cannot make direction conversion machine" << LogIO::EXCEPTION;
1501 - }
1574 + mFrame_p.resetPosition(rowAntenna1Position);
1502 1575
1503 - } else {
1504 - mFrame_p.resetEpoch(epoch);
1505 - if (lastAntID_p != vb.antenna1()(row)) {
1506 - if (lastAntID_p == -1) {
1507 - // antenna ID is unset
1508 - logIO_p << LogIO::DEBUGGING
1509 - << "update antenna position for conversion: new MS ID " << msId_p
1510 - << ", antenna ID " << vb.antenna1()(row) << LogIO::POST;
1511 - } else {
1512 - logIO_p << LogIO::DEBUGGING
1513 - << "update antenna position for conversion: MS ID " << msId_p
1514 - << ", last antenna ID " << lastAntID_p
1515 - << ", new antenna ID " << vb.antenna1()(row) << LogIO::POST;
1516 - }
1517 - MPosition pos;
1518 - lastAntID_p = vb.antenna1()(row);
1519 - pos = vb.msColumns().antenna().positionMeas()(lastAntID_p);
1520 - mFrame_p.resetPosition(pos);
1576 + // Remember antenna id for next call,
1577 + // which may be done using a different VisBuffer ...
1578 + lastAntID_p = rowAntenna1;
1579 + }
1521 1580 }
1522 - }
1523 1581
1524 - if (!nullPointingTable) {
1525 - if (dointerp) {
1526 - MDirection newdir = directionMeas(act_mspc, pointIndex, vb.time()(row));
1527 - worldPosMeas = (*pointingToImage)(newdir);
1528 - //Vector<Double> newdirv = newdir.getAngle("rad").getValue();
1529 - //cerr<<"dir0="<<newdirv(0)<<endl;
1530 -
1531 - //fprintf(pfile,"%.8f %.8f \n", newdirv(0), newdirv(1));
1532 - //printf("%lf %lf \n", newdirv(0), newdirv(1));
1582 + // 5. First: interpolate pointing direction if needed,
1583 + // Then: convert the result to image's reference frame
1584 + if (havePointings) {
1585 + if (dointerp) {
1586 + MDirection newdir = directionMeas(act_mspc, pointingIndex, rowTime);
1587 + worldPosMeas = (*pointingToImage)(newdir);
1588 +
1589 + // Debug stuff
1590 + //Vector<Double> newdirv = newdir.getAngle("rad").getValue();
1591 + //cerr<<"dir0="<<newdirv(0)<<endl;
1592 +
1593 + //fprintf(pfile,"%.8f %.8f \n", newdirv(0), newdirv(1));
1594 + //printf("%lf %lf \n", newdirv(0), newdirv(1));
1595 + } else {
1596 + worldPosMeas = (*pointingToImage)(directionMeas(act_mspc, pointingIndex));
1597 + }
1533 1598 } else {
1534 - worldPosMeas = (*pointingToImage)(directionMeas(act_mspc, pointIndex));
1599 + // Without pointings, this converts the direction of the phase center
1600 + worldPosMeas = (*pointingToImage)(vb.direction1()(row));
1535 1601 }
1536 - } else {
1537 - worldPosMeas = (*pointingToImage)(vb.direction1()(row));
1538 - }
1539 1602
1540 - Bool result = directionCoord.toPixel(xyPos, worldPosMeas);
1541 - if (!result) {
1542 - logIO_p << "Failed to find a pixel for pointing direction of "
1543 - << MVTime(worldPosMeas.getValue().getLong("rad")).string(MVTime::TIME) << ", " << MVAngle(worldPosMeas.getValue().getLat("rad")).string(MVAngle::ANGLE) << LogIO::WARN << LogIO::POST;
1544 - return false;
1545 - }
1603 + // 6. Convert world position coordinates to image pixel coordinates
1604 + Bool result = directionCoord.toPixel(xyPos, worldPosMeas);
1605 + if (!result) {
1606 + logIO_p << "Failed to find a pixel for pointing direction of "
1607 + << MVTime(worldPosMeas.getValue().getLong("rad")).string(MVTime::TIME)
1608 + << ", " << MVAngle(worldPosMeas.getValue().getLat("rad")).string(MVAngle::ANGLE)
1609 + << LogIO::WARN << LogIO::POST;
1610 + return false;
1611 + }
1546 1612
1547 - if ((pointingDirCol_p == "SOURCE_OFFSET") || (pointingDirCol_p == "POINTING_OFFSET")) {
1548 - //there is no sense to track in offset coordinates...hopefully the
1549 - //user set the image coords right
1550 - fixMovingSource_p = false;
1551 - }
1552 - if (fixMovingSource_p) {
1553 - if (xyPosMovingOrig_p.nelements() < 2) {
1554 - directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1613 + // 7. Handle moving sources
1614 + if ((pointingDirCol_p == "SOURCE_OFFSET") || (pointingDirCol_p == "POINTING_OFFSET")) {
1615 + // It makes no sense to track in offset coordinates...
1616 + // hopefully the user sets the image coords right
1617 + fixMovingSource_p = false;
1555 1618 }
1556 - //via HADEC or AZEL for parallax of near sources
1557 - MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1558 - MDirection tmphadec=MDirection::Convert(movingDir_p, outref1)();
1559 - MDirection actSourceDir=(*pointingToImage)(tmphadec);
1560 - Vector<Double> actPix;
1561 - directionCoord.toPixel(actPix, actSourceDir);
1562 1619
1563 - //cout << row << " scan " << vb.scan()(row) << "xyPos " << xyPos << " xyposmovorig " << xyPosMovingOrig_p << " actPix " << actPix << endl;
1620 + if (fixMovingSource_p) {
1621 + if (xyPosMovingOrig_p.nelements() < 2) {
1622 + directionCoord.toPixel(xyPosMovingOrig_p, firstMovingDir_p);
1623 + }
1624 + //via HADEC or AZEL for parallax of near sources
1625 + MDirection::Ref outref1(MDirection::AZEL, mFrame_p);
1626 + MDirection tmphadec = MDirection::Convert(movingDir_p, outref1)();
1627 + MDirection actSourceDir = (*pointingToImage)(tmphadec);
1628 + Vector<Double> actPix;
1629 + directionCoord.toPixel(actPix, actSourceDir);
1564 1630
1565 - xyPos=xyPos+xyPosMovingOrig_p-actPix;
1566 - }
1631 + //cout << row << " scan " << vb.scan()(row) << "xyPos " << xyPos
1632 + // << " xyposmovorig " << xyPosMovingOrig_p << " actPix " << actPix << endl;
1567 1633
1568 - return result;
1569 - // Convert to pixel coordinates
1634 + xyPos = xyPos + xyPosMovingOrig_p - actPix;
1635 + }
1636 +
1637 + return result;
1570 1638 }
1571 1639
1572 1640 MDirection SDGrid::directionMeas(const MSPointingColumns& mspc, const Int& index){
1573 1641 if (pointingDirCol_p == "TARGET") {
1574 1642 return mspc.targetMeas(index);
1575 1643 } else if (pointingDirCol_p == "POINTING_OFFSET") {
1576 1644 if (!mspc.pointingOffsetMeasCol().isNull()) {
1577 1645 return mspc.pointingOffsetMeas(index);
1578 1646 }
1579 1647 cerr << "No PONTING_OFFSET column in POINTING table" << endl;

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut