bug 1436:

Created a W-projection gridder that applies the corrections

CEP/Imager/LofarGridder/CMakeLists.txt

 # $Id$
 
-lofar_package(LofarGridder 0.1 DEPENDS MWImager)
+#lofar_package(LofarGridder 0.1 DEPENDS MWImager)
+#lofar_package(LofarGridder 0.1 DEPENDS BBSKernel)
 
 lofar_find_package(AskapSoft REQUIRED)
-lofar_find_package(Casacore COMPONENTS images REQUIRED)
+lofar_find_package(LOFAR COMPONENTS bbskernel parmdb mwcommon blob common REQUIRED)
+lofar_find_package(Casacore COMPONENTS ms images REQUIRED)
+lofar_find_package(Boost REQUIRED)
 
 add_subdirectory(include/LofarGridder)
 add_subdirectory(src)

CEP/Imager/LofarGridder/include/LofarGridder/LofarGridder.h

-//# LofarGridder.h: Test visibility gridder.
+//# LofarGridder.h: Gridder for LOFAR data correcting for DD effects
 //#
 //# Copyright (C) 2009
 //# ASTRON (Netherlands Institute for Radio Astronomy)
@@ -22,23 +22,26 @@
 //#
 //# @author Ger van Diepen <diepen at astron dot nl>
 
-#ifndef LOFAR_LOFARGRIDDER_TESTGRIDDER_H
-#define LOFAR_LOFARGRIDDER_TESTGRIDDER_H
+#ifndef LOFAR_LOFARGRIDDER_LOFARGRIDDER_H
+#define LOFAR_LOFARGRIDDER_LOFARGRIDDER_H
 
 #include <gridding/TableVisGridder.h>
 #include <Common/ParameterSet.h>
+#include <Common/lofar_vector.h>
+#include <Common/lofar_string.h>
 
 namespace LOFAR
 {
-  // @brief Gridder to test dynamic loading of gridders
+  // @brief Gridder for LOFAR data correcting for DD effects
   //
   // @ingroup testgridder
-  class LofarGridder : public askap::synthesis::TableVisGridder
+
+  class LofarGridder : public askap::synthesis::IVisGridder
   {
   public:
 
-    // Standard two dimensional box gridding
-    LofarGridder();
+    // Construct from the given parset.
+    explicit LofarGridder (const ParameterSet&);
 
     // Clone this Gridder
     virtual askap::synthesis::IVisGridder::ShPtr clone();
@@ -55,19 +58,89 @@ namespace LOFAR
     // @brief Register the gridder create function with its name.
     static void registerGridder();
 
+    // @brief Initialise the gridding
+    // @param axes axes specifications
+    // @param shape Shape of output image: u,v,pol,chan
+    // @param dopsf Make the psf?
+    virtual void initialiseGrid (const askap::scimath::Axes& axes,
+                                 const casa::IPosition& shape,
+                                 const bool dopsf=true);
+
+    // @brief Grid the visibility data.
+    // @param acc const data accessor to work with
+    // @note a non-const adapter is created behind the scene. If no on-the-fly
+    // visibility correction is performed, this adapter is equivalent to the
+    // original const data accessor.
+    virtual void grid (askap::synthesis::IConstDataAccessor& acc);
+      
+    // Form the final output image
+    // @param out Output double precision image or PSF
+    virtual void finaliseGrid (casa::Array<double>& out);
+
+    // @brief Calculate weights image
+    // @details Form the sum of the convolution function squared, 
+    // multiplied by the weights for each different convolution 
+    // function. This is used in the evaluation of the position
+    // dependent sensitivity
+    // @param out Output double precision sum of weights images
+    virtual void finaliseWeights (casa::Array<double>& out);
+
+    // @brief Initialise the degridding
+    // @param axes axes specifications
+    // @param image Input image: cube: u,v,pol,chan
+    virtual void initialiseDegrid (const askap::scimath::Axes& axes,
+                                   const casa::Array<double>& image);
+
+    // @brief Make context-dependant changes to the gridder behaviour
+    // @param context context
+    virtual void customiseForContext (casa::String context);
+      
+    // @brief assign weights
+    // @param viswt shared pointer to visibility weights
+    virtual void initVisWeights (askap::synthesis::IVisWeights::ShPtr viswt);
+      
+    // @brief Degrid the visibility data.
+    // @param[in] acc non-const data accessor to work with  
+    virtual void degrid (askap::synthesis::IDataAccessor& acc);
+
+    // @brief Finalise degridding.
+    virtual void finaliseDegrid();
+
     // @brief Initialise the indices
     // @param[in] acc const data accessor to work with
-    virtual void initIndices(const askap::synthesis::IConstDataAccessor& acc);
+    ///    virtual void initIndices (const askap::synthesis::IConstDataAccessor& acc);
 
     // @brief Correct for gridding convolution function
     // @param image image to be corrected
-    virtual void correctConvolution(casa::Array<double>& image);
+    ///    virtual void correctConvolution(casa::Array<double>& image);
 				
   protected:
     // Initialize convolution function
     // @param[in] acc const data accessor to work with
-    virtual void initConvolutionFunction
-    (const askap::synthesis::IConstDataAccessor& acc);
+    ///    virtual void initConvolutionFunction
+    ///    (const askap::synthesis::IConstDataAccessor& acc);
+
+  private:
+    // Initialize the corrections.
+    void initCorrections (const askap::synthesis::IConstDataAccessor& acc);
+
+
+    //# Data members.
+    // Is the gridder correction initialized?
+    bool itsInitialized;
+    // Nr of frequency channels to average.
+    uint itsFreqAvg;
+    // Nr of time slots to average.
+    uint itsTimeAvg;
+    // The effects to correct for.
+    vector<string> itsCorrect;
+    // The ParmDB to use.
+    string itsParmDBName;
+    // The gridder to use after correcting and averaging.
+    //# Unfortunatelty we need to use TableVisGridder because the
+    //# function getImageCntre is available only there.
+    askap::synthesis::IVisGridder::ShPtr itsIGridder;
+    askap::synthesis::TableVisGridder* itsGridder;
   };
 
 } //# end namespace

