Changeset 2671 for trunk/src/STGrid.cpp

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.

File:

• trunk/src/STGrid.cpp (modified) (9 diffs)

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" ) {