1 /***************************************************************************
2 * pixel.cpp is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 #include <algorithm>
21 #include "mgl2/canvas.h"
22 #include "mgl2/thread.h"
23 #if MGL_HAVE_OMP
24 #include <omp.h>
25 #endif
26
27 //-----------------------------------------------------------------------------
pxl_combine(long id,long n,const void *)28 void mglCanvas::pxl_combine(long id, long n, const void *)
29 {
30 #if !MGL_HAVE_PTHREAD
31 #pragma omp parallel for
32 #endif
33 for(long i=id;i<n;i+=mglNumThr)
34 {
35 unsigned char *cc = C+12*i, c[4], *b=GB+4*i, *g=G4+4*i;
36 c[0]=b[0]; c[1]=b[1]; c[2]=b[2]; c[3]=b[3];
37 combine(c,cc+8); combine(c,cc+4); combine(c,cc);
38 g[0]=c[0]; g[1]=c[1]; g[2]=c[2]; g[3]=c[3];
39 }
40 }
41 //-----------------------------------------------------------------------------
pxl_memcpy(long id,long n,const void *)42 void mglCanvas::pxl_memcpy(long id, long n, const void *)
43 {
44 #if !MGL_HAVE_PTHREAD
45 #pragma omp parallel for
46 #endif
47 for(long i=id;i<n;i+=mglNumThr)
48 { unsigned char *g=G4+4*i, *c=C+12*i;
49 g[0]=c[0]; g[1]=c[1]; g[2]=c[2]; g[3]=c[3]; }
50 }
51 //-----------------------------------------------------------------------------
pxl_backgr(long id,long n,const void *)52 void mglCanvas::pxl_backgr(long id, long n, const void *)
53 {
54 #if !MGL_HAVE_PTHREAD
55 #pragma omp parallel for
56 #endif
57 for(long i=id;i<n;i+=mglNumThr)
58 { unsigned char *b=GB+4*i, c[4]={BDef[0],BDef[1],BDef[2],BDef[3]}, *g=G+3*i;
59 combine(c,G4+4*i); g[0]=c[0]; g[1]=c[1]; g[2]=c[2]; }
60 }
61 //-----------------------------------------------------------------------------
pxl_transform(long id,long n,const void *)62 void mglCanvas::pxl_transform(long id, long n, const void *)
63 {
64 const float *b = Bp.b;
65 const float dx = -Bp.x*Width/2, dy = -Bp.y*Height/2, dz = Depth/2.;
66 #if !MGL_HAVE_PTHREAD
67 #pragma omp parallel for
68 #endif
69 for(long i=id;i<n;i+=mglNumThr)
70 {
71 mglPnt &p=Pnt[i];
72 if(p.sub>=0)
73 {
74 float x = p.xx-Width/2., y = p.yy-Height/2., z = p.zz-Depth/2.;
75 p.x = b[0]*x + b[1]*y + b[2]*z + dx;
76 p.y = b[3]*x + b[4]*y + b[5]*z + dy;
77 p.z = b[6]*x + b[7]*y + b[8]*z + dz;
78 float d = get_persp(Bp.pf,p.z,Depth);
79 p.x = Width/2. + d*p.x; p.y = Height/2. + d*p.y;
80 }
81 }
82 }
83 //-----------------------------------------------------------------------------
pxl_setz_adv(long id,long n,const void *)84 void mglCanvas::pxl_setz_adv(long id, long n, const void *)
85 {
86 #if !MGL_HAVE_PTHREAD
87 #pragma omp parallel for
88 #endif
89 for(long i=id;i<n;i+=mglNumThr)
90 {
91 mglPrim &q=Prm[i]; q.z = Pnt[q.n1].z;
92 if(q.type==3) q.z = (q.z + Pnt[q.n2].z + Pnt[q.n3].z + Pnt[q.n4].z)/4;
93 else if(q.type==2) q.z = (q.z + Pnt[q.n2].z + Pnt[q.n3].z)/3;
94 else if(q.type==1) q.z = (q.z + Pnt[q.n2].z)/2;
95 }
96 }
97 //-----------------------------------------------------------------------------
pxl_pntcol(long id,long n,const void *)98 void mglCanvas::pxl_pntcol(long id, long n, const void *)
99 {
100 #if !MGL_HAVE_PTHREAD
101 #pragma omp parallel for
102 #endif
103 for(long i=id;i<n;i+=mglNumThr)
104 { mglRGBA c; col2int(Pnt[i],c.r,HighId-1); pnt_col[i]=c.c; }
105 }
106 //-----------------------------------------------------------------------------
pxl_setz(long id,long n,const void *)107 void mglCanvas::pxl_setz(long id, long n, const void *)
108 {
109 #if !MGL_HAVE_PTHREAD
110 #pragma omp parallel for
111 #endif
112 for(long i=id;i<n;i+=mglNumThr)
113 { mglPrim &q=Prm[i]; q.z = Pnt[q.n1].z; }
114 }
115 //-----------------------------------------------------------------------------
pxl_primdr(long id,long,const void *)116 void mglCanvas::pxl_primdr(long id, long , const void *)
117 {
118 #define Q 4 // should be >= sqrt(2*num_thr) ???
119 const int nx=Q,ny=Q; // TODO find dependence on Q for 1, 2, 4, 8 threads. Try to select optimal
120 #if !MGL_HAVE_PTHREAD
121 #pragma omp parallel for
122 #endif
123 for(long i=id;i<nx*ny;i+=mglNumThr)
124 {
125 mglDrawReg d; d.set(this,nx,ny,i);
126 if(Quality&MGL_DRAW_NORM) for(size_t k=0;k<Prm.size();k++)
127 {
128 if(Stop) break;
129 const mglPrim &p=GetPrm(k); d.copy(p);
130 long n1=p.n1, n2=p.n2, n3=p.n3, n4=p.n4;
131 switch(p.type)
132 {
133 case 3: quad_draw(Pnt[n1],Pnt[n2],Pnt[n3],Pnt[n4],&d); break;
134 case 1: line_draw(Pnt[n1],Pnt[n2],&d); break;
135 case 4: glyph_draw(p,&d); break;
136 case 0: mark_draw(Pnt[n1],n4,p.s,&d); break;
137 case 2: trig_draw(Pnt[n1],Pnt[n2],Pnt[n3],true,&d); break;
138 }
139 }
140 else if(Quality&MGL_DRAW_FAST) for(size_t k=0;k<Prm.size();k++)
141 {
142 if(Stop) break;
143 const mglPrim &p=GetPrm(k); d.copy(p);
144 long n1=p.n1, n2=p.n2, n3=p.n3, n4=p.n4;
145 switch(p.type)
146 {
147 case 3: trig_draw(Pnt[n1],Pnt[n2],Pnt[n4],true,&d);
148 trig_draw(Pnt[n1],Pnt[n3],Pnt[n4],true,&d); break;
149 case 1: line_draw(Pnt[n1],Pnt[n2],&d); break;
150 case 4: glyph_draw(p,&d); break;
151 case 0: mark_draw(Pnt[n1],n4,p.s,&d); break;
152 case 2: trig_draw(Pnt[n1],Pnt[n2],Pnt[n3],true,&d); break;
153 }
154 }
155 else for(size_t k=0;k<Prm.size();k++)
156 {
157 if(Stop) break;
158 const mglPrim &p=GetPrm(k); d.copy(p);
159 long n1=p.n1, n2=p.n2, n3=p.n3, n4=p.n4;
160 switch(p.type)
161 {
162 case 3: fast_draw(Pnt[n1],Pnt[n4],&d); fast_draw(Pnt[n2],Pnt[n3],&d); break;
163 case 1: fast_draw(Pnt[n1],Pnt[n2],&d); break;
164 case 4: glyph_draw(p,&d); break;
165 case 0: mark_draw(Pnt[n1],n4,p.s,&d); break;
166 case 2: fast_draw(Pnt[n1],Pnt[n2],&d); fast_draw(Pnt[n1],Pnt[n3],&d);
167 fast_draw(Pnt[n2],Pnt[n3],&d); break;
168 }
169 }
170 }
171 }
172 //-----------------------------------------------------------------------------
pxl_dotsdr(long id,long n,const void *)173 void mglCanvas::pxl_dotsdr(long id, long n, const void *)
174 {
175 const float *b = Bp.b;
176 const float dx = -Bp.x*Width/2, dy = -Bp.y*Height/2, dz = Depth/2.;
177 #if !MGL_HAVE_PTHREAD
178 #pragma omp parallel for
179 #endif
180 for(long i=id;i<n;i+=mglNumThr)
181 {
182 unsigned char r[4]={0,0,0,255};
183 const mglPnt &p=Pnt[i];
184 if(p.sub<0) continue;
185 float x = p.xx-Width/2., y = p.yy-Height/2., z = p.zz-Depth/2.,xx,yy,zz;
186 xx = b[0]*x + b[1]*y + b[2]*z + dx;
187 yy = b[3]*x + b[4]*y + b[5]*z + dy;
188 zz = b[6]*x + b[7]*y + b[8]*z + dz;
189 float d = get_persp(Bp.pf,zz,Depth);
190 xx = Width/2. + d*xx; yy = Height/2. + d*yy;
191
192 r[0] = (unsigned char)(255*p.r);
193 r[1] = (unsigned char)(255*p.g);
194 r[2] = (unsigned char)(255*p.b);
195 long i0=long(xx)+Width*(Height-1-long(yy));
196 unsigned char *c = C+12*i0;
197 if(i0>=0 && i0<Width*Height && zz>Z[3*i0])
198 { Z[3*i0]=z; c[0]=r[0]; c[1]=r[1]; c[2]=r[2]; c[3]=r[3]; OI[i0]=-1; }
199 }
200 }
201 //-----------------------------------------------------------------------------
pxl_other(long id,long n,const void * p)202 void mglCanvas::pxl_other(long id, long n, const void *p)
203 {
204 const mglCanvas *gr = (const mglCanvas *)p;
205 if(Quality&MGL_DRAW_NORM)
206 #if !MGL_HAVE_PTHREAD
207 #pragma omp parallel for
208 #endif
209 for(long k=id;k<n;k+=mglNumThr)
210 {
211 long i = k%Width, j = Height-1-(k/Width);
212 pnt_plot(i,j,gr->Z[3*k+2],gr->C+12*k+8,gr->OI[k]);
213 pnt_plot(i,j,gr->Z[3*k+1],gr->C+12*k+4,gr->OI[k]);
214 pnt_plot(i,j,gr->Z[3*k],gr->C+12*k,gr->OI[k]);
215 }
216 else
217 #if !MGL_HAVE_PTHREAD
218 #pragma omp parallel for
219 #endif
220 for(long k=id;k<n;k+=mglNumThr)
221 {
222 long i = k%Width, j = Height-1-(k/Width);
223 pnt_plot(i,j,gr->Z[3*k],gr->C+12*k,gr->OI[k]);
224 }
225 }
226 //-----------------------------------------------------------------------------
pnt_plot(long x,long y,mreal z,const unsigned char ci[4],int obj_id)227 void mglCanvas::pnt_plot(long x,long y,mreal z,const unsigned char ci[4], int obj_id)
228 {
229 if(ci[3])
230 {
231 long i0=x+Width*(Height-1-y);
232 unsigned char *cc = C+12*i0, c[4];
233 c[0]=ci[0]; c[1]=ci[1]; c[2]=ci[2]; c[3]=ci[3];
234 float *zz = Z+3*i0, zf = FogDist*(z/Depth-0.5-FogDz);
235 // try to remove double transparency near vertexes
236 if(fabs(z-zz[0])<1 && OI[i0]==obj_id && abs(cc[0]-ci[0])+abs(cc[1]-ci[1])+abs(cc[2]-ci[2])<5)
237 {
238 if(cc[3]<ci[3])
239 { cc[0]=c[0]; cc[1]=c[1]; cc[2]=c[2]; cc[3]=c[3]; }
240 return;
241 }
242 if(zf<0) // add fog
243 {
244 int d = int(255.f-255.f*exp(5.f*zf));
245 unsigned char cb[4] = {BDef[0], BDef[1], BDef[2], (unsigned char)d};
246 if(d==255) return;
247 combine(c,cb);
248 }
249 if(Quality&MGL_DRAW_NORM)
250 {
251 if(z>=zz[1]) // shift point on slice down and paste new point
252 {
253 zz[2] = zz[1]; combine(cc+8,cc+4);
254 if(z>=zz[0])
255 { zz[1] = zz[0]; zz[0] = z; OI[i0]=obj_id;
256 cc[4]=cc[0]; cc[0]=c[0]; cc[5]=cc[1]; cc[1]=c[1];
257 cc[6]=cc[2]; cc[2]=c[2]; cc[7]=cc[3]; cc[3]=c[3]; }
258 else
259 { zz[1] = z; cc[4]=c[0]; cc[5]=c[1]; cc[6]=c[2]; cc[7]=c[3]; }
260 }
261 else
262 {
263 if(z>=zz[2]) // shift point on slice down and paste new point
264 { zz[2] = z; combine(cc+8,c); }
265 else // point below the background
266 { combine(c,cc+8); cc[8]=c[0]; cc[9]=c[1]; cc[10]=c[2];cc[11]=c[3]; }
267 }
268 }
269 if(Quality&MGL_DRAW_FAST)
270 {
271 if(z>=zz[0]) // point upper the background
272 { zz[0]=z; combine(cc,c); OI[i0]=obj_id; }
273 else
274 { combine(c,cc); cc[6]=cc[2]; cc[2]=c[2]; cc[7]=cc[3]; cc[3]=c[3]; }
275 }
276 else
277 {
278 if(z>=zz[0]) // point upper the background
279 { zz[0]=z; cc[0]=c[0]; cc[1]=c[1]; cc[2]=c[2]; cc[3]=c[3]; OI[i0]=obj_id; }
280 }
281 }
282 }
283 //-----------------------------------------------------------------------------
mexp(float x)284 inline float mexp(float x) // exp(-x) ~ 1/(1+x+x^2/2+x^3/4+x^5/40)
285 { return 1.f/(1.f+x*(1.f+x*x/2.f)); }
286 //{ return 1/(1+x*(1+x/2*(1+x/2*(1+x*x/10)))); }
287 //{ return exp(-x); }
col2int(const mglPnt & p,unsigned char * r,int obj_id) const288 void mglCanvas::col2int(const mglPnt &p,unsigned char *r, int obj_id) const
289 {
290 // if(!r) return r; // NOTE r must be provided!
291 if(p.a<=0) { r[0]=r[1]=r[2]=r[3]=0; return; }
292 float b0=0,b1=0,b2=0, ar,ag,ab,dif;
293 const size_t nl = p.sub>=0?p.sub:-1-p.sub;
294 const bool glob = !get(MGL_LOCAL_LIGHT);
295 ar = ag = ab = glob?AmbBr:Sub[nl].AmbBr;
296 dif = glob?DifBr:Sub[nl].DifBr;
297 const mglLight *gll = glob?light:Sub[nl].light;
298
299 if(mgl_isnum(p.u+p.v+p.w))
300 {
301 for(long i=0;i<10;i++)
302 {
303 const mglLight &ll=gll[i];
304 if(!ll.n) continue;
305 if(mgl_isnan(ll.q.x)) // source at infinity
306 {
307 const mglPoint &lp = ll.p;
308 float nn = 2*(p.u*lp.x+p.v*lp.y+p.w*lp.z) / (p.u*p.u+p.v*p.v+p.w*p.w+1e-6f);
309 float d0 = lp.x - p.u*nn;
310 float d1 = lp.y - p.v*nn;
311 float d2 = lp.z - p.w*nn;
312 nn = 1 + d2/sqrt(d0*d0+d1*d1+d2*d2+1e-6f);
313
314 // nn = exp(-ll.a*nn)*ll.b*2;
315 nn = mexp(ll.a*nn)*ll.b*2.f;
316 const mglColor &lc = ll.c;
317 b0 += nn*lc.r;
318 b1 += nn*lc.g;
319 b2 += nn*lc.b;
320 }
321 else // diffuse and specular light
322 {
323 const mglPoint &lp=ll.p, &lq=ll.q;
324 float d0 = lq.x-p.x; // direction to light source
325 float d1 = lq.y-p.y;
326 float d2 = lq.z-p.z;
327 float nn = 1+(d0*lp.x+d1*lp.y+d2*lp.z)/sqrt(d0*d0+d1*d1+d2*d2+1e-6f);
328 // float bb = exp(-3*ll.a*nn); nn = bb*dif*2;
329 float bb = mexp(3*ll.a*nn); nn = bb*dif*2;
330 const mglColor &lc = ll.c;
331 ar += nn*lc.r;
332 ag += nn*lc.g;
333 ab += nn*lc.b;
334
335 nn = 2*(p.u*d0+p.v*d1+p.w*d2) / (p.u*p.u+p.v*p.v+p.w*p.w+1e-6f);
336 d0 -= p.u*nn; d1 -= p.v*nn; d2 -= p.w*nn;
337 nn = 1 + d2/sqrt(d0*d0+d1*d1+d2*d2+1e-6f);
338
339 // nn = exp(-ll.a*nn)*bb*ll.b*2;
340 nn = mexp(ll.a*nn)*bb*ll.b*2.f;
341 b0 += nn*lc.r;
342 b1 += nn*lc.g;
343 b2 += nn*lc.b;
344 }
345 }
346 b0 += (ar>1 ? 1:ar)*p.r; // diffuse light
347 b1 += (ag>1 ? 1:ag)*p.g;
348 b2 += (ab>1 ? 1:ab)*p.b;
349 b0 = b0<1 ? b0 : 1; // normalize components
350 b1 = b1<1 ? b1 : 1;
351 b2 = b2<1 ? b2 : 1;
352 }
353 else
354 { b0=p.r; b1=p.g; b2=p.b; }
355 // try to highlight faces
356 if(obj_id==HighId) { b0*=0.7; b1*=0.7; b2*=0.7; }
357 r[0] = (unsigned char)(255*b0);
358 r[1] = (unsigned char)(255*b1);
359 r[2] = (unsigned char)(255*b2);
360 // r[3] = get(MGL_ENABLE_ALPHA) ? (unsigned char)(255*p.a) : 255;
361 r[3] = (unsigned char)((Quality&MGL_DRAW_NORM)?255*p.a:255);
362 // return r;
363 }
364 //-----------------------------------------------------------------------------
365 /// color mixing: color c1 is under color c2 !!!
combine(unsigned char * c1,const unsigned char * c2) const366 void mglCanvas::combine(unsigned char *c1, const unsigned char *c2) const
367 {
368 if(c2[3])
369 {
370 const unsigned a1=c1[3], a2=c2[3];
371 if((Flag&3)==0)
372 {
373 unsigned b1=255-a2;
374 c1[0] = (c1[0]*b1 + c2[0]*a2)/256;
375 c1[1] = (c1[1]*b1 + c2[1]*a2)/256;
376 c1[2] = (c1[2]*b1 + c2[2]*a2)/256;
377 c1[3] = (unsigned char)(a2+a1*b1/255);
378 }
379 else if((Flag&3)==1)
380 {
381 c1[0] = (unsigned char)((255-a1*(255-c1[0])/256)*(255-a2*(255-c2[0])/256)/256);
382 c1[1] = (unsigned char)((255-a1*(255-c1[1])/256)*(255-a2*(255-c2[1])/256)/256);
383 c1[2] = (unsigned char)((255-a1*(255-c1[2])/256)*(255-a2*(255-c2[2])/256)/256);
384 c1[3] = 255;
385 }
386 else if((Flag&3)==2)
387 {
388 unsigned b1,b2,b3;
389 b1 = (c1[0]*a1 + c2[0]*a2)/255; c1[0] = b1<255 ? b1 : 255;
390 b2 = (c1[1]*a1 + c2[1]*a2)/255; c1[1] = b2<255 ? b2 : 255;
391 b3 = (c1[2]*a1 + c2[2]*a2)/255; c1[2] = b3<255 ? b3 : 255;
392 c1[3] = 255;
393 }
394 }
395 }
396 //-----------------------------------------------------------------------------
397 // it looks as MSV=4 is optimal for speed-vs-quality
398 #define MSV 4
visible(long i,long j,const unsigned char m[8],mreal pw,int a)399 int inline visible(long i, long j, const unsigned char m[8], mreal pw, int a) // Check if pixel visible
400 {
401 float c = mgl_cos[(a+360)%360], s = mgl_cos[(a+450)%360];
402 int ii,jj,ss=0;
403 #if MSV==4
404 ii = int(0.5+((i-0.33)*c+(j-0.33)*s)/pw)&7; jj = int(0.5+((j-0.33)*c-(i-0.33)*s)/pw)&7;
405 if(m[jj] & (1L<<ii)) ss++;
406 ii = int(0.5+((i-0.33)*c+(j+0.33)*s)/pw)&7; jj = int(0.5+((j+0.33)*c-(i-0.33)*s)/pw)&7;
407 if(m[jj] & (1L<<ii)) ss++;
408 ii = int(0.5+((i+0.33)*c+(j-0.33)*s)/pw)&7; jj = int(0.5+((j-0.33)*c-(i+0.33)*s)/pw)&7;
409 if(m[jj] & (1L<<ii)) ss++;
410 ii = int(0.5+((i+0.33)*c+(j+0.33)*s)/pw)&7; jj = int(0.5+((j+0.33)*c-(i+0.33)*s)/pw)&7;
411 if(m[jj] & (1L<<ii)) ss++;
412 #elif MSV==9
413 ii = int(0.5+((i-0.25)*c+(j-0.25)*s)/pw)&7; jj = int(0.5+((j-0.25)*c-(i-0.25)*s)/pw)&7;
414 if(m[jj] & (1L<<ii)) ss++;
415 ii = int(0.5+((i-0.25)*c+j*s)/pw)&7; jj = int(0.5+(j*c-(i-0.25)*s)/pw)&7;
416 if(m[jj] & (1L<<ii)) ss++;
417 ii = int(0.5+((i-0.25)*c+(j+0.25)*s)/pw)&7; jj = int(0.5+((j+0.25)*c-(i-0.25)*s)/pw)&7;
418 if(m[jj] & (1L<<ii)) ss++;
419
420 ii = int(0.5+(i*c+(j-0.25)*s)/pw)&7; jj = int(0.5+((j-0.25)*c-i*s)/pw)&7;
421 if(m[jj] & (1L<<ii)) ss++;
422 ii = int(0.5+(i*c+j*s)/pw)&7; jj = int(0.5+(j*c-i*s)/pw)&7;
423 if(m[jj] & (1L<<ii)) ss++;
424 ii = int(0.5+(i*c+(j+0.25)*s)/pw)&7; jj = int(0.5+((j+0.25)*c-i*s)/pw)&7;
425 if(m[jj] & (1L<<ii)) ss++;
426
427 ii = int(0.5+((i+0.25)*c+(j-0.25)*s)/pw)&7; jj = int(0.5+((j-0.25)*c-(i+0.25)*s)/pw)&7;
428 if(m[jj] & (1L<<ii)) ss++;
429 ii = int(0.5+((i+0.25)*c+j*s)/pw)&7; jj = int(0.5+(j*c-(i+0.25)*s)/pw)&7;
430 if(m[jj] & (1L<<ii)) ss++;
431 ii = int(0.5+((i+0.25)*c+(j+0.25)*s)/pw)&7; jj = int(0.5+((j+0.25)*c-(i+0.25)*s)/pw)&7;
432 if(m[jj] & (1L<<ii)) ss++;
433 #elif MSV==1
434 ii = int(0.5+(i*c+j*s)/pw)&7; jj = int(0.5+(j*c-i*s)/pw)&7;
435 if(m[jj] & (1L<<ii)) ss++;
436 #endif
437 return ss;
438 }
439 //-----------------------------------------------------------------------------
440 /* Bilinear interpolation r(u,v) = r0 + (r1-r0)*u + (r2-r0)*v + (r3+r0-r1-r2)*u*v
441 is used (where r is one of {x,y,z,R,G,B,A}. Variables u,v are determined
442 for each point (x,y) and selected one pair which 0<u<1 and 0<v<1.*/
quad_draw(const mglPnt & p1,const mglPnt & p2,const mglPnt & p3,const mglPnt & p4,const mglDrawReg * d)443 void mglCanvas::quad_draw(const mglPnt &p1, const mglPnt &p2, const mglPnt &p3, const mglPnt &p4, const mglDrawReg *d)
444 {
445 if(Quality&MGL_DRAW_LMEM)
446 {
447 if(!(Quality&3))
448 { fast_draw(p1,p4,d); fast_draw(p2,p3,d); return; }
449 if(!(Quality&MGL_DRAW_NORM))
450 { trig_draw(p1,p2,p4,true,d); trig_draw(p1,p3,p4,true,d); return; }
451 }
452 unsigned char r[4];
453 long y1,x1,y2,x2;
454 mglPnt d1(p2-p1), d2(p3-p1), d3(p4+p1-p2-p3);
455
456 if(d1.x==0 && d1.y==0) { trig_draw(p1,p3,p4,true,d); return; }
457 if(d2.x==0 && d2.y==0) { trig_draw(p1,p2,p4,true,d); return; }
458
459 x1 = long(mgl_min(mgl_min(p1.x,p2.x), mgl_min(p3.x,p4.x))); // bounding box
460 y1 = long(mgl_min(mgl_min(p1.y,p2.y), mgl_min(p3.y,p4.y)));
461 x2 = long(mgl_max(mgl_max(p1.x,p2.x), mgl_max(p3.x,p4.x)));
462 y2 = long(mgl_max(mgl_max(p1.y,p2.y), mgl_max(p3.y,p4.y)));
463 x1=mgl_imax(x1,d->x1); x2=mgl_imin(x2,d->x2);
464 y1=mgl_imax(y1,d->y1); y2=mgl_imin(y2,d->y2);
465 if(x1>x2 || y1>y2) return;
466
467 const float dd = d1.x*d2.y-d1.y*d2.x;
468 const float dsx =-4*(d2.y*d3.x - d2.x*d3.y)*d1.y;
469 const float dsy = 4*(d2.y*d3.x - d2.x*d3.y)*d1.x;
470
471 mglPoint n1(mglPoint(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z)^mglPoint(p3.x-p1.x,p3.y-p1.y,p3.z-p1.z));
472 mglPoint n2(mglPoint(p2.x-p4.x,p2.y-p4.y,p2.z-p4.z)^mglPoint(p3.x-p4.x,p3.y-p4.y,p3.z-p4.z));
473 mglPoint nr((n1.x+n2.x)*0.5,(n1.y+n2.y)*0.5,(n1.z+n2.z)*0.5);
474
475 const float x0 = p1.x, y0 = p1.y;
476 const int oi = d->ObjId, ang=d->angle;
477 const mreal pw = d->PenWidth;
478 const uint64_t pd = d->PDef;
479
480 mglPnt tmp(p1+d1+d2+d3), pp(p1);
481 if(mgl_isnan(tmp.u) && mgl_isnum(tmp.v))
482 { pp.u = nr.x; pp.v = nr.y; pp.w = nr.z;
483 d1.u=d1.v=d1.w=d2.u=d2.v=d2.w=d3.u=d3.v=d3.w=0; }
484
485 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
486 {
487 int ms = (pd==MGL_SOLID_MASK) ? 4 : visible(i,j,d->m, pw,ang);
488 if(ms!=0)
489 {
490 float xx = (i-x0), yy = (j-y0), s;
491 s = dsx*xx + dsy*yy + (dd+d3.y*xx-d3.x*yy)*(dd+d3.y*xx-d3.x*yy);
492 if(s>=0)
493 {
494 s = sqrt(s);
495 float qu = d3.x*yy - d3.y*xx + dd + s;
496 float qv = d3.y*xx - d3.x*yy + dd + s;
497 float u = 2.f*(d2.y*xx - d2.x*yy)/qu;
498 float v = 2.f*(d1.x*yy - d1.y*xx)/qv;
499 if(u*(1.f-u)<0.f || v*(1.f-v)<0.f) // first root bad
500 {
501 qu = d3.x*yy - d3.y*xx + dd - s;
502 qv = d3.y*xx - d3.x*yy + dd - s;
503 // u = v = -1.f;
504 u = 2.f*(d2.y*xx - d2.x*yy)/qu; v = 2.f*(d1.x*yy - d1.y*xx)/qv;
505 if(u*(1.f-u)<0.f || v*(1.f-v)<0.f) continue; // second root bad
506 }
507 mglPnt p(pp+d1*u+d2*v+d3*(u*v)); col2int(p,r,oi);
508 r[3] = ms*r[3]/MSV;
509 if(r[3]) pnt_plot(i,j,p.z,r,oi);
510 }
511 }
512 }
513 }
514 //-----------------------------------------------------------------------------
515 /* Linear interpolation r(u,v) = r0 + (r1-r0)*u + (r2-r0)*v is used, where r is
516 one of {x,y,z,R,G,B,A}. Variables u,v are determined for each point (x,y).
517 Point plotted is u>0 and v>0 and u+v<1.*/
trig_draw(const mglPnt & p1,const mglPnt & p2,const mglPnt & p3,bool anorm,const mglDrawReg * d)518 void mglCanvas::trig_draw(const mglPnt &p1, const mglPnt &p2, const mglPnt &p3, bool anorm, const mglDrawReg *d)
519 {
520 if(!(Quality&3) && anorm)
521 { fast_draw(p1,p2,d); fast_draw(p1,p3,d); fast_draw(p2,p3,d); return; }
522 unsigned char r[4];
523 long y1,x1,y2,x2;
524 mglPnt d1(p2-p1), d2(p3-p1);
525
526 const float tmp = d2.x*d1.y - d1.x*d2.y;
527 if(fabs(tmp)<1e-5) return; // points lies on the same line
528 const float dyv =-d1.x/tmp, dxv = d1.y/tmp;
529 const float dyu = d2.x/tmp, dxu =-d2.y/tmp;
530
531 x1 = long(mgl_min(p1.x<p2.x?p1.x:p2.x, p3.x)); // bounding box
532 y1 = long(mgl_min(p1.y<p2.y?p1.y:p2.y, p3.y));
533 x2 = long(mgl_max(p1.x>p2.x?p1.x:p2.x, p3.x));
534 y2 = long(mgl_max(p1.y>p2.y?p1.y:p2.y, p3.y));
535 x1=x1>d->x1?x1:d->x1; x2=x2<d->x2?x2:d->x2;
536 y1=y1>d->y1?y1:d->y1; y2=y2<d->y2?y2:d->y2;
537 if(x1>x2 || y1>y2) return;
538 // default normale
539 const mglPoint nr(mglPoint(p2.x-p1.x,p2.y-p1.y,p2.z-p1.z)^mglPoint(p3.x-p1.x,p3.y-p1.y,p3.z-p1.z));
540 const float x0 = p1.x, y0 = p1.y;
541 // provide additional height to be well visible on the surfaces
542 const float dz = anorm? 0 : (Width>2 ? 1 : 1e-5*Width);
543 const int oi = d->ObjId, ang=d->angle;
544 const mreal pw = d->PenWidth;
545 const uint64_t pd = d->PDef;
546
547 mglPnt tp(p1+d1+d2), pp(p1);
548 if(mgl_isnan(tp.u) && mgl_isnum(tp.v))
549 { pp.u = nr.x; pp.v = nr.y; pp.w = nr.z;
550 d1.u=d1.v=d1.w=d2.u=d2.v=d2.w=0; }
551
552 if(Quality&MGL_DRAW_NORM) for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
553 {
554 int ms = (pd==MGL_SOLID_MASK) ? 4 : visible(i,j,d->m, pw,ang);
555 if(ms!=0)
556 {
557 float xx = (i-x0), yy = (j-y0);
558 float u = dxu*xx+dyu*yy, v = dxv*xx+dyv*yy;
559 if(u<0 || v<0 || u+v>1) continue;
560 mglPnt p(pp+d1*u+d2*v); col2int(p,r,oi);
561 r[3] = ms*r[3]/MSV;
562 if(r[3]) pnt_plot(i,j,p.z+dz,r,oi);
563 }
564 }
565 else
566 {
567 col2int(p1,r,oi);
568 float zz = p1.z+dz, dz1=d1.z, dz2=d2.z;
569 unsigned char ra = r[3];
570 if(ra) for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
571 {
572 int ms = (pd==MGL_SOLID_MASK) ? 4 : visible(i,j,d->m, pw,ang);
573 if(ms!=0)
574 {
575 float xx = (i-x0), yy = (j-y0);
576 float u = dxu*xx+dyu*yy, v = dxv*xx+dyv*yy;
577 if(u<0 || v<0 || u+v>1) continue;
578 r[3] = ms*ra/MSV;
579 pnt_plot(i,j,zz+dz1*u+dz2*v,r,oi);
580 }
581 }
582 }
583 }
584 //-----------------------------------------------------------------------------
mgl_sline(unsigned char c,float x)585 inline unsigned char mgl_sline(unsigned char c,float x)
586 { x*=x/2; return (unsigned char)(c/(1+x+x*x/5)); }
line_draw(const mglPnt & p1,const mglPnt & p2,const mglDrawReg * dr)587 void mglCanvas::line_draw(const mglPnt &p1, const mglPnt &p2, const mglDrawReg *dr)
588 {
589 if((Quality&3)==MGL_DRAW_WIRE) { fast_draw(p1,p2,dr); return; } // previously was <2. This may slightly slow down for Quality=1
590 unsigned char r[4];
591 long y1,x1,y2,x2;
592
593 const float dz = Width>2 ? 1 : 1e-5*Width; // provide additional height to be well visible on the surfaces
594 const int oi = dr->ObjId;
595 const float pw=dr->PenWidth*(oi==HighId?2:1), dpw=pen_delta*(oi==HighId?2:3);
596 const mglPnt d(p2-p1);
597 bool hor = fabs(d.x)>fabs(d.y);
598
599 x1 = long(p1.x<p2.x?p1.x:p2.x); y1 = long(p1.y<p2.y?p1.y:p2.y); // bounding box
600 x2 = long(p1.x>p2.x?p1.x:p2.x); y2 = long(p1.y>p2.y?p1.y:p2.y);
601 x1 -= pw+10/dpw; x2 += pw+10/dpw;
602 y1 -= pw+10/dpw; y2 += pw+10/dpw;
603 x1=x1>dr->x1?x1:dr->x1; x2=x2<dr->x2?x2:dr->x2;
604 y1=y1>dr->y1?y1:dr->y1; y2=y2<dr->y2?y2:dr->y2;
605 const float dd = hypot(d.x, d.y);
606 if(x1>x2 || y1>y2 || dd<1e-5) return;
607
608 const float dxv = d.y/dd, dyv =-d.x/dd;
609 const float dxu = d.x/dd, dyu = d.y/dd;
610
611 const uint64_t pd = dr->PDef;
612 const mreal pp = dr->pPos, V = (pw-1)*(pw-1)/4, S = (1-pw)/2;
613 if(hor) for(long i=x1;i<=x2;i++)
614 {
615 y1 = int(p1.y+d.y*(i-p1.x)/d.x - pw - 10/dpw);
616 y2 = int(p1.y+d.y*(i-p1.x)/d.x + pw + 10/dpw);
617 y1=y1>dr->y1?y1:dr->y1; y2=y2<dr->y2?y2:dr->y2;
618 for(long j=y1;j<=y2;j++)
619 {
620 float xx = (i-p1.x), yy = (j-p1.y);
621 float u = dxu*xx+dyu*yy, v = dxv*xx+dyv*yy; v = v*v;
622 if(u<0) { v += 16*u*u; u=0; }
623 else if(u>dd) { v += 16*(u-dd)*(u-dd); u=dd; }
624 if( pd & ((uint64_t)1<<(long(pp+u/pw)%16)) )
625 {
626 mglPnt p(p1+d*(u/dd)); col2int(p,r,oi);
627 r[3] = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
628 if(r[3]) pnt_plot(i,j,p.z+dz,r,oi);
629 }
630 }
631 }
632 else for(long j=y1;j<=y2;j++)
633 {
634 x1 = int(p1.x+d.x*(j-p1.y)/d.y - pw - 10/dpw);
635 x2 = int(p1.x+d.x*(j-p1.y)/d.y + pw + 10/dpw);
636 x1=x1>dr->x1?x1:dr->x1; x2=x2<dr->x2?x2:dr->x2;
637 for(long i=x1;i<=x2;i++)
638 {
639 float xx = (i-p1.x), yy = (j-p1.y);
640 float u = dxu*xx+dyu*yy, v = dxv*xx+dyv*yy; v = v*v;
641 if(u<0) { v += 16*u*u; u=0; }
642 else if(u>dd) { v += 16*(u-dd)*(u-dd); u=dd; }
643 if( pd & ((uint64_t)1<<(long(pp+u/pw)%16)) )
644 {
645 mglPnt p(p1+d*(u/dd)); col2int(p,r,oi);
646 r[3] = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
647 if(r[3]) pnt_plot(i,j,p.z+dz,r,oi);
648 }
649 }
650 }
651 }
652 //-----------------------------------------------------------------------------
pnt_fast(long x,long y,mreal z,const unsigned char ci[4],int obj_id)653 void mglCanvas::pnt_fast(long x,long y,mreal z,const unsigned char ci[4], int obj_id)
654 {
655 long i0=x+Width*(Height-1-y);
656 unsigned char *cc = C+12*i0;
657 if(ci[3]!=0 && z>Z[3*i0]) // point upper the background
658 { Z[3*i0]=z; OI[i0]=obj_id;
659 cc[0]=ci[0]; cc[1]=ci[1]; cc[2]=ci[2]; cc[3]=ci[3]; }
660 }
661 //-----------------------------------------------------------------------------
fast_draw(const mglPnt & p1,const mglPnt & p2,const mglDrawReg * dr)662 void mglCanvas::fast_draw(const mglPnt &p1, const mglPnt &p2, const mglDrawReg *dr)
663 {
664 if(p1.x==p2.x && p1.y==p2.y) return;
665 const mglPnt d(p2-p1);
666 const int oi = dr->ObjId;
667 unsigned char r[4]; col2int(p1,r,oi);
668 long y1,x1,y2,x2;
669
670 const bool hor = fabs(d.x)>fabs(d.y);
671
672 x1 = long(p1.x<p2.x?p1.x:p2.x); y1 = long(p1.y<p2.y?p1.y:p2.y); // bounding box
673 x2 = long(p1.x>p2.x?p1.x:p2.x); y2 = long(p1.y>p2.y?p1.y:p2.y);
674 x1=x1>dr->x1?x1:dr->x1; x2=x2<dr->x2?x2:dr->x2;
675 y1=y1>dr->y1?y1:dr->y1; y2=y2<dr->y2?y2:dr->y2;
676 if(x1>x2 || y1>y2) return;
677 const float dz = Width>2 ? 1 : 1e-5*Width; // provide additional height to be well visible on the surfaces
678
679 if(hor) for(long i=x1;i<=x2;i++)
680 {
681 long c = long(p1.y+d.y*(i-p1.x)/d.x);
682 if(c>=y1 && c<=y2)
683 pnt_fast(i, c, p1.z+d.z*(i-p1.x)/d.x+dz, r,oi);
684 }
685 else for(long i=y1;i<=y2;i++)
686 {
687 long c = long(p1.x+d.x*(i-p1.y)/d.y);
688 if(c>=x1 && c<=x2)
689 pnt_fast(c, i, p1.z+d.z*(i-p1.y)/d.y+dz, r,oi);
690 }
691 }
692 //-----------------------------------------------------------------------------
pnt_draw(const mglPnt & p,const mglDrawReg * dr)693 void mglCanvas::pnt_draw(const mglPnt &p, const mglDrawReg *dr)
694 {
695 // if(k<0 || !dr) return;
696 const int oi = dr->ObjId;
697 const float pw=(oi==HighId?6:3)*dr->PenWidth,dpw=(oi==HighId?2:3)*pen_delta;
698 unsigned char cs[4], cc;
699 col2int(p,cs,oi); cc = cs[3];
700 if(cc==0) return;
701 const long s = long(pw+10/dpw+fabs(pw));
702 const long i1=mgl_max(-s,dr->x1-p.x),i2=mgl_min(s,dr->x2-p.x);
703 const long j1=mgl_max(-s,dr->y1-p.y),j2=mgl_min(s,dr->y2-p.y);
704 const mreal V = (pw-1)*(pw-1)/4, S = (1-pw)/2;
705 if(!(Quality&3)) for(long j=j1;j<=j2;j++) for(long i=i1;i<=i2;i++) // fast draw
706 {
707 float v = i*i+j*j;
708 if(v>1+V) continue;
709 if(cs[3]) pnt_plot(p.x+i,p.y+j,p.z,cs,oi);
710 }
711 else for(long j=j1;j<=j2;j++) for(long i=i1;i<=i2;i++)
712 {
713 float v = i*i+j*j;
714 cs[3] = v<V ? cc : mgl_sline(cc,dpw*(sqrt(v)+S));
715 if(cs[3]) pnt_plot(p.x+i,p.y+j,p.z,cs,oi);
716 }
717 }
718 //-----------------------------------------------------------------------------
mark_draw(const mglPnt & q,char type,mreal size,mglDrawReg * dr)719 void mglCanvas::mark_draw(const mglPnt &q, char type, mreal size, mglDrawReg *dr)
720 {
721 const int oi = dr->ObjId;
722 unsigned char cs[4]; col2int(q,cs,oi);
723 const unsigned char ca = cs[3];// = size>0 ? 255 : 255*q.t;
724 const mreal ss=(strchr("xsSoO",type)?1:1.1)*fabs(size), dpw=(oi==HighId?2:3)*pen_delta;
725 mreal PW=1;
726
727 if(type=='.' || ss==0)
728 {
729 PW = 3*(ss?ss:sqrt(font_factor/400));
730 if(oi==HighId) PW *= 2;
731 const mreal pw = PW;
732 mreal dd = pw+10/dpw;
733 long x1 = long(q.x-dd), y1 = long(q.y-dd); // bounding box
734 long x2 = long(q.x+dd), y2 = long(q.y+dd);
735 x1=x1>dr->x1?x1:dr->x1; x2=x2<dr->x2?x2:dr->x2;
736 y1=y1>dr->y1?y1:dr->y1; y2=y2<dr->y2?y2:dr->y2;
737 if(x1>x2 || y1>y2) return;
738 const float V=(pw-1)*(pw-1)/4,S=(1-pw)/2;
739
740 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
741 {
742 float dx=i-q.x, dy=j-q.y, v=dx*dx+dy*dy;
743 int sum = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
744 cs[3] = ca*sum/255;
745 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
746 }
747 }
748 else
749 {
750 dr->PDef = MGL_SOLID_MASK; dr->angle = 0;
751 PW = dr->PenWidth*sqrt(fabs(50*size));
752 if(PW<1) PW=1;
753 if(oi==HighId) PW *= 2;
754 const mreal pw = PW;
755
756 mreal dd = ss+pw+10/dpw;
757 long x1 = long(q.x-dd), y1 = long(q.y-dd); // bounding box
758 long x2 = long(q.x+dd), y2 = long(q.y+dd);
759 x1=x1>dr->x1?x1:dr->x1; x2=x2<dr->x2?x2:dr->x2;
760 y1=y1>dr->y1?y1:dr->y1; y2=y2<dr->y2?y2:dr->y2;
761 if(x1>x2 || y1>y2) return;
762 const float V=(pw-1)*(pw-1)/4,S=(1-pw)/2;
763
764 switch(type)
765 {
766 case 'P':
767 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
768 {
769 float dx=i-q.x, dy=j-q.y, v,u;
770 int sum=0;
771 u = fabs(dy)-ss; v = (dx-ss)*(dx-ss)+(u<0?0:u*u);
772 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
773 u = fabs(dy)-ss; v = (dx+ss)*(dx+ss)+(u<0?0:u*u);
774 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
775 u = fabs(dy)-ss; v = dx*dx+(u<0?0:u*u);
776 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
777 u = fabs(dx)-ss; v = (dy-ss)*(dy-ss)+(u<0?0:u*u);
778 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
779 u = fabs(dx)-ss; v = (dy+ss)*(dy+ss)+(u<0?0:u*u);
780 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
781 u = fabs(dx)-ss; v = dy*dy+(u<0?0:u*u);
782 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
783 sum = sum>255?255:sum; cs[3] = ca*sum/255;
784 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
785 }
786 break;
787 case '+':
788 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
789 {
790 float dx=i-q.x, dy=j-q.y, v,u;
791 int sum=0;
792 u = fabs(dy)-ss; v = dx*dx+(u<0?0:u*u);
793 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
794 u = fabs(dx)-ss; v = dy*dy+(u<0?0:u*u);
795 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
796 sum = sum>255?255:sum; cs[3] = ca*sum/255;
797 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
798 }
799 break;
800 case 'X':
801 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
802 {
803 float dx=i-q.x, dy=j-q.y, v,u;
804 int sum=0;
805 u = fabs(dy)-ss; v = (dx-ss)*(dx-ss)+(u<0?0:u*u);
806 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
807 u = fabs(dy)-ss; v = (dx+ss)*(dx+ss)+(u<0?0:u*u);
808 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
809 u = fabs(dx)-ss; v = (dy-ss)*(dy-ss)+(u<0?0:u*u);
810 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
811 u = fabs(dx)-ss; v = (dy+ss)*(dy+ss)+(u<0?0:u*u);
812 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
813
814 u = fabs(dx+dy)-2*ss; v = dx-dy; v = v*v+(u<0?0:u*u);
815 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
816 u = fabs(dx-dy)-2*ss; v = dx+dy; v = v*v+(u<0?0:u*u);
817 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
818
819 sum = sum>255?255:sum; cs[3] = ca*sum/255;
820 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
821 }
822 break;
823 case 'x':
824 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
825 {
826 float dx=i-q.x, dy=j-q.y, v,u;
827 int sum=0;
828 u = fabs(dx+dy)-2*ss; v = dx-dy; v = v*v+(u<0?0:u*u);
829 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
830 u = fabs(dx-dy)-2*ss; v = dx+dy; v = v*v+(u<0?0:u*u);
831 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
832 sum = sum>255?255:sum; cs[3] = ca*sum/255;
833 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
834 }
835 break;
836 case 'S':
837 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
838 {
839 float dx=i-q.x, dy=j-q.y, v,u;
840 u = fabs(dy)-ss; if(u<0) u=0;
841 v = fabs(dx)-ss; if(v<0) v=0; v = u*u+v*v;
842 int sum = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
843 cs[3] = ca*sum/255;
844 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
845 }
846 break;
847 case 's':
848 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
849 {
850 float dx=i-q.x, dy=j-q.y, v,u;
851 int sum=0;
852 u = fabs(dy)-ss; v = (dx-ss)*(dx-ss)+(u<0?0:u*u);
853 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
854 u = fabs(dy)-ss; v = (dx+ss)*(dx+ss)+(u<0?0:u*u);
855 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
856 u = fabs(dx)-ss; v = (dy-ss)*(dy-ss)+(u<0?0:u*u);
857 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
858 u = fabs(dx)-ss; v = (dy+ss)*(dy+ss)+(u<0?0:u*u);
859 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
860 sum = sum>255?255:sum; cs[3] = ca*sum/255;
861 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
862 }
863 break;
864 case 'D':
865 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
866 {
867 float dx=i-q.x, dy=j-q.y, v,u;
868 u = fabs(dx-dy)-ss; if(u<0) u=0;
869 v = fabs(dx+dy)-ss; if(v<0) v=0; v = u*u+v*v;
870 int sum = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
871 cs[3] = ca*sum/255;
872 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
873 }
874 break;
875 case 'd':
876 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
877 {
878 float dx=i-q.x, dy=j-q.y, v,u;
879 int sum=0;
880 u = fabs(dx+dy)-ss; v = (dx-dy-ss)*(dx-dy-ss)+(u<0?0:u*u);
881 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
882 u = fabs(dx+dy)-ss; v = (dx-dy+ss)*(dx-dy+ss)+(u<0?0:u*u);
883 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
884 u = fabs(dx-dy)-ss; v = (dx+dy-ss)*(dx+dy-ss)+(u<0?0:u*u);
885 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
886 u = fabs(dx-dy)-ss; v = (dx+dy+ss)*(dx+dy+ss)+(u<0?0:u*u);
887 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
888 sum = sum>255?255:sum; cs[3] = ca*sum/255;
889 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
890 }
891 break;
892 case 'Y':
893 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
894 {
895 float dx=i-q.x, dy=j-q.y, v,u;
896 int sum=0;
897 u = fabs(dy+ss/2)-ss/2; v = dx*dx+(u<0?0:u*u);
898 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
899 u = fabs(0.87*dx+0.5*dy-ss/2)-ss/2; v = (0.5*dx-0.87*dy)*(0.5*dx-0.87*dy)+(u<0?0:u*u);
900 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
901 u = fabs(-0.87*dx+0.5*dy-ss/2)-ss/2; v = (0.5*dx+0.87*dy)*(0.5*dx+0.87*dy)+(u<0?0:u*u);
902 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
903 sum = sum>255?255:sum; cs[3] = ca*sum/255;
904 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
905 }
906 break;
907 case '*':
908 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
909 {
910 float dx=i-q.x, dy=j-q.y, v,u;
911 int sum=0;
912 u = fabs(dy)-ss; v = dx*dx+(u<0?0:u*u);
913 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
914 u = fabs(0.87*dx+0.5*dy)-ss; v = (0.5*dx-0.87*dy)*(0.5*dx-0.87*dy)+(u<0?0:u*u);
915 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
916 u = fabs(-0.87*dx+0.5*dy)-ss; v = (0.5*dx+0.87*dy)*(0.5*dx+0.87*dy)+(u<0?0:u*u);
917 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
918 sum = sum>255?255:sum; cs[3] = ca*sum/255;
919 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
920 }
921 break;
922 case 'T':
923 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
924 {
925 float dx=i-q.x, dy=j-q.y, v,u;
926 u=dy/1.5+ss/3; v=(dx+ss-u)/2;
927 if(u>0 && v>0 && u+v<ss) cs[3]=ca;
928 else
929 {
930 int sum=0;
931 u = fabs(dx)-ss; v = dy+ss/2; v = v*v+(u<0?0:u*u);
932 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
933 u = fabs(0.55*dx+0.83*dy)-0.9*ss; v = 0.83*dx-0.55*dy+0.55*ss; v = v*v+(u<0?0:u*u);
934 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
935 u = fabs(0.55*dx-0.83*dy)-0.9*ss; v = 0.83*dx+0.55*dy-0.55*ss; v = v*v+(u<0?0:u*u);
936 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
937 sum = sum>255?255:sum; cs[3] = ca*sum/255;
938 }
939 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
940 }
941 break;
942 case '^':
943 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
944 {
945 float dx=i-q.x, dy=j-q.y, v,u;
946 int sum=0;
947 u = fabs(dx)-ss; v = dy+ss/2; v = v*v+(u<0?0:u*u);
948 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
949 u = fabs(0.55*dx+0.83*dy)-0.9*ss; v = 0.83*dx-0.55*dy+0.55*ss; v = v*v+(u<0?0:u*u);
950 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
951 u = fabs(0.55*dx-0.83*dy)-0.9*ss; v = 0.83*dx+0.55*dy-0.55*ss; v = v*v+(u<0?0:u*u);
952 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
953 sum = sum>255?255:sum; cs[3] = ca*sum/255;
954 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
955 }
956 break;
957 case 'V':
958 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
959 {
960 float dx=i-q.x, dy=j-q.y, v,u;
961 u=-dy/1.5+ss/3; v=(dx+ss-u)/2;
962 if(u>0 && v>0 && u+v<ss) cs[3]=ca;
963 else
964 {
965 int sum=0;
966 u = fabs(dx)-ss; v = dy-ss/2; v = v*v+(u<0?0:u*u);
967 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
968 u = fabs(0.55*dx+0.83*dy)-0.9*ss; v = 0.83*dx-0.55*dy-0.55*ss; v = v*v+(u<0?0:u*u);
969 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
970 u = fabs(0.55*dx-0.83*dy)-0.9*ss; v = 0.83*dx+0.55*dy+0.55*ss; v = v*v+(u<0?0:u*u);
971 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
972 sum = sum>255?255:sum; cs[3] = ca*sum/255;
973 }
974 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
975 }
976 break;
977 case 'v':
978 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
979 {
980 float dx=i-q.x, dy=j-q.y, v,u;
981 int sum=0;
982 u = fabs(dx)-ss; v = dy-ss/2; v = v*v+(u<0?0:u*u);
983 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
984 u = fabs(0.55*dx+0.83*dy)-0.9*ss; v = 0.83*dx-0.55*dy-0.55*ss; v = v*v+(u<0?0:u*u);
985 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
986 u = fabs(0.55*dx-0.83*dy)-0.9*ss; v = 0.83*dx+0.55*dy+0.55*ss; v = v*v+(u<0?0:u*u);
987 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
988 sum = sum>255?255:sum; cs[3] = ca*sum/255;
989 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
990 }
991 break;
992 case 'L':
993 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
994 {
995 float dx=i-q.x, dy=j-q.y, v,u;
996 u=-dx/1.5+ss/3; v=(dy+ss-u)/2;
997 if(u>0 && v>0 && u+v<ss) cs[3]=ca;
998 else
999 {
1000 int sum=0;
1001 u = fabs(dy)-ss; v = dx-ss/2; v = v*v+(u<0?0:u*u);
1002 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1003 u = fabs(0.55*dy+0.83*dx)-0.9*ss; v = 0.83*dy-0.55*dx-0.55*ss; v = v*v+(u<0?0:u*u);
1004 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1005 u = fabs(0.55*dy-0.83*dx)-0.9*ss; v = 0.83*dy+0.55*dx+0.55*ss; v = v*v+(u<0?0:u*u);
1006 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1007 sum = sum>255?255:sum; cs[3] = ca*sum/255;
1008 }
1009 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1010 }
1011 break;
1012 case '<':
1013 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1014 {
1015 float dx=i-q.x, dy=j-q.y, v,u;
1016 int sum=0;
1017 u = fabs(dy)-ss; v = dx-ss/2; v = v*v+(u<0?0:u*u);
1018 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1019 u = fabs(0.55*dy+0.83*dx)-0.9*ss; v = 0.83*dy-0.55*dx-0.55*ss; v = v*v+(u<0?0:u*u);
1020 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1021 u = fabs(0.55*dy-0.83*dx)-0.9*ss; v = 0.83*dy+0.55*dx+0.55*ss; v = v*v+(u<0?0:u*u);
1022 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1023 sum = sum>255?255:sum; cs[3] = ca*sum/255;
1024 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1025 }
1026 break;
1027 case 'R':
1028 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1029 {
1030 float dx=i-q.x, dy=j-q.y, v,u;
1031 u=dx/1.5+ss/3; v=(dy+ss-u)/2;
1032 if(u>0 && v>0 && u+v<ss) cs[3]=ca;
1033 else
1034 {
1035 int sum=0;
1036 u = fabs(dy)-ss; v = dx+ss/2; v = v*v+(u<0?0:u*u);
1037 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1038 u = fabs(0.55*dy+0.83*dx)-0.9*ss; v = 0.83*dy-0.55*dx+0.55*ss; v = v*v+(u<0?0:u*u);
1039 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1040 u = fabs(0.55*dy-0.83*dx)-0.9*ss; v = 0.83*dy+0.55*dx-0.55*ss; v = v*v+(u<0?0:u*u);
1041 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1042 sum = sum>255?255:sum; cs[3] = ca*sum/255;
1043 }
1044 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1045 }
1046 break;
1047 case '>':
1048 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1049 {
1050 float dx=i-q.x, dy=j-q.y, v,u;
1051 int sum=0;
1052 u = fabs(dy)-ss; v = dx+ss/2; v = v*v+(u<0?0:u*u);
1053 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1054 u = fabs(0.55*dy+0.83*dx)-0.9*ss; v = 0.83*dy-0.55*dx+0.55*ss; v = v*v+(u<0?0:u*u);
1055 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1056 u = fabs(0.55*dy-0.83*dx)-0.9*ss; v = 0.83*dy+0.55*dx-0.55*ss; v = v*v+(u<0?0:u*u);
1057 sum += v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1058 sum = sum>255?255:sum; cs[3] = ca*sum/255;
1059 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1060 }
1061 break;
1062 case 'O':
1063 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1064 {
1065 float dx=i-q.x, dy=j-q.y, v;
1066 v = hypot(dx,dy)-ss; v=v<0?0:v*v;
1067 int sum = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1068 cs[3] = ca*sum/255;
1069 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1070 }
1071 break;
1072 case 'o':
1073 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1074 {
1075 float dx=i-q.x, dy=j-q.y, v;
1076 v = hypot(dx,dy)-ss; v=v*v;
1077 int sum = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1078 cs[3] = ca*sum/255;
1079 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1080 }
1081 break;
1082 case 'C':
1083 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1084 {
1085 float dx=i-q.x, dy=j-q.y, v;
1086 v = hypot(dx,dy)-ss; v=v*v;
1087 int sum = v<V ? 255 : mgl_sline(255,dpw*(sqrt(v)+S));
1088 v = dx*dx+dy*dy;
1089 sum += v<(2*pw-1)*(2*pw-1)/4 ? 255 : mgl_sline(255,dpw*(sqrt(v)+(1-2*pw)/2));
1090 sum = sum>255?255:sum; cs[3] = ca*sum/255;
1091 if(cs[3]) pnt_plot(i,j,q.z+1,cs,oi);
1092 }
1093 break;
1094 }
1095 }
1096 }
1097 //-----------------------------------------------------------------------------
1098 // scale direction for new view/zoom
GetGlyphPhi(const mglPnt & q,float phi)1099 float mglCanvas::GetGlyphPhi(const mglPnt &q, float phi)
1100 {
1101 float x,y,z,ll;
1102 if(q.sub<0)
1103 { x = q.u; y = q.v; z = q.w; }
1104 else
1105 {
1106 x = Bp.b[0]*q.u + Bp.b[1]*q.v + Bp.b[2]*q.w;
1107 y = Bp.b[3]*q.u + Bp.b[4]*q.v + Bp.b[5]*q.w;
1108 z = Bp.b[6]*q.u + Bp.b[7]*q.v + Bp.b[8]*q.w;
1109
1110 float dv= get_persp(Bp.pf,q.z,Depth);
1111 float c = get_pfact(Bp.pf,Depth);
1112 x += (q.x-Width/2)*z*c*dv;
1113 y += (q.y-Height/2)*z*c*dv;
1114 }
1115 ll = x*x+y*y;
1116 if(ll < 1e-10) return NAN;
1117 if(ll==ll && phi<1e4)
1118 {
1119 phi = -atan2(y,x)*180/M_PI;
1120 // if(fabs(phi)>90) phi+=180; // NOTE this is 2nd part of rotation changes (see also text_plot())
1121 }
1122 else phi=0;
1123 return phi;
1124 }
1125 //-----------------------------------------------------------------------------
glyph_draw(const mglPrim & P,mglDrawReg * d)1126 void mglCanvas::glyph_draw(const mglPrim &P, mglDrawReg *d)
1127 {
1128 float phi = GetGlyphPhi(Pnt[P.n2],P.w);
1129 if(mgl_isnan(phi)) return;
1130 if(d) { d->PDef = MGL_SOLID_MASK; d->angle = 0; d->PenWidth=(P.n3&4)?1.2:0.8; }
1131
1132 mglPnt p=Pnt[P.n1]; p.a=1;
1133 mreal fact = get_persp(Bp.pf,p.z,Depth);
1134 mreal pf=p.sub<0?1:sqrt((Bp.b[0]*Bp.b[0]+Bp.b[1]*Bp.b[1]+Bp.b[3]*Bp.b[3]+Bp.b[4]*Bp.b[4])/2)*fact;
1135 mreal size=P.s, f = P.p*pf*P.s;
1136 p.u *= pf*size; p.v *= pf*size;
1137
1138 const mglGlyph &g = Glf[P.n4];
1139 if(P.n3&8)
1140 {
1141 if(!(P.n3&4)) glyph_line(phi,p,f,true, d);
1142 glyph_line(phi,p,f,false, d);
1143 }
1144 else
1145 {
1146 if(!(P.n3&4)) glyph_fill(phi,p,f,g, d);
1147 glyph_wire(phi,p,f,g, d);
1148 }
1149 }
1150 //-----------------------------------------------------------------------------
mgl_addpnts(mreal x1,mreal y1,mreal x2,mreal y2,std::vector<mreal> * b)1151 void static mgl_addpnts(mreal x1,mreal y1,mreal x2,mreal y2, std::vector<mreal> *b)
1152 {
1153 // if(x1>x2) { mreal t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t; }
1154 if(y1<y2) for(int i=long(y1);i<=long(y2)+1;i++)
1155 {
1156 mreal d = (i-y1)/(y2-y1);
1157 if(d>=0 && d<=1) b[i].push_back(x1+d*(x2-x1));
1158 }
1159 else for(int i=long(y2);i<=long(y1)+1;i++)
1160 {
1161 // mreal xx1 = x1+(x2-x1)*(i-y1)/(y2-y1);
1162 // mreal xx2 = x1+(x2-x1)*(i+1-y1)/(y2-y1);
1163 // if(xx1>xx2) { mreal t=xx1; xx1=xx2; xx2=t; }
1164 // if(xx1<x1) xx1=x1;
1165 // if(xx2>x2) xx2=x2;
1166 // if(i>y1 && i<y2)
1167 // {
1168 // b[i].push_back(xx1);
1169 // if(xx2>=xx1+1) { b[i].push_back(xx2); b[i].push_back(xx2); }
1170 // }
1171 mreal d = (i-y1)/(y2-y1);
1172 if(d>=0 && d<=1) b[i].push_back(x1+d*(x2-x1));
1173 }
1174 }
glyph_fill(mreal phi,const mglPnt & pp,mreal f,const mglGlyph & g,const mglDrawReg * d)1175 void mglCanvas::glyph_fill(mreal phi, const mglPnt &pp, mreal f, const mglGlyph &g, const mglDrawReg *d)
1176 {
1177 if(g.trig && g.nt>0) // slow but look very nice :(
1178 {
1179 const mreal co=cos(phi*M_PI/180), si=sin(phi*M_PI/180);
1180 mglPnt q0=pp, q1=pp, q2=pp;
1181 q0.u=q0.v=q1.u=q1.v=q2.u=q2.v=NAN;
1182 for(long ik=0;ik<g.nt;ik++)
1183 {
1184 long ii = 6*ik; mreal x, y;
1185 x = pp.u+g.trig[ii]*f; y = pp.v+g.trig[ii+1]*f;
1186 q0.x = pp.x+(x*co+y*si)/2; q0.y = pp.y+(y*co-x*si)/2; ii+=2;
1187 x = pp.u+g.trig[ii]*f; y = pp.v+g.trig[ii+1]*f;
1188 q1.x = pp.x+(x*co+y*si)/2; q1.y = pp.y+(y*co-x*si)/2; ii+=2;
1189 x = pp.u+g.trig[ii]*f; y = pp.v+g.trig[ii+1]*f;
1190 q2.x = pp.x+(x*co+y*si)/2; q2.y = pp.y+(y*co-x*si)/2;
1191 trig_draw(q0,q1,q2,false,d);
1192 }
1193 return;
1194 }
1195 if(!g.line || g.nl<=0) return;
1196 const mreal co=cos(phi*M_PI/180), si=sin(phi*M_PI/180);
1197
1198 mreal x1 = 1e10, x2=-1e10, y1=1e10, y2=-1e10;
1199 for(long i=0;i<g.nl;i++) // find sizes of glyph
1200 {
1201 long ii=2*i;
1202 if(g.line[ii]==0x3fff && g.line[ii+1]==0x3fff) continue;
1203 mreal x = pp.u + g.line[ii]*f, y = pp.v + g.line[ii+1]*f;
1204 mreal xx = pp.x+(x*co+y*si)/2, yy = pp.y+(y*co-x*si)/2;
1205 if(xx<x1) x1=xx;
1206 if(xx>x2) x2=xx;
1207 if(yy<y1) y1=yy;
1208 if(yy>y2) y2=yy;
1209 }
1210 x1-=2; x2+=2; y1-=2; y2+=2;
1211 long w = long(x2-x1+1), h = long(y2-y1+1), il=0;
1212 long x0=long(x1), y0=long(y1),i1=1,i2=w-2,j1=1,j2=h-2;
1213 if(d) // apply mglDrawReg
1214 {
1215 if(x0+i1<d->x1) i1 = d->x1-x0;
1216 if(x0+i2>d->x2) i2 = d->x2-x0;
1217 if(y0+j1<d->y1) j1 = d->y1-y0;
1218 if(y0+j2>d->y2) j2 = d->y2-y0;
1219 }
1220 else
1221 {
1222 if(x0+i1<0) i1 = -x0;
1223 if(x0+i2>Width) i2 = Width-x0;
1224 if(y0+j1<0) j1 = -y0;
1225 if(y0+j2>Height)j2 = Height-y0;
1226 }
1227 if(i1>=i2 || j1>=j2) return;
1228
1229 std::vector<mreal> *b = new std::vector<mreal>[h];
1230 const float dz = Width>2 ? 1 : 1e-5*Width; // provide additional height to be well visible on the surfaces
1231 const int oi = d?d->ObjId:-1;
1232 unsigned char r[4]; col2int(pp,r,oi);
1233 if(r[3]) for(long i=0;i<g.nl;i++) // add bounding points
1234 {
1235 long ii=2*i;
1236 mreal x = pp.u + g.line[ii]*f, y = pp.v + g.line[ii+1]*f;
1237 mreal xx1 = pp.x+(x*co+y*si)/2-x1, yy1 = pp.y+(y*co-x*si)/2-y1, xx2, yy2;
1238 if(g.line[ii]==0x3fff && g.line[ii+1]==0x3fff) // line breakthrough
1239 { il = i+1; continue; }
1240 else if(i==g.nl-1 || (g.line[ii+2]==0x3fff && g.line[ii+3]==0x3fff)) // enclose the circle
1241 {
1242 ii=2*il; x = pp.u + g.line[ii]*f; y = pp.v + g.line[ii+1]*f;
1243 xx2 = pp.x+(x*co+y*si)/2-x1; yy2 = pp.y+(y*co-x*si)/2-y1;
1244 }
1245 else // ordinary line
1246 {
1247 ii+=2; x = pp.u + g.line[ii]*f; y = pp.v + g.line[ii+1]*f;
1248 xx2 = pp.x+(x*co+y*si)/2-x1; yy2 = pp.y+(y*co-x*si)/2-y1;
1249 }
1250 mgl_addpnts(xx1,yy1,xx2,yy2,b);
1251 // draw boundary lines in any case ???
1252 if(fabs(xx2-xx1)>fabs(yy2-yy1)) // horizontal line
1253 {
1254 mreal d = (yy2-yy1)/(xx2-xx1), a = yy1-d*xx1+0.5;
1255 if(xx1>xx2) { mreal t=xx1; xx1=xx2; xx2=t; }
1256 for(long k=xx1;k<=xx2;k++)
1257 {
1258 long ii = long(k), jj = long(a+d*k);
1259 if(ii>=i1 && ii<=i2 && jj>=j1 && jj<=j2) pnt_plot(x0+ii,y0+jj,pp.z+dz,r,oi);
1260 }
1261 }
1262 else // vertical line
1263 {
1264 mreal d = (xx2-xx1)/(yy2-yy1), a = xx1-d*yy1+0.5;
1265 if(yy1>yy2) { mreal t=yy1; yy1=yy2; yy2=t; }
1266 for(long k=yy1;k<=yy2;k++)
1267 {
1268 long jj = long(k), ii = long(a+d*k);
1269 if(ii>=i1 && ii<=i2 && jj>=j1 && jj<=j2) pnt_plot(x0+ii,y0+jj,pp.z+dz,r,oi);
1270 }
1271 }
1272 }
1273 // TODO add smoothing -- if 3 neighbors >0 => set 1; if 3 neighbors=0 => set 0 ???
1274 for(long j=j1;j<=j2;j++) // draw glyph
1275 {
1276 if(b[j].size()<2) continue;
1277 std::sort(b[j].begin(),b[j].end());
1278 for(size_t k=0;k<b[j].size();k+=2)
1279 {
1280 long ii1 = long(b[j][k]+0.5), ii2=long(b[j][k+1]+0.5);
1281 // if(ii1==ii2 && b[j].size()%2==1) { k++; ii2=long(b[j][k+1]+0.5); }
1282 if(ii1<i1) ii1=i1;
1283 if(ii2>i2) ii2=i2;
1284 for(long i=ii1;i<=ii2;i++) pnt_plot(x0+i,y0+j,pp.z+dz,r,oi);
1285 }
1286 }
1287 delete []b;
1288 }
1289 //-----------------------------------------------------------------------------
glyph_wire(mreal phi,const mglPnt & pp,mreal f,const mglGlyph & g,const mglDrawReg * d)1290 void mglCanvas::glyph_wire(mreal phi, const mglPnt &pp, mreal f, const mglGlyph &g, const mglDrawReg *d)
1291 {
1292 if(!g.line || g.nl<=0) return;
1293 long il=0;
1294 const mreal co=cos(phi*M_PI/180), si=sin(phi*M_PI/180);
1295 mglPnt q0=pp, q1=pp; q0.u=q0.v=q1.u=q1.v=NAN;
1296 mglPoint p1,p2;
1297 for(long ik=0;ik<g.nl;ik++)
1298 {
1299 long ii = 2*ik;
1300 if(g.line[ii]==0x3fff && g.line[ii+1]==0x3fff) // line breakthrough
1301 { il = ik+1; continue; }
1302 else if(ik==g.nl-1 || (g.line[ii+2]==0x3fff && g.line[ii+3]==0x3fff))
1303 { // enclose the circle
1304 mreal x,y;
1305 x = pp.u+g.line[ii]*f; y = pp.v+g.line[ii+1]*f;
1306 q0.x = pp.x+(x*co+y*si)/2; q0.y = pp.y+(y*co-x*si)/2; ii=2*il;
1307 x = pp.u+g.line[ii]*f; y = pp.v+g.line[ii+1]*f;
1308 q1.x = pp.x+(x*co+y*si)/2; q1.y = pp.y+(y*co-x*si)/2;
1309 line_draw(q0,q1,d);
1310 }
1311 else
1312 { // normal line
1313 mreal x,y;
1314 x = pp.u+g.line[ii]*f; y = pp.v+g.line[ii+1]*f;
1315 q0.x = pp.x+(x*co+y*si)/2; q0.y = pp.y+(y*co-x*si)/2; ii+=2;
1316 x = pp.u+g.line[ii]*f; y = pp.v+g.line[ii+1]*f;
1317 q1.x = pp.x+(x*co+y*si)/2; q1.y = pp.y+(y*co-x*si)/2;
1318 line_draw(q0,q1,d);
1319 }
1320 }
1321 }
1322 //-----------------------------------------------------------------------------
glyph_line(mreal phi,const mglPnt & pp,mreal f,bool solid,const mglDrawReg * d)1323 void mglCanvas::glyph_line(mreal phi, const mglPnt &pp, mreal f, bool solid, const mglDrawReg *d)
1324 {
1325 const mreal co=cos(phi*M_PI/180), si=sin(phi*M_PI/180);
1326 mglPnt q0=pp,q1=pp,q2=pp,q3=pp;
1327 q0.u=q0.v=q1.u=q1.v=q2.u=q2.v=q3.u=q3.v=NAN;
1328
1329 mreal dy = 0.004,x,y;
1330 x=pp.u; y=pp.v-dy; q0.x=pp.x+(x*co+y*si)/2; q0.y=pp.y+(y*co-x*si)/2;
1331 x=pp.u+f; y=pp.v-dy; q1.x=pp.x+(x*co+y*si)/2; q1.y=pp.y+(y*co-x*si)/2;
1332 x=pp.u; y=pp.v+dy; q2.x=pp.x+(x*co+y*si)/2; q2.y=pp.y+(y*co-x*si)/2;
1333 x=pp.u+f; y=pp.v+dy; q3.x=pp.x+(x*co+y*si)/2; q3.y=pp.y+(y*co-x*si)/2;
1334
1335 if(solid) quad_draw(q0,q1,q3,q2,d);
1336 else
1337 {
1338 line_draw(q0,q1,d); line_draw(q2,q1,d);
1339 line_draw(q0,q3,d); line_draw(q2,q3,d);
1340 }
1341 }
1342 //-----------------------------------------------------------------------------
setPp(mglPnt & q,const mglPoint & p)1343 long mglCanvas::setPp(mglPnt &q, const mglPoint &p)
1344 {
1345 q.xx=q.x=p.x; q.yy=q.y=p.y; q.zz=q.z=p.z;
1346 long k;
1347 #pragma omp critical(pnt)
1348 {k=Pnt.size(); MGL_PUSH(Pnt,q,mutexPnt);}
1349 return k;
1350 }
1351 //-----------------------------------------------------------------------------
arrow_draw(long n1,long n2,char st,float ll)1352 void mglCanvas::arrow_draw(long n1, long n2, char st, float ll)
1353 {
1354 const mglPnt &p1=Pnt[n1], &p2=Pnt[n2];
1355 mglPnt q=p1; //q.u=q.v=q.w=0;
1356
1357 mglPoint kl(p1.x-p2.x,p1.y-p2.y,p1.z-p2.z), kt, p0(p1.x,p1.y,p1.z), p;
1358 mreal d = hypot(kl.x,kl.y);
1359 if(d==0) return;
1360 kl /= d; kt = !kl;
1361 kl *= ll; kt *= ll;
1362
1363 Reserve(8);
1364 long k1,k2,k3,k4;
1365
1366 switch(st) // S,D -- cube, T -- sq.pyramid, I -- square, O -- sphere???, A,K,V -- cone???
1367 {
1368 case 'I':
1369 k1=setPp(q,p0+kt); k2=setPp(q,p0-kt); line_plot(k1,k2); break;
1370 case 'D':
1371 k1=setPp(q,p0+kl); k2=setPp(q,p0-kl); k3=setPp(q,p0+kt); k4=setPp(q,p0-kt);
1372 trig_plot(k1,k2,k3); trig_plot(k1,k2,k4); break;
1373 case 'S':
1374 k1=setPp(q,p0+kl+kt); k2=setPp(q,p0+kl-kt);
1375 k3=setPp(q,p0-kl-kt); k4=setPp(q,p0-kl+kt);
1376 quad_plot(k1,k2,k4,k3); break;
1377 case 'X':
1378 k1=setPp(q,p0+kl+kt); k2=setPp(q,p0+kl-kt);
1379 k3=setPp(q,p0-kl-kt); k4=setPp(q,p0-kl+kt);
1380 line_plot(k1,k3); line_plot(k2,k4); break;
1381 case 'T':
1382 k1=setPp(q,p0-kl+kt); k2=setPp(q,p0-kl-kt); k3=setPp(q,p0+kl);
1383 trig_plot(k1,k2,k3); break;
1384 case 'K':
1385 k1=setPp(q,p0+kt); k2=setPp(q,p0-kt); line_plot(k1,k2);
1386 case 'A':
1387 k1=setPp(q,p0-2.*kl+kt); k2=setPp(q,p0-2.*kl-kt); k3=setPp(q,p0-1.5*kl);
1388 trig_plot(n1,k3,k1); trig_plot(n1,k3,k2); break;
1389 case 'V':
1390 k1=setPp(q,p0+2.*kl+kt); k2=setPp(q,p0+2.*kl-kt); k3=setPp(q,p0+1.5*kl);
1391 trig_plot(n1,k3,k1); trig_plot(n1,k3,k2); break;
1392 case 'O': // let draw icosahedron
1393 {
1394 const int n = 12; k1=setPp(q,p0+kl);
1395 for(int i=1;i<=n;i++)
1396 {
1397 mreal u = 2*i*M_PI/n;
1398 k2 = k1; k1 = setPp(q,p0+kl*cos(u)+kt*sin(u));
1399 trig_plot(n1,k1,k2);
1400 }
1401 break;
1402 }
1403 }
1404 }
1405 //-----------------------------------------------------------------------------
arrow_plot_3d(long n1,long n2,char st,float ll)1406 void mglCanvas::arrow_plot_3d(long n1, long n2, char st, float ll)
1407 {
1408 const mglPnt &p1=Pnt[n1], &p2=Pnt[n2];
1409 mglPnt q=p1; //q.u=q.v=q.w=0;
1410
1411 mglPoint kl(p1.x-p2.x,p1.y-p2.y,p1.z-p2.z), kt, kz, p0(p1.x,p1.y,p1.z), p;
1412 if(kl.norm()==0) return;
1413 kl.Normalize(); kt = !kl; kz = kl^kt;
1414 kl *= ll; kt *= ll; kz *= ll;
1415
1416 Reserve(8);
1417 long k1,k2,k3,k4,k5, k6,k7,k8;
1418
1419 switch(st) // S,D -- cube, T -- sq.pyramid, I -- square, O -- sphere???, A,K,V -- cone???
1420 {
1421 case 'I':
1422 k1=setPp(q,p0+kt); k2=setPp(q,p0+kz);
1423 k3=setPp(q,p0-kt); k4=setPp(q,p0-kz);
1424 quad_plot(k1,k2,k4,k3); break;
1425 case 'D':
1426 k1=setPp(q,p0+kl); k2=setPp(q,p0-kl); k5=k3=setPp(q,p0+kt);
1427 k4=setPp(q,p0+kz); trig_plot(k1,k3,k4); trig_plot(k2,k3,k4); k3=k4;
1428 k4=setPp(q,p0-kt); trig_plot(k1,k3,k4); trig_plot(k2,k3,k4); k3=k4;
1429 k4=setPp(q,p0-kz); trig_plot(k1,k3,k4); trig_plot(k2,k3,k4); k3=k4;
1430 trig_plot(k1,k3,k5); trig_plot(k2,k3,k5); break;
1431 case 'S':
1432 k1=setPp(q,p0+kl+kt); k2=setPp(q,p0+kl+kz); k3=setPp(q,p0+kl-kt); k4=setPp(q,p0+kl-kz);
1433 k5=setPp(q,p0-kl+kt); k6=setPp(q,p0-kl+kz); k7=setPp(q,p0-kl-kt); k8=setPp(q,p0-kl-kz);
1434 quad_plot(k1,k2,k4,k3); quad_plot(k1,k2,k5,k6); quad_plot(k3,k2,k7,k6);
1435 quad_plot(k1,k4,k5,k8); quad_plot(k3,k4,k7,k8); quad_plot(k5,k6,k8,k7); break;
1436 case 'X':
1437 k1=setPp(q,p0+kl+kt); k2=setPp(q,p0+kl+kz); k3=setPp(q,p0+kl-kt); k4=setPp(q,p0+kl-kz);
1438 k5=setPp(q,p0-kl+kt); k6=setPp(q,p0-kl+kz); k7=setPp(q,p0-kl-kt); k8=setPp(q,p0-kl-kz);
1439 line_plot(k1,k7); line_plot(k2,k8); line_plot(k3,k5); line_plot(k4,k6); break;
1440 case 'T':
1441 k1=setPp(q,p0-kl+kt); k2=setPp(q,p0-kl+kz); k3=setPp(q,p0-kl-kt);
1442 k4=setPp(q,p0-kl-kz); k5=setPp(q,p0+kl);
1443 trig_plot(k1,k2,k5); trig_plot(k2,k3,k5);
1444 trig_plot(k3,k4,k5); trig_plot(k1,k4,k5); break;
1445 case 'K':
1446 k1=setPp(q,p0+kt); k2=setPp(q,p0+kz);
1447 k3=setPp(q,p0-kt); k4=setPp(q,p0-kz); quad_plot(k1,k2,k4,k3);
1448 case 'A':
1449 k1=setPp(q,p0-2.*kl+kt); k2=setPp(q,p0-2.*kl+kz); k3=setPp(q,p0-2.*kl-kt);
1450 k4=setPp(q,p0-2.*kl-kz); k5=setPp(q,p0-1.5*kl);
1451 trig_plot(n1,k5,k1); trig_plot(n1,k5,k2);
1452 trig_plot(n1,k5,k3); trig_plot(n1,k5,k4); break;
1453 case 'V':
1454 k1=setPp(q,p0+2.*kl+kt); k2=setPp(q,p0+2.*kl+kz); k3=setPp(q,p0+2.*kl-kt);
1455 k4=setPp(q,p0+2.*kl-kz); k5=setPp(q,p0+1.5*kl);
1456 trig_plot(n1,k5,k1); trig_plot(n1,k5,k2);
1457 trig_plot(n1,k5,k3); trig_plot(n1,k5,k4); break;
1458 case 'O': // let draw icosahedron
1459 {
1460 const int n = 12, m = n/2; Reserve(n*m);
1461 long *nn=new long[2*n], n1=setPp(q,p0+kl), n2=setPp(q,p0-kl);
1462 mreal u,v,rr;
1463 for(long i=0;i<m;i++) for(long j=0;j<n;j++)
1464 {
1465 if(i>0 && i<m-1)
1466 {
1467 u = i*M_PI/(m-1.); v = 2*M_PI*j/(n-1.)-1; rr = sin(u);
1468 nn[j+n]=nn[j]; nn[j]=setPp(q,p0+kl*cos(u)+kt*rr*cos(v)+kz*rr*sin(v));
1469 }
1470 else if(i==0) nn[j] = n1;
1471 else if(i==m-1) { nn[j+n]=nn[j]; nn[j]=n2; }
1472 if(i*j>0) quad_plot(nn[j-1], nn[j], nn[j+n-1], nn[j+n]);
1473 }
1474 delete []nn; break;
1475 }
1476 }
1477 }
1478 //-----------------------------------------------------------------------------
quad_vis(const mglPnt & p1,const mglPnt & p2,const mglPnt & p3,const mglPnt & p4) const1479 bool mglCanvas::quad_vis(const mglPnt &p1, const mglPnt &p2, const mglPnt &p3, const mglPnt &p4) const
1480 {
1481 long y1,x1,y2,x2;
1482 mglPnt d1(p2-p1), d2(p3-p1), d3(p4+p1-p2-p3);
1483
1484 if(d1.x==0 && d1.y==0) return trig_vis(p1,p3,p4);
1485 if(d2.x==0 && d2.y==0) return trig_vis(p1,p2,p4);
1486
1487 x1 = long(mgl_min(mgl_min(p1.x,p2.x), mgl_min(p3.x,p4.x))); // bounding box
1488 y1 = long(mgl_min(mgl_min(p1.y,p2.y), mgl_min(p3.y,p4.y)));
1489 x2 = long(mgl_max(mgl_max(p1.x,p2.x), mgl_max(p3.x,p4.x)));
1490 y2 = long(mgl_max(mgl_max(p1.y,p2.y), mgl_max(p3.y,p4.y)));
1491 x1=mgl_imax(x1,0); x2=mgl_imin(x2,Width);
1492 y1=mgl_imax(y1,0); y2=mgl_imin(y2,Height);
1493 // if(x1>x2 || y1>y2) return;
1494
1495 const float dd = d1.x*d2.y-d1.y*d2.x;
1496 const float dsx =-4*(d2.y*d3.x - d2.x*d3.y)*d1.y;
1497 const float dsy = 4*(d2.y*d3.x - d2.x*d3.y)*d1.x;
1498
1499 const float x0 = p1.x, y0 = p1.y;
1500 bool vis = false;
1501 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1502 {
1503 float xx = (i-x0), yy = (j-y0), s;
1504 s = dsx*xx + dsy*yy + (dd+d3.y*xx-d3.x*yy)*(dd+d3.y*xx-d3.x*yy);
1505 if(s>=0)
1506 {
1507 s = sqrt(s);
1508 float qu = d3.x*yy - d3.y*xx + dd + s;
1509 float qv = d3.y*xx - d3.x*yy + dd + s;
1510 float u = 2.f*(d2.y*xx - d2.x*yy)/qu;
1511 float v = 2.f*(d1.x*yy - d1.y*xx)/qv;
1512 if(u*(1.f-u)<0.f || v*(1.f-v)<0.f) // first root bad
1513 {
1514 qu = d3.x*yy - d3.y*xx + dd - s;
1515 qv = d3.y*xx - d3.x*yy + dd - s;
1516 // u = v = -1.f;
1517 u = 2.f*(d2.y*xx - d2.x*yy)/qu; v = 2.f*(d1.x*yy - d1.y*xx)/qv;
1518 if(u*(1.f-u)<0.f || v*(1.f-v)<0.f) continue; // second root bad
1519 }
1520 float zz = p1.z+d1.z*u+d2.z*v+d3.z*(u*v);
1521 if(zz>=Z[3*(i+Width*(Height-1-j))]-2) vis=true;
1522 }
1523 }
1524 return vis;
1525 }
1526 //-----------------------------------------------------------------------------
trig_vis(const mglPnt & p1,const mglPnt & p2,const mglPnt & p3) const1527 bool mglCanvas::trig_vis(const mglPnt &p1, const mglPnt &p2, const mglPnt &p3) const
1528 {
1529 long y1,x1,y2,x2;
1530 const mglPnt d1(p2-p1), d2(p3-p1);
1531
1532 const float tmp = d2.x*d1.y - d1.x*d2.y;
1533 if(fabs(tmp)<1e-5) return false; // points lies on the same line
1534 const float dyv =-d1.x/tmp, dxv = d1.y/tmp;
1535 const float dyu = d2.x/tmp, dxu =-d2.y/tmp;
1536
1537 x1 = long(mgl_min(p1.x<p2.x?p1.x:p2.x, p3.x)); // bounding box
1538 y1 = long(mgl_min(p1.y<p2.y?p1.y:p2.y, p3.y));
1539 x2 = long(mgl_max(p1.x>p2.x?p1.x:p2.x, p3.x));
1540 y2 = long(mgl_max(p1.y>p2.y?p1.y:p2.y, p3.y));
1541 x1=x1>0?x1:0; x2=x2<Width?x2:Width;
1542 y1=y1>0?y1:0; y2=y2<Height?y2:Height;
1543 // if(x1>x2 || y1>y2) return;
1544 // default normale
1545 const float x0 = p1.x, y0 = p1.y;
1546 bool vis=false;
1547 // provide additional height to be well visible on the surfaces
1548 for(long j=y1;j<=y2;j++) for(long i=x1;i<=x2;i++)
1549 {
1550 float xx = (i-x0), yy = (j-y0);
1551 float u = dxu*xx+dyu*yy, v = dxv*xx+dyv*yy;
1552 if(u<0 || v<0 || u+v>1) continue;
1553 float zz = p1.z+d1.z*u+d2.z*v;
1554 if(zz>=Z[3*(i+Width*(Height-1-j))]-2) vis=true;
1555 }
1556 return vis;
1557 }
1558 //-----------------------------------------------------------------------------
1559