Changeset 2671


Ignore:
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:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/asapgrid.py

    r2669 r2671  
    135135        self.gridder._defineimage( nx, ny, cellx, celly, center )
    136136
    137     def setFunc( self, func='box', width=-1 ):
     137    def setFunc( self, func='box', width=-1, gwidth="", jwidth="" ):
    138138        """
    139139        Set convolution function. Possible options are 'box' (Box-car,
    … …  
    151151                 choose pre-defined value for each convolution function.
    152152        """
    153         self.gridder._setfunc( func, width )
     153        fname = func.upper()
     154        if fname == 'GAUSS' or fname == 'GJINC':
     155            gw = str(gwidth)
     156            jw = str(jwidth)
     157            w = str(width)
     158            if w[0] == '-': w = ''
     159            #self.gridder._setfunc(fname, -1, w, gw, jw)
     160            self.gridder._setfunc(fname,convsupport=-1,gwidth=gw,jwidth=jw,truncate=w)
     161        else:
     162            self.gridder._setfunc( func, convsupport=width )
    154163
    155164    def setWeight( self, weightType='uniform' ):
    … …  
    223232        asaplog.push('plot: elapsed time %s sec'%(t1-t0))
    224233        asaplog.post('DEBUG','asapgrid.plot')
     234
     235    def plotFunc(self, clear=True):
     236        """
     237        Support function to see the shape of current grid function.
     238
     239        clear -- clear panel if True. Default is True.
     240        """
     241        pl.figure(11)
     242        if clear:
     243            pl.clf()
     244        f = self.gridder._getfunc()
     245        convsampling = 100
     246        a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
     247        pl.plot(a,f,'.-')
     248        pl.xlabel('pixel')
     249        pl.ylabel('convFunc')
    225250       
    226251class asapgrid2:
    … …  
    367392                 choose pre-defined value for each convolution function.
    368393        """
    369         self.gridder._setfunc( func, width )
     394        fname = func.upper()
     395        if fname == 'GAUSS' or fname == 'GJINC':
     396            gw = str(gwidth)
     397            jw = str(jwidth)
     398            w = str(width)
     399            if w[0] == '-': w = ''
     400            #self.gridder._setfunc(fname, -1, w, gw, jw)
     401            self.gridder._setfunc(fname,convsupport=-1,gwidth=gw,jwidth=jw,truncate=w)
     402        else:
     403            self.gridder._setfunc( func, convsupport=width )
    370404
    371405    def setWeight( self, weightType='uniform' ):
    … …  
    408442        tp = 0 if rcParams['scantable.storage']=='memory' else 1
    409443        return scantable( self.gridder._get( tp ), average=False )
     444   
     445    def plotFunc(self, clear=True):
     446        """
     447        Support function to see the shape of current grid function.
     448
     449        clear -- clear panel if True. Default is True.
     450        """
     451        pl.figure(11)
     452        if clear:
     453            pl.clf()
     454        f = self.gridder._getfunc()
     455        convsampling = 100
     456        a = numpy.arange(0,len(f)/convsampling,1./convsampling,dtype=float)
     457        pl.plot(a,f,'.-')
     458        pl.xlabel('pixel')
     459        pl.ylabel('convFunc')
    410460
    411461class _SDGridPlotter:

  • trunk/src/STGrid.cpp

    r2669 r2671  
    8383  convSupport_ = -1 ;
    8484  userSupport_ = -1 ;
     85  truncate_ = "";
     86  gwidth_ = "";
     87  jwidth_ = "";
    8588  convSampling_ = 100 ;
    8689  nprocessed_ = 0 ;
    … …  
    143146 
    144147void STGrid::setFunc( string convType,
    145                       int convSupport )
     148                      int convSupport,
     149                      string truncate,
     150                      string gwidth,
     151                      string jwidth )
    146152{
    147153  convType_ = String( convType ) ;
    148154  convType_.upcase() ;
    149155  userSupport_ = (Int)convSupport ;
     156  truncate_ = String( truncate );
     157  gwidth_ = String( gwidth );
     158  jwidth_ = String( jwidth );
    150159}
    151160
    … …  
    429438  }
    430439
     440  // Warn if gauss or gjinc gridding with non-square cells
     441  if ((cellx_ != celly_) && (convType_=="GAUSS"||convType_=="GJINC")) {
     442    os << LogIO::WARN
     443       << "The " << convType_ << " gridding doesn't support non-square grid." << endl
     444       << "Result may be wrong." << LogIO::POST;
     445  }
     446
    431447  // grid parameter
    432448  os << LogIO::DEBUGGING ;
    … …  
    442458  os << "   center = " << center_ << endl ;
    443459  os << "   weighting = " << wtype_ << endl ;
    444   os << "   convfunc = " << convType_ << " with support " << convSupport_ << endl ;
     460  os << "   convfunc = " << convType_ << endl;
     461  if (convType_ == "GAUSS") {
     462    os << "      gwidth = " << gwidth_ << endl;
     463    os << "      truncate = " << truncate_ << endl;
     464  }
     465  else if (convType_ == "GJINC") {
     466    os << "      gwidth = " << gwidth_ << endl;
     467    os << "      jwidth = " << jwidth_ << endl;
     468    os << "      truncate = " << truncate_ << endl;
     469  }
     470  else {
     471    os << "      support = " << convSupport_ << endl;
     472  }
    445473  os << "   doclip = " << (doclip_?"True":"False") << endl ;
    446474  os << "----------" << LogIO::POST ;
    … …  
    15671595#if defined(NEED_UNDERSCORES)
    15681596#define grdsf grdsf_
     1597#define grdgauss grdgauss_
     1598#define grdjinc1 grdjinc1_
    15691599#endif
    15701600extern "C" {
    1571    void grdsf(Double*, Double*);
     1601  void grdsf(Double*, Double*);
     1602  void grdgauss(Double*, Double*, Double*);
     1603  void grdjinc1(Double*, Double*, Int*, Double*);
    15721604}
    15731605void STGrid::spheroidalFunc( Vector<Float> &convFunc )
    … …  
    15821614}
    15831615
    1584 void STGrid::gaussFunc( Vector<Float> &convFunc )
     1616void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate )
    15851617{
    15861618  convFunc = 0.0 ;
    1587   // HWHM of the Gaussian is convSupport_ / 4
    1588   // To take into account Gaussian tail, kernel cutoff is set to 4 * HWHM
    1589   Int len = convSampling_ * convSupport_ ;
    1590   Double hwhm = len * 0.25 ;
     1619  Int len = (Int)(truncate*Double(convSampling_)+0.5);
     1620  Double out, val;
    15911621  for ( Int i = 0 ; i < len ; i++ ) {
    1592     Double val = Double(i) / hwhm ;
    1593     convFunc(i) = exp( -log(2)*val*val ) ;
     1622    val = Double(i) / Double(convSampling_) ;
     1623    grdgauss(&hwhm, &val, &out);
     1624    convFunc(i) = out;
     1625  }
     1626}
     1627
     1628void STGrid::gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate )
     1629{
     1630  convFunc = 0.0;
     1631  Double out1, out2, val;
     1632  Int normalize = 1;
     1633  if  (truncate >= 0.0) {
     1634    Int len = (Int)(truncate*Double(convSampling_)+0.5);
     1635    for (Int i = 0 ; i < len ; i++) {
     1636      val = Double(i) / Double(convSampling_);
     1637      grdgauss(&hwhm, &val, &out1);
     1638      grdjinc1(&c, &val, &normalize, &out2);
     1639      convFunc(i) = out1 * out2;
     1640    }
     1641  }
     1642  else {
     1643    Int len = convFunc.nelements();
     1644    for (Int i = 0 ; i < len ; i++) {
     1645      val = Double(i) / Double(convSampling_);
     1646      grdjinc1(&c, &val, &normalize, &out2);
     1647      if (out2 <= 0.0) {
     1648        LogIO os(LogOrigin("STGrid","gjincFunc",WHERE));
     1649        os << LogIO::DEBUG1 << "convFunc is automatically truncated at radius " << val << LogIO::POST;
     1650        break;
     1651      }
     1652      grdgauss(&hwhm, &val, &out1);
     1653      convFunc(i) = out1 * out2;
     1654    }
    15941655  }
    15951656}
    … …  
    16001661}
    16011662
     1663vector<float> STGrid::getConvFunc()
     1664{
     1665  LogIO os(LogOrigin("STGrid","getConvFunc",WHERE));
     1666  Vector<Float> convFunc;
     1667  vector<float> out;
     1668
     1669  if (cellx_ <= 0.0 || celly_ <= 0.0) {
     1670    selectData();
     1671    setupGrid();
     1672  }
     1673
     1674  if (convType_ == "BOX" || convType_ == "SF") {
     1675    setConvFunc(convFunc);
     1676  }
     1677  else if (convType_ == "GAUSS") {
     1678    Quantum<Double> q1,q2;
     1679    readQuantity(q1,gwidth_);
     1680    readQuantity(q2,truncate_);
     1681//     if (celly_ <= 0.0
     1682//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
     1683//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel"))) {
     1684//       throw AipsError("You have to call defineImage to get correct convFunc");
     1685//     }
     1686    setConvFunc(convFunc);
     1687  }
     1688  else if (convType_ == "GJINC") {
     1689    Quantum<Double> q1,q2,q3;
     1690    readQuantity(q1,gwidth_);
     1691    readQuantity(q2,truncate_);
     1692    readQuantity(q3,jwidth_);
     1693//     if (celly_ <= 0.0
     1694//         && ((!q1.getUnit().empty()&&q1.getUnit()!="pixel") ||
     1695//             (!q2.getUnit().empty()&&q2.getUnit()!="pixel") ||
     1696//             (!q3.getUnit().empty()&&q3.getUnit()!="pixel"))) {
     1697//       throw AipsError("You have to call defineImage to get correct convFunc");
     1698//     }
     1699    setConvFunc(convFunc);
     1700  }
     1701  else if (convType_ == "PB") {
     1702    throw AipsError("Grid function PB is not available");
     1703  }
     1704  else {
     1705    throw AipsError("Unknown grid function: "+convType_);
     1706  }
     1707
     1708  convFunc.tovector(out);
     1709  return out;
     1710}
     1711
    16021712void STGrid::setConvFunc( Vector<Float> &convFunc )
    16031713{
     1714  LogIO os(LogOrigin("STGrid","setConvFunc",WHERE));
    16041715  convSupport_ = userSupport_ ;
    16051716  if ( convType_ == "BOX" ) {
    … …  
    16091720    convFunc.resize( convSize ) ;
    16101721    boxFunc( convFunc, convSize ) ;
     1722    os << LogIO::DEBUGGING
     1723       << "convType_ = " << convType_ << endl
     1724       << "convSupport_ = " << convSupport_ << LogIO::POST;
    16111725  }
    16121726  else if ( convType_ == "SF" ) {
    … …  
    16161730    convFunc.resize( convSize ) ;
    16171731    spheroidalFunc( convFunc ) ;
     1732    os << LogIO::DEBUGGING
     1733       << "convType_ = " << convType_ << endl
     1734       << "convSupport_ = " << convSupport_ << LogIO::POST;
    16181735  }
    16191736  else if ( convType_ == "GAUSS" ) {
    1620     // to take into account Gaussian tail
    1621     if ( convSupport_ < 0 )
    1622       convSupport_ = 4 ; // 1 * 4
    1623     else {
    1624       convSupport_ = userSupport_ * 4 ;
    1625     }
    1626     Int convSize = convSampling_ * ( 2 * convSupport_ + 2 ) ;
     1737    // determine pixel gwidth
     1738    // default is HWHM corresponding to b = 1.0 (Mangum et al. 2007)
     1739    Double pixelGW = sqrt(log(2.0));
     1740    Quantum<Double> q ;
     1741    if (!gwidth_.empty()) {
     1742      readQuantity( q, gwidth_ );
     1743      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
     1744        pixelGW = q.getValue();
     1745      }
     1746      else {
     1747        pixelGW = q.getValue("rad")/celly_;
     1748      }
     1749    }
     1750    // determine truncation radius
     1751    // default is 3 * HWHM
     1752    Double truncate = 3.0 * pixelGW;
     1753    if (!truncate_.empty()) {
     1754      readQuantity( q, truncate_ );
     1755      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
     1756        truncate = q.getValue();
     1757      }
     1758      else {
     1759        truncate = q.getValue("rad")/celly_;
     1760      }
     1761    }     
     1762    convSupport_ = (Int)(truncate+0.5);
     1763    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
    16271764    convFunc.resize( convSize ) ;
    1628     gaussFunc( convFunc ) ;
     1765    gaussFunc( convFunc, pixelGW, truncate ) ;
     1766    os << LogIO::DEBUGGING
     1767       << "convType_ = " << convType_ << endl
     1768       << "convSupport_ = " << convSupport_ << endl
     1769       << "truncate_ = " << truncate << endl
     1770       << "gwidth_ = " << pixelGW << "pixel" << LogIO::POST;
     1771  }
     1772  else if ( convType_ == "GJINC" ) {
     1773    // determine pixel gwidth
     1774    // default is HWHM corresponding to b = 2.52 (Mangum et al. 2007)
     1775    Double pixelGW = sqrt(log(2.0)) * 2.52;
     1776    Quantum<Double> q ;
     1777    if (!gwidth_.empty()) {
     1778      readQuantity( q, gwidth_ );
     1779      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
     1780        pixelGW = q.getValue();
     1781      }
     1782      else {
     1783        pixelGW = q.getValue("rad")/celly_;
     1784      }
     1785    }
     1786    // determine pixel c
     1787    // default is c = 1.55 (Mangum et al. 2007)
     1788    Double pixelJW = 1.55;
     1789    if (!jwidth_.empty()) {
     1790      readQuantity( q, jwidth_ );
     1791      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
     1792        pixelJW = q.getValue();
     1793      }
     1794      else {
     1795        pixelJW = q.getValue("rad")/celly_;
     1796      }
     1797    }     
     1798    // determine truncation radius
     1799    // default is -1.0 (truncate at first null)
     1800    Double truncate = -1.0;
     1801    if (!truncate_.empty()) {
     1802      readQuantity( q, truncate_ );
     1803      if ( q.getUnit().empty() || q.getUnit()=="pixel" ) {
     1804        truncate = q.getValue();
     1805      }
     1806      else {
     1807        truncate = q.getValue("rad")/celly_;
     1808      }
     1809    }     
     1810    convSupport_ = (truncate >= 0.0) ? (Int)(truncate+0.5) : (Int)(2*pixelJW+0.5);
     1811    Int convSize = convSampling_ * ( 2*convSupport_ + 2 ) ;
     1812    convFunc.resize( convSize ) ;
     1813    gjincFunc( convFunc, pixelGW, pixelJW, truncate ) ;
     1814    os << LogIO::DEBUGGING
     1815       << "convType_ = " << convType_ << endl
     1816       << "convSupport_ = " << convSupport_ << endl
     1817       << "truncate_ = " << truncate << endl
     1818       << "gwidth_ = " << pixelGW << "pixel" << endl
     1819       << "jwidth_ = " << pixelJW << "pixel" << LogIO::POST;
    16291820  }
    16301821  else if ( convType_ == "PB" ) {

  • trunk/src/STGrid.h

    r2669 r2671  
    6161                    string scenter="" ) ;
    6262  void setFunc( string convtype="box",
    63                 int convsupport=-1 ) ;
     63                int convsupport=-1,
     64                string truncate="",
     65                string gwidth="",
     66                string jwidth="" ) ;
    6467
    6568  void setWeight( const string wType="uniform" ) ;
    … …  
    7174 
    7275  string saveData( string outfile="" ) ;
     76
     77  // support function to know how grid function looks like
     78  vector<float> getConvFunc();
    7379
    7480//private:
    … …  
    138144  void boxFunc( Vector<Float> &convFunc, Int &convSize ) ;
    139145  void spheroidalFunc( Vector<Float> &convFunc ) ;
    140   void gaussFunc( Vector<Float> &convFunc ) ;
     146  void gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate ) ;
     147  void gjincFunc( Vector<Float> &convFunc, Double hwhm, Double c, Double truncate );
    141148  void pbFunc( Vector<Float> &convFunc ) ;
    142149  void setConvFunc( Vector<Float> &convFunc ) ;
    … …  
    232239  Int convSupport_ ;
    233240  Int userSupport_ ;
     241  String gwidth_;
     242  String jwidth_;
     243  String truncate_;
    234244  Int convSampling_ ;
    235245  Vector<uInt> pollist_ ;

  • trunk/src/python_STGrid.cpp

    r2593 r2671  
    2424    .def("_setpollist", &STGrid::setPolList)
    2525    .def("_setscanlist", &STGrid::setScanList)
    26     .def("_defineimage", &STGrid::defineImage)
    27     .def("_setfunc", &STGrid::setFunc)
     26    .def("_defineimage", &STGrid::defineImage,
     27         (boost::python::arg("nx")=-1,
     28          boost::python::arg("ny")=-1,
     29          boost::python::arg("scellx")="",
     30          boost::python::arg("scelly")="",
     31          boost::python::arg("scenter")=""))
     32    .def("_setfunc", &STGrid::setFunc,
     33         (boost::python::arg("convsupport")=-1,
     34          boost::python::arg("truncate")="",
     35          boost::python::arg("gwidth")="",
     36          boost::python::arg("jwidth")=""))
    2837    .def("_grid", &STGrid::grid)
    2938    .def("_setin", &STGrid::setFileIn)
    … …  
    3342    .def("_disableclip", &STGrid::disableClip)
    3443    .def("_save", &STGrid::saveData)
     44    .def("_getfunc", &STGrid::getConvFunc)
    3545    ;
    3646
    … …  
    4252    .def("_setpollist", &STGrid2::setPolList)
    4353    .def("_setscanlist", &STGrid2::setScanList)
    44     .def("_defineimage", &STGrid2::defineImage)
    45     .def("_setfunc", &STGrid2::setFunc)
     54    .def("_defineimage", &STGrid2::defineImage,
     55         (boost::python::arg("nx")=-1,
     56          boost::python::arg("ny")=-1,
     57          boost::python::arg("scellx")="",
     58          boost::python::arg("scelly")="",
     59          boost::python::arg("scenter")=""))
     60    .def("_setfunc", &STGrid2::setFunc,
     61         (boost::python::arg("convsupport")=-1,
     62          boost::python::arg("truncate")="",
     63          boost::python::arg("gwidth")="",
     64          boost::python::arg("jwidth")=""))
    4665    .def("_grid", &STGrid2::grid)
    4766    .def("_setin", &STGrid2::setScantable)
    … …  
    5170    .def("_disableclip", &STGrid2::disableClip)
    5271    .def("_get", &STGrid2::getResultAsScantable)
     72    .def("_getfunc", &STGrid2::getConvFunc)
    5373    ;
    5474   
Note: See TracChangeset for help on using the changeset viewer.