Task #1892: Added a rudimentary zenith imager (not working yet)

.gitattributes

@@ -317,6 +317,7 @@ CEP/DP3/AOFlagger/include/AOFlagger/imaging/fourproductcorrelatortester.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/imaging/model.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/imaging/observatorium.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/imaging/uvimager.h -text
+CEP/DP3/AOFlagger/include/AOFlagger/imaging/zenithimager.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/msio/arraycolumniterator.h -text
 CEP/DP3/AOFlagger/include/AOFlagger/msio/baselinematrixloader.h -text
@@ -592,6 +593,8 @@ CEP/DP3/AOFlagger/src/gui/tfstatoptionwindow.cpp -text
 CEP/DP3/AOFlagger/src/imaging/fourproductcorrelatortester.cpp -text
 CEP/DP3/AOFlagger/src/imaging/model.cpp -text
 CEP/DP3/AOFlagger/src/imaging/uvimager.cpp -text
+CEP/DP3/AOFlagger/src/imaging/zenithimager.cpp -text
+CEP/DP3/AOFlagger/src/imgzenith.cpp -text
 CEP/DP3/AOFlagger/src/ms2uv.cpp -text
 CEP/DP3/AOFlagger/src/msinfo.cpp -text
 CEP/DP3/AOFlagger/src/msio/baselinematrixloader.cpp -text

CEP/DP3/AOFlagger/include/AOFlagger/imaging/zenithimager.h

