关于多值隐函数绘图问题

forcal

有关于多值隐函数绘图的问题向大家请教！

以这个函数为例（x取-8～8，y取-8～8）：f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;

可以求得x为任意值时对应的多个y值。绘图时，将每一组值(x,y)绘制为一个点，点数多了就会形成曲线，如下图：

但如何将相邻的两个点连成线呢？就是说如何判断两个点是相邻的呢？

例如，以下代码求得了x＝1.1及x＝1.2时对应的y值；x＝1.1的y值有6个，x＝1.2的y值也有6个，如何判断哪些点该连成线呢？

1. !using("fcopt");
   
2. f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;
   
3. isolve[HFor("f"), optrange,-8,8, optthis,1, optpara,1.1];
   
4. isolve[HFor("f"), optrange,-8,8, optthis,1, optpara,1.2];

复制代码

1.1                       -5.413569621185505        0.
1.1                       -1.52029233271966         0.
1.1                       -0.8671650244626143       8.881784197001252e-016
1.1                       0.8671650244626143        8.881784197001252e-016
1.1                       1.52029233271966          0.
1.1                       5.413569621185504         0.
6.
1.2                       -5.274704716901752        0.
1.2                       -1.731609595222745        0.
1.2                       -0.9268324331975657       0.
1.2                       0.9268324331975657        0.
1.2                       1.731609595222745         0.
1.2                       5.274704716901752         0.
6.

wayne

1# forcal
点画多了,越稠密,就越连续了,
================

或者求出一定数目的点,然后线性插值,样条插值,...............

wayne

不知道你具体是在问什么

是关于画图引擎的渲染技巧吗
还是就题论题,只是画出该函数的 点图,然后标记一条直线?
===================================

前者,你可以参考gnuplot
如果是后者,你可以把方程化成极坐标的形式,就清晰了

forcal

不知道你具体是在问什么

是关于画图引擎的渲染技巧吗
还是就题论题,只是画出该函数的 点图,然后标记一条直线?
===================================

前者,你可以参考gnuplot
如果是后者,你可以把方程化成极 ...
wayne 发表于 2011-1-13 17:08

画出函数的点图是比较容易的，我想把点连成一条折线（曲线）。
例如matlab等画这个图的曲线是连续的，但并不是点图。

＝＝＝＝＝＝＝＝＝
有些隐函数方程化成极坐标的形式是比较困难的。

wayne

4# forcal

例如matlab等画这个图的曲线是连续的，但并不是点图。

那你做一个测试,试试:
先产生一个点序列，100个 ，,比如   {x,sin(x)}
然后shuffle一下，把这100对的顺序全部搞乱，越乱越好，再plot一下，看看，..................

wayne

4# forcal
如果你说的是contour plot,
那个实现起来要复杂一些, 你用google搜搜 contour plot Implementation

三天后我会回来的

wayne

4# forcal
算了,再给你补一贴吧:

用OpenGL 画你给的那个函数:


代码:

1. /* generates contours using marching squares */
   
2. 
   
3. /* region size */
   
4. 
   
5. #define X_MAX 8.0
   
6. #define Y_MAX 8.0
   
7. #define X_MIN -8.0
   
8. #define Y_MIN -8.0
   
9. 
   
10. /* number of cells */
    
11. 
    
12. #define N_X 500
    
13. #define N_Y 500
    
14. 
    
15. /* contour value */
    
16. 
    
17. #define THRESHOLD 0.0
    
18. 
    
19. #include <GL/glut.h>
    
20. 
    
21. void display()
    
22. {
    
23.    double f(double,double);
    
24.    int cell(double, double, double, double);
    
25.    void lines(int, int, int, double, double, double, double);
    
26. 
    
27.    double data[N_X][N_Y];
    
28.    int i,j;
    
29.    int c;
    
30. 
    
31.    glClear(GL_COLOR_BUFFER_BIT);
    
32. /* form data array from function */
    
33. 
    
34.     for(i=0;i<N_X;i++) for (j=0;j<N_Y;j++)
    
35.        data[i][j]=f(X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0), Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0));
    
36. 
    
37.     /* process each cell */
    
38. 
    
39.     for(i=0;i<N_X;i++) for (j=0;j<N_Y;j++)
    
40.     {
    
41.        c=cell(data[i][j], data[i+1][j], data[i+1][j+1], data[i][j+1]);
    
42.        lines(c,i,j,data[i][j], data[i+1][j], data[i+1][j+1], data[i][j+1]);
    
43.     }
    
