Changeset 2671

Timestamp:

10/18/12 18:36:25 (12 years ago)

Author:

Takeshi Nakazato

Message:

New Development: No

JIRA Issue: Yes CAS-4429

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Implemented gauss and gjinc gridding. Added some arguments to
pass necessary parameters for those grid functions.

Also, defined helper function that plots current grid function
against radius in pixel.

Location:

trunk

Files:

• python/asapgrid.py (modified) (5 diffs)
• src/STGrid.cpp (modified) (9 diffs)
• src/STGrid.h (modified) (4 diffs)
• src/python_STGrid.cpp (modified) (4 diffs)

trunk/python/asapgrid.py

@@ -135 +135 @@
         self.gridder._defineimage( nx, ny, cellx, celly, center )
 
-    def setFunc( self, func='box', width=-1 ):
+    def setFunc( self, func='box', width=-1, gwidth="", jwidth="" ):
         """
         Set convolution function. Possible options are 'box' (Box-car,
@@ -151 +151 @@
                  choose pre-defined value for each convolution function.
         """
-        self.gridder._setfunc( func, width )
+        fname = func.upper()
+        if fname == 'GAUSS' or fname == 'GJINC':
+            gw = str(gwidth)
+            jw = str(jwidth)
+            w = str(width)
+            if w[0] == '-': w = ''
+            #self.gridder._setfunc(fname, -1, w, gw, jw)
+            self.gridder._setfunc(fname,convsupport=-1,gwidth=gw,jwidth=jw,truncate=w)
+        else:
+            self.gridder._setfunc( func, convsupport=width )
 
     def setWeight( self, weightType='uniform' ):
@@ -223 +232 @@
         asaplog.push('plot: elapsed time %s sec'%(t1-t0))
         asaplog.post('DEBUG','asapgrid.plot')
+
+    def plotFunc(self, clear=True):
+        """
+        Support function to see the shape of current grid function.
+
+        clear -- clear panel if True. Default is True.
+        """
+        pl.figure(11)
+        if clear:
+            pl.clf()
+        f = self.gridder._getfunc()
+        convsampling = 100
+        a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
+        pl.plot(a,f,'.-')
+        pl.xlabel('pixel')
+        pl.ylabel('convFunc')
         
 class asapgrid2:
