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 int m=n > 0 ? f[0].length : 0; 249 pen[][] p=new pen[n][m]; 250 real step=(Max == Min) ? 0.0 : (palette.length-1)/(Max-Min); 251 for(int i=0; i < n; ++i) { 252 real[] fi=f[i]; 253 p[i]=sequence(new pen(int j) {return palette[round((fi[j]-Min)*step)];},m); 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(); 273 274void palette(picture pic=currentpicture, Label L="", bounds bounds, 275 pair initial, pair final, axis axis=Right, pen[] palette, 276 pen p=currentpen, paletteticks ticks=PaletteTicks, 277 bool copy=true, bool antialias=false) 278{ 279 real initialz=pic.scale.z.T(bounds.min); 280 real finalz=pic.scale.z.T(bounds.max); 281 bounds mz=autoscale(initialz,finalz,pic.scale.z.scale); 282 283 axisT axis; 284 axis(pic,axis); 285 real angle=degrees(axis.align.dir); 286 287 initial=Scale(pic,initial); 288 final=Scale(pic,final); 289 290 pair lambda=final-initial; 291 bool vertical=(floor((angle+45)/90) % 2 == 0); 292 pair perp,par; 293 294 if(vertical) {perp=E; par=N;} else {perp=N; par=E;} 295 296 path g=(final-dot(lambda,par)*par)--final; 297 path g2=initial--final-dot(lambda,perp)*perp; 298 299 if(sgn(dot(lambda,perp)*dot(axis.align.dir,perp)) == -1) { 300 path tmp=g; 301 g=g2; 302 g2=tmp; 303 } 304 305 if(copy) palette=copy(palette); 306 Label L=L.copy(); 307 if(L.defaultposition) L.position(0.5); 308 L.align(axis.align); 309 L.p(p); 310 if(vertical && L.defaulttransform) { 311 frame f; 312 add(f,Label(L.s,(0,0),L.p)); 313 if(length(max(f)-min(f)) > ylabelwidth*fontsize(L.p)) 314 L.transform(rotate(90)); 315 } 316 real[][] pdata={sequence(palette.length)}; 317 318 transform T; 319 pair Tinitial,Tfinal; 320 if(vertical) { 321 T=swap; 322 Tinitial=T*initial; 323 Tfinal=T*final; 324 } else { 325 Tinitial=initial; 326 Tfinal=final; 327 } 328 329 pic.add(new void(frame f, transform t) { 330 _image(f,pdata,Tinitial,Tfinal,palette,t*T,copy=false, 331 antialias=antialias); 332 },true); 333 334 ticklocate locate=ticklocate(initialz,finalz,pic.scale.z,mz.min,mz.max); 335 axis(pic,L,g,g2,p,ticks(sgn(axis.side.x*dot(lambda,par))),locate,mz.divisor, 336 true); 337 338 pic.add(new void(frame f, transform t) { 339 pair Z0=t*initial; 340 pair Z1=t*final; 341 draw(f,Z0--(Z0.x,Z1.y)--Z1--(Z1.x,Z0.y)--cycle,p); 342 },true); 343 344 pic.addBox(initial,final); 345} 346 347// A grayscale palette 348pen[] Grayscale(int NColors=256) 349{ 350 real ninv=1.0/(NColors-1.0); 351 return sequence(new pen(int i) {return gray(i*ninv);},NColors); 352} 353 354// A color wheel palette 355pen[] Wheel(int NColors=32766) 356{ 357 if(settings.gray) return Grayscale(NColors); 358 359 int nintervals=6; 360 int n=-quotient(NColors,-nintervals); 361 362 pen[] Palette; 363 if(n == 0) return Palette; 364 365 Palette=new pen[n*nintervals]; 366 real ninv=1.0/n; 367 368 for(int i=0; i < n; ++i) { 369 real ininv=i*ninv; 370 real ininv1=1.0-ininv; 371 Palette[i]=rgb(1.0,0.0,ininv); 372 Palette[n+i]=rgb(ininv1,0.0,1.0); 373 Palette[2n+i]=rgb(0.0,ininv,1.0); 374 Palette[3n+i]=rgb(0.0,1.0,ininv1); 375 Palette[4n+i]=rgb(ininv,1.0,0.0); 376 Palette[5n+i]=rgb(1.0,ininv1,0.0); 377 } 378 return Palette; 379} 380 381// A rainbow palette 382pen[] Rainbow(int NColors=32766) 383{ 384 if(settings.gray) return Grayscale(NColors); 385 386 int offset=1; 387 int nintervals=5; 388 int n=-quotient(NColors-1,-nintervals); 389 390 pen[] Palette; 391 if(n == 0) return Palette; 392 393 Palette=new pen[n*nintervals+offset]; 394 real ninv=1.0/n; 395 396 for(int i=0; i < n; ++i) { 397 real ininv=i*ninv; 398 real ininv1=1.0-ininv; 399 Palette[i]=rgb(ininv1,0.0,1.0); 400 Palette[n+i]=rgb(0.0,ininv,1.0); 401 Palette[2n+i]=rgb(0.0,1.0,ininv1); 402 Palette[3n+i]=rgb(ininv,1.0,0.0); 403 Palette[4n+i]=rgb(1.0,ininv1,0.0); 404 } 405 Palette[4n+n]=rgb(1.0,0.0,0.0); 406 407 return Palette; 408} 409 410private pen[] BWRainbow(int NColors, bool two) 411{ 412 if(settings.gray) return Grayscale(NColors); 413 414 int offset=1; 415 int nintervals=6; 416 int divisor=3; 417 418 if(two) nintervals += 6; 419 420 int num=NColors-offset; 421 int n=-quotient(num,-nintervals*divisor)*divisor; 422 NColors=n*nintervals+offset; 423 424 pen[] Palette; 425 if(n == 0) return Palette; 426 427 Palette=new pen[NColors]; 428 real ninv=1.0/n; 429 430 int k=0; 431 432 if(two) { 433 for(int i=0; i < n; ++i) { 434 real ininv=i*ninv; 435 real ininv1=1.0-ininv; 436 Palette[i]=rgb(ininv1,0.0,1.0); 437 Palette[n+i]=rgb(0.0,ininv,1.0); 438 Palette[2n+i]=rgb(0.0,1.0,ininv1); 439 Palette[3n+i]=rgb(ininv,1.0,0.0); 440 Palette[4n+i]=rgb(1.0,ininv1,0.0); 441 Palette[5n+i]=rgb(1.0,0.0,ininv); 442 } 443 k += 6n; 444 } 445 446 if(two) 447 for(int i=0; i < n; ++i) 448 Palette[k+i]=rgb(1.0-i*ninv,0.0,1.0); 449 else { 450 int n3=-quotient(n,-3); 451 int n23=2*n3; 452 real third=n3*ninv; 453 real twothirds=n23*ninv; 454 for(int i=0; i < n3; ++i) { 455 real ininv=i*ninv; 456 Palette[k+i]=rgb(ininv,0.0,ininv); 457 Palette[k+n3+i]=rgb(third,0.0,third+ininv); 458 Palette[k+n23+i]=rgb(third-ininv,0.0,twothirds+ininv); 459 } 460 } 461 k += n; 462 463 for(int i=0; i < n; ++i) { 464 real ininv=i*ninv; 465 real ininv1=1.0-ininv; 466 Palette[k+i]=rgb(0.0,ininv,1.0); 467 Palette[k+n+i]=rgb(0.0,1.0,ininv1); 468 Palette[k+2n+i]=rgb(ininv,1.0,0.0); 469 Palette[k+3n+i]=rgb(1.0,ininv1,0.0); 470 Palette[k+4n+i]=rgb(1.0,ininv,ininv); 471 } 472 Palette[k+5n]=rgb(1.0,1.0,1.0); 473 474 return Palette; 475} 476 477// Quantize palette to exactly n values 478pen[] quantize(pen[] Palette, int n) 479{ 480 if(Palette.length == 0) abort("cannot quantize empty palette"); 481 if(n <= 1) abort("palette must contain at least two pens"); 482 real step=(Palette.length-1)/(n-1); 483 return sequence(new pen(int i) { 484 return Palette[round(i*step)]; 485 },n); 486} 487 488// A rainbow palette tapering off to black/white at the spectrum ends, 489pen[] BWRainbow(int NColors=32761) 490{ 491 return BWRainbow(NColors,false); 492} 493 494// A double rainbow palette tapering off to black/white at the spectrum ends, 495// with a linearly scaled intensity. 496pen[] BWRainbow2(int NColors=32761) 497{ 498 pen[] Palette=BWRainbow(NColors,true); 499 int n=Palette.length; 500 real ninv=1.0/n; 501 for(int i=0; i < n; ++i) 502 Palette[i]=i*ninv*Palette[i]; 503 return Palette; 504} 505 506//A palette varying linearly over the specified array of pens, using 507// NColors in each interpolation interval. 508pen[] Gradient(int NColors=256 ... pen[] p) 509{ 510 pen[] P; 511 if(p.length < 2) abort("at least 2 colors must be specified"); 512 real step=NColors > 1 ? (1/(NColors-1)) : 1; 513 for(int i=0; i < p.length-1; ++i) { 514 pen begin=p[i]; 515 pen end=p[i+1]; 516 P.append(sequence(new pen(int j) { 517 return interp(begin,end,j*step); 518 },NColors)); 519 } 520 return P; 521} 522 523pen[] cmyk(pen[] Palette) 524{ 525 int n=Palette.length; 526 for(int i=0; i < n; ++i) 527 Palette[i]=cmyk(Palette[i]); 528 return Palette; 529} 530