CEP/Imager/LofarGridder/include/LofarGridder/LofarWProjectGridder.h

+//# LofarWProjectGridder.h: Gridder for LOFAR data correcting for DD effects
+//#
+//# Copyright (C) 2009
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+//#
+//# @author Ger van Diepen <diepen at astron dot nl>
+
+#ifndef LOFAR_LOFARGRIDDER_LOFARWPROJECTGRIDDER_H
+#define LOFAR_LOFARGRIDDER_LOFARWPROJECTGRIDDER_H
+
+//# LOFAR includes
+#include <BBSKernel/StationResponse.h>
+#include <Common/ParameterSet.h>
+#include <Common/lofar_vector.h>
+#include <Common/lofar_string.h>
+
+//# askap includes
+#include <gridding/WProjectVisGridder.h>
+
+//# casacore includes
+#include <ms/MeasurementSets/MeasurementSet.h>
+#include <casa/Quanta/MVDirection.h>
+
+namespace LOFAR
+{
+  // @brief Gridder for LOFAR data correcting for DD effects
+  //
+  // @ingroup testgridder
+
+  class LofarWProjectGridder : public askap::synthesis::WProjectVisGridder
+  {
+  public:
+
+    // Construct from the given parset.
+    explicit LofarWProjectGridder (const ParameterSet&);
+
+    // Clone this Gridder
+    virtual askap::synthesis::IVisGridder::ShPtr clone();
+
+    virtual ~LofarWProjectGridder();
+
+    // @brief Function to create the gridder from a parset.
+    // This function will be registered in the gridder registry.
+    static askap::synthesis::IVisGridder::ShPtr makeGridder
+    (const ParameterSet&);
+
+    // @brief Return the (unique) name of the gridder.
+    static const std::string& gridderName();
+
+    // @brief Register the gridder create function with its name.
+    static void registerGridder();
+
+  private:
+    // Correct the visibilities before gridding or after degridding.
+    void correctVisibilities (askap::synthesis::IDataAccessor& acc,
+                              bool forward);
+
+    // Initialize the corrections.
+    void initCorrections (const askap::synthesis::IConstDataAccessor& acc);
+
+    // Form an Instrument object.
+    BBS::Instrument makeInstrument (const casa::MeasurementSet& ms);
+
+
+    //# Data members.
+    // Is the gridder correction initialized?
+    bool itsInitialized;
+    // The ParmDB to use.
+    string itsParmDBName;
+    // The config file to use.
+    string itsConfigName;
+    string itsConfigPath;
+    // The StationResponse object to calculate the corrections.
+    BBS::StationResponse::Ptr itsResponse;
+    // Time and frequency interval in MS.
+    double itsTimeInterval;
+    double itsFreqInterval;
+    double itsStartFreq;
+    // The last facet center used.
+    casa::MVDirection itsLastFacetCenter;
+  };
+
+} //# end namespace
+
+#endif

