关于多值隐函数绘图问题

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 */
3. /* region size */
5. #define X_MAX 8.0
6. #define Y_MAX 8.0
7. #define X_MIN -8.0
8. #define Y_MIN -8.0
10. /* number of cells */
12. #define N_X 500
13. #define N_Y 500
15. /* contour value */
17. #define THRESHOLD 0.0
19. #include <GL/glut.h>
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);
27. double data[N_X][N_Y];
28. int i,j;
29. int c;
31. glClear(GL_COLOR_BUFFER_BIT);
32. /* form data array from function */
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));
37. /* process each cell */
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. }
47. /* define function f(x,y) */
49. double f(double x, double y)
50. {
52. double a=36;
54. /* Ovals of Cassini */
56. return (x*x+y*y)*(x*x+y*y)*(x*x+y*y)-a*(x*x-y*y)*(x*x-y*y);
57. }
60. /* define cell vertices */
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. }
73. /* draw line segments for each case */
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. }
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. }
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();
167. }
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. }
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. }
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. }

复制代码

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 绘隐函数图原理能大致说一下吗？看不太懂。