C++ Mask2DPtr::SetValue方法代码示例

void StatisticalFlagger::EnlargeFlags(Mask2DPtr mask, size_t timeSize, size_t frequencySize)
{
	Mask2DCPtr old = Mask2D::CreateCopy(mask);
	for(size_t y=0;y<mask->Height();++y)
	{
		size_t top, bottom;
		if(y > frequencySize)
			top = y - frequencySize;
		else
			top = 0;
		if(y + frequencySize < mask->Height() - 1)
			bottom = y + frequencySize;
		else
			bottom = mask->Height() - 1;
		
		for(size_t x=0;x<mask->Width();++x)
		{
			size_t left, right;
			if(x > timeSize)
				left = x - timeSize;
			else
				left = 0;
			if(x + timeSize < mask->Width() - 1)
				right = x + timeSize;
			else
				right = mask->Width() - 1;
			
			if(SquareContainsFlag(old, left, top, right, bottom))
				mask->SetValue(x, y, true);
		}
	}
}

void StatisticalFlagger::FlagFrequency(Mask2DPtr mask, size_t y)
{
	for(size_t x=0;x<mask->Width();++x)
	{
		mask->SetValue(x, y, true);
	}
}

std::pair<TimeFrequencyData,TimeFrequencyMetaDataPtr> RSPReader::ReadSingleBeamlet(unsigned long timestepStart, unsigned long timestepEnd, unsigned beamletCount, unsigned beamletIndex)
{
	std::pair<TimeFrequencyData,TimeFrequencyMetaDataPtr> data = ReadAllBeamlets(timestepStart, timestepEnd, beamletCount);
	
	const unsigned width = timestepEnd - timestepStart;
	Image2DPtr realX = Image2D::CreateZeroImagePtr(width, 1);
	Image2DPtr imaginaryX = Image2D::CreateZeroImagePtr(width, 1);
	Image2DPtr realY = Image2D::CreateZeroImagePtr(width, 1);
	Image2DPtr imaginaryY = Image2D::CreateZeroImagePtr(width, 1);
	Mask2DPtr mask = Mask2D::CreateUnsetMaskPtr(width, 1);
	
	TimeFrequencyData allX = data.first.Make(Polarization::XX);
	TimeFrequencyData allY = data.first.Make(Polarization::YY);
	Image2DCPtr xr = allX.GetRealPart();
	Image2DCPtr xi = allX.GetImaginaryPart();
	Image2DCPtr yr = allY.GetRealPart();
	Image2DCPtr yi = allY.GetImaginaryPart();
	Mask2DCPtr maskWithBeamlets = data.first.GetSingleMask();
	
	for(unsigned x=0;x<width;++x)
	{
		realX->SetValue(x, 0, xr->Value(x, beamletIndex));
		imaginaryX->SetValue(x, 0, xi->Value(x, beamletIndex));
		realY->SetValue(x, 0, yr->Value(x, beamletIndex));
		imaginaryY->SetValue(x, 0, yi->Value(x, beamletIndex));
		mask->SetValue(x, 0, maskWithBeamlets->Value(x, beamletIndex));
	}
	data.first = TimeFrequencyData(Polarization::XX, realX, imaginaryX, Polarization::YY, realY, imaginaryY);
	data.first.SetGlobalMask(mask);
	BandInfo band = data.second->Band();
	band.channels[0] = data.second->Band().channels[beamletIndex];
	band.channels.resize(1);
	data.second->SetBand(band);
	return data;
}

void StatisticalFlagger::FlagTime(Mask2DPtr mask, size_t x)
{
	for(size_t y=0;y<mask->Height();++y)
	{
		mask->SetValue(x, y, true);
	}
}

