diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h b/CEP/DP3/AOFlagger/include/AOFlagger/gui/mswindow.h index 598b4922604cd74148b52ba2db50aade5c0d3f11..5a317bbb21a6d8d7356606cef7dd9003130aea29 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 a5ab688666441c0f810f76a3f99aa4053b87327f..e53baabc98eae0bf5d1e16978458ad4a18865186 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 f3644835d07f8ba3797852f8b896bbe37aab8119..9ce213ac20dc7017b1571c756870edb692e92686 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 c58933213738a6440b896ffadb376676ae85cd79..84874e310847d4948223e7d413a995addc281fd5 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 88b98d5d46139e7908b4a84295efaf011737f262..9b9ad7481e1a4cafa8889a563ff492debe8208db 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 6b27e112e4c2f58450056c056b6a8e87d0c4f02c..b88c9e7443bf941db4fd703353cc7aa618a4506e 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 740125ef62652a04dcf24724a71b206adef87756..251662fdf8cfb09bb7a8b3b302c0b752d7128b1d 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 558b9b5f29bf8bc495864f3dce3ee7ce1160f239..a99b72648b129512ee87c718635993327655a58f 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 0dfb2d5cf11032af5f04f70678f82fa3ce8833fd..1c154d973bec927228e1be379c28c1b01b5e8016 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 60c9e3d8ce1bfd4e3242249f67184359f152cd3c..3143dcdd24e37f07ef9069b4b67ca011e20a2592 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();