From 124e3302028f0b0f6cbfd56a011d92c67d04564d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Offringa?= <offringa@astron.nl>
Date: Wed, 14 Mar 2012 10:41:00 +0000
Subject: [PATCH] Task #1892: Added a rudimentary zenith imager (not working
yet)
---
.gitattributes | 3 +
.../include/AOFlagger/imaging/zenithimager.h | 56 ++++++
.../include/AOFlagger/msio/antennainfo.h | 19 +-
.../include/AOFlagger/msio/measurementset.h | 3 +-
CEP/DP3/AOFlagger/src/CMakeLists.txt | 1 +
.../AOFlagger/src/imaging/zenithimager.cpp | 120 ++++++++++++
CEP/DP3/AOFlagger/src/imgzenith.cpp | 173 ++++++++++++++++++
CEP/DP3/AOFlagger/src/msinfo.cpp | 2 +-
CEP/DP3/AOFlagger/src/msio/measurementset.cpp | 4 +-
9 files changed, 375 insertions(+), 6 deletions(-)
create mode 100644 CEP/DP3/AOFlagger/include/AOFlagger/imaging/zenithimager.h
create mode 100644 CEP/DP3/AOFlagger/src/imaging/zenithimager.cpp
create mode 100644 CEP/DP3/AOFlagger/src/imgzenith.cpp
diff --git a/.gitattributes b/.gitattributes
index 73df04b9d8b..ca09c7cf98c 100644
--- a/.gitattributes
+++ b/.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
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/imaging/zenithimager.h b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/zenithimager.h
new file mode 100644
index 00000000000..d42e59bfa02
--- /dev/null
+++ b/CEP/DP3/AOFlagger/include/AOFlagger/imaging/zenithimager.h
@@ -0,0 +1,56 @@
+/***************************************************************************
+ * 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
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h b/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h
index 9baf25ee4c2..c021bc987d9 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/msio/antennainfo.h
+++ b/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
{
diff --git a/CEP/DP3/AOFlagger/include/AOFlagger/msio/measurementset.h b/CEP/DP3/AOFlagger/include/AOFlagger/msio/measurementset.h
index 11d4c3a9881..9921e349e86 100644
--- a/CEP/DP3/AOFlagger/include/AOFlagger/msio/measurementset.h
+++ b/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; }
diff --git a/CEP/DP3/AOFlagger/src/CMakeLists.txt b/CEP/DP3/AOFlagger/src/CMakeLists.txt
index 74c6960061c..c85b2fd76b1 100644
--- a/CEP/DP3/AOFlagger/src/CMakeLists.txt
+++ b/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})
diff --git a/CEP/DP3/AOFlagger/src/imaging/zenithimager.cpp b/CEP/DP3/AOFlagger/src/imaging/zenithimager.cpp
new file mode 100644
index 00000000000..7ca135d2530
--- /dev/null
+++ b/CEP/DP3/AOFlagger/src/imaging/zenithimager.cpp
@@ -0,0 +1,120 @@
+/***************************************************************************
+ * 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);
+}
diff --git a/CEP/DP3/AOFlagger/src/imgzenith.cpp b/CEP/DP3/AOFlagger/src/imgzenith.cpp
new file mode 100644
index 00000000000..328018d4ca8
--- /dev/null
+++ b/CEP/DP3/AOFlagger/src/imgzenith.cpp
@@ -0,0 +1,173 @@
+/***************************************************************************
+ * 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";
+ }
+ }
+}
diff --git a/CEP/DP3/AOFlagger/src/msinfo.cpp b/CEP/DP3/AOFlagger/src/msinfo.cpp
index b3160e59599..48f146d74d3 100644
--- a/CEP/DP3/AOFlagger/src/msinfo.cpp
+++ b/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 {
diff --git a/CEP/DP3/AOFlagger/src/msio/measurementset.cpp b/CEP/DP3/AOFlagger/src/msio/measurementset.cpp
index 747867dbeb1..88ccfaf3f80 100644
--- a/CEP/DP3/AOFlagger/src/msio/measurementset.cpp
+++ b/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");
--
GitLab