void Morphology::floodFill(Mask2DCPtr mask, SegmentedImagePtr output, Mask2DPtr *matrices, size_t x, size_t y, size_t z, size_t value, int **hCounts, int **vCounts)
{
	std::stack<MorphologyPoint3D> points;
	MorphologyPoint3D startPoint;
	startPoint.x = x;
	startPoint.y = y;
	startPoint.z = z;
	points.push(startPoint);
	do {
		MorphologyPoint3D p = points.top();
		points.pop();
		if(mask->Value(p.x, p.y))
		{
			if(output->Value(p.x, p.y) == 0)
			{
				output->SetValue(p.x, p.y, value);
			} else {
				// now we need to decide whether to change this sample to the new segment or not
				if(hCounts[p.y][p.x] < vCounts[p.y][p.x] && p.z == 2)
					output->SetValue(p.x, p.y, value);
			}
		}
		Mask2DPtr matrix = matrices[p.z];
		matrix->SetValue(p.x, p.y, false);
		if((p.z == 0 || p.z == 2) && matrices[1]->Value(p.x,p.y))
		{
			MorphologyPoint3D newP;
			newP.x = p.x; newP.y = p.y; newP.z = 1;
			points.push(newP);
		}
		if(p.x > 0 && matrix->Value(p.x-1,p.y))
		{
			MorphologyPoint3D newP;
			newP.x = p.x-1; newP.y = p.y; newP.z = p.z;
			points.push(newP);
		}
		if(p.x < mask->Width()-1 && matrix->Value(p.x+1,p.y))
		{
			MorphologyPoint3D newP;
			newP.x = p.x+1; newP.y = p.y; newP.z = p.z; newP.z = p.z;
			points.push(newP);
		}
		if(p.y > 0 && matrix->Value(p.x,p.y-1))
		{
			MorphologyPoint3D newP;
			newP.x = p.x; newP.y = p.y-1; newP.z = p.z;
			points.push(newP);
		}
		if(p.y < mask->Height()-1 && matrix->Value(p.x,p.y+1))
		{
			MorphologyPoint3D newP;
			newP.x = p.x; newP.y = p.y+1; newP.z = p.z;
			points.push(newP);
		}
	} while(points.size() != 0);
}

void StatisticalFlagger::DilateFlagsHorizontally(Mask2DPtr mask, size_t timeSize)
{
	if(timeSize != 0)
	{
		Mask2DPtr destination = Mask2D::CreateUnsetMaskPtr(mask->Width(), mask->Height());
		if(timeSize > mask->Width()) timeSize = mask->Width();
		const int intSize = (int) timeSize;
		
		for(size_t y=0;y<mask->Height();++y)
		{
			int dist = intSize + 1;
			for(size_t x=0;x<timeSize;++x)
			{
				if(mask->Value(x, y))
					dist = - intSize;
				dist++;
			}
			for(size_t x=0;x<mask->Width() - timeSize;++x)
			{
				if(mask->Value(x + timeSize, y))
					dist = -intSize;
				if(dist <= intSize)
				{
					destination->SetValue(x, y, true);
					dist++;
				} else {
					destination->SetValue(x, y, false);
				}
			}
			for(size_t x=mask->Width() - timeSize;x<mask->Width();++x)
			{
				if(dist <= intSize)
				{
					destination->SetValue(x, y, true);
					dist++;
				} else {
					destination->SetValue(x, y, false);
				}
			}
		}
		mask->Swap(destination);
	}
}

void StatisticalFlagger::ApplyMarksInFrequency(Mask2DPtr mask, int **flagMarks)
{
	for(size_t x=0;x<mask->Width();++x)
	{
		int startedCount = 0;
		for(size_t y=0;y<mask->Height();++y)
		{
			startedCount += flagMarks[y][x];
			if(startedCount > 0)
				mask->SetValue(x, y, true);
		}
	}
}