+/***************************************************************************
+ *   Copyright (C) 2012 by A.R. Offringa                                   *
+ *   offringa@astro.rug.nl                                                 *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#ifndef ZENITH_IMAGER_H
+#define ZENITH_IMAGER_H
+
+#include <complex>
+
+#include <AOFlagger/msio/image2d.h>
+
+class ZenithImager
+{
+	public:
+		ZenithImager();
+		
+		void Initialize(size_t resolution);
+		
+		void Clear();
+		
+		void Add(const class BandInfo &band, const std::complex<float> *samples, const bool *isRFI, double u, double v, double w, double phaseRotation);
+		
+		void FourierTransform(Image2DPtr &real, Image2DPtr &imaginary);
+		
+		Image2DCPtr UVReal() const { return _real; }
+		Image2DCPtr UVImaginary() const { return _imaginary; }
+	private:
+		void add(const class BandInfo &band, double r, double i, double u, double v, double w, double phaseRotation, double wavelength);
+		
+		static long double Wavelength(long double frequency)
+		{
+			return 299792458.0L / frequency;
+		}
+		
+		size_t _resolution, _totalCount, _outsideCount;
+		
+		Image2DPtr _real, _imaginary, _weights;
+};
+
+#endif

CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h

@@ -41,7 +41,7 @@ struct EarthPosition : public Serializable {
 		s.width(16);
 		s.precision(16);
 		s << x << "," << y << "," << z << " (alt " << sqrtl(x*x+y*y+z*z) << "), or "
-		<< "N" << Lattitude()*180/M_PI << " E" << Longitude()*180/M_PI;
+		<< "N" << Latitude()*180/M_PI << " E" << Longitude()*180/M_PI;
 		return s.str();
 	}
 	EarthPosition FromITRS(long double x, long double y, long double z);
@@ -50,8 +50,18 @@ struct EarthPosition : public Serializable {
 	{
 		return atan2l(y, x);
 	}
+	
+	long double LongitudeL() const
+	{
+		return atan2l(y, x);
+	}
 
-	double Lattitude() const
+	double Latitude() const
+	{
+		return atan2l(z, sqrtl((long double) x*x + y*y));
+	}
+	
+	long double LatitudeL() const
 	{
 		return atan2l(z, sqrtl((long double) x*x + y*y));
 	}
@@ -60,6 +70,11 @@ struct EarthPosition : public Serializable {
 	{
 		return sqrtl((long double) x*x+y*y+z*z);
 	}
+	
+	double AltitudeL() const
+	{
+		return sqrtl((long double) x*x+y*y+z*z);
+	}
 
 	virtual void Serialize(std::ostream &stream) const
 	{

CEP/DP3/AOFlagger/include/AOFlagger/msio/measurementset.h

@@ -126,11 +126,12 @@ class MeasurementSet {
 		}
 		size_t GetPolarizationCount();
 		static size_t GetPolarizationCount(const std::string &filename);
+		static struct BandInfo GetBandInfo(const std::string &filename, unsigned bandIndex);
 		size_t AntennaCount();
 		size_t FieldCount();
 		size_t BandCount();
 		struct AntennaInfo GetAntennaInfo(unsigned antennaId);
-		struct BandInfo GetBandInfo(unsigned bandIndex);
+		struct BandInfo GetBandInfo(unsigned bandIndex) {return GetBandInfo(_location, bandIndex);}
 		struct FieldInfo GetFieldInfo(unsigned fieldIndex);
 		void DataMerge(const MeasurementSet &source);
 		std::string Location() const throw() { return _location; }

CEP/DP3/AOFlagger/src/CMakeLists.txt

@@ -248,6 +248,7 @@ lofar_add_bin_program(msinfo msinfo.cpp)
 lofar_add_bin_program(ns2bbs ns2bbs.cpp)
 lofar_add_bin_program(colormapper colormapper.cpp)
 lofar_add_bin_program(versionaoflagger versionaoflagger.cc Package__Version.cc)
+lofar_add_bin_program(imgzenith imgzenith.cpp imaging/zenithimager.cpp)
 
 if(BOOST_ASIO_H_FOUND AND SIGCXX_FOUND)
   lofar_add_bin_program(aoquality aoquality.cpp ${REMOTE_FILES})

CEP/DP3/AOFlagger/src/imaging/zenithimager.cpp

+/***************************************************************************
+ *   Copyright (C) 2012 by A.R. Offringa                                   *
+ *   offringa@astro.rug.nl                                                 *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#include <AOFlagger/imaging/zenithimager.h>
+
+#include <AOFlagger/msio/antennainfo.h>
+
+#include <AOFlagger/util/ffttools.h>
+
+ZenithImager::ZenithImager() : _resolution(0), _totalCount(0), _outsideCount(0)
+{
+}
+
+void ZenithImager::Initialize(size_t resolution)
+{
+	_resolution = resolution;
+	_real = Image2D::CreateZeroImagePtr(resolution, resolution);
+	_imaginary = Image2D::CreateZeroImagePtr(resolution, resolution);
+	_weights = Image2D::CreateZeroImagePtr(resolution, resolution);
+}
+
+void ZenithImager::Clear()
+{
+	_real->SetAll(0.0);
+	_imaginary->SetAll(0.0);
+	_weights->SetAll(0.0);
+	_totalCount = 0;
+	_outsideCount = 0;
+}
+
+void ZenithImager::Add(const BandInfo &band, const std::complex<float> *samples, const bool *isRFI, double u, double v, double w, double phaseRotation)
+{
+	size_t n = band.channelCount;
+	for(size_t f=0;f<n;++f)
+	{
+		if(!isRFI[f])
+		{
+			const double wavelength = Wavelength(band.channels[f].frequencyHz);
+			const double r = samples[f].real(), i = samples[f].imag();
+			add(band, r, i, u, v, w, phaseRotation, wavelength);
+		}
+	}
+}
+
+void ZenithImager::add(const BandInfo &band, double r, double i, double u, double v, double w, double phaseRotation, double wavelength)
+{
+	const double factor = 1.0 / wavelength;
+	u *= factor;
+	v *= factor;
+	//w *= factor;
+	phaseRotation *= factor * 2 * M_PI;
+	
+	const double centre = (double) (_resolution + 1.0) * 0.5;
+	
+	int uPos = (int) round(u + centre);
+	int vPos = (int) round(v + centre);
+	
+	if(uPos >= 0 && vPos >= 0 && uPos < (int) _resolution && vPos < (int) _resolution)
+	{
+		const double sinR = sin(phaseRotation), cosR = cos(phaseRotation);
+		const double rotatedR = r * cosR - i * sinR;
+		const double rotatedI = r * sinR + i * cosR;
+		
+		_real->AddValue(uPos, vPos, rotatedR);
+		_imaginary->AddValue(uPos, vPos, rotatedI);
+		_weights->AddValue(uPos, vPos, 1.0);
+		
+		int uPos2 = (int) round(centre - u);
+		int vPos2 = (int) round(centre - v);
+		if(uPos2 >= 0 && vPos2 >= 0 && uPos2 < (int) _resolution && vPos2 < (int) _resolution)
+		{
+			_real->AddValue(uPos2, vPos2, rotatedR);
+			_imaginary->AddValue(uPos2, vPos2, -rotatedI);
+			_weights->AddValue(uPos2, vPos2, 1.0);
+		}
+	} else {
+		_outsideCount++;
+	}
+	
+	++_totalCount;
+}
+
+void ZenithImager::FourierTransform(Image2DPtr &real, Image2DPtr &imaginary)
+{
+	std::cout << "Performing FT of " << _totalCount << " samples, " << _outsideCount << " fell outside uv area... " << std::flush;
+	
+	double normFactor = _weights->Sum() / ((num_t) _real->Height() * _real->Width());
+	for(size_t y=0;y<_real->Height();++y) {
+		for(size_t x=0;x<_real->Width();++x) {
+			num_t weight = _weights->Value(x, y);
+			if(weight != 0.0)
+			{
+				_real->SetValue(x, y, _real->Value(x, y) * normFactor / weight);
+				_imaginary->SetValue(x, y, _imaginary->Value(x, y) * normFactor / weight);
+				_weights->SetValue(x, y, 1.0);
+			} 
+		}
+	}
+	
+	real = Image2D::CreateUnsetImagePtr(_resolution, _resolution);
+	imaginary = Image2D::CreateUnsetImagePtr(_resolution, _resolution);
+	FFTTools::CreateFFTImage(*_real, *_imaginary, *real, *imaginary);
+}

CEP/DP3/AOFlagger/src/imgzenith.cpp

+/***************************************************************************
+ *   Copyright (C) 2011 by A.R. Offringa                                   *
+ *   offringa@astro.rug.nl                                                 *
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ *   This program is distributed in the hope that it will be useful,       *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
+ *   GNU General Public License for more details.                          *
+ *                                                                         *
+ *   You should have received a copy of the GNU General Public License     *
+ *   along with this program; if not, write to the                         *
+ *   Free Software Foundation, Inc.,                                       *
+ *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
+ ***************************************************************************/
+
+#include <tables/Tables/Table.h>
+
+#include <measures/Measures/UVWMachine.h>
+#include <measures/Measures/MEpoch.h>
+
+#include <ms/MeasurementSets/MeasurementSet.h>
+
+#include <AOFlagger/msio/antennainfo.h>
+#include <AOFlagger/msio/measurementset.h>
+
+#include <AOFlagger/imaging/zenithimager.h>
+#include <AOFlagger/msio/pngfile.h>
+
+const casa::Unit radUnit("rad");
+const casa::Unit dayUnit("d");
+
+void repoint(casa::MDirection &phaseDirection, casa::MPosition &position, const casa::MEpoch &obstime, double &u, double &v, double &w, double &phaseRotation)
+{
+	//MPosition location(MVPosition(Quantity(1, "km"), Quantity(150, "deg"), Quantity(20, "deg")), MPosition::WGS84);
+	
+	//casa::MEpoch obstime(casa::Quantity(t, dayUnit), casa::MEpoch::UTC);
+	//std::cout << "Time=" << obstime.getValue() << '\n';
+	
+	casa::MeasFrame timeAndLocation(obstime, position);
+	
+	casa::MDirection::Ref refApparent(casa::MDirection::APP, timeAndLocation);
+	
+	// Calculate zenith
+	casa::MDirection outDirection(casa::Quantity(M_PI, radUnit), casa::Quantity(0.0, radUnit), refApparent);
+	//std::cout << "Out=" << outDirection.getValue() << '\n';
+	
+	// Construct a CASA UVW converter
+	casa::UVWMachine uvwConverter(outDirection, phaseDirection);
+	casa::Vector<double> uvwVector(3);
+	uvwVector[0] = u;
+	uvwVector[1] = v;
+	uvwVector[2] = w;
+	//std::cout << "In: " << u << ',' << v << ',' << w << '\n';
+	//std::cout << "Phase shift: " << uvwConverter.getPhase(uvwVector) << '\n';
+	//uvwConverter.convertUVW(uvwVector); // getPhase already does that!
+	u = uvwVector[0];
+	v = uvwVector[1];
+	w = uvwVector[2];
+	//std::cout << "Out: " << u << ',' << v << ',' << w << '\n';
+	//std::cout << "Phase centre: " << uvwConverter.phaseCenter().getValue() << '\n';
+}
+
+int main(int argc, char *argv[])
+{
+	std::string filename(argv[1]);
+	size_t integrationSteps;
+	if(argc >= 3)
+		integrationSteps = atoi(argv[2]);
+	else
+		integrationSteps = 60;
+	
+	std::cout << "Opening " << filename << "...\n";
+	
+	const unsigned polarizationCount = MeasurementSet::GetPolarizationCount(filename);
+	const BandInfo band = MeasurementSet::GetBandInfo(filename, 0);
+	
+	const unsigned antennaIndex = 0;
+	
+	casa::MeasurementSet table(argv[1]);
+	casa::MEpoch::ROScalarColumn timeColumn(table, "TIME");
+	casa::ROArrayColumn<double> uvwColumn(table, "UVW");
+	casa::ROScalarColumn<int> ant1Column(table, "ANTENNA1");
+	casa::ROScalarColumn<int> ant2Column(table, "ANTENNA2");
+	casa::ROArrayColumn<casa::Complex> dataColumn(table, "DATA");
+	casa::ROArrayColumn<bool> flagColumn(table, "FLAG");
+	
+	casa::Table antennaTable(table.antenna());
+	casa::MPosition::ROScalarColumn antPositionColumn(antennaTable, "POSITION");
+	casa::ROScalarColumn<casa::String> antNameColumn(antennaTable, "NAME");
+	casa::MPosition position  = antPositionColumn(antennaIndex);
+	std::cout << "Imaging zenith of antenna " << antNameColumn(antennaIndex)
+		<< ", pos=" << position.getValue() << '\n';
+		
+	casa::Table fieldTable(table.field());
+	casa::MDirection::ROArrayColumn phaseDirColumn(fieldTable, "PHASE_DIR");
+	casa::MDirection phaseDirection = *phaseDirColumn(0).begin();
+	std::cout << "Phase direction: " << phaseDirection.getValue() << '\n';
+
+	std::complex<float> *samples[polarizationCount];
+	bool *isRFI[polarizationCount];
+	for(unsigned p = 0; p < polarizationCount; ++p)
+	{
+		isRFI[p] = new bool[band.channelCount];
+		samples[p] = new std::complex<float>[band.channelCount];
+	}
+
+	unsigned row = 0;
+	ZenithImager imager;
+	imager.Initialize(2048);
+	size_t timeStep = 0;
+	while(row<table.nrow())
+	{
+		const casa::MEpoch t = timeColumn(row);
+		do
+		{
+			if(ant1Column(row) != ant2Column(row))
+			{
+				casa::Array<double> uvw = uvwColumn(row);
+				casa::Array<double>::const_iterator uvw_i = uvw.begin();
+				double u = *uvw_i; ++uvw_i;
+				double v = *uvw_i; ++uvw_i;
+				double w = *uvw_i;
+				double phaseRotation;
+				repoint(phaseDirection, position, t, u, v, w, phaseRotation);
+				
+				const casa::Array<casa::Complex> dataArray = dataColumn(row);
+				const casa::Array<bool> flagArray = flagColumn(row);
+				
+				casa::Array<casa::Complex>::const_iterator dataIter = dataArray.begin();
+				casa::Array<bool>::const_iterator flagIter = flagArray.begin();
+				
+				for(unsigned channel = 0; channel<band.channelCount; ++channel)
+				{
+					for(unsigned p = 0; p < polarizationCount; ++p)
+					{
+						samples[p][channel] = *dataIter;
+						isRFI[p][channel] = *flagIter;
+						
+						++dataIter;
+						++flagIter;
+					}
+				}
+				
+				imager.Add(band, samples[0], isRFI[0], u, v, w, phaseRotation);
+			}
+			++row;
+		} while(row<table.nrow() && timeColumn(row).getValue() == t.getValue());
+		
+		timeStep++;
+		if(timeStep % integrationSteps == 0)
+		{
+			Image2DPtr real, imaginary;
+			imager.FourierTransform(real, imaginary);
+			
+			std::stringstream s;
+			if(timeStep < 10000) s << '0';
+			if(timeStep < 1000) s << '0';
+			if(timeStep < 100) s << '0';
+			if(timeStep < 10) s << '0';
+			s << timeStep << ".png";
+			std::cout << "Saving " << s.str() << "... " << std::flush;
+			PngFile::Save(*real, std::string("zen/") + s.str());
+			PngFile::Save(*imager.UVReal(), std::string("zen-uv/") + s.str());
+			imager.Clear();
+			std::cout << "Done.\n";
+		}
+	}
+}