CEP/Imager/LofarGridder/src/CMakeLists.txt

@@ -3,7 +3,7 @@
 include(LofarPackageVersion)
 
 lofar_add_library(lofargridder
-  Package__Version.cc Register.cc LofarGridder.cc
+  Package__Version.cc Register.cc LofarWProjectGridder.cc
   )
 
 lofar_add_bin_program(versionlofargridder versionlofargridder.cc)

CEP/Imager/LofarGridder/src/LofarGridder.cc

-//# TestGridder.cc: Test visibility gridder.
+//# LofarGridder.cc: Gridder for LOFAR data correcting for DD effects
 //#
 //# Copyright (C) 2009
 //# ASTRON (Netherlands Institute for Radio Astronomy)
@@ -22,17 +22,39 @@
 //#
 //# @author Ger van Diepen <diepen at astron dot nl>
 
+#include <lofar_config.h>
+#include <ParmDB/Grid.h>
 
+//# ASKAP includes
 #include <LofarGridder/LofarGridder.h>
 #include <gridding/VisGridderFactory.h>
+#include <dataaccess/TableConstDataAccessor.h>
+#include <dataaccess/TableConstDataIterator.h>
+#include <dataaccess/OnDemandBufferDataAccessor.h>
 
+//#casacore includes
+#include <tables/Tables/TableRecord.h>
+
+using namespace LOFAR::BBS;
+using namespace casa;
+using namespace askap;
 using namespace askap::synthesis;
 
 namespace LOFAR
 {
 
-  LofarGridder::LofarGridder()
-  {}
+  LofarGridder::LofarGridder (const ParameterSet& parset)
+    : itsInitialized (false)
+  {
+    itsTimeAvg = parset.getInt ("average.timestep", 1);
+    itsFreqAvg = parset.getInt ("average.freqstep", 1);
+    itsCorrect = parset.getStringVector ("correct");
+    itsParmDBName = parset.getString ("correct.parmdb");
+    string gridderName = parset.getString ("name");
+    itsIGridder = VisGridderFactory::createGridder (gridderName, parset);
+    itsGridder = dynamic_cast<TableVisGridder*>(itsIGridder.get());
+    ASSERT (itsGridder != 0);
+  }
 
   LofarGridder::~LofarGridder()
   {}
@@ -43,10 +65,10 @@ namespace LOFAR
     return IVisGridder::ShPtr(new LofarGridder(*this));
   }
 
-  IVisGridder::ShPtr LofarGridder::makeGridder (const ParameterSet&)
+  IVisGridder::ShPtr LofarGridder::makeGridder (const ParameterSet& parset)
   {
     std::cout << "LofarGridder::makeGridder" << std::endl;
-    return IVisGridder::ShPtr(new LofarGridder());
+    return IVisGridder::ShPtr(new LofarGridder(parset));
   }
 
   const std::string& LofarGridder::gridderName()
@@ -60,24 +82,103 @@ namespace LOFAR
     VisGridderFactory::registerGridder (gridderName(), &makeGridder);
   }
 