void ThresholdMitigater::VerticalSumThresholdLarge(Image2DCPtr input, Mask2DPtr mask, num_t threshold)
{
	Mask2DPtr maskCopy = Mask2D::CreateCopy(mask);
	const size_t width = mask->Width(), height = mask->Height();
	if(Length <= height)
	{
		for(size_t x=0;x<width;++x)
		{
			num_t sum = 0.0;
			size_t count = 0, yTop, yBottom;

			for(yBottom=0;yBottom<Length-1;++yBottom)
			{
				if(!mask->Value(x, yBottom))
				{
					sum += input->Value(x, yBottom);
					++count;
				}
			}

			yTop = 0;
			while(yBottom < height)
			{
				// add the sample at the bottom
				if(!mask->Value(x, yBottom))
				{
					sum += input->Value(x, yBottom);
					++count;
				}
				// Check
				if(count>0 && fabs(sum/count) > threshold)
				{
					for(size_t i=0;i<Length;++i)
						maskCopy->SetValue(x, yTop + i, true);
				}
				// subtract the sample at the top
				if(!mask->Value(x, yTop))
				{
					sum -= input->Value(x, yTop);
					--count;
				}
				++yTop;
				++yBottom;
			}
		}
	}
	mask->Swap(maskCopy);
}

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

  void BHFitsImageSet::loadImageData(TimeFrequencyData &data, const TimeFrequencyMetaDataPtr &metaData, const BHFitsImageSetIndex &index)
  {
		std::vector<num_t> buffer(_width * _height);
		_file->ReadCurrentImageData(0, &buffer[0], _width * _height);
		
		int
			rangeStart = _timeRanges[index._imageIndex].start, 
			rangeEnd = _timeRanges[index._imageIndex].end;
		Image2DPtr image = Image2D::CreateZeroImagePtr(rangeEnd-rangeStart, _height);

		std::vector<num_t>::const_iterator bufferPtr = buffer.begin() + _height*rangeStart;
		for(int x=rangeStart; x!=rangeEnd; ++x)
		{
			for(int y=0; y!=_height; ++y)
			{
				image->SetValue(x-rangeStart, y, *bufferPtr);
				++bufferPtr;
			}
		}
		data = TimeFrequencyData(TimeFrequencyData::AmplitudePart, SinglePolarisation, image);

		try {
			FitsFile flagFile(flagFilePath());
			flagFile.Open(FitsFile::ReadOnlyMode);
			flagFile.ReadCurrentImageData(0, &buffer[0], _width * _height);
			bufferPtr = buffer.begin() + _height*rangeStart;
			Mask2DPtr mask = Mask2D::CreateUnsetMaskPtr(rangeEnd-rangeStart, _height);
			for(int x=rangeStart; x!=rangeEnd; ++x)
			{
				for(int y=0; y!=_height; ++y)
				{
					bool flag = false;
					if(*bufferPtr == 0.0)
						flag = false;
					else if(*bufferPtr == 1.0)
						flag = true;
					else std::runtime_error("Expecting a flag file with only ones and zeros, but this file contained other values.");
					mask->SetValue(x-rangeStart, y, flag);
					++bufferPtr;
				}
			}
			data.SetGlobalMask(mask);
		} catch(std::exception &)
		{
			// Flag file could not be read; probably does not exist. Ignore this, flags will be initialized to false.
		}

		double
			frequencyDelta = _file->GetDoubleKeywordValue("CDELT1"),
			timeDelta = _file->GetDoubleKeywordValue("CDELT2");
		BandInfo band;
		for(int ch=0; ch!=_height; ++ch)
		{
			ChannelInfo channel;
			channel.frequencyHz = ch * frequencyDelta * 1000000.0;
			band.channels.push_back(channel);
		}
		metaData->SetBand(band);

		const int rangeWidth = rangeEnd-rangeStart;
		std::vector<double> observationTimes(rangeWidth);
		for(int t=0; t!=rangeWidth; ++t)
			observationTimes[t] = (t + rangeStart) * timeDelta;
		metaData->SetObservationTimes(observationTimes);

		AntennaInfo antennaInfo;
		antennaInfo.id = 0;
		antennaInfo.name = RangeName(index._imageIndex);
		antennaInfo.diameter = 0.0;
		antennaInfo.mount = "Unknown";
		antennaInfo.station = GetTelescopeName();
		metaData->SetAntenna1(antennaInfo);
		metaData->SetAntenna2(antennaInfo);
  }

