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