-  void LofarGridder::initIndices(const IConstDataAccessor&) 
+  void LofarGridder::initialiseGrid(const scimath::Axes& axes,
+                                    const casa::IPosition& shape,
+                                    const bool dopsf)
+  {
+    itsGridder->initialiseGrid (axes, shape, dopsf);
+  }
+
+  void LofarGridder::grid(IConstDataAccessor& acc)
+  {
+    if (! itsInitialized) {
+      initCorrections (acc);
+    }
+    // do correction, uvw rotation, and averaging here
+    // For time being only correction is done.
+    // It is possible to do the following:
+    ///  Array<Complex> data (acc.visibility());
+    // which makes a reference copy of the array, thus changs it later.
+    // This is an optimization we won't make yet.
+    // Can use PolConverter::isLinear() to check if data are XX,XY,YX,YY
+    OnDemandBufferDataAccessor acc2(acc);
+    Array<Complex>& data = acc2.rwVisibility();
+    itsGridder->grid (acc2);
+    // Get facet center in J2000.
+    casa::MVDirection center = itsGridder->getImageCentre();
+    double time = acc.time();    // time in sec (do I need exposure?)
+    const Vector<double>& freq = acc.frequency();
+    Axis::ShPtr timeAxis(new RegularAxis ());
+    Axis::ShPtr freqAxis(new RegularAxis ());
+  }
+      
+  void LofarGridder::finaliseGrid(casa::Array<double>& out)
+  {
+    // If averaging is done, it has to grid the last timeslots.
+    itsGridder->finaliseGrid (out);
+  }
+
+  void LofarGridder::finaliseWeights(casa::Array<double>& out)
+  {
+    itsGridder->finaliseWeights (out);
+  }
+
+  void LofarGridder::initialiseDegrid(const scimath::Axes& axes,
+                                      const casa::Array<double>& image)
+  {
+    itsGridder->initialiseDegrid (axes, image);
+  }
+
+  void LofarGridder::customiseForContext(casa::String context)
+  {
+    itsGridder->customiseForContext (context);
+  }
+      
+  void LofarGridder::initVisWeights(IVisWeights::ShPtr viswt)
   {
+    itsGridder->initVisWeights (viswt);
   }
       
-  void LofarGridder::initConvolutionFunction(const IConstDataAccessor&)
+  void LofarGridder::degrid(IDataAccessor& acc)
   {
-    itsSupport=0;
-    itsOverSample=1;
-    itsCSize=2*(itsSupport+1)*itsOverSample+1; // 3
-    itsCCenter=(itsCSize-1)/2; // 1
-    itsConvFunc.resize(1);
-    itsConvFunc[0].resize(itsCSize, itsCSize); // 3, 3, 1
-    itsConvFunc[0].set(0.0);
-    itsConvFunc[0](itsCCenter,itsCCenter)=1.0; // 1,1,0 = 1
+    if (! itsInitialized) {
+      initCorrections (acc);
     }
+    // do correction, uvw rotation, and averaging here
+    // For time being only correction is done.
+    OnDemandBufferDataAccessor acc2(acc);
+    Array<Complex>& data = acc2.rwVisibility();
+    itsGridder->degrid (acc2);
+  }
+
+  void LofarGridder::finaliseDegrid()
+  {
+    itsGridder->finaliseDegrid();
+  }
+
+//   void LofarGridder::initIndices(const IConstDataAccessor& acc) 
+//   {
+//     itsGridder->initIndices (acc);
+//   }
+
+//   void LofarGridder::initConvolutionFunction(const IConstDataAccessor& acc)
+//   {
+//     itsGridder->initConvolutionFunction (acc);
+//   }
+    
+//   void LofarGridder::correctConvolution(casa::Array<double>& image)
+//   {
+//     itsGridder->correctConvolution (image);
+//   }
 
-  void LofarGridder::correctConvolution(casa::Array<double>& image)
+  void LofarGridder::initCorrections (const IConstDataAccessor& acc)
   {
+    const TableConstDataAccessor& tacc =
+      dynamic_cast<const TableConstDataAccessor&>(acc);
+    const TableConstDataIterator& titer = tacc.iterator();
+    Table ms = titer.table();
+    Table antTab(ms.keywordSet().asTable("ANTENNA"));
+    itsInitialized = true;
   }
 
 } //# end namespace

