C++ Image2DCPtr::Width方法代码示例

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

Image2DPtr HighPassFilter::ApplyLowPass(const Image2DCPtr &image, const Mask2DCPtr &mask)
{
	initializeKernel();
	Image2DPtr
		outputImage = Image2D::CreateUnsetImagePtr(image->Width(), image->Height()),
		weights = Image2D::CreateUnsetImagePtr(image->Width(), image->Height());
	setFlaggedValuesToZeroAndMakeWeightsSSE(image, outputImage, mask, weights);
	applyLowPassSSE(outputImage);
	applyLowPassSSE(weights);
	elementWiseDivideSSE(outputImage, weights);
	weights.reset();
	return outputImage;
}

void Compress::Write(std::ofstream &stream, Image2DCPtr image, Mask2DCPtr mask)
{
	const num_t
		max = ThresholdTools::MaxValue(image, mask),
		min = ThresholdTools::MinValue(image, mask),
		mid = (min + max) / 2.0;
	const num_t normalizeFactor = (num_t) ((2<<22) + ((2<<22)-1)) / (max - min);
	const uint32_t
		width = image->Width(),
		height = image->Height();
	const char mode = 0;

	stream.write(reinterpret_cast<const char*>(&max), sizeof(max));
	stream.write(reinterpret_cast<const char*>(&min), sizeof(min));
	stream.write(reinterpret_cast<const char*>(&width), sizeof(width));
	stream.write(reinterpret_cast<const char*>(&height), sizeof(height));
	stream.write(&mode, sizeof(mode));

	for(unsigned y=0;y<height;++y)
	{
		for(unsigned x=0;x<width;++x)
		{
			if(!mask->Value(x, y))
			{
				int32_t value = (int32_t) round((image->Value(x, y) - mid) * normalizeFactor);
				stream.write(reinterpret_cast<char*>(&value)+1, 3);
			}
		}
	}
}

void RFIGuiController::PlotPowerTime()
{
	if(IsImageLoaded())
	{
		Plot2D &plot = _plotManager->NewPlot2D("Power over time");
		plot.SetLogarithmicYAxis(true);

		TimeFrequencyData activeData = ActiveData();
		Image2DCPtr image = activeData.GetSingleImage();
		Mask2DPtr mask =
			Mask2D::CreateSetMaskPtr<false>(image->Width(), image->Height());
		Plot2DPointSet &totalPlot = plot.StartLine("Total");
		RFIPlots::MakePowerTimePlot(totalPlot, image, mask, MetaData());

		mask = Mask2D::CreateCopy(activeData.GetSingleMask());
		if(!mask->AllFalse())
		{
			Plot2DPointSet &uncontaminatedPlot = plot.StartLine("Uncontaminated");
			RFIPlots::MakePowerTimePlot(uncontaminatedPlot, image, mask, MetaData());
	
			mask->Invert();
			Plot2DPointSet &rfiPlot = plot.StartLine("RFI");
			RFIPlots::MakePowerTimePlot(rfiPlot, image, mask, MetaData());
		}

		_plotManager->Update();
	}
}

void RFIGuiController::PlotPowerSNR()
{
	Image2DCPtr
		image = ActiveData().GetSingleImage(),
		model = RevisedData().GetSingleImage();
	if(IsImageLoaded())
	{
		Plot2D &plot = _plotManager->NewPlot2D("SNR spectrum");
		plot.SetLogarithmicYAxis(true);

		Mask2DPtr mask =
			Mask2D::CreateSetMaskPtr<false>(image->Width(), image->Height());
		Plot2DPointSet &totalPlot = plot.StartLine("Total");
		RFIPlots::MakeSNRSpectrumPlot(totalPlot, image, model, mask);

		mask = Mask2D::CreateCopy(ActiveData().GetSingleMask());
		if(!mask->AllFalse())
		{
			Plot2DPointSet &uncontaminatedPlot = plot.StartLine("Uncontaminated");
			RFIPlots::MakeSNRSpectrumPlot(uncontaminatedPlot, image, model, mask);
	
			mask->Invert();
			Plot2DPointSet &rfiPlot = plot.StartLine("RFI");
			RFIPlots::MakeSNRSpectrumPlot(rfiPlot, image, model, mask);
		}

		_plotManager->Update();
	}
}

