Commits

Federico Montesino Pouzols authored 92999617da4 Merge 
Merge remote-tracking branch 'origin/master' into bugfix/CAS-11397
937df492483
master

code/mstransform/MSTransform/MSTransformManager.cc

Modified
17 17 //# You should have received a copy of the GNU Lesser General Public
18 18 //# License along with this library; if not, write to the Free Software
19 19 //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 20 //# MA 02111-1307 USA
21 21 //# $Id: $
22 22 ​
23 23 #include <mstransform/MSTransform/MSTransformManager.h>
24 24 ​
25 25 #include <mstransform/TVI/PolAverageTVI.h>
26 26 ​
27 +#include <limits>
27 28 ​
28 29 using namespace casacore;
29 30 namespace casa { //# NAMESPACE CASA - BEGIN
30 31 ​
31 32 /////////////////////////////////////////////
32 33 ////// MS Transform Framework utilities /////
33 34 /////////////////////////////////////////////
34 35 namespace MSTransformations
35 36 {
36 37 Double wtToSigma(Double weight)
1644 1645 }
1645 1646 ​
1646 1647 propagateWeights(propagateWeights_p);
1647 1648 setChannelAverageKernel(weightmode_p);
1648 1649 setSmoothingKernel(smoothmode_p);
1649 1650 ​
1650 1651 ​
1651 1652 // Set Regridding kernel
1652 1653 if (fftShiftEnabled_p)
1653 1654 {
1654 - if (combinespws_p)
1655 - {
1656 - regridCoreComplex_p = &MSTransformManager::interpol1Dfftshift;
1657 - regridCoreFloat_p = &MSTransformManager::interpol1Dfftshift;
1658 - }
1659 - else
1660 - {
1661 - regridCoreComplex_p = &MSTransformManager::fftshift;
1662 - regridCoreFloat_p = &MSTransformManager::fftshift;
1663 - }
1655 + regridCoreComplex_p = &MSTransformManager::interpol1Dfftshift;
1656 + regridCoreFloat_p = &MSTransformManager::interpol1Dfftshift;
1664 1657 }
1665 1658 else
1666 1659 {
1667 1660 regridCoreComplex_p = &MSTransformManager::interpol1D;
1668 1661 regridCoreFloat_p = &MSTransformManager::interpol1D;
1669 1662 }
1670 1663 ​
1671 1664 //// Determine the frequency transformation methods to use ////
1672 1665 ​
1673 1666 // Drop channels with non-uniform width when doing only channel average
5833 5826 }
5834 5827 ​
5835 5828 Double newCentralFrequencyBeforeRegriddingAtReferenceTime = inputOutputSpwMap_p[spwIndex].first.CHAN_FREQ_aux[centralChan];
5836 5829 Double absoluteShift = newCentralFrequencyBeforeRegriddingAtCurrentTime - newCentralFrequencyBeforeRegriddingAtReferenceTime;
5837 5830 ​
5838 5831 Double chanWidth = inputOutputSpwMap_p[spwIndex].second.CHAN_FREQ[1] - inputOutputSpwMap_p[spwIndex].second.CHAN_FREQ[0];
5839 5832 Double bandwidth = inputOutputSpwMap_p[spwIndex].second.CHAN_FREQ[inputOutputSpwMap_p[spwIndex].second.NUM_CHAN-1] - inputOutputSpwMap_p[spwIndex].second.CHAN_FREQ[0];
5840 5833 bandwidth += chanWidth;
5841 5834 ​
5842 5835 fftShift_p = - absoluteShift / bandwidth;
5836 +​
5837 + ostringstream current;
5838 + current << setprecision(numeric_limits<double>::max_digits10)
5839 + << newCentralFrequencyBeforeRegriddingAtCurrentTime;
5840 + ostringstream reference;
5841 + reference << setprecision(numeric_limits<double>::max_digits10)
5842 + << newCentralFrequencyBeforeRegriddingAtReferenceTime;
5843 + logger_p << LogIO::NORMAL << LogOrigin("MSTransformManager",__FUNCTION__)
5844 + << "Using fftshift interpolation. The absolute shift is the "
5845 + << "new central frequency at current (input SPW) time - new "
5846 + << "central frequency "
5847 + << "at reference (output SPW) time\nAbsolute shift: "
5848 + << current
5849 + << " - " << reference
5850 + << " = " << absoluteShift << ", bandwidth " << bandwidth
5851 + << ", relative shift: " << fftShift_p << LogIO::POST;
5852 +​
5843 5853 }
5844 5854 else
5845 5855 {
5846 5856 for(uInt chan_idx=0; chan_idx<inputOutputSpwMap_p[spwIndex].first.CHAN_FREQ.size(); chan_idx++)
5847 5857 {
5848 5858 inputOutputSpwMap_p[spwIndex].first.CHAN_FREQ_aux[chan_idx] =
5849 5859 freqTransEngine_p(inputOutputSpwMap_p[spwIndex].first.CHAN_FREQ[chan_idx]).get(MSTransformations::Hz).getValue();
5850 5860 }
5851 5861 ​
5852 5862 /*
7116 7126 writeCube(inputDataCube,outputDataCol,rowRef);
7117 7127 if (outputFlagCol != NULL)
7118 7128 {
7119 7129 writeCube(vb->flagCube(),*outputFlagCol,rowRef);
7120 7130 }
7121 7131 ​
7122 7132 return;
7123 7133 }
7124 7134 ​
7125 7135 // -----------------------------------------------------------------------
7126 -//
7136 +// combine - for combinespws=True
7127 7137 // -----------------------------------------------------------------------
7128 7138 template <class T> void MSTransformManager::combineCubeOfData( vi::VisBuffer2 *vb,
7129 7139 RefRows &rowRef,
7130 7140 const Cube<T> &inputDataCube,
7131 7141 ArrayColumn<T> &outputDataCol,
7132 7142 ArrayColumn<Bool> *outputFlagCol,
7133 7143 const Cube<Float> &inputWeightCube)
7134 7144 {
7135 7145 // Write flag column too?
7136 7146 if (outputFlagCol != NULL)
7506 7516 // -----------------------------------------------------------------------
7507 7517 template <class T> void MSTransformManager::transformAndWriteCubeOfData( Int inputSpw,
7508 7518 RefRows &rowRef,
7509 7519 const Cube<T> &inputDataCube,
7510 7520 const Cube<Bool> &inputFlagsCube,
7511 7521 const Cube<Float> &inputWeightsCube,
7512 7522 IPosition &outputPlaneShape,
7513 7523 ArrayColumn<T> &outputDataCol,
7514 7524 ArrayColumn<Bool> *outputFlagCol)
7515 7525 {
7526 + logger_p << LogIO::DEBUG1 << LogOrigin("MSTransformManager",__FUNCTION__)
7527 + << "Shape of input data cube: " << inputDataCube.shape()
7528 + << ", output plane shape: " << outputPlaneShape
7529 + << LogIO::POST;
7530 +​
7516 7531 // Write flag column too?
7517 7532 if (outputFlagCol != NULL)
7518 7533 {
7519 7534 writeOutputFlagsPlaneSlices_p = &MSTransformManager::writeOutputFlagsPlaneSlices;
7520 7535 writeOutputFlagsPlaneReshapedSlices_p = &MSTransformManager::writeOutputFlagsPlaneReshapedSlices;
7521 7536 writeOutputFlagsPlane_p = &MSTransformManager::writeOutputFlagsPlane;
7522 7537 }
7523 7538 else
7524 7539 {
7525 7540 writeOutputFlagsPlaneSlices_p = &MSTransformManager::dontWriteOutputFlagsPlaneSlices;
7614 7629 }
7615 7630 ​
7616 7631 // -----------------------------------------------------------------------
7617 7632 //
7618 7633 // -----------------------------------------------------------------------
7619 7634 void MSTransformManager::setWeightsPlaneByReference( uInt inputRow,
7620 7635 const Cube<Float> &inputWeightsCube,
7621 7636 Matrix<Float> &inputWeightsPlane)
7622 7637 {
7623 7638 inputWeightsPlane = inputWeightsCube.xyPlane(inputRow);
7624 -​
7625 - return;
7626 7639 }
7627 7640 ​
7628 7641 // -----------------------------------------------------------------------
7629 7642 //
7630 7643 // -----------------------------------------------------------------------
7631 7644 template <class T> void MSTransformManager::transformAndWritePlaneOfData( Int inputSpw,
7632 7645 uInt row,
7633 7646 Matrix<T> &inputDataPlane,
7634 7647 Matrix<Bool> &inputFlagsPlane,
7635 7648 Matrix<Float> &inputWeightsPlane,
7662 7675 // Fill output stripes by reference
7663 7676 outputDataStripe.reference(outputDataPlane.row(corrIndex));
7664 7677 outputFlagsStripe.reference(outputFlagsPlane.row(corrIndex));
7665 7678 ​
7666 7679 transformStripeOfData(inputSpw,inputDataStripe,inputFlagsStripe,
7667 7680 inputWeightsStripe,outputDataStripe,outputFlagsStripe);
7668 7681 }
7669 7682 ​
7670 7683 // Write output planes
7671 7684 writeOutputPlanes(row,outputDataPlane,outputFlagsPlane,outputDataCol,*outputFlagCol);
7672 -​
7673 - return;
7674 7685 }
7675 7686 ​
7676 7687 // -----------------------------------------------------------------------
7677 7688 //
7678 7689 // -----------------------------------------------------------------------
7679 7690 void MSTransformManager::writeOutputPlanes( uInt row,
7680 7691 Matrix<Complex> &outputDataPlane,
7681 7692 Matrix<Bool> &outputFlagsPlane,
7682 7693 ArrayColumn<Complex> &outputDataCol,
7683 7694 ArrayColumn<Bool> &outputFlagCol)
7684 7695 {
7685 7696 (*this.*writeOutputPlanesComplex_p)(row,outputDataPlane,outputFlagsPlane,outputDataCol,outputFlagCol);
7686 - return;
7687 7697 }
7688 7698 ​
7689 7699 // -----------------------------------------------------------------------
7690 7700 //
7691 7701 // -----------------------------------------------------------------------
7692 7702 void MSTransformManager::writeOutputPlanes( uInt row,
7693 7703 Matrix<Float> &outputDataPlane,
7694 7704 Matrix<Bool> &outputFlagsPlane,
7695 7705 ArrayColumn<Float> &outputDataCol,
7696 7706 ArrayColumn<Bool> &outputFlagCol)
7697 7707 {
7698 7708 (*this.*writeOutputPlanesFloat_p)(row,outputDataPlane,outputFlagsPlane,outputDataCol,outputFlagCol);
7699 - return;
7700 7709 }
7701 7710 ​
7702 7711 // -----------------------------------------------------------------------
7703 7712 //
7704 7713 // -----------------------------------------------------------------------
7705 7714 void MSTransformManager::setOutputbuffer(Cube<Complex> *& dataBufferPointer,Cube<Bool> *& flagBufferPointer)
7706 7715 {
7707 7716 switch (dataBuffer_p)
7708 7717 {
7709 7718 case MSTransformations::visCube:
7912 7921 template <class T> void MSTransformManager::writeOutputPlanesInBlock( uInt row,
7913 7922 Matrix<T> &outputDataPlane,
7914 7923 Matrix<Bool> &outputFlagsPlane,
7915 7924 ArrayColumn<T> &outputDataCol,
7916 7925 ArrayColumn<Bool> &outputFlagCol)
7917 7926 {
7918 7927 IPosition outputPlaneShape = outputDataPlane.shape();
7919 7928 outputDataCol.setShape(row,outputPlaneShape);
7920 7929 outputDataCol.put(row, outputDataPlane);
7921 7930 (*this.*writeOutputFlagsPlane_p)(outputFlagsPlane,outputFlagCol, outputPlaneShape, row);
7922 -​
7923 - return;
7924 7931 }
7925 7932 ​
7926 7933 ​
7927 7934 // -----------------------------------------------------------------------
7928 7935 //
7929 7936 // -----------------------------------------------------------------------
7930 7937 void MSTransformManager::writeOutputFlagsPlane( Matrix<Bool> &outputPlane,
7931 7938 ArrayColumn<Bool> &outputCol,
7932 7939 IPosition &outputPlaneShape,
7933 7940 uInt &outputRow)
7934 7941 {
7935 7942 outputCol.setShape(outputRow,outputPlaneShape);
7936 7943 outputCol.put(outputRow, outputPlane);
7937 - return;
7938 7944 }
7939 7945 ​
7940 7946 // -----------------------------------------------------------------------
7941 7947 //
7942 7948 // -----------------------------------------------------------------------
7943 7949 template <class T> void MSTransformManager::writeOutputPlanesInSlices( uInt row,
7944 7950 Matrix<T> &outputDataPlane,
7945 7951 Matrix<Bool> &outputFlagsPlane,
7946 7952 ArrayColumn<T> &outputDataCol,
7947 7953 ArrayColumn<Bool> &outputFlagCol)
7961 7967 (*this.*writeOutputFlagsPlaneSlices_p)( outputFlagsPlane,outputFlagCol,
7962 7968 sliceX,sliceY,outputPlaneShape_i,outRow);
7963 7969 }
7964 7970 ​
7965 7971 uInt outRow = row+spw_i;
7966 7972 IPosition outputPlaneShape_tail(2,nCorrs,tailOfChansforLastSpw_p);
7967 7973 Slice sliceY(chansPerOutputSpw_p*spw_i,tailOfChansforLastSpw_p);
7968 7974 writeOutputPlaneReshapedSlices(outputDataPlane,outputDataCol,sliceX,sliceY,outputPlaneShape_tail,outRow);
7969 7975 (*this.*writeOutputFlagsPlaneReshapedSlices_p)( outputFlagsPlane,outputFlagCol,
7970 7976 sliceX,sliceY,outputPlaneShape_tail,outRow);
7971 -​
7972 - return;
7973 7977 }
7974 7978 ​
7975 7979 // -----------------------------------------------------------------------
7976 7980 //
7977 7981 // -----------------------------------------------------------------------
7978 7982 void MSTransformManager::writeOutputFlagsPlaneSlices( Matrix<Bool> &outputPlane,
7979 7983 ArrayColumn<Bool> &outputCol,
7980 7984 Slice &sliceX,
7981 7985 Slice &sliceY,
7982 7986 IPosition &outputPlaneShape,
7983 7987 uInt &outputRow)
7984 7988 {
7985 7989 writeOutputPlaneSlices(outputPlane,outputCol,sliceX,sliceY,outputPlaneShape,outputRow);
7986 - return;
7987 7990 }
7988 7991 ​
7989 7992 // -----------------------------------------------------------------------
7990 7993 //
7991 7994 // -----------------------------------------------------------------------
7992 7995 void MSTransformManager::writeOutputFlagsPlaneReshapedSlices( Matrix<Bool> &outputPlane,
7993 7996 ArrayColumn<Bool> &outputCol,
7994 7997 Slice &sliceX,
7995 7998 Slice &sliceY,
7996 7999 IPosition &outputPlaneShape,
7997 8000 uInt &outputRow)
7998 8001 {
7999 8002 writeOutputPlaneReshapedSlices(outputPlane,outputCol,sliceX,sliceY,outputPlaneShape,outputRow);
8000 - return;
8001 8003 }
8002 8004 ​
8003 8005 // -----------------------------------------------------------------------
8004 8006 //
8005 8007 // -----------------------------------------------------------------------
8006 8008 template <class T> void MSTransformManager::writeOutputPlaneSlices( Matrix<T> &outputPlane,
8007 8009 ArrayColumn<T> &outputCol,
8008 8010 Slice &sliceX,
8009 8011 Slice &sliceY,
8010 8012 IPosition &outputPlaneShape,
8011 8013 uInt &outputRow)
8012 8014 {
8013 8015 Matrix<T> outputPlane_i = outputPlane(sliceX,sliceY);
8014 8016 outputCol.setShape(outputRow,outputPlaneShape);
8015 8017 outputCol.put(outputRow, outputPlane_i);
8016 - return;
8017 8018 }
8018 8019 ​
8019 8020 // -----------------------------------------------------------------------
8020 8021 //
8021 8022 // -----------------------------------------------------------------------
8022 8023 template <class T> void MSTransformManager::writeOutputPlaneReshapedSlices( Matrix<T> &outputPlane,
8023 8024 ArrayColumn<T> &outputCol,
8024 8025 Slice &sliceX,
8025 8026 Slice &sliceY,
8026 8027 IPosition &outputPlaneShape,
8027 8028 uInt &outputRow)
8028 8029 {
8029 8030 Matrix<T> outputPlane_i = outputPlane(sliceX,sliceY);
8030 8031 outputPlane_i.resize(outputPlaneShape,true);
8031 8032 outputCol.setShape(outputRow,outputPlaneShape);
8032 8033 outputCol.put(outputRow, outputPlane_i);
8033 - return;
8034 8034 }
8035 8035 ​
8036 8036 // -----------------------------------------------------------------------
8037 8037 //
8038 8038 // -----------------------------------------------------------------------
8039 8039 void MSTransformManager::setWeightStripeByReference( uInt corrIndex,
8040 8040 Matrix<Float> &inputWeightsPlane,
8041 8041 Vector<Float> &inputWeightsStripe)
8042 8042 {
8043 8043 inputWeightsStripe.reference(inputWeightsPlane.row(corrIndex));
8044 - return;
8045 8044 }
8046 8045 ​
8047 8046 // -----------------------------------------------------------------------
8048 8047 //
8049 8048 // -----------------------------------------------------------------------
8050 8049 void MSTransformManager::transformStripeOfData(Int inputSpw,
8051 8050 const Vector<Complex> &inputDataStripe,
8052 8051 const Vector<Bool> &inputFlagsStripe,
8053 8052 const Vector<Float> &inputWeightsStripe,
8054 8053 Vector<Complex> &outputDataStripe,
8055 8054 Vector<Bool> &outputFlagsStripe)
8056 8055 {
8057 - (*this.*transformStripeOfDataComplex_p)( inputSpw,inputDataStripe,inputFlagsStripe,
8058 - inputWeightsStripe,outputDataStripe,outputFlagsStripe);
8059 - return;
8056 + auto shapeBefore = outputDataStripe.shape();
8057 + (*this.*transformStripeOfDataComplex_p)(inputSpw, inputDataStripe, inputFlagsStripe,
8058 + inputWeightsStripe, outputDataStripe,
8059 + outputFlagsStripe);
8060 + auto shapeAfter = outputDataStripe.shape();
8061 + if (shapeAfter != shapeBefore) {
8062 + ostringstream msg;
8063 + msg << "Shape of output complex data stripe changed after applying "
8064 + << "transformation. Output shape expected before transformation: "
8065 + << shapeBefore
8066 + << ". Output shape produced by transformation: " << shapeAfter;
8067 + logger_p << LogIO::DEBUG1 << LogOrigin("MSTransformManager",__FUNCTION__)
8068 + << LogIO::POST;
8069 + throw AipsError(msg.str());
8070 + }
8060 8071 }
8061 8072 ​
8062 8073 // -----------------------------------------------------------------------
8063 8074 //
8064 8075 // -----------------------------------------------------------------------
8065 8076 void MSTransformManager::transformStripeOfData(Int inputSpw,
8066 8077 const Vector<Float> &inputDataStripe,
8067 8078 const Vector<Bool> &inputFlagsStripe,
8068 8079 const Vector<Float> &inputWeightsStripe,
8069 8080 Vector<Float> &outputDataStripe,
8070 8081 Vector<Bool> &outputFlagsStripe)
8071 8082 {
8072 8083 (*this.*transformStripeOfDataFloat_p)( inputSpw,inputDataStripe,inputFlagsStripe,inputWeightsStripe,
8073 8084 outputDataStripe,outputFlagsStripe);
8074 - return;
8075 8085 }
8076 8086 ​
8077 8087 // -----------------------------------------------------------------------
8078 8088 //
8079 8089 // -----------------------------------------------------------------------
8080 8090 template <class T> void MSTransformManager::average(Int inputSpw,
8081 8091 const Vector<T> &inputDataStripe,
8082 8092 const Vector<Bool> &inputFlagsStripe,
8083 8093 const Vector<Float> &inputWeightsStripe,
8084 8094 Vector<T> &outputDataStripe,
8703 8713 Vector<Bool> &outputFlagsStripe)
8704 8714 {
8705 8715 ​
8706 8716 regridCore( inputSpw,
8707 8717 inputDataStripe,
8708 8718 inputFlagsStripe,
8709 8719 inputWeightsStripe,
8710 8720 outputDataStripe,
8711 8721 outputFlagsStripe);
8712 8722 ​
8713 - return;
8714 8723 }
8715 8724 ​
8716 8725 // -----------------------------------------------------------------------
8717 8726 //
8718 8727 // -----------------------------------------------------------------------
8719 8728 void MSTransformManager::regridCore(Int inputSpw,
8720 8729 const Vector<Complex> &inputDataStripe,
8721 8730 const Vector<Bool> &inputFlagsStripe,
8722 8731 const Vector<Float> &inputWeightsStripe,
8723 8732 Vector<Complex> &outputDataStripe,
8724 8733 Vector<Bool> &outputFlagsStripe)
8725 8734 {
8726 8735 ​
8727 8736 (*this.*regridCoreComplex_p)( inputSpw,
8728 8737 inputDataStripe,
8729 8738 inputFlagsStripe,
8730 8739 inputWeightsStripe,
8731 8740 outputDataStripe,
8732 8741 outputFlagsStripe);
8733 - return;
8734 8742 }
8735 8743 ​
8736 8744 // -----------------------------------------------------------------------
8737 8745 //
8738 8746 // -----------------------------------------------------------------------
8739 8747 void MSTransformManager::regridCore(Int inputSpw,
8740 8748 const Vector<Float> &inputDataStripe,
8741 8749 const Vector<Bool> &inputFlagsStripe,
8742 8750 const Vector<Float> &inputWeightsStripe,
8743 8751 Vector<Float> &outputDataStripe,
8744 8752 Vector<Bool> &outputFlagsStripe)
8745 8753 {
8746 8754 (*this.*regridCoreFloat_p)( inputSpw,
8747 8755 inputDataStripe,
8748 8756 inputFlagsStripe,
8749 8757 inputWeightsStripe,
8750 8758 outputDataStripe,
8751 8759 outputFlagsStripe);
8752 - return;
8753 8760 }
8754 8761 ​
8755 8762 // -----------------------------------------------------------------------
8756 8763 //
8757 8764 // -----------------------------------------------------------------------
8758 8765 void MSTransformManager::fftshift(Int ,
8759 8766 const Vector<Complex> &inputDataStripe,
8760 8767 const Vector<Bool> &inputFlagsStripe,
8761 8768 const Vector<Float> &,
8762 8769 Vector<Complex> &outputDataStripe,
8763 8770 Vector<Bool> &outputFlagsStripe)
8764 8771 {
8765 8772 fFFTServer_p.fftshift(outputDataStripe,
8766 8773 outputFlagsStripe,
8767 8774 (const Vector<Complex>)inputDataStripe,
8768 8775 (const Vector<Bool>)inputFlagsStripe,
8769 8776 (const uInt)0, // In vectors axis 0 is the only dimension
8770 8777 (const Double)fftShift_p,
8771 8778 false, // A good data point has its flag set to false
8772 8779 false);
8773 - return;
8774 8780 }
8775 8781 ​
8776 8782 // -----------------------------------------------------------------------
8777 8783 //
8778 8784 // -----------------------------------------------------------------------
8779 8785 void MSTransformManager::fftshift(Int ,
8780 8786 const Vector<Float> &inputDataStripe,
8781 8787 const Vector<Bool> &inputFlagsStripe,
8782 8788 const Vector<Float> &,
8783 8789 Vector<Float> &outputDataStripe,
8784 8790 Vector<Bool> &outputFlagsStripe)
8785 8791 {
8786 8792 fFFTServer_p.fftshift(outputDataStripe,
8787 8793 outputFlagsStripe,
8788 8794 (const Vector<Float>)inputDataStripe,
8789 8795 (const Vector<Bool>)inputFlagsStripe,
8790 8796 (const uInt)0, // In vectors axis 0 is the only dimension
8791 8797 (const Double)fftShift_p,
8792 8798 false); // A good data point has its flag set to false
8793 - return;
8794 8799 }
8795 8800 ​
8796 8801 // -----------------------------------------------------------------------
8797 8802 //
8798 8803 // -----------------------------------------------------------------------
8799 8804 template <class T> void MSTransformManager::interpol1D(Int inputSpw,
8800 8805 const Vector<T> &inputDataStripe,
8801 8806 const Vector<Bool> &inputFlagsStripe,
8802 8807 const Vector<Float> &,
8803 8808 Vector<T> &outputDataStripe,
8880 8885 continue;
8881 8886 ​
8882 8887 outputDataStripe[outIdx] = (outputDataStripe[outIdx] * outWeights[outIdx] +
8883 8888 intermDataStripe[mapIdx]) /
8884 8889 (1. + outWeights[outIdx]);
8885 8890 outWeights[outIdx] += 1;
8886 8891 outputFlagsStripe[outIdx] |= intermFlagsStripe[mapIdx];
8887 8892 }
8888 8893 }
8889 8894 ​
8890 -// -----------------------------------------------------------------------
8891 -//
8892 -// -----------------------------------------------------------------------
8895 +// ------------------------------------------------------------------------
8896 +// casacore::fftshift does not interpolate, it needs interpolation+fftshift
8897 +// ------------------------------------------------------------------------
8893 8898 template <class T> void MSTransformManager::interpol1Dfftshift(Int inputSpw,
8894 8899 const Vector<T> &inputDataStripe,
8895 8900 const Vector<Bool> &inputFlagsStripe,
8896 8901 const Vector<Float> &inputWeightsStripe,
8897 8902 Vector<T> &outputDataStripe,
8898 8903 Vector<Bool> &outputFlagsStripe)
8899 8904 {
8900 - Vector<T> regriddedDataStripe(inputDataStripe.shape(),T());
8901 - Vector<Bool> regriddedFlagsStripe(inputFlagsStripe.shape(),false);
8905 + Vector<T> regriddedDataStripe(outputDataStripe.shape(),T());
8906 + Vector<Bool> regriddedFlagsStripe(outputFlagsStripe.shape(),false);
8902 8907 ​
8903 - // This linear interpolation provides an uniform grid (pre-condition to apply fftshift)
8904 - interpol1D(inputSpw,inputDataStripe,inputFlagsStripe,inputWeightsStripe,regriddedDataStripe,regriddedFlagsStripe);
8908 + // This linear interpolation provides a uniform grid (pre-condition to apply fftshift)
8909 + interpol1D(inputSpw,inputDataStripe,inputFlagsStripe,inputWeightsStripe,regriddedDataStripe,regriddedFlagsStripe);
8905 8910 ​
8906 - // fftshift takes care of time
8907 - fftshift(inputSpw,regriddedDataStripe,regriddedFlagsStripe,inputWeightsStripe,outputDataStripe,outputFlagsStripe);
8908 -​
8909 - return;
8911 + // fftshift takes care of time
8912 + fftshift(inputSpw,regriddedDataStripe,regriddedFlagsStripe,inputWeightsStripe,outputDataStripe,outputFlagsStripe);
8910 8913 }
8911 8914 ​
8912 8915 // -----------------------------------------------------------------------
8913 8916 //
8914 8917 // -----------------------------------------------------------------------
8915 8918 template <class T> void MSTransformManager::averageRegrid(Int inputSpw,
8916 8919 const Vector<T> &inputDataStripe,
8917 8920 const Vector<Bool> &inputFlagsStripe,
8918 8921 const Vector<Float> &inputWeightsStripe,
8919 8922 Vector<T> &outputDataStripe,

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

Add shortcut