yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #00743
[svn] r1559 - in trunk: . gui/py lib/smoothing
Author: eudoxos
Date: 2008-10-26 22:47:44 +0100 (Sun, 26 Oct 2008)
New Revision: 1559
Modified:
trunk/SConstruct
trunk/gui/py/_utils.cpp
trunk/lib/smoothing/WeightedAverage2d.cpp
trunk/lib/smoothing/WeightedAverage2d.hpp
Log:
!! IMPORTANT !! Introduced new dependency on "python-numpy" package
(numpy/ndarrayobject.h header), which is now checked for by scons.
Tell me if this is not OK!
1. Fix many things in the weighted average smoother, add polygon clipping.
2. Add utils.pointInsidePolygon to rapidly compue what it says; this is what
needs numpy.
Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct 2008-10-25 19:03:17 UTC (rev 1558)
+++ trunk/SConstruct 2008-10-26 21:47:44 UTC (rev 1559)
@@ -324,9 +324,10 @@
if not ok: featureNotOK('log4cxx')
env.Append(CPPDEFINES=['LOG4CXX'])
if 'python' in env['features']:
- ok=(conf.CheckPython() and conf.CheckIPython() # not needed now: and conf.CheckScientificPython()
- and (conf.CheckLibWithHeader('boost_python-mt','boost/python.hpp','c++','boost::python::scope();',autoadd=1)
- or conf.CheckLibWithHeader('boost_python','boost/python.hpp','c++','boost::python::scope();',autoadd=1)))
+ ok=(conf.CheckPython()
+ and conf.CheckIPython() # not needed now: and conf.CheckScientificPython()
+ and CheckLib_maybeMT(conf,'boost_python','boost/python.hpp','c++','boost::python::scope();')
+ and conf.CheckCXXHeader(['Python.h','numpy/ndarrayobject.h'],'<>'))
if not ok: featureNotOK('python')
env.Append(CPPDEFINES=['EMBED_PYTHON'])
if env['useMiniWm3']: env.Append(LIBS='miniWm3',CPPDEFINES=['MINIWM3'])
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2008-10-25 19:03:17 UTC (rev 1558)
+++ trunk/gui/py/_utils.cpp 2008-10-26 21:47:44 UTC (rev 1559)
@@ -6,6 +6,11 @@
#include<yade/pkg-dem/SpheresContactGeometry.hpp>
#include<yade/pkg-common/NormalShearInteractions.hpp>
#include<cmath>
+
+#include<numpy/ndarrayobject.h>
+
+
+
using namespace boost::python;
#ifdef LOG4CXX
@@ -199,6 +204,39 @@
return ret;
}
+/* Tell us whether a point lies in polygon given by array of points.
+ * @param xy is the point that is being tested
+ * @param vertices is Numeric.array (or list or tuple) of vertices of the polygon.
+ * Every row of the array is x and y coordinate, numer of rows is >= 3 (triangle).
+ *
+ * Copying the algorithm from http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
+ * is gratefully acknowledged:
+ *
+ * License to Use:
+ * Copyright (c) 1970-2003, Wm. Randolph Franklin
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
+ * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimers.
+ * 2. Redistributions in binary form must reproduce the above copyright notice in the documentation and/or other materials provided with the distribution.
+ * 3. The name of W. Randolph Franklin may not be used to endorse or promote products derived from this Software without specific prior written permission.
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ * http://numpy.scipy.org/numpydoc/numpy-13.html told me how to use Numeric.array from c
+ */
+bool pointInsidePolygon(python::tuple xy, python::object vertices){
+ Real testx=python::extract<double>(xy[0])(),testy=python::extract<double>(xy[1])();
+ char** vertData; int rows, cols; PyArrayObject* vert=(PyArrayObject*)vertices.ptr();
+ int result=PyArray_As2D((PyObject**)&vert /* is replaced */ ,&vertData,&rows,&cols,PyArray_DOUBLE);
+ if(result!=0) throw invalid_argument("Unable to cast vertices to 2d array");
+ if(cols!=2 || rows<3) throw invalid_argument("Vertices must have 2 columns (x and y) and at least 3 rows.");
+ int i /*current node*/, j/*previous node*/; bool inside=false;
+ for(i=0,j=rows-1; i<rows; j=i++){
+ double vx_i=*(double*)(vert->data+i*vert->strides[0]), vy_i=*(double*)(vert->data+i*vert->strides[0]+vert->strides[1]), vx_j=*(double*)(vert->data+j*vert->strides[0]), vy_j=*(double*)(vert->data+j*vert->strides[0]+vert->strides[1]);
+ if (((vy_i>testy)!=(vy_j>testy)) && (testx < (vx_j-vx_i) * (testy-vy_i) / (vy_j-vy_i) + vx_i) ) inside=!inside;
+ }
+ Py_DECREF(vert);
+ return inside;
+}
+
/* Project 3d point into 2d using spiral projection along given axis;
* the returned tuple is
*
@@ -218,7 +256,6 @@
Real hRef=dH_dTheta*(theta-theta0);
long period;
Real h=Shop::periodicWrap(pt[axis],hRef-.5*Mathr::PI*dH_dTheta,hRef+.5*Mathr::PI*dH_dTheta,&period);
- if(abs(h-hRef)>0.005) cerr<<"@@@ "<<h-hRef<<" "<<h<<" "<<hRef<<" "<<pt[axis]<<" ("<<hRef-.5*Mathr::PI*dH_dTheta<<","<<hRef+.5*Mathr::PI*dH_dTheta<<") "<<theta<<" "<<endl;
return python::make_tuple(python::make_tuple(r,h-hRef),theta);
}
BOOST_PYTHON_FUNCTION_OVERLOADS(spiralProject_overloads,spiralProject,2,4);
@@ -229,6 +266,9 @@
BOOST_PYTHON_FUNCTION_OVERLOADS(unbalancedForce_overloads,Shop::unbalancedForce,0,1);
BOOST_PYTHON_MODULE(_utils){
+ // http://numpy.scipy.org/numpydoc/numpy-13.html mentions this must be done in module init, otherwise we will crash
+ import_array();
+
def("PWaveTimeStep",PWaveTimeStep);
def("aabbExtrema",aabbExtrema,aabbExtrema_overloads(args("cutoff","centers")));
def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor")));
@@ -243,6 +283,7 @@
def("sumBexMoments",sumBexMoments);
def("createInteraction",Shop__createExplicitInteraction);
def("spiralProject",spiralProject,spiralProject_overloads(args("axis","theta0")));
+ def("pointInsidePolygon",pointInsidePolygon);
}
Modified: trunk/lib/smoothing/WeightedAverage2d.cpp
===================================================================
--- trunk/lib/smoothing/WeightedAverage2d.cpp 2008-10-25 19:03:17 UTC (rev 1558)
+++ trunk/lib/smoothing/WeightedAverage2d.cpp 2008-10-26 21:47:44 UTC (rev 1559)
@@ -1,5 +1,19 @@
#include"WeightedAverage2d.hpp"
+/* Tell whether point is inside polygon
+ *
+ * See _utils.cpp: pointInsidePolygon for docs and license.
+ */
+bool pyGaussAverage::pointInsidePolygon(const Vector2r& pt, const vector<Vector2r>& vertices){
+ int i /*current node*/, j/*previous node*/; bool inside=false; int rows=(int)vertices.size();
+ const Real& testx=pt[0],testy=pt[1];
+ for(i=0,j=rows-1; i<rows; j=i++){
+ const Real& vx_i=vertices[i][0], vy_i=vertices[i][1], vx_j=vertices[j][0], vy_j=vertices[j][1];
+ if (((vy_i>testy)!=(vy_j>testy)) && (testx < (vx_j-vx_i) * (testy-vy_i) / (vy_j-vy_i) + vx_i) ) inside=!inside;
+ }
+ return inside;
+}
+
BOOST_PYTHON_MODULE(WeightedAverage2d)
{
boost::python::class_<pyGaussAverage>("GaussAverage",python::init<python::tuple,python::tuple,python::tuple,Real>())
@@ -7,6 +21,7 @@
.def("avg",&pyGaussAverage::avg)
.add_property("stDev",&pyGaussAverage::stDev_get,&pyGaussAverage::stDev_set)
.add_property("relThreshold",&pyGaussAverage::relThreshold_get,&pyGaussAverage::relThreshold_set)
+ .add_property("clips",&pyGaussAverage::clips_get,&pyGaussAverage::clips_set)
;
};
Modified: trunk/lib/smoothing/WeightedAverage2d.hpp
===================================================================
--- trunk/lib/smoothing/WeightedAverage2d.hpp 2008-10-25 19:03:17 UTC (rev 1558)
+++ trunk/lib/smoothing/WeightedAverage2d.hpp 2008-10-26 21:47:44 UTC (rev 1559)
@@ -119,23 +119,23 @@
virtual Real getWeight(const Vector2r& meanPt, const T& e){
Vector2r pos=getPosition(e);
Real rSq=pow(meanPt[0]-pos[0],2)+pow(meanPt[1]-pos[1],2);
- if(rSq>relThreshold*stDev) return 0.;
- return (1./sqrt(2*M_PI*stDev))*exp(-rSq/(2*stDev));
+ if(rSq>relThreshold*stDev*stDev) return 0.;
+ return (1./sqrt(2*M_PI*stDev*stDev))*exp(-rSq/(2*stDev*stDev));
}
virtual vector<Vector2i> filterCells(const Vector2r& refPt){return WeightedAverage<T,Tvalue>::grid->circleFilter(refPt,stDev*relThreshold);}
};
-/* Class for doing template specialization of gaussian kernel average on SGKA_dataPt and for testing */
-struct dataPt{
+/* Class for doing template specialization of gaussian kernel average on SGKA_Scalar2d and for testing */
+struct Scalar2d{
Vector2r pos;
- int val;
+ Real val;
};
/* Final specialization; it only needs to define getValue and getPosition -- these functions contain knowledge about the element class itself */
-struct SGKA_dataPt: public SymmGaussianKernelAverage<dataPt,Real>{
- SGKA_dataPt(const shared_ptr<GridContainer<dataPt> >& _grid, Real _stDev, Real _relThreshold=3): SymmGaussianKernelAverage<dataPt,Real>(_grid,_stDev,_relThreshold){}
- virtual Real getValue(const dataPt& dp){return (Real)dp.val;}
- virtual Vector2r getPosition(const dataPt& dp){return dp.pos;}
+struct SGKA_Scalar2d: public SymmGaussianKernelAverage<Scalar2d,Real>{
+ SGKA_Scalar2d(const shared_ptr<GridContainer<Scalar2d> >& _grid, Real _stDev, Real _relThreshold=3): SymmGaussianKernelAverage<Scalar2d,Real>(_grid,_stDev,_relThreshold){}
+ virtual Real getValue(const Scalar2d& dp){return (Real)dp.val;}
+ virtual Vector2r getPosition(const Scalar2d& dp){return dp.pos;}
};
/* simplified interface for python:
@@ -148,19 +148,50 @@
*
* */
class pyGaussAverage{
- //struct dataPt{Vector2r pos; Real val;};
+ //struct Scalar2d{Vector2r pos; Real val;};
Vector2r tuple2vec2r(const python::tuple& t){return Vector2r(python::extract<Real>(t[0])(),python::extract<Real>(t[1])());}
Vector2i tuple2vec2i(const python::tuple& t){return Vector2i(python::extract<int>(t[0])(),python::extract<int>(t[1])());}
- shared_ptr<SGKA_dataPt> sgka;
+ shared_ptr<SGKA_Scalar2d> sgka;
+ struct Poly2d{vector<Vector2r> vertices; bool inclusive;};
+ vector<Poly2d> clips;
public:
pyGaussAverage(python::tuple lo, python::tuple hi, python::tuple nCells, Real stDev){
- shared_ptr<GridContainer<dataPt> > g(new GridContainer<dataPt>(tuple2vec2r(lo),tuple2vec2r(hi),tuple2vec2i(nCells)));
- sgka=shared_ptr<SGKA_dataPt>(new SGKA_dataPt(g,stDev));
+ shared_ptr<GridContainer<Scalar2d> > g(new GridContainer<Scalar2d>(tuple2vec2r(lo),tuple2vec2r(hi),tuple2vec2i(nCells)));
+ sgka=shared_ptr<SGKA_Scalar2d>(new SGKA_Scalar2d(g,stDev));
}
- void addPt(Real val, python::tuple pos){dataPt d; d.pos=tuple2vec2r(pos); d.val=val; sgka->grid->add(d,d.pos); }
- Real avg(python::tuple pt){return sgka->computeAverage(tuple2vec2r(pt));}
-
+ bool pointInsidePolygon(const Vector2r&,const vector<Vector2r>&);
+ bool ptIsClipped(const Vector2r& pt){
+ FOREACH(const Poly2d& poly, clips){
+ bool inside=pointInsidePolygon(pt,poly.vertices);
+ if((inside && !poly.inclusive) || (!inside && poly.inclusive)) return true;
+ }
+ return false;
+ }
+ bool addPt(Real val, python::tuple pos){Scalar2d d; d.pos=tuple2vec2r(pos); if(ptIsClipped(d.pos)) return false; d.val=val; sgka->grid->add(d,d.pos); return true; }
+ Real avg(python::tuple _pt){Vector2r pt=tuple2vec2r(_pt); if(ptIsClipped(pt)) return std::numeric_limits<Real>::quiet_NaN(); return sgka->computeAverage(pt);}
Real stDev_get(){return sgka->stDev;} void stDev_set(Real s){sgka->stDev=s;}
Real relThreshold_get(){return sgka->relThreshold;} void relThreshold_set(Real rt){sgka->relThreshold=rt;}
+ python::list clips_get(){
+ python::list ret;
+ FOREACH(const Poly2d& poly, clips){
+ python::list vertices;
+ FOREACH(const Vector2r& v, poly.vertices) vertices.append(python::make_tuple(v[0],v[1]));
+ ret.append(python::make_tuple(vertices,poly.inclusive));
+ }
+ return ret;
+ }
+ void clips_set(python::list l){
+ /* [ ( [(x1,y1),(x2,y2),…], true), … ] */
+ clips.clear();
+ for(int i=0; i<python::len(l); i++){
+ python::tuple polyDesc=python::extract<python::tuple>(l[i])();
+ python::list coords=python::extract<python::list>(polyDesc[0]);
+ Poly2d poly; poly.inclusive=python::extract<bool>(polyDesc[1]);
+ for(int j=0; j<python::len(coords); j++){
+ poly.vertices.push_back(tuple2vec2r(python::extract<python::tuple>(coords[j])));
+ }
+ clips.push_back(poly);
+ }
+ }
};
