Menu
▾
▴
[Plplot-cvs] SF.net SVN: plplot:[11594] trunk/examples
|
From: <ai...@us...> - 2011-03-04 21:03:07
|
Revision: 11594
http://plplot.svn.sourceforge.net/plplot/?rev=11594&view=rev
Author: airwin
Date: 2011-03-04 21:03:00 +0000 (Fri, 04 Mar 2011)
Log Message:
-----------
Update C and python versions of example 27 so that fills work correctly
both for -eofill (even-odd fill) and with that option dropped
(non-zero winding rule fill).
All the previous fill issues for self-intersecting boundaries were an
artifact of the incorrect angle termination condition (via windings).
The correct windings must be calculated such that _both_ phi and phiw
are the minimum possible integer multiple of 2 pi. It is easy to show
that this condition leads to windings = params[1]/gcd(params[0],
params[1]) (with appropriate adjustments to convert to integers and
use absolute values).
With this change to the maximum phi and phiw, the results look
beautiful for both -eofill and not.
I ascribe the previous bad filling results for the -eofill case as due
to the boundaries of certain areas being duplicated an even number of
times (which reverses the decision about whether a given area will be
filled or not for the even-odd fill rule). I ascribe the previous
slightly bad fill results for the non-zero winding number fill rule
to the last large boundary segment that was required to close the
boundary for the previous case.
A side benefit of no longer using an excessive number of windings is
that NPTS can be divided by 10 and yet still have a smooth result.
I also changed the xmin, xmax, ymin, ymax logic so that it worked
reliably for all cases.
For the python case I also took advantage of numpy array capabilities.
The resulting C and python results are consistent. Please propagate
these windings, NPT, and min/max changes to all other languages.
Modified Paths:
--------------
trunk/examples/c/x27c.c
trunk/examples/python/xw27.py
Modified: trunk/examples/c/x27c.c
===================================================================
--- trunk/examples/c/x27c.c 2011-03-03 22:46:14 UTC (rev 11593)
+++ trunk/examples/c/x27c.c 2011-03-04 21:03:00 UTC (rev 11594)
@@ -28,6 +28,7 @@
void cycloid( void );
void spiro( PLFLT data[], int fill );
+PLINT gcd (PLINT a, PLINT b);
//--------------------------------------------------------------------------
// main
@@ -41,6 +42,10 @@
main( int argc, const char *argv[] )
{
// R, r, p, N
+ // R and r should be integers to give correct termination of the
+ // angle loop using gcd.
+ // N.B. N is just a place holder since it is no longer used
+ // (because we now have proper termination of the angle loop).
PLFLT params[9][4] = {
21.0, 7.0, 7.0, 3.0, // Deltoid
21.0, 7.0, 10.0, 3.0,
@@ -113,7 +118,25 @@
}
//--------------------------------------------------------------------------
+// Calculate greatest common divisor following pseudo-code for the
+// Euclidian algorithm at http://en.wikipedia.org/wiki/Euclidean_algorithm
+PLINT gcd (PLINT a, PLINT b)
+{
+ PLINT t;
+ a = abs(a);
+ b = abs(b);
+ while ( b != 0 )
+ {
+ t = b;
+ b = a % b;
+ a = t;
+ }
+ return a;
+}
+
+//--------------------------------------------------------------------------
+
void
cycloid( void )
{
@@ -125,7 +148,7 @@
void
spiro( PLFLT params[], int fill )
{
-#define NPNT 20000
+#define NPNT 2000
static PLFLT xcoord[NPNT + 1];
static PLFLT ycoord[NPNT + 1];
@@ -139,26 +162,29 @@
PLFLT xmax;
PLFLT ymin;
PLFLT ymax;
- PLFLT scale;
// Fill the coordinates
- windings = (int) params[3];
+ // Proper termination of the angle loop very near the beginning
+ // point, see
+ // http://mathforum.org/mathimages/index.php/Hypotrochoid.
+ windings = (PLINT) abs(params[1])/gcd((PLINT) params[0], (PLINT) params[1]);
steps = NPNT / windings;
- dphi = 8.0 * acos( -1.0 ) / (PLFLT) steps;
+ dphi = 2.0 * PI / (PLFLT) steps;
- xmin = 0.0; // This initialisation is safe!
- xmax = 0.0;
- ymin = 0.0;
- ymax = 0.0;
-
for ( i = 0; i <= windings * steps; i++ )
{
phi = (PLFLT) i * dphi;
phiw = ( params[0] - params[1] ) / params[1] * phi;
xcoord[i] = ( params[0] - params[1] ) * cos( phi ) + params[2] * cos( phiw );
ycoord[i] = ( params[0] - params[1] ) * sin( phi ) - params[2] * sin( phiw );
-
+ if ( i == 0)
+ {
+ xmin = xcoord[i];
+ xmax = xcoord[i];
+ ymin = ycoord[i];
+ ymax = ycoord[i];
+ }
if ( xmin > xcoord[i] )
xmin = xcoord[i];
if ( xmax < xcoord[i] )
@@ -169,18 +195,10 @@
ymax = ycoord[i];
}
- if ( xmax - xmin > ymax - ymin )
- {
- scale = xmax - xmin;
- }
- else
- {
- scale = ymax - ymin;
- }
- xmin = -0.65 * scale;
- xmax = 0.65 * scale;
- ymin = -0.65 * scale;
- ymax = 0.65 * scale;
+ xmin -= 0.15 * (xmax - xmin);
+ xmax += 0.15 * (xmax - xmin);
+ ymin -= 0.15 * (ymax - ymin);
+ ymax += 0.15 * (ymax - ymin);
plwind( xmin, xmax, ymin, ymax );
Modified: trunk/examples/python/xw27.py
===================================================================
--- trunk/examples/python/xw27.py 2011-03-03 22:46:14 UTC (rev 11593)
+++ trunk/examples/python/xw27.py 2011-03-04 21:03:00 UTC (rev 11594)
@@ -23,6 +23,7 @@
#
from plplot_py_demos import *
+import types
# main
#
@@ -32,7 +33,11 @@
def main():
- # R, r, p, N
+ # R, r, p, N
+ # R and r should be integers to give correct termination of the
+ # angle loop using gcd.
+ # N.B. N is just a place holder since it is no longer used
+ # (because we now have proper termination of the angle loop).
params = [ [21.0, 7.0, 7.0, 3.0], # Deltoid
[21.0, 7.0, 10.0, 3.0],
[21.0, -7.0, 10.0, 3.0],
@@ -74,48 +79,41 @@
plvpor( 0.0, 1.0, 0.0, 1.0 )
spiro( params[i], 1 )
+def gcd(a, b):
+ if not (type(a) is types.IntType and type(b) is types.IntType):
+ raise RuntimeError, "gcd arguments must be integers"
+ a = abs(a);
+ b = abs(b);
+ while(b != 0):
+ t = b
+ b = a % b
+ a = t
+ return a
def spiro(params, fill):
# Fill the coordinates
- NPNT = 20000
+ NPNT = 2000
- windings = int(params[3])
+ # Proper termination of the angle loop very near the beginning
+ # point, see
+ # http://mathforum.org/mathimages/index.php/Hypotrochoid.
+ windings = int(abs(params[1])/gcd(int(params[0]), int(params[1])))
steps = int(NPNT/windings)
- dphi = 8.0*arccos(-1.0)/steps
+ dphi = 2.0*pi/float(steps)
- xmin = 0.0 # This initialisation is safe!
- xmax = 0.0
- ymin = 0.0
- ymax = 0.0
-
- # Add 0. to convert to real array for Numeric. numpy does not require this.
- xcoord = 0. + zeros(windings*steps+1)
- ycoord = 0. + zeros(windings*steps+1)
-
- for i in range(windings*steps+1) :
- phi = i * dphi
- phiw = (params[0]-params[1])/params[1]*phi
- xcoord[i] = (params[0]-params[1])*cos(phi) + params[2]*cos(phiw)
- ycoord[i] = (params[0]-params[1])*sin(phi) - params[2]*sin(phiw)
-
- if ( xmin > xcoord[i] ) :
- xmin = xcoord[i]
- if ( xmax < xcoord[i] ) :
- xmax = xcoord[i]
- if ( ymin > ycoord[i] ) :
- ymin = ycoord[i]
- if ( ymax < ycoord[i] ) :
- ymax = ycoord[i]
-
- if ( xmax-xmin > ymax-ymin ) :
- scale = xmax - xmin
- else :
- scale = ymax - ymin
+ phi = arange(windings*steps+1)*dphi
+ phiw = (params[0]-params[1])/params[1]*phi
+ xcoord = (params[0]-params[1])*cos(phi) + params[2]*cos(phiw)
+ ycoord = (params[0]-params[1])*sin(phi) - params[2]*sin(phiw)
+ xmin = min(xcoord)
+ xmax = max(xcoord)
+ ymin = min(ycoord)
+ ymax = max(ycoord)
- xmin = - 0.65 * scale
- xmax = 0.65 * scale
- ymin = - 0.65 * scale
- ymax = 0.65 * scale
+ xmin -= 0.15 * (xmax - xmin)
+ xmax += 0.15 * (xmax - xmin)
+ ymin -= 0.15 * (ymax - ymin)
+ ymax += 0.15 * (ymax - ymin)
plwind( xmin, xmax, ymin, ymax )
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
