- Timestamp:
- 10/18/12 18:36:25 (12 years ago)
- Location:
- trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/python/asapgrid.py
r2669 r2671 135 135 self.gridder._defineimage( nx, ny, cellx, celly, center ) 136 136 137 def setFunc( self, func='box', width=-1 ):137 def setFunc( self, func='box', width=-1, gwidth="", jwidth="" ): 138 138 """ 139 139 Set convolution function. Possible options are 'box' (Box-car, … … 151 151 choose pre-defined value for each convolution function. 152 152 """ 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 ) 154 163 155 164 def setWeight( self, weightType='uniform' ): … … 223 232 asaplog.push('plot: elapsed time %s sec'%(t1-t0)) 224 233 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') 225 250 226 251 class asapgrid2: … … 367 392 choose pre-defined value for each convolution function. 368 393 """ 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 ) 370 404 371 405 def setWeight( self, weightType='uniform' ): … … 408 442 tp = 0 if rcParams['scantable.storage']=='memory' else 1 409 443 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') 410 460 411 461 class _SDGridPlotter: -
trunk/src/STGrid.cpp
r2669 r2671 83 83 convSupport_ = -1 ; 84 84 userSupport_ = -1 ; 85 truncate_ = ""; 86 gwidth_ = ""; 87 jwidth_ = ""; 85 88 convSampling_ = 100 ; 86 89 nprocessed_ = 0 ; … … 143 146 144 147 void STGrid::setFunc( string convType, 145 int convSupport ) 148 int convSupport, 149 string truncate, 150 string gwidth, 151 string jwidth ) 146 152 { 147 153 convType_ = String( convType ) ; 148 154 convType_.upcase() ; 149 155 userSupport_ = (Int)convSupport ; 156 truncate_ = String( truncate ); 157 gwidth_ = String( gwidth ); 158 jwidth_ = String( jwidth ); 150 159 } 151 160 … … 429 438 } 430 439 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 431 447 // grid parameter 432 448 os << LogIO::DEBUGGING ; … … 442 458 os << " center = " << center_ << endl ; 443 459 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 } 445 473 os << " doclip = " << (doclip_?"True":"False") << endl ; 446 474 os << "----------" << LogIO::POST ; … … 1567 1595 #if defined(NEED_UNDERSCORES) 1568 1596 #define grdsf grdsf_ 1597 #define grdgauss grdgauss_ 1598 #define grdjinc1 grdjinc1_ 1569 1599 #endif 1570 1600 extern "C" { 1571 void grdsf(Double*, Double*); 1601 void grdsf(Double*, Double*); 1602 void grdgauss(Double*, Double*, Double*); 1603 void grdjinc1(Double*, Double*, Int*, Double*); 1572 1604 } 1573 1605 void STGrid::spheroidalFunc( Vector<Float> &convFunc ) … … 1582 1614 } 1583 1615 1584 void STGrid::gaussFunc( Vector<Float> &convFunc )1616 void STGrid::gaussFunc( Vector<Float> &convFunc, Double hwhm, Double truncate ) 1585 1617 { 1586 1618 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; 1591 1621 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 1628 void 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 } 1594 1655 } 1595 1656 } … … 1600 1661 } 1601 1662 1663 vector<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 1602 1712 void STGrid::setConvFunc( Vector<Float> &convFunc ) 1603 1713 { 1714 LogIO os(LogOrigin("STGrid","setConvFunc",WHERE)); 1604 1715 convSupport_ = userSupport_ ; 1605 1716 if ( convType_ == "BOX" ) { … … 1609 1720 convFunc.resize( convSize ) ; 1610 1721 boxFunc( convFunc, convSize ) ; 1722 os << LogIO::DEBUGGING 1723 << "convType_ = " << convType_ << endl 1724 << "convSupport_ = " << convSupport_ << LogIO::POST; 1611 1725 } 1612 1726 else if ( convType_ == "SF" ) { … … 1616 1730 convFunc.resize( convSize ) ; 1617 1731 spheroidalFunc( convFunc ) ; 1732 os << LogIO::DEBUGGING 1733 << "convType_ = " << convType_ << endl 1734 << "convSupport_ = " << convSupport_ << LogIO::POST; 1618 1735 } 1619 1736 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 ) ; 1627 1764 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; 1629 1820 } 1630 1821 else if ( convType_ == "PB" ) { -
trunk/src/STGrid.h
r2669 r2671 61 61 string scenter="" ) ; 62 62 void setFunc( string convtype="box", 63 int convsupport=-1 ) ; 63 int convsupport=-1, 64 string truncate="", 65 string gwidth="", 66 string jwidth="" ) ; 64 67 65 68 void setWeight( const string wType="uniform" ) ; … … 71 74 72 75 string saveData( string outfile="" ) ; 76 77 // support function to know how grid function looks like 78 vector<float> getConvFunc(); 73 79 74 80 //private: … … 138 144 void boxFunc( Vector<Float> &convFunc, Int &convSize ) ; 139 145 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 ); 141 148 void pbFunc( Vector<Float> &convFunc ) ; 142 149 void setConvFunc( Vector<Float> &convFunc ) ; … … 232 239 Int convSupport_ ; 233 240 Int userSupport_ ; 241 String gwidth_; 242 String jwidth_; 243 String truncate_; 234 244 Int convSampling_ ; 235 245 Vector<uInt> pollist_ ; -
trunk/src/python_STGrid.cpp
r2593 r2671 24 24 .def("_setpollist", &STGrid::setPolList) 25 25 .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")="")) 28 37 .def("_grid", &STGrid::grid) 29 38 .def("_setin", &STGrid::setFileIn) … … 33 42 .def("_disableclip", &STGrid::disableClip) 34 43 .def("_save", &STGrid::saveData) 44 .def("_getfunc", &STGrid::getConvFunc) 35 45 ; 36 46 … … 42 52 .def("_setpollist", &STGrid2::setPolList) 43 53 .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")="")) 46 65 .def("_grid", &STGrid2::grid) 47 66 .def("_setin", &STGrid2::setScantable) … … 51 70 .def("_disableclip", &STGrid2::disableClip) 52 71 .def("_get", &STGrid2::getResultAsScantable) 72 .def("_getfunc", &STGrid2::getConvFunc) 53 73 ; 54 74
Note:
See TracChangeset
for help on using the changeset viewer.