From a0c4c343e3c48d7970843a2ce1d03d9153381f29 Mon Sep 17 00:00:00 2001
From: Adriaan Renting <renting@astron.nl>
Date: Thu, 15 Feb 2007 16:24:02 +0000
Subject: [PATCH] BugId: 1028 Made some changes so it actually compiles now.
Needs some more work to actually produce the right results.
---
.../CS1_DFTImager/ImagerProcessControl.h | 18 +-
Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc | 705 ++++++-----------
Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc~ | 706 ++++++------------
Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h | 51 +-
Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h~ | 51 +-
Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc | 101 ++-
Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc~ | 103 ++-
Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h | 84 ++-
Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h~ | 86 ++-
.../CS1_DFTImager/src/ImagerProcessControl.cc | 21 +-
.../src/ImagerProcessControl.cc~ | 25 +-
Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc | 82 +-
Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc~ | 82 +-
Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h | 11 +-
Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h~ | 11 +-
Appl/CEP/CS1/CS1_DFTImager/src/Makefile.am | 2 +
Appl/CEP/CS1/CS1_DFTImager/src/Makefile.in | 6 +-
17 files changed, 856 insertions(+), 1289 deletions(-)
diff --git a/Appl/CEP/CS1/CS1_DFTImager/include/CS1_DFTImager/ImagerProcessControl.h b/Appl/CEP/CS1/CS1_DFTImager/include/CS1_DFTImager/ImagerProcessControl.h
index 515b29dde28..12a139157e0 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/include/CS1_DFTImager/ImagerProcessControl.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/include/CS1_DFTImager/ImagerProcessControl.h
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by Adriaan Renting *
+ * Copyright (C) 2007 by Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -17,8 +17,8 @@
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
-#ifndef LOFARFLAGGERPROCESSCONTROL_H
-#define LOFARFLAGGERPROCESSCONTROL_H
+#ifndef IMAGER_PROCESSCONTROL_H
+#define IMAGER_PROCESSCONTROL_H
#include <PLC/ProcessControl.h>
#include <APS/ParameterSet.h>
@@ -27,7 +27,7 @@
@author Adriaan Renting
*/
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
@@ -40,13 +40,13 @@ namespace LOFAR
std::string itsMSName;
std::string itsImageName;
int itsResolution;
- MS_File* intMS;
- Image_File* intImage;
+ MS_File* itsMS;
+ Image_File* itsImage;
DFTImager* itsImager;
public:
- FlaggerProcessControl(void);
+ ImagerProcessControl(void);
- ~FlaggerProcessControl(void);
+ ~ImagerProcessControl(void);
// \name Command to control the processes.
// There are a dozen commands that can be sent to a application process
// to control its flow. The return values for these command are:<br>
@@ -78,7 +78,7 @@ namespace LOFAR
tribool snapshot(const std::string&);
std::string askInfo(const std::string&);
- }; //class FlaggerProcessControl
+ }; //class ImagerProcessControl
} //namespace CS1
}; //namespace LOFAR
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc
index 5d7adeb3f92..b0c92c6db62 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -21,10 +21,6 @@
#include <tables/Tables/TableIter.h>
#include "DFTImager.h"
#include <casa/Quanta/MVEpoch.h>
-#include <images/Images/PagedImage.h>
-#include <images/Images/ImageInfo.h>
-#include <images/Images/SubImage.h>
-#include <images/Images/ImageUtilities.h>
using namespace casa;
@@ -32,7 +28,7 @@ enum CorrelationTypes {None=0,I=1,Q=2,U=3,V=4,RR=5,RL=6,LR=7,LL=8,XX=9,XY=10,YX=
//===============>>> itoa <<<===============
//ANSI C++ doesn't seem to have a decent function for this, or I'm not aware of it. Need to rename it to IntToStr(), to avoid confusion
-std::string itoa(int value, int base=10)
+std::string itoa(int value, int base=10)
{ // found on http://www.jb.man.ac.uk/~slowe/cpp/itoa.html
//maybe copyright Robert Jan Schaper, Ray-Yuan Sheu, Rodrigo de Salvo Braz, Wes Garland and John Maloney
enum { kMaxDigits = 35 };
@@ -52,515 +48,248 @@ std::string itoa(int value, int base=10)
return buf;
}
-//===============>>> ComplexMedianFlagger::ComplexMedianFlagger <<<===============
-/* initialize some meta data and get the datastorage the right size. */
-ComplexMedianFlagger::ComplexMedianFlagger(MS_File* InputMSfile,
- int InputWindowSize,
- bool UseOnlyXpolarizations,
- double InputMinThreshold,
- double InputMaxThreshold)
+namespace LOFAR
{
- MSfile = InputMSfile;
- WindowSize = InputWindowSize;
- MinThreshold = InputMinThreshold;
- MaxThreshold = InputMaxThreshold;
- NumAntennae = (*MSfile).itsNumAntennae;
- NumPairs = (*MSfile).itsNumPairs;
- NumBands = (*MSfile).itsNumBands;
- NumChannels = (*MSfile).itsNumChannels;
- NumPolarizations = (*MSfile).itsNumPolarizations;
- NoiseLevel = (*MSfile).itsNoiseLevel;
- NumTimeslots = (*MSfile).itsNumTimeslots;
- AntennaNames = (*MSfile).itsAntennaNames;
-
- PairsIndex.resize(NumPairs);
- Statistics = Cube<int>(NumBands, NumAntennae, NumAntennae, 0);
- PolarizationsToCheck.resize(NumPolarizations);
- DeterminePolarizationsToCheck(UseOnlyXpolarizations);
-
- TimeslotData.resize(NumPairs*NumBands);
- for (int i = 0; i < NumPairs*NumBands; i++)
- { TimeslotData[i].resize(NumPolarizations, NumChannels, WindowSize);
- }
-
- int index = 0;
- for (int i = 0; i < NumAntennae; i++)
- { for(int j = i; j < NumAntennae; j++)
- { PairsIndex[index] = pairii(i, j);
- BaselineIndex[pairii(i, j)] = index++;
- }
- }
- ComputeBaselineLengths();
-}
-
-//===============>>> ComplexMedianFlagger::~ComplexMedianFlagger <<<===============
-
-ComplexMedianFlagger::~ComplexMedianFlagger()
-{
-}
-
-//===============>>> ComplexMedianFlagger::DetermineCorrelationsToCheck <<<===============
-/* create a list of polarizations we want to check, maybe we only want to to XY, YX */
-void ComplexMedianFlagger::DeterminePolarizationsToCheck(bool UseOnlyXpolarizations)
-{
- if (UseOnlyXpolarizations)
+ namespace CS1
{
- bool noCorrError = true;
- for (int i = 0; i < NumPolarizations; i++)
- {
- switch((*MSfile).itsPolarizations[i])
- {
- case None:
- case I:
- case RR:
- case LL:
- case XX:
- case YY:
- if (UseOnlyXpolarizations)
- PolarizationsToCheck[i] = false;
- break;
- case Q:
- case U:
- case V:
- case RL:
- case LR:
- case XY:
- case YX:
- noCorrError = false;
- if (UseOnlyXpolarizations)
- PolarizationsToCheck[i] = true;
- break;
- }
- }
- if (noCorrError)
+ //===============>>> DFTImager::DFTImager <<<===============
+ /* initialize some meta data and get the datastorage the right size. */
+ DFTImager::DFTImager(MS_File* InputMSfile,
+ Image_File*InputImageFile)
{
- cout << "There are no crosspolarizations to flag!";
- exit(1);
- }
- }
- else
- {
- for (int i = 0; i < NumPolarizations; i++)
- { PolarizationsToCheck[i] = true;
- }
- }
-}
+ MSfile = InputMSfile;
+ Imagefile = InputImageFile;
+ NumAntennae = (*MSfile).itsNumAntennae;
+ NumPairs = (*MSfile).itsNumPairs;
+ NumBands = (*MSfile).itsNumBands;
+ NumChannels = (*MSfile).itsNumChannels;
+ NumPolarizations = (*MSfile).itsNumPolarizations;
+ NumTimeslots = (*MSfile).itsNumTimeslots;
+ AntennaNames = (*MSfile).itsAntennaNames;
+ Resolution = 41;
-//===============>>> ComplexMedianFlagger::ComputeBaselineLengths <<<===============
-/* compute baseline lengths, and determine the longest one.*/
-void ComplexMedianFlagger::ComputeBaselineLengths()
-{
- MaxBaselineLength = 0.0;
- BaselineLengths.resize(NumPairs);
- //Antenna positions
- MSAntenna antenna = (*MSfile).antenna();
- ROArrayColumn<Double> position(antenna, "POSITION");
- for (int i = 0; i < NumAntennae; i++ )
- {
- for (int j = i; j < NumAntennae; j++)
- {
- Vector<Double> p(position(i) - position(j));
- double temp = sqrt(p(0)*p(0) + p(1)*p(1) + p(2)*p(2));
- BaselineLengths[BaselineIndex[pairii(i, j)]] = temp;
- if (temp > MaxBaselineLength && temp < 3000000) //radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243
- { MaxBaselineLength = temp; // non-existent antenna's can have position (0,0,0)
- }
- }
- }
-}
+ PairsIndex.resize(NumPairs);
+ // Statistics = Cube<int>(NumBands, NumAntennae, NumAntennae, 0);
+ // PolarizationsToCheck.resize(NumPolarizations);
+ // DeterminePolarizationsToCheck(UseOnlyXpolarizations);
-//===============>>> ComplexMedianFlagger::FlagDataOrBaselines <<<===============
-/*This function outputs the gathered statistics.*/
-void ComplexMedianFlagger::ProcessStatistics()
-{
- vector<int> bands(NumBands);
- vector<int> antennae(NumAntennae);
- unsigned int namelength = 6;
- for(int i = 0; i < NumAntennae; i++)
- {
- if (namelength < AntennaNames[i].size())
- { namelength = AntennaNames[i].size();
- }
- }
- for (int i = 0; i < NumBands; i++)
- {
- cout << "Band: " << i+1 << endl;
- cout << string(namelength+1,' ');
- for(int j = 0; j < NumAntennae; j++)
- {
- string out = AntennaNames[j];
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl;
- for(int j = 0; j < NumAntennae; j++)
- {
- string out = AntennaNames[j];
- out.resize(namelength+1,' ');
- cout << out;
- for(int k = 0; k < NumAntennae; k++)
- {
- if (k < j) //We print a complete array, but we have inspected only those where k >= j
- {
- int val = 100 * Statistics(i,k,j) / (NumChannels * NumTimeslots * NumPolarizations);
- bands[i] += val;
- antennae[j] += val;
- antennae[k] += val;
- string out = itoa(val) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- else
- {
- int val = 100 * Statistics(i,j,k) / (NumChannels * NumTimeslots * NumPolarizations);
- bands[i] += val;
- antennae[j] += val;
- antennae[k] += val;
- string out = itoa(val) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- /*{ cout << rms << "Faulty baseline detected: Antenna " << j+1
- << ":" << k+1 << " SpectralWindow: "<< i+1 << endl;
- }
- }*/
+ TimeslotData.resize(NumPairs*NumBands);
+ for (int i = 0; i < NumPairs*NumBands; i++)
+ { TimeslotData[i].resize(NumPolarizations, NumChannels, WindowSize);
}
- cout << endl;
- }
- }
- cout << "Bands (flagged %): ";
- for (int i = 0; i < NumBands; i++)
- {
- string out = string("IVC") + itoa(i);
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl << " ";
- for (int i = 0; i < NumBands; i++)
- {
- string out = itoa(bands[i] / (NumAntennae*NumAntennae)) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl << "Antennae (flagged %): " ;
- for(int j = 0; j < NumAntennae; j++)
- {
- string out = AntennaNames[j];
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl << " ";
- for (int i = 0; i < NumAntennae; i++)
- {
- string out = itoa(antennae[i] / (NumBands*NumAntennae*2)) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl;
-}
-//===============>>> ComplexMedianFlagger::FlagTimeslot <<<===============
-/* This function inspects each visibility in a cetain baseline-band
-and flags on complexe distance, then determines to flag the entire baseline-band
-based on the RMS of the points it didn't flag.*/
-bool ComplexMedianFlagger::FlagBaselineBand(Matrix<Bool>* Flags,
- Cube<Complex>* Timeslots,
- vector<double>* rms,
- int* flagCounter,
- double FlagThreshold,
- int Position,
- bool flagRMS)
-{
- Vector<Float> Reals(WindowSize);
- Vector<Float> Imags(WindowSize);
- vector<double> RMS(NumPolarizations);
- vector<int> RMSCounter(NumPolarizations);
- bool FlagCompleteRow = true;
- int flagcount = 0;
- for (int i = NumChannels-1; i >= 0; i--)
- {
- bool FlagAllPolarizations = false;
- for (int j = NumPolarizations-1; j >= 0; j--)
- { //we need to loop twice, once to determine FlagAllCorrelations
- if (!FlagAllPolarizations && PolarizationsToCheck[j])
- {
- for (int k = 0; k < WindowSize; k++)
- { //This might be faster in some other way ?
- Reals[k] = (*Timeslots)(j, i, k).real();
- Imags[k] = (*Timeslots)(j, i, k).imag();
+
+ int index = 0;
+ for (int i = 0; i < NumAntennae; i++)
+ { for(int j = i; j < NumAntennae; j++)
+ { PairsIndex[index] = pairii(i, j);
+ BaselineIndex[pairii(i, j)] = index++;
}
- float real = medianInPlace(Reals, false);
- float imag = medianInPlace(Imags, false);
- FlagAllPolarizations |= FlagThreshold < abs((*Timeslots)(j, i, Position)
- - Complex(real, imag));
- }
- }
- for (int j = NumPolarizations-1; j >= 0; j--)
- { //the second loop we set the flags or calculate RMS
- if (FlagAllPolarizations)
- { (*Flags)(j, i) = true;
- flagcount++;
- }
- else
- {
- FlagCompleteRow = false;
- (*Flags)(j, i) = false; //could check existing flag ?
- double value = abs((*Timeslots)(j, i, Position));
- RMS[j] += value*value;
- (RMSCounter[j])++;
}
+ ComputeBaselineLengths();
}
- }
- //these need to be separated out into a different function for clarity
- if (flagcount > 0.9 * NumChannels * NumPolarizations) //more as 90% bad
- {
- FlagCompleteRow = true;
- (*Flags) = true;
- flagcount = NumChannels * NumPolarizations;
- }
- if (!FlagCompleteRow && flagRMS) //RMSCounter > 0
- {
- for (int j = NumPolarizations-1; j >= 0; j--)
- {
- (*rms)[j] = sqrt( RMS[j] / RMSCounter[j] );
-// if (PolarizationsToCheck[j] && ((*rms)[j] > 1.5 || (*rms)[j] <= 0.5 * NoiseLevel)) //we think this RMS value is non-physical
- if (PolarizationsToCheck[j] && ((*rms)[j] > 25 || (*rms)[j] <= 0.5 * NoiseLevel)) //we think this RMS value is non-physical
- {
- FlagCompleteRow = true;
- (*Flags) = true;
- flagcount = NumChannels * NumPolarizations;
- }
- }
- }
- (*flagCounter) += flagcount;
- return FlagCompleteRow;
-}
-//===============>>> ComplexMedianFlagger::FlagBaseline <<<===============
-/* This function iterates over baseline and band and uses FlagBaselineBand() to determine
- for each one if it needs to be flagged. It treats the autocorrelations separately,
- to detect entire malfunctioning telescopes. Finally it writes the flags.
-*/
-void ComplexMedianFlagger::FlagTimeslot(TableIterator* flag_iter,
- bool flagDatapoints,
- bool flagRMS,
- bool ExistingFlags,
- Int Position)
-{
- Cube<bool> Flags(NumPolarizations, NumChannels, NumPairs*NumBands, false);
- vector<bool> FlagCompleteRow(NumPairs*NumBands);
- vector<int> stats(NumPairs*NumBands);
- for (int i = 0; i < NumBands; i++)
- {
- vector<double> rms(NumPolarizations);
- vector< vector<double> > AutoCorrRMS(NumPolarizations, NumAntennae);
- for (int j = 0; j < NumPolarizations; j++) // this needs to be written more elegantly
- {
- if ((*MSfile).itsPolarizations[j] == YX)
- { PolarizationsToCheck[j] = false; // We're not checking YX because WSRT writes 0.0 into them
- }
- }
- for(int j = 0; j < NumAntennae; j++)
- { //check auto correlations
- int index = i*NumPairs + BaselineIndex[pairii(j, j)];
- Matrix<Bool> flags = Flags.xyPlane(index);
- FlagCompleteRow[index] = FlagBaselineBand(&flags,
- &(TimeslotData[index]),
- &rms, &(stats[index]),
- 3000 * NoiseLevel,
- Position, flagRMS);
- for(int k = 0; k < NumPolarizations; k++)
- { AutoCorrRMS[k][j] = rms[k];
- }
- }
- for(int k = 0; k < NumPolarizations; k++) //check on median
+
+ //===============>>> DFTImager::~DFTImager <<<===============
+
+ DFTImager::~DFTImager()
{
- if (PolarizationsToCheck[k])
- {
- double median = medianInPlace(Vector<double>(AutoCorrRMS[k]), false);
- for(int j = 0; j < NumAntennae; j++)
- {
- int index = i*NumPairs + BaselineIndex[pairii(j, j)];
- if (AutoCorrRMS[k][j] < median/10 || AutoCorrRMS[k][j] > median*10)
- { FlagCompleteRow[index] = true;
- }
- }
- }
}
- for (int j = 0; j < NumPolarizations; j++) // this needs to be written more elegantly
+
+ //===============>>> DFTImager::DetermineCorrelationsToCheck <<<===============
+ /* create a list of polarizations we want to check, maybe we only want to to XY, YX */
+ // void DFTImager::DeterminePolarizationsToCheck(bool UseOnlyXpolarizations)
+ // {
+ // if (UseOnlyXpolarizations)
+ // {
+ // bool noCorrError = true;
+ // for (int i = 0; i < NumPolarizations; i++)
+ // {
+ // switch((*MSfile).itsPolarizations[i])
+ // {
+ // case None:
+ // case I:
+ // case RR:
+ // case LL:
+ // case XX:
+ // case YY:
+ // if (UseOnlyXpolarizations)
+ // PolarizationsToCheck[i] = false;
+ // break;
+ // case Q:
+ // case U:
+ // case V:
+ // case RL:
+ // case LR:
+ // case XY:
+ // case YX:
+ // noCorrError = false;
+ // if (UseOnlyXpolarizations)
+ // PolarizationsToCheck[i] = true;
+ // break;
+ // }
+ // }
+ // if (noCorrError)
+ // {
+ // cout << "There are no crosspolarizations to flag!";
+ // exit(1);
+ // }
+ // }
+ // else
+ // {
+ // for (int i = 0; i < NumPolarizations; i++)
+ // { PolarizationsToCheck[i] = true;
+ // }
+ // }
+ // }
+
+ //===============>>> DFTImager::ComputeBaselineLengths <<<===============
+ /* compute baseline lengths, and determine the longest one.*/
+ void DFTImager::ComputeBaselineLengths()
{
- if ((*MSfile).itsPolarizations[j] == YX)
- { PolarizationsToCheck[j] = true; // We're not checking YX
- }
- }
- // this code coud be separated into different functions
- vector< vector<double> > CrossCorrRMS(NumPolarizations, NumPairs - NumAntennae);
- for(int j = 0; j < NumAntennae; j++)
- {
- for(int k = j+1; k < NumAntennae; k++)
- {
- int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- Matrix<Bool> flags = Flags.xyPlane(index);
- if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
- { //we skip the non-existent telescopes
- FlagCompleteRow[index] = FlagBaselineBand(&flags,
- &(TimeslotData[index]),
- &rms, &(stats[index]),
- (MinThreshold + (MaxThreshold - MinThreshold)
- * BaselineLengths[BaselineIndex[pairii(j, k)]]
- / MaxBaselineLength
- ) * NoiseLevel,
- Position, flagRMS);
- for(int l = 0; l < NumPolarizations; l++)
- { CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] = rms[l]; //a very complex way of indexing, could be simplified?
- }
- if ( ( FlagCompleteRow[i*NumPairs + BaselineIndex[pairii(j, j)]]
- ||FlagCompleteRow[i*NumPairs + BaselineIndex[pairii(k, k)]])
- &&flagRMS && !FlagCompleteRow[index])
- { //We flag it if the autocorrelation is flagged
- FlagCompleteRow[index] = true;
+ MaxBaselineLength = 0.0;
+ BaselineLengths.resize(NumPairs);
+ //Antenna positions
+ MSAntenna antennas = (*MSfile).antennas();
+ ROArrayColumn<Double> position(antennas, "POSITION");
+ for (int i = 0; i < NumAntennae; i++ )
+ {
+ for (int j = i; j < NumAntennae; j++)
+ {
+ Vector<Double> p(position(i) - position(j));
+ double temp = sqrt(p(0)*p(0) + p(1)*p(1) + p(2)*p(2));
+ BaselineLengths[BaselineIndex[pairii(i, j)]] = temp;
+ if (temp > MaxBaselineLength && temp < 3000000) //radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243
+ { MaxBaselineLength = temp; // non-existent antenna's can have position (0,0,0)
}
- }
+ }
}
- }
- for(int l = 0; l < NumPolarizations; l++) // check on median
+ }
+
+ //===============>>> DFTImager::ImageTimeslot <<<===============
+ /* This function iterates over baseline and band and uses FlagBaselineBand() to determine
+ for each one if it needs to be flagged. It treats the autocorrelations separately,
+ to detect entire malfunctioning telescopes. Finally it writes the flags.
+ */
+ void DFTImager::ImageTimeslot()
{
- if (PolarizationsToCheck[l])
- {
- double median = medianInPlace(Vector<double>(CrossCorrRMS[l]), false);
+/* Cube<bool> Image(NumPolarizations, NumChannels, NumPairs*NumBands, false);
+ vector<bool> test(NumPairs*NumBands);
+ for (int i = 0; i < NumBands; i++)
+ {
+ vector< vector<double> > CrossCorrRMS(NumPolarizations, NumPairs - NumAntennae);
for(int j = 0; j < NumAntennae; j++)
- {
+ {
for(int k = j+1; k < NumAntennae; k++)
- {
+ {
int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- if (!FlagCompleteRow[index])
- { if( CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] < median/10 //a very complex way of indexing, could be simplified?
- || CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] > median*10)
- { FlagCompleteRow[index] = true;
- }
+ Matrix<Bool> image = Image.xyPlane(index);
+ if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
+ { //we skip the non-existent telescopes
+ test = ImageBaselineBand(&image,
+ &(TimeslotData[index]),
+ Position);
}
- }
+ }
}
- }
- }
- //this code could be separated in to different functions
- for(int j = 0; j < NumAntennae; j++) //write the data
- {
- for(int k = j; k < NumAntennae; k++)
- {
- int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- Matrix<Bool> flags = Flags.xyPlane(index);
- if (FlagCompleteRow[index])
+ //this code could be separated in to different functions
+ for(int j = 0; j < NumAntennae; j++) //write the data
{
- Matrix<Bool> flags = Flags.xyPlane(index);
- flags = true;
- Statistics(i,j,k) += NumChannels * NumPolarizations;
+ for(int k = j; k < NumAntennae; k++)
+ {
+ int index = i*NumPairs + BaselineIndex[pairii(j, k)];
+ Matrix<bool> image = Image.xyPlane(index);
+ (*Imagefile).WriteImageData(&image);
+ }
}
- else
- { Statistics(i,j,k) += stats[index];
- }
- (*MSfile).WriteDataPointFlags(flag_iter, &flags, FlagCompleteRow[index], ExistingFlags);
- (*flag_iter)++;
- }
+ }*/
}
- }
-}
-//===============>>> ComplexMedianFlagger::UpdateTimeslotData <<<===============
-/* This function reads the visibility data for one timeslot and checks if
- a mosaicing mode is active*/
-bool ComplexMedianFlagger::UpdateTimeslotData(vector<int>* OldFields,
- vector<int>* OldBands,
- int* TimeCounter,
- Table* TimeslotTable,
- double* Time)
-{
- int rowcount = (*TimeslotTable).nrow();
- ROTableVector<Int> antenna1((*TimeslotTable), "ANTENNA1");
- ROTableVector<Int> antenna2((*TimeslotTable), "ANTENNA2");
- ROTableVector<Int> bandnr ((*TimeslotTable), "DATA_DESC_ID");
- ROTableVector<Int> fieldid ((*TimeslotTable), "FIELD_ID");
- ROTableVector<Double> time ((*TimeslotTable), "TIME_CENTROID");//for testing purposes
- ROArrayColumn<Complex> data ((*TimeslotTable), "DATA");
- Cube<Complex> Data(NumPolarizations, NumChannels, rowcount);
+ //===============>>> DFTImager::UpdateTimeslotData <<<===============
+ /* This function reads the visibility data for one timeslot and checks if
+ a mosaicing mode is active*/
+ /* The datastructure Timeslotdata is rather complex because it flattens
+ half-filled antenna x antenna x bands matrix into a vector of NumPairs X bands length. */
+ bool DFTImager::UpdateTimeslotData(vector<int>* OldFields,
+ vector<int>* OldBands,
+ int* TimeCounter,
+ Table* TimeslotTable,
+ double* Time)
+ {
+ int rowcount = (*TimeslotTable).nrow();
+ ROTableVector<Int> antenna1((*TimeslotTable), "ANTENNA1");
+ ROTableVector<Int> antenna2((*TimeslotTable), "ANTENNA2");
+ ROTableVector<Int> bandnr ((*TimeslotTable), "DATA_DESC_ID");
+ ROTableVector<Int> fieldid ((*TimeslotTable), "FIELD_ID");
+ ROTableVector<Double> time ((*TimeslotTable), "TIME_CENTROID");//for testing purposes
+ ROArrayColumn<Complex> data ((*TimeslotTable), "DATA");
+ Cube<Complex> Data(NumPolarizations, NumChannels, rowcount);
- data.getColumn(Data); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
- (*Time) = time(0);//for testing purposes, might be useful in the future
-
- bool NewField_or_Frequency = false;
-
- for (int i = 0; i < rowcount; i++)
- {
- int bi = BaselineIndex[pairii(antenna1(i), antenna2(i))];
- int field = fieldid(i);
- int band = bandnr(i);
- int index = (band % NumBands) * NumPairs + bi;
- if ((*TimeCounter) > WindowSize - 1)
- {
- if (field != (*OldFields)[bi]) //pointing mosaicing
- { NewField_or_Frequency = true;
- }
- if (band != (*OldBands)[index]) //frequency mosaicing
- { NewField_or_Frequency = true;
- }
- }
- (*OldFields)[bi] = field;
- (*OldBands)[index] = band;
-
- TimeslotData[index].xyPlane((*TimeCounter) % WindowSize) = Data.xyPlane(i); //TimeslotData is only WindowSize size!
- }
- (*TimeCounter)++; //we want to reset if we detect a gap in DATA_DESC_ID or FIELD. Maybe TIME too???
- return NewField_or_Frequency; //assuming everybody is changing field or frequency at the same time!
-}
+ data.getColumn(Data); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
+ (*Time) = time(0);//for testing purposes, might be useful in the future
-//===============>>> ComplexMedianFlagger::FlagDataOrBaselines <<<===============
-/* This function iterates over the data per timeslot and uses Flagtimeslot()
- to actually flag datapoints (if flagDatapoints), and entire baselines (if flagRMS)*/
-void ComplexMedianFlagger::FlagDataOrBaselines(bool flagDatapoints,
- bool flagRMS,
- bool ExistingFlags)
-{
- TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
- TableIterator flag_iter = (*MSfile).TimeAntennaIterator();
- int TimeCounter = 0;
- double Time = 0.0;//for testing purposes
- vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
- vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
-
- int step = NumTimeslots / 10 + 1; //not exact but it'll do
- int row = 0;
-
- while (!timeslot_iter.pastEnd())
- {
- Table TimeslotTable = timeslot_iter.table();
- bool NewFieldorFreq = UpdateTimeslotData(&OldFields,
- &OldBands,
- &TimeCounter,
- &TimeslotTable,
- &Time);
- //cout << "Processing: " << MVTime(Time/(24*3600)).string(MVTime::YMD) << endl; //for testing purposes
-
- if (TimeCounter == WindowSize - 1)
- { //We have filled WindowSize timeslots and need to flag the first WindowSize/2 timeslots
- for (int position = 0; position < WindowSize/2; position++)
- { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags, position);
+ bool NewField_or_Frequency = false;
+
+ for (int i = 0; i < rowcount; i++)
+ {
+ int bi = BaselineIndex[pairii(antenna1(i), antenna2(i))];
+ int field = fieldid(i);
+ int band = bandnr(i);
+ int index = (band % NumBands) * NumPairs + bi;
+ if ((*TimeCounter) > WindowSize - 1)
+ {
+ if (field != (*OldFields)[bi]) //pointing mosaicing
+ { NewField_or_Frequency = true;
+ }
+ if (band != (*OldBands)[index]) //frequency mosaicing
+ { NewField_or_Frequency = true;
+ }
+ }
+ (*OldFields)[bi] = field;
+ (*OldBands)[index] = band;
+
+ TimeslotData[index].xyPlane((*TimeCounter) % WindowSize) = Data.xyPlane(i); //TimeslotData is only WindowSize size!
}
+ (*TimeCounter)++; //we want to reset if we detect a gap in DATA_DESC_ID or FIELD. Maybe TIME too???
+ return NewField_or_Frequency; //assuming everybody is changing field or frequency at the same time!
}
- if (TimeCounter > WindowSize - 1) //nothing special just falg a timeslot
- { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags,
- (TimeCounter + WindowSize/2) % WindowSize);
- }
- timeslot_iter++;
- if (row++ % step == 0) // to tell the user how much % we have processed,
- { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
- }
- if (timeslot_iter.pastEnd() || NewFieldorFreq)
- { //We still need to flag the last WindowSize/2 timeslots
- for (int position = WindowSize/2 + 1; position < WindowSize; position++)
- { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags,
- (TimeCounter + position) % WindowSize);
+
+ //===============>>> DFTImager::FlagDataOrBaselines <<<===============
+ /* This function iterates over the data per timeslot and uses Flagtimeslot()
+ to actually flag datapoints (if flagDatapoints), and entire baselines (if flagRMS)*/
+ void DFTImager::MakeImage(int Resolution)
+ {
+ TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
+ int TimeCounter = 0;
+ double Time = 0.0;//for testing purposes
+ vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
+ vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
+ vector< Cube<float> > Image(NumBands);
+
+ int step = NumTimeslots / 10 + 1; //not exact but it'll do
+ int row = 0;
+
+ while (!timeslot_iter.pastEnd())
+ {
+ Table TimeslotTable = timeslot_iter.table();
+ bool NewFieldorFreq = UpdateTimeslotData(&OldFields,
+ &OldBands,
+ &TimeCounter,
+ &TimeslotTable,
+ &Time);
+ if (NewFieldorFreq)
+ {
+ cout << "Error: current version can not handle multiple fields/frequqncies in dataset." << endl;
+ break;
+ }
+ cout << "Processing: " << MVTime(Time/(24*3600)).string(MVTime::YMD) << endl; //for testing purposes
+
+ ImageTimeslot();
+ timeslot_iter++;
+ if (row++ % step == 0) // to tell the user how much % we have processed,
+ { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
+ }
}
- TimeCounter = 0; //reset because we have changed to a new Field or frequency
+ Imagefile->WriteImage(&Image);
}
- }
- ProcessStatistics(); //is there a baseline that should be flagged?
-}
-//===============>>> ComplexMedianFlagger <<<===============
+ //===============>>> DFTImager <<<===============
+ } //namespace CS1
+}; //namespace LOFAR
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc~ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc~
index 02f185171ad..b164325ffff 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.cc~
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -21,19 +21,14 @@
#include <tables/Tables/TableIter.h>
#include "DFTImager.h"
#include <casa/Quanta/MVEpoch.h>
-#include <images/Images/PagedImage.h>
-#include <images/Images/ImageInfo.h>
-#include <images/Images/SubImage.h>
-#include <images/Images/ImageUtilities.h>
-using namespace WSRT; //should be nfra ?
using namespace casa;
enum CorrelationTypes {None=0,I=1,Q=2,U=3,V=4,RR=5,RL=6,LR=7,LL=8,XX=9,XY=10,YX=11,YY=12}; //found somewhere in AIPS++, don't remember where
//===============>>> itoa <<<===============
//ANSI C++ doesn't seem to have a decent function for this, or I'm not aware of it. Need to rename it to IntToStr(), to avoid confusion
-std::string itoa(int value, int base=10)
+std::string itoa(int value, int base=10)
{ // found on http://www.jb.man.ac.uk/~slowe/cpp/itoa.html
//maybe copyright Robert Jan Schaper, Ray-Yuan Sheu, Rodrigo de Salvo Braz, Wes Garland and John Maloney
enum { kMaxDigits = 35 };
@@ -53,515 +48,248 @@ std::string itoa(int value, int base=10)
return buf;
}
-//===============>>> ComplexMedianFlagger::ComplexMedianFlagger <<<===============
-/* initialize some meta data and get the datastorage the right size. */
-ComplexMedianFlagger::ComplexMedianFlagger(MS_File* InputMSfile,
- int InputWindowSize,
- bool UseOnlyXpolarizations,
- double InputMinThreshold,
- double InputMaxThreshold)
+namespace LOFAR
{
- MSfile = InputMSfile;
- WindowSize = InputWindowSize;
- MinThreshold = InputMinThreshold;
- MaxThreshold = InputMaxThreshold;
- NumAntennae = (*MSfile).itsNumAntennae;
- NumPairs = (*MSfile).itsNumPairs;
- NumBands = (*MSfile).itsNumBands;
- NumChannels = (*MSfile).itsNumChannels;
- NumPolarizations = (*MSfile).itsNumPolarizations;
- NoiseLevel = (*MSfile).itsNoiseLevel;
- NumTimeslots = (*MSfile).itsNumTimeslots;
- AntennaNames = (*MSfile).itsAntennaNames;
-
- PairsIndex.resize(NumPairs);
- Statistics = Cube<int>(NumBands, NumAntennae, NumAntennae, 0);
- PolarizationsToCheck.resize(NumPolarizations);
- DeterminePolarizationsToCheck(UseOnlyXpolarizations);
-
- TimeslotData.resize(NumPairs*NumBands);
- for (int i = 0; i < NumPairs*NumBands; i++)
- { TimeslotData[i].resize(NumPolarizations, NumChannels, WindowSize);
- }
-
- int index = 0;
- for (int i = 0; i < NumAntennae; i++)
- { for(int j = i; j < NumAntennae; j++)
- { PairsIndex[index] = pairii(i, j);
- BaselineIndex[pairii(i, j)] = index++;
- }
- }
- ComputeBaselineLengths();
-}
-
-//===============>>> ComplexMedianFlagger::~ComplexMedianFlagger <<<===============
-
-ComplexMedianFlagger::~ComplexMedianFlagger()
-{
-}
-
-//===============>>> ComplexMedianFlagger::DetermineCorrelationsToCheck <<<===============
-/* create a list of polarizations we want to check, maybe we only want to to XY, YX */
-void ComplexMedianFlagger::DeterminePolarizationsToCheck(bool UseOnlyXpolarizations)
-{
- if (UseOnlyXpolarizations)
+ namespace CS1
{
- bool noCorrError = true;
- for (int i = 0; i < NumPolarizations; i++)
- {
- switch((*MSfile).itsPolarizations[i])
- {
- case None:
- case I:
- case RR:
- case LL:
- case XX:
- case YY:
- if (UseOnlyXpolarizations)
- PolarizationsToCheck[i] = false;
- break;
- case Q:
- case U:
- case V:
- case RL:
- case LR:
- case XY:
- case YX:
- noCorrError = false;
- if (UseOnlyXpolarizations)
- PolarizationsToCheck[i] = true;
- break;
- }
- }
- if (noCorrError)
+ //===============>>> DFTImager::DFTImager <<<===============
+ /* initialize some meta data and get the datastorage the right size. */
+ DFTImager::DFTImager(MS_File* InputMSfile,
+ Image_File*InputImageFile)
{
- cout << "There are no crosspolarizations to flag!";
- exit(1);
- }
- }
- else
- {
- for (int i = 0; i < NumPolarizations; i++)
- { PolarizationsToCheck[i] = true;
- }
- }
-}
+ MSfile = InputMSfile;
+ Imagefile = InputImageFile;
+ NumAntennae = (*MSfile).itsNumAntennae;
+ NumPairs = (*MSfile).itsNumPairs;
+ NumBands = (*MSfile).itsNumBands;
+ NumChannels = (*MSfile).itsNumChannels;
+ NumPolarizations = (*MSfile).itsNumPolarizations;
+ NumTimeslots = (*MSfile).itsNumTimeslots;
+ AntennaNames = (*MSfile).itsAntennaNames;
+ Resolution = 41;
-//===============>>> ComplexMedianFlagger::ComputeBaselineLengths <<<===============
-/* compute baseline lengths, and determine the longest one.*/
-void ComplexMedianFlagger::ComputeBaselineLengths()
-{
- MaxBaselineLength = 0.0;
- BaselineLengths.resize(NumPairs);
- //Antenna positions
- MSAntenna antenna = (*MSfile).antenna();
- ROArrayColumn<Double> position(antenna, "POSITION");
- for (int i = 0; i < NumAntennae; i++ )
- {
- for (int j = i; j < NumAntennae; j++)
- {
- Vector<Double> p(position(i) - position(j));
- double temp = sqrt(p(0)*p(0) + p(1)*p(1) + p(2)*p(2));
- BaselineLengths[BaselineIndex[pairii(i, j)]] = temp;
- if (temp > MaxBaselineLength && temp < 3000000) //radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243
- { MaxBaselineLength = temp; // non-existent antenna's can have position (0,0,0)
- }
- }
- }
-}
+ PairsIndex.resize(NumPairs);
+ // Statistics = Cube<int>(NumBands, NumAntennae, NumAntennae, 0);
+ // PolarizationsToCheck.resize(NumPolarizations);
+ // DeterminePolarizationsToCheck(UseOnlyXpolarizations);
-//===============>>> ComplexMedianFlagger::FlagDataOrBaselines <<<===============
-/*This function outputs the gathered statistics.*/
-void ComplexMedianFlagger::ProcessStatistics()
-{
- vector<int> bands(NumBands);
- vector<int> antennae(NumAntennae);
- unsigned int namelength = 6;
- for(int i = 0; i < NumAntennae; i++)
- {
- if (namelength < AntennaNames[i].size())
- { namelength = AntennaNames[i].size();
- }
- }
- for (int i = 0; i < NumBands; i++)
- {
- cout << "Band: " << i+1 << endl;
- cout << string(namelength+1,' ');
- for(int j = 0; j < NumAntennae; j++)
- {
- string out = AntennaNames[j];
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl;
- for(int j = 0; j < NumAntennae; j++)
- {
- string out = AntennaNames[j];
- out.resize(namelength+1,' ');
- cout << out;
- for(int k = 0; k < NumAntennae; k++)
- {
- if (k < j) //We print a complete array, but we have inspected only those where k >= j
- {
- int val = 100 * Statistics(i,k,j) / (NumChannels * NumTimeslots * NumPolarizations);
- bands[i] += val;
- antennae[j] += val;
- antennae[k] += val;
- string out = itoa(val) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- else
- {
- int val = 100 * Statistics(i,j,k) / (NumChannels * NumTimeslots * NumPolarizations);
- bands[i] += val;
- antennae[j] += val;
- antennae[k] += val;
- string out = itoa(val) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- /*{ cout << rms << "Faulty baseline detected: Antenna " << j+1
- << ":" << k+1 << " SpectralWindow: "<< i+1 << endl;
- }
- }*/
+ TimeslotData.resize(NumPairs*NumBands);
+ for (int i = 0; i < NumPairs*NumBands; i++)
+ { TimeslotData[i].resize(NumPolarizations, NumChannels, WindowSize);
}
- cout << endl;
- }
- }
- cout << "Bands (flagged %): ";
- for (int i = 0; i < NumBands; i++)
- {
- string out = string("IVC") + itoa(i);
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl << " ";
- for (int i = 0; i < NumBands; i++)
- {
- string out = itoa(bands[i] / (NumAntennae*NumAntennae)) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl << "Antennae (flagged %): " ;
- for(int j = 0; j < NumAntennae; j++)
- {
- string out = AntennaNames[j];
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl << " ";
- for (int i = 0; i < NumAntennae; i++)
- {
- string out = itoa(antennae[i] / (NumBands*NumAntennae*2)) + "%";
- out.resize(namelength+1,' ');
- cout << out;
- }
- cout << endl;
-}
-//===============>>> ComplexMedianFlagger::FlagTimeslot <<<===============
-/* This function inspects each visibility in a cetain baseline-band
-and flags on complexe distance, then determines to flag the entire baseline-band
-based on the RMS of the points it didn't flag.*/
-bool ComplexMedianFlagger::FlagBaselineBand(Matrix<Bool>* Flags,
- Cube<Complex>* Timeslots,
- vector<double>* rms,
- int* flagCounter,
- double FlagThreshold,
- int Position,
- bool flagRMS)
-{
- Vector<Float> Reals(WindowSize);
- Vector<Float> Imags(WindowSize);
- vector<double> RMS(NumPolarizations);
- vector<int> RMSCounter(NumPolarizations);
- bool FlagCompleteRow = true;
- int flagcount = 0;
- for (int i = NumChannels-1; i >= 0; i--)
- {
- bool FlagAllPolarizations = false;
- for (int j = NumPolarizations-1; j >= 0; j--)
- { //we need to loop twice, once to determine FlagAllCorrelations
- if (!FlagAllPolarizations && PolarizationsToCheck[j])
- {
- for (int k = 0; k < WindowSize; k++)
- { //This might be faster in some other way ?
- Reals[k] = (*Timeslots)(j, i, k).real();
- Imags[k] = (*Timeslots)(j, i, k).imag();
+
+ int index = 0;
+ for (int i = 0; i < NumAntennae; i++)
+ { for(int j = i; j < NumAntennae; j++)
+ { PairsIndex[index] = pairii(i, j);
+ BaselineIndex[pairii(i, j)] = index++;
}
- float real = medianInPlace(Reals, false);
- float imag = medianInPlace(Imags, false);
- FlagAllPolarizations |= FlagThreshold < abs((*Timeslots)(j, i, Position)
- - Complex(real, imag));
- }
- }
- for (int j = NumPolarizations-1; j >= 0; j--)
- { //the second loop we set the flags or calculate RMS
- if (FlagAllPolarizations)
- { (*Flags)(j, i) = true;
- flagcount++;
- }
- else
- {
- FlagCompleteRow = false;
- (*Flags)(j, i) = false; //could check existing flag ?
- double value = abs((*Timeslots)(j, i, Position));
- RMS[j] += value*value;
- (RMSCounter[j])++;
}
+ ComputeBaselineLengths();
}
- }
- //these need to be separated out into a different function for clarity
- if (flagcount > 0.9 * NumChannels * NumPolarizations) //more as 90% bad
- {
- FlagCompleteRow = true;
- (*Flags) = true;
- flagcount = NumChannels * NumPolarizations;
- }
- if (!FlagCompleteRow && flagRMS) //RMSCounter > 0
- {
- for (int j = NumPolarizations-1; j >= 0; j--)
- {
- (*rms)[j] = sqrt( RMS[j] / RMSCounter[j] );
-// if (PolarizationsToCheck[j] && ((*rms)[j] > 1.5 || (*rms)[j] <= 0.5 * NoiseLevel)) //we think this RMS value is non-physical
- if (PolarizationsToCheck[j] && ((*rms)[j] > 25 || (*rms)[j] <= 0.5 * NoiseLevel)) //we think this RMS value is non-physical
- {
- FlagCompleteRow = true;
- (*Flags) = true;
- flagcount = NumChannels * NumPolarizations;
- }
- }
- }
- (*flagCounter) += flagcount;
- return FlagCompleteRow;
-}
-//===============>>> ComplexMedianFlagger::FlagBaseline <<<===============
-/* This function iterates over baseline and band and uses FlagBaselineBand() to determine
- for each one if it needs to be flagged. It treats the autocorrelations separately,
- to detect entire malfunctioning telescopes. Finally it writes the flags.
-*/
-void ComplexMedianFlagger::FlagTimeslot(TableIterator* flag_iter,
- bool flagDatapoints,
- bool flagRMS,
- bool ExistingFlags,
- Int Position)
-{
- Cube<bool> Flags(NumPolarizations, NumChannels, NumPairs*NumBands, false);
- vector<bool> FlagCompleteRow(NumPairs*NumBands);
- vector<int> stats(NumPairs*NumBands);
- for (int i = 0; i < NumBands; i++)
- {
- vector<double> rms(NumPolarizations);
- vector< vector<double> > AutoCorrRMS(NumPolarizations, NumAntennae);
- for (int j = 0; j < NumPolarizations; j++) // this needs to be written more elegantly
- {
- if ((*MSfile).itsPolarizations[j] == YX)
- { PolarizationsToCheck[j] = false; // We're not checking YX because WSRT writes 0.0 into them
- }
- }
- for(int j = 0; j < NumAntennae; j++)
- { //check auto correlations
- int index = i*NumPairs + BaselineIndex[pairii(j, j)];
- Matrix<Bool> flags = Flags.xyPlane(index);
- FlagCompleteRow[index] = FlagBaselineBand(&flags,
- &(TimeslotData[index]),
- &rms, &(stats[index]),
- 3000 * NoiseLevel,
- Position, flagRMS);
- for(int k = 0; k < NumPolarizations; k++)
- { AutoCorrRMS[k][j] = rms[k];
- }
- }
- for(int k = 0; k < NumPolarizations; k++) //check on median
+
+ //===============>>> DFTImager::~DFTImager <<<===============
+
+ DFTImager::~DFTImager()
{
- if (PolarizationsToCheck[k])
- {
- double median = medianInPlace(Vector<double>(AutoCorrRMS[k]), false);
- for(int j = 0; j < NumAntennae; j++)
- {
- int index = i*NumPairs + BaselineIndex[pairii(j, j)];
- if (AutoCorrRMS[k][j] < median/10 || AutoCorrRMS[k][j] > median*10)
- { FlagCompleteRow[index] = true;
- }
- }
- }
}
- for (int j = 0; j < NumPolarizations; j++) // this needs to be written more elegantly
+
+ //===============>>> DFTImager::DetermineCorrelationsToCheck <<<===============
+ /* create a list of polarizations we want to check, maybe we only want to to XY, YX */
+ // void DFTImager::DeterminePolarizationsToCheck(bool UseOnlyXpolarizations)
+ // {
+ // if (UseOnlyXpolarizations)
+ // {
+ // bool noCorrError = true;
+ // for (int i = 0; i < NumPolarizations; i++)
+ // {
+ // switch((*MSfile).itsPolarizations[i])
+ // {
+ // case None:
+ // case I:
+ // case RR:
+ // case LL:
+ // case XX:
+ // case YY:
+ // if (UseOnlyXpolarizations)
+ // PolarizationsToCheck[i] = false;
+ // break;
+ // case Q:
+ // case U:
+ // case V:
+ // case RL:
+ // case LR:
+ // case XY:
+ // case YX:
+ // noCorrError = false;
+ // if (UseOnlyXpolarizations)
+ // PolarizationsToCheck[i] = true;
+ // break;
+ // }
+ // }
+ // if (noCorrError)
+ // {
+ // cout << "There are no crosspolarizations to flag!";
+ // exit(1);
+ // }
+ // }
+ // else
+ // {
+ // for (int i = 0; i < NumPolarizations; i++)
+ // { PolarizationsToCheck[i] = true;
+ // }
+ // }
+ // }
+
+ //===============>>> DFTImager::ComputeBaselineLengths <<<===============
+ /* compute baseline lengths, and determine the longest one.*/
+ void DFTImager::ComputeBaselineLengths()
{
- if ((*MSfile).itsPolarizations[j] == YX)
- { PolarizationsToCheck[j] = true; // We're not checking YX
- }
- }
- // this code coud be separated into different functions
- vector< vector<double> > CrossCorrRMS(NumPolarizations, NumPairs - NumAntennae);
- for(int j = 0; j < NumAntennae; j++)
- {
- for(int k = j+1; k < NumAntennae; k++)
- {
- int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- Matrix<Bool> flags = Flags.xyPlane(index);
- if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
- { //we skip the non-existent telescopes
- FlagCompleteRow[index] = FlagBaselineBand(&flags,
- &(TimeslotData[index]),
- &rms, &(stats[index]),
- (MinThreshold + (MaxThreshold - MinThreshold)
- * BaselineLengths[BaselineIndex[pairii(j, k)]]
- / MaxBaselineLength
- ) * NoiseLevel,
- Position, flagRMS);
- for(int l = 0; l < NumPolarizations; l++)
- { CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] = rms[l]; //a very complex way of indexing, could be simplified?
- }
- if ( ( FlagCompleteRow[i*NumPairs + BaselineIndex[pairii(j, j)]]
- ||FlagCompleteRow[i*NumPairs + BaselineIndex[pairii(k, k)]])
- &&flagRMS && !FlagCompleteRow[index])
- { //We flag it if the autocorrelation is flagged
- FlagCompleteRow[index] = true;
+ MaxBaselineLength = 0.0;
+ BaselineLengths.resize(NumPairs);
+ //Antenna positions
+ MSAntenna antennas = (*MSfile).antennas();
+ ROArrayColumn<Double> position(antennas, "POSITION");
+ for (int i = 0; i < NumAntennae; i++ )
+ {
+ for (int j = i; j < NumAntennae; j++)
+ {
+ Vector<Double> p(position(i) - position(j));
+ double temp = sqrt(p(0)*p(0) + p(1)*p(1) + p(2)*p(2));
+ BaselineLengths[BaselineIndex[pairii(i, j)]] = temp;
+ if (temp > MaxBaselineLength && temp < 3000000) //radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243
+ { MaxBaselineLength = temp; // non-existent antenna's can have position (0,0,0)
}
- }
+ }
}
- }
- for(int l = 0; l < NumPolarizations; l++) // check on median
+ }
+
+ //===============>>> DFTImager::ImageTimeslot <<<===============
+ /* This function iterates over baseline and band and uses FlagBaselineBand() to determine
+ for each one if it needs to be flagged. It treats the autocorrelations separately,
+ to detect entire malfunctioning telescopes. Finally it writes the flags.
+ */
+/* void DFTImager::ImageTimeslot()
{
- if (PolarizationsToCheck[l])
- {
- double median = medianInPlace(Vector<double>(CrossCorrRMS[l]), false);
+ Cube<bool> Image(NumPolarizations, NumChannels, NumPairs*NumBands, false);
+ vector<bool> test(NumPairs*NumBands);
+ for (int i = 0; i < NumBands; i++)
+ {
+ vector< vector<double> > CrossCorrRMS(NumPolarizations, NumPairs - NumAntennae);
for(int j = 0; j < NumAntennae; j++)
- {
+ {
for(int k = j+1; k < NumAntennae; k++)
- {
+ {
int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- if (!FlagCompleteRow[index])
- { if( CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] < median/10 //a very complex way of indexing, could be simplified?
- || CrossCorrRMS[l][j*NumAntennae + k - ((j+1)*(j+2))/2] > median*10)
- { FlagCompleteRow[index] = true;
- }
+ Matrix<Bool> image = Image.xyPlane(index);
+ if ((BaselineLengths[BaselineIndex[pairii(j, k)]] < 3000000))//radius of the Earth in meters? WSRT sometimes has fake telescopes at 3854243 m
+ { //we skip the non-existent telescopes
+ test = ImageBaselineBand(&image,
+ &(TimeslotData[index]),
+ Position);
}
- }
+ }
}
- }
- }
- //this code could be separated in to different functions
- for(int j = 0; j < NumAntennae; j++) //write the data
- {
- for(int k = j; k < NumAntennae; k++)
- {
- int index = i*NumPairs + BaselineIndex[pairii(j, k)];
- Matrix<Bool> flags = Flags.xyPlane(index);
- if (FlagCompleteRow[index])
+ //this code could be separated in to different functions
+ for(int j = 0; j < NumAntennae; j++) //write the data
{
- Matrix<Bool> flags = Flags.xyPlane(index);
- flags = true;
- Statistics(i,j,k) += NumChannels * NumPolarizations;
+ for(int k = j; k < NumAntennae; k++)
+ {
+ int index = i*NumPairs + BaselineIndex[pairii(j, k)];
+ Matrix<bool> image = Image.xyPlane(index);
+ (*Imagefile).WriteImageData(&image);
+ }
}
- else
- { Statistics(i,j,k) += stats[index];
- }
- (*MSfile).WriteDataPointFlags(flag_iter, &flags, FlagCompleteRow[index], ExistingFlags);
- (*flag_iter)++;
}
- }
- }
-}
+ }*/
-//===============>>> ComplexMedianFlagger::UpdateTimeslotData <<<===============
-/* This function reads the visibility data for one timeslot and checks if
- a mosaicing mode is active*/
-bool ComplexMedianFlagger::UpdateTimeslotData(vector<int>* OldFields,
- vector<int>* OldBands,
- int* TimeCounter,
- Table* TimeslotTable,
- double* Time)
-{
- int rowcount = (*TimeslotTable).nrow();
- ROTableVector<Int> antenna1((*TimeslotTable), "ANTENNA1");
- ROTableVector<Int> antenna2((*TimeslotTable), "ANTENNA2");
- ROTableVector<Int> bandnr ((*TimeslotTable), "DATA_DESC_ID");
- ROTableVector<Int> fieldid ((*TimeslotTable), "FIELD_ID");
- ROTableVector<Double> time ((*TimeslotTable), "TIME_CENTROID");//for testing purposes
- ROArrayColumn<Complex> data ((*TimeslotTable), "DATA");
- Cube<Complex> Data(NumPolarizations, NumChannels, rowcount);
+ //===============>>> DFTImager::UpdateTimeslotData <<<===============
+ /* This function reads the visibility data for one timeslot and checks if
+ a mosaicing mode is active*/
+ /* The datastructure Timeslotdata is rather complex because it flattens
+ half-filled antenna x antenna x bands matrix into a vector of NumPairs X bands length. */
+ bool DFTImager::UpdateTimeslotData(vector<int>* OldFields,
+ vector<int>* OldBands,
+ int* TimeCounter,
+ Table* TimeslotTable,
+ double* Time)
+ {
+ int rowcount = (*TimeslotTable).nrow();
+ ROTableVector<Int> antenna1((*TimeslotTable), "ANTENNA1");
+ ROTableVector<Int> antenna2((*TimeslotTable), "ANTENNA2");
+ ROTableVector<Int> bandnr ((*TimeslotTable), "DATA_DESC_ID");
+ ROTableVector<Int> fieldid ((*TimeslotTable), "FIELD_ID");
+ ROTableVector<Double> time ((*TimeslotTable), "TIME_CENTROID");//for testing purposes
+ ROArrayColumn<Complex> data ((*TimeslotTable), "DATA");
+ Cube<Complex> Data(NumPolarizations, NumChannels, rowcount);
- data.getColumn(Data); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
- (*Time) = time(0);//for testing purposes, might be useful in the future
-
- bool NewField_or_Frequency = false;
-
- for (int i = 0; i < rowcount; i++)
- {
- int bi = BaselineIndex[pairii(antenna1(i), antenna2(i))];
- int field = fieldid(i);
- int band = bandnr(i);
- int index = (band % NumBands) * NumPairs + bi;
- if ((*TimeCounter) > WindowSize - 1)
- {
- if (field != (*OldFields)[bi]) //pointing mosaicing
- { NewField_or_Frequency = true;
- }
- if (band != (*OldBands)[index]) //frequency mosaicing
- { NewField_or_Frequency = true;
- }
- }
- (*OldFields)[bi] = field;
- (*OldBands)[index] = band;
-
- TimeslotData[index].xyPlane((*TimeCounter) % WindowSize) = Data.xyPlane(i); //TimeslotData is only WindowSize size!
- }
- (*TimeCounter)++; //we want to reset if we detect a gap in DATA_DESC_ID or FIELD. Maybe TIME too???
- return NewField_or_Frequency; //assuming everybody is changing field or frequency at the same time!
-}
+ data.getColumn(Data); //We're not checking Data.nrow() Data.ncolumn(), assuming all data is the same size.
+ (*Time) = time(0);//for testing purposes, might be useful in the future
-//===============>>> ComplexMedianFlagger::FlagDataOrBaselines <<<===============
-/* This function iterates over the data per timeslot and uses Flagtimeslot()
- to actually flag datapoints (if flagDatapoints), and entire baselines (if flagRMS)*/
-void ComplexMedianFlagger::FlagDataOrBaselines(bool flagDatapoints,
- bool flagRMS,
- bool ExistingFlags)
-{
- TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
- TableIterator flag_iter = (*MSfile).TimeAntennaIterator();
- int TimeCounter = 0;
- double Time = 0.0;//for testing purposes
- vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
- vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
-
- int step = NumTimeslots / 10 + 1; //not exact but it'll do
- int row = 0;
-
- while (!timeslot_iter.pastEnd())
- {
- Table TimeslotTable = timeslot_iter.table();
- bool NewFieldorFreq = UpdateTimeslotData(&OldFields,
- &OldBands,
- &TimeCounter,
- &TimeslotTable,
- &Time);
- //cout << "Processing: " << MVTime(Time/(24*3600)).string(MVTime::YMD) << endl; //for testing purposes
-
- if (TimeCounter == WindowSize - 1)
- { //We have filled WindowSize timeslots and need to flag the first WindowSize/2 timeslots
- for (int position = 0; position < WindowSize/2; position++)
- { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags, position);
+ bool NewField_or_Frequency = false;
+
+ for (int i = 0; i < rowcount; i++)
+ {
+ int bi = BaselineIndex[pairii(antenna1(i), antenna2(i))];
+ int field = fieldid(i);
+ int band = bandnr(i);
+ int index = (band % NumBands) * NumPairs + bi;
+ if ((*TimeCounter) > WindowSize - 1)
+ {
+ if (field != (*OldFields)[bi]) //pointing mosaicing
+ { NewField_or_Frequency = true;
+ }
+ if (band != (*OldBands)[index]) //frequency mosaicing
+ { NewField_or_Frequency = true;
+ }
+ }
+ (*OldFields)[bi] = field;
+ (*OldBands)[index] = band;
+
+ TimeslotData[index].xyPlane((*TimeCounter) % WindowSize) = Data.xyPlane(i); //TimeslotData is only WindowSize size!
}
+ (*TimeCounter)++; //we want to reset if we detect a gap in DATA_DESC_ID or FIELD. Maybe TIME too???
+ return NewField_or_Frequency; //assuming everybody is changing field or frequency at the same time!
}
- if (TimeCounter > WindowSize - 1) //nothing special just falg a timeslot
- { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags,
- (TimeCounter + WindowSize/2) % WindowSize);
- }
- timeslot_iter++;
- if (row++ % step == 0) // to tell the user how much % we have processed,
- { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
- }
- if (timeslot_iter.pastEnd() || NewFieldorFreq)
- { //We still need to flag the last WindowSize/2 timeslots
- for (int position = WindowSize/2 + 1; position < WindowSize; position++)
- { FlagTimeslot(&flag_iter, flagDatapoints, flagRMS, ExistingFlags,
- (TimeCounter + position) % WindowSize);
+
+ //===============>>> DFTImager::FlagDataOrBaselines <<<===============
+ /* This function iterates over the data per timeslot and uses Flagtimeslot()
+ to actually flag datapoints (if flagDatapoints), and entire baselines (if flagRMS)*/
+ void DFTImager::MakeImage(int Resolution)
+ {
+ TableIterator timeslot_iter = (*MSfile).TimeslotIterator();
+ int TimeCounter = 0;
+ double Time = 0.0;//for testing purposes
+ vector<int> OldFields(NumPairs); //to check on multipointing and mosaicing
+ vector<int> OldBands(NumPairs*NumBands); //to check on multifrequency and freq mosaicing
+ vector< Cube<float> > Image(NumBands);
+
+ int step = NumTimeslots / 10 + 1; //not exact but it'll do
+ int row = 0;
+
+ while (!timeslot_iter.pastEnd())
+ {
+ Table TimeslotTable = timeslot_iter.table();
+ bool NewFieldorFreq = UpdateTimeslotData(&OldFields,
+ &OldBands,
+ &TimeCounter,
+ &TimeslotTable,
+ &Time);
+ if (NewFieldorFreq)
+ {
+ cout << "Error: current version can not handle multiple fields/frequqncies in dataset." << endl;
+ break;
+ }
+ cout << "Processing: " << MVTime(Time/(24*3600)).string(MVTime::YMD) << endl; //for testing purposes
+
+ ImageTimeslot();
+ timeslot_iter++;
+ if (row++ % step == 0) // to tell the user how much % we have processed,
+ { cout << 10*(row/step) << "%" << endl; //not very accurate for low numbers of timeslots, but it'll do for now
+ }
}
- TimeCounter = 0; //reset because we have changed to a new Field or frequency
+ Imagefile->WriteImage(&Image);
}
- }
- ProcessStatistics(); //is there a baseline that should be flagged?
-}
-//===============>>> ComplexMedianFlagger <<<===============
+ //===============>>> DFTImager <<<===============
+ } //namespace CS1
+}; //namespace LOFAR
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h
index 4c6eb889b73..6767aa4b078 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -17,8 +17,8 @@
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
-#ifndef __FLAGGER_COMPLEXMEDIANFLAGGER_H__
-#define __FLAGGER_COMPLEXMEDIANFLAGGER_H__
+#ifndef __IMAGER_DFTIMAGER_H__
+#define __IMAGER_DFTIMAGER_H__
#include <casa/Arrays.h>
#include <utility>
@@ -41,31 +41,50 @@ namespace LOFAR
using casa::Table;
using std::list;
using std::vector;
-
+
typedef pair<int, int> pairii;
-
+
class DFTImager
{
public:
- DFTImager(MS_File* MSfile,
- PagedImage Imagefile,
- int Resolution);
+ DFTImager(MS_File* MSfile,
+ Image_File* Imagefile);
~DFTImager();
-
- void MakeImage();
-
+
+ void MakeImage(int Resolution);
+
protected:
+ int NumAntennae;
+ int NumPairs;
+ int NumBands;
int NumChannels;
- MS_File* MSfile;
-
+ int NumPolarizations;
+ int WindowSize;
+ int NumTimeslots;
+ double MinThreshold;
+ double MaxThreshold;
+ double MaxBaselineLength;
+ vector<double> BaselineLengths;
+ vector<pairii> PairsIndex;
+ map<pairii, int> BaselineIndex;
+ vector< Cube<Complex> > TimeslotData;
+ vector<casa::String> AntennaNames;
+ Cube< int > Statistics;
+ int Resolution;
+
+ MS_File* MSfile;
+ Image_File* Imagefile;
+
+ private:
+ void ComputeBaselineLengths();
bool UpdateTimeslotData(vector<int>* OldFields,
vector<int>* OldBands,
int* TimeCounter,
Table* TimeslotTable,
double* Time);
- private:
- }; // ComplexMedianFlagger
+ void ImageTimeslot(void);
+ }; // DFTImager
}; // namespace CS1
}; // namespace LOFAR
-#endif // __FLAGGER_COMPLEXMEDIANFLAGGER_H__
+#endif // __IMAGER_DFTIMAGER_H__
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h~ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h~
index 4c6eb889b73..a5abade15aa 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/DFTImager.h~
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -17,8 +17,8 @@
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
-#ifndef __FLAGGER_COMPLEXMEDIANFLAGGER_H__
-#define __FLAGGER_COMPLEXMEDIANFLAGGER_H__
+#ifndef __IMAGER_DFTIMAGER_H__
+#define __IMAGER_DFTIMAGER_H__
#include <casa/Arrays.h>
#include <utility>
@@ -41,31 +41,50 @@ namespace LOFAR
using casa::Table;
using std::list;
using std::vector;
-
+
typedef pair<int, int> pairii;
-
+
class DFTImager
{
public:
- DFTImager(MS_File* MSfile,
- PagedImage Imagefile,
- int Resolution);
+ DFTImager(MS_File* MSfile,
+ Image_File* Imagefile);
~DFTImager();
-
- void MakeImage();
-
+
+ void MakeImage(int Resolution);
+
protected:
+ int NumAntennae;
+ int NumPairs;
+ int NumBands;
int NumChannels;
- MS_File* MSfile;
-
+ int NumPolarizations;
+ int WindowSize;
+ int NumTimeslots;
+ double MinThreshold;
+ double MaxThreshold;
+ double MaxBaselineLength;
+ vector<double> BaselineLengths;
+ vector<pairii> PairsIndex;
+ map<pairii, int> BaselineIndex;
+ vector< Cube<Complex> > TimeslotData;
+ vector<casa::String> AntennaNames;
+ Cube< int > Statistics;
+ int Resolution;
+
+ MS_File* MSfile;
+ Image_File* Imagefile;
+
+ private:
+ void ComputeBaselineLengths();
bool UpdateTimeslotData(vector<int>* OldFields,
vector<int>* OldBands,
int* TimeCounter,
Table* TimeslotTable,
double* Time);
- private:
- }; // ComplexMedianFlagger
+ void ImageTimeslot();
+ }; // DFTImager
}; // namespace CS1
}; // namespace LOFAR
-#endif // __FLAGGER_COMPLEXMEDIANFLAGGER_H__
+#endif // __IMAGER_DFTIMAGER_H__
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc
index c99df156210..b20362ffdab 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc
@@ -28,6 +28,9 @@
#include<casa/aips.h>
#include <images/Images/PagedImage.h>
+#include <images/Images/ImageInfo.h>
+#include <images/Images/SubImage.h>
+#include <images/Images/ImageUtilities.h>
//#include <trial/ImgCrdSys/ImageCoordinate.h>
#include <casa/Arrays/Vector.h>
@@ -41,10 +44,47 @@
#include <casa/iostream.h>
#include <casa/namespace.h>
-main()
+#include "Image_File.h"
+
+
+using namespace casa;
+
+namespace LOFAR
{
- cout << "untested" << endl;
-}
+ namespace CS1
+ {
+
+ //===============>>> Image_file::Image_file <<<===============
+
+ Image_File::Image_File(const std::string& imagename)
+ {
+ ImageName = imagename;
+ }
+
+ //===============>>> Image_file::~Image_file <<<===============
+
+ Image_File::~Image_File()
+ {
+ delete Image;
+ }
+
+ //===============>>> Image_file::init <<<===============
+
+ void Image_File::init(casa::IPosition shape, size_t length)
+ {
+// Image = new casa::PagedImage(shape, NULL, Table::Update);
+ }
+
+ //===============>>> Image_file::~WriteImage <<<===============
+
+ void Image_File::WriteImage(std::vector< casa::Cube<float> >* image)
+ {
+ casa::Cube<float> Bal = (*image)[0];
+ init(Bal.shape(), image->size());
+ Image->putSlice(Bal, casa::IPosition(3, 0, 0, 0), casa::IPosition(3, 0, 0, 0));
+ }
+
+
/*
//#predeclarations
@@ -85,7 +125,7 @@ int main(int argc, char **argv)
cout << "OK" << endl;
} catch (AipsError x) {
cerr << "Caught Exception: " << x.getMesg() << endl;
- }
+ }
return 0;
}
@@ -141,7 +181,7 @@ void describeImage(const String &message, const PagedImage<Float> &image)
cout << " shape: "
<< image.coordinates().imageShape() << endl;
// * /
- cout << " ok: "<< image.ok()<< endl;
+ cout << " ok: "<< image.ok()<< endl;
}
// build a set of ImageCoordinates in three dimensions
@@ -161,7 +201,7 @@ ImageCoordinate build3Dcoords()
// Let's create an EarthPosition with full description of all parameters.
// We need a type - GEOCENTRIC seems good.
EarthPosition::Type theType(EarthPosition::GEOCENTRIC);
- // We need the time and date of the observation...
+ // We need the time and date of the observation...
Double julianDate = 2449376.0;
// and we can add on the UT.
julianDate += 16.52/24.0;
@@ -173,11 +213,11 @@ ImageCoordinate build3Dcoords()
ourPosition(1) = 34.1562394;
// geocentric radius (in meters) goes in the last field.
ourPosition(2) = 6372139.592;
- // then use these to build our EarthPosition.
+ // then use these to build our EarthPosition.
EarthPosition myObs(theType, julianDate, ourPosition);
// we need to know a rotation - here it is zero.
Double myRot = 0;
- // we need to know where the spherical position is to be on our 2-d
+ // we need to know where the spherical position is to be on our 2-d
// projection i.e. what pixel is associated with my object's position?
Vector<Double> thePixel(2);
thePixel(0) = 55.0;
@@ -190,10 +230,10 @@ ImageCoordinate build3Dcoords()
// Now we can make fruit of our labor - the ProjectedPosition itself.
ProjectedPosition myMapping(myMethod, myVectorType, myEpoch, myCoords,
myObs, myRot, thePixel, theBinning);
-
+
// we need to know what the units are of the measured value
MeasuredValue::Type myValueUnit(MeasuredValue::RADIO_VELOCITY);
- // we need to know what the above units are in reference
+ // we need to know what the above units are in reference
// here it is the velocity of the earth.
ReferenceValue myRefValue(ReferenceValue::VELOCITY, 9.56e+03);
// we need to know the value of the measurement
@@ -204,7 +244,7 @@ ImageCoordinate build3Dcoords()
Double myValuePos(27.3);
// Now we may construct the LinearAxis itself.
LinearAxis myLinearAxis(myValueUnit,myValue,myRefValue,myBin,myValuePos);
-
+
ImageCoordinate coords;
coords.addAxis(myMapping);
coords.addAxis(myLinearAxis);
@@ -215,20 +255,20 @@ ImageCoordinate build3Dcoords()
void createStandardImageOnDisk(const String &filename)
{
if(verbose_) cout<< "-- createStandardImageOnDisk --"<< endl;
-
+
ImageCoordinate coords(build3Dcoords());
IPosition imageShape(3, 50,47,31);
-
+
PagedImage<Float> standard(imageShape, coords, filename);
standard.set(TEST_PIXEL_VALUE);
-
+
if(verbose_) describeImage("standard image", standard);
}
// these are the constructors to test:
// PagedImage(const IPosition &shape, const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
-// PagedImage(const IPosition &shape, const Array<T> &array,
+// PagedImage(const IPosition &shape, const Array<T> &array,
// const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
// PagedImage(const String &filename, uInt rowNumber);
@@ -242,12 +282,12 @@ void testConstructors()
ImageCoordinate coords(build3Dcoords());
PagedImage<Float> image1(shape, coords, NAME0);
- if (verbose_)
+ if (verbose_)
describeImage("shape, coords, filename ctor: ", image1);
PagedImage<Float> image2(shape, coords, NAME1, True);
-
- if (verbose_)
+
+ if (verbose_)
describeImage("array, array, coords, filename, masking? ctor: ", image2);
PagedImage<Float> image3(STANDARD);
@@ -266,10 +306,10 @@ void testConstructors()
throw(AipsError("file-constructed image not ok"));
delete image4;
}
-
+
deleteFile(NAME0);
deleteFile(NAME1);
-}
+}
// test some of the data member manipulation functions
void testSetMemberFunctions()
@@ -279,7 +319,7 @@ void testSetMemberFunctions()
ImageCoordinate coords;
{
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
-
+
if(verbose_) describeImage("shape, coords, filename ctor: ", image6);
if(image6.ok())
throw(AipsError("shape-constructed image6 is ok, but shouldn't be!"));
@@ -288,20 +328,20 @@ void testSetMemberFunctions()
coords = build3Dcoords();
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
if(!image6.ok())
- throw(AipsError("shape-constructed image6 - is not ok!"));
+ throw(AipsError("shape-constructed image6 - is not ok!"));
image6.rename(NAME2);
if(verbose_) describeImage("image6, after rename", image6);
if(!image6.ok())
throw(AipsError("image6 after rename - not ok"));
-
+
image6.setCoordinateInfo(build3Dcoords());
if(verbose_) describeImage("image6, after setCoordinates", image6);
if(!image6.ok())
throw(AipsError("image6 after setCoordinates - not ok"));
}
-
+
deleteFile(NAME2);
-}
+}
//
void testArrayofImagesAndAssignmentOperator()
@@ -320,7 +360,7 @@ void testArrayofImagesAndAssignmentOperator()
filenames [7] = NAME7;
filenames [8] = NAME8;
filenames [9] = NAME9;
- {
+ {
Image<Float> images [max];
for(uInt i=0; i< max; i++) {
if(verbose_) cout<< "-- image array "<< i<< " --"<< endl;
@@ -328,14 +368,14 @@ void testArrayofImagesAndAssignmentOperator()
images [i] = Image<Float>(IPosition(2,10,10));
images [i].setCoordinateInfo(MinimalCoords());
images [i].setName(filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape constructed image, setCoords, setName in array",
images [i]);
} // if i is odd
else {
images [i] =
Image<Float>(IPosition(i+1,5), MinimalCoords(), filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape,coords,name constructed image in array,",
images [i]);
} // else
@@ -343,7 +383,7 @@ void testArrayofImagesAndAssignmentOperator()
throw(AipsError("one of array of images not ok"));
} // for i
} // scope
-
+
if(verbose_) cout<< "-- about to delete image files --"<< endl;
for(uInt i=0; i< max; i++)
@@ -351,3 +391,6 @@ void testArrayofImagesAndAssignmentOperator()
// * /
}
*/
+ //===============>>> Image_File <<<===============
+ }
+}
\ No newline at end of file
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc~ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc~
index 6de948f2231..8f7f6936338 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.cc~
@@ -28,7 +28,10 @@
#include<casa/aips.h>
#include <images/Images/PagedImage.h>
-///#include <trial/ImgCrdSys/ImageCoordinate.h>
+#include <images/Images/ImageInfo.h>
+#include <images/Images/SubImage.h>
+#include <images/Images/ImageUtilities.h>
+//#include <trial/ImgCrdSys/ImageCoordinate.h>
#include <casa/Arrays/Vector.h>
#include <casa/Arrays/Cube.h>
@@ -41,10 +44,47 @@
#include <casa/iostream.h>
#include <casa/namespace.h>
-main()
+#include "Image_File.h"
+
+
+using namespace casa;
+
+namespace LOFAR
{
- cout << "untested" << endl;
-}
+ namespace CS1
+ {
+
+ //===============>>> Image_file::Image_file <<<===============
+
+ Image_File::Image_File(const std::string& imagename)
+ {
+ ImageName = imagename;
+ }
+
+ //===============>>> Image_file::~Image_file <<<===============
+
+ Image_File::~Image_File()
+ {
+ delete Image;
+ }
+
+ //===============>>> Image_file::init <<<===============
+
+ void Image_File::init(casa::IPosition shape, size_t length)
+ {
+// Image = new casa::PagedImage(shape, NULL, Table::Update);
+ }
+
+ //===============>>> Image_file::~WriteImage <<<===============
+
+ void Image_File::WriteImage(std::vector< casa::Cube<float> >* image)
+ {
+ casa::Cube<float> Bal = (*image)[0];
+ init(Bal.shape(), image->size());
+ Image->putSlice((*image[0]), casa::IPosition(3, 0, 0, 0), casa::IPosition(3, 0, 0, 0));
+ }
+
+
/*
//#predeclarations
@@ -85,7 +125,7 @@ int main(int argc, char **argv)
cout << "OK" << endl;
} catch (AipsError x) {
cerr << "Caught Exception: " << x.getMesg() << endl;
- }
+ }
return 0;
}
@@ -141,7 +181,7 @@ void describeImage(const String &message, const PagedImage<Float> &image)
cout << " shape: "
<< image.coordinates().imageShape() << endl;
// * /
- cout << " ok: "<< image.ok()<< endl;
+ cout << " ok: "<< image.ok()<< endl;
}
// build a set of ImageCoordinates in three dimensions
@@ -161,7 +201,7 @@ ImageCoordinate build3Dcoords()
// Let's create an EarthPosition with full description of all parameters.
// We need a type - GEOCENTRIC seems good.
EarthPosition::Type theType(EarthPosition::GEOCENTRIC);
- // We need the time and date of the observation...
+ // We need the time and date of the observation...
Double julianDate = 2449376.0;
// and we can add on the UT.
julianDate += 16.52/24.0;
@@ -173,11 +213,11 @@ ImageCoordinate build3Dcoords()
ourPosition(1) = 34.1562394;
// geocentric radius (in meters) goes in the last field.
ourPosition(2) = 6372139.592;
- // then use these to build our EarthPosition.
+ // then use these to build our EarthPosition.
EarthPosition myObs(theType, julianDate, ourPosition);
// we need to know a rotation - here it is zero.
Double myRot = 0;
- // we need to know where the spherical position is to be on our 2-d
+ // we need to know where the spherical position is to be on our 2-d
// projection i.e. what pixel is associated with my object's position?
Vector<Double> thePixel(2);
thePixel(0) = 55.0;
@@ -190,10 +230,10 @@ ImageCoordinate build3Dcoords()
// Now we can make fruit of our labor - the ProjectedPosition itself.
ProjectedPosition myMapping(myMethod, myVectorType, myEpoch, myCoords,
myObs, myRot, thePixel, theBinning);
-
+
// we need to know what the units are of the measured value
MeasuredValue::Type myValueUnit(MeasuredValue::RADIO_VELOCITY);
- // we need to know what the above units are in reference
+ // we need to know what the above units are in reference
// here it is the velocity of the earth.
ReferenceValue myRefValue(ReferenceValue::VELOCITY, 9.56e+03);
// we need to know the value of the measurement
@@ -204,7 +244,7 @@ ImageCoordinate build3Dcoords()
Double myValuePos(27.3);
// Now we may construct the LinearAxis itself.
LinearAxis myLinearAxis(myValueUnit,myValue,myRefValue,myBin,myValuePos);
-
+
ImageCoordinate coords;
coords.addAxis(myMapping);
coords.addAxis(myLinearAxis);
@@ -215,20 +255,20 @@ ImageCoordinate build3Dcoords()
void createStandardImageOnDisk(const String &filename)
{
if(verbose_) cout<< "-- createStandardImageOnDisk --"<< endl;
-
+
ImageCoordinate coords(build3Dcoords());
IPosition imageShape(3, 50,47,31);
-
+
PagedImage<Float> standard(imageShape, coords, filename);
standard.set(TEST_PIXEL_VALUE);
-
+
if(verbose_) describeImage("standard image", standard);
}
// these are the constructors to test:
// PagedImage(const IPosition &shape, const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
-// PagedImage(const IPosition &shape, const Array<T> &array,
+// PagedImage(const IPosition &shape, const Array<T> &array,
// const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
// PagedImage(const String &filename, uInt rowNumber);
@@ -242,12 +282,12 @@ void testConstructors()
ImageCoordinate coords(build3Dcoords());
PagedImage<Float> image1(shape, coords, NAME0);
- if (verbose_)
+ if (verbose_)
describeImage("shape, coords, filename ctor: ", image1);
PagedImage<Float> image2(shape, coords, NAME1, True);
-
- if (verbose_)
+
+ if (verbose_)
describeImage("array, array, coords, filename, masking? ctor: ", image2);
PagedImage<Float> image3(STANDARD);
@@ -266,10 +306,10 @@ void testConstructors()
throw(AipsError("file-constructed image not ok"));
delete image4;
}
-
+
deleteFile(NAME0);
deleteFile(NAME1);
-}
+}
// test some of the data member manipulation functions
void testSetMemberFunctions()
@@ -279,7 +319,7 @@ void testSetMemberFunctions()
ImageCoordinate coords;
{
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
-
+
if(verbose_) describeImage("shape, coords, filename ctor: ", image6);
if(image6.ok())
throw(AipsError("shape-constructed image6 is ok, but shouldn't be!"));
@@ -288,20 +328,20 @@ void testSetMemberFunctions()
coords = build3Dcoords();
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
if(!image6.ok())
- throw(AipsError("shape-constructed image6 - is not ok!"));
+ throw(AipsError("shape-constructed image6 - is not ok!"));
image6.rename(NAME2);
if(verbose_) describeImage("image6, after rename", image6);
if(!image6.ok())
throw(AipsError("image6 after rename - not ok"));
-
+
image6.setCoordinateInfo(build3Dcoords());
if(verbose_) describeImage("image6, after setCoordinates", image6);
if(!image6.ok())
throw(AipsError("image6 after setCoordinates - not ok"));
}
-
+
deleteFile(NAME2);
-}
+}
//
void testArrayofImagesAndAssignmentOperator()
@@ -320,7 +360,7 @@ void testArrayofImagesAndAssignmentOperator()
filenames [7] = NAME7;
filenames [8] = NAME8;
filenames [9] = NAME9;
- {
+ {
Image<Float> images [max];
for(uInt i=0; i< max; i++) {
if(verbose_) cout<< "-- image array "<< i<< " --"<< endl;
@@ -328,14 +368,14 @@ void testArrayofImagesAndAssignmentOperator()
images [i] = Image<Float>(IPosition(2,10,10));
images [i].setCoordinateInfo(MinimalCoords());
images [i].setName(filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape constructed image, setCoords, setName in array",
images [i]);
} // if i is odd
else {
images [i] =
Image<Float>(IPosition(i+1,5), MinimalCoords(), filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape,coords,name constructed image in array,",
images [i]);
} // else
@@ -343,7 +383,7 @@ void testArrayofImagesAndAssignmentOperator()
throw(AipsError("one of array of images not ok"));
} // for i
} // scope
-
+
if(verbose_) cout<< "-- about to delete image files --"<< endl;
for(uInt i=0; i< max; i++)
@@ -351,3 +391,6 @@ void testArrayofImagesAndAssignmentOperator()
// * /
}
*/
+ //===============>>> Image_File <<<===============
+ }
+}
\ No newline at end of file
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h
index ec4f2c6d6bd..876236ca933 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h
@@ -24,22 +24,25 @@
#ifndef __IMAGER_IMAGER_FILE_H__
#define __IMAGER_IMAGER_FILE_H__
+#include <vector>
+#include <string>
+
#include<casa/aips.h>
#include <images/Images/PagedImage.h>
-///#include <trial/ImgCrdSys/ImageCoordinate.h>
+//#include <trial/ImgCrdSys/ImageCoordinate.h>
-#include <casa/Arrays/Vector.h>
+//#include <casa/Arrays/Vector.h>
#include <casa/Arrays/Cube.h>
-#include <casa/Exceptions/Error.h>
+//#include <casa/Exceptions/Error.h>
#include <casa/Arrays/IPosition.h>
-#include <tables/Tables.h>
-#include <casa/BasicSL/String.h>
+//#include <tables/Tables.h>
+//#include <casa/BasicSL/String.h>
-#include <casa/stdlib.h>
-#include <casa/iostream.h>
+//#include <casa/stdlib.h>
+//#include <casa/iostream.h>
-#include <casa/namespace.h>
-namespace LOFAR
+//#include <casa/namespace.h>
+namespace LOFAR
{
namespace CS1
{
@@ -48,14 +51,15 @@ namespace LOFAR
public:
Image_File(const std::string& imagename);
~Image_File();
-
+ void WriteImage(vector< casa::Cube<float> >* image);
+
protected:
- string ImageName;
- casa::PagedImage* Image;
- void init();
+ std::string ImageName;
+ casa::PagedImage<float>* Image;
+ void init(casa::IPosition shape, size_t length);
private:
}; // Image_File
- } //namespace CS1
+ }; //namespace CS1
}; //namespace LOFAR
/*
@@ -97,7 +101,7 @@ int main(int argc, char **argv)
cout << "OK" << endl;
} catch (AipsError x) {
cerr << "Caught Exception: " << x.getMesg() << endl;
- }
+ }
return 0;
}
@@ -153,7 +157,7 @@ void describeImage(const String &message, const PagedImage<Float> &image)
cout << " shape: "
<< image.coordinates().imageShape() << endl;
// * /
- cout << " ok: "<< image.ok()<< endl;
+ cout << " ok: "<< image.ok()<< endl;
}
// build a set of ImageCoordinates in three dimensions
@@ -173,7 +177,7 @@ ImageCoordinate build3Dcoords()
// Let's create an EarthPosition with full description of all parameters.
// We need a type - GEOCENTRIC seems good.
EarthPosition::Type theType(EarthPosition::GEOCENTRIC);
- // We need the time and date of the observation...
+ // We need the time and date of the observation...
Double julianDate = 2449376.0;
// and we can add on the UT.
julianDate += 16.52/24.0;
@@ -185,11 +189,11 @@ ImageCoordinate build3Dcoords()
ourPosition(1) = 34.1562394;
// geocentric radius (in meters) goes in the last field.
ourPosition(2) = 6372139.592;
- // then use these to build our EarthPosition.
+ // then use these to build our EarthPosition.
EarthPosition myObs(theType, julianDate, ourPosition);
// we need to know a rotation - here it is zero.
Double myRot = 0;
- // we need to know where the spherical position is to be on our 2-d
+ // we need to know where the spherical position is to be on our 2-d
// projection i.e. what pixel is associated with my object's position?
Vector<Double> thePixel(2);
thePixel(0) = 55.0;
@@ -202,10 +206,10 @@ ImageCoordinate build3Dcoords()
// Now we can make fruit of our labor - the ProjectedPosition itself.
ProjectedPosition myMapping(myMethod, myVectorType, myEpoch, myCoords,
myObs, myRot, thePixel, theBinning);
-
+
// we need to know what the units are of the measured value
MeasuredValue::Type myValueUnit(MeasuredValue::RADIO_VELOCITY);
- // we need to know what the above units are in reference
+ // we need to know what the above units are in reference
// here it is the velocity of the earth.
ReferenceValue myRefValue(ReferenceValue::VELOCITY, 9.56e+03);
// we need to know the value of the measurement
@@ -216,7 +220,7 @@ ImageCoordinate build3Dcoords()
Double myValuePos(27.3);
// Now we may construct the LinearAxis itself.
LinearAxis myLinearAxis(myValueUnit,myValue,myRefValue,myBin,myValuePos);
-
+
ImageCoordinate coords;
coords.addAxis(myMapping);
coords.addAxis(myLinearAxis);
@@ -227,20 +231,20 @@ ImageCoordinate build3Dcoords()
void createStandardImageOnDisk(const String &filename)
{
if(verbose_) cout<< "-- createStandardImageOnDisk --"<< endl;
-
+
ImageCoordinate coords(build3Dcoords());
IPosition imageShape(3, 50,47,31);
-
+
PagedImage<Float> standard(imageShape, coords, filename);
standard.set(TEST_PIXEL_VALUE);
-
+
if(verbose_) describeImage("standard image", standard);
}
// these are the constructors to test:
// PagedImage(const IPosition &shape, const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
-// PagedImage(const IPosition &shape, const Array<T> &array,
+// PagedImage(const IPosition &shape, const Array<T> &array,
// const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
// PagedImage(const String &filename, uInt rowNumber);
@@ -254,12 +258,12 @@ void testConstructors()
ImageCoordinate coords(build3Dcoords());
PagedImage<Float> image1(shape, coords, NAME0);
- if (verbose_)
+ if (verbose_)
describeImage("shape, coords, filename ctor: ", image1);
PagedImage<Float> image2(shape, coords, NAME1, True);
-
- if (verbose_)
+
+ if (verbose_)
describeImage("array, array, coords, filename, masking? ctor: ", image2);
PagedImage<Float> image3(STANDARD);
@@ -278,10 +282,10 @@ void testConstructors()
throw(AipsError("file-constructed image not ok"));
delete image4;
}
-
+
deleteFile(NAME0);
deleteFile(NAME1);
-}
+}
// test some of the data member manipulation functions
void testSetMemberFunctions()
@@ -291,7 +295,7 @@ void testSetMemberFunctions()
ImageCoordinate coords;
{
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
-
+
if(verbose_) describeImage("shape, coords, filename ctor: ", image6);
if(image6.ok())
throw(AipsError("shape-constructed image6 is ok, but shouldn't be!"));
@@ -300,20 +304,20 @@ void testSetMemberFunctions()
coords = build3Dcoords();
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
if(!image6.ok())
- throw(AipsError("shape-constructed image6 - is not ok!"));
+ throw(AipsError("shape-constructed image6 - is not ok!"));
image6.rename(NAME2);
if(verbose_) describeImage("image6, after rename", image6);
if(!image6.ok())
throw(AipsError("image6 after rename - not ok"));
-
+
image6.setCoordinateInfo(build3Dcoords());
if(verbose_) describeImage("image6, after setCoordinates", image6);
if(!image6.ok())
throw(AipsError("image6 after setCoordinates - not ok"));
}
-
+
deleteFile(NAME2);
-}
+}
//
void testArrayofImagesAndAssignmentOperator()
@@ -332,7 +336,7 @@ void testArrayofImagesAndAssignmentOperator()
filenames [7] = NAME7;
filenames [8] = NAME8;
filenames [9] = NAME9;
- {
+ {
Image<Float> images [max];
for(uInt i=0; i< max; i++) {
if(verbose_) cout<< "-- image array "<< i<< " --"<< endl;
@@ -340,14 +344,14 @@ void testArrayofImagesAndAssignmentOperator()
images [i] = Image<Float>(IPosition(2,10,10));
images [i].setCoordinateInfo(MinimalCoords());
images [i].setName(filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape constructed image, setCoords, setName in array",
images [i]);
} // if i is odd
else {
images [i] =
Image<Float>(IPosition(i+1,5), MinimalCoords(), filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape,coords,name constructed image in array,",
images [i]);
} // else
@@ -355,7 +359,7 @@ void testArrayofImagesAndAssignmentOperator()
throw(AipsError("one of array of images not ok"));
} // for i
} // scope
-
+
if(verbose_) cout<< "-- about to delete image files --"<< endl;
for(uInt i=0; i< max; i++)
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h~ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h~
index ab749b09e67..c3391201e61 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Image_File.h~
@@ -24,38 +24,42 @@
#ifndef __IMAGER_IMAGER_FILE_H__
#define __IMAGER_IMAGER_FILE_H__
+#include <vector>
+#include <string>
+
#include<casa/aips.h>
#include <images/Images/PagedImage.h>
-///#include <trial/ImgCrdSys/ImageCoordinate.h>
+//#include <trial/ImgCrdSys/ImageCoordinate.h>
-#include <casa/Arrays/Vector.h>
+//#include <casa/Arrays/Vector.h>
#include <casa/Arrays/Cube.h>
-#include <casa/Exceptions/Error.h>
+//#include <casa/Exceptions/Error.h>
#include <casa/Arrays/IPosition.h>
-#include <tables/Tables.h>
-#include <casa/BasicSL/String.h>
+//#include <tables/Tables.h>
+//#include <casa/BasicSL/String.h>
-#include <casa/stdlib.h>
-#include <casa/iostream.h>
+//#include <casa/stdlib.h>
+//#include <casa/iostream.h>
-#include <casa/namespace.h>
-namespace LOFAR
+//#include <casa/namespace.h>
+namespace LOFAR
{
namespace CS1
{
- class MS_File
+ class Image_File
{
public:
Image_File(const std::string& imagename);
~Image_File();
-
+ void WriteImage(vector< casa::Cube<float> >* image);
+
protected:
- string ImageName;
- casa::PagedImage* Image;
- void init();
+ std::string ImageName;
+ casa::PagedImage<float>* Image;
+ void init(casa::IPosition shape);
private:
}; // Image_File
- } //namespace CS1
+ }; //namespace CS1
}; //namespace LOFAR
/*
@@ -97,7 +101,7 @@ int main(int argc, char **argv)
cout << "OK" << endl;
} catch (AipsError x) {
cerr << "Caught Exception: " << x.getMesg() << endl;
- }
+ }
return 0;
}
@@ -153,7 +157,7 @@ void describeImage(const String &message, const PagedImage<Float> &image)
cout << " shape: "
<< image.coordinates().imageShape() << endl;
// * /
- cout << " ok: "<< image.ok()<< endl;
+ cout << " ok: "<< image.ok()<< endl;
}
// build a set of ImageCoordinates in three dimensions
@@ -173,7 +177,7 @@ ImageCoordinate build3Dcoords()
// Let's create an EarthPosition with full description of all parameters.
// We need a type - GEOCENTRIC seems good.
EarthPosition::Type theType(EarthPosition::GEOCENTRIC);
- // We need the time and date of the observation...
+ // We need the time and date of the observation...
Double julianDate = 2449376.0;
// and we can add on the UT.
julianDate += 16.52/24.0;
@@ -185,11 +189,11 @@ ImageCoordinate build3Dcoords()
ourPosition(1) = 34.1562394;
// geocentric radius (in meters) goes in the last field.
ourPosition(2) = 6372139.592;
- // then use these to build our EarthPosition.
+ // then use these to build our EarthPosition.
EarthPosition myObs(theType, julianDate, ourPosition);
// we need to know a rotation - here it is zero.
Double myRot = 0;
- // we need to know where the spherical position is to be on our 2-d
+ // we need to know where the spherical position is to be on our 2-d
// projection i.e. what pixel is associated with my object's position?
Vector<Double> thePixel(2);
thePixel(0) = 55.0;
@@ -202,10 +206,10 @@ ImageCoordinate build3Dcoords()
// Now we can make fruit of our labor - the ProjectedPosition itself.
ProjectedPosition myMapping(myMethod, myVectorType, myEpoch, myCoords,
myObs, myRot, thePixel, theBinning);
-
+
// we need to know what the units are of the measured value
MeasuredValue::Type myValueUnit(MeasuredValue::RADIO_VELOCITY);
- // we need to know what the above units are in reference
+ // we need to know what the above units are in reference
// here it is the velocity of the earth.
ReferenceValue myRefValue(ReferenceValue::VELOCITY, 9.56e+03);
// we need to know the value of the measurement
@@ -216,7 +220,7 @@ ImageCoordinate build3Dcoords()
Double myValuePos(27.3);
// Now we may construct the LinearAxis itself.
LinearAxis myLinearAxis(myValueUnit,myValue,myRefValue,myBin,myValuePos);
-
+
ImageCoordinate coords;
coords.addAxis(myMapping);
coords.addAxis(myLinearAxis);
@@ -227,20 +231,20 @@ ImageCoordinate build3Dcoords()
void createStandardImageOnDisk(const String &filename)
{
if(verbose_) cout<< "-- createStandardImageOnDisk --"<< endl;
-
+
ImageCoordinate coords(build3Dcoords());
IPosition imageShape(3, 50,47,31);
-
+
PagedImage<Float> standard(imageShape, coords, filename);
standard.set(TEST_PIXEL_VALUE);
-
+
if(verbose_) describeImage("standard image", standard);
}
// these are the constructors to test:
// PagedImage(const IPosition &shape, const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
-// PagedImage(const IPosition &shape, const Array<T> &array,
+// PagedImage(const IPosition &shape, const Array<T> &array,
// const MinimalCoords &coordinateInfo,
// const String &nameOfNewFile, uInt rowNumber);
// PagedImage(const String &filename, uInt rowNumber);
@@ -254,12 +258,12 @@ void testConstructors()
ImageCoordinate coords(build3Dcoords());
PagedImage<Float> image1(shape, coords, NAME0);
- if (verbose_)
+ if (verbose_)
describeImage("shape, coords, filename ctor: ", image1);
PagedImage<Float> image2(shape, coords, NAME1, True);
-
- if (verbose_)
+
+ if (verbose_)
describeImage("array, array, coords, filename, masking? ctor: ", image2);
PagedImage<Float> image3(STANDARD);
@@ -278,10 +282,10 @@ void testConstructors()
throw(AipsError("file-constructed image not ok"));
delete image4;
}
-
+
deleteFile(NAME0);
deleteFile(NAME1);
-}
+}
// test some of the data member manipulation functions
void testSetMemberFunctions()
@@ -291,7 +295,7 @@ void testSetMemberFunctions()
ImageCoordinate coords;
{
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
-
+
if(verbose_) describeImage("shape, coords, filename ctor: ", image6);
if(image6.ok())
throw(AipsError("shape-constructed image6 is ok, but shouldn't be!"));
@@ -300,20 +304,20 @@ void testSetMemberFunctions()
coords = build3Dcoords();
PagedImage<Float> image6(IPosition(3,10,10,4), coords, NAME2);
if(!image6.ok())
- throw(AipsError("shape-constructed image6 - is not ok!"));
+ throw(AipsError("shape-constructed image6 - is not ok!"));
image6.rename(NAME2);
if(verbose_) describeImage("image6, after rename", image6);
if(!image6.ok())
throw(AipsError("image6 after rename - not ok"));
-
+
image6.setCoordinateInfo(build3Dcoords());
if(verbose_) describeImage("image6, after setCoordinates", image6);
if(!image6.ok())
throw(AipsError("image6 after setCoordinates - not ok"));
}
-
+
deleteFile(NAME2);
-}
+}
//
void testArrayofImagesAndAssignmentOperator()
@@ -332,7 +336,7 @@ void testArrayofImagesAndAssignmentOperator()
filenames [7] = NAME7;
filenames [8] = NAME8;
filenames [9] = NAME9;
- {
+ {
Image<Float> images [max];
for(uInt i=0; i< max; i++) {
if(verbose_) cout<< "-- image array "<< i<< " --"<< endl;
@@ -340,14 +344,14 @@ void testArrayofImagesAndAssignmentOperator()
images [i] = Image<Float>(IPosition(2,10,10));
images [i].setCoordinateInfo(MinimalCoords());
images [i].setName(filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape constructed image, setCoords, setName in array",
images [i]);
} // if i is odd
else {
images [i] =
Image<Float>(IPosition(i+1,5), MinimalCoords(), filenames [i]);
- if(verbose_)
+ if(verbose_)
describeImage("shape,coords,name constructed image in array,",
images [i]);
} // else
@@ -355,7 +359,7 @@ void testArrayofImagesAndAssignmentOperator()
throw(AipsError("one of array of images not ok"));
} // for i
} // scope
-
+
if(verbose_) cout<< "-- about to delete image files --"<< endl;
for(uInt i=0; i< max; i++)
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc b/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc
index 19535e217ac..cd9cb1e80b1 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -31,7 +31,7 @@
#define IMAGER_VERSION "0.10"
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
@@ -39,7 +39,7 @@ namespace LOFAR
ImagerProcessControl::ImagerProcessControl()
: ProcessControl()
{
- ItsMS = NULL;
+ itsMS = NULL;
itsImage = NULL;
itsImager = NULL;
}
@@ -62,8 +62,8 @@ namespace LOFAR
tribool ImagerProcessControl::define()
{
LOFAR::ACC::APS::ParameterSet* ParamSet = LOFAR::ACC::APS::globalParameterSet();
- itsMS = ParamSet->getString("ms");
- itsImager = ParamSet->getString("image");
+ itsMSName = ParamSet->getString("ms");
+ itsImageName = ParamSet->getString("image");
itsResolution = ParamSet->getInt32("resolution");
return true;
}
@@ -73,9 +73,7 @@ namespace LOFAR
{
try{
std::cout << "Runnning Imager please wait..." << std::endl;
- itsImager->MakeImage(itsMS,
- itsImage,
- itsResolution);
+ itsImager->MakeImage(itsResolution);
}
catch(casa::AipsError& err)
{
@@ -97,8 +95,9 @@ namespace LOFAR
string("This is experimental software, please report errors or requests to renting@astron.nl\n") +
string("Documentation can be found at: www.astron.nl/~renting\n");
cout << itsMS << endl;
- itsMS = new WSRT::MS_File(itsMS);
- itsImager = new DFTImager (itsMS, itsResolution);
+ itsMS = new MS_File(itsMSName);
+ itsImage = new Image_File(itsImageName);
+ itsImager = new DFTImager (itsMS, itsImage);
}
catch(casa::AipsError& err)
{
@@ -115,7 +114,7 @@ namespace LOFAR
//===============>>> ImagerProcessControl::quit <<<================================
tribool ImagerProcessControl::quit()
- {
+ {
if (itsMS)
{
delete itsMS;
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc~ b/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc~
index 4e6d97248a8..026dc2ae71f 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/ImagerProcessControl.cc~
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -31,7 +31,7 @@
#define IMAGER_VERSION "0.10"
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
@@ -39,7 +39,7 @@ namespace LOFAR
ImagerProcessControl::ImagerProcessControl()
: ProcessControl()
{
- ItsMS = NULL;
+ itsMS = NULL;
itsImage = NULL;
itsImager = NULL;
}
@@ -60,22 +60,20 @@ namespace LOFAR
//===============>>> ImagerProcessControl::define <<<==============================
tribool ImagerProcessControl::define()
- {
+ {
LOFAR::ACC::APS::ParameterSet* ParamSet = LOFAR::ACC::APS::globalParameterSet();
- itsMS = ParamSet->getString("ms");
- itsImager = ParamSet->getString("image");
+ itsMSName = ParamSet->getString("ms");
+ itsImageName = ParamSet->getString("image");
itsResolution = ParamSet->getInt32("resolution");
return true;
}
//===============>>> ImagerProcessControl::run <<<=================================
tribool ImagerProcessControl::run()
- {
+ {
try{
std::cout << "Runnning Imager please wait..." << std::endl;
- itsImager->MakeImage(itsMS,
- itsImage,
- itsResolution);
+ itsImager->MakeImage();
}
catch(casa::AipsError& err)
{
@@ -97,8 +95,9 @@ namespace LOFAR
string("This is experimental software, please report errors or requests to renting@astron.nl\n") +
string("Documentation can be found at: www.astron.nl/~renting\n");
cout << itsMS << endl;
- itsMS = new WSRT::MS_File(itsMS);
- itsImager = new DFTImager (itsMS, itsResolution);
+ itsMS = new MS_File(itsMSName);
+ itsImage = new Image_File(itsImageName);
+ itsImager = new DFTImager (itsMS, itsImage, itsResolution);
}
catch(casa::AipsError& err)
{
@@ -115,7 +114,7 @@ namespace LOFAR
//===============>>> ImagerProcessControl::quit <<<================================
tribool ImagerProcessControl::quit()
- {
+ {
if (itsMS)
{
delete itsMS;
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc
index 0ba609c0605..60272f89f46 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -28,55 +28,55 @@
using namespace casa;
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
-
+
//===============>>> MS_File::MS_File <<<===============
-
+
MS_File::MS_File(const std::string& msname)
{
MSName = msname;
MS = new MeasurementSet(MSName, Table::Update);
init();
}
-
+
//===============>>> MS_File::~MS_File <<<===============
-
+
MS_File::~MS_File()
{
delete MS;
}
-
+
//===============>>> MS_File::init <<<===============
-
+
void MS_File::init()
{
//Number of samples
- itsNumSamples = (*MS).nrow();
+ itsNumSamples = (*MS).nrow();
std::cout << "NumSamples" << itsNumSamples << std::endl;
//Number of Fields
MSField fields = (*MS).field();
itsNumFields = fields.nrow();
std::cout << "NumFields" << itsNumFields << std::endl;
-
- //Number of Antennae
+
+ //Number of Antennae
MSAntenna antennae = (*MS).antenna();
itsNumAntennae = antennae.nrow();
std::cout << "NumAntennae" << itsNumAntennae << std::endl;
-
+
//Antenna Names
ROScalarColumn<String> ANT_NAME_col(antennae, "NAME");
Vector<String> ant_names = ANT_NAME_col.getColumn();
ant_names.tovector(itsAntennaNames);
-
+
//Number of channels in the Band
MSSpectralWindow spectral_window = (*MS).spectralWindow();
ROScalarColumn<Int> NUM_CHAN_col(spectral_window, "NUM_CHAN");
itsNumChannels = NUM_CHAN_col(0);
std::cout << "NumChannels" << itsNumChannels << std::endl;
-
+
//Number of polarizations
MSPolarization polarization = (*MS).polarization();
ROScalarColumn<Int> NUM_CORR_col(polarization, "NUM_CORR");
@@ -85,69 +85,57 @@ namespace LOFAR
itsPolarizations.resize(itsNumPolarizations);
CORR_TYPE_col.get(0, itsPolarizations);
std::cout << "NumPolarizations" << itsNumPolarizations << std::endl;
-
+
//calculate theoretical noise level
ROScalarColumn<Double> EXPOSURE_col((*MS), "EXPOSURE");
Double exposure = EXPOSURE_col(0);
-
+
ROScalarColumn<Double> TOTAL_BANDWIDTH_col(spectral_window, "TOTAL_BANDWIDTH");
Double bandwidth = TOTAL_BANDWIDTH_col(0) / itsNumChannels;
-
+
itsNoiseLevel = 1.0 / sqrt(bandwidth * exposure);
std::cout << "Noiselevel" << itsNoiseLevel << std::endl;
-
+
//calculate number of timeslots
ROScalarColumn<Double> INTERVAL_col((*MS), "INTERVAL");
Double interval = INTERVAL_col(0);
-
+
//Number of timeslots
ROScalarColumn<Double> TIME_CENTROID_col((*MS), "TIME_CENTROID");
Double firstdate = TIME_CENTROID_col(0);
Double lastdate = TIME_CENTROID_col(itsNumSamples-1);
std::cout << "interval" << interval << std::endl;
-
+
itsNumTimeslots = (int)((lastdate-firstdate)/interval) + 1;
std::cout << "Numtimeslots" << itsNumTimeslots << std::endl;
-
+
//calculate number of baselines.
itsNumPairs = (itsNumAntennae) * (itsNumAntennae + 1) / 2; //Triangular numbers formula
std::cout << "NumPairs" << itsNumPairs << std::endl;
-
+
//calculate number of Bands
itsNumBands = itsNumSamples / (itsNumPairs * itsNumTimeslots);
std::cout << "NumBands" << itsNumBands << std::endl;
-
+
}
-
+
//===============>>> MS_File::antennas <<<===============
-
- MSAntenna MS_File::antenna()
- {
- return (*MS).antennas();
+
+ MSAntenna MS_File::antennas()
+ {
+ return (*MS).antenna();
}
-
+
//===============>>> MS_File::BaselineIterator <<<===============
-
+ /* Returns a table of records for one timeslot at a time*/
TableIterator MS_File::TimeslotIterator()
- {
+ {
Block<String> ms_iteration_variables(1);
ms_iteration_variables[0] = "TIME_CENTROID";
-
- return TableIterator((*MS), ms_iteration_variables);
- }
-
- //===============>>> MS_File::BaselineIterator <<<===============
-
- TableIterator MS_File::TimeAntennaIterator()
- {
- Block<String> ms_iteration_variables(4);
- ms_iteration_variables[0] = "TIME_CENTROID";
- ms_iteration_variables[1] = "DATA_DESC_ID";
- ms_iteration_variables[2] = "ANTENNA1";
- ms_iteration_variables[3] = "ANTENNA2";
-
+
return TableIterator((*MS), ms_iteration_variables);
- }
+ }
+
//===============>>> MS_File <<<===============
}
-}
\ No newline at end of file
+}
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc~ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc~
index 214db27983c..c00f61aed4b 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.cc~
@@ -1,5 +1,5 @@
/***************************************************************************
- * Copyright (C) 2006 by ASTRON, Adriaan Renting *
+ * Copyright (C) 2007 by ASTRON, Adriaan Renting *
* renting@astron.nl *
* *
* This program is free software; you can redistribute it and/or modify *
@@ -28,55 +28,55 @@
using namespace casa;
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
-
+
//===============>>> MS_File::MS_File <<<===============
-
+
MS_File::MS_File(const std::string& msname)
{
MSName = msname;
MS = new MeasurementSet(MSName, Table::Update);
init();
}
-
+
//===============>>> MS_File::~MS_File <<<===============
-
+
MS_File::~MS_File()
{
delete MS;
}
-
+
//===============>>> MS_File::init <<<===============
-
+
void MS_File::init()
{
//Number of samples
- itsNumSamples = (*MS).nrow();
+ itsNumSamples = (*MS).nrow();
std::cout << "NumSamples" << itsNumSamples << std::endl;
//Number of Fields
MSField fields = (*MS).field();
itsNumFields = fields.nrow();
std::cout << "NumFields" << itsNumFields << std::endl;
-
- //Number of Antennae
+
+ //Number of Antennae
MSAntenna antennae = (*MS).antenna();
itsNumAntennae = antennae.nrow();
std::cout << "NumAntennae" << itsNumAntennae << std::endl;
-
+
//Antenna Names
ROScalarColumn<String> ANT_NAME_col(antennae, "NAME");
Vector<String> ant_names = ANT_NAME_col.getColumn();
ant_names.tovector(itsAntennaNames);
-
+
//Number of channels in the Band
MSSpectralWindow spectral_window = (*MS).spectralWindow();
ROScalarColumn<Int> NUM_CHAN_col(spectral_window, "NUM_CHAN");
itsNumChannels = NUM_CHAN_col(0);
std::cout << "NumChannels" << itsNumChannels << std::endl;
-
+
//Number of polarizations
MSPolarization polarization = (*MS).polarization();
ROScalarColumn<Int> NUM_CORR_col(polarization, "NUM_CORR");
@@ -85,69 +85,57 @@ namespace LOFAR
itsPolarizations.resize(itsNumPolarizations);
CORR_TYPE_col.get(0, itsPolarizations);
std::cout << "NumPolarizations" << itsNumPolarizations << std::endl;
-
+
//calculate theoretical noise level
ROScalarColumn<Double> EXPOSURE_col((*MS), "EXPOSURE");
Double exposure = EXPOSURE_col(0);
-
+
ROScalarColumn<Double> TOTAL_BANDWIDTH_col(spectral_window, "TOTAL_BANDWIDTH");
Double bandwidth = TOTAL_BANDWIDTH_col(0) / itsNumChannels;
-
+
itsNoiseLevel = 1.0 / sqrt(bandwidth * exposure);
std::cout << "Noiselevel" << itsNoiseLevel << std::endl;
-
+
//calculate number of timeslots
ROScalarColumn<Double> INTERVAL_col((*MS), "INTERVAL");
Double interval = INTERVAL_col(0);
-
+
//Number of timeslots
ROScalarColumn<Double> TIME_CENTROID_col((*MS), "TIME_CENTROID");
Double firstdate = TIME_CENTROID_col(0);
Double lastdate = TIME_CENTROID_col(itsNumSamples-1);
std::cout << "interval" << interval << std::endl;
-
+
itsNumTimeslots = (int)((lastdate-firstdate)/interval) + 1;
std::cout << "Numtimeslots" << itsNumTimeslots << std::endl;
-
+
//calculate number of baselines.
itsNumPairs = (itsNumAntennae) * (itsNumAntennae + 1) / 2; //Triangular numbers formula
std::cout << "NumPairs" << itsNumPairs << std::endl;
-
+
//calculate number of Bands
itsNumBands = itsNumSamples / (itsNumPairs * itsNumTimeslots);
std::cout << "NumBands" << itsNumBands << std::endl;
-
+
}
-
- //===============>>> MS_File::BaselineIterator <<<===============
-
- MSAntenna MS_File::antenna()
- {
+
+ //===============>>> MS_File::antennas <<<===============
+
+ MSAntenna MS_File::antennas()
+ {
return (*MS).antenna();
}
-
+
//===============>>> MS_File::BaselineIterator <<<===============
-
+ /* Returns one timeslot at a time*/
TableIterator MS_File::TimeslotIterator()
- {
+ {
Block<String> ms_iteration_variables(1);
ms_iteration_variables[0] = "TIME_CENTROID";
-
- return TableIterator((*MS), ms_iteration_variables);
- }
-
- //===============>>> MS_File::BaselineIterator <<<===============
-
- TableIterator MS_File::TimeAntennaIterator()
- {
- Block<String> ms_iteration_variables(4);
- ms_iteration_variables[0] = "TIME_CENTROID";
- ms_iteration_variables[1] = "DATA_DESC_ID";
- ms_iteration_variables[2] = "ANTENNA1";
- ms_iteration_variables[3] = "ANTENNA2";
-
+
return TableIterator((*MS), ms_iteration_variables);
- }
+ }
+
//===============>>> MS_File <<<===============
}
-}
\ No newline at end of file
+}
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h
index 24129773290..c1a530418cc 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h
@@ -26,7 +26,7 @@
#include <tables/Tables.h>
#include <tables/Tables/TableIter.h>
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
@@ -35,7 +35,7 @@ namespace LOFAR
public:
MS_File(const std::string& msname);
~MS_File();
-
+
int itsNumSamples;
int itsNumAntennae;
int itsNumFields;
@@ -47,18 +47,17 @@ namespace LOFAR
double itsNoiseLevel;
std::vector<casa::String> itsAntennaNames;
casa::Vector<casa::Int> itsPolarizations;
-
+
casa::TableIterator TimeslotIterator();
- casa::TableIterator TimeAntennaIterator();
casa::MSAntenna antennas();
-
+
protected:
string MSName;
casa::MeasurementSet* MS;
void init();
private:
}; // MS_File
- } //namespace CS1
+ }; //namespace CS1
}; //namespace LOFAR
#endif // __IMAGER_MS_FILE_H__
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h~ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h~
index e7eb09c9c7c..55b1c7469e6 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h~
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/MS_File.h~
@@ -26,7 +26,7 @@
#include <tables/Tables.h>
#include <tables/Tables/TableIter.h>
-namespace LOFAR
+namespace LOFAR
{
namespace CS1
{
@@ -35,7 +35,7 @@ namespace LOFAR
public:
MS_File(const std::string& msname);
~MS_File();
-
+
int itsNumSamples;
int itsNumAntennae;
int itsNumFields;
@@ -47,18 +47,17 @@ namespace LOFAR
double itsNoiseLevel;
std::vector<casa::String> itsAntennaNames;
casa::Vector<casa::Int> itsPolarizations;
-
+
casa::TableIterator TimeslotIterator();
- casa::TableIterator TimeAntennaIterator();
casa::MSAntenna antennas();
-
+
protected:
string MSName;
casa::MeasurementSet* MS;
void init();
private:
}; // MS_File
- } //namespace CS1_Flagger
+ } //namespace CS1
}; //namespace LOFAR
#endif // __IMAGER_MS_FILE_H__
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.am b/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.am
index 2c4377bfd77..0708a9188a7 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.am
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.am
@@ -5,6 +5,8 @@ CS1_DFTImager_SOURCES = Imager.cc \
DFTImager.h \
MS_File.cc \
MS_File.h \
+ Image_File.cc \
+ Image_File.h \
ImagerProcessControl.cc
CS1_DFTImager_LDADD =
CS1_DFTImager_DEPENDENCIES = $(LOFAR_DEPEND)
diff --git a/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.in b/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.in
index 5b00dacb3a6..452d1c49b62 100644
--- a/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.in
+++ b/Appl/CEP/CS1/CS1_DFTImager/src/Makefile.in
@@ -63,7 +63,8 @@ am__installdirs = "$(DESTDIR)$(bindir)"
binPROGRAMS_INSTALL = $(INSTALL_PROGRAM)
PROGRAMS = $(bin_PROGRAMS)
am_CS1_DFTImager_OBJECTS = Imager.$(OBJEXT) DFTImager.$(OBJEXT) \
- MS_File.$(OBJEXT) ImagerProcessControl.$(OBJEXT)
+ MS_File.$(OBJEXT) Image_File.$(OBJEXT) \
+ ImagerProcessControl.$(OBJEXT)
CS1_DFTImager_OBJECTS = $(am_CS1_DFTImager_OBJECTS)
DEFAULT_INCLUDES = -I. -I$(srcdir) -I$(top_builddir)
depcomp = $(SHELL) $(top_srcdir)/depcomp
@@ -231,6 +232,8 @@ CS1_DFTImager_SOURCES = Imager.cc \
DFTImager.h \
MS_File.cc \
MS_File.h \
+ Image_File.cc \
+ Image_File.h \
ImagerProcessControl.cc
CS1_DFTImager_LDADD =
@@ -339,6 +342,7 @@ distclean-compile:
-rm -f *.tab.c
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/DFTImager.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Image_File.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/Imager.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ImagerProcessControl.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/MS_File.Po@am__quote@
--
GitLab