From 98edd33d68c2e0edab723e4877cd37b04e0a7fbb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@astron.nl>
Date: Mon, 17 Jan 2011 16:48:54 +0000
Subject: [PATCH] Bug 1491: directional filtering
---
.../include/AOFlagger/gui/mswindow.h | 7 ++-
.../AOFlagger/gui/timefrequencywidget.h | 1 +
.../AOFlagger/msio/timefrequencymetadata.h | 16 ++++++
.../AOFlagger/rfi/strategy/cutareaaction.h | 30 ++++++++++
.../rfi/strategy/fouriertransformaction.h | 11 +++-
.../include/AOFlagger/rfi/uvprojection.h | 11 +---
.../include/AOFlagger/util/ffttools.h | 1 +
.../src/gui/complexplaneplotwindow.cpp | 6 --
CEP/DP3/AOFlagger/src/gui/mswindow.cpp | 5 ++
CEP/DP3/AOFlagger/src/util/ffttools.cpp | 56 +++++++++++++++++++
10 files changed, 127 insertions(+), 17 deletions(-)
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h
index 598b4922604..5a317bbb21a 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h
+++ b/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;
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/timefrequencywidget.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/timefrequencywidget.h
index a5ab6886664..e53baabc98e 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/gui/timefrequencywidget.h
+++ b/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; }
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencymetadata.h b/CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencymetadata.h
index f3644835d07..9ce213ac20d 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/msio/timefrequencymetadata.h
+++ b/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();
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/cutareaaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/cutareaaction.h
index c5893321373..84874e31084 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/cutareaaction.h
+++ b/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);
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/fouriertransformaction.h b/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/fouriertransformaction.h
index 88b98d5d461..9b9ad7481e1 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/strategy/fouriertransformaction.h
+++ b/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
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/uvprojection.h b/CEP/DP3/AOFlagger/include/AOFlagger/rfi/uvprojection.h
index 6b27e112e4c..b88c9e7443b 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/rfi/uvprojection.h
+++ b/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);
- }
}
}
}
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/util/ffttools.h b/CEP/DP3/AOFlagger/include/AOFlagger/util/ffttools.h
index 740125ef626..251662fdf8c 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/util/ffttools.h
+++ b/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:
diff --git a/CEP/DP3/AOFlagger/src/gui/complexplaneplotwindow.cpp b/CEP/DP3/AOFlagger/src/gui/complexplaneplotwindow.cpp
index 558b9b5f29b..a99b72648b1 100644
--- a/CEP/DP3/AOFlagger/src/gui/complexplaneplotwindow.cpp
+++ b/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;
diff --git a/CEP/DP3/AOFlagger/src/gui/mswindow.cpp b/CEP/DP3/AOFlagger/src/gui/mswindow.cpp
index 0dfb2d5cf11..1c154d973be 100644
--- a/CEP/DP3/AOFlagger/src/gui/mswindow.cpp
+++ b/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'>"
diff --git a/CEP/DP3/AOFlagger/src/util/ffttools.cpp b/CEP/DP3/AOFlagger/src/util/ffttools.cpp
index 60c9e3d8ce1..3143dcdd24e 100644
--- a/CEP/DP3/AOFlagger/src/util/ffttools.cpp
+++ b/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();
--
GitLab