Bug 1491: directional filtering

CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h

@@ -115,6 +115,11 @@ class MSWindow : public Gtk::Window {
 		void onToggleFlags();
 		void onToggleMap();
 		void onToggleImage();
+		void onWinsorizedColorsPressed()
+		{
+			_timeFrequencyWidget.SetWindorizedColors(_winsorizedColorsButton->get_active());
+			_timeFrequencyWidget.Update();
+		}
 		void onCompress();
 		void onQuit() { hide(); }
 		void onActionFileOpen();
@@ -206,7 +211,7 @@ class MSWindow : public Gtk::Window {
 		Glib::RefPtr<Gtk::ToggleAction>
 			_originalFlagsButton, _altFlagsButton,
 			_originalImageButton, _backgroundImageButton, _diffImageButton,
-			_timeGraphButton;
+			_winsorizedColorsButton, _timeGraphButton;
 		Glib::RefPtr<Gtk::RadioAction>
 			_mapLogButton, _mapBWButton, _mapColorButton,
 			_gaussianTestSetsButton, _rayleighTestSetsButton;

CEP/DP3/AOFlagger/include/AOFlagger/gui/timefrequencywidget.h

@@ -43,6 +43,7 @@ class TimeFrequencyWidget : public Gtk::DrawingArea {
 		void SetShowOriginalFlagging(bool newValue) { _showOriginalFlagging = newValue; }
 		void SetShowAlternativeFlagging(bool newValue) { _showAlternativeFlagging = newValue; }
 		void SetColorMap(TFMap colorMap) { _colorMap = colorMap; }
+		void SetWindorizedColors(bool value) { _winsorizedStretch = value; }
 		void Update(); 
 		void AddAlternativeFlagging(Mask2DCPtr mask);
 		Image2DCPtr Image() { return _image; }

CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencymetadata.h

@@ -45,6 +45,22 @@ class TimeFrequencyMetaData
 			_uvw(0)
 		{
 		}
+		TimeFrequencyMetaData(const TimeFrequencyMetaData &source)
+			: _antenna1(0), _antenna2(0), _band(0), _field(0), _observationTimes(0), _uvw(0)
+		{
+			if(source._antenna1 != 0)
+				_antenna1 = new AntennaInfo(*source._antenna1);
+			if(source._antenna2 != 0)
+				_antenna2 = new AntennaInfo(*source._antenna2);
+			if(source._band != 0)
+				_band = new BandInfo(*source._band);
+			if(source._field != 0)
+				_field = new FieldInfo(*source._field);
+			if(source._observationTimes != 0)
+				_observationTimes = new std::vector<double>(*source._observationTimes);
+			if(source._uvw != 0)
+				_uvw = new std::vector<class UVW>(*source._uvw);
+		}
 		~TimeFrequencyMetaData()
 		{
 			ClearAntenna1();

CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/cutareaaction.h

@@ -21,6 +21,7 @@
 #define CUTAREAACTION_H
 
 #include <AOFlagger/msio/timefrequencydata.h>
+#include <AOFlagger/msio/timefrequencymetadata.h>
 
 #include "actionblock.h"
 #include "artifactset.h"
@@ -60,11 +61,14 @@ namespace rfiStrategy {
 				TimeFrequencyData oldOriginal(artifacts.OriginalData());
 				TimeFrequencyData oldRevised(artifacts.RevisedData());
 				TimeFrequencyData oldContaminated(artifacts.ContaminatedData());
+				TimeFrequencyMetaDataCPtr oldMetaData = artifacts.MetaData();
 
 				Cut(artifacts.OriginalData());
 				Cut(artifacts.RevisedData());
 				Cut(artifacts.ContaminatedData());
 
+				artifacts.SetMetaData(Cut(oldMetaData));
+
 				ActionBlock::Perform(artifacts, progress);
 				
 				PlaceBack(artifacts.OriginalData(), oldOriginal);
@@ -74,6 +78,7 @@ namespace rfiStrategy {
 				artifacts.SetOriginalData(oldOriginal);
 				artifacts.SetRevisedData(oldRevised);
 				artifacts.SetContaminatedData(oldContaminated);
+				artifacts.SetMetaData(oldMetaData);
 			}
 
 			virtual ActionType Type() const { return CutAreaActionType; }
@@ -84,6 +89,31 @@ namespace rfiStrategy {
 				size_t endChannel = data.ImageHeight() - _bottomChannels;
 				data.Trim(_startTimeSteps, _topChannels, endTime, endChannel);
 			}
+			TimeFrequencyMetaDataPtr Cut(TimeFrequencyMetaDataCPtr &metaData)
+			{
+				TimeFrequencyMetaData *newMetaData = new TimeFrequencyMetaData(*metaData);
+				if(newMetaData->HasUVW())
+				{
+					std::vector<UVW> newUVW;
+					newUVW.insert(newUVW.begin(), newMetaData->UVW().begin()+_startTimeSteps, newMetaData->UVW().end()-(_startTimeSteps+_endTimeSteps));
+					newMetaData->SetUVW(newUVW);
+				}
+				if(newMetaData->HasObservationTimes())
+				{
+					std::vector<double> times;
+					times.insert(times.begin(), newMetaData->ObservationTimes().begin() +_startTimeSteps, newMetaData->ObservationTimes().end()-(_startTimeSteps+_endTimeSteps));
+					newMetaData->SetObservationTimes(times);
+				}
+				return TimeFrequencyMetaDataPtr(newMetaData);
+				if(newMetaData->HasBand())
+				{
+					BandInfo band(newMetaData->Band());
+					band.channelCount -= _topChannels - _bottomChannels;
+					band.channels.erase(band.channels.end() - (_topChannels+_bottomChannels), band.channels.end());
+					band.channels.erase(band.channels.begin(), band.channels.begin() + _topChannels);
+					newMetaData->SetBand(band);
+				}
+			}
 			void PlaceBack(class TimeFrequencyData &cuttedData, class TimeFrequencyData &oldData)
 			{
 				oldData.CopyFrom(cuttedData, _startTimeSteps, _topChannels);

CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/fouriertransformaction.h

@@ -35,7 +35,7 @@ namespace rfiStrategy {
 	class FourierTransformAction : public Action
 	{
 		public:
-			FourierTransformAction() : Action(), _inverse(false)
+			FourierTransformAction() : Action(), _inverse(false), _dynamic(false), _sections(64)
 			{
 			}
 			virtual std::string Description()
@@ -60,12 +60,19 @@ namespace rfiStrategy {
 				Image2DPtr
 					real = Image2D::CreateCopy(data.GetImage(0)),
 					imaginary = Image2D::CreateCopy(data.GetImage(1));
-				FFTTools::CreateHorizontalFFTImage(*real, *imaginary, _inverse);
+				if(_dynamic)
+				{
+					FFTTools::CreateDynamicHorizontalFFTImage(real, imaginary, _sections, _inverse);
+				} else {
+					FFTTools::CreateHorizontalFFTImage(*real, *imaginary, _inverse);
+				}
 				data.SetImage(0, real);
 				data.SetImage(1, imaginary);
 			}
 
 			bool _inverse;
+			bool _dynamic;
+			unsigned _sections;
 	};
 
 } // namespace

CEP/DP3/AOFlagger/include/AOFlagger/rfi/uvprojection.h

@@ -152,15 +152,10 @@ class UVProjection
 					rangeEnd = inputWidth - rangeStart;
 				for(unsigned x=0;x<inputWidth;++x)
 				{
-					if(x > rangeStart && x < rangeEnd)
-					{
-						if(weights->Value(x, y) != 0.0)
-							destination->SetValue(x, y, destination->Value(x, y) / weights->Value(x, y));
-						else
-							AOLogger::Warn << "UV projection did not fill entire range\n";
-					} else {
+					if(x > rangeStart && x < rangeEnd && weights->Value(x, y) != 0.0)
+						destination->SetValue(x, y, destination->Value(x, y) / weights->Value(x, y));
+					else
 						destination->SetValue(x, y, 0.0);
-					}
 				}
 			}
 		}

CEP/DP3/AOFlagger/include/AOFlagger/util/ffttools.h

@@ -81,6 +81,7 @@ class FFTTools{
 		static void Multiply(Image2D &leftReal, Image2D &leftImaginary, const Image2D &rightReal, const Image2D &rightImaginary);
 		static void Sqrt(Image2D &image);
 		static void CreateHorizontalFFTImage(Image2D &real, Image2D &imaginary, bool inverse=false);
+		static void CreateDynamicHorizontalFFTImage(Image2DPtr real, Image2DPtr imaginary, unsigned sections, bool inverse=false);
 		static Image2DPtr AngularTransform(Image2DCPtr input);
 		static void FFT(SampleRowPtr realRow, SampleRowPtr imaginaryRow);
 	private:

CEP/DP3/AOFlagger/src/gui/complexplaneplotwindow.cpp

@@ -212,17 +212,12 @@ void ComplexPlanePlotWindow::onPlotPressed()
 				else
 					plot.StartLine("Fit (real)");
 				size_t middleY = (2*y + avgSize) / 2;
-				double timeStart = _observationTimes[x];
 				double deltaTime;
 				if(_observationTimes.size()>1)
 					deltaTime = _observationTimes[1] - _observationTimes[0];
 				else
 					deltaTime = 1.0;
-				double timeEnd = _observationTimes[x+length-1]+deltaTime;
-				long double frequency = _msWindow.TimeFrequencyMetaData()->Band().channels[middleY].frequencyHz;
 				Baseline baseline(_msWindow.TimeFrequencyMetaData()->Antenna1(), _msWindow.TimeFrequencyMetaData()->Antenna2());
-				long double delayRA = _msWindow.TimeFrequencyMetaData()->Field().delayDirectionRA;
-				long double delayDec = _msWindow.TimeFrequencyMetaData()->Field().delayDirectionDec;
 				long double fringeCount =
 					UVImager::GetFringeCount(x, x+length, middleY, _msWindow.TimeFrequencyMetaData());
 				long double fringeFrequency = fringeCount / length;
@@ -350,7 +345,6 @@ void ComplexPlanePlotWindow::setDetailsLabel()
 		deltaTime = _observationTimes[1] - _observationTimes[0];
 	else
 		deltaTime = 1.0;
-	double timeEnd = _observationTimes[x+length-1]+deltaTime;
 	long double frequency = metaData->Band().channels[middleY].frequencyHz;
 	Baseline baseline(metaData->Antenna1(), metaData->Antenna2());
 	long double delayRA = metaData->Field().delayDirectionRA;

CEP/DP3/AOFlagger/src/gui/mswindow.cpp

@@ -514,6 +514,9 @@ void MSWindow::createToolbar()
 	_timeGraphButton = Gtk::ToggleAction::create("TimeGraph", "Time graph");
 	_timeGraphButton->set_active(false); 
 	_actionGroup->add(_timeGraphButton, sigc::mem_fun(*this, &MSWindow::onTimeGraphButtonPressed) );
+	_winsorizedColorsButton = Gtk::ToggleAction::create("WinsorizedColors", "Winsorized colors");
+	_winsorizedColorsButton->set_active(true); 
+	_actionGroup->add(_winsorizedColorsButton, sigc::mem_fun(*this, &MSWindow::onWinsorizedColorsPressed) );
 
 	_actionGroup->add( Gtk::Action::create("PlotDist", "Plot _distribution"),
   sigc::mem_fun(*this, &MSWindow::onPlotDistPressed) );
@@ -681,6 +684,8 @@ void MSWindow::createToolbar()
     "      <menuitem action='MapLog'/>"
     "      <menuitem action='MapColor'/>"
     "      <separator/>"
+    "      <menuitem action='WinsorizedColors'/>"
+    "      <separator/>"
     "      <menuitem action='TimeGraph'/>"
 	  "    </menu>"
 	  "    <menu action='MenuPlot'>"

CEP/DP3/AOFlagger/src/util/ffttools.cpp

@@ -302,6 +302,62 @@ void FFTTools::CreateHorizontalFFTImage(Image2D &real, Image2D &imaginary, bool
 	fftw_free(in);
 }
 
+void FFTTools::CreateDynamicHorizontalFFTImage(Image2DPtr real, Image2DPtr imaginary, unsigned sections, bool inverse)
+{
+	const size_t width = real->Width();
+	if(real->Height() == 0 || width == 0) return;
+	SampleRowPtr
+		realRow = SampleRow::CreateFromRowSum(real, 0, real->Height()),
+		imaginaryRow = SampleRow::CreateFromRowSum(imaginary, 0, imaginary->Height());
+
+	Image2DPtr
+		destReal = Image2D::CreateEmptyImagePtr(real->Width(), real->Height()),
+		destImag = Image2D::CreateEmptyImagePtr(real->Width(), real->Height());
+	
+	unsigned long n_in = width;
+	fftw_complex *in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n_in);
+	fftw_complex *out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n_in);
+
+	int sign = -1;
+	if(inverse)
+		sign = 1;
+
+	for(unsigned sec=0;sec<sections;++sec)
+	{
+		const unsigned
+			secStart = width * sec / (sections + 1),
+			secEnd = width * (sec + 2) / (sections + 1);
+
+		for(unsigned x=secStart;x<secEnd;++x) {
+			in[x-secStart][0] = realRow->Value(x);
+			in[x-secStart][1] = imaginaryRow->Value(x);
+		}
+
+		fftw_plan plan = fftw_plan_dft_1d(secEnd - secStart, in, out, sign, FFTW_ESTIMATE);
+		fftw_execute(plan);
+		fftw_destroy_plan(plan);
+
+		size_t maxF = secEnd - secStart;
+		if(maxF > destReal->Height()) maxF = destReal->Height();
+		unsigned xEnd = width*(sec+1)/sections;
+		for(unsigned long x=width*sec/sections;x<xEnd;++x) {
+			for(unsigned long y=0;y<maxF;++y) {
+					destReal->SetValue(x, y, out[y][0]);
+					destImag->SetValue(x, y, out[y][1]);
+				}
+			for(unsigned long y=maxF;y<destReal->Height();++y)
+			{
+				destReal->SetValue(x, y, 0.0);
+				destImag->SetValue(x, y, 0.0);
+			}
+		}
+	}
+	fftw_free(out);
+	fftw_free(in);
+	real->SetValues(destReal);
+	imaginary->SetValues(destImag);
+}
+
 Image2DPtr FFTTools::AngularTransform(Image2DCPtr image)
 {
 	size_t minDim = image->Width() > image->Height() ? image->Height() : image->Width();