std::pair<TimeFrequencyData,TimeFrequencyMetaDataPtr> RSPReader::ReadChannelBeamlet(unsigned long timestepStart, unsigned long timestepEnd, unsigned beamletCount, unsigned beamletIndex)
{
	const unsigned width = timestepEnd - timestepStart;
	
	std::pair<TimeFrequencyData,TimeFrequencyMetaDataPtr> data = ReadSingleBeamlet(timestepStart*(unsigned long) 256, timestepEnd*(unsigned long) 256, beamletCount, beamletIndex);

	TimeFrequencyData allX = data.first.Make(Polarization::XX);
	TimeFrequencyData allY = data.first.Make(Polarization::YY);
	Image2DCPtr xr = allX.GetRealPart();
	Image2DCPtr xi = allX.GetImaginaryPart();
	Image2DCPtr yr = allY.GetRealPart();
	Image2DCPtr yi = allY.GetImaginaryPart();
	Mask2DCPtr mask = data.first.GetSingleMask();
	
	Image2DPtr
		outXR = Image2D::CreateUnsetImagePtr(width, 256),
		outXI = Image2D::CreateUnsetImagePtr(width, 256),
		outYR = Image2D::CreateUnsetImagePtr(width, 256),
		outYI = Image2D::CreateUnsetImagePtr(width, 256);
	Mask2DPtr
		outMask = Mask2D::CreateUnsetMaskPtr(width, 256);
	
	std::vector<double> observationTimes;
	for(unsigned long timestep = 0;timestep < timestepEnd-timestepStart;++timestep)
	{
		unsigned long timestepIndex = timestep * 256;
		SampleRow
			realX = SampleRow::MakeFromRow(xr.get(), timestepIndex, 256, 0),
			imaginaryX = SampleRow::MakeFromRow(xi.get(), timestepIndex, 256, 0),
			realY = SampleRow::MakeFromRow(yr.get(), timestepIndex, 256, 0),
			imaginaryY = SampleRow::MakeFromRow(yi.get(), timestepIndex, 256, 0);
		
		FFTTools::FFT(realX, imaginaryX);
		FFTTools::FFT(realY, imaginaryY);
		
		realX.SetVerticalImageValues(outXR.get(), timestep);
		imaginaryX.SetVerticalImageValues(outXI.get(), timestep);
		realY.SetVerticalImageValues(outYR.get(), timestep);
		imaginaryY.SetVerticalImageValues(outYI.get(), timestep);
		
		observationTimes.push_back(data.second->ObservationTimes()[timestepIndex + 256/2]);

		size_t validValues = 0;
		for(unsigned y=0;y<256;++y)
		{
			if(!mask->Value(timestepIndex + y, 0))
				++validValues;
		}
		for(unsigned y=0;y<256;++y)
		{
			outMask->SetValue(timestep, y , validValues == 0);
		}
	}
	
	data.first = TimeFrequencyData(Polarization::XX, outXR, outXI, Polarization::YY, outYR, outYI);
	data.first.SetGlobalMask(outMask);
	BandInfo band = data.second->Band();
	band.channels.clear();
	for(unsigned i=0;i<256;++i)
	{
		ChannelInfo channel;
		channel.frequencyHz = i+1;
		channel.frequencyIndex = i;
		band.channels.push_back(channel);
	}
	data.second->SetBand(band);
	data.second->SetObservationTimes(observationTimes);
	return data;
}

