1private import graph; 2 3private transform swap=(0,0,0,1,1,0); 4 5typedef bounds range(picture pic, real min, real max); 6 7range Range(bool automin=false, real min=-infinity, 8 bool automax=false, real max=infinity) 9{ 10 return new bounds(picture pic, real dmin, real dmax) { 11 // autoscale routine finds reasonable limits 12 bounds mz=autoscale(pic.scale.z.T(dmin), 13 pic.scale.z.T(dmax), 14 pic.scale.z.scale); 15 // If automin/max, use autoscale result, else 16 // if min/max is finite, use specified value, else 17 // use minimum/maximum data value 18 real pmin=automin ? pic.scale.z.Tinv(mz.min) : (finite(min) ? min : dmin); 19 real pmax=automax ? pic.scale.z.Tinv(mz.max) : (finite(max) ? max : dmax); 20 return bounds(pmin,pmax); 21 }; 22} 23 24range Automatic=Range(true,true); 25range Full=Range(); 26 27void image(frame f, real[][] data, pair initial, pair final, pen[] palette, 28 bool transpose=(initial.x < final.x && initial.y < final.y), 29 transform t=identity(), bool copy=true, bool antialias=false) 30{ 31 transform T=transpose ? swap : identity(); 32 _image(f,copy ? copy(data) : data,T*initial,T*final,palette,t*T,copy=false, 33 antialias=antialias); 34} 35 36void image(frame f, pen[][] data, pair initial, pair final, 37 bool transpose=(initial.x < final.x && initial.y < final.y), 38 transform t=identity(), bool copy=true, bool antialias=false) 39{ 40 transform T=transpose ? swap : identity(); 41 _image(f,copy ? copy(data) : data,T*initial,T*final,t*T,copy=false, 42 antialias=antialias); 43} 44 45// Reduce color palette to approximate range of data relative to "display" 46// range => errors of 1/palette.length in resulting color space. 47pen[] adjust(picture pic, real min, real max, real rmin, real rmax, 48 pen[] palette) 49{ 50 real dmin=pic.scale.z.T(min); 51 real dmax=pic.scale.z.T(max); 52 real delta=rmax-rmin; 53 if(delta > 0) { 54 real factor=palette.length/delta; 55 int minindex=floor(factor*(dmin-rmin)); 56 if(minindex < 0) minindex=0; 57 int maxindex=ceil(factor*(dmax-rmin)); 58 if(maxindex > palette.length) maxindex=palette.length; 59 if(minindex > 0 || maxindex < palette.length) 60 return palette[minindex:maxindex]; 61 } 62 return palette; 63} 64 65private real[] sequencereal; 66 67bounds image(picture pic=currentpicture, real[][] f, range range=Full, 68 pair initial, pair final, pen[] palette, 69 bool transpose=(initial.x < final.x && initial.y < final.y), 70 bool copy=true, bool antialias=false) 71{ 72 if(copy) f=copy(f); 73 if(copy) palette=copy(palette); 74 75 real m=min(f); 76 real M=max(f); 77 bounds bounds=range(pic,m,M); 78 real rmin=pic.scale.z.T(bounds.min); 79 real rmax=pic.scale.z.T(bounds.max); 80 palette=adjust(pic,m,M,rmin,rmax,palette); 81 82 // Crop data to allowed range and scale 83 if(range != Full || pic.scale.z.scale.T != identity || 84 pic.scale.z.postscale.T != identity) { 85 scalefcn T=pic.scale.z.T; 86 real m=bounds.min; 87 real M=bounds.max; 88 for(int i=0; i < f.length; ++i) 89 f[i]=map(new real(real x) {return T(min(max(x,m),M));},f[i]); 90 } 91 92 initial=Scale(pic,initial); 93 final=Scale(pic,final); 94 95 pic.addBox(initial,final); 96 97 transform T; 98 if(transpose) { 99 T=swap; 100 initial=T*initial; 101 final=T*final; 102 } 103 104 pic.add(new void(frame F, transform t) { 105 _image(F,f,initial,final,palette,t*T,copy=false,antialias=antialias); 106 },true); 107 return bounds; // Return bounds used for color space 108} 109 110bounds image(picture pic=currentpicture, real f(real, real), 111 range range=Full, pair initial, pair final, 112 int nx=ngraph, int ny=nx, pen[] palette, bool antialias=false) 113{ 114 // Generate data, taking scaling into account 115 real xmin=pic.scale.x.T(initial.x); 116 real xmax=pic.scale.x.T(final.x); 117 real ymin=pic.scale.y.T(initial.y); 118 real ymax=pic.scale.y.T(final.y); 119 real[][] data=new real[ny][nx]; 120 for(int j=0; j < ny; ++j) { 121 real y=pic.scale.y.Tinv(interp(ymin,ymax,(j+0.5)/ny)); 122 scalefcn Tinv=pic.scale.x.Tinv; 123 // Take center point of each bin 124 data[j]=sequence(new real(int i) { 125 return f(Tinv(interp(xmin,xmax,(i+0.5)/nx)),y); 126 },nx); 127 } 128 return image(pic,data,range,initial,final,palette,transpose=false, 129 copy=false,antialias=antialias); 130} 131 132void image(picture pic=currentpicture, pen[][] data, pair initial, pair final, 133 bool transpose=(initial.x < final.x && initial.y < final.y), 134 bool copy=true, bool antialias=false) 135{ 136 if(copy) data=copy(data); 137 138 initial=Scale(pic,initial); 139 final=Scale(pic,final); 140 141 pic.addBox(initial,final); 142 143 transform T; 144 if(transpose) { 145 T=swap; 146 initial=T*initial; 147 final=T*final; 148 } 149 150 pic.add(new void(frame F, transform t) { 151 _image(F,data,initial,final,t*T,copy=false,antialias=antialias); 152 },true); 153} 154 155void image(picture pic=currentpicture, pen f(int, int), int width, int height, 156 pair initial, pair final, 157 bool transpose=(initial.x < final.x && initial.y < final.y), 158 bool antialias=false) 159{ 160 initial=Scale(pic,initial); 161 final=Scale(pic,final); 162 163 pic.addBox(initial,final); 164 165 transform T; 166 if(transpose) { 167 T=swap; 168 int temp=width; 169 width=height; 170 height=temp; 171 initial=T*initial; 172 final=T*final; 173 } 174 175 pic.add(new void(frame F, transform t) { 176 _image(F,f,width,height,initial,final,t*T,antialias=antialias); 177 },true); 178} 179 180bounds image(picture pic=currentpicture, pair[] z, real[] f, 181 range range=Full, pen[] palette) 182{ 183 if(z.length != f.length) 184 abort("z and f arrays have different lengths"); 185 186 real m=min(f); 187 real M=max(f); 188 bounds bounds=range(pic,m,M); 189 real rmin=pic.scale.z.T(bounds.min); 190 real rmax=pic.scale.z.T(bounds.max); 191 192 palette=adjust(pic,m,M,rmin,rmax,palette); 193 rmin=max(rmin,m); 194 rmax=min(rmax,M); 195 196 // Crop data to allowed range and scale 197 if(range != Full || pic.scale.z.scale.T != identity || 198 pic.scale.z.postscale.T != identity) { 199 scalefcn T=pic.scale.z.T; 200 real m=bounds.min; 201 real M=bounds.max; 202 f=map(new real(real x) {return T(min(max(x,m),M));},f); 203 } 204 205 int[] edges={0,0,1}; 206 int N=palette.length-1; 207 208 int[][] trn=triangulate(z); 209 real step=rmax == rmin ? 0.0 : N/(rmax-rmin); 210 for(int i=0; i < trn.length; ++i) { 211 int[] trni=trn[i]; 212 int i0=trni[0], i1=trni[1], i2=trni[2]; 213 pen color(int i) {return palette[round((f[i]-rmin)*step)];} 214 gouraudshade(pic,z[i0]--z[i1]--z[i2]--cycle, 215 new pen[] {color(i0),color(i1),color(i2)},edges); 216 } 217 return bounds; // Return bounds used for color space 218} 219 220bounds image(picture pic=currentpicture, real[] x, real[] y, real[] f, 221 range range=Full, pen[] palette) 222{ 223 int n=x.length; 224 if(n != y.length) 225 abort("x and y arrays have different lengths"); 226 227 pair[] z=sequence(new pair(int i) {return (x[i],y[i]);},n); 228 return image(pic,z,f,range,palette); 229} 230 231// Construct a pen[] array from f using the specified palette. 232pen[] palette(real[] f, pen[] palette) 233{ 234 real Min=min(f); 235 real Max=max(f); 236 if(palette.length == 0) return new pen[]; 237 real step=Max == Min ? 0.0 : (palette.length-1)/(Max-Min); 238 return sequence(new pen(int i) {return palette[round((f[i]-Min)*step)];}, 239 f.length); 240} 241 242// Construct a pen[][] array from f using the specified palette. 243pen[][] palette(real[][] f, pen[] palette) 244{ 245 real Min=min(f); 246 real Max=max(f); 247 int n=f.length; 248 pen[][] p=new pen[n][]; 249 real step=(Max == Min) ? 0.0 : (palette.length-1)/(Max-Min); 250 for(int i=0; i < n; ++i) { 251 real[] fi=f[i]; 252 p[i]=sequence(new pen(int j) {return palette[round((fi[j]-Min)*step)];}, 253 f[i].length); 254 } 255 return p; 256} 257 258typedef ticks paletteticks(int sign=-1); 259 260paletteticks PaletteTicks(Label format="", ticklabel ticklabel=null, 261 bool beginlabel=true, bool endlabel=true, 262 int N=0, int n=0, real Step=0, real step=0, 263 pen pTick=nullpen, pen ptick=nullpen) 264{ 265 return new ticks(int sign=-1) { 266 format.align(sign > 0 ? RightSide : LeftSide); 267 return Ticks(sign,format,ticklabel,beginlabel,endlabel,N,n,Step,step, 268 true,true,extend=true,pTick,ptick); 269 }; 270} 271 272paletteticks PaletteTicks=PaletteTicks(); 273paletteticks NoTicks=new ticks(int sign=-1) {return NoTicks;}; 274 275void palette(picture pic=currentpicture, Label L="", bounds bounds, 276 pair initial, pair final, axis axis=Right, pen[] palette, 277 pen p=currentpen, paletteticks ticks=PaletteTicks, 278 bool copy=true, bool antialias=false) 279{ 280 real initialz=pic.scale.z.T(bounds.min); 281 real finalz=pic.scale.z.T(bounds.max); 282 bounds mz=autoscale(initialz,finalz,pic.scale.z.scale); 283 284 axisT axis; 285 axis(pic,axis); 286 real angle=degrees(axis.align.dir); 287 288 initial=Scale(pic,initial); 289 final=Scale(pic,final); 290 291 pair lambda=final-initial; 292 bool vertical=(floor((angle+45)/90) % 2 == 0); 293 pair perp,par; 294 295 if(vertical) {perp=E; par=N;} else {perp=N; par=E;} 296 297 path g=(final-dot(lambda,par)*par)--final; 298 path g2=initial--final-dot(lambda,perp)*perp; 299 300 if(sgn(dot(lambda,perp)*dot(axis.align.dir,perp)) == -1) { 301 path tmp=g; 302 g=g2; 303 g2=tmp; 304 } 305 306 if(copy) palette=copy(palette); 307 Label L=L.copy(); 308 if(L.defaultposition) L.position(0.5); 309 L.align(axis.align); 310 L.p(p); 311 if(vertical && L.defaulttransform) { 312 frame f; 313 add(f,Label(L.s,(0,0),L.p)); 314 if(length(max(f)-min(f)) > ylabelwidth*fontsize(L.p)) 315 L.transform(rotate(90)); 316 } 317 real[][] pdata={sequence(palette.length)}; 318 319 transform T; 320 pair Tinitial,Tfinal; 321 if(vertical) { 322 T=swap; 323 Tinitial=T*initial; 324 Tfinal=T*final; 325 } else { 326 Tinitial=initial; 327 Tfinal=final; 328 } 329 330 pic.add(new void(frame f, transform t) { 331 _image(f,pdata,Tinitial,Tfinal,palette,t*T,copy=false, 332 antialias=antialias); 333 },true); 334 335 ticklocate locate=ticklocate(initialz,finalz,pic.scale.z,mz.min,mz.max); 336 axis(pic,L,g,g2,p,ticks(sgn(axis.side.x*dot(lambda,par))),locate,mz.divisor, 337 true); 338 339 pic.add(new void(frame f, transform t) { 340 pair Z0=t*initial; 341 pair Z1=t*final; 342 draw(f,Z0--(Z0.x,Z1.y)--Z1--(Z1.x,Z0.y)--cycle,p); 343 },true); 344 345 pic.addBox(initial,final); 346} 347 348// A grayscale palette 349pen[] Grayscale(int NColors=256) 350{ 351 real ninv=1.0/(NColors-1.0); 352 return sequence(new pen(int i) {return gray(i*ninv);},NColors); 353} 354 355// A color wheel palette 356pen[] Wheel(int NColors=32766) 357{ 358 if(settings.gray) return Grayscale(NColors); 359 360 int nintervals=6; 361 if(NColors <= nintervals) NColors=nintervals+1; 362 int n=-quotient(NColors,-nintervals); 363 364 pen[] Palette; 365 366 Palette=new pen[n*nintervals]; 367 real ninv=1.0/n; 368 369 for(int i=0; i < n; ++i) { 370 real ininv=i*ninv; 371 real ininv1=1.0-ininv; 372 Palette[i]=rgb(1.0,0.0,ininv); 373 Palette[n+i]=rgb(ininv1,0.0,1.0); 374 Palette[2n+i]=rgb(0.0,ininv,1.0); 375 Palette[3n+i]=rgb(0.0,1.0,ininv1); 376 Palette[4n+i]=rgb(ininv,1.0,0.0); 377 Palette[5n+i]=rgb(1.0,ininv1,0.0); 378 } 379 return Palette; 380} 381 382// A rainbow palette 383pen[] Rainbow(int NColors=32766) 384{ 385 if(settings.gray) return Grayscale(NColors); 386 387 int offset=1; 388 int nintervals=5; 389 if(NColors <= nintervals) NColors=nintervals+1; 390 int n=-quotient(NColors-1,-nintervals); 391 392 pen[] Palette; 393 394 Palette=new pen[n*nintervals+offset]; 395 real ninv=1.0/n; 396 397 for(int i=0; i < n; ++i) { 398 real ininv=i*ninv; 399 real ininv1=1.0-ininv; 400 Palette[i]=rgb(ininv1,0.0,1.0); 401 Palette[n+i]=rgb(0.0,ininv,1.0); 402 Palette[2n+i]=rgb(0.0,1.0,ininv1); 403 Palette[3n+i]=rgb(ininv,1.0,0.0); 404 Palette[4n+i]=rgb(1.0,ininv1,0.0); 405 } 406 Palette[4n+n]=rgb(1.0,0.0,0.0); 407 408 return Palette; 409} 410 411private pen[] BWRainbow(int NColors, bool two) 412{ 413 if(settings.gray) return Grayscale(NColors); 414 415 int offset=1; 416 int nintervals=6; 417 int divisor=3; 418 419 if(two) nintervals += 6; 420 421 int Nintervals=nintervals*divisor; 422 if(NColors <= Nintervals) NColors=Nintervals+1; 423 int num=NColors-offset; 424 int n=-quotient(num,-Nintervals)*divisor; 425 NColors=n*nintervals+offset; 426 427 pen[] Palette; 428 429 Palette=new pen[NColors]; 430 real ninv=1.0/n; 431 432 int k=0; 433 434 if(two) { 435 for(int i=0; i < n; ++i) { 436 real ininv=i*ninv; 437 real ininv1=1.0-ininv; 438 Palette[i]=rgb(ininv1,0.0,1.0); 439 Palette[n+i]=rgb(0.0,ininv,1.0); 440 Palette[2n+i]=rgb(0.0,1.0,ininv1); 441 Palette[3n+i]=rgb(ininv,1.0,0.0); 442 Palette[4n+i]=rgb(1.0,ininv1,0.0); 443 Palette[5n+i]=rgb(1.0,0.0,ininv); 444 } 445 k += 6n; 446 } 447 448 if(two) 449 for(int i=0; i < n; ++i) 450 Palette[k+i]=rgb(1.0-i*ninv,0.0,1.0); 451 else { 452 int n3=-quotient(n,-3); 453 int n23=2*n3; 454 real third=n3*ninv; 455 real twothirds=n23*ninv; 456 for(int i=0; i < n3; ++i) { 457 real ininv=i*ninv; 458 Palette[k+i]=rgb(ininv,0.0,ininv); 459 Palette[k+n3+i]=rgb(third,0.0,third+ininv); 460 Palette[k+n23+i]=rgb(third-ininv,0.0,twothirds+ininv); 461 } 462 } 463 k += n; 464 465 for(int i=0; i < n; ++i) { 466 real ininv=i*ninv; 467 real ininv1=1.0-ininv; 468 Palette[k+i]=rgb(0.0,ininv,1.0); 469 Palette[k+n+i]=rgb(0.0,1.0,ininv1); 470 Palette[k+2n+i]=rgb(ininv,1.0,0.0); 471 Palette[k+3n+i]=rgb(1.0,ininv1,0.0); 472 Palette[k+4n+i]=rgb(1.0,ininv,ininv); 473 } 474 Palette[k+5n]=rgb(1.0,1.0,1.0); 475 476 return Palette; 477} 478 479// Quantize palette to exactly n values 480pen[] quantize(pen[] Palette, int n) 481{ 482 if(Palette.length == 0) abort("cannot quantize empty palette"); 483 if(n <= 1) abort("palette must contain at least two pens"); 484 real step=(Palette.length-1)/(n-1); 485 return sequence(new pen(int i) { 486 return Palette[round(i*step)]; 487 },n); 488} 489 490// A rainbow palette tapering off to black/white at the spectrum ends, 491pen[] BWRainbow(int NColors=32761) 492{ 493 return BWRainbow(NColors,false); 494} 495 496// A double rainbow palette tapering off to black/white at the spectrum ends, 497// with a linearly scaled intensity. 498pen[] BWRainbow2(int NColors=32761) 499{ 500 pen[] Palette=BWRainbow(NColors,true); 501 int n=Palette.length; 502 real ninv=1.0/n; 503 for(int i=0; i < n; ++i) 504 Palette[i]=i*ninv*Palette[i]; 505 return Palette; 506} 507 508//A palette varying linearly over the specified array of pens, using 509// NColors in each interpolation interval. 510pen[] Gradient(int NColors=256 ... pen[] p) 511{ 512 pen[] P; 513 if(p.length < 2) abort("at least 2 colors must be specified"); 514 real step=NColors > 1 ? (1/(NColors-1)) : 1; 515 for(int i=0; i < p.length-1; ++i) { 516 pen begin=p[i]; 517 pen end=p[i+1]; 518 P.append(sequence(new pen(int j) { 519 return interp(begin,end,j*step); 520 },NColors)); 521 } 522 return P; 523} 524 525pen[] cmyk(pen[] Palette) 526{ 527 int n=Palette.length; 528 for(int i=0; i < n; ++i) 529 Palette[i]=cmyk(Palette[i]); 530 return Palette; 531} 532