void RFIGuiController::PlotPowerRMS()
{
	if(IsImageLoaded())
	{
		Plot2D &plot = _plotManager->NewPlot2D("Spectrum RMS");
		plot.SetLogarithmicYAxis(true);

		TimeFrequencyData activeData = ActiveData();
		Image2DCPtr image = activeData.GetSingleImage();
		Mask2DPtr mask =
			Mask2D::CreateSetMaskPtr<false>(image->Width(), image->Height());
		Plot2DPointSet &beforeSet = plot.StartLine("Before");
		RFIPlots::MakeRMSSpectrumPlot(beforeSet, image, mask);

		mask = Mask2D::CreateCopy(activeData.GetSingleMask());
		if(!mask->AllFalse())
		{
			Plot2DPointSet &afterSet = plot.StartLine("After");
			RFIPlots::MakeRMSSpectrumPlot(afterSet, image, mask);
	
			//mask->Invert();
			//Plot2DPointSet &rfiSet = plot.StartLine("RFI");
			//RFIPlots::MakeRMSSpectrumPlot(rfiSet, _timeFrequencyWidget.Image(), mask);
		}

		_plotManager->Update();
	}
}

	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 ThresholdMitigater::HorizontalSumThreshold(Image2DCPtr input, Mask2DPtr mask, num_t threshold)
{
	if(Length <= input->Width())
	{
		size_t width = input->Width()-Length+1; 
		for(size_t y=0;y<input->Height();++y) {
			for(size_t x=0;x<width;++x) {
				num_t sum = 0.0;
				size_t count = 0;
				for(size_t i=0;i<Length;++i) {
					if(!mask->Value(x+i, y)) {
						sum += input->Value(x+i, y);
						count++;
					}
				}
				if(count>0 && fabs(sum/count) > threshold) {
					for(size_t i=0;i<Length;++i)
						mask->SetValue(x + i, y, true);
				}
			}
		}
	}
}

void HighPassFilter::setFlaggedValuesToZeroAndMakeWeights(const Image2DCPtr &inputImage, const Image2DPtr &outputImage, const Mask2DCPtr &inputMask, const Image2DPtr &weightsOutput)
{
	const size_t width = inputImage->Width();
	for(size_t y=0;y<inputImage->Height();++y)
	{
		for(size_t x=0;x<width;++x)
		{
			if(inputMask->Value(x, y) || !isfinite(inputImage->Value(x, y)))
			{
				outputImage->SetValue(x, y, 0.0);
				weightsOutput->SetValue(x, y, 0.0);
			} else {
				outputImage->SetValue(x, y, inputImage->Value(x, y));
				weightsOutput->SetValue(x, y, 1.0);
			}
		}
	}
}

void ThresholdMitigater::VerticalVarThreshold(Image2DCPtr input, Mask2DPtr mask, size_t length, num_t threshold)
{
	size_t height = input->Height()-length+1; 
	for(size_t y=0;y<height;++y) {
		for(size_t x=0;x<input->Width();++x) {
			bool flag = true;
			for(size_t i=0;i<length;++i) {
				if(input->Value(x, y+i) <= threshold && input->Value(x, y+i) >= -threshold) {
					flag = false;
					break;
				}
			}
			if(flag) {
				for(size_t i=0;i<length;++i)
					mask->SetValue(x, y + i, true);
			}
		}
	}
}

void ThresholdMitigater::HorizontalVarThreshold(Image2DCPtr input, Mask2DPtr mask, size_t length, num_t threshold)
{
	size_t width = input->Width()-length+1;
	for(size_t y=0;y<input->Height();++y) {
		for(size_t x=0;x<width;++x) {
			bool flag = true;
			for(size_t i=0;i<length;++i) {
				if(input->Value(x+i, y) < threshold && input->Value(x+i, y) > -threshold) {
					flag = false;
					break;
				}
			}
			if(flag) {
				for(size_t i=0;i<length;++i)
					mask->SetValue(x + i, y, true);
			}
		}
	}
}