std::pair<TimeFrequencyData,TimeFrequencyMetaDataPtr> RSPReader::ReadAllBeamlets(unsigned long timestepStart, unsigned long timestepEnd, unsigned beamletCount)
{
	const unsigned width = timestepEnd - timestepStart;
	Image2DPtr realX = Image2D::CreateZeroImagePtr(width, beamletCount);
	Image2DPtr imaginaryX = Image2D::CreateZeroImagePtr(width, beamletCount);
	Image2DPtr realY = Image2D::CreateZeroImagePtr(width, beamletCount);
	Image2DPtr imaginaryY = Image2D::CreateZeroImagePtr(width, beamletCount);
	Mask2DPtr mask = Mask2D::CreateSetMaskPtr<true>(width, beamletCount);
	
	std::ifstream file(_rawFile.c_str(), std::ios_base::binary | std::ios_base::in);
	size_t frame = 0;
	std::set<short> stations;
	
	TimeFrequencyMetaDataPtr metaData = TimeFrequencyMetaDataPtr(new TimeFrequencyMetaData());
	BandInfo band;
	for(size_t i=0;i<beamletCount;++i)
	{
		ChannelInfo channel;
		channel.frequencyHz = i+1;
		channel.frequencyIndex = i;
		band.channels.push_back(channel);
	}
	metaData->SetBand(band);
	
	std::vector<double> observationTimes;
	
	// Read a header and determine the reading start position
	// Because timestepStart might fall within a block, the 
	RCPApplicationHeader firstHeader;
	firstHeader.Read(file);
	const unsigned long bytesPerFrame = beamletCount * firstHeader.nofBlocks * RCPBeamletData::SIZE + RCPApplicationHeader::SIZE;
	const unsigned long startFrame = timestepStart / (unsigned long) firstHeader.nofBlocks;
	const unsigned long startByte = startFrame * bytesPerFrame;
	const unsigned long offsetFromStart = timestepStart - (startFrame * firstHeader.nofBlocks);
	//Logger::Debug << "Seeking to " << startByte << " (timestepStart=" << timestepStart << ", offsetFromStart=" << offsetFromStart << ", startFrame=" << startFrame << ",bytesPerFrame=" << bytesPerFrame << ")\n";
	file.seekg(startByte, std::ios_base::beg);
	
	// Read the frames
	unsigned long x=0;
	while(x < width + offsetFromStart && file.good()) {
		RCPApplicationHeader header;
		header.Read(file);
		if(header.versionId != 2)
		{
			std::stringstream s;
			s << "Corrupted header found in frame " << frame << "!";
			throw std::runtime_error(s.str());
		}
		if(stations.count(header.stationId)==0)
		{
			stations.insert(header.stationId);
			AntennaInfo antenna;
			std::stringstream s;
			s << "LOFAR station with index " << header.stationId;
			antenna.name = s.str();
			metaData->SetAntenna1(antenna);
			metaData->SetAntenna2(antenna);
		}
		for(size_t j=0;j<beamletCount;++j)
		{
			for(size_t i=0;i<header.nofBlocks;++i)
			{
				RCPBeamletData data;
				data.Read(file);
				if(i + x < width + offsetFromStart && i + x >= offsetFromStart)
				{
					const unsigned long pos = i + x - offsetFromStart;
					realX->SetValue(pos, j, data.xr);
					imaginaryX->SetValue(pos, j, data.xi);
					realY->SetValue(pos, j, data.yr);
					imaginaryY->SetValue(pos, j, data.yi);
					mask->SetValue(pos, j, false);
				}
			}
		}
		x += header.nofBlocks;
		++frame;
	}
	//Logger::Debug << "Read " << frame << " frames.\n";
	
	for(unsigned long i=0;i<width;++i)
	{
		const unsigned long pos = i + timestepStart;
		const double time =
			(double) pos * (double) STATION_INTEGRATION_STEPS / (double) _clockSpeed;
		observationTimes.push_back(time);
	}
	
	metaData->SetObservationTimes(observationTimes);
	
	std::pair<TimeFrequencyData,TimeFrequencyMetaDataPtr> data;
	data.first = TimeFrequencyData(Polarization::XX, realX, imaginaryX, Polarization::YY, realY, imaginaryY);
	data.first.SetGlobalMask(mask);
	data.second = metaData;
	return data;
}