diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a82b59ca88fca9a5a4bac3a2d3975be24f48fbd..8332a903bb8c5865480c0f14f360e35eebe2c6bd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -238,7 +238,6 @@ set(PLOT_FILES set(GUI_FILES rfigui/controllers/imagecomparisoncontroller.cpp rfigui/controllers/rfiguicontroller.cpp - rfigui/complexplaneplotwindow.cpp rfigui/gotowindow.cpp rfigui/imagepropertieswindow.cpp rfigui/opendialog.cpp @@ -308,10 +307,8 @@ set(ALGORITHMS_FILES algorithms/baselineselector.cpp algorithms/baselinetimeplaneimager.cpp algorithms/combinatorialthresholder.cpp - algorithms/fringestoppingfitter.cpp algorithms/fringetestcreater.cpp algorithms/highpassfilter.cpp - algorithms/localfitmethod.cpp algorithms/morphology.cpp algorithms/sinusfitter.cpp algorithms/siroperator.cpp diff --git a/algorithms/fringestoppingfitter.cpp b/algorithms/fringestoppingfitter.cpp deleted file mode 100644 index 2890da4c8847052b00c37d8230d4bdc8d66c236f..0000000000000000000000000000000000000000 --- a/algorithms/fringestoppingfitter.cpp +++ /dev/null @@ -1,424 +0,0 @@ -#include "../structures/antennainfo.h" - -#include "../imaging/uvimager.h" - -#include "../util/logger.h" -#include "../util/rng.h" - -#include "../algorithms/fringestoppingfitter.h" -#include "../algorithms/sinusfitter.h" - -namespace algorithms { - -FringeStoppingFitter::FringeStoppingFitter() - : _originalData(nullptr), - _fringesToConsider(1.0), - _minWindowSize(32), - _maxWindowSize(128), - _returnFittedValue(false), - _returnMeanValue(false), - _fringeFit(true), - _newPhaseCentreDec(M_PInl * 0.5), - _newPhaseCentreRA(0.0) {} - -FringeStoppingFitter::~FringeStoppingFitter() {} - -void FringeStoppingFitter::PerformFit(unsigned taskNumber) { - if (_fringeFit) { - PerformDynamicFrequencyFitOnOneChannel(taskNumber, _maxWindowSize); - } else { - const size_t x = taskNumber; - - if (_fitChannelsIndividually) { - for (size_t y = 0; y < _originalData->ImageHeight(); ++y) { - num_t r, i; - CalculateFitValue(*_originalData->GetRealPart(), - *_originalData->GetImaginaryPart(), x, y, 1, r, i); - - _realBackground->SetValue(x, y, r); - _imaginaryBackground->SetValue(x, y, i); - } - } else { - num_t r, i; - CalculateFitValue(*_originalData->GetRealPart(), - *_originalData->GetImaginaryPart(), x, 0, - _originalData->ImageHeight(), r, i); - - for (size_t y = 0; y < _originalData->ImageHeight(); ++y) { - _realBackground->SetValue(x, y, r); - _imaginaryBackground->SetValue(x, y, i); - } - } - } -} - -void FringeStoppingFitter::PerformStaticFrequencyFitOnOneChannel(unsigned y) { - if (_fitChannelsIndividually) { - for (size_t x = 0; x < _originalData->ImageWidth(); ++x) { - num_t r, i; - CalculateFitValue(*_originalData->GetRealPart(), - *_originalData->GetImaginaryPart(), x, y, 1, r, i); - - _realBackground->SetValue(x, y, r); - _imaginaryBackground->SetValue(x, y, i); - } - } else { - for (size_t x = 0; x < _originalData->ImageWidth(); ++x) { - num_t r, i; - CalculateFitValue(*_originalData->GetRealPart(), - *_originalData->GetImaginaryPart(), x, 0, - _originalData->ImageHeight(), r, i); - - _realBackground->SetValue(x, y, r); - _imaginaryBackground->SetValue(x, y, i); - } - } -} - -void FringeStoppingFitter::PerformFringeStop() { - Image2DCPtr real = _originalData->GetRealPart(), - imaginary = _originalData->GetImaginaryPart(); - - const Baseline baseline(*_antenna1Info, *_antenna2Info); - - for (size_t x = 0; x < real->Width(); ++x) { - for (size_t y = 0; y < real->Height(); ++y) { - const num_t r = real->Value(x, y); - num_t i = imaginary->Value(x, y); - const num_t fringeCount = UVImager::GetFringeCount(0, x, y, _metaData); - - num_t cosfreq = cosn(fringeCount * 2.0L * M_PIn), - sinfreq = sinn(fringeCount * 2.0L * M_PIn); - - const num_t newR = r * cosfreq - i * sinfreq; - i = r * sinfreq + i * cosfreq; - - _realBackground->SetValue(x, y, newR); - _imaginaryBackground->SetValue(x, y, i); - } - } -} - -num_t FringeStoppingFitter::CalculateMaskedAverage(const Image2D& image, - size_t x, size_t yFrom, - size_t yLength) { - num_t average = 0.0L; - size_t count = 0; - for (size_t i = yFrom; i < yFrom + yLength; ++i) { - if (!_originalMask->Value(x, i) && std::isfinite(image.Value(x, i))) { - average += image.Value(x, i); - ++count; - } - } - return average / (long double)count; -} - -num_t FringeStoppingFitter::CalculateUnmaskedAverage(const Image2D& image, - size_t x, size_t yFrom, - size_t yLength) { - num_t average = 0.0L; - size_t count = 0; - for (size_t i = yFrom; i < yFrom + yLength; ++i) { - if (std::isfinite(image.Value(x, i))) { - average += image.Value(x, i); - ++count; - } - } - return average / (num_t)count; -} - -void FringeStoppingFitter::CalculateFitValue(const Image2D& real, - const Image2D& imaginary, size_t x, - size_t yFrom, size_t yLength, - num_t& rValue, num_t& iValue) { - size_t windowWidth; - const size_t yMid = yFrom + yLength / 2; - const num_t estimatedFrequency = GetFringeFrequency(x, yMid); - windowWidth = (size_t)ceil(_fringesToConsider / estimatedFrequency); - if (windowWidth > _maxWindowSize) windowWidth = _maxWindowSize; - if (windowWidth < _minWindowSize) windowWidth = _minWindowSize; - - size_t xLeft, xRight; - xRight = x + windowWidth / 2; - if (x >= (windowWidth / 2)) { - xLeft = x - (windowWidth / 2); - } else { - xLeft = 0; - xRight += (windowWidth / 2) - x; - } - if (xRight > real.Width()) { - if (xLeft > xRight - real.Width()) - xLeft -= xRight - real.Width(); - else - xLeft = 0; - xRight = real.Width(); - } - - const num_t fringeFrequency = GetFringeFrequency(x, yMid); - num_t* dataT = new num_t[xRight - xLeft]; - num_t* dataR = new num_t[xRight - xLeft]; - num_t* dataI = new num_t[xRight - xLeft]; - - size_t index = 0; - for (size_t i = xLeft; i < xRight; ++i) { - dataR[index] = CalculateMaskedAverage(real, i, yFrom, yLength); - dataI[index] = CalculateMaskedAverage(imaginary, i, yFrom, yLength); - if (std::isfinite(dataR[index]) && std::isfinite(dataI[index])) { - dataT[index] = i; - ++index; - } - } - if (index == 0) { - for (size_t i = xLeft; i < xRight; ++i) { - dataR[index] = CalculateUnmaskedAverage(real, i, yFrom, yLength); - dataI[index] = CalculateUnmaskedAverage(imaginary, i, yFrom, yLength); - if (std::isfinite(dataR[index]) && std::isfinite(dataI[index])) { - dataT[index] = i; - ++index; - } - } - } - - SinusFitter fitter; - num_t phaseR, phaseI, amplitude; - - fitter.FindPhaseAndAmplitudeComplex(phaseR, amplitude, dataR, dataI, dataT, - index, fringeFrequency); - phaseI = fmodn(phaseR + M_PIn * 0.5L, 2.0L * M_PIn); - const num_t rfiValueR = - SinusFitter::Value(phaseR, amplitude, x, fringeFrequency, 0.0); - const num_t rfiValueI = - SinusFitter::Value(phaseI, amplitude, x, fringeFrequency, 0.0); - - if (_returnFittedValue) { - const num_t meanR = fitter.FindMean(phaseR, amplitude, dataR, dataT, index, - fringeFrequency); - const num_t meanI = fitter.FindMean(phaseI, amplitude, dataI, dataT, index, - fringeFrequency); - rValue = rfiValueR + meanR; - iValue = rfiValueI + meanI; - } else if (_returnMeanValue) { - rValue = fitter.FindMean(phaseR, amplitude, dataR, dataT, index, - fringeFrequency); - iValue = fitter.FindMean(phaseI, amplitude, dataI, dataT, index, - fringeFrequency); - } else { - const num_t observedValueR = - CalculateMaskedAverage(real, x, yFrom, yLength); - const num_t observedValueI = - CalculateMaskedAverage(imaginary, x, yFrom, yLength); - rValue = observedValueR - rfiValueR; - iValue = observedValueI - rfiValueI; - } - - delete[] dataR; - delete[] dataI; - delete[] dataT; -} - -num_t FringeStoppingFitter::GetFringeFrequency(size_t x, size_t y) { - double deltaTime; - if (_observationTimes->size() > 1) - deltaTime = (*_observationTimes)[1] - (*_observationTimes)[0]; - else - deltaTime = 1.0; - const num_t observeFreq = _bandInfo->channels[y].frequencyHz; - const Baseline baseline(*_antenna1Info, *_antenna2Info); - const num_t delayRA = _fieldInfo->delayDirectionRA; - const num_t delayDec = _fieldInfo->delayDirectionDec; - return deltaTime * UVImager::GetFringeStopFrequency(x, baseline, delayRA, - delayDec, observeFreq, - _metaData); -} - -void FringeStoppingFitter::GetRFIValue(num_t& r, num_t& i, int x, int y, - num_t rfiPhase, num_t rfiStrength) { - const numl_t earthRotation = UVImager::TimeToEarthLattitude(x, _metaData); - const Baseline baseline = _metaData->Baseline(); - const numl_t newWPos = UVImager::GetWPosition( - _newPhaseCentreDec, _newPhaseCentreRA, - _metaData->Band().channels[y].frequencyHz, earthRotation, - baseline.DeltaX(), baseline.DeltaY()); - const num_t rotations = - UVImager::GetFringeCount(0, x, y, _metaData) - newWPos; - - r = cosn(rotations * 2.0 * M_PIn + rfiPhase) * rfiStrength; - i = -sinn(rotations * 2.0 * M_PIn + rfiPhase) * rfiStrength; -} - -void FringeStoppingFitter::GetMeanValue(num_t& rMean, num_t& iMean, num_t phase, - num_t amplitude, const SampleRow& real, - const SampleRow& imaginary, - unsigned xStart, unsigned xEnd, - unsigned y) { - rMean = 0.0; - iMean = 0.0; - for (unsigned t = xStart; t < xEnd; ++t) { - num_t r, i; - GetRFIValue(r, i, t, y, phase, amplitude); - rMean += real.Value(t) - r; - iMean += imaginary.Value(t) - i; - } - rMean /= (num_t)(xEnd - xStart); - iMean /= (num_t)(xEnd - xStart); -} - -void FringeStoppingFitter::MinimizeRFIFitError(num_t& phase, num_t& amplitude, - const SampleRow& real, - const SampleRow& imaginary, - unsigned xStart, unsigned xEnd, - unsigned y) const throw() { - // calculate 1/N * \sum_x v(t) e^{2 i \pi \tau_g(t)}, where \tau_g(t) is the - // number of phase rotations because of the geometric delay as function of - // time t. - - num_t sumR = 0.0, sumI = 0.0; - size_t n = 0; - - const num_t dx = - _metaData->Antenna2().position.x - _metaData->Antenna1().position.x; - const num_t dy = - _metaData->Antenna2().position.y - _metaData->Antenna1().position.y; - - for (unsigned t = xStart; t < xEnd; ++t) { - const num_t vR = real.Value(t); - const num_t vI = imaginary.Value(t); - - if (std::isfinite(vR) && std::isfinite(vI)) { - const num_t tRotation = - UVImager::TimeToEarthLattitude(_metaData->ObservationTimes()[t]); - const num_t tauge = UVImager::GetFringeCount(0, t, y, _metaData); - const num_t taugeNew = UVImager::GetWPosition( - _newPhaseCentreDec, _newPhaseCentreRA, - _metaData->Band().channels[y].frequencyHz, tRotation, dx, dy); - const num_t phaseShift = tauge - taugeNew; - - sumR += vR * cosn(-2.0 * M_PIn * phaseShift); - sumR += vI * sinn(-2.0 * M_PIn * phaseShift); - - sumI += vR * sinn(-2.0 * M_PIn * phaseShift); - sumI -= vI * cosn(-2.0 * M_PIn * phaseShift); - ++n; - } - } - - sumR /= (num_t)n; - sumI /= (num_t)n; - - phase = SinusFitter::Phase(sumR, sumI); - amplitude = sqrtn(sumR * sumR + sumI * sumI); -} - -void FringeStoppingFitter::PerformDynamicFrequencyFitOnOneChannel(unsigned y) { - SampleRow real = - SampleRow::MakeFromRow(_originalData->GetRealPart().get(), y), - imaginary = SampleRow::MakeFromRow( - _originalData->GetImaginaryPart().get(), y); - PerformDynamicFrequencyFitOnOneRow(real, imaginary, y); -} - -void FringeStoppingFitter::PerformDynamicFrequencyFitOnOneRow( - const SampleRow& real, const SampleRow& imaginary, unsigned y) { - num_t phase, strength; - MinimizeRFIFitError(phase, strength, real, imaginary, 0, - _originalData->ImageWidth(), y); - Logger::Debug << "Amplitude found: " << strength << " phase found: " << phase - << '\n'; - for (size_t x = 0; x < _originalData->ImageWidth(); ++x) { - num_t rfiR, rfiI; - GetRFIValue(rfiR, rfiI, x, y, phase, strength); - if (_returnFittedValue) { - num_t rMean, iMean; - GetMeanValue(rMean, iMean, phase, strength, real, imaginary, 0, - _originalData->ImageWidth(), y); - _realBackground->SetValue(x, y, rfiR + rMean); - _imaginaryBackground->SetValue(x, y, rfiI + iMean); - } else { - _realBackground->SetValue(x, y, rfiR); - _imaginaryBackground->SetValue(x, y, rfiI); - } - } -} - -void FringeStoppingFitter::PerformDynamicFrequencyFitOnOneChannel( - unsigned y, unsigned windowSize) { - SampleRow real = SampleRow::MakeFromRowWithMissings( - _originalData->GetRealPart().get(), _originalMask.get(), y), - imaginary = SampleRow::MakeFromRowWithMissings( - _originalData->GetImaginaryPart().get(), _originalMask.get(), - y); - PerformDynamicFrequencyFitOnOneRow(real, imaginary, y, windowSize); -} - -void FringeStoppingFitter::PerformDynamicFrequencyFitOnOneRow( - const SampleRow& real, const SampleRow& imaginary, unsigned y, - unsigned windowSize) { - const unsigned halfWindowSize = windowSize / 2; - for (size_t x = 0; x < real.Size(); ++x) { - size_t windowStart, windowEnd; - if (x > halfWindowSize) - windowStart = x - halfWindowSize; - else - windowStart = 0; - if (x + halfWindowSize < real.Size()) - windowEnd = x + halfWindowSize; - else - windowEnd = real.Size(); - num_t windowPhase, windowStrength; - MinimizeRFIFitError(windowPhase, windowStrength, real, imaginary, - windowStart, windowEnd, y); - - num_t rfiR, rfiI; - GetRFIValue(rfiR, rfiI, x, y, windowPhase, windowStrength); - if (_returnFittedValue) { - num_t rMean, iMean; - GetMeanValue(rMean, iMean, windowPhase, windowStrength, real, imaginary, - windowStart, windowEnd, y); - _realBackground->SetValue(x, y, rfiR + rMean); - _imaginaryBackground->SetValue(x, y, rfiI + iMean); - } else { - _realBackground->SetValue(x, y, rfiR); - _imaginaryBackground->SetValue(x, y, rfiI); - } - } -} - -void FringeStoppingFitter::PerformDynamicFrequencyFit() { - for (size_t y = 0; y < _originalData->ImageHeight(); ++y) { - PerformDynamicFrequencyFitOnOneChannel(y); - } -} - -void FringeStoppingFitter::PerformDynamicFrequencyFit(unsigned yStart, - unsigned yEnd, - unsigned windowSize) { - SampleRow real = SampleRow::MakeFromRowSum(_originalData->GetRealPart().get(), - yStart, yEnd), - imaginary = SampleRow::MakeFromRowSum( - _originalData->GetImaginaryPart().get(), yStart, yEnd); - PerformDynamicFrequencyFitOnOneRow(real, imaginary, (yStart + yEnd) / 2, - windowSize); -} - -void FringeStoppingFitter::PerformDynamicFrequencyFit(unsigned yStart, - unsigned yEnd) { - SampleRow real = SampleRow::MakeFromRowSum(_originalData->GetRealPart().get(), - yStart, yEnd), - imaginary = SampleRow::MakeFromRowSum( - _originalData->GetImaginaryPart().get(), yStart, yEnd); - PerformDynamicFrequencyFitOnOneRow(real, imaginary, (yStart + yEnd) / 2); -} - -num_t FringeStoppingFitter::GetAmplitude(unsigned yStart, unsigned yEnd) { - const unsigned y = (yStart + yEnd) / 2; - SampleRow real = SampleRow::MakeFromRowSum(_originalData->GetRealPart().get(), - yStart, yEnd), - imaginary = SampleRow::MakeFromRowSum( - _originalData->GetImaginaryPart().get(), yStart, yEnd); - num_t phase, amplitude; - MinimizeRFIFitError(phase, amplitude, real, imaginary, 0, - _originalData->ImageWidth(), y); - return amplitude; -} - -} // namespace algorithms diff --git a/algorithms/fringestoppingfitter.h b/algorithms/fringestoppingfitter.h deleted file mode 100644 index d1f5207e60b39ff797bb3785dcf0308b562189ae..0000000000000000000000000000000000000000 --- a/algorithms/fringestoppingfitter.h +++ /dev/null @@ -1,129 +0,0 @@ -#ifndef FRINGESTOPPINGFITTER_H -#define FRINGESTOPPINGFITTER_H - -#include <vector> - -#include "../algorithms/surfacefitmethod.h" - -#include "../structures/samplerow.h" -#include "../structures/timefrequencymetadata.h" - -namespace algorithms { - -class FringeStoppingFitter final : public SurfaceFitMethod { - public: - FringeStoppingFitter(); - virtual ~FringeStoppingFitter(); - - void SetMetaData(TimeFrequencyMetaDataCPtr metaData) { - _metaData = metaData; - _fieldInfo = &metaData->Field(); - _bandInfo = &metaData->Band(); - _antenna1Info = &metaData->Antenna1(); - _antenna2Info = &metaData->Antenna2(); - _observationTimes = &metaData->ObservationTimes(); - } - void Initialize(const TimeFrequencyData& input) override { - _originalData = &input; - _realBackground = Image2D::CreateZeroImagePtr(_originalData->ImageWidth(), - _originalData->ImageHeight()); - _imaginaryBackground = Image2D::CreateZeroImagePtr( - _originalData->ImageWidth(), _originalData->ImageHeight()); - _originalMask = input.GetSingleMask(); - } - unsigned int TaskCount() const { - return _fringeFit ? _originalData->ImageHeight() - : _originalData->ImageWidth(); - } - void PerformFit(unsigned taskNumber) override; - void PerformStaticFrequencyFitOnOneChannel(unsigned y); - void PerformFringeStop(); - class TimeFrequencyData Background() const { - return TimeFrequencyData(aocommon::Polarization::StokesI, _realBackground, - _imaginaryBackground); - } - void SetFringesToConsider(long double fringesToConsider) { - _fringesToConsider = fringesToConsider; - } - void SetMinWindowSize(size_t minWindowSize) { - _minWindowSize = minWindowSize; - } - void SetMaxWindowSize(size_t maxWindowSize) { - _maxWindowSize = maxWindowSize; - } - void SetFitChannelsIndividually(bool fitChannelsIndividually) { - _fitChannelsIndividually = fitChannelsIndividually; - } - void SetReturnFittedValue(bool returnFittedValue) { - _returnFittedValue = returnFittedValue; - } - void SetReturnMeanValue(bool returnMeanValue) { - _returnMeanValue = returnMeanValue; - } - void PerformDynamicFrequencyFit(); - void PerformDynamicFrequencyFitOnOneChannel(unsigned y); - void PerformDynamicFrequencyFitOnOneChannel(unsigned y, unsigned windowSize); - void PerformDynamicFrequencyFit(unsigned yStart, unsigned yEnd); - void PerformDynamicFrequencyFit(unsigned yStart, unsigned yEnd, - unsigned windowSize); - num_t GetAmplitude(unsigned yStart, unsigned yEnd); - - numl_t NewPhaseCentreRA() const { return _newPhaseCentreRA; } - void SetNewPhaseCentreRA(long double newPhaseCentreRA) { - _newPhaseCentreRA = newPhaseCentreRA; - } - - numl_t NewPhaseCentreDec() const { return _newPhaseCentreDec; } - void SetNewPhaseCentreDec(long double newPhaseCentreDec) { - _newPhaseCentreDec = newPhaseCentreDec; - } - - private: - num_t CalculateFitValue(const Image2D& image, size_t y); - inline num_t CalculateMaskedAverage(const Image2D& image, size_t x, - size_t yFrom, size_t yLength); - inline num_t CalculateUnmaskedAverage(const Image2D& image, size_t x, - size_t yFrom, size_t yLength); - void CalculateFitValue(const Image2D& real, const Image2D& imaginary, - size_t x, size_t yFrom, size_t yLength, num_t& rValue, - num_t& iValue); - num_t GetFringeFrequency(size_t x, size_t y); - - void GetRFIValue(num_t& r, num_t& i, int x, int y, num_t rfiPhase, - num_t rfiStrength); - void GetMeanValue(num_t& rMean, num_t& iMean, num_t phase, num_t amplitude, - const SampleRow& real, const SampleRow& imaginary, - unsigned xStart, unsigned xEnd, unsigned y); - void MinimizeRFIFitError(num_t& phase, num_t& amplitude, - const SampleRow& real, const SampleRow& imaginary, - unsigned xStart, unsigned xEnd, unsigned y) const - throw(); - - void PerformDynamicFrequencyFitOnOneRow(const SampleRow& real, - const SampleRow& imaginary, - unsigned y); - void PerformDynamicFrequencyFitOnOneRow(const SampleRow& real, - const SampleRow& imaginary, - unsigned y, unsigned windowSize); - - Mask2DCPtr _originalMask; - const class TimeFrequencyData* _originalData; - - Image2DPtr _realBackground, _imaginaryBackground; - - TimeFrequencyMetaDataCPtr _metaData; - const class FieldInfo* _fieldInfo; - const class BandInfo* _bandInfo; - const class AntennaInfo *_antenna1Info, *_antenna2Info; - const std::vector<double>* _observationTimes; - num_t _fringesToConsider; - size_t _minWindowSize, _maxWindowSize; - bool _fitChannelsIndividually; - bool _returnFittedValue, _returnMeanValue; - bool _fringeFit; - numl_t _newPhaseCentreDec, _newPhaseCentreRA; -}; - -} // namespace algorithms - -#endif diff --git a/algorithms/localfitmethod.cpp b/algorithms/localfitmethod.cpp deleted file mode 100644 index 1cacc4a1835352d0279dec276fe4be062e12ebe6..0000000000000000000000000000000000000000 --- a/algorithms/localfitmethod.cpp +++ /dev/null @@ -1,305 +0,0 @@ -#include "localfitmethod.h" - -#include "../structures/image2d.h" -#include "../msio/pngfile.h" - -#include "../util/rng.h" - -#include "thresholdtools.h" - -#include <cmath> -#include <iostream> -#include <vector> - -namespace algorithms { - -LocalFitMethod::LocalFitMethod() : _weights(nullptr) {} - -LocalFitMethod::~LocalFitMethod() { ClearWeights(); } - -void LocalFitMethod::Initialize(const TimeFrequencyData& input) { - ClearWeights(); - _original = input.GetSingleImage(); - _background2D = - Image2D::CreateZeroImagePtr(_original->Width(), _original->Height()); - _background = TimeFrequencyData(input.ComplexRepresentation(), - input.Polarizations()[0], _background2D); - _mask = input.GetSingleMask(); - if (_hSquareSize * 2 > _original->Width()) - _hSquareSize = _original->Width() / 2; - if (_vSquareSize * 2 > _original->Height()) - _vSquareSize = _original->Height() / 2; - switch (_method) { - case None: - case Average: - case Median: - case Minimum: - break; - case GaussianWeightedAverage: - case FastGaussianWeightedAverage: - InitializeGaussianWeights(); - break; - } -} - -void LocalFitMethod::ClearWeights() { - if (_weights) { - for (unsigned y = 0; y < _vSquareSize * 2 + 1; ++y) delete[] _weights[y]; - delete[] _weights; - _weights = nullptr; - } -} - -void LocalFitMethod::InitializeGaussianWeights() { - _weights = new num_t*[_vSquareSize * 2 + 1]; - for (int y = 0; y < (int)(_vSquareSize * 2 + 1); ++y) { - _weights[y] = new num_t[_hSquareSize * 2 + 1]; - for (int x = 0; x < (int)(_hSquareSize * 2 + 1); ++x) { - _weights[y][x] = - RNG::EvaluateGaussian2D(x - (int)_hSquareSize, y - (int)_vSquareSize, - _hKernelSize, _vKernelSize); - } - } -} - -unsigned LocalFitMethod::TaskCount() { - if (_method == FastGaussianWeightedAverage) - return 1; - else - return _original->Height(); -} - -void LocalFitMethod::PerformFit(unsigned taskNumber) { - if (_mask == nullptr) throw std::runtime_error("Mask has not been set!"); - if (_method == FastGaussianWeightedAverage) { - CalculateWeightedAverageFast(); - } else { - const unsigned y = taskNumber; - for (unsigned x = 0; x < _original->Width(); ++x) - _background2D->SetValue(x, y, CalculateBackgroundValue(x, y)); - } -} - -long double LocalFitMethod::CalculateBackgroundValue(unsigned x, unsigned y) { - ThreadLocal local; - local.image = this; - local.currentX = x; - local.currentY = y; - - if (local.currentY >= _vSquareSize) - local.startY = local.currentY - _vSquareSize; - else - local.startY = 0; - local.endY = local.currentY + _vSquareSize; - if (local.endY >= _original->Height()) local.endY = _original->Height() - 1; - - if (local.currentX >= _hSquareSize) - local.startX = local.currentX - _hSquareSize; - else - local.startX = 0; - local.endX = local.currentX + _hSquareSize; - if (local.endX >= _original->Width()) local.endX = _original->Width() - 1; - local.emptyWindows = 0; - - switch (_method) { - case None: - case FastGaussianWeightedAverage: - return 0.0; - case Median: - return CalculateMedian(x, y, local); - case Minimum: - return CalculateMinimum(x, y, local); - case Average: - return CalculateAverage(x, y, local); - case GaussianWeightedAverage: - return CalculateWeightedAverage(x, y, local); - default: - throw std::runtime_error( - "The LocalFitMethod was not initialized before a fit was executed."); - } -} - -long double LocalFitMethod::CalculateAverage(unsigned x, unsigned y, - ThreadLocal& local) { - long double sum = 0.0; - unsigned long count = 0; - for (unsigned yi = local.startY; yi <= local.endY; ++yi) { - for (unsigned xi = local.startX; xi <= local.endX; ++xi) { - if (!_mask->Value(xi, yi) && std::isfinite(_original->Value(xi, yi))) { - sum += _original->Value(xi, yi); - ++count; - } - } - } - if (count != 0) - return sum / (long double)count; - else - return _original->Value(x, y); -} - -long double LocalFitMethod::CalculateMedian(unsigned x, unsigned y, - ThreadLocal& local) { - // unsigned maxSize = (local.endY-local.startY)*(local.endX-local.startX); - std::vector<long double> orderList; - unsigned long count = 0; - for (unsigned yi = local.startY; yi <= local.endY; ++yi) { - for (unsigned xi = local.startX; xi <= local.endX; ++xi) { - if (!_mask->Value(xi, yi) && std::isfinite(_original->Value(xi, yi))) { - orderList.push_back(_original->Value(xi, yi)); - ++count; - } - } - } - if (count == 0) { - local.emptyWindows++; - return _original->Value(x, y); - } else if (count % 2 == 1) { - std::nth_element(orderList.begin(), orderList.begin() + count / 2, - orderList.end()); - const long double mOdd = orderList[count / 2]; - return mOdd; - } else { - std::nth_element(orderList.begin(), orderList.begin() + count / 2, - orderList.end()); - const long double mOdd = orderList[count / 2]; - std::nth_element(orderList.begin(), orderList.begin() + (count / 2 - 1), - orderList.end()); - const long double mEven = orderList[count / 2 - 1]; - return (mOdd + mEven) * 0.5L; - } -} - -long double LocalFitMethod::CalculateMinimum(unsigned x, unsigned y, - ThreadLocal& local) { - long double minimum = 1e100; - unsigned long count = 0; - for (unsigned yi = local.startY; yi <= local.endY; ++yi) { - for (unsigned xi = local.startX; xi <= local.endX; ++xi) { - if (!_mask->Value(xi, yi) && std::isfinite(_original->Value(xi, yi)) && - _original->Value(xi, yi) < minimum) { - minimum = _original->Value(xi, yi); - ++count; - } - } - } - if (count == 0) - return _original->Value(x, y); - else - return minimum; -} - -long double LocalFitMethod::CalculateWeightedAverage(unsigned x, unsigned y, - ThreadLocal& local) { - long double sum = 0.0; - long double totalWeight = 0.0; - for (unsigned j = local.startY; j <= local.endY; ++j) { - for (unsigned i = local.startX; i <= local.endX; ++i) { - if (!_mask->Value(i, j) && std::isfinite(_original->Value(i, j))) { - const long double weight = - _weights[j - y + _vSquareSize][i - x + _hSquareSize]; - sum += _original->Value(i, j) * weight; - totalWeight += weight; - } - } - } - if (totalWeight != 0.0) { - return sum / totalWeight; - } else { - sum = 0.0; - totalWeight = 0.0; - for (unsigned j = local.startY; j <= local.endY; ++j) { - for (unsigned i = local.startX; i <= local.endX; ++i) { - if (std::isfinite(_original->Value(i, j))) { - const long double weight = - _weights[j - y + _vSquareSize][i - x + _hSquareSize]; - sum += _original->Value(i, j) * weight; - totalWeight += weight; - } - } - } - if (totalWeight != 0.0) { - return sum / totalWeight; - } else { - sum = 0.0; - totalWeight = 0.0; - for (unsigned j = local.startY; j <= local.endY; ++j) { - for (unsigned i = local.startX; i <= local.endX; ++i) { - const long double weight = - _weights[j - y + _vSquareSize][i - x + _hSquareSize]; - sum += _original->Value(i, j) * weight; - totalWeight += weight; - } - } - return sum / totalWeight; - } - } -} - -void LocalFitMethod::CalculateWeightedAverageFast() { - *_background2D = *_original; - ThresholdTools::SetFlaggedValuesToZero(_background2D.get(), _mask.get()); - PerformGaussianConvolution(_background2D); - const Image2DPtr flagWeights = CreateFlagWeightsMatrix(); - PerformGaussianConvolution(flagWeights); - ElementWiseDivide(_background2D, flagWeights); -} - -void LocalFitMethod::ElementWiseDivide(Image2DPtr leftHand, - Image2DCPtr rightHand) { - for (unsigned y = 0; y < leftHand->Height(); ++y) { - for (unsigned x = 0; x < leftHand->Width(); ++x) { - if (rightHand->Value(x, y) == 0.0) - leftHand->SetValue(x, y, 0.0); - else - leftHand->SetValue(x, y, - leftHand->Value(x, y) / rightHand->Value(x, y)); - } - } -} - -Image2DPtr LocalFitMethod::CreateFlagWeightsMatrix() { - Image2DPtr image = - Image2D::CreateUnsetImagePtr(_mask->Width(), _mask->Height()); - for (unsigned y = 0; y < image->Height(); ++y) { - for (unsigned x = 0; x < image->Width(); ++x) { - if (!_mask->Value(x, y) && std::isfinite(_original->Value(x, y))) - image->SetValue(x, y, 1.0); - else - image->SetValue(x, y, 0.0); - } - } - return image; -} - -void LocalFitMethod::PerformGaussianConvolution(Image2DPtr input) { - // Guassian convolution can be separated in two 1D convolution - // because of properties of the 2D Gaussian function. - const Image2DPtr temp = - Image2D::CreateZeroImagePtr(input->Width(), input->Height()); - for (int i = -_hSquareSize; i <= (int)_hSquareSize; ++i) { - const num_t gaus = _weights[_vSquareSize][i + _hSquareSize]; - for (unsigned y = 0; y < input->Height(); ++y) { - const unsigned xStart = i >= 0 ? 0 : -i; - const unsigned xEnd = i <= 0 ? input->Width() : input->Width() - i; - for (unsigned x = xStart; x < xEnd; ++x) { - if (std::isfinite(input->Value(x + i, y))) - temp->AddValue(x, y, input->Value(x + i, y) * gaus); - } - } - } - - input->SetAll(0.0); - for (int j = -_vSquareSize; j <= (int)_vSquareSize; ++j) { - const num_t gaus = _weights[j + _vSquareSize][_hSquareSize]; - const unsigned yStart = j >= 0 ? 0 : -j; - const unsigned yEnd = j <= 0 ? input->Height() : input->Height() - j; - for (unsigned y = yStart; y < yEnd; ++y) { - for (unsigned x = 0; x < input->Width(); ++x) { - if (std::isfinite(temp->Value(x, y + j))) - input->AddValue(x, y, temp->Value(x, y + j) * gaus); - } - } - } -} - -} // namespace algorithms diff --git a/algorithms/localfitmethod.h b/algorithms/localfitmethod.h deleted file mode 100644 index 6f5848dfefdabd778759e4e2c87ba77c3e9ea5b0..0000000000000000000000000000000000000000 --- a/algorithms/localfitmethod.h +++ /dev/null @@ -1,104 +0,0 @@ -#ifndef LocalFitMethod_H -#define LocalFitMethod_H - -#include <string> -#include <mutex> - -#include "../structures/image2d.h" -#include "../structures/mask2d.h" -#include "../structures/timefrequencydata.h" - -#include "surfacefitmethod.h" - -namespace algorithms { - -class LocalFitMethod final : public SurfaceFitMethod { - public: - enum Method { - None, - Average, - GaussianWeightedAverage, - FastGaussianWeightedAverage, - Median, - Minimum - }; - LocalFitMethod(); - ~LocalFitMethod(); - void SetToAverage(unsigned hSquareSize, unsigned vSquareSize) { - ClearWeights(); - _hSquareSize = hSquareSize; - _vSquareSize = vSquareSize; - _method = Average; - } - void SetToWeightedAverage(unsigned hSquareSize, unsigned vSquareSize, - long double hKernelSize, long double vKernelSize) { - ClearWeights(); - _hSquareSize = hSquareSize; - _vSquareSize = vSquareSize; - _method = FastGaussianWeightedAverage; - _hKernelSize = hKernelSize; - _vKernelSize = vKernelSize; - } - void SetToMedianFilter(unsigned hSquareSize, unsigned vSquareSize) { - ClearWeights(); - _hSquareSize = hSquareSize; - _vSquareSize = vSquareSize; - _method = Median; - } - void SetToMinimumFilter(unsigned hSquareSize, unsigned vSquareSize) { - ClearWeights(); - _hSquareSize = hSquareSize; - _vSquareSize = vSquareSize; - _method = Minimum; - } - void SetToNone() { - ClearWeights(); - _method = None; - } - void SetParameters(unsigned hSquareSize, unsigned vSquareSize, - enum Method method) { - ClearWeights(); - _hSquareSize = hSquareSize; - _vSquareSize = vSquareSize; - _method = method; - } - virtual void Initialize(const TimeFrequencyData& input) final override; - //[[ deprecated("Trying to make surfacemethod go away") ]] - unsigned TaskCount(); - virtual void PerformFit(unsigned taskNumber) final override; - TimeFrequencyData Background() const { return _background; } - - private: - struct ThreadLocal { - LocalFitMethod* image; - unsigned currentX, currentY; - unsigned startX, startY, endX, endY; - size_t emptyWindows; - }; - long double CalculateBackgroundValue(unsigned x, unsigned y); - long double FitBackground(unsigned x, unsigned y, ThreadLocal& local); - long double CalculateAverage(unsigned x, unsigned y, ThreadLocal& local); - long double CalculateMedian(unsigned x, unsigned y, ThreadLocal& local); - long double CalculateMinimum(unsigned x, unsigned y, ThreadLocal& local); - long double CalculateWeightedAverage(unsigned x, unsigned y, - ThreadLocal& local); - void ClearWeights(); - void InitializeGaussianWeights(); - void PerformGaussianConvolution(Image2DPtr input); - void CalculateWeightedAverageFast(); - Image2DPtr CreateFlagWeightsMatrix(); - void ElementWiseDivide(Image2DPtr leftHand, Image2DCPtr rightHand); - - Image2DCPtr _original; - TimeFrequencyData _background; - Image2DPtr _background2D; - Mask2DCPtr _mask; - unsigned _hSquareSize, _vSquareSize; - num_t** _weights; - long double _hKernelSize, _vKernelSize; - enum Method _method; -}; - -} // namespace algorithms - -#endif diff --git a/algorithms/surfacefitmethod.h b/algorithms/surfacefitmethod.h deleted file mode 100644 index f4ab22ebd33e0f1597d1714704a44d467ca49f10..0000000000000000000000000000000000000000 --- a/algorithms/surfacefitmethod.h +++ /dev/null @@ -1,20 +0,0 @@ -#ifndef IMAGE_H -#define IMAGE_H - -#include "../structures/image2d.h" -#include "../structures/timefrequencydata.h" - -#include "../util/ffttools.h" - -namespace algorithms { - -class SurfaceFitMethod { - public: - virtual void Initialize(const TimeFrequencyData& input) = 0; - virtual void PerformFit(unsigned taskNumber) = 0; - virtual ~SurfaceFitMethod() {} -}; - -} // namespace algorithms - -#endif diff --git a/algorithms/svdmitigater.h b/algorithms/svdmitigater.h index f03e6f3d61b37d5e884a14f9ec15b9f6b24fbdc2..3d6f23eaf7609b4d6ef29d929b05702ae7bf0e5b 100644 --- a/algorithms/svdmitigater.h +++ b/algorithms/svdmitigater.h @@ -4,8 +4,7 @@ #include <iostream> #include "../structures/image2d.h" - -#include "surfacefitmethod.h" +#include "../structures/timefrequencydata.h" // Needs to be included LAST #include "../util/f2c.h" @@ -14,12 +13,12 @@ class XYPlot; namespace algorithms { -class SVDMitigater final : public SurfaceFitMethod { +class SVDMitigater final { public: SVDMitigater(); ~SVDMitigater(); - void Initialize(const TimeFrequencyData& data) override { + void Initialize(const TimeFrequencyData& data) { Clear(); _data = data; _iteration = 0; @@ -28,7 +27,7 @@ class SVDMitigater final : public SurfaceFitMethod { _iteration++; RemoveSingularValues(_removeCount); } - void PerformFit(unsigned) override { PerformFit(); } + void PerformFit(unsigned) { PerformFit(); } void RemoveSingularValues(unsigned singularValueCount) { if (!IsDecomposed()) Decompose(); diff --git a/algorithms/testsetgenerator.cpp b/algorithms/testsetgenerator.cpp index 4651e92698316fffe1e8b8e842e425b4dc634f2a..6be91acb5d327f1469a602687c0b11c3e57f2185 100644 --- a/algorithms/testsetgenerator.cpp +++ b/algorithms/testsetgenerator.cpp @@ -4,6 +4,9 @@ #include <vector> #include <utility> +#include "combinatorialthresholder.h" +#include "highpassfilter.h" + #include "../structures/image2d.h" #include "../msio/pngfile.h" @@ -15,9 +18,6 @@ #include "../imaging/model.h" #include "../imaging/observatorium.h" -#include "combinatorialthresholder.h" -#include "localfitmethod.h" - namespace algorithms { void TestSetGenerator::AddSpectralLine(Image2D& data, Mask2D& rfi, @@ -521,22 +521,13 @@ void TestSetGenerator::SetModelData(Image2D& image, unsigned sources, void TestSetGenerator::SubtractBackground(Image2D& image) { const Mask2DPtr zero = Mask2D::CreateSetMaskPtr<false>(image.Width(), image.Height()); - LocalFitMethod fittedImage; - fittedImage.SetToWeightedAverage(20, 40, 7.5, 15.0); + HighPassFilter filter; + filter.SetHKernelSigmaSq(7.5 * 7.5); + filter.SetVKernelSigmaSq(15.0 * 15.0); + filter.SetHWindowSize(20); + filter.SetVWindowSize(40); const Image2DPtr imagePtr = Image2D::MakePtr(image); - TimeFrequencyData data(TimeFrequencyData::AmplitudePart, - aocommon::Polarization::StokesI, imagePtr); - data.SetGlobalMask(zero); - fittedImage.Initialize(data); - for (unsigned i = 0; i < fittedImage.TaskCount(); ++i) - fittedImage.PerformFit(i); - image = - Image2D::MakeFromDiff(image, *fittedImage.Background().GetSingleImage()); - for (unsigned y = 0; y < image.Height(); ++y) { - for (unsigned x = 0; x < image.Width(); ++x) { - image.AddValue(x, y, 1.0); - } - } + image = *filter.ApplyHighPass(imagePtr, zero); } Image2D TestSetGenerator::sampleRFIDistribution(unsigned width, unsigned height, diff --git a/rfigui/complexplaneplotwindow.cpp b/rfigui/complexplaneplotwindow.cpp deleted file mode 100644 index 18c249d4976321a50c333f3ed52809c5cd7faeee..0000000000000000000000000000000000000000 --- a/rfigui/complexplaneplotwindow.cpp +++ /dev/null @@ -1,384 +0,0 @@ -#include "complexplaneplotwindow.h" - -#include <set> -#include <sstream> - -#include <gtkmm/messagedialog.h> - -#include "../structures/antennainfo.h" -#include "../structures/date.h" -#include "../structures/image2d.h" - -#include "../util/rfiplots.h" - -#include "../algorithms/combinatorialthresholder.h" -#include "../algorithms/fringestoppingfitter.h" - -#include "../plot/plotmanager.h" - -#include "../imaging/uvimager.h" - -#include "rfiguiwindow.h" - -using algorithms::FringeStoppingFitter; - -ComplexPlanePlotWindow::ComplexPlanePlotWindow(RFIGuiWindow& rfiGuiWindow, - PlotManager& plotManager) - : _rfiGuiWindow(rfiGuiWindow), - _plotManager(plotManager), - _detailsFrame("Location details"), - _detailsLabel(), - _xPositionLabel("time start position:"), - _yPositionLabel("frequency start position:"), - _lengthLabel("time length:"), - _ySumLengthLabel("Frequency averaging size:"), - _xPositionScale(Gtk::ORIENTATION_HORIZONTAL), - _yPositionScale(Gtk::ORIENTATION_HORIZONTAL), - _lengthScale(Gtk::ORIENTATION_HORIZONTAL), - _ySumLengthScale(Gtk::ORIENTATION_HORIZONTAL), - _realVersusImaginaryButton("Real versus imaginary"), - _timeVersusRealButton("Time versus real"), - _allValuesButton("All values"), - _unmaskedValuesButton("Unmasked values"), - _maskedValuesButton("Masked values"), - _fittedValuesButton("Fitted values (constant freq)"), - _individualSampleFitButton("Channel fitted values (varying freq)"), - _fringeFitButton("Fringe fitted (varying freq)"), - _dynamicFringeFitButton("Dynamic fringe fitted (varying freq+amp)"), - _plotButton("Plot"), - _xMax(_rfiGuiWindow.GetOriginalData().ImageWidth()), - _yMax(_rfiGuiWindow.GetOriginalData().ImageHeight()) { - _detailsFrame.add(_detailsBox); - _detailsBox.show(); - - _detailsBox.pack_start(_detailsLabel); - _detailsLabel.set_line_wrap(true); - _detailsLabel.set_max_width_chars(40); - _detailsLabel.show(); - - _mainBox.pack_start(_detailsFrame); - _detailsFrame.show(); - - _mainBox.pack_start(_xPositionLabel); - _xPositionLabel.show(); - _xPositionScale.set_range(0.0, _rfiGuiWindow.GetOriginalData().ImageWidth()); - _xPositionScale.signal_value_changed().connect( - sigc::mem_fun(*this, &ComplexPlanePlotWindow::onTimeStartChanged)); - _mainBox.pack_start(_xPositionScale); - _xPositionScale.show(); - - _mainBox.pack_start(_yPositionLabel); - _yPositionLabel.show(); - _yPositionScale.set_range(0.0, _rfiGuiWindow.GetOriginalData().ImageHeight()); - _yPositionScale.signal_value_changed().connect( - sigc::mem_fun(*this, &ComplexPlanePlotWindow::onFreqChanged)); - _mainBox.pack_start(_yPositionScale); - _yPositionScale.show(); - - _mainBox.pack_start(_lengthLabel); - _lengthLabel.show(); - _lengthScale.set_range(1.0, _rfiGuiWindow.GetOriginalData().ImageWidth()); - _lengthScale.signal_value_changed().connect( - sigc::mem_fun(*this, &ComplexPlanePlotWindow::onTimeDurationChanged)); - _mainBox.pack_start(_lengthScale); - _lengthScale.show(); - - _mainBox.pack_start(_ySumLengthLabel); - _ySumLengthLabel.show(); - _ySumLengthScale.set_range(1.0, - _rfiGuiWindow.GetOriginalData().ImageHeight()); - _ySumLengthScale.signal_value_changed().connect( - sigc::mem_fun(*this, &ComplexPlanePlotWindow::onFreqSizeChanged)); - _mainBox.pack_start(_ySumLengthScale); - _ySumLengthScale.show(); - - _mainBox.pack_start(_realVersusImaginaryButton); - _realVersusImaginaryButton.show(); - Gtk::RadioButtonGroup group = _realVersusImaginaryButton.get_group(); - _timeVersusRealButton.set_group(group); - _realVersusImaginaryButton.set_active(true); - _mainBox.pack_start(_timeVersusRealButton); - _timeVersusRealButton.show(); - - _mainBox.pack_start(_allValuesButton); - _allValuesButton.set_active(true); - _allValuesButton.show(); - - _mainBox.pack_start(_maskedValuesButton); - _maskedValuesButton.show(); - - _mainBox.pack_start(_unmaskedValuesButton); - _unmaskedValuesButton.show(); - - _mainBox.pack_start(_fittedValuesButton); - _fittedValuesButton.set_active(true); - _fittedValuesButton.show(); - - _mainBox.pack_start(_individualSampleFitButton); - _individualSampleFitButton.show(); - - _mainBox.pack_start(_fringeFitButton); - _fringeFitButton.show(); - - _mainBox.pack_start(_dynamicFringeFitButton); - _dynamicFringeFitButton.show(); - - _plotButton.signal_clicked().connect( - sigc::mem_fun(*this, &ComplexPlanePlotWindow::onPlotPressed)); - _buttonBox.pack_start(_plotButton); - _plotButton.show(); - - _mainBox.pack_start(_buttonBox); - _buttonBox.show(); - - add(_mainBox); - _mainBox.show(); - - _observationTimes = _rfiGuiWindow.SelectedMetaData()->ObservationTimes(); - - setDetailsLabel(); -} - -ComplexPlanePlotWindow::~ComplexPlanePlotWindow() {} - -void ComplexPlanePlotWindow::onPlotPressed() { - if (_rfiGuiWindow.HasImage()) { - try { - XYPlot& plot = _plotManager.NewPlot2D("Complex plane"); - const size_t x = (size_t)_xPositionScale.get_value(); - const size_t y = (size_t)_yPositionScale.get_value(); - const size_t length = (size_t)_lengthScale.get_value(); - const size_t avgSize = (size_t)_ySumLengthScale.get_value(); - const bool realVersusImaginary = _realVersusImaginaryButton.get_active(); - const TimeFrequencyData& data = _rfiGuiWindow.GetActiveData(); - - if (_allValuesButton.get_active()) { - XYPointSet* pointSet; - if (realVersusImaginary) - pointSet = &plot.StartLine("Data", XYPointSet::DrawPoints); - else - pointSet = &plot.StartLine("Data (real)", XYPointSet::DrawPoints); - const Mask2DPtr mask = - Mask2D::CreateSetMaskPtr<false>(_rfiGuiWindow.AltMask()->Width(), - _rfiGuiWindow.AltMask()->Height()); - RFIPlots::MakeComplexPlanePlot(*pointSet, data, x, length, y, avgSize, - mask, realVersusImaginary, false); - - if (!realVersusImaginary) { - pointSet = - &plot.StartLine("Data (imaginary)", XYPointSet::DrawPoints); - RFIPlots::MakeComplexPlanePlot(*pointSet, data, x, length, y, avgSize, - mask, realVersusImaginary, true); - } - } - - if (_unmaskedValuesButton.get_active()) { - XYPointSet* pointSet = &plot.StartLine("Without RFI"); - RFIPlots::MakeComplexPlanePlot(*pointSet, data, x, length, y, avgSize, - _rfiGuiWindow.AltMask(), - realVersusImaginary, false); - if (!realVersusImaginary) { - pointSet = &plot.StartLine("Without RFI (I)"); - RFIPlots::MakeComplexPlanePlot(*pointSet, data, x, length, y, avgSize, - _rfiGuiWindow.AltMask(), - realVersusImaginary, true); - } - } - - if (_maskedValuesButton.get_active()) { - XYPointSet* pointSet = &plot.StartLine("Only RFI"); - const Mask2DPtr mask = Mask2D::MakePtr(*_rfiGuiWindow.AltMask()); - mask->Invert(); - RFIPlots::MakeComplexPlanePlot(*pointSet, data, x, length, y, avgSize, - mask, realVersusImaginary, false); - if (!realVersusImaginary) { - pointSet = &plot.StartLine("Only RFI (I)"); - RFIPlots::MakeComplexPlanePlot(*pointSet, data, x, length, y, avgSize, - mask, realVersusImaginary, true); - } - } - - if (_fittedValuesButton.get_active()) { - XYPointSet* pointSet; - if (realVersusImaginary) - pointSet = &plot.StartLine("Fit"); - else - pointSet = &plot.StartLine("Fit (real)"); - const size_t middleY = (2 * y + avgSize) / 2; - const Baseline baseline(_rfiGuiWindow.SelectedMetaData()->Antenna1(), - _rfiGuiWindow.SelectedMetaData()->Antenna2()); - const long double fringeCount = UVImager::GetFringeCount( - x, x + length, middleY, _rfiGuiWindow.SelectedMetaData()); - const long double fringeFrequency = fringeCount / length; - const Mask2DPtr mask = - Mask2D::CreateSetMaskPtr<false>(_rfiGuiWindow.AltMask()->Width(), - _rfiGuiWindow.AltMask()->Height()); - RFIPlots::MakeFittedComplexPlot(*pointSet, data, x, length, y, avgSize, - mask, fringeFrequency, - realVersusImaginary, false); - if (!realVersusImaginary) { - pointSet = &plot.StartLine("Fit (imaginary)"); - RFIPlots::MakeFittedComplexPlot(*pointSet, data, x, length, y, - avgSize, mask, fringeFrequency, - realVersusImaginary, true); - } - } - - if (_individualSampleFitButton.get_active()) { - FringeStoppingFitter fitter; - fitter.Initialize(data); - fitter.SetFitChannelsIndividually(true); - fitter.SetFringesToConsider(1.0L); - fitter.SetMaxWindowSize(256); - fitter.SetReturnFittedValue(true); - fitter.SetReturnMeanValue(false); - - fitter.SetMetaData(_rfiGuiWindow.SelectedMetaData()); - fitter.PerformStaticFrequencyFitOnOneChannel(y); - - XYPointSet* pointSet; - if (realVersusImaginary) - pointSet = &plot.StartLine("Fit"); - else - pointSet = &plot.StartLine("Fit (real)"); - const Mask2DPtr mask = - Mask2D::CreateSetMaskPtr<false>(_rfiGuiWindow.AltMask()->Width(), - _rfiGuiWindow.AltMask()->Height()); - RFIPlots::MakeComplexPlanePlot(*pointSet, fitter.Background(), x, - length, y, avgSize, mask, - realVersusImaginary, false); - - fitter.SetReturnFittedValue(false); - fitter.SetReturnMeanValue(true); - if (!realVersusImaginary) { - pointSet = &plot.StartLine("Fit (imaginary)"); - RFIPlots::MakeComplexPlanePlot(*pointSet, fitter.Background(), x, - length, y, avgSize, mask, - realVersusImaginary, true); - } - - fitter.PerformStaticFrequencyFitOnOneChannel(y); - - pointSet = &plot.StartLine("Center"); - RFIPlots::MakeComplexPlanePlot(*pointSet, fitter.Background(), x, - length, y, avgSize, mask, - realVersusImaginary, false); - - if (!realVersusImaginary) { - pointSet = &plot.StartLine("Center (I)"); - RFIPlots::MakeComplexPlanePlot(*pointSet, fitter.Background(), x, - length, y, avgSize, mask, - realVersusImaginary, true); - } - } - - if (_fringeFitButton.get_active() || - _dynamicFringeFitButton.get_active()) { - /*FringeStoppingFitter fitter; - Image2DPtr zero = Image2D::CreateZeroImagePtr(data.ImageWidth(), - data.ImageHeight()); Image2DPtr ones = - Image2D::CreateZeroImagePtr(data.ImageWidth(), data.ImageHeight()); - for(size_t yi=0;yi<ones->Height();++yi) - for(size_t xi=0;xi<ones->Width();++xi) - ones->SetValue(xi, yi, 1.0L); - TimeFrequencyData data(StokesIPolarisation, ones, zero); - fitter.Initialize(data); - fitter.SetFitChannelsIndividually(true); - - fitter.SetMetaData(_rfiGuiWindow.TimeFrequencyMetaData()); - fitter.PerformFringeStop(); - - plot.StartLine("Fringe rotation"); - Mask2DPtr mask = - Mask2D::CreateSetMaskPtr<false>(_rfiGuiWindow.AltMask()->Width(), - _rfiGuiWindow.AltMask()->Height()); RFIPlots::MakeComplexPlanePlot(plot, - fitter.Background(), x, length, y, avgSize, mask, realVersusImaginary, - false); - - if(!realVersusImaginary) - { - plot.StartLine("Fringe rotation (I)"); - RFIPlots::MakeComplexPlanePlot(plot, fitter.Background(), x, - length, y, avgSize, mask, realVersusImaginary, true); - }*/ - - FringeStoppingFitter fitter; - fitter.Initialize(data); - - fitter.SetMetaData(_rfiGuiWindow.SelectedMetaData()); - // fitter.PerformFringeStop(); - fitter.SetReturnFittedValue(true); - if (_dynamicFringeFitButton.get_active()) - fitter.PerformDynamicFrequencyFit(y, y + avgSize, 200); - else - fitter.PerformDynamicFrequencyFit(y, y + avgSize); - - XYPointSet* pointSet; - if (realVersusImaginary) - pointSet = &plot.StartLine("Fit"); - else - pointSet = &plot.StartLine("Fit (real)"); - const Mask2DPtr mask = - Mask2D::CreateSetMaskPtr<false>(_rfiGuiWindow.AltMask()->Width(), - _rfiGuiWindow.AltMask()->Height()); - RFIPlots::MakeComplexPlanePlot(*pointSet, fitter.Background(), x, - length, y, avgSize, mask, - realVersusImaginary, false); - - if (!realVersusImaginary) { - pointSet = &plot.StartLine("Fit (imaginary)"); - RFIPlots::MakeComplexPlanePlot(*pointSet, fitter.Background(), x, - length, y, avgSize, mask, - realVersusImaginary, true); - } - } - - _plotManager.Update(); - - } catch (std::exception& e) { - Gtk::MessageDialog dialog(*this, e.what(), false, Gtk::MESSAGE_ERROR); - dialog.run(); - } - } -} - -void ComplexPlanePlotWindow::setDetailsLabel() { - const size_t x = (size_t)_xPositionScale.get_value(); - const size_t y = (size_t)_yPositionScale.get_value(); - const size_t length = (size_t)_lengthScale.get_value(); - const size_t avgSize = (size_t)_ySumLengthScale.get_value(); - const size_t middleY = (2 * y + avgSize) / 2; - const TimeFrequencyMetaDataCPtr metaData = _rfiGuiWindow.SelectedMetaData(); - - const double timeStart = _observationTimes[x]; - double deltaTime; - if (_observationTimes.size() > 1) - deltaTime = _observationTimes[1] - _observationTimes[0]; - else - deltaTime = 1.0; - const long double frequency = metaData->Band().channels[middleY].frequencyHz; - const Baseline baseline(metaData->Antenna1(), metaData->Antenna2()); - const long double delayRA = metaData->Field().delayDirectionRA; - const long double delayDec = metaData->Field().delayDirectionDec; - const long double intFringeFreq = - UVImager::GetFringeCount(x, x + length, y, metaData); - const long double midFringeFreq = UVImager::GetFringeStopFrequency( - (x * 2 + length) / 2, baseline, delayRA, delayDec, frequency, metaData); - - std::stringstream s; - s << "Start time: " << Date::AipsMJDToString(timeStart) << std::endl - << "Frequency: " << frequency / 1000000.0L << "Mhz" << std::endl - << "Baseline: " << baseline.Distance() << "m" << std::endl - << "Delay direction: " << delayRA << "RA, " << delayDec << "dec." - << std::endl - << "(=" << RightAscension::ToString(delayRA) << " RA, " - << Declination::ToString(delayDec) << " dec.)" << std::endl - << "Mid fringe stopping freq: " << midFringeFreq << "(Hz)" << std::endl - << "Fringe count: " << intFringeFreq << std::endl - << "Fringe length: " << 1.0L / intFringeFreq << "(s)" << std::endl - << "Time step: " << deltaTime << "(s)" << std::endl - << "Samples/fringe: " << (1.0L / (deltaTime * intFringeFreq)) << std::endl - << "Fringes in domain: " << intFringeFreq << std::endl; - - _detailsLabel.set_text(s.str()); -} diff --git a/rfigui/complexplaneplotwindow.h b/rfigui/complexplaneplotwindow.h deleted file mode 100644 index 735f70a9ebd16e046f4fdf77297b22db8b504aba..0000000000000000000000000000000000000000 --- a/rfigui/complexplaneplotwindow.h +++ /dev/null @@ -1,76 +0,0 @@ -#ifndef COMPLEXPLANEPLOTWINDOW_H -#define COMPLEXPLANEPLOTWINDOW_H - -#include <vector> - -#include <gtkmm/box.h> -#include <gtkmm/button.h> -#include <gtkmm/buttonbox.h> -#include <gtkmm/checkbutton.h> -#include <gtkmm/frame.h> -#include <gtkmm/label.h> -#include <gtkmm/radiobutton.h> -#include <gtkmm/scale.h> -#include <gtkmm/window.h> - -class ComplexPlanePlotWindow : public Gtk::Window { - public: - ComplexPlanePlotWindow(class RFIGuiWindow& rfiGuiWindow, - class PlotManager& plotManager); - ~ComplexPlanePlotWindow(); - - private: - size_t XStart() const throw() { return (size_t)_xPositionScale.get_value(); } - size_t XLength() const throw() { return (size_t)_lengthScale.get_value(); } - - size_t YStart() const throw() { return (size_t)_yPositionScale.get_value(); } - size_t YLength() const throw() { - return (size_t)_ySumLengthScale.get_value(); - } - - void onPlotPressed(); - void onTimeStartChanged() { - if (XStart() + XLength() > _xMax) - _lengthScale.set_value(_xMax - XStart()); - else - setDetailsLabel(); - } - void onTimeDurationChanged() { - if (XStart() + XLength() > _xMax) - _xPositionScale.set_value(_xMax - XLength()); - else - setDetailsLabel(); - } - void onFreqChanged() { - if (YStart() + YLength() > _yMax) - _ySumLengthScale.set_value(_yMax - YStart()); - else - setDetailsLabel(); - } - void onFreqSizeChanged() { - if (YStart() + YLength() > _yMax) - _yPositionScale.set_value(_yMax - YLength()); - else - setDetailsLabel(); - } - void setDetailsLabel(); - - class RFIGuiWindow& _rfiGuiWindow; - class PlotManager& _plotManager; - Gtk::Frame _detailsFrame; - Gtk::VBox _mainBox, _detailsBox; - Gtk::Label _detailsLabel, _xPositionLabel, _yPositionLabel, _lengthLabel, - _ySumLengthLabel; - Gtk::Scale _xPositionScale, _yPositionScale, _lengthScale, _ySumLengthScale; - Gtk::RadioButton _realVersusImaginaryButton, _timeVersusRealButton; - Gtk::CheckButton _allValuesButton, _unmaskedValuesButton, _maskedValuesButton, - _fittedValuesButton, _individualSampleFitButton, _fringeFitButton, - _dynamicFringeFitButton; - - Gtk::Button _plotButton; - Gtk::ButtonBox _buttonBox; - std::vector<double> _observationTimes; - size_t _xMax, _yMax; -}; - -#endif diff --git a/rfigui/controllers/rfiguicontroller.cpp b/rfigui/controllers/rfiguicontroller.cpp index 12c442d70cd1ce01a90b753d76980f7d7579aba5..c428cbf53ea99cc133dc0e95ef2759dda4e3b904 100644 --- a/rfigui/controllers/rfiguicontroller.cpp +++ b/rfigui/controllers/rfiguicontroller.cpp @@ -1,6 +1,5 @@ #include "rfiguicontroller.h" -#include "../../algorithms/fringestoppingfitter.h" #include "../../algorithms/svdmitigater.h" #include "../../algorithms/testsetgenerator.h" diff --git a/rfigui/rfiguimenu.cpp b/rfigui/rfiguimenu.cpp index 1bfbc94d1d7d4fcd165090defa13fd4b177ede83..0661b0fad50f48f97590fc2eb60ac3dc86ed9b73 100644 --- a/rfigui/rfiguimenu.cpp +++ b/rfigui/rfiguimenu.cpp @@ -184,8 +184,6 @@ void RFIGuiMenu::makePlotMenu() { "Plot _distribution"); addItem(_menuPlot, _miPlotLogLogDistribution, OnPlotLogLogDistPressed, "Plot _log-log dist"); - addItem(_menuPlot, _miPlotComplexPlane, OnPlotComplexPlanePressed, - "Plot _complex plane"); addItem(_menuPlot, _miPlotSingularValues, OnPlotSingularValuesPressed, "Plot _singular values"); } diff --git a/rfigui/rfiguimenu.h b/rfigui/rfiguimenu.h index 529ecbd9aa7d200eda0a46fb9fc13bd3ee529d46..d212cc1f7e199558ef4ba347abf0819c359295f7 100644 --- a/rfigui/rfiguimenu.h +++ b/rfigui/rfiguimenu.h @@ -49,7 +49,6 @@ class RFIGuiMenu { // Plot sigc::signal<void> OnPlotDistPressed; sigc::signal<void> OnPlotLogLogDistPressed; - sigc::signal<void> OnPlotComplexPlanePressed; sigc::signal<void> OnPlotMeanSpectrumPressed; sigc::signal<void> OnPlotSumSpectrumPressed; sigc::signal<void> OnPlotPowerSpectrumPressed; @@ -338,7 +337,7 @@ class RFIGuiMenu { Gtk::MenuItem _miPlotTime, _miPlotFrequency; Gtk::Menu _menuPlotTime, _menuPlotFrequency; Gtk::MenuItem _miPlotDistribution, _miPlotLogLogDistribution, - _miPlotComplexPlane, _miPlotMeanSpectrum; + _miPlotMeanSpectrum; Gtk::MenuItem _miPlotSumSpectrum, _miPlotPowerSpectrum, _miPlotFrequencyScatter, _miPlotTimeMean; Gtk::MenuItem _miPlotTimeScatter, _miPlotSingularValues; diff --git a/rfigui/rfiguiwindow.cpp b/rfigui/rfiguiwindow.cpp index 2fbcebd84f2d5e969161ad9979e51b4d75dd080a..1356db479071fa2046a6f35c22838528ceabe7c9 100644 --- a/rfigui/rfiguiwindow.cpp +++ b/rfigui/rfiguiwindow.cpp @@ -29,7 +29,6 @@ #include "controllers/imagecomparisoncontroller.h" #include "controllers/rfiguicontroller.h" -#include "complexplaneplotwindow.h" #include "gotowindow.h" #include "histogramwindow.h" #include "imagepropertieswindow.h" @@ -126,8 +125,6 @@ RFIGuiWindow::RFIGuiWindow(RFIGuiController* controller, // Plot _menu->OnPlotDistPressed.connect([&]() { onPlotDistPressed(); }); _menu->OnPlotLogLogDistPressed.connect([&]() { onPlotLogLogDistPressed(); }); - _menu->OnPlotComplexPlanePressed.connect( - [&]() { onPlotComplexPlanePressed(); }); _menu->OnPlotMeanSpectrumPressed.connect( [&]() { onPlotMeanSpectrumPressed(); }); _menu->OnPlotSumSpectrumPressed.connect( @@ -781,14 +778,6 @@ void RFIGuiWindow::onPlotDistPressed() { _controller->PlotDist(); } void RFIGuiWindow::onPlotLogLogDistPressed() { _controller->PlotLogLogDist(); } -void RFIGuiWindow::onPlotComplexPlanePressed() { - if (HasImage()) { - _plotComplexPlaneWindow.reset( - new ComplexPlanePlotWindow(*this, _controller->PlotManager())); - _plotComplexPlaneWindow->show(); - } -} - void RFIGuiWindow::onPlotPowerSpectrumPressed() { _controller->PlotPowerSpectrum(); } diff --git a/rfigui/rfiguiwindow.h b/rfigui/rfiguiwindow.h index 5e1d5ee538335e52806fe869f8b634c5a76ab4ae..6a7b1fd8cd5ec3862eeb09d5f488e3c2a950155b 100644 --- a/rfigui/rfiguiwindow.h +++ b/rfigui/rfiguiwindow.h @@ -192,7 +192,6 @@ class RFIGuiWindow : public Gtk::Window { void onShowStats(); void onPlotDistPressed(); void onPlotLogLogDistPressed(); - void onPlotComplexPlanePressed(); void onPlotMeanSpectrumPressed(); void onPlotSumSpectrumPressed(); void onPlotPowerSpectrumPressed(); @@ -245,8 +244,7 @@ class RFIGuiWindow : public Gtk::Window { std::unique_ptr<class HistogramWindow> _histogramWindow; std::unique_ptr<class PlotWindow> _plotWindow; - std::unique_ptr<Gtk::Window> _gotoWindow, _plotComplexPlaneWindow, - _imagePropertiesWindow; + std::unique_ptr<Gtk::Window> _gotoWindow, _imagePropertiesWindow; std::unique_ptr<class ProgressWindow> _progressWindow; std::unique_ptr<class RFIGuiMenu> _menu; diff --git a/test/algorithms/highpassfiltertest.cpp b/test/algorithms/highpassfiltertest.cpp index 657e6ca42d80c013e0fbf827069f8a884ce50f9b..756c235d9547794d74f4003f3ab316805c9b8f6d 100644 --- a/test/algorithms/highpassfiltertest.cpp +++ b/test/algorithms/highpassfiltertest.cpp @@ -3,85 +3,14 @@ #include "../../structures/mask2d.h" #include "../../structures/image2d.h" -#include "../../algorithms/localfitmethod.h" #include "../../algorithms/highpassfilter.h" #include <boost/test/unit_test.hpp> -using aocommon::Polarization; - using algorithms::HighPassFilter; -using algorithms::LocalFitMethod; BOOST_AUTO_TEST_SUITE(high_pass_filter, *boost::unit_test::label("algorithms")) -BOOST_AUTO_TEST_CASE(filter) { - const size_t width = 99, height = 99; - Image2DPtr testImage = Image2D::CreateZeroImagePtr(width, height); - testImage->SetValue(10, 10, 1.0); - testImage->SetValue(15, 15, 2.0); - testImage->SetValue(20, 20, 0.5); - - // Fitting - LocalFitMethod fitMethod; - TimeFrequencyData data(TimeFrequencyData::AmplitudePart, - Polarization::StokesI, - Image2DCPtr(new Image2D(*testImage))); - fitMethod.SetToWeightedAverage(10, 20, 2.5, 5.0); - fitMethod.Initialize(data); - for (size_t i = 0; i < fitMethod.TaskCount(); ++i) fitMethod.PerformFit(i); - Image2DPtr fitResult(new Image2D(Image2D::MakeFromDiff( - *testImage, *fitMethod.Background().GetSingleImage()))); - - // High-pass filter - HighPassFilter filter; - Image2DPtr filterResult(new Image2D(*testImage)); - filter.SetHWindowSize(21); - filter.SetVWindowSize(41); - filter.SetHKernelSigmaSq(2.5); - filter.SetVKernelSigmaSq(5.0); - filterResult = filter.ApplyHighPass( - filterResult, Mask2D::CreateSetMaskPtr<false>(width, height)); - - ImageAsserter::AssertEqual(filterResult, fitResult); -} - -BOOST_AUTO_TEST_CASE(small_image_filter) { - const size_t width = 8, height = 8; - Image2DPtr testImage = Image2D::CreateZeroImagePtr(width, height); - testImage->SetValue(1, 1, 1.0); - - // Fitting - LocalFitMethod fitMethod; - TimeFrequencyData data(TimeFrequencyData::AmplitudePart, - Polarization::StokesI, - Image2DPtr(new Image2D(*testImage))); - fitMethod.SetToWeightedAverage(10, 20, 2.5, 5.0); - fitMethod.Initialize(data); - for (size_t i = 0; i < fitMethod.TaskCount(); ++i) fitMethod.PerformFit(i); - Image2DCPtr fitResult(new Image2D(Image2D::MakeFromDiff( - *testImage, *fitMethod.Background().GetSingleImage()))); - - // High-pass filter - HighPassFilter filter; - Image2DPtr filterResult(new Image2D(*testImage)); - filter.SetHWindowSize(21); - filter.SetVWindowSize(41); - filter.SetHKernelSigmaSq(2.5); - filter.SetVKernelSigmaSq(5.0); - filterResult = filter.ApplyHighPass( - filterResult, Mask2D::CreateSetMaskPtr<false>(width, height)); - - BOOST_CHECK(true); // avoid warnings by boost - - // This test will fail, but the high-pass filter is actually slightly better - // than the older "fitter" -- it will keep the kernel as large as possible, - // while the sliding window fit can be one off. The test is still good to - // guard for out of bounds errors. - // ImageAsserter::AssertEqual(filterResult, fitResult, "Convolution with - // kernel that is larger than the image"); -} - BOOST_AUTO_TEST_CASE(filter_with_mask) { const size_t width = 8, height = 8; Image2DPtr image = Image2D::CreateZeroImagePtr(width, height); diff --git a/test/experiments/highpassfilterexperiment.cpp b/test/experiments/highpassfilterexperiment.cpp index 61cebeeebacfd139c4edf30f3cbc7ec08c1f6c08..79f273899799f2ae02c18d9286f89476f7c1e2a5 100644 --- a/test/experiments/highpassfilterexperiment.cpp +++ b/test/experiments/highpassfilterexperiment.cpp @@ -1,17 +1,13 @@ #include <boost/test/unit_test.hpp> #include "../../algorithms/highpassfilter.h" -#include "../../algorithms/localfitmethod.h" #include "../../util/rng.h" #include "../../util/stopwatch.h" #include <iostream> -using aocommon::Polarization; - using algorithms::HighPassFilter; -using algorithms::LocalFitMethod; BOOST_AUTO_TEST_SUITE(high_pass_filter_experiment, *boost::unit_test::label("experiment") * @@ -45,21 +41,6 @@ static void InitializeFlagged(Image2DPtr& image, Mask2DPtr& mask) { } } -BOOST_AUTO_TEST_CASE(time_fitting) { - Image2DPtr image; - Mask2DPtr mask; - Initialize(image, mask); - - LocalFitMethod fitMethod; - TimeFrequencyData data(TimeFrequencyData::AmplitudePart, - Polarization::StokesI, image); - fitMethod.SetToWeightedAverage(10, 20, 2.5, 5.0); - fitMethod.Initialize(data); - Stopwatch watch(true); - for (size_t i = 0; i < fitMethod.TaskCount(); ++i) fitMethod.PerformFit(i); - std::cout << " time token: " << watch.ToString() << ' '; -} - BOOST_AUTO_TEST_CASE(time_high_pass_filter) { Image2DPtr image; Mask2DPtr mask; @@ -75,21 +56,6 @@ BOOST_AUTO_TEST_CASE(time_high_pass_filter) { std::cout << " time token: " << watch.ToString() << ' '; } -BOOST_AUTO_TEST_CASE(time_flagged_fitting) { - Image2DPtr image; - Mask2DPtr mask; - InitializeFlagged(image, mask); - - LocalFitMethod fitMethod; - TimeFrequencyData data(TimeFrequencyData::AmplitudePart, - Polarization::StokesI, image); - fitMethod.SetToWeightedAverage(10, 20, 2.5, 5.0); - fitMethod.Initialize(data); - Stopwatch watch(true); - for (size_t i = 0; i < fitMethod.TaskCount(); ++i) fitMethod.PerformFit(i); - std::cout << " time token: " << watch.ToString() << ' '; -} - BOOST_AUTO_TEST_CASE(time_flagged_high_pass_filter) { Image2DPtr image; Mask2DPtr mask;