44.     glFlush();
    
45. }
    
46. 
    
47. /* define function f(x,y)      */
    
48. 
    
49. double f(double x, double y)
    
50. {
    
51. 
    
52.    double a=36;
    
53. 
    
54.    /* Ovals of Cassini  */
    
55. 
    
56.    return (x*x+y*y)*(x*x+y*y)*(x*x+y*y)-a*(x*x-y*y)*(x*x-y*y);
    
57. }
    
58. 
    
59. 
    
60. /* define cell vertices */
    
61. 
    
62. int cell(double a, double b, double c , double d)
    
63. {
    
64.     int n=0;
    
65.     if(a>THRESHOLD) n+=1;
    
66.     if(b>THRESHOLD) n+=8;
    
67.     if(c>THRESHOLD) n+=4;
    
68.     if(d>THRESHOLD) n+=2;
    
69.     return n;
    
70. }
    
71. 
    
72. 
    
73. /* draw line segments for each case */
    
74. 
    
75. void lines(int num, int i, int j, double a, double b, double c, double d)
    
76. {
    
77.      void draw_one(int, int, int, double, double, double, double);
    
78.      void draw_adjacent(int, int, int, double, double, double, double);
    
79.      void draw_opposite(int, int, int, double, double, double, double);
    
80.      switch(num)
    
81.      {
    
82.      case 1: case 2: case 4: case 7: case 8: case 11: case 13: case 14:
    
83.         draw_one(num, i,j,a,b,c,d);
    
84.         break;
    
85.      case 3: case 6: case 9: case 12:
    
86.         draw_adjacent(num,i,j,a,b,c,d);
    
87.         break;
    
88.      case 5: case 10:
    
89.         draw_opposite(num, i,j,a,b,c,d);
    
90.         break;
    
91.      case 0: case 15: break;
    
92.      }
    
93. }
    
94. 
    
95. void draw_one(int num, int i, int j, double a, double b, double c, double d)
    
96. {
    
97.   double x1, y1, x2, y2;
    
98.   double ox, oy;
    
99.   double dx, dy;
    
100.   dx=(X_MAX-(X_MIN))/(N_X-1.0);
     
101.   dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
     
102.   ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
     
103.   oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
     
104.     switch(num)
     
105.     {
     
106.     case 1: case 14:
     
107.       x1=ox;
     
108.       y1=oy+dy*(THRESHOLD-a)/(d-a);
     
109.       x2=ox+dx*(THRESHOLD-a)/(b-a);
     
110.       y2=oy;
     
111.       break;
     
112.     case 2: case 13:
     
113.       x1=ox;
     
114.       y1=oy+dy*(THRESHOLD-a)/(d-a);
     
115.       x2=ox+dx*(THRESHOLD-d)/(c-d);
     
116.       y2=oy+dy;
     
117.       break;
     
118.     case 4: case 11:
     
119.       x1=ox+dx*(THRESHOLD-d)/(c-d);
     
120.       y1=oy+dy;
     
121.       x2=ox+dx;
     
122.       y2=oy+dy*(THRESHOLD-b)/(c-b);
     
123.       break;
     
124.     case 7: case 8:
     
125.       x1=ox+dx*(THRESHOLD-a)/(b-a);
     
126.       y1=oy;
     
127.       x2=ox+dx;
     
128.       y2=oy+dy*(THRESHOLD-b)/(c-b);
     
129.       break;
     
130.     }
     
131.   glBegin(GL_LINES);
     
132.     glVertex2d(x1, y1);
     
133.     glVertex2d(x2, y2);
     
134.   glEnd();
     
135. }
     
136. 
     
137. void draw_adjacent(int num, int i, int j, double a, double b,
     
138. double c, double d)
     
139. {
     
140.   double x1, y1, x2, y2;
     
141.   double ox, oy;
     
142.   double dx, dy;
     
143.   dx=(X_MAX-(X_MIN))/(N_X-1.0);
     
144.   dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
     
145.   ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
     
146.   oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
     
147.   switch(num)
     
148.   {
     
149.     case 3: case 12:
     
150.       x1=ox+dx*(THRESHOLD-a)/(b-a);
     
151.       y1=oy;
     
152.       x2=ox+dx*(THRESHOLD-d)/(c-d);
     
153.       y2=oy+dy;
     
154.       break;
     
155.     case 6: case 9:
     
156.       x1=ox;
     
157.       y1=oy+dy*(THRESHOLD-a)/(d-a);
     
158.       x2=ox+dx;
     
159.       y2=oy+dy*(THRESHOLD-b)/(c-b);
     
160.       break;
     
161.   }
     
162.   glBegin(GL_LINES);
     
163.     glVertex2d(x1, y1);
     
164.     glVertex2d(x2, y2);
     
165.   glEnd();
     
166. 
     
167. }
     