CEP/DP3/AOFlagger/src/msinfo.cpp

@@ -90,7 +90,7 @@ int main(int argc, char *argv[])
 			cout.width(16);
 			cout.precision(16);
 			const AntennaInfo antenna = set.GetAntennaInfo(i);
-			cout << antenna.id << '\t' << antenna.position.Lattitude()*180.0/M_PI << '\t' << antenna.position.Longitude()*180.0/M_PI << '\t' << antenna.position.Altitude() << '\n';
+			cout << antenna.id << '\t' << antenna.position.Latitude()*180.0/M_PI << '\t' << antenna.position.Longitude()*180.0/M_PI << '\t' << antenna.position.Altitude() << '\n';
 		}
 	}
 	else {

CEP/DP3/AOFlagger/src/msio/measurementset.cpp

@@ -198,10 +198,10 @@ struct AntennaInfo MeasurementSet::GetAntennaInfo(unsigned antennaId)
 	return info;
 }
 
-struct BandInfo MeasurementSet::GetBandInfo(unsigned bandIndex)
+BandInfo MeasurementSet::GetBandInfo(const std::string &filename, unsigned bandIndex)
 {
 	BandInfo band;
-	casa::MeasurementSet ms(_location);
+	casa::MeasurementSet ms(filename);
 	casa::Table spectralWindowTable = ms.spectralWindow();
 	casa::ROScalarColumn<int> numChanCol(spectralWindowTable, "NUM_CHAN");
 	casa::ROArrayColumn<double> frequencyCol(spectralWindowTable, "CHAN_FREQ");