C++ ArtifactSet::ContaminatedData方法代码示例

/**
 * Automatic selection selects all timesteps which RMS is higher than some value relative to the stddev of
 * all timesteps.
 */
void TimeSelectionAction::AutomaticSelection(ArtifactSet &artifacts)
{
	Image2DCPtr image = artifacts.ContaminatedData().GetSingleImage();
	SampleRowPtr timesteps = SampleRow::CreateEmpty(image->Width());
	Mask2DPtr mask = Mask2D::CreateCopy(artifacts.ContaminatedData().GetSingleMask());
	for(size_t x=0;x<image->Width();++x)
	{
		SampleRowPtr row = SampleRow::CreateFromColumnWithMissings(image, mask, x);
		timesteps->SetValue(x, row->RMSWithMissings());
	}
	bool change;
	MedianWindow<num_t>::SubtractMedian(timesteps, 512);
	do {
		num_t median = 0.0;
		num_t stddev = timesteps->StdDevWithMissings(0.0);
		change = false;
		for(size_t x=0;x<timesteps->Size();++x)
		{
			if(!timesteps->ValueIsMissing(x) && (timesteps->Value(x) - median > stddev * _threshold || median - timesteps->Value(x) > stddev * _threshold))
			{
				mask->SetAllVertically<true>(x);
				timesteps->SetValueMissing(x);
				change = true;
			}
		}
	} while(change);
	artifacts.ContaminatedData().SetGlobalMask(mask);
}

	void SlidingWindowFitAction::Perform(ArtifactSet &artifacts, class ProgressListener &listener)
	{
		LocalFitMethod method;
		switch(_parameters.method)
		{
			case SlidingWindowFitParameters::None:
				method.SetToNone();
				break;
			case SlidingWindowFitParameters::Average:
				method.SetToAverage(
					_parameters.timeDirectionWindowSize,
					_parameters.frequencyDirectionWindowSize);
				break;
			case SlidingWindowFitParameters::GaussianWeightedAverage:
				method.SetToWeightedAverage(
					_parameters.timeDirectionWindowSize,
					_parameters.frequencyDirectionWindowSize,
					_parameters.timeDirectionKernelSize,
					_parameters.frequencyDirectionKernelSize);
				break;
			case SlidingWindowFitParameters::Median:
				method.SetToMedianFilter(
					_parameters.timeDirectionWindowSize,
					_parameters.frequencyDirectionWindowSize);
				break;
			case SlidingWindowFitParameters::Minimum:
				method.SetToMinimumFilter(
					_parameters.timeDirectionWindowSize,
					_parameters.frequencyDirectionWindowSize);
				break;
		}

		method.Initialize(artifacts.ContaminatedData());
		
		size_t taskCount = method.TaskCount();
		for(size_t i=0;i<taskCount;++i)
		{
			method.PerformFit(i);
			listener.OnProgress(*this, i+1, taskCount);
		}
		TimeFrequencyData newRevisedData = method.Background();
		newRevisedData.SetMask(artifacts.RevisedData());

		TimeFrequencyData *contaminatedData =
			TimeFrequencyData::CreateTFDataFromDiff(artifacts.ContaminatedData(), newRevisedData);
		contaminatedData->SetMask(artifacts.ContaminatedData());

		artifacts.SetRevisedData(newRevisedData);
		artifacts.SetContaminatedData(*contaminatedData);

		delete contaminatedData;
	}

	void PlotAction::plotPolarizationFlagCounts(ArtifactSet &artifacts)
	{
		if(artifacts.PolarizationStatistics() == 0)
			throw BadUsageException("No polarization statistics in the artifact set");

		TimeFrequencyData &data = artifacts.ContaminatedData();
		artifacts.PolarizationStatistics()->Add(data);
	}

	void PlotAction::plotFrequencyPower(ArtifactSet &artifacts)
	{
		if(artifacts.FrequencyPowerPlot() == 0)
			throw BadUsageException("No frequency power plot in the artifact set");

		TimeFrequencyData &data = artifacts.ContaminatedData();
		TimeFrequencyMetaDataCPtr meta = artifacts.MetaData();
		artifacts.FrequencyPowerPlot()->Add(data, meta);
	}

	void PlotAction::plotAntennaFlagCounts(ArtifactSet &artifacts)
	{
		if(artifacts.AntennaFlagCountPlot() == 0)
			throw BadUsageException("No antenna flag count plot in the artifact set");

		TimeFrequencyData &data = artifacts.ContaminatedData();
		TimeFrequencyMetaDataCPtr meta = artifacts.MetaData();
		artifacts.AntennaFlagCountPlot()->Add(data, meta);
	}

	void PlotAction::plotSpectrumPerBaseline(ArtifactSet &artifacts)
	{
		if(artifacts.FrequencyPowerPlot() == 0)
			throw BadUsageException("No frequency power plot in the artifact set");

		TimeFrequencyData &data = artifacts.ContaminatedData();
		TimeFrequencyMetaDataCPtr meta = artifacts.MetaData();
		artifacts.FrequencyPowerPlot()->SetLogYAxis(_logYAxis);
		artifacts.FrequencyPowerPlot()->StartNewLine(meta->Antenna1().name + " x " + meta->Antenna2().name);
		artifacts.FrequencyPowerPlot()->Add(data, meta);
	}

	void ImagerAction::Perform(ArtifactSet &artifacts, ProgressListener &progress)
	{
		boost::mutex::scoped_lock lock(_imagerMutex);
		UVImager *imager = artifacts.Imager();
		if(imager == 0)
			throw BadUsageException("No imager available to create image.");
		TimeFrequencyData &data = artifacts.ContaminatedData();
		TimeFrequencyMetaDataCPtr metaData = artifacts.MetaData();
		if(data.PolarisationCount() > 1)
		{
			TimeFrequencyData *tmp = data.CreateTFData(StokesIPolarisation);
			data = *tmp;
			delete tmp;
		}
		
		bool btPlaneImager = true;
		if(btPlaneImager)
		{
			typedef double ImagerNumeric;
			BaselineTimePlaneImager<ImagerNumeric> btImager;
			BandInfo band = metaData->Band();
			Image2DCPtr
				inputReal = data.GetRealPart(),
				inputImag = data.GetImaginaryPart();
			Mask2DCPtr mask = data.GetSingleMask();
			size_t width = inputReal->Width();
			
			for(size_t t=0;t!=width;++t)
			{
				UVW uvw = metaData->UVW()[t];
				size_t channelCount = inputReal->Height();
				std::vector<std::complex<ImagerNumeric> > data(channelCount);
				for(size_t ch=0;ch!=channelCount;++ch) {
					if(mask->Value(t, ch))
						data[ch] = std::complex<ImagerNumeric>(0.0, 0.0);
					else
						data[ch] = std::complex<ImagerNumeric>(inputReal->Value(t, ch), inputImag->Value(t, ch));
				}
				
				btImager.Image(uvw.u, uvw.v, uvw.w, band.channels[0].frequencyHz, band.channels[1].frequencyHz-band.channels[0].frequencyHz, channelCount, &(data[0]), imager->FTReal());
			}
		} else {
			progress.OnStartTask(*this, 0, 1, "Imaging baseline");
			for(size_t y=0;y<data.ImageHeight();++y)
			{
				imager->Image(data, metaData, y);
				progress.OnProgress(*this, y, data.ImageHeight());
			}
			progress.OnEndTask(*this);
		}
	}

	void SVDAction::Perform(ArtifactSet &artifacts, class ProgressListener &listener)
	{
		SVDMitigater mitigater;
		mitigater.Initialize(artifacts.ContaminatedData());
		mitigater.SetRemoveCount(_singularValueCount);
		for(size_t i=0;i<mitigater.TaskCount();++i)
		{
			mitigater.PerformFit(i);
			listener.OnProgress(*this, i+1, mitigater.TaskCount());
		}

		TimeFrequencyData newRevisedData = mitigater.Background();
		newRevisedData.SetMask(artifacts.RevisedData());

		TimeFrequencyData *contaminatedData =
			TimeFrequencyData::CreateTFDataFromDiff(artifacts.ContaminatedData(), newRevisedData);
		contaminatedData->SetMask(artifacts.ContaminatedData());

		artifacts.SetRevisedData(newRevisedData);
		artifacts.SetContaminatedData(*contaminatedData);

		delete contaminatedData;
	}

	void PlotAction::plotAntennaFlagCounts(ArtifactSet &artifacts)
	{
		if(artifacts.AntennaFlagCountPlot() == 0)
			throw BadUsageException("No antenna flag count plot in the artifact set");
		
		if(artifacts.HasMetaData() && artifacts.MetaData()->HasAntenna1() && artifacts.MetaData()->HasAntenna2())
		{
			TimeFrequencyData &data = artifacts.ContaminatedData();
			TimeFrequencyMetaDataCPtr meta = artifacts.MetaData();
			artifacts.AntennaFlagCountPlot()->Add(data, meta);
		} else {
			AOLogger::Warn << "The strategy contains an action that makes an antenna plot, but the image set did not provide meta data.\n"
				"Plot will not be made.\n";
		}
	}

	void PlotAction::plotBaselineRMS(ArtifactSet &artifacts)
	{
		if(artifacts.PolarizationStatistics() == 0)
			throw BadUsageException("No polarization statistics in the artifact set");

		TimeFrequencyData &data = artifacts.ContaminatedData();
		TimeFrequencyMetaDataCPtr metaData = artifacts.MetaData();
		double rms = 0.0;
		for(unsigned i=0;i<data.PolarisationCount();++i)
		{
			TimeFrequencyData *polarisation = data.CreateTFDataFromPolarisationIndex(i);
			Mask2DCPtr mask = polarisation->GetSingleMask();
			for(unsigned j=0;j<polarisation->ImageCount();++j)
			{
				Image2DCPtr image = polarisation->GetImage(j);
				rms += ThresholdTools::RMS(image, mask);
			}
			delete polarisation;
		}
		rms /= data.PolarisationCount();
		;
		AOLogger::Info << "RMS of " << metaData->Antenna1().name << " x " << metaData->Antenna2().name << ": "
			<< rms << '\n';
	}

	void SpatialCompositionAction::Perform(ArtifactSet &artifacts, ProgressListener &progress)
	{
		size_t imageCount = artifacts.ContaminatedData().ImageCount();
		std::vector<Image2DPtr> images(imageCount);
		for(size_t p=0;p<imageCount;++p)
			images[p] = Image2D::CreateZeroImagePtr(artifacts.ContaminatedData().ImageWidth(), artifacts.ContaminatedData().ImageHeight());

		std::string filename = artifacts.ImageSet()->File();
		SpatialMSImageSet set(filename);
		ImageSetIndex *index = set.StartIndex();
		size_t progressStep = 0, totalProgress = artifacts.ContaminatedData().ImageWidth() * artifacts.ContaminatedData().ImageHeight()/256;
		while(index->IsValid())
		{
			TimeFrequencyData *data = set.LoadData(*index);
			SpatialMatrixMetaData metaData(set.SpatialMetaData(*index));
			for(size_t p=0;p!=imageCount;++p)
			{
				switch(_operation)
				{
					case SumCrossCorrelationsOperation:
						images[p]->SetValue(metaData.TimeIndex(), metaData.ChannelIndex(), sumCrossCorrelations(data->GetImage(p)));
						break;
					case SumAutoCorrelationsOperation:
						images[p]->SetValue(metaData.TimeIndex(), metaData.ChannelIndex(), sumAutoCorrelations(data->GetImage(p)));
						break;
					case EigenvalueDecompositionOperation: {
						num_t value = eigenvalue(data->GetImage(p), data->GetImage(p+1));
						images[p]->SetValue(metaData.TimeIndex(), metaData.ChannelIndex(), value);
						images[p+1]->SetValue(metaData.TimeIndex(), metaData.ChannelIndex(), 0.0);
						++p;
						} break;
					case EigenvalueRemovalOperation: {
						std::pair<num_t, num_t> value = removeEigenvalue(data->GetImage(p), data->GetImage(p+1));
						images[p]->SetValue(metaData.TimeIndex(), metaData.ChannelIndex(), value.first);
						images[p+1]->SetValue(metaData.TimeIndex(), metaData.ChannelIndex(), value.second);
						++p;
						} break;
				}
			}
			delete data;
			index->Next();
			++progressStep;
			progress.OnProgress(*this, progressStep/256, totalProgress);
		}
		delete index;

		TimeFrequencyData newRevisedData = artifacts.RevisedData();
		for(size_t p=0;p<imageCount;++p)
			newRevisedData.SetImage(p, images[p]);
		
		newRevisedData.SetMask(artifacts.RevisedData());

		TimeFrequencyData *contaminatedData =
			TimeFrequencyData::CreateTFDataFromDiff(artifacts.ContaminatedData(), newRevisedData);
		contaminatedData->SetMask(artifacts.ContaminatedData());

		artifacts.SetRevisedData(newRevisedData);
		artifacts.SetContaminatedData(*contaminatedData);

		delete contaminatedData;

	}