CEP/Imager/LofarGridder/src/LofarWProjectGridder.cc

+//# LofarWProjectGridder.cc: WProjection gridder for LOFAR data correcting for DD effects
+//#
+//# Copyright (C) 2010
+//# ASTRON (Netherlands Institute for Radio Astronomy)
+//# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//#
+//# This file is part of the LOFAR software suite.
+//# The LOFAR software suite is free software: you can redistribute it and/or
+//# modify it under the terms of the GNU General Public License as published
+//# by the Free Software Foundation, either version 3 of the License, or
+//# (at your option) any later version.
+//#
+//# The LOFAR software suite is distributed in the hope that it will be useful,
+//# but WITHOUT ANY WARRANTY; without even the implied warranty of
+//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//# GNU General Public License for more details.
+//#
+//# You should have received a copy of the GNU General Public License along
+//# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//#
+//# $Id$
+//#
+//# @author Ger van Diepen <diepen at astron dot nl>
+
+#include <lofar_config.h>
+#include <ParmDB/Grid.h>
+
+//# ASKAP includes
+#include <LofarGridder/LofarWProjectGridder.h>
+#include <gridding/VisGridderFactory.h>
+#include <dataaccess/TableConstDataAccessor.h>
+#include <dataaccess/TableConstDataIterator.h>
+
+//# casacore includes
+#include <tables/Tables/TableRecord.h>
+#include <ms/MeasurementSets/MSColumns.h>
+#include <ms/MeasurementSets/MSSpWindowColumns.h>
+#include <ms/MeasurementSets/MSAntennaColumns.h>
+#include <ms/MeasurementSets/MSObsColumns.h>
+#include <measures/Measures/MeasTable.h>
+#include <casa/OS/Path.h>
+
+using namespace LOFAR::BBS;
+using namespace casa;
+using namespace askap;
+using namespace askap::synthesis;
+
+namespace LOFAR
+{
+
+  LofarWProjectGridder::LofarWProjectGridder (const ParameterSet& parset)
+    : WProjectVisGridder (*(dynamic_cast<const WProjectVisGridder*>
+                            (WProjectVisGridder::createGridder(parset).get()))),
+      itsInitialized (false)
+  {
+    // Currently only beam can be corrected for.
+    ///    itsTimeAvg = parset.getInt ("average.timestep", 1);
+    ///    itsFreqAvg = parset.getInt ("average.freqstep", 1);
+    ///    itsCorrect = parset.getStringVector ("correct");
+    itsParmDBName = parset.getString ("ParmDB.Instrument", "");
+    itsConfigName = parset.getString ("Beam.StationConfig.Name");
+    itsConfigPath = parset.getString ("Beam.StationConfig.Path");
+  }
+
+  LofarWProjectGridder::~LofarWProjectGridder()
+  {}
+
+  // Clone a copy of this Gridder
+  IVisGridder::ShPtr LofarWProjectGridder::clone()
+  {
+    return IVisGridder::ShPtr(new LofarWProjectGridder(*this));
+  }
+
+  IVisGridder::ShPtr LofarWProjectGridder::makeGridder
+  (const ParameterSet& parset)
+  {
+    std::cout << "LofarWProjectGridder::makeGridder" << std::endl;
+    return IVisGridder::ShPtr(new LofarWProjectGridder(parset));
+  }
+
+  const std::string& LofarWProjectGridder::gridderName()
+  {
+    static std::string name("LofarGridder.WProject");
+    return name;
+  }
+
+  void LofarWProjectGridder::registerGridder()
+  {
+    VisGridderFactory::registerGridder (gridderName(), &makeGridder);
+  }
+
+  void LofarWProjectGridder::correctVisibilities (IDataAccessor& acc,
+                                                  bool forward)
+  {
+    if (!itsInitialized) {
+      initCorrections (acc);
+    }
+    // If changed, set the facet center as direction of interest.
+    MVDirection facetCenter = getImageCentre();
+    if (facetCenter != itsLastFacetCenter) {
+      itsLastFacetCenter = facetCenter;
+      itsResponse->setDirection (MDirection(facetCenter, MDirection::J2000));
+    }
+    // Get read/write access to the visibilities. It makes a copy if needed.
+    Array<Complex>& visData = acc.rwVisibility();
+    ASSERT (visData.contiguousStorage());
+    Complex* data = visData.data();
+    uInt nfreq = visData.shape()[1];
+    double time = acc.time();
+    // Set the axes for the BBS request.
+    Axis::ShPtr freqAxis
+      (new RegularAxis(itsStartFreq-itsFreqInterval/2, itsFreqInterval, nfreq));
+    Axis::ShPtr timeAxis
+      (new RegularAxis(time-itsTimeInterval/2, itsTimeInterval, 1));
+    Grid grid(timeAxis, freqAxis);
+    itsResponse->setEvalGrid (grid);
+    const casa::Vector<uInt>& ant1 = acc.antenna1();
+    const casa::Vector<uInt>& ant2 = acc.antenna2();
+    for (uInt i=0; i<ant1.size(); ++i) {
+      JonesMatrix::View j1 = itsResponse->evaluate (ant1[i]);
+      JonesMatrix::View j2 = itsResponse->evaluate (ant2[i]);
+      for (uInt j=0; j<nfreq; ++j) {
+        dcomplex cl0 = j1(j,0).getDComplex(0,0);
+        dcomplex cl1 = j1(j,0).getDComplex(1,0);
+        dcomplex cl2 = j1(j,0).getDComplex(0,1);
+        dcomplex cl3 = j1(j,0).getDComplex(1,1);
+        dcomplex cr0 = conj (j2(j,0).getDComplex(0,0));
+        dcomplex cr1 = conj (j2(j,0).getDComplex(1,0));
+        dcomplex cr2 = conj (j2(j,0).getDComplex(0,1));
+        dcomplex cr3 = conj (j2(j,0).getDComplex(1,1));
+        dcomplex d0 = data[0];
+        dcomplex d1 = data[1];
+        dcomplex d2 = data[2];
+        dcomplex d3 = data[3];
+        if (forward) {
+          // For gridding use inverse of responses.
+          //   data = inv(res(i)) * data * inv(conjtranspose(res(j)))
+          // Note that inv(conjtranspose(A)) == conjtranspose(inv(A))
+          // inv(a b) == ( d -b)  / (ad-bc)
+          //    (c d)    (-c  a)
+          dcomplex tmp0 = cl3*d0 - cl1*d2;
+          dcomplex tmp1 = cl3*d1 - cl1*d3;
+          dcomplex tmp2 = cl0*d2 - cl2*d0;
+          dcomplex tmp3 = cl0*d3 - cl2*d1;
+          dcomplex factor = 1. / ((cl0*cl3 - cl1*cl2) * (cr0*cr3 - cr1*cr2));
+          data[0] = factor * (tmp0*cr3 - tmp1*cr1);
+          data[1] = factor * (tmp1*cr0 - tmp0*cr2);
+          data[2] = factor * (tmp2*cr3 - tmp3*cr1);
+          data[3] = factor * (tmp3*cr0 - tmp2*cr2);
+        } else {
+          // Degridding; data = res(i) * data * conjtranspose(res(j))
+          dcomplex tmp0 = cl0*d0 + cl1*d2;
+          dcomplex tmp1 = cl0*d1 + cl1*d3;
+          dcomplex tmp2 = cl2*d0 + cl3*d2;
+          dcomplex tmp3 = cl2*d1 + cl3*d3;
+          data[0] = tmp0*cr0 + tmp1*cr1;
+          data[1] = tmp0*cr2 + tmp1*cr3;
+          data[2] = tmp2*cr0 + tmp3*cr1;
+          data[3] = tmp2*cr2 + tmp3*cr3;
+        }
+        data += 4;
+      }
+    }
+  }
+
+  void LofarWProjectGridder::initCorrections (const IConstDataAccessor& acc)
+  {
+    // Ensure there are polarizations XX,XY,YX,YY.
+    ASSERT (acc.stokes().size() == 4);
+    ASSERT (acc.stokes()[0] == Stokes::XX);
+    ASSERT (acc.stokes()[1] == Stokes::XY);
+    ASSERT (acc.stokes()[2] == Stokes::YX);
+    ASSERT (acc.stokes()[3] == Stokes::YY);
+    const TableConstDataAccessor& tacc =
+      dynamic_cast<const TableConstDataAccessor&>(acc);
+    const TableConstDataIterator& titer = tacc.iterator();
+    MeasurementSet ms(titer.table());
+    // Get the time interval.
+    ROMSColumns msCol(ms);
+    // Get the reference frequency and width.
+    // It assumes that the interval is constant for this band.
+    ROMSSpWindowColumns spwCols(ms.spectralWindow());
+    double refFreq = spwCols.refFrequency()(0);
+    itsFreqInterval = spwCols.chanWidth()(0).data()[0];
+    casa::Vector<double> freqs = spwCols.chanFreq()(0);
+    itsStartFreq = freqs[0];
+    for (uInt i=1; i<freqs.size(); ++i) {
+      ASSERTSTR (casa::near(freqs[i], itsStartFreq + i*itsFreqInterval),
+                 "Frequency channels are not regularly spaced");
+    }
+    // Form the StationResponse object.
+    itsResponse = StationResponse::Ptr (new StationResponse
+                                        (makeInstrument(ms),
+                                         itsConfigName, Path(itsConfigPath),
+                                         refFreq));
+    // Set the pointing direction (for beamforming).
+    // Use first value of MDirection array in first row in FIELD subtable.
+    ROMSFieldColumns fieldCols(ms.field());
+    itsResponse->setPointing (fieldCols.delayDirMeasCol()(0).data()[0]);
+    itsTimeInterval = msCol.interval()(0);
+    // Initialization is done.
+    itsInitialized = true;
+  }
+
+  BBS::Instrument LofarWProjectGridder::makeInstrument
+  (const MeasurementSet& ms)
+  {
+    // Get all station positions.
+    ROMSAntennaColumns antCols(ms.antenna());
+    vector<Station> stations;
+    stations.reserve (ms.antenna().nrow());
+    for (uInt i=0; i<ms.antenna().nrow(); ++i) {
+      stations.push_back (Station(antCols.name()(i),
+                                  antCols.positionMeas()(i)));
+    }
+    // Find observatory position.
+    // If not found, set it to the position of the middle station.
+    MPosition arrayPos;
+    Bool fndObs = False;
+    if (ms.observation().nrow() > 0) {
+      ROMSObservationColumns obsCols(ms.observation());
+      String telescope = obsCols.telescopeName()(0);
+      fndObs = MeasTable::Observatory (arrayPos, telescope);
+    }
+    if (!fndObs  &&  stations.size() > 0) {
+      arrayPos = stations[stations.size()/2].position();
+    }
+    return BBS::Instrument ("LOFAR", arrayPos, stations);
+  }
+
+} //# end namespace

CEP/Imager/LofarGridder/src/Register.cc

@@ -23,10 +23,11 @@
 //# @author Ger van Diepen <diepen at astron dot nl>
 
 //# Includes
+#include <lofar_config.h>
 #include <LofarGridder/Register.h>
-#include <LofarGridder/LofarGridder.h>
+#include <LofarGridder/LofarWProjectGridder.h>
 
 void register_lofargridder()
 {
-  LOFAR::LofarGridder::registerGridder();
+  LOFAR::LofarWProjectGridder::registerGridder();
 }

CEP/Imager/LofarGridder/test/tLofarGridder.cc

@@ -37,7 +37,7 @@ int main (int argc, const char** argv)
 
     // Create a LofarGridder instance.
     ParameterSet parset;
-    parset.add ("gridder", "LofarGridder");
+    parset.add ("gridder", "LofarGridder.WProject");
     IVisGridder::ShPtr gridder = VisGridderFactory::make (parset);
 
   } catch (const std::exception& x) {