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