diff --git a/.gitattributes b/.gitattributes
index 3c6a54f7234036675b51fb28b33d22e33ed64dd2..7a6d4861f42afbfbfd90f2fe2953b13f36805a13 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -398,6 +398,7 @@ CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/frequencyselectionaction.h
CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/fringestopaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/imageraction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/iterationaction.h -text
+CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/plotaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/quickcalibrateaction.h -text
CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/rawappenderaction.h -text
@@ -642,6 +643,7 @@ CEP/DP3/AOFlagger/src/strategy/actions/foreachmsaction.cpp -text
CEP/DP3/AOFlagger/src/strategy/actions/frequencyselectionaction.cpp -text
CEP/DP3/AOFlagger/src/strategy/actions/fringestopaction.cpp -text
CEP/DP3/AOFlagger/src/strategy/actions/imageraction.cpp -text
+CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp -text
CEP/DP3/AOFlagger/src/strategy/actions/plotaction.cpp -text
CEP/DP3/AOFlagger/src/strategy/actions/slidingwindowfitaction.cpp -text
CEP/DP3/AOFlagger/src/strategy/actions/spatialcompositionaction.cpp -text
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/quality/defaultstatistics.h b/CEP/DP3/AOFlagger/include/AOFlagger/quality/defaultstatistics.h
index 76e8c424fd4d1ebda043c97a05121b8699f5e279..b001c011c97db5809f0d83080e080ba7be31c716 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/quality/defaultstatistics.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/quality/defaultstatistics.h
@@ -101,6 +101,25 @@ class DefaultStatistics : public Serializable
return *this;
}
+ DefaultStatistics ToSinglePolarization() const
+ {
+ if(_polarizationCount == 1)
+ return *this;
+
+ DefaultStatistics singlePol(1);
+ for(unsigned p=0;p<_polarizationCount;++p)
+ {
+ singlePol.rfiCount[0] += rfiCount[p];
+ singlePol.count[0] += count[p];
+ singlePol.sum[0] += sum[p];
+ singlePol.sumP2[0] += sumP2[p];
+ singlePol.dCount[0] += dCount[p];
+ singlePol.dSum[0] += dSum[p];
+ singlePol.dSumP2[0] += dSumP2[p];
+ }
+ return singlePol;
+ }
+
virtual void Serialize(std::ostream &stream) const
{
SerializeToUInt32(stream, _polarizationCount);
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/quality/statisticscollection.h b/CEP/DP3/AOFlagger/include/AOFlagger/quality/statisticscollection.h
index 54db7169a4c400d3d7cd51c82039729ab694b6de..7ee511e0c18dcf6b740797009af83899c70dee7c 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/quality/statisticscollection.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/quality/statisticscollection.h
@@ -127,6 +127,11 @@ class StatisticsCollection : public Serializable
loadBaseline<false>(qualityData);
}
+ void LoadTimeStatisticsOnly(QualityTablesFormatter &qualityData)
+ {
+ loadTime<false>(qualityData);
+ }
+
void Add(QualityTablesFormatter &qualityData)
{
loadTime<true>(qualityData);
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/action.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/action.h
index a1854242d107efe270293dfaa1f423b1601a1fb3..7dd743ed2c45ffe1d54e48e6b160e8006044232d 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/action.h
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/action.h
@@ -52,6 +52,7 @@ namespace rfiStrategy {
FringeStopActionType,
ImagerActionType,
IterationBlockType,
+ NormalizeVarianceActionType,
PlotActionType,
QuickCalibrateActionType,
RawAppenderActionType,
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h
new file mode 100644
index 0000000000000000000000000000000000000000..25c8f1bed93ea3a23897c95a806877dd1a9b59ee
--- /dev/null
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/strategy/actions/normalizevarianceaction.h
@@ -0,0 +1,59 @@
+/***************************************************************************
+ * Copyright (C) 2008 by A.R. Offringa *
+ * offringa@astro.rug.nl *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#ifndef RFISTRATEGYNORMALIZEVARIANCEACTION_H
+#define RFISTRATEGYNORMALIZEVARIANCEACTION_H
+
+#include <map>
+
+#include <AOFlagger/strategy/actions/action.h>
+
+#include <AOFlagger/strategy/control/artifactset.h>
+
+namespace rfiStrategy {
+
+ /**
+ @author A.R. Offringa <offringa@astro.rug.nl>
+ */
+ class NormalizeVarianceAction : public Action {
+ public:
+ NormalizeVarianceAction() :
+ _medianFilterSizeInS(60.0*15.0)
+ {
+ }
+ virtual ~NormalizeVarianceAction()
+ {
+ }
+ virtual std::string Description()
+ {
+ return "Normalize variance over time";
+ }
+ virtual void Perform(ArtifactSet &artifacts, ProgressListener &progress);
+
+ virtual ActionType Type() const { return NormalizeVarianceActionType; }
+
+ private:
+ double _medianFilterSizeInS;
+
+ void correctData(std::vector<Image2DPtr> &data, size_t timeStep, double stddev);
+ void correctDataUpTo(std::vector<Image2DPtr> &data, size_t &dataTimeIndex, double time, const std::vector<double> &observationTimes, double stddev);
+ };
+}
+
+#endif
diff --git a/CEP/DP3/AOFlagger/src/CMakeLists.txt b/CEP/DP3/AOFlagger/src/CMakeLists.txt
index c85b2fd76b1dc0684399cdca2d2ae64ec6887db6..693248aa5b6303c82f1f4c77c473408106a099da 100644
--- a/CEP/DP3/AOFlagger/src/CMakeLists.txt
+++ b/CEP/DP3/AOFlagger/src/CMakeLists.txt
@@ -116,6 +116,7 @@ set(STRATEGY_ACTION_FILES
strategy/actions/frequencyselectionaction.cpp
strategy/actions/fringestopaction.cpp
strategy/actions/imageraction.cpp
+ strategy/actions/normalizevarianceaction.cpp
strategy/actions/plotaction.cpp
strategy/actions/slidingwindowfitaction.cpp
strategy/actions/spatialcompositionaction.cpp
diff --git a/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp b/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..06d4b38c8b36fadb1ced4956e5b5f945a12c1e39
--- /dev/null
+++ b/CEP/DP3/AOFlagger/src/strategy/actions/normalizevarianceaction.cpp
@@ -0,0 +1,151 @@
+/***************************************************************************
+ * Copyright (C) 2008 by A.R. Offringa *
+ * offringa@astro.rug.nl *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+#include <AOFlagger/strategy/actions/normalizevarianceaction.h>
+
+#include <map>
+
+#include <AOFlagger/util/progresslistener.h>
+
+#include <AOFlagger/quality/defaultstatistics.h>
+#include <AOFlagger/quality/qualitytablesformatter.h>
+#include <AOFlagger/quality/statisticalvalue.h>
+#include <AOFlagger/quality/statisticscollection.h>
+
+#include <AOFlagger/strategy/algorithms/medianwindow.h>
+#include <AOFlagger/strategy/imagesets/msimageset.h>
+#include <AOFlagger/quality/statisticsderivator.h>
+
+namespace rfiStrategy {
+
+void NormalizeVarianceAction::Perform(ArtifactSet &artifacts, ProgressListener &progress)
+{
+ if(!artifacts.HasImageSet())
+ throw std::runtime_error("Normalize variance called without image set");
+ ImageSet *imageSet = artifacts.ImageSet();
+ MSImageSet *msImageSet = dynamic_cast<MSImageSet*>(imageSet);
+ if(msImageSet == 0)
+ throw std::runtime_error("Normalize variance actions needs measurement set");
+ std::string filename = msImageSet->Reader()->Set().Location();
+ QualityTablesFormatter qtables(filename);
+ StatisticsCollection statCollection(msImageSet->Reader()->Set().GetPolarizationCount());
+ statCollection.LoadTimeStatisticsOnly(qtables);
+ statCollection.IntegrateTimeToOneChannel();
+
+ TimeFrequencyData &original = artifacts.OriginalData();
+ const std::vector<double> &observationTimes = artifacts.MetaData()->ObservationTimes();
+ size_t
+ width = original.ImageWidth();
+
+ std::vector<Image2DPtr> data;
+ for(unsigned img=0;img<original.ImageCount();++img)
+ data.push_back(Image2D::CreateCopy(original.GetImage(img)));
+
+ // Calculate all stddevs
+ const std::map<double, DefaultStatistics> &statMap = statCollection.TimeStatistics();
+ std::map<double, double> stddevs;
+ for(std::map<double, DefaultStatistics>::const_iterator i = statMap.begin();
+ i != statMap.end(); ++i)
+ {
+ double stddev = StatisticsDerivator::GetStatisticAmplitude(
+ QualityTablesFormatter::DStandardDeviationStatistic,
+ i->second.ToSinglePolarization(), 0);
+ stddevs.insert(std::pair<double, double>(i->first, stddev));
+ }
+
+ // Add the first half of the window
+ const double halfWindowTime = _medianFilterSizeInS * 0.5;
+ MedianWindow<double> window;
+ std::map<double, double>::const_iterator windowRightSideIterator = stddevs.begin();
+ const double startTime = windowRightSideIterator->first;
+ do {
+ window.Add(windowRightSideIterator->second);
+ ++windowRightSideIterator;
+ } while(windowRightSideIterator->first - startTime < halfWindowTime);
+
+ // Add the second half, and start correcting the data
+ size_t dataTimeIndex = 0;
+ while(windowRightSideIterator != stddevs.end() &&
+ windowRightSideIterator->first - startTime < _medianFilterSizeInS)
+ {
+ correctDataUpTo(data, dataTimeIndex, windowRightSideIterator->first, observationTimes, window.Median());
+ window.Add(windowRightSideIterator->second);
+ ++windowRightSideIterator;
+ }
+
+ // Slide window until right side hits end
+ std::map<double, double>::const_iterator windowLeftSideIterator = stddevs.begin();
+ const double endTime = stddevs.rbegin()->first;
+ while(windowRightSideIterator != stddevs.end() && windowRightSideIterator->first < endTime)
+ {
+ correctDataUpTo(data, dataTimeIndex, windowRightSideIterator->first, observationTimes, window.Median());
+
+ window.Add(windowRightSideIterator->second);
+ window.Remove(windowLeftSideIterator->second);
+
+ ++windowRightSideIterator;
+ ++windowLeftSideIterator;
+ }
+
+ // Slide until window center hits end
+ while(windowRightSideIterator->first + halfWindowTime < endTime)
+ {
+ correctDataUpTo(data, dataTimeIndex, windowLeftSideIterator->first + _medianFilterSizeInS, observationTimes, window.Median());
+ window.Remove(windowLeftSideIterator->second);
+ ++windowLeftSideIterator;
+ }
+
+ while(dataTimeIndex < width)
+ {
+ correctData(data, dataTimeIndex, window.Median());
+ ++dataTimeIndex;
+ }
+
+ // Replace images
+ for(unsigned img=0;img<original.ImageCount();++img)
+ original.SetImage(img, data[img]);
+}
+
+void NormalizeVarianceAction::correctDataUpTo(std::vector<Image2DPtr> &data, size_t &dataTimeIndex, double rightSideTime, const std::vector<double> &observationTimes, double stddev)
+{
+ size_t width = (*data.begin())->Width();
+ double halfWindowWidth = _medianFilterSizeInS*0.5;
+ while(dataTimeIndex < width &&
+ observationTimes[dataTimeIndex] + halfWindowWidth < rightSideTime)
+ {
+ correctData(data, dataTimeIndex, stddev);
+ ++dataTimeIndex;
+ }
+}
+
+void NormalizeVarianceAction::correctData(std::vector<Image2DPtr> &data, size_t timeStep, double stddev)
+{
+ AOLogger::Debug << "Scaling timestep " << timeStep << " with " << stddev << '\n';
+ num_t oneOverStddev = 1.0 / stddev;
+
+ for(std::vector<Image2DPtr>::iterator i=data.begin();i!=data.end();++i)
+ {
+ Image2DPtr image = *i;
+
+ for(unsigned y=0;y<image->Height();++y)
+ image->SetValue(timeStep, y, image->Value(timeStep, y) * oneOverStddev);
+ }
+}
+
+} // end of namespace
diff --git a/CEP/DP3/AOFlagger/src/strategy/control/actionfactory.cpp b/CEP/DP3/AOFlagger/src/strategy/control/actionfactory.cpp
index d57c131eeafb8f746903bc1e5b1cfcbf606edd26..dc925153ee6ad0abcb04c21c9bd4e9454152dcd3 100644
--- a/CEP/DP3/AOFlagger/src/strategy/control/actionfactory.cpp
+++ b/CEP/DP3/AOFlagger/src/strategy/control/actionfactory.cpp
@@ -42,6 +42,7 @@
#include <AOFlagger/strategy/actions/fringestopaction.h>
#include <AOFlagger/strategy/actions/imageraction.h>
#include <AOFlagger/strategy/actions/iterationaction.h>
+#include <AOFlagger/strategy/actions/normalizevarianceaction.h>
#include <AOFlagger/strategy/actions/plotaction.h>
#include <AOFlagger/strategy/actions/quickcalibrateaction.h>
#include <AOFlagger/strategy/actions/rawappenderaction.h>
@@ -85,6 +86,7 @@ const std::vector<std::string> ActionFactory::GetActionList()
list.push_back("Fringe stopping recovery");
list.push_back("Image");
list.push_back("Iteration");
+ list.push_back("Normalize variance");
list.push_back("Phase adapter");
list.push_back("Plot");
list.push_back("Quickly calibrate");
@@ -149,6 +151,8 @@ Action *ActionFactory::CreateAction(const std::string &action)
return new ImagerAction();
else if(action == "Iteration")
return new IterationBlock();
+ else if(action == "Normalize variance")
+ return new NormalizeVarianceAction();
else if(action == "Phase adapter")
return new Adapter();
else if(action == "Plot")