168. 
     
169. void draw_opposite(int num, int i, int j, double a, double b,
     
170. double c, double d)
     
171. {
     
172.   double x1,y1,x2,y2,x3,y3,x4,y4;
     
173.   double ox, oy;
     
174.   double dx, dy;
     
175.   dx=(X_MAX-(X_MIN))/(N_X-1.0);
     
176.   dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
     
177.   ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
     
178.   oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
     
179.   switch(num)
     
180.   {
     
181.     case5:
     
182.       x1=ox;
     
183.       y1=oy+dy*(THRESHOLD-a)/(d-a);
     
184.       x2=ox+dx*(THRESHOLD-a)/(b-a);
     
185.       y2=oy;
     
186.       x3=ox+dx*(THRESHOLD-d)/(c-d);
     
187.       y3=oy+dy;
     
188.       x4=ox+dx;
     
189.       y4=oy+dy*(THRESHOLD-b)/(c-b);
     
190.       break;
     
191.     case 10:
     
192.       x1=ox;
     
193.       y1=oy+dy*(THRESHOLD-a)/(d-a);
     
194.       x2=ox+dx*(THRESHOLD-d)/(c-d);
     
195.       y2=oy+dy;
     
196.       x3=ox+dy*(THRESHOLD-a)/(b-a);
     
197.       y3=oy;
     
198.       x4=ox+dx;
     
199.       y4=oy+dy*(THRESHOLD-b)/(c-b);
     
200.       break;
     
201.   }
     
202.   glBegin(GL_LINES);
     
203.     glVertex2d(x1, y1);
     
204.     glVertex2d(x2, y2);
     
205.     glVertex2d(x3, y3);
     
206.     glVertex2d(x4, y4);
     
207.   glEnd();
     
208. }
     
209. 
     
210. void myReshape(int w, int h)
     
211. {
     
212.     glViewport(0, 0, w, h);
     
213.     glMatrixMode(GL_PROJECTION);
     
214.     glLoadIdentity();
     
215.     if (w <= h) gluOrtho2D(X_MIN, X_MAX, Y_MIN * (GLfloat) h / (GLfloat) w, Y_MAX* (GLfloat) h / (GLfloat) w);
     
216.     else gluOrtho2D(X_MIN * (GLfloat) w / (GLfloat) h, X_MAX * (GLfloat) w / (GLfloat) h, Y_MIN, Y_MAX);
     
217.     glMatrixMode(GL_MODELVIEW);
     
218. }
     
219. 
     
220. int main(int argc, char **argv)
     
221. {
     
222.    glutInit(&argc, argv);
     
223.    glutInitWindowSize(500, 500);
     
224.    glutCreateWindow("contour plot");
     
225.    glutReshapeFunc(myReshape);
     
226.    glutDisplayFunc(display);
     
227.    glClearColor(0.0,0.0,0.0,1.0);
     
228.    glColor3f(1.0,1.0,1.0);
     
229.    glutMainLoop();
     
230. }
     
231. 
     

复制代码

forcal

4# forcal


那你做一个测试,试试:
先产生一个点序列，100个 ，,比如   {x,sin(x)}
然后shuffle一下，把这100对的顺序全部搞乱，越乱越好，再plot一下，看看，..................
wayne 发表于 2011-1-13 18:56

我的意思是MATLAB提供了一个ezplot函数绘制隐函数图形，似乎不是用描点方式绘图的，它是如何绘制的？

对于单值函数f(x)=... ...；绘制曲线图是很容易的，将所有点按x的次序连起来即可。

但对于隐函数f(x,y)=0，一个x对应多个y值，该如何绘制曲线图？

forcal

4# forcal
如果你说的是contour plot,
那个实现起来要复杂一些, 你用google搜搜 contour plot Implementation

三天后我会回来的
wayne 发表于 2011-1-13 19:10

英文不是很好，等你回来，呵呵。

forcal

4# forcal
算了,再给你补一贴吧:

用OpenGL 画你给的那个函数:
wayne 发表于 2011-1-13 20:24

这个绘图原理与我所说的不一样，我原想先解出所有要绘制的点，然后连成线，看来比较复杂。

这个OpenGL 绘隐函数图原理能大致说一下吗？看不太懂。