void HighPassFilter::setFlaggedValuesToZeroAndMakeWeightsSSE(const Image2DCPtr &inputImage, const Image2DPtr &outputImage, const Mask2DCPtr &inputMask, const Image2DPtr &weightsOutput)
{
	const size_t width = inputImage->Width();
	const __m128i zero4i = _mm_set_epi32(0, 0, 0, 0);
	const __m128 zero4 = _mm_set_ps(0.0, 0.0, 0.0, 0.0);
	const __m128 one4 = _mm_set_ps(1.0, 1.0, 1.0, 1.0);
	for(size_t y=0;y<inputImage->Height();++y)
	{
		const bool *rowPtr = inputMask->ValuePtr(0, y);
		const float *inputPtr = inputImage->ValuePtr(0, y);
		float *outputPtr = outputImage->ValuePtr(0, y);
		float *weightsPtr = weightsOutput->ValuePtr(0, y);
		const float *end = inputPtr + width;
		while(inputPtr < end)
		{
			
			// Assign each integer to one bool in the mask
			// Convert false to 0xFFFFFFFF and true to 0
			__m128 conditionMask = _mm_castsi128_ps(
				_mm_cmpeq_epi32(_mm_set_epi32(rowPtr[3] || !isfinite(inputPtr[3]), rowPtr[2] || !isfinite(inputPtr[2]),
																			rowPtr[1] || !isfinite(inputPtr[1]), rowPtr[0] || !isfinite(inputPtr[0])),
												zero4i));
			
			_mm_store_ps(weightsPtr, _mm_or_ps(
				_mm_and_ps(conditionMask, one4),
				_mm_andnot_ps(conditionMask, zero4)
			));
			_mm_store_ps(outputPtr, _mm_or_ps(
				_mm_and_ps(conditionMask, _mm_load_ps(inputPtr)),
				_mm_andnot_ps(conditionMask, zero4)
			));
			
			rowPtr += 4;
			outputPtr += 4;
			inputPtr += 4;
			weightsPtr += 4;
		}
	}
}

void Compress::WriteSubtractFrequencies(std::ofstream &stream, Image2DCPtr image, Mask2DCPtr mask)
{
	const num_t
		max = ThresholdTools::MaxValue(image, mask),
		min = ThresholdTools::MinValue(image, mask);
	const num_t normalizeFactor = (num_t) ((2<<22) + ((2<<22)-1)) / (max - min);
	//const num_t normalizeFactor = 256.0;
	const uint32_t
		width = image->Width(),
		height = image->Height();
	const char mode = 1;

	stream.write(reinterpret_cast<const char*>(&max), sizeof(max));
	stream.write(reinterpret_cast<const char*>(&min), sizeof(min));
	stream.write(reinterpret_cast<const char*>(&width), sizeof(width));
	stream.write(reinterpret_cast<const char*>(&height), sizeof(height));
	stream.write(&mode, sizeof(mode));

	std::vector<int32_t> basis(width);
	for(size_t x=0;x<width;++x)
	{
		SampleRowPtr row = SampleRow::CreateFromColumn(image, x);
		basis[x] = (int32_t) round(row->Median() * normalizeFactor);
	}
	stream.write(reinterpret_cast<char*>(&basis[0]), sizeof(basis));

	for(unsigned y=0;y<height;++y)
	{
		for(unsigned x=0;x<width;++x)
		{
			if(!mask->Value(x, y))
			{
				int32_t value = (int32_t) (round(image->Value(x, y) * normalizeFactor) - basis[x]);
				stream.write(reinterpret_cast<char*>(&value)+1, 3);
			}
		}
	}
}

void ImagePlaneWindow::onPlotVertically()
{
	if(_heatMapPlot.HasImage())
	{
		Plot plot("Image-vertical-axis.pdf");
		plot.SetXAxisText("Declination index");
		plot.SetYAxisText("Amplitude");
		//plot.SetLogScale(false, true);
		plot.StartLine();
		Image2DCPtr image = _heatMapPlot.Image();
		for(size_t y=0;y<image->Height();++y)
		{
			num_t sum = 0.0;
			for(size_t x=0;x<image->Width();++x)
			{
				sum += image->Value(x, y);
			}
			plot.PushDataPoint(y, sum);
		}
		plot.Close();
		plot.Show();
	}
}

void FrequencyPowerPlot::Add(class TimeFrequencyData &data, TimeFrequencyMetaDataCPtr meta)
{
	Image2DCPtr image = data.GetSingleImage();
	Mask2DCPtr mask = data.GetSingleMask();
	for(size_t y=0;y<image->Height();++y)
	{
		double frequency = meta->Band().channels[y].frequencyHz;
		size_t count = 0;
		long double value = 0.0L;

		for(size_t x=0;x<image->Width();++x)
		{
			if(!mask->Value(x, y))
			{
				++count;
				value += image->Value(x, y);
			}
		}
		MapItem &item = _values[frequency];
		item.total += value;
		item.count += count;
	}
}