C++ ArtifactSet类代码示例

/**
 * 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);
}

ArtifactSet RuleNode::changedInputArtifacts(const ArtifactSet &allCompatibleInputs,
                                            const ArtifactSet &explicitlyDependsOn,
                                            const ArtifactSet &auxiliaryInputs) const
{
    ArtifactSet changedInputArtifacts;
    if (explicitlyDependsOn != m_oldExplicitlyDependsOn)
        return allCompatibleInputs;
    if (!m_needsToConsiderChangedInputs)
        return changedInputArtifacts;

    for (Artifact * const artifact : explicitlyDependsOn) {
        if (artifact->timestamp() > m_lastApplicationTime)
            return allCompatibleInputs;
    }
    if (auxiliaryInputs != m_oldAuxiliaryInputs)
        return allCompatibleInputs;
    for (Artifact * const artifact : auxiliaryInputs) {
        if (artifact->timestamp() > m_lastApplicationTime)
            return allCompatibleInputs;
    }
    for (Artifact * const artifact : allCompatibleInputs) {
        if (artifact->timestamp() > m_lastApplicationTime)
            changedInputArtifacts.insert(artifact);
    }
    return changedInputArtifacts;
}

	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::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::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 TimeConvolutionAction::PerformFFTSincOperation(ArtifactSet &artifacts, Image2DPtr real, Image2DPtr imag) const
{
	fftw_complex
		*fftIn = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * real->Width()),
		*fftOut = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * real->Width());
	
	// FFTW plan routines are not thread safe, so lock.
	boost::mutex::scoped_lock lock(artifacts.IOMutex());
	fftw_plan
		fftPlanForward = fftw_plan_dft_1d(real->Width(), fftIn, fftOut, FFTW_FORWARD, FFTW_MEASURE),
		fftPlanBackward = fftw_plan_dft_1d(real->Width(), fftIn, fftOut, FFTW_BACKWARD, FFTW_MEASURE);
	lock.unlock();
	
	const size_t width = real->Width();

	const BandInfo band = artifacts.MetaData()->Band();
	for(unsigned y=0;y<real->Height();++y)
	{
		const numl_t sincScale = ActualSincScaleInSamples(artifacts, band.channels[y].frequencyHz);
		const numl_t limitFrequency = (numl_t) width / sincScale;
		if(y == real->Height()/2)
		{
			AOLogger::Debug << "Horizontal sinc scale: " << sincScale << " (filter scale: " << Angle::ToString(ActualSincScaleAsRaDecDist(artifacts, band.channels[y].frequencyHz)) << ")\n";
		}
		if(sincScale > 1.0)
		{
			for(unsigned x=0;x<width;++x)
			{
				fftIn[x][0] = real->Value(x, y);
				fftIn[x][1] = imag->Value(x, y);
			}
			
			fftw_execute_dft(fftPlanForward, fftIn, fftOut);
			size_t filterIndexSize = (limitFrequency > 1.0) ? (size_t) ceil(limitFrequency/2.0) : 1;
			// Remove the high frequencies [filterIndexSize : n-filterIndexSize]
			for(size_t f=filterIndexSize;f<width - filterIndexSize;++f)
			{
				fftOut[f][0] = 0.0;
				fftOut[f][1] = 0.0;
			}
			fftw_execute_dft(fftPlanBackward, fftOut, fftIn);
			
			const double n = width;
			for(unsigned x=0;x<width;++x)
			{
				real->SetValue(x, y, fftIn[x][0] / n);
				imag->SetValue(x, y, fftIn[x][1] / n);
			}
		}
	}
	fftw_free(fftIn);
	fftw_free(fftOut);
}

	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 Transformer::setupInputs(QScriptValue targetScriptValue, const ArtifactSet &inputs,
        const QString &defaultModuleName)
{
    const auto scriptEngine = static_cast<ScriptEngine *>(targetScriptValue.engine());
    QScriptValue scriptValue = translateInOutputs(scriptEngine, inputs, defaultModuleName);
    targetScriptValue.setProperty(StringConstants::inputsVar(), scriptValue);
    QScriptValue inputScriptValue;
    if (inputs.size() == 1) {
        Artifact *input = *inputs.cbegin();
        const FileTags &fileTags = input->fileTags();
        QBS_ASSERT(!fileTags.empty(), return);
        QScriptValue inputsForFileTag = scriptValue.property(fileTags.cbegin()->toString());
        inputScriptValue = inputsForFileTag.property(0);
    }

QScriptValue Transformer::translateInOutputs(ScriptEngine *scriptEngine,
                                             const ArtifactSet &artifacts,
                                             const QString &defaultModuleName)
{
    using TagArtifactsMap = QMap<QString, QList<Artifact*>>;
    TagArtifactsMap tagArtifactsMap;
    for (Artifact *artifact : artifacts)
        for (const FileTag &fileTag : artifact->fileTags())
            tagArtifactsMap[fileTag.toString()].push_back(artifact);
    for (TagArtifactsMap::Iterator it = tagArtifactsMap.begin(); it != tagArtifactsMap.end(); ++it)
        std::sort(it.value().begin(), it.value().end(), compareByFilePath);

    QScriptValue jsTagFiles = scriptEngine->newObject();
    for (TagArtifactsMap::const_iterator tag = tagArtifactsMap.constBegin(); tag != tagArtifactsMap.constEnd(); ++tag) {
        const QList<Artifact*> &artifacts = tag.value();
        QScriptValue jsFileConfig = scriptEngine->newArray(artifacts.size());
        int i = 0;
        for (Artifact * const artifact : artifacts) {
            jsFileConfig.setProperty(i++, translateFileConfig(scriptEngine, artifact,
                                                              defaultModuleName));
        }
        jsTagFiles.setProperty(tag.key(), jsFileConfig);
    }

    return jsTagFiles;
}

	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 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 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 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 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 ForEachBaselineAction::Perform(ArtifactSet &artifacts, ProgressListener &progress)
	{
		if(!artifacts.HasImageSet())
		{
			progress.OnStartTask(*this, 0, 1, "For each baseline (no image set)");
			progress.OnEndTask(*this);
			AOLogger::Warn <<
				"I executed a ForEachBaselineAction without an active imageset: something is\n"
				"likely wrong. Check your strategy and the input files.\n";
		} else if(_selection == Current)
		{
			ActionBlock::Perform(artifacts, progress);
		} else
		{
			ImageSet *imageSet = artifacts.ImageSet();
			MSImageSet *msImageSet = dynamic_cast<MSImageSet*>(imageSet);
			if(msImageSet != 0)
			{
				// Check memory usage
				ImageSetIndex *tempIndex = msImageSet->StartIndex();
				size_t timeStepCount = msImageSet->ObservationTimesVector(*tempIndex).size();
				delete tempIndex;
				size_t channelCount = msImageSet->GetBandInfo(0).channels.size();
				size_t estMemorySizePerThread = 8/*bp complex*/ * 4 /*polarizations*/ * timeStepCount * channelCount * 3 /* approx copies of the data that will be made in memory*/;
				AOLogger::Debug << "Estimate of memory each thread will use: " << estMemorySizePerThread/(1024*1024) << " MB.\n";
				size_t compThreadCount = _threadCount;
				if(compThreadCount > 0) --compThreadCount;
				if(estMemorySizePerThread * compThreadCount > 12ul*1024ul*1024ul*1024ul)
				{
					size_t maxThreads = (12ul * 1024ul * 1024ul * 1024ul) / estMemorySizePerThread;
					if(maxThreads < 1) maxThreads = 1;
					AOLogger::Warn <<
						"WARNING WARNING WARNING WARNING WARNING WARNING WARNING WARNING!\n"
						"This measurement set is TOO LARGE to be processed with " << _threadCount << " threads!\n" <<
						_threadCount << " threads would require " << ((estMemorySizePerThread*compThreadCount)/(1024*1024)) << " MB of memory approximately.\n"
						"Number of threads that will actually be used: " << maxThreads << "\n"
						"This might hurt performance a lot!\n\n";
					_threadCount = maxThreads;
				}
			}
			if(!_antennaeToSkip.empty())
			{
				AOLogger::Debug << "The following antenna's will be skipped: ";
				for(std::set<size_t>::const_iterator i=_antennaeToSkip.begin();i!=_antennaeToSkip.end(); ++i)
					AOLogger::Debug << (*i) << ' ';
				AOLogger::Debug <<'\n';
			}
			if(!_antennaeToInclude.empty())
			{
				AOLogger::Debug << "Only the following antenna's will be included: ";
				for(std::set<size_t>::const_iterator i=_antennaeToInclude.begin();i!=_antennaeToInclude.end(); ++i)
					AOLogger::Debug << (*i) << ' ';
				AOLogger::Debug <<'\n';
			}

			if(artifacts.MetaData() != 0)
			{
				_hasInitAntennae = true;
				if(artifacts.MetaData()->HasAntenna1())
					_initAntenna1 = artifacts.MetaData()->Antenna1();
				else
					_hasInitAntennae = false;
				if(artifacts.MetaData()->HasAntenna2())
					_initAntenna2 = artifacts.MetaData()->Antenna2();
				else
					_hasInitAntennae = false;
			}
			_artifacts = &artifacts;
			
			_initPartIndex = 0;
			_finishedBaselines = false;
			_baselineCount = 0;
			_baselineProgress = 0;
			_nextIndex = 0;
			
			// Count the baselines that are to be processed
			ImageSetIndex *iteratorIndex = imageSet->StartIndex();
			while(iteratorIndex->IsValid())
			{
				if(IsBaselineSelected(*iteratorIndex))
					++_baselineCount;
				iteratorIndex->Next();
			}
			delete iteratorIndex;
			AOLogger::Debug << "Will process " << _baselineCount << " baselines.\n";
			
			// Initialize thread data and threads
			_loopIndex = imageSet->StartIndex();
			_progressTaskNo = new int[_threadCount];
			_progressTaskCount = new int[_threadCount];
			progress.OnStartTask(*this, 0, 1, "Initializing");

			boost::thread_group threadGroup;
			ReaderFunction reader(*this);
			threadGroup.create_thread(reader);
			
			size_t mathThreads = mathThreadCount();
			for(unsigned i=0;i<mathThreads;++i)
			{
				PerformFunction function(*this, progress, i);
//.........这里部分代码省略.........