@@ -367 +392 @@
                  choose pre-defined value for each convolution function.
         """
-        self.gridder._setfunc( func, width )
+        fname = func.upper()
+        if fname == 'GAUSS' or fname == 'GJINC':
+            gw = str(gwidth)
+            jw = str(jwidth)
+            w = str(width)
+            if w[0] == '-': w = ''
+            #self.gridder._setfunc(fname, -1, w, gw, jw)
+            self.gridder._setfunc(fname,convsupport=-1,gwidth=gw,jwidth=jw,truncate=w)
+        else:
+            self.gridder._setfunc( func, convsupport=width )
 
     def setWeight( self, weightType='uniform' ):
@@ -408 +442 @@
         tp = 0 if rcParams['scantable.storage']=='memory' else 1
         return scantable( self.gridder._get( tp ), average=False )
+    
+    def plotFunc(self, clear=True):
+        """
+        Support function to see the shape of current grid function.
+
+        clear -- clear panel if True. Default is True.
+        """
+        pl.figure(11)
+        if clear:
+            pl.clf()
+        f = self.gridder._getfunc()
+        convsampling = 100
+        a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
+        pl.plot(a,f,'.-')
+        pl.xlabel('pixel')
+        pl.ylabel('convFunc')
 
 class _SDGridPlotter:

trunk/src/STGrid.cpp

@@ -83 +83 @@
   convSupport_ = -1 ;
   userSupport_ = -1 ;
+  truncate_ = "";
+  gwidth_ = "";
+  jwidth_ = "";
   convSampling_ = 100 ;
   nprocessed_ = 0 ;
@@ -143 +146 @@
   
 void STGrid::setFunc( string convType,
-                      int convSupport ) 
+                      int convSupport,
+                      string truncate,
+                      string gwidth,
+                      string jwidth ) 
 {
   convType_ = String( convType ) ;
   convType_.upcase() ;
   userSupport_ = (Int)convSupport ;
+  truncate_ = String( truncate );
+  gwidth_ = String( gwidth );
+  jwidth_ = String( jwidth );
 }
 
@@ -429 +438 @@
   }
 
+  // Warn if gauss or gjinc gridding with non-square cells
+  if ((cellx_ != celly_) && (convType_=="GAUSS"||convType_=="GJINC")) {
+    os << LogIO::WARN 
+       << "The " << convType_ << " gridding doesn't support non-square grid." << endl
+       << "Result may be wrong." << LogIO::POST;
+  }
+
   // grid parameter
   os << LogIO::DEBUGGING ;
@@ -442 +458 @@
   os << "   center = " << center_ << endl ;
   os << "   weighting = " << wtype_ << endl ;
-  os << "   convfunc = " << convType_ << " with support " << convSupport_ << endl ; 
+  os << "   convfunc = " << convType_ << endl;
+  if (convType_ == "GAUSS") {
+    os << "      gwidth = " << gwidth_ << endl;
+    os << "      truncate = " << truncate_ << endl;
+  }
+  else if (convType_ == "GJINC") {
+    os << "      gwidth = " << gwidth_ << endl;
+    os << "      jwidth = " << jwidth_ << endl;
+    os << "      truncate = " << truncate_ << endl;
+  }
+  else {
+    os << "      support = " << convSupport_ << endl;
+  }
   os << "   doclip = " << (doclip_?"True":"False") << endl ;
   os << "----------" << LogIO::POST ;
@@ -1567 +1595 @@
 #if defined(NEED_UNDERSCORES)
 #define grdsf grdsf_
+#define grdgauss grdgauss_
+#define grdjinc1 grdjinc1_
 #endif
 extern "C" { 
-   void grdsf(Double*, Double*);
+  void grdsf(Double*, Double*);
+  void grdgauss(Double*, Double*, Double*); 
+  void grdjinc1(Double*, Double*, Int*, Double*);
 }
 void STGrid::spheroidalFunc( Vector<Float> &convFunc ) 
@@ -1582 +1614 @@
 }
 
-void STGrid::gaussFunc( Vector<Float> &convFunc ) 
+void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate ) 
 {
   convFunc = 0.0 ;
-  // HWHM of the Gaussian is convSupport_ / 4
-  // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
-  Int len = convSampling_ * convSupport_ ;
-  Double hwhm = len * 0.25 ;
+  Int len = (Int)(truncate*Double(convSampling_)+0.5);
+  Double out, val;
   for ( Int i = 0 ; i < len ; i++ ) {
-    Double val = Double(i) / hwhm ;
-    convFunc(i) = exp( -log(2)*val*val ) ;
+    val = Double(i) / Double(convSampling_) ;
+    grdgauss(&hwhm, &val, &out);
+    convFunc(i) = out;
+  }
+}
+
+void STGrid::gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate )
+{
+  convFunc = 0.0;
+  Double out1, out2, val;
+  Int normalize = 1;
+  if  (truncate >= 0.0) { 
+    Int len = (Int)(truncate*Double(convSampling_)+0.5);
+    for (Int i = 0 ; i < len ; i++) {
+      val = Double(i) / Double(convSampling_);
+      grdgauss(&hwhm, &val, &out1);
+      grdjinc1(&c, &val, &normalize, &out2);
+      convFunc(i) = out1 * out2;
+    }
+  }
+  else {
+    Int len = convFunc.nelements();
+    for (Int i = 0 ; i < len ; i++) {
+      val = Double(i) / Double(convSampling_);
+      grdjinc1(&c, &val, &normalize, &out2);
+      if (out2 <= 0.0) {
+        LogIO os(LogOrigin("STGrid","gjincFunc",WHERE));
+        os << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << val << LogIO::POST;
+        break;
+      }
+      grdgauss(&hwhm, &val, &out1);
+      convFunc(i) = out1 * out2;
+    }
   }
 }
@@ -1600 +1661 @@
 }
 
+vector<float> STGrid::getConvFunc()
+{
+  LogIO os(LogOrigin("STGrid","getConvFunc",WHERE));
+  Vector<Float> convFunc;
+  vector<float> out;
+
+  if (cellx_ <= 0.0 || celly_ <= 0.0) {
+    selectData();
+    setupGrid();
+  }
+
+  if (convType_ == "BOX" || convType_ == "SF") {
+    setConvFunc(convFunc);
+  }
+  else if (convType_ == "GAUSS") {
+    Quantum<Double> q1,q2;
+    readQuantity(q1,gwidth_);
+    readQuantity(q2,truncate_);
+//     if (celly_ <= 0.0 
+//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
+//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel"))) {
+//       throw AipsError("You have to call defineImage to get correct convFunc");
+//     }
+    setConvFunc(convFunc);
+  }
+  else if (convType_ == "GJINC") {
+    Quantum<Double> q1,q2,q3;
+    readQuantity(q1,gwidth_);
+    readQuantity(q2,truncate_);
+    readQuantity(q3,jwidth_);
+//     if (celly_ <= 0.0 
+//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
+//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel") ||
+//             (!q3.getUnit().empty()&&q3.getUnit()!="pixel"))) {
+//       throw AipsError("You have to call defineImage to get correct convFunc");
+//     }
+    setConvFunc(convFunc);
+  }
+  else if (convType_ == "PB") {
+    throw AipsError("Grid function PB is not available");
+  }
+  else {
+    throw AipsError("Unknown grid function: "+convType_);
+  }
+
+  convFunc.tovector(out);
+  return out;
+}
+
 void STGrid::setConvFunc( Vector<Float> &convFunc )
 {
+  LogIO os(LogOrigin("STGrid","setConvFunc",WHERE));
   convSupport_ = userSupport_ ;
   if ( convType_ == "BOX" ) {
@@ -1609 +1720 @@
     convFunc.resize( convSize ) ;
     boxFunc( convFunc, convSize ) ;
+    os << LogIO::DEBUGGING
+       << "convType_ = " << convType_ << endl
+       << "convSupport_ = " << convSupport_ << LogIO::POST;
   }
   else if ( convType_ == "SF" ) {
@@ -1616 +1730 @@
     convFunc.resize( convSize ) ;
     spheroidalFunc( convFunc ) ;
+    os << LogIO::DEBUGGING
+       << "convType_ = " << convType_ << endl
+       << "convSupport_ = " << convSupport_ << LogIO::POST;
   }
   else if ( convType_ == "GAUSS" ) {
-    // to take into account Gaussian tail
-    if ( convSupport_ < 0 )
-      convSupport_ = 4 ; // 1 * 4
-    else {
-      convSupport_ = userSupport_ * 4 ;
-    }
-    Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
+    // determine pixel gwidth
+    // default is HWHM corresponding to b = 1.0 (Mangum et al. 2007)
+    Double pixelGW = sqrt(log(2.0)); 
+    Quantum<Double> q ;
+    if (!gwidth_.empty()) {
+      readQuantity( q, gwidth_ );
+      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
+        pixelGW = q.getValue();
+      }
+      else {
+        pixelGW = q.getValue("rad")/celly_;
+      }
+    }
+    // determine truncation radius
+    // default is 3 * HWHM
+    Double truncate = 3.0 * pixelGW; 
+    if (!truncate_.empty()) {
+      readQuantity( q, truncate_ );
+      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
+        truncate = q.getValue();
+      }
+      else {
+        truncate = q.getValue("rad")/celly_;
+      }
+    }      
+    convSupport_ = (Int)(truncate+0.5);
+    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
     convFunc.resize( convSize ) ;
-    gaussFunc( convFunc ) ;
+    gaussFunc( convFunc, pixelGW, truncate ) ;
+    os << LogIO::DEBUGGING
+       << "convType_ = " << convType_ << endl
+       << "convSupport_ = " << convSupport_ << endl
+       << "truncate_ = " << truncate << endl
+       << "gwidth_ = " << pixelGW << "pixel" << LogIO::POST;
+  }
+  else if ( convType_ == "GJINC" ) {
+    // determine pixel gwidth
+    // default is HWHM corresponding to b = 2.52 (Mangum et al. 2007)
+    Double pixelGW = sqrt(log(2.0)) * 2.52; 
+    Quantum<Double> q ;
+    if (!gwidth_.empty()) {
+      readQuantity( q, gwidth_ );
+      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
+        pixelGW = q.getValue();
+      }
+      else {
+        pixelGW = q.getValue("rad")/celly_;
+      }
+    }
+    // determine pixel c
+    // default is c = 1.55 (Mangum et al. 2007)
+    Double pixelJW = 1.55;
+    if (!jwidth_.empty()) {
+      readQuantity( q, jwidth_ );
+      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
+        pixelJW = q.getValue();
+      }
+      else {
+        pixelJW = q.getValue("rad")/celly_;
+      }
+    }      
+    // determine truncation radius
+    // default is -1.0 (truncate at first null)
+    Double truncate = -1.0;
+    if (!truncate_.empty()) {
+      readQuantity( q, truncate_ );
+      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
+        truncate = q.getValue();
+      }
+      else {
+        truncate = q.getValue("rad")/celly_;
+      }
+    }      
+    convSupport_ = (truncate >= 0.0) ? (Int)(truncate+0.5) : (Int)(2*pixelJW+0.5);
+    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
+    convFunc.resize( convSize ) ;
+    gjincFunc( convFunc, pixelGW, pixelJW, truncate ) ;
+    os << LogIO::DEBUGGING
+       << "convType_ = " << convType_ << endl
+       << "convSupport_ = " << convSupport_ << endl
+       << "truncate_ = " << truncate << endl
+       << "gwidth_ = " << pixelGW << "pixel" << endl
+       << "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
   }
   else if ( convType_ == "PB" ) {

trunk/src/STGrid.h

@@ -61 +61 @@
                     string scenter="" ) ;
   void setFunc( string convtype="box",
-                int convsupport=-1 ) ;
+                int convsupport=-1,
+                string truncate="",
+                string gwidth="",
+                string jwidth="" ) ;
 
   void setWeight( const string wType="uniform" ) ;
@@ -71 +74 @@
   
   string saveData( string outfile="" ) ;
+
+  // support function to know how grid function looks like
+  vector<float> getConvFunc();
 
 //private:
@@ -138 +144 @@
   void boxFunc( Vector<Float> &convFunc, Int &convSize ) ;
   void spheroidalFunc( Vector<Float> &convFunc ) ;
-  void gaussFunc( Vector<Float> &convFunc ) ;
+  void gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate ) ;
+  void gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate );
   void pbFunc( Vector<Float> &convFunc ) ;
   void setConvFunc( Vector<Float> &convFunc ) ;
@@ -232 +239 @@
   Int convSupport_ ;
   Int userSupport_ ;
+  String gwidth_;
+  String jwidth_;
+  String truncate_;
   Int convSampling_ ;
   Vector<uInt> pollist_ ;

trunk/src/python_STGrid.cpp

@@ -24 +24 @@
     .def("_setpollist", &STGrid::setPolList)
     .def("_setscanlist", &STGrid::setScanList)
-    .def("_defineimage", &STGrid::defineImage)
-    .def("_setfunc", &STGrid::setFunc)
+    .def("_defineimage", &STGrid::defineImage,
+         (boost::python::arg("nx")=-1,
+          boost::python::arg("ny")=-1,
+          boost::python::arg("scellx")="",
+          boost::python::arg("scelly")="",
+          boost::python::arg("scenter")=""))
+    .def("_setfunc", &STGrid::setFunc,
+         (boost::python::arg("convsupport")=-1,
+          boost::python::arg("truncate")="",
+          boost::python::arg("gwidth")="",
+          boost::python::arg("jwidth")=""))
     .def("_grid", &STGrid::grid)
     .def("_setin", &STGrid::setFileIn)
@@ -33 +42 @@
     .def("_disableclip", &STGrid::disableClip) 
     .def("_save", &STGrid::saveData)
+    .def("_getfunc", &STGrid::getConvFunc)
     ;
 
@@ -42 +52 @@
     .def("_setpollist", &STGrid2::setPolList)
     .def("_setscanlist", &STGrid2::setScanList)
-    .def("_defineimage", &STGrid2::defineImage)
-    .def("_setfunc", &STGrid2::setFunc)
+    .def("_defineimage", &STGrid2::defineImage,
+         (boost::python::arg("nx")=-1,
+          boost::python::arg("ny")=-1,
+          boost::python::arg("scellx")="",
+          boost::python::arg("scelly")="",
+          boost::python::arg("scenter")=""))
+    .def("_setfunc", &STGrid2::setFunc,
+         (boost::python::arg("convsupport")=-1,
+          boost::python::arg("truncate")="",
+          boost::python::arg("gwidth")="",
+          boost::python::arg("jwidth")=""))
     .def("_grid", &STGrid2::grid)
     .def("_setin", &STGrid2::setScantable)
@@ -51 +70 @@
     .def("_disableclip", &STGrid2::disableClip) 
     .def("_get", &STGrid2::getResultAsScantable)
+    .def("_getfunc", &STGrid2::getConvFunc)
     ;