1 #define COLORDEPTH 255
2 #define DEPTHEXTEN 32760
3
4 typedef double point[3];
5
6 typedef void (*BFunction)();
7 typedef double (*ZFunction)(double x, double y);
8 typedef void (*SFunction)(point p, double x, double y);
9 typedef void (*CFunction)(point p, double x);
10
11 const double CONV_RAD=M_PI/180;
12 int do_init = 1, orient = 1;
13 int width, height;
14
15 double x_1,x_2,y_1,y_2,z_1,z_2;
16 double ratio_x, ratio_y, ratio_z;
17 double light_x, light_y, light_z;
18 double m_rot[12]={1,0,0,0,0,1,0,0,0,0,1,0};
19
20 short int *z_depth = NULL;
21 short int *luminosity = NULL;
22
max(double x,double y)23 double max(double x, double y)
24 {
25 if(x>=y) return x; else return y;
26 }
27
min(double x,double y)28 double min(double x, double y)
29 {
30 if(x>=y) return x; else return y;
31 }
32
setcolor(char * pix,char * val)33 void setcolor(char *pix, char *val)
34 {
35 memcpy(pix, val, 3);
36 }
37
create_3d_buffer()38 void create_3d_buffer()
39 {
40 z_depth = (short int *)malloc(width*height*sizeof(short int));
41 luminosity = (short int *)malloc(width*height*sizeof(short int));
42 }
43
free_3d_buffer()44 void free_3d_buffer()
45 {
46 if (z_depth) free(z_depth);
47 if (luminosity) free(luminosity);
48 z_depth = NULL;
49 luminosity = NULL;
50 }
51
init_3d_buffer()52 void init_3d_buffer()
53 {
54 int i, j, k;
55 double light;
56
57 if (do_init) {
58 for (j=0; j<height; j++)
59 for (i=0; i<width; i++) {
60 k = i+j*width;
61 z_depth[k] = -DEPTHEXTEN;
62 luminosity[k] = 0;
63 }
64
65 light = sqrt(light_x*light_x + light_y*light_y + light_z*light_z);
66 if (light>0) {
67 light_x = light_x/light;
68 light_y = light_y/light;
69 light_z = light_z/light;
70 }
71 ratio_x = (x_2-x_1)/(double)width;
72 ratio_y = (y_2-y_1)/(double)height;
73 if(ratio_x==0 || ratio_y==0) exit(1);
74 ratio_z = z_2-z_1;
75 if (ratio_z<=0) exit(1); else ratio_z = ((double)DEPTHEXTEN)/ratio_z;
76 do_init = 0;
77 }
78 }
79
translation(double a,double b,double c)80 void translation(double a, double b, double c)
81 {
82 m_rot[3] = a;
83 m_rot[7] = b;
84 m_rot[11] = c;
85 }
86
scale(double a)87 void scale(double a)
88 {
89 int i;
90 for (i=0; i<=11; i++) m_rot[i] *= a;
91 }
92
rotation_angles(double theta,double phi,double psi)93 void rotation_angles(double theta, double phi, double psi)
94 {
95 double ca,cb,cc, sa,sb,sc;
96
97 ca=theta*CONV_RAD;
98 cb=phi*CONV_RAD;
99 cc=psi*CONV_RAD;
100
101 sa=sin(ca); ca=cos(ca);
102 sb=sin(cb); cb=cos(cb);
103 sc=sin(cc); cc=cos(cc);
104
105 m_rot[0] = ca*cb;
106 m_rot[1] = -sa*cc-ca*sb*sc;
107 m_rot[2] = sa*sc-ca*sb*cc;
108
109 m_rot[4] = sa*cb;
110 m_rot[5] = ca*cc-sa*sb*sc;
111 m_rot[6] = -ca*sc-sa*sb*cc;
112
113 m_rot[8] = sb;
114 m_rot[9] = cb*sc;
115 m_rot[10]= cb*cc;
116 }
117
rotate(point p,point q)118 void rotate(point p, point q)
119 {
120 int i,j;
121 for (i=0; i<=2; i++) {
122 j=4*i;
123 q[i] = m_rot[j]*p[0]+m_rot[j+1]*p[1]+m_rot[j+2]*p[2]+m_rot[j+3];
124 }
125 }
126
setpoint(point p,point q)127 void setpoint(point p, point q)
128 {
129 int i;
130 for(i=0; i<=2; i++) p[i] = q[i];
131 }
132
mesh(point * p,int h,int i,int j)133 void mesh(point *p, int h, int i, int j)
134 {
135 point a[3], b, c;
136 double ca,cb,cc,cd, s,t, adds,addt, x, y, z, z0,z1,z2;
137 double delta_x, delta_y, delta_z, delta;
138 int u1, u2, v1, v2, col;
139 int k,l,m,n;
140
141 init_3d_buffer();
142 rotate(p[h],a[0]);
143 rotate(p[i],a[1]);
144 rotate(p[j],a[2]);
145
146 u1=width; u2=-1;
147 v1=height; v2=-1;
148
149 for(k=0; k<=2; k++) {
150 l=(int)((a[k][0]-x_1)/ratio_x);
151 if(l<u1) u1=l;
152 if(l>u2) u2=l;
153 l=(int)((y_2-a[k][1])/ratio_y);
154 if(l<v1) v1=l;
155 if(l>v2) v2=l;
156 }
157
158 if(u1<0) u1=0;
159 if(u2>=width) u2=width-1;
160 if(v1<0) v1=0;
161 if(v2>=height) v2=height-1;
162 if(u1>u2 || v1>v2) return;
163 z0 = a[0][2];
164 z1 = a[1][2]-z0;
165 z2 = a[2][2]-z0;
166 z0 = (z0-z_1)*ratio_z;
167 z1 = z1*ratio_z;
168 z2 = z2*ratio_z;
169 for (k=0; k<=2; k++) {
170 b[k] = a[1][k] - a[0][k];
171 c[k] = a[2][k] - a[0][k];
172 }
173 delta_x = b[1]*c[2] - b[2]*c[1];
174 delta_y = -b[0]*c[2] + b[2]*c[0];
175 delta_z = b[0]*c[1] - b[1]*c[0];
176
177 if (fabs(delta_z)<1E-5) return;
178
179 delta = 1/delta_z;
180
181 ca = c[1] * delta;
182 cb = -c[0] * delta;
183 cc = -b[1] * delta;
184 cd = b[0] * delta;
185
186 adds = ca*ratio_x;
187 addt = cc*ratio_x;
188
189 delta = -1/(sqrt(delta_x*delta_x + delta_y*delta_y + delta_z*delta_z));
190 delta_x = delta_x * delta;
191 delta_y = delta_y * delta;
192 delta_z = delta_z * delta;
193
194 if (orient<0) {
195 /*
196 for (k=0; k<=2; k++)
197 c[k] = (a[0][k] + a[1][k] + a[2][k])/3.0 - m_rot[3+4*k];
198 delta = c[0] * delta_x + c[1] * delta_y + c[2] * delta_z;
199 if (delta>0) printf("+"); else printf("-");
200 if (delta<0) {
201 delta_x = -delta_x;
202 delta_y = -delta_y;
203 delta_z = -delta_z;
204 }
205 */
206 delta_x = -delta_x;
207 delta_y = -delta_y;
208 delta_z = -delta_z;
209 }
210 col = 32767 *
211 (1.000001 + light_x * delta_x + light_y*delta_y + light_z * delta_z);
212
213 for(l=v1; l<=v2; l++) {
214 x = x_1 + u1*ratio_x - a[0][0];
215 y = y_2 - l*ratio_y - a[0][1];
216 s = ca*x + cb*y;
217 t = cc*x + cd*y;
218 for(k=u1; k<=u2; k++) {
219 if (s>=-0.0001 && t>=-0.0001 && s+t<=1.0001) {
220 m = (int) (z0 + s*z1 + t*z2);
221 if (m<-DEPTHEXTEN) m=-DEPTHEXTEN;
222 if (m>DEPTHEXTEN) m=DEPTHEXTEN;
223 n = k+l*width;
224 if (m>z_depth[n]) {
225 z_depth[n] = m;
226 luminosity[n] = col;
227 }
228 }
229 s = s + adds;
230 t = t + addt;
231 }
232 }
233 }
234
build(BFunction bfunct)235 void build(BFunction bfunct)
236 {
237 init_3d_buffer();
238 bfunct();
239 }
240
light(ZFunction zfunct)241 void light(ZFunction zfunct)
242 {
243 int i, j;
244 int k, l;
245 double x,y;
246
247 init_3d_buffer();
248
249 for(j=0; j<height; j++) {
250 y = y_2-ratio_y*j;
251 for (i=0; i<width; i++) {
252 x = ratio_x*i+x_1;
253 k = (int)((zfunct(x,y)-z_1)*ratio_z);
254 if (k<0) k=0;
255 if (k>DEPTHEXTEN) k=DEPTHEXTEN;
256 l = i+j*width;
257 if (k>z_depth[l]) z_depth[l] = k;
258 }
259 }
260 }
261
surface(SFunction sfunct,double xi,double xs,int nx,double yi,double ys,int ny)262 void surface(SFunction sfunct, double xi, double xs, int nx,
263 double yi, double ys, int ny)
264 {
265 int i, j, l;
266 double hx,hy;
267 point p[4];
268
269 if (nx<=0 || ny<=0) return;
270
271 init_3d_buffer();
272
273 hx = (xs-xi)/nx;
274 hy = (ys-yi)/ny;
275
276 for (j=0; j<ny; j++) {
277 for (l=0; l<=1; l++)
278 sfunct(p[l],xi,yi+(j+l)*hy);
279
280 for (i=1; i<=nx; i++) {
281 setpoint(p[2],p[0]); setpoint(p[3],p[1]);
282 for (l=0; l<=1; l++)
283 sfunct(p[l], xi+i*hx,yi+(j+l)*hy);
284 mesh(p,0,1,2); mesh(p,1,3,2);
285 }
286 }
287 }
288
segment(point * p)289 void segment(point *p)
290 {
291 int i, u, x0, y0, x1, y1, z, delta;
292 int nx, ny, dx, dy, dxp, dyp, min, max;
293 double z0, dz;
294 point a[2];
295
296 init_3d_buffer();
297
298 for (i=0; i<=1; i++)
299 rotate(p[i], a[i]);
300
301 x0 = (int)((a[0][0] - x_1)/ratio_x);
302 y0 = (int)((a[0][1] - y_1)/ratio_y);
303 x1 = (int)((a[1][0] - x_1)/ratio_x);
304 y1 = (int)((a[1][1] - y_1)/ratio_y);
305 dz = ((double)DEPTHEXTEN)/(z_2-z_1);
306 z0 = (a[0][2] - z_1)*dz;
307 dz = (a[1][2] - a[0][2])*dz;
308
309 nx = x1-x0;
310 dxp = 1;
311 if (nx<0) {
312 dxp = -1;
313 nx = -nx;
314 }
315 ny = y1-y0;
316 dyp = 1;
317 if (ny<0) {
318 dyp = -1;
319 ny = -ny;
320 }
321 if (nx>ny) {
322 min = ny;
323 max = nx;
324 dx = dxp;
325 dy = 0;
326 } else {
327 min = nx;
328 max = ny;
329 dx = 0;
330 dy = dyp;
331 }
332
333 dz = dz/((double)max);
334
335 delta = max/2 - max;
336 for (i=0; i<=max; i++) {
337 u = x0+width*y0;
338 if (x0>=0 && x0<width && y0>=0 && y0<height) {
339 u = x0+width*y0;
340 z = (int) (z0 + ((double)i) * dz);
341 if (z>z_depth[u]) {
342 luminosity[u] = 1;
343 z_depth[u] = z;
344 }
345 }
346 delta += min;
347 if (delta>=max/2) {
348 delta -= max;
349 x0 += dxp;
350 y0 += dyp;
351 } else {
352 x0 += dx;
353 y0 += dy;
354 }
355 }
356 }
357
curve(CFunction cfunct,double ti,double ts,int n)358 void curve(CFunction cfunct, double ti, double ts, int n)
359 {
360 int i, l;
361 double h;
362 point p[2];
363
364 if (n<=0) return;
365
366 init_3d_buffer();
367
368 h = (ts-ti)/n;
369
370 for (i=0; i<=n; i++) {
371 if (i>=1) setpoint(p[0], p[1]);
372 cfunct(p[1], ti+i*h);
373 if (i>=1) segment(p);
374 }
375 }
376
377
378