1 /* ----------------------------------------------------------------------
2    SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3    http://sparta.sandia.gov
4    Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5    Sandia National Laboratories
6 
7    Copyright (2014) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14 
15 /* ----------------------------------------------------------------------
16    Contributing author: Nathan Fabian (Sandia)
17 ------------------------------------------------------------------------- */
18 
19 #include "mpi.h"
20 #include "math.h"
21 #include "ctype.h"
22 #include "stdlib.h"
23 #include "string.h"
24 #include "image.h"
25 #include "math_extra.h"
26 #include "update.h"
27 #include "random_mars.h"
28 #include "random_knuth.h"
29 #include "math_const.h"
30 #include "error.h"
31 #include "memory.h"
32 
33 #ifdef SPARTA_JPEG
34 #include "jpeglib.h"
35 #endif
36 
37 #ifdef SPARTA_PNG
38 #include <png.h>
39 #include <zlib.h>
40 #include <setjmp.h>
41 #include "version.h"
42 #endif
43 
44 using namespace SPARTA_NS;
45 using namespace MathConst;
46 
47 #define NCOLORS 140
48 #define NELEMENTS 109
49 #define BIG 1.0e20
50 #define EPSILON 1.0e-6
51 
52 enum{NUMERIC,MINVALUE,MAXVALUE};
53 enum{CONTINUOUS,DISCRETE,SEQUENTIAL};
54 enum{ABSOLUTE,FRACTIONAL};
55 enum{NO,YES};
56 
57 /* ---------------------------------------------------------------------- */
58 
Image(SPARTA * sparta,int nmap_caller)59 Image::Image(SPARTA *sparta, int nmap_caller) : Pointers(sparta)
60 {
61   MPI_Comm_rank(world,&me);
62   MPI_Comm_size(world,&nprocs);
63 
64   // defaults for 3d viz
65 
66   width = height = 512;
67   theta = 60.0 * MY_PI/180.0;
68   phi = 30.0 * MY_PI/180.0;
69   zoom = 1.0;
70   persp = 0.0;
71   shiny = 1.0;
72   ssao = NO;
73 
74   up[0] = 0.0;
75   up[1] = 0.0;
76   up[2] = 1.0;
77 
78   // colors
79 
80   ncolors = 0;
81   username = NULL;
82   userrgb = NULL;
83 
84   boxcolor = color2rgb("yellow");
85   background[0] = background[1] = background[2] = 0;
86 
87   // define nmap colormaps, all with default settings
88 
89   nmap = nmap_caller;
90   maps = new ColorMap*[nmap];
91   for (int i = 0; i < nmap; i++)
92     maps[i] = new ColorMap(sparta,this);
93 
94   // static parameters
95 
96   FOV = MY_PI/6.0;              // 30 degrees
97   ambientColor[0] = 0.0;
98   ambientColor[1] = 0.0;
99   ambientColor[2] = 0.0;
100 
101   keyLightPhi = -MY_PI4;        // -45 degrees
102   keyLightTheta = MY_PI/6.0;    // 30 degrees
103   keyLightColor[0] = 0.9;
104   keyLightColor[1] = 0.9;
105   keyLightColor[2] = 0.9;
106 
107   fillLightPhi = MY_PI/6.0;     // 30 degrees
108   fillLightTheta = 0;
109   fillLightColor[0] = 0.45;
110   fillLightColor[1] = 0.45;
111   fillLightColor[2] = 0.45;
112 
113   backLightPhi = MY_PI;         // 180 degrees
114   backLightTheta = MY_PI/12.0;  // 15 degrees
115   backLightColor[0] = 0.9;
116   backLightColor[1] = 0.9;
117   backLightColor[2] = 0.9;
118 
119   random = NULL;
120 
121   // MPI_Gatherv vectors
122 
123   recvcounts = NULL;
124   displs = NULL;
125 }
126 
127 /* ---------------------------------------------------------------------- */
128 
~Image()129 Image::~Image()
130 {
131   for (int i = 0; i < nmap; i++) delete maps[i];
132   delete [] maps;
133 
134   for (int i = 0; i < ncolors; i++) delete [] username[i];
135   memory->sfree(username);
136   memory->destroy(userrgb);
137 
138   memory->destroy(depthBuffer);
139   memory->destroy(surfaceBuffer);
140   memory->destroy(imageBuffer);
141   memory->destroy(depthcopy);
142   memory->destroy(surfacecopy);
143   memory->destroy(rgbcopy);
144 
145   if (random) delete random;
146 
147   memory->destroy(recvcounts);
148   memory->destroy(displs);
149 }
150 
151 /* ----------------------------------------------------------------------
152    allocate image and depth buffers
153    called after image size is set
154 ------------------------------------------------------------------------- */
155 
buffers()156 void Image::buffers()
157 {
158   npixels = width * height;
159   memory->create(depthBuffer,npixels,"image:depthBuffer");
160   memory->create(surfaceBuffer,2*npixels,"image:surfaceBuffer");
161   memory->create(imageBuffer,3*npixels,"image:imageBuffer");
162   memory->create(depthcopy,npixels,"image:depthcopy");
163   memory->create(surfacecopy,2*npixels,"image:surfacecopy");
164   memory->create(rgbcopy,3*npixels,"image:rgbcopy");
165 }
166 
167 /* ----------------------------------------------------------------------
168    reset view parameters
169    called once if view is STATIC
170    called before every render if view is DYNAMIC
171 ------------------------------------------------------------------------- */
172 
view_params(double boxxlo,double boxxhi,double boxylo,double boxyhi,double boxzlo,double boxzhi)173 void Image::view_params(double boxxlo, double boxxhi, double boxylo,
174 			double boxyhi, double boxzlo, double boxzhi)
175 {
176   // camDir points at the camera, view direction = -camDir
177 
178   camDir[0] = sin(theta)*cos(phi);
179   camDir[1] = sin(theta)*sin(phi);
180   camDir[2] = cos(theta);
181 
182   // normalize up vector
183 
184   if (up[0] == 0.0 && up[1] == 0.0 && up[2] == 0.0)
185     error->all(FLERR,"Invalid image up vector");
186   MathExtra::norm3(up);
187 
188   // adjust camDir by epsilon if camDir and up are parallel
189   // do this by tweaking view direction, not up direction
190   // try to insure continuous images as changing view passes thru up
191   // sufficient to handle common cases where theta = 0 or 180 is degenerate?
192 
193   double dot = MathExtra::dot3(up,camDir);
194   if (fabs(dot) > 1.0-EPSILON) {
195     if (theta == 0.0) {
196       camDir[0] = sin(EPSILON)*cos(phi);
197       camDir[1] = sin(EPSILON)*sin(phi);
198       camDir[2] = cos(EPSILON);
199     } else if (theta == MY_PI) {
200       camDir[0] = sin(theta-EPSILON)*cos(phi);
201       camDir[1] = sin(theta-EPSILON)*sin(phi);
202       camDir[2] = cos(theta-EPSILON);
203     } else {
204       camDir[0] = sin(theta+EPSILON)*cos(phi);
205       camDir[1] = sin(theta+EPSILON)*sin(phi);
206       camDir[2] = cos(theta+EPSILON);
207     }
208   }
209 
210   // camUp = camDir x (Up x camDir)
211 
212   MathExtra::cross3(up,camDir,camRight);
213   MathExtra::norm3(camRight);
214   MathExtra::cross3(camDir,camRight,camUp);
215   if (camUp[0] == 0.0 && camUp[1] == 0.0 && camUp[2] == 0.0)
216     error->all(FLERR,"Invalid image up vector");
217   MathExtra::norm3(camUp);
218 
219   // zdist = camera distance = function of zoom & bounding box
220   // camPos = camera position = function of camDir and zdist
221 
222   double delx = 2.0*(boxxhi-boxxlo);
223   double dely = 2.0*(boxyhi-boxylo);
224   double delz = 2.0*(boxzhi-boxzlo);
225   double maxdel = MAX(delx,dely);
226   maxdel = MAX(maxdel,delz);
227 
228   zdist = maxdel;
229   zdist /= tan(FOV);
230   zdist += 0.5 * (delx*camDir[0] + dely*camDir[1] + delz*camDir[2]);
231   zdist /= zoom;
232 
233   camPos[0] = camDir[0] * zdist;
234   camPos[1] = camDir[1] * zdist;
235   camPos[2] = camDir[2] * zdist;
236 
237   // light directions in terms of -camDir = z
238 
239   keyLightDir[0] = cos(keyLightTheta) * sin(keyLightPhi);
240   keyLightDir[1] = sin(keyLightTheta);
241   keyLightDir[2] = cos(keyLightTheta) * cos(keyLightPhi);
242 
243   fillLightDir[0] = cos(fillLightTheta) * sin(fillLightPhi);
244   fillLightDir[1] = sin(fillLightTheta);
245   fillLightDir[2] = cos(fillLightTheta) * cos(fillLightPhi);
246 
247   backLightDir[0] = cos(backLightTheta) * sin(backLightPhi);
248   backLightDir[1] = sin(backLightTheta);
249   backLightDir[2] = cos(backLightTheta) * cos(backLightPhi);
250 
251   keyHalfDir[0] = 0 + keyLightDir[0];
252   keyHalfDir[1] = 0 + keyLightDir[1];
253   keyHalfDir[2] = 1 + keyLightDir[2];
254   MathExtra::norm3(keyHalfDir);
255 
256   // adjust shinyness of the reflection
257 
258   specularHardness = 16.0 * shiny;
259   specularIntensity = shiny;
260 
261   // adjust strength of the SSAO
262 
263   if (ssao) {
264     if (!random) {
265       random = new RanKnuth(update->ranmaster->uniform());
266       double seed = update->ranmaster->uniform();
267       random->reset(seed,me,100);
268     }
269     SSAORadius = maxdel * 0.05 * ssaoint;
270     SSAOSamples = static_cast<int> (8.0 + 32.0*ssaoint);
271     SSAOJitter = MY_PI / 12;
272     ambientColor[0] = 0.5;
273     ambientColor[1] = 0.5;
274     ambientColor[2] = 0.5;
275   }
276 
277   // param for rasterizing spheres
278 
279   tanPerPixel = -(maxdel / (double) height);
280 }
281 
282 /* ----------------------------------------------------------------------
283    initialize image to background color and depth buffer
284    no need to init surfaceBuffer, since will be based on depth
285 ------------------------------------------------------------------------- */
286 
clear()287 void Image::clear()
288 {
289   int red = background[0];
290   int green = background[1];
291   int blue = background[2];
292 
293   int ix,iy;
294   for (iy = 0; iy < height; iy ++)
295     for (ix = 0; ix < width; ix ++) {
296       imageBuffer[iy * width * 3 + ix * 3 + 0] = red;
297       imageBuffer[iy * width * 3 + ix * 3 + 1] = green;
298       imageBuffer[iy * width * 3 + ix * 3 + 2] = blue;
299       depthBuffer[iy * width + ix] = -1;
300     }
301 }
302 
303 /* ----------------------------------------------------------------------
304    merge image from each processor into one composite image
305    done pixel by pixel, respecting depth buffer
306    hi procs send to lo procs, cascading down logarithmically
307 ------------------------------------------------------------------------- */
308 
merge()309 void Image::merge()
310 {
311   MPI_Request requests[3];
312   MPI_Status statuses[3];
313 
314   int nhalf = 1;
315   while (nhalf < nprocs) nhalf *= 2;
316   nhalf /= 2;
317 
318   while (nhalf) {
319     if (me < nhalf && me+nhalf < nprocs) {
320       MPI_Irecv(rgbcopy,npixels*3,MPI_BYTE,me+nhalf,0,world,&requests[0]);
321       MPI_Irecv(depthcopy,npixels,MPI_DOUBLE,me+nhalf,0,world,&requests[1]);
322       if (ssao)
323 	MPI_Irecv(surfacecopy,npixels*2,MPI_DOUBLE,
324 		  me+nhalf,0,world,&requests[2]);
325       if (ssao) MPI_Waitall(3,requests,statuses);
326       else MPI_Waitall(2,requests,statuses);
327 
328       for (int i = 0; i < npixels; i++) {
329         if (depthBuffer[i] < 0 || (depthcopy[i] >= 0 &&
330 				   depthcopy[i] < depthBuffer[i])) {
331           depthBuffer[i] = depthcopy[i];
332           imageBuffer[i*3+0] = rgbcopy[i*3+0];
333           imageBuffer[i*3+1] = rgbcopy[i*3+1];
334           imageBuffer[i*3+2] = rgbcopy[i*3+2];
335           if (ssao) {
336             surfaceBuffer[i*2+0] = surfacecopy[i*2+0];
337             surfaceBuffer[i*2+1] = surfacecopy[i*2+1];
338           }
339         }
340       }
341 
342     } else if (me >= nhalf && me < 2*nhalf) {
343       MPI_Send(imageBuffer,npixels*3,MPI_BYTE,me-nhalf,0,world);
344       MPI_Send(depthBuffer,npixels,MPI_DOUBLE,me-nhalf,0,world);
345       if (ssao) MPI_Send(surfaceBuffer,npixels*2,MPI_DOUBLE,me-nhalf,0,world);
346     }
347 
348     nhalf /= 2;
349   }
350 
351   // extra SSAO enhancement
352   // bcast full image to all procs
353   // each works on subset of pixels
354   // MPI_Gather() result back to proc 0
355   // use Gatherv() if subset of pixels is not the same size on every proc
356 
357   if (ssao) {
358     MPI_Bcast(imageBuffer,npixels*3,MPI_BYTE,0,world);
359     MPI_Bcast(surfaceBuffer,npixels*2,MPI_DOUBLE,0,world);
360     MPI_Bcast(depthBuffer,npixels,MPI_DOUBLE,0,world);
361     compute_SSAO();
362 
363     int pixelstart = 3 * static_cast<int> (1.0*me/nprocs * npixels);
364     int pixelstop = 3 * static_cast<int> (1.0*(me+1)/nprocs * npixels);
365     int mypixels = pixelstop - pixelstart;
366 
367     if (npixels % nprocs == 0) {
368       MPI_Gather(imageBuffer+pixelstart,mypixels,MPI_BYTE,
369                  rgbcopy,mypixels,MPI_BYTE,0,world);
370 
371     } else {
372       if (recvcounts == NULL) {
373         memory->create(recvcounts,nprocs,"image:recvcounts");
374         memory->create(displs,nprocs,"image:displs");
375         MPI_Allgather(&mypixels,1,MPI_INT,recvcounts,1,MPI_INT,world);
376         displs[0] = 0;
377         for (int i = 1; i < nprocs; i++)
378           displs[i] = displs[i-1] + recvcounts[i-1];
379       }
380 
381       MPI_Gatherv(imageBuffer+pixelstart,mypixels,MPI_BYTE,
382                   rgbcopy,recvcounts,displs,MPI_BYTE,0,world);
383     }
384 
385     writeBuffer = rgbcopy;
386   } else {
387     writeBuffer = imageBuffer;
388   }
389 }
390 
391 /* ----------------------------------------------------------------------
392    draw a line as a cylinder
393 ------------------------------------------------------------------------- */
394 
draw_line(double * x1,double * x2,double * color,double diameter)395 void Image::draw_line(double *x1, double *x2, double *color, double diameter)
396 {
397   draw_cylinder(x1,x2,color,diameter,3);
398 }
399 
400 /* ----------------------------------------------------------------------
401    draw outline of a 3d box as 12 cylinders
402 ------------------------------------------------------------------------- */
403 
draw_box(double (* corners)[3],double * color,double diameter)404 void Image::draw_box(double (*corners)[3], double *color, double diameter)
405 {
406   draw_cylinder(corners[0],corners[1],color,diameter,3);
407   draw_cylinder(corners[2],corners[3],color,diameter,3);
408   draw_cylinder(corners[0],corners[2],color,diameter,3);
409   draw_cylinder(corners[1],corners[3],color,diameter,3);
410   draw_cylinder(corners[0],corners[4],color,diameter,3);
411   draw_cylinder(corners[1],corners[5],color,diameter,3);
412   draw_cylinder(corners[2],corners[6],color,diameter,3);
413   draw_cylinder(corners[3],corners[7],color,diameter,3);
414   draw_cylinder(corners[4],corners[5],color,diameter,3);
415   draw_cylinder(corners[6],corners[7],color,diameter,3);
416   draw_cylinder(corners[4],corners[6],color,diameter,3);
417   draw_cylinder(corners[5],corners[7],color,diameter,3);
418 }
419 
420 /* ----------------------------------------------------------------------
421    draw outline of a 2d rectangle as 4 cylinders
422 ------------------------------------------------------------------------- */
423 
draw_box2d(double (* corners)[3],double * color,double diameter)424 void Image::draw_box2d(double (*corners)[3], double *color, double diameter)
425 {
426   draw_cylinder(corners[0],corners[1],color,diameter,3);
427   draw_cylinder(corners[2],corners[3],color,diameter,3);
428   draw_cylinder(corners[0],corners[2],color,diameter,3);
429   draw_cylinder(corners[1],corners[3],color,diameter,3);
430 }
431 
432 /* ----------------------------------------------------------------------
433    draw XYZ axes in red/green/blue
434    axes = 4 end points
435 ------------------------------------------------------------------------- */
436 
draw_axes(double (* axes)[3],double diameter)437 void Image::draw_axes(double (*axes)[3], double diameter)
438 {
439   draw_cylinder(axes[0],axes[1],color2rgb("red"),diameter,3);
440   draw_cylinder(axes[0],axes[2],color2rgb("green"),diameter,3);
441   draw_cylinder(axes[0],axes[3],color2rgb("blue"),diameter,3);
442 }
443 
444 /* ----------------------------------------------------------------------
445    draw sphere at x with surfaceColor and diameter
446    render pixel by pixel onto image plane with depth buffering
447 ------------------------------------------------------------------------- */
448 
draw_sphere(double * x,double * surfaceColor,double diameter)449 void Image::draw_sphere(double *x, double *surfaceColor, double diameter)
450 {
451   int ix,iy;
452   double projRad;
453   double xlocal[3],surface[3];
454   double depth;
455 
456   xlocal[0] = x[0] - xctr;
457   xlocal[1] = x[1] - yctr;
458   xlocal[2] = x[2] - zctr;
459 
460   double xmap = MathExtra::dot3(camRight,xlocal);
461   double ymap = MathExtra::dot3(camUp,xlocal);
462   double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(xlocal,camDir);
463 
464   double radius = 0.5*diameter;
465   double radsq = radius*radius;
466   double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
467     -tanPerPixel / zoom;
468   double pixelRadiusFull = radius / pixelWidth;
469   int pixelRadius = static_cast<int> (pixelRadiusFull + 0.5) + 1;
470 
471   double xf = xmap / pixelWidth;
472   double yf = ymap / pixelWidth;
473   int xc = static_cast<int> (xf);
474   int yc = static_cast<int> (yf);
475   double width_error = xf - xc;
476   double height_error = yf - yc;
477 
478   // shift 0,0 to screen center (vs lower left)
479 
480   xc += width / 2;
481   yc += height / 2;
482 
483   for (iy = yc - pixelRadius; iy <= yc + pixelRadius; iy++) {
484     for (ix = xc - pixelRadius; ix <= xc + pixelRadius; ix++) {
485       if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;
486 
487       surface[1] = ((iy - yc) - height_error) * pixelWidth;
488       surface[0] = ((ix - xc) - width_error) * pixelWidth;
489       projRad = surface[0]*surface[0] + surface[1]*surface[1];
490 
491       // outside the sphere in the projected image
492 
493       if (projRad > radsq) continue;
494       surface[2] = sqrt(radsq - projRad);
495       depth = dist - surface[2];
496 
497       surface[0] /= radius;
498       surface[1] /= radius;
499       surface[2] /= radius;
500 
501       draw_pixel (ix, iy, depth, surface, surfaceColor);
502     }
503   }
504 }
505 
506 /* ----------------------------------------------------------------------
507    draw axis oriented brick at x with surfaceColor and diameter[3] in extent
508    render pixel by pixel onto image plane with depth buffering
509 ------------------------------------------------------------------------- */
510 
draw_brick(double * x,double * surfaceColor,double * diameter)511 void Image::draw_brick(double *x, double *surfaceColor, double *diameter)
512 {
513   double xlocal[3],surface[3],normal[3];
514   double t,tdir[3];
515   double depth;
516 
517   xlocal[0] = x[0] - xctr;
518   xlocal[1] = x[1] - yctr;
519   xlocal[2] = x[2] - zctr;
520 
521   double xmap = MathExtra::dot3(camRight,xlocal);
522   double ymap = MathExtra::dot3(camUp,xlocal);
523   double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(xlocal,camDir);
524 
525   double radius[3];
526   radius[0] = 0.5*diameter[0];
527   radius[1] = 0.5*diameter[1];
528   radius[2] = 0.5*diameter[2];
529 
530   double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
531     -tanPerPixel / zoom;
532 
533   double halfWidth = MAX(diameter[0],diameter[1]);
534   halfWidth = MAX(halfWidth,diameter[2]);
535   double pixelHalfWidthFull = halfWidth / pixelWidth;
536   int pixelHalfWidth = static_cast<int> (pixelHalfWidthFull + 0.5);
537 
538   double xf = xmap / pixelWidth;
539   double yf = ymap / pixelWidth;
540   int xc = static_cast<int> (xf);
541   int yc = static_cast<int> (yf);
542   double width_error = xf - xc;
543   double height_error = yf - yc;
544 
545   // shift 0,0 to screen center (vs lower left)
546 
547   xc += width / 2;
548   yc += height / 2;
549 
550   for (int iy = yc - pixelHalfWidth; iy <= yc + pixelHalfWidth; iy ++) {
551     for (int ix = xc - pixelHalfWidth; ix <= xc + pixelHalfWidth; ix ++) {
552       if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;
553 
554       double sy = ((iy - yc) - height_error) * pixelWidth;
555       double sx = ((ix - xc) - width_error) * pixelWidth;
556       surface[0] = camRight[0] * sx + camUp[0] * sy;
557       surface[1] = camRight[1] * sx + camUp[1] * sy;
558       surface[2] = camRight[2] * sx + camUp[2] * sy;
559 
560       // iterate through each of the 6 axis-oriented planes of the box
561       // only render up to 3 which are facing the camera
562       // these checks short circuit a dot product, testing for > 0
563 
564       for (int dim = 0; dim < 3; dim ++) {
565         if (camDir[dim] > 0) {          // positive faces camera
566           t = (radius[dim] - surface[dim]) / camDir[dim];
567           normal[0] = camRight[dim];
568           normal[1] = camUp[dim];
569           normal[2] = camDir[dim];
570         } else if (camDir[dim] < 0) {   // negative faces camera
571           t = -(radius[dim] + surface[dim]) / camDir[dim];
572           normal[0] = -camRight[dim];
573           normal[1] = -camUp[dim];
574           normal[2] = -camDir[dim];
575         }
576         if (camDir[dim] != 0) {
577           tdir[0] = camDir[0] * t;
578           tdir[1] = camDir[1] * t;
579           tdir[2] = camDir[2] * t;
580 
581           bool xin = ((surface[0]+tdir[0]) >= -radius[0]) &&
582 	    ((surface[0]+tdir[0]) <= radius[0]);
583           bool yin = ((surface[1]+tdir[1]) >= -radius[1]) &&
584 	    ((surface[1]+tdir[1]) <= radius[1]);
585           bool zin = ((surface[2]+tdir[2]) >= -radius[2]) &&
586 	    ((surface[2]+tdir[2]) <= radius[2]);
587 
588           switch (dim) {
589 	  case 0:
590 	    if (yin & zin) {
591 	      depth = dist - t;
592 	      draw_pixel (ix, iy, depth, normal, surfaceColor);
593 	    }
594 	    break;
595 	  case 1:
596 	    if (xin & zin) {
597 	      depth = dist - t;
598 	      draw_pixel (ix, iy, depth, normal, surfaceColor);
599 	    }
600 	    break;
601 	  case 2:
602 	    if (xin & yin) {
603 	      depth = dist - t;
604 	      draw_pixel (ix, iy, depth, normal, surfaceColor);
605 	    }
606 	    break;
607           }
608         }
609       }
610     }
611   }
612 }
613 
614 /* ----------------------------------------------------------------------
615    draw cylinder from x to y with surfaceColor and diameter
616    render pixel by pixel onto image plane with depth buffering
617    if sflag = 0, draw no end spheres
618    if sflag = 1, draw 1st end sphere
619    if sflag = 2, draw 2nd end sphere
620    if sflag = 3, draw both end spheres
621 ------------------------------------------------------------------------- */
622 
draw_cylinder(double * x,double * y,double * surfaceColor,double diameter,int sflag)623 void Image::draw_cylinder(double *x, double *y,
624 			  double *surfaceColor, double diameter, int sflag)
625 {
626   double surface[3], normal[3];
627   double mid[3],xaxis[3],yaxis[3],zaxis[3];
628   double camLDir[3], camLRight[3], camLUp[3];
629   double zmin, zmax;
630 
631   if (sflag % 2) draw_sphere(x,surfaceColor,diameter);
632   if (sflag/2) draw_sphere(y,surfaceColor,diameter);
633 
634   double radius = 0.5*diameter;
635   double radsq = radius*radius;
636 
637   zaxis[0] = y[0] - x[0];
638   zaxis[1] = y[1] - x[1];
639   zaxis[2] = y[2] - x[2];
640 
641   double rasterWidth = fabs(MathExtra::dot3(zaxis, camRight)) + diameter;
642   double rasterHeight = fabs(MathExtra::dot3(zaxis, camUp)) + diameter;
643 
644   mid[0] = (y[0] + x[0]) * 0.5 - xctr;
645   mid[1] = (y[1] + x[1]) * 0.5 - yctr;
646   mid[2] = (y[2] + x[2]) * 0.5 - zctr;
647 
648   double len = MathExtra::len3(zaxis);
649   MathExtra::scale3(1.0/len,zaxis);
650   len *= 0.5;
651   zmax = len;
652   zmin = -len;
653 
654   double xmap = MathExtra::dot3(camRight,mid);
655   double ymap = MathExtra::dot3(camUp,mid);
656   double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(mid,camDir);
657 
658   double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
659     -tanPerPixel / zoom;
660 
661   double xf = xmap / pixelWidth;
662   double yf = ymap / pixelWidth;
663   int xc = static_cast<int> (xf);
664   int yc = static_cast<int> (yf);
665   double width_error = xf - xc;
666   double height_error = yf - yc;
667 
668   // shift 0,0 to screen center (vs lower left)
669 
670   xc += width / 2;
671   yc += height / 2;
672 
673   double pixelHalfWidthFull = (rasterWidth * 0.5) / pixelWidth;
674   double pixelHalfHeightFull = (rasterHeight * 0.5) / pixelWidth;
675   int pixelHalfWidth = static_cast<int> (pixelHalfWidthFull + 0.5);
676   int pixelHalfHeight = static_cast<int> (pixelHalfHeightFull + 0.5);
677 
678   if (zaxis[0] == camDir[0] && zaxis[1] == camDir[1] && zaxis[2] == camDir[2])
679     return;
680   if (zaxis[0] == -camDir[0] && zaxis[1] == -camDir[1] &&
681       zaxis[2] == -camDir[2]) return;
682 
683   MathExtra::cross3(zaxis,camDir,yaxis);
684   MathExtra::norm3(yaxis);
685   MathExtra::cross3(yaxis,zaxis,xaxis);
686   MathExtra::norm3(xaxis);
687 
688   camLDir[0] = MathExtra::dot3(camDir,xaxis);
689   camLDir[1] = 0.0;
690   camLDir[2] = MathExtra::dot3(camDir,zaxis);
691 
692   camLRight[0] = MathExtra::dot3(camRight,xaxis);
693   camLRight[1] = MathExtra::dot3(camRight,yaxis);
694   camLRight[2] = MathExtra::dot3(camRight,zaxis);
695   MathExtra::norm3(camLRight);
696 
697   camLUp[0] = MathExtra::dot3(camUp,xaxis);
698   camLUp[1] = MathExtra::dot3(camUp,yaxis);
699   camLUp[2] = MathExtra::dot3(camUp,zaxis);
700   MathExtra::norm3(camLUp);
701 
702   double a = camLDir[0] * camLDir[0];
703 
704   for (int iy = yc - pixelHalfHeight; iy <= yc + pixelHalfHeight; iy ++) {
705     for (int ix = xc - pixelHalfWidth; ix <= xc + pixelHalfWidth; ix ++) {
706       if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;
707 
708       double sy = ((iy - yc) - height_error) * pixelWidth;
709       double sx = ((ix - xc) - width_error) * pixelWidth;
710       surface[0] = camLRight[0] * sx + camLUp[0] * sy;
711       surface[1] = camLRight[1] * sx + camLUp[1] * sy;
712       surface[2] = camLRight[2] * sx + camLUp[2] * sy;
713 
714       double b = 2 * camLDir[0] * surface[0];
715       double c = surface[0] * surface[0] + surface[1] * surface[1] - radsq;
716 
717       double partial = b*b - 4*a*c;
718       if (partial < 0) continue;
719       partial = sqrt (partial);
720 
721       double t = (-b + partial) / (2*a);
722       double t2 = (-b - partial) / (2*a);
723       if (t2 > t) { t = t2; }
724 
725       surface[0] += t * camLDir[0];
726       surface[1] += t * camLDir[1];
727       surface[2] += t * camLDir[2];
728 
729       if (surface[2] > zmax || surface[2] < zmin) continue;
730 
731       // convert surface into the surface normal
732 
733       normal[0] = surface[0] / radius;
734       normal[1] = surface[1] / radius;
735       normal[2] = 0.0;
736 
737       // in camera space
738 
739       surface[0] = MathExtra::dot3 (normal, camLRight);
740       surface[1] = MathExtra::dot3 (normal, camLUp);
741       surface[2] = MathExtra::dot3 (normal, camLDir);
742 
743       double depth = dist - t;
744       draw_pixel (ix, iy, depth, surface, surfaceColor);
745     }
746   }
747 }
748 
749 /* ----------------------------------------------------------------------
750    draw triangle with 3 corner points x,y,z and surfaceColor
751 ------------------------------------------------------------------------- */
752 
draw_triangle(double * x,double * y,double * z,double * surfaceColor)753 void Image::draw_triangle(double *x, double *y, double *z, double *surfaceColor)
754 {
755   double d1[3], d1len, d2[3], d2len, normal[3], invndotd;
756   double xlocal[3], ylocal[3], zlocal[3];
757   double surface[3];
758   double depth;
759 
760   xlocal[0] = x[0] - xctr;
761   xlocal[1] = x[1] - yctr;
762   xlocal[2] = x[2] - zctr;
763   ylocal[0] = y[0] - xctr;
764   ylocal[1] = y[1] - yctr;
765   ylocal[2] = y[2] - zctr;
766   zlocal[0] = z[0] - xctr;
767   zlocal[1] = z[1] - yctr;
768   zlocal[2] = z[2] - zctr;
769 
770   MathExtra::sub3 (xlocal, ylocal, d1);
771   d1len = MathExtra::len3 (d1);
772   MathExtra::scale3 (1.0 / d1len, d1);
773   MathExtra::sub3 (zlocal, ylocal, d2);
774   d2len = MathExtra::len3 (d2);
775   MathExtra::scale3 (1.0 / d2len, d2);
776 
777   MathExtra::cross3 (d1, d2, normal);
778   MathExtra::norm3 (normal);
779 
780   // ----------------
781   // old code
782 
783   //invndotd = 1.0 / MathExtra::dot3(normal, camDir);
784 
785   // invalid triangle (parallel)
786 
787   //if (invndotd == 0) return;
788 
789   // ----------------
790   // new code
791 
792   double ndotd = MathExtra::dot3(normal,camDir);
793   if (ndotd >= 0.0) return;
794   invndotd = 1.0 / ndotd;
795 
796   // ----------------
797 
798   double r[3],u[3];
799 
800   r[0] = MathExtra::dot3(camRight,xlocal);
801   r[1] = MathExtra::dot3(camRight,ylocal);
802   r[2] = MathExtra::dot3(camRight,zlocal);
803 
804   u[0] = MathExtra::dot3(camUp,xlocal);
805   u[1] = MathExtra::dot3(camUp,ylocal);
806   u[2] = MathExtra::dot3(camUp,zlocal);
807 
808   double rasterLeft = r[0] - MIN(r[0],MIN(r[1],r[2]));
809   double rasterRight = MAX(r[0],MAX(r[1],r[2])) - r[0];
810   double rasterDown = u[0] - MIN(u[0],MIN(u[1],u[2]));
811   double rasterUp = MAX(u[0],MAX(u[1],u[2])) - u[0];
812 
813   double xmap = MathExtra::dot3(camRight,xlocal);
814   double ymap = MathExtra::dot3(camUp,xlocal);
815   double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(xlocal,camDir);
816 
817   double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
818     -tanPerPixel / zoom;
819 
820   double xf = xmap / pixelWidth;
821   double yf = ymap / pixelWidth;
822   int xc = static_cast<int> (xf);
823   int yc = static_cast<int> (yf);
824   double width_error = xf - xc;
825   double height_error = yf - yc;
826 
827   // shift 0,0 to screen center (vs lower left)
828 
829   xc += width / 2;
830   yc += height / 2;
831 
832   double pixelLeftFull = rasterLeft / pixelWidth;
833   double pixelRightFull = rasterRight / pixelWidth;
834   double pixelDownFull = rasterDown / pixelWidth;
835   double pixelUpFull = rasterUp / pixelWidth;
836 
837   //old
838   //int pixelLeft = static_cast<int> (pixelLeftFull + 0.5);
839   //int pixelRight = static_cast<int> (pixelRightFull + 0.5);
840   //int pixelDown = static_cast<int> (pixelDownFull + 0.5);
841   //int pixelUp = static_cast<int> (pixelUpFull + 0.5);
842 
843   // new
844   int pixelLeft = static_cast<int> (pixelLeftFull + 1.0);
845   int pixelRight = static_cast<int> (pixelRightFull + 1.0);
846   int pixelDown = static_cast<int> (pixelDownFull + 1.0);
847   int pixelUp = static_cast<int> (pixelUpFull + 1.0);
848 
849   for (int iy = yc - pixelDown; iy <= yc + pixelUp; iy ++) {
850     for (int ix = xc - pixelLeft; ix <= xc + pixelRight; ix ++) {
851       if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;
852 
853       double sy = ((iy - yc) - height_error) * pixelWidth;
854       double sx = ((ix - xc) - width_error) * pixelWidth;
855       surface[0] = camRight[0] * sx + camUp[0] * sy;
856       surface[1] = camRight[1] * sx + camUp[1] * sy;
857       surface[2] = camRight[2] * sx + camUp[2] * sy;
858 
859       double t = -MathExtra::dot3(normal,surface) * invndotd;
860 
861       // test inside the triangle
862 
863       double p[3];
864       p[0] = xlocal[0] + surface[0] + camDir[0] * t;
865       p[1] = xlocal[1] + surface[1] + camDir[1] * t;
866       p[2] = xlocal[2] + surface[2] + camDir[2] * t;
867 
868       double s1[3], s2[3], s3[3];
869       double c1[3], c2[3];
870 
871       MathExtra::sub3 (zlocal, xlocal, s1);
872       MathExtra::sub3 (ylocal, xlocal, s2);
873       MathExtra::sub3 (p, xlocal, s3);
874       MathExtra::cross3 (s1, s2, c1);
875       MathExtra::cross3 (s1, s3, c2);
876       if (MathExtra::dot3 (c1, c2) <= 0) continue;
877 
878       MathExtra::sub3 (xlocal, ylocal, s1);
879       MathExtra::sub3 (zlocal, ylocal, s2);
880       MathExtra::sub3 (p, ylocal, s3);
881       MathExtra::cross3 (s1, s2, c1);
882       MathExtra::cross3 (s1, s3, c2);
883       if (MathExtra::dot3 (c1, c2) <= 0) continue;
884 
885       MathExtra::sub3 (ylocal, zlocal, s1);
886       MathExtra::sub3 (xlocal, zlocal, s2);
887       MathExtra::sub3 (p, zlocal, s3);
888       MathExtra::cross3 (s1, s2, c1);
889       MathExtra::cross3 (s1, s3, c2);
890       if (MathExtra::dot3 (c1, c2) <= 0) continue;
891 
892       double cNormal[3];
893       cNormal[0] = MathExtra::dot3(camRight, normal);
894       cNormal[1] = MathExtra::dot3(camUp, normal);
895       cNormal[2] = MathExtra::dot3(camDir, normal);
896 
897       depth = dist - t;
898       draw_pixel(ix,iy,depth,cNormal,surfaceColor);
899     }
900   }
901 }
902 
903 /* ---------------------------------------------------------------------- */
904 
draw_pixel(int ix,int iy,double depth,double * surface,double * surfaceColor)905 void Image::draw_pixel(int ix, int iy, double depth,
906 			   double *surface, double *surfaceColor)
907 {
908   double diffuseKey,diffuseFill,diffuseBack,specularKey;
909   if (depth < 0 || (depthBuffer[ix + iy*width] >= 0 &&
910 		    depth >= depthBuffer[ix + iy*width])) return;
911   depthBuffer[ix + iy*width] = depth;
912 
913   // store only the tangent relative to the camera normal (0,0,-1)
914 
915   surfaceBuffer[0 + ix * 2 + iy*width * 2] = surface[1];
916   surfaceBuffer[1 + ix * 2 + iy*width * 2] = -surface[0];
917 
918   diffuseKey = saturate(MathExtra::dot3(surface, keyLightDir));
919   diffuseFill = saturate(MathExtra::dot3(surface, fillLightDir));
920   diffuseBack = saturate(MathExtra::dot3(surface, backLightDir));
921   specularKey = pow(saturate(MathExtra::dot3(surface, keyHalfDir)),
922 		    specularHardness) * specularIntensity;
923 
924   double c[3];
925   c[0] = surfaceColor[0] * ambientColor[0];
926   c[1] = surfaceColor[1] * ambientColor[1];
927   c[2] = surfaceColor[2] * ambientColor[2];
928 
929   c[0] += surfaceColor[0] * keyLightColor[0] * diffuseKey;
930   c[1] += surfaceColor[1] * keyLightColor[1] * diffuseKey;
931   c[2] += surfaceColor[2] * keyLightColor[2] * diffuseKey;
932 
933   c[0] += keyLightColor[0] * specularKey;
934   c[1] += keyLightColor[1] * specularKey;
935   c[2] += keyLightColor[2] * specularKey;
936 
937   c[0] += surfaceColor[0] * fillLightColor[0] * diffuseFill;
938   c[1] += surfaceColor[1] * fillLightColor[1] * diffuseFill;
939   c[2] += surfaceColor[2] * fillLightColor[2] * diffuseFill;
940 
941   c[0] += surfaceColor[0] * backLightColor[0] * diffuseBack;
942   c[1] += surfaceColor[1] * backLightColor[1] * diffuseBack;
943   c[2] += surfaceColor[2] * backLightColor[2] * diffuseBack;
944 
945   c[0] = saturate(c[0]);
946   c[1] = saturate(c[1]);
947   c[2] = saturate(c[2]);
948 
949   imageBuffer[0 + ix*3 + iy*width*3] = static_cast<int>(c[0] * 255.0);
950   imageBuffer[1 + ix*3 + iy*width*3] = static_cast<int>(c[1] * 255.0);
951   imageBuffer[2 + ix*3 + iy*width*3] = static_cast<int>(c[2] * 255.0);
952 }
953 
954 /* ---------------------------------------------------------------------- */
955 
compute_SSAO()956 void Image::compute_SSAO()
957 {
958   // used for rasterizing the spheres
959 
960   double delTheta = 2.0*MY_PI / SSAOSamples;
961 
962   // typical neighborhood value for shading
963 
964   double pixelWidth = (tanPerPixel > 0) ? tanPerPixel :
965 	-tanPerPixel / zoom;
966   int pixelRadius = (int) trunc (SSAORadius / pixelWidth + 0.5);
967 
968   // each proc is assigned a subset of contiguous pixels from the full image
969   // pixels are contiguous in x (columns within a row), then by row
970   // index = pixels from 0 to npixel-1
971   // x = column # from 0 to width-1
972   // y = row # from 0 to height-1
973 
974   int pixelstart = static_cast<int> (1.0*me/nprocs * npixels);
975   int pixelstop = static_cast<int> (1.0*(me+1)/nprocs * npixels);
976 
977   for (int index = pixelstart; index < pixelstop; index++) {
978     int x = index % width;
979     int y = index / width;
980 
981     double cdepth = depthBuffer[index];
982     if (cdepth < 0) { continue; }
983 
984     double sx = surfaceBuffer[index * 2 + 0];
985     double sy = surfaceBuffer[index * 2 + 1];
986     double sin_t = -sqrt(sx*sx + sy*sy);
987 
988     // DEBUG - remove randomness so image is same on any proc count
989     //double mytheta = 0.5 * SSAOJitter;
990     double mytheta = random->uniform() * SSAOJitter;
991     double ao = 0.0;
992 
993     for (int s = 0; s < SSAOSamples; s ++) {
994       double hx = cos(mytheta);
995       double hy = sin(mytheta);
996       mytheta += delTheta;
997 
998       // multiply by z cross surface tangent
999       // so that dot (aka cos) works here
1000 
1001       double scaled_sin_t = sin_t * (hx*sy + hy*sx);
1002 
1003       // Bresenham's line algorithm to march over depthBuffer
1004 
1005       int dx = static_cast<int> (hx * pixelRadius);
1006       int dy = static_cast<int> (hy * pixelRadius);
1007       int ex = x + dx;
1008       if (ex < 0) { ex = 0; } if (ex >= width) { ex = width - 1; }
1009       int ey = y + dy;
1010       if (ey < 0) { ey = 0; } if (ey >= height) { ey = height - 1; }
1011       double delta;
1012       int small, large;
1013       double lenIncr;
1014       if (fabs(hx) > fabs(hy)) {
1015         small = (hx > 0) ? 1 : -1;
1016         large = (hy > 0) ? width : -width;
1017         delta = fabs(hy / hx);
1018       } else {
1019         small = (hy > 0) ? width : -width;
1020         large = (hx > 0) ? 1 : -1;
1021         delta = fabs(hx / hy);
1022       }
1023       lenIncr = sqrt (1 + delta * delta) * pixelWidth;
1024 
1025       // initialize with one step
1026       // because the center point doesn't need testing
1027 
1028       int end = ex + ey * width;
1029       int ind = index + small;
1030       double len = lenIncr;
1031       double err = delta;
1032       if (err >= 1.0) {
1033         ind += large;
1034         err -= 1.0;
1035       }
1036 
1037       double minPeak = -1;
1038       double peakLen = 0.0;
1039       int stepsTaken = 1;
1040       while ((small > 0 && ind <= end) || (small < 0 && ind >= end)) {
1041         if (ind < 0 || ind >= (width*height)) {
1042           break;
1043         }
1044 
1045         // cdepth - depthBuffer B/C we want it in the negative z direction
1046 
1047         if (minPeak < 0 || (depthBuffer[ind] >= 0 &&
1048                             depthBuffer[ind] < minPeak)) {
1049           minPeak = depthBuffer[ind];
1050           peakLen = len;
1051         }
1052         ind += small;
1053         len += lenIncr;
1054         err += delta;
1055         if (err >= 1.0) {
1056           ind += large;
1057           err -= 1.0;
1058         }
1059         stepsTaken ++;
1060       }
1061 
1062       if (peakLen > 0) {
1063         double h = atan ((cdepth - minPeak) / peakLen);
1064         ao += saturate(sin (h) - scaled_sin_t);
1065       } else {
1066         ao += saturate(-scaled_sin_t);
1067       }
1068     }
1069     ao /= (double)SSAOSamples;
1070 
1071     double c[3];
1072     c[0] = (double) (*(unsigned char *) &imageBuffer[index * 3 + 0]);
1073     c[1] = (double) (*(unsigned char *) &imageBuffer[index * 3 + 1]);
1074     c[2] = (double) (*(unsigned char *) &imageBuffer[index * 3 + 2]);
1075     c[0] *= (1.0 - ao);
1076     c[1] *= (1.0 - ao);
1077     c[2] *= (1.0 - ao);
1078     imageBuffer[index * 3 + 0] = (int) c[0];
1079     imageBuffer[index * 3 + 1] = (int) c[1];
1080     imageBuffer[index * 3 + 2] = (int) c[2];
1081   }
1082 }
1083 
1084 /* ---------------------------------------------------------------------- */
1085 
1086 #ifdef SPARTA_JPEG
write_JPG(FILE * fp)1087 void Image::write_JPG(FILE *fp)
1088 {
1089   struct jpeg_compress_struct cinfo;
1090   struct jpeg_error_mgr jerr;
1091   JSAMPROW row_pointer;
1092 
1093   cinfo.err = jpeg_std_error(&jerr);
1094   jpeg_create_compress(&cinfo);
1095   jpeg_stdio_dest(&cinfo,fp);
1096   cinfo.image_width = width;
1097   cinfo.image_height = height;
1098   cinfo.input_components = 3;
1099   cinfo.in_color_space = JCS_RGB;
1100 
1101   jpeg_set_defaults(&cinfo);
1102   jpeg_set_quality(&cinfo,85,TRUE);
1103   jpeg_start_compress(&cinfo,TRUE);
1104 
1105   while (cinfo.next_scanline < cinfo.image_height) {
1106     row_pointer = (JSAMPROW)
1107       &writeBuffer[(cinfo.image_height - 1 - cinfo.next_scanline) * 3 * width];
1108     jpeg_write_scanlines(&cinfo,&row_pointer,1);
1109   }
1110 
1111   jpeg_finish_compress(&cinfo);
1112   jpeg_destroy_compress(&cinfo);
1113 }
1114 #else
write_JPG(FILE *)1115 void Image::write_JPG(FILE*)
1116 {
1117 }
1118 #endif
1119 
1120 /* ---------------------------------------------------------------------- */
1121 
1122 #ifdef SPARTA_PNG
write_PNG(FILE * fp)1123 void Image::write_PNG(FILE *fp)
1124 {
1125   png_structp png_ptr;
1126   png_infop info_ptr;
1127 
1128   png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1129 
1130   if (!png_ptr) return;
1131 
1132   info_ptr = png_create_info_struct(png_ptr);
1133   if (!info_ptr) {
1134     png_destroy_write_struct(&png_ptr, NULL);
1135     return;
1136   }
1137 
1138   if (setjmp(png_jmpbuf(png_ptr))) {
1139     png_destroy_write_struct(&png_ptr, &info_ptr);
1140     return;
1141   }
1142 
1143   png_init_io(png_ptr, fp);
1144   png_set_compression_level(png_ptr,Z_BEST_COMPRESSION);
1145   png_set_IHDR(png_ptr,info_ptr,width,height,8,PNG_COLOR_TYPE_RGB,
1146     PNG_INTERLACE_NONE,PNG_COMPRESSION_TYPE_DEFAULT,PNG_FILTER_TYPE_DEFAULT);
1147 
1148   png_text text_ptr[2];
1149   memset(text_ptr,0,2*sizeof(png_text));
1150 
1151   char key0[]  = "Software";
1152   char text0[] = "SPARTA " SPARTA_VERSION;
1153   char key1[]  = "Description";
1154   char text1[] = "Dump image snapshot";
1155   text_ptr[0].key = key0;
1156   text_ptr[0].text = text0;
1157   text_ptr[1].key = key1;
1158   text_ptr[1].text = text1;
1159   text_ptr[0].compression = PNG_TEXT_COMPRESSION_NONE;
1160   text_ptr[1].compression = PNG_TEXT_COMPRESSION_NONE;
1161 
1162   png_set_text(png_ptr,info_ptr,text_ptr,1);
1163   png_write_info(png_ptr,info_ptr);
1164 
1165   png_bytep row_pointers[height];
1166   for (int i=0; i < height; ++i)
1167     row_pointers[i] = (png_bytep) &writeBuffer[(height-i-1)*3*width];
1168 
1169   png_write_image(png_ptr, row_pointers);
1170   png_write_end(png_ptr, info_ptr);
1171 
1172   png_destroy_write_struct(&png_ptr, &info_ptr);
1173 }
1174 #else
write_PNG(FILE *)1175 void Image::write_PNG(FILE *)
1176 {
1177 }
1178 #endif
1179 
1180 /* ---------------------------------------------------------------------- */
1181 
write_PPM(FILE * fp)1182 void Image::write_PPM(FILE *fp)
1183 {
1184   fprintf(fp,"P6\n%d %d\n255\n",width,height);
1185 
1186   int y;
1187   for (y = height-1; y >= 0; y--)
1188     fwrite(&writeBuffer[y*width*3],3,width,fp);
1189 }
1190 
1191 /* ----------------------------------------------------------------------
1192    return static/dynamic status of color map index
1193 ------------------------------------------------------------------------- */
1194 
map_dynamic(int index)1195 int Image::map_dynamic(int index)
1196 {
1197   return maps[index]->dynamic;
1198 }
1199 
1200 /* ----------------------------------------------------------------------
1201    redefine properties of the color map index
1202    return 1 if any error in args, else return 0
1203 ------------------------------------------------------------------------- */
1204 
map_reset(int index,int narg,char ** arg)1205 int Image::map_reset(int index, int narg, char **arg)
1206 {
1207   return maps[index]->reset(narg,arg);
1208 }
1209 
1210 /* ----------------------------------------------------------------------
1211    set min/max bounds of dynamic color map index
1212 ------------------------------------------------------------------------- */
1213 
map_minmax(int index,double mindynamic,double maxdynamic)1214 int Image::map_minmax(int index, double mindynamic, double maxdynamic)
1215 {
1216   return maps[index]->minmax(mindynamic,maxdynamic);
1217 }
1218 
1219 /* ----------------------------------------------------------------------
1220    return 3-vector color corresponding to value from color map index
1221 ------------------------------------------------------------------------- */
1222 
map_value2color(int index,double value)1223 double *Image::map_value2color(int index, double value)
1224 {
1225   return maps[index]->value2color(value);
1226 }
1227 
1228 /* ----------------------------------------------------------------------
1229    add a new color to username and userrgb
1230    redefine RGB values in userrgb if name already exists
1231    return 1 if RGB values are invalid, else return 0
1232 ------------------------------------------------------------------------- */
1233 
addcolor(char * name,double r,double g,double b)1234 int Image::addcolor(char *name, double r, double g, double b)
1235 {
1236   int icolor;
1237   for (icolor = 0; icolor < ncolors; icolor++)
1238     if (strcmp(name,username[icolor]) == 0) break;
1239 
1240   if (icolor == ncolors) {
1241     username = (char **)
1242       memory->srealloc(username,(ncolors+1)*sizeof(char *),"image:username");
1243     memory->grow(userrgb,ncolors+1,3,"image:userrgb");
1244     ncolors++;
1245   }
1246 
1247   int n = strlen(name) + 1;
1248   username[icolor] = new char[n];
1249   strcpy(username[icolor],name);
1250 
1251   if (r < 0.0 || r > 1.0 || g < 0.0 || g > 1.0 || b < 0.0 || b > 1.0)
1252     return 1;
1253 
1254   userrgb[icolor][0] = r;
1255   userrgb[icolor][1] = g;
1256   userrgb[icolor][2] = b;
1257 
1258   return 0;
1259 }
1260 
1261 /* ----------------------------------------------------------------------
1262    if index > 0, return ptr to index-1 color from rgb
1263    if index < 0, return ptr to -index-1 color from userrgb
1264    if index = 0, search the 2 lists of color names for the string color
1265    search user-defined color names first, then the list of NCOLORS names
1266    return a pointer to the 3 floating point RGB values or NULL if didn't find
1267 ------------------------------------------------------------------------- */
1268 
color2rgb(const char * color,int index)1269 double *Image::color2rgb(const char *color, int index)
1270 {
1271   static const char *name[NCOLORS] = {
1272     "aliceblue",
1273     "antiquewhite",
1274     "aqua",
1275     "aquamarine",
1276     "azure",
1277     "beige",
1278     "bisque",
1279     "black",
1280     "blanchedalmond",
1281     "blue",
1282     "blueviolet",
1283     "brown",
1284     "burlywood",
1285     "cadetblue",
1286     "chartreuse",
1287     "chocolate",
1288     "coral",
1289     "cornflowerblue",
1290     "cornsilk",
1291     "crimson",
1292     "cyan",
1293     "darkblue",
1294     "darkcyan",
1295     "darkgoldenrod",
1296     "darkgray",
1297     "darkgreen",
1298     "darkkhaki",
1299     "darkmagenta",
1300     "darkolivegreen",
1301     "darkorange",
1302     "darkorchid",
1303     "darkred",
1304     "darksalmon",
1305     "darkseagreen",
1306     "darkslateblue",
1307     "darkslategray",
1308     "darkturquoise",
1309     "darkviolet",
1310     "deeppink",
1311     "deepskyblue",
1312     "dimgray",
1313     "dodgerblue",
1314     "firebrick",
1315     "floralwhite",
1316     "forestgreen",
1317     "fuchsia",
1318     "gainsboro",
1319     "ghostwhite",
1320     "gold",
1321     "goldenrod",
1322     "gray",
1323     "green",
1324     "greenyellow",
1325     "honeydew",
1326     "hotpink",
1327     "indianred",
1328     "indigo",
1329     "ivory",
1330     "khaki",
1331     "lavender",
1332     "lavenderblush",
1333     "lawngreen",
1334     "lemonchiffon",
1335     "lightblue",
1336     "lightcoral",
1337     "lightcyan",
1338     "lightgoldenrodyellow",
1339     "lightgreen",
1340     "lightgrey",
1341     "lightpink",
1342     "lightsalmon",
1343     "lightseagreen",
1344     "lightskyblue",
1345     "lightslategray",
1346     "lightsteelblue",
1347     "lightyellow",
1348     "lime",
1349     "limegreen",
1350     "linen",
1351     "magenta",
1352     "maroon",
1353     "mediumaquamarine",
1354     "mediumblue",
1355     "mediumorchid",
1356     "mediumpurple",
1357     "mediumseagreen",
1358     "mediumslateblue",
1359     "mediumspringgreen",
1360     "mediumturquoise",
1361     "mediumvioletred",
1362     "midnightblue",
1363     "mintcream",
1364     "mistyrose",
1365     "moccasin",
1366     "navajowhite",
1367     "navy",
1368     "oldlace",
1369     "olive",
1370     "olivedrab",
1371     "orange",
1372     "orangered",
1373     "orchid",
1374     "palegoldenrod",
1375     "palegreen",
1376     "paleturquoise",
1377     "palevioletred",
1378     "papayawhip",
1379     "peachpuff",
1380     "peru",
1381     "pink",
1382     "plum",
1383     "powderblue",
1384     "purple",
1385     "red",
1386     "rosybrown",
1387     "royalblue",
1388     "saddlebrown",
1389     "salmon",
1390     "sandybrown",
1391     "seagreen",
1392     "seashell",
1393     "sienna",
1394     "silver",
1395     "skyblue",
1396     "slateblue",
1397     "slategray",
1398     "snow",
1399     "springgreen",
1400     "steelblue",
1401     "tan",
1402     "teal",
1403     "thistle",
1404     "tomato",
1405     "turquoise",
1406     "violet",
1407     "wheat",
1408     "white",
1409     "whitesmoke",
1410     "yellow",
1411     "yellowgreen"
1412   };
1413 
1414   static double rgb[NCOLORS][3] = {
1415     {240/255.0, 248/255.0, 255/255.0},
1416     {250/255.0, 235/255.0, 215/255.0},
1417     {0/255.0, 255/255.0, 255/255.0},
1418     {127/255.0, 255/255.0, 212/255.0},
1419     {240/255.0, 255/255.0, 255/255.0},
1420     {245/255.0, 245/255.0, 220/255.0},
1421     {255/255.0, 228/255.0, 196/255.0},
1422     {0/255.0, 0/255.0, 0/255.0},
1423     {255/255.0, 255/255.0, 205/255.0},
1424     {0/255.0, 0/255.0, 255/255.0},
1425     {138/255.0, 43/255.0, 226/255.0},
1426     {165/255.0, 42/255.0, 42/255.0},
1427     {222/255.0, 184/255.0, 135/255.0},
1428     {95/255.0, 158/255.0, 160/255.0},
1429     {127/255.0, 255/255.0, 0/255.0},
1430     {210/255.0, 105/255.0, 30/255.0},
1431     {255/255.0, 127/255.0, 80/255.0},
1432     {100/255.0, 149/255.0, 237/255.0},
1433     {255/255.0, 248/255.0, 220/255.0},
1434     {220/255.0, 20/255.0, 60/255.0},
1435     {0/255.0, 255/255.0, 255/255.0},
1436     {0/255.0, 0/255.0, 139/255.0},
1437     {0/255.0, 139/255.0, 139/255.0},
1438     {184/255.0, 134/255.0, 11/255.0},
1439     {169/255.0, 169/255.0, 169/255.0},
1440     {0/255.0, 100/255.0, 0/255.0},
1441     {189/255.0, 183/255.0, 107/255.0},
1442     {139/255.0, 0/255.0, 139/255.0},
1443     {85/255.0, 107/255.0, 47/255.0},
1444     {255/255.0, 140/255.0, 0/255.0},
1445     {153/255.0, 50/255.0, 204/255.0},
1446     {139/255.0, 0/255.0, 0/255.0},
1447     {233/255.0, 150/255.0, 122/255.0},
1448     {143/255.0, 188/255.0, 143/255.0},
1449     {72/255.0, 61/255.0, 139/255.0},
1450     {47/255.0, 79/255.0, 79/255.0},
1451     {0/255.0, 206/255.0, 209/255.0},
1452     {148/255.0, 0/255.0, 211/255.0},
1453     {255/255.0, 20/255.0, 147/255.0},
1454     {0/255.0, 191/255.0, 255/255.0},
1455     {105/255.0, 105/255.0, 105/255.0},
1456     {30/255.0, 144/255.0, 255/255.0},
1457     {178/255.0, 34/255.0, 34/255.0},
1458     {255/255.0, 250/255.0, 240/255.0},
1459     {34/255.0, 139/255.0, 34/255.0},
1460     {255/255.0, 0/255.0, 255/255.0},
1461     {220/255.0, 220/255.0, 220/255.0},
1462     {248/255.0, 248/255.0, 255/255.0},
1463     {255/255.0, 215/255.0, 0/255.0},
1464     {218/255.0, 165/255.0, 32/255.0},
1465     {128/255.0, 128/255.0, 128/255.0},
1466     {0/255.0, 128/255.0, 0/255.0},
1467     {173/255.0, 255/255.0, 47/255.0},
1468     {240/255.0, 255/255.0, 240/255.0},
1469     {255/255.0, 105/255.0, 180/255.0},
1470     {205/255.0, 92/255.0, 92/255.0},
1471     {75/255.0, 0/255.0, 130/255.0},
1472     {255/255.0, 240/255.0, 240/255.0},
1473     {240/255.0, 230/255.0, 140/255.0},
1474     {230/255.0, 230/255.0, 250/255.0},
1475     {255/255.0, 240/255.0, 245/255.0},
1476     {124/255.0, 252/255.0, 0/255.0},
1477     {255/255.0, 250/255.0, 205/255.0},
1478     {173/255.0, 216/255.0, 230/255.0},
1479     {240/255.0, 128/255.0, 128/255.0},
1480     {224/255.0, 255/255.0, 255/255.0},
1481     {250/255.0, 250/255.0, 210/255.0},
1482     {144/255.0, 238/255.0, 144/255.0},
1483     {211/255.0, 211/255.0, 211/255.0},
1484     {255/255.0, 182/255.0, 193/255.0},
1485     {255/255.0, 160/255.0, 122/255.0},
1486     {32/255.0, 178/255.0, 170/255.0},
1487     {135/255.0, 206/255.0, 250/255.0},
1488     {119/255.0, 136/255.0, 153/255.0},
1489     {176/255.0, 196/255.0, 222/255.0},
1490     {255/255.0, 255/255.0, 224/255.0},
1491     {0/255.0, 255/255.0, 0/255.0},
1492     {50/255.0, 205/255.0, 50/255.0},
1493     {250/255.0, 240/255.0, 230/255.0},
1494     {255/255.0, 0/255.0, 255/255.0},
1495     {128/255.0, 0/255.0, 0/255.0},
1496     {102/255.0, 205/255.0, 170/255.0},
1497     {0/255.0, 0/255.0, 205/255.0},
1498     {186/255.0, 85/255.0, 211/255.0},
1499     {147/255.0, 112/255.0, 219/255.0},
1500     {60/255.0, 179/255.0, 113/255.0},
1501     {123/255.0, 104/255.0, 238/255.0},
1502     {0/255.0, 250/255.0, 154/255.0},
1503     {72/255.0, 209/255.0, 204/255.0},
1504     {199/255.0, 21/255.0, 133/255.0},
1505     {25/255.0, 25/255.0, 112/255.0},
1506     {245/255.0, 255/255.0, 250/255.0},
1507     {255/255.0, 228/255.0, 225/255.0},
1508     {255/255.0, 228/255.0, 181/255.0},
1509     {255/255.0, 222/255.0, 173/255.0},
1510     {0/255.0, 0/255.0, 128/255.0},
1511     {253/255.0, 245/255.0, 230/255.0},
1512     {128/255.0, 128/255.0, 0/255.0},
1513     {107/255.0, 142/255.0, 35/255.0},
1514     {255/255.0, 165/255.0, 0/255.0},
1515     {255/255.0, 69/255.0, 0/255.0},
1516     {218/255.0, 112/255.0, 214/255.0},
1517     {238/255.0, 232/255.0, 170/255.0},
1518     {152/255.0, 251/255.0, 152/255.0},
1519     {175/255.0, 238/255.0, 238/255.0},
1520     {219/255.0, 112/255.0, 147/255.0},
1521     {255/255.0, 239/255.0, 213/255.0},
1522     {255/255.0, 239/255.0, 213/255.0},
1523     {205/255.0, 133/255.0, 63/255.0},
1524     {255/255.0, 192/255.0, 203/255.0},
1525     {221/255.0, 160/255.0, 221/255.0},
1526     {176/255.0, 224/255.0, 230/255.0},
1527     {128/255.0, 0/255.0, 128/255.0},
1528     {255/255.0, 0/255.0, 0/255.0},
1529     {188/255.0, 143/255.0, 143/255.0},
1530     {65/255.0, 105/255.0, 225/255.0},
1531     {139/255.0, 69/255.0, 19/255.0},
1532     {250/255.0, 128/255.0, 114/255.0},
1533     {244/255.0, 164/255.0, 96/255.0},
1534     {46/255.0, 139/255.0, 87/255.0},
1535     {255/255.0, 245/255.0, 238/255.0},
1536     {160/255.0, 82/255.0, 45/255.0},
1537     {192/255.0, 192/255.0, 192/255.0},
1538     {135/255.0, 206/255.0, 235/255.0},
1539     {106/255.0, 90/255.0, 205/255.0},
1540     {112/255.0, 128/255.0, 144/255.0},
1541     {255/255.0, 250/255.0, 250/255.0},
1542     {0/255.0, 255/255.0, 127/255.0},
1543     {70/255.0, 130/255.0, 180/255.0},
1544     {210/255.0, 180/255.0, 140/255.0},
1545     {0/255.0, 128/255.0, 128/255.0},
1546     {216/255.0, 191/255.0, 216/255.0},
1547     {253/255.0, 99/255.0, 71/255.0},
1548     {64/255.0, 224/255.0, 208/255.0},
1549     {238/255.0, 130/255.0, 238/255.0},
1550     {245/255.0, 222/255.0, 179/255.0},
1551     {255/255.0, 255/255.0, 255/255.0},
1552     {245/255.0, 245/255.0, 245/255.0},
1553     {255/255.0, 255/255.0, 0/255.0},
1554     {154/255.0, 205/255.0, 50/255.0}
1555   };
1556 
1557   if (index > 0) {
1558     if (index > NCOLORS) return NULL;
1559     return rgb[index-1];
1560   }
1561   if (index < 0) {
1562     if (-index > ncolors) return NULL;
1563     return userrgb[-index-1];
1564   }
1565 
1566   for (int i = 0; i < ncolors; i++)
1567     if (strcmp(color,username[i]) == 0) return userrgb[i];
1568   for (int i = 0; i < NCOLORS; i++)
1569     if (strcmp(color,name[i]) == 0) return rgb[i];
1570   return NULL;
1571 }
1572 
1573 /* ----------------------------------------------------------------------
1574    return number of default colors
1575 ------------------------------------------------------------------------- */
1576 
default_colors()1577 int Image::default_colors()
1578 {
1579   return NCOLORS;
1580 }
1581 
1582 /* ----------------------------------------------------------------------
1583    search the list of element names for the string element
1584    return a pointer to the 3 floating point RGB values
1585    this list is used by AtomEye and is taken from its Mendeleyev.c file
1586 ------------------------------------------------------------------------- */
1587 
element2color(char * element)1588 double *Image::element2color(char *element)
1589 {
1590   static const char *name[NELEMENTS] = {
1591     "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne",
1592     "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",
1593     "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
1594     "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr",
1595     "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
1596     "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
1597     "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
1598     "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
1599     "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
1600     "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
1601     "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt"
1602   };
1603 
1604   static double rgb[NELEMENTS][3] = {
1605     {0.8, 0.8, 0.8},
1606     {0.6431372549, 0.6666666667, 0.6784313725},
1607     {0.7, 0.7, 0.7},
1608     {0.6431372549, 0.6666666667, 0.6784313725},
1609     {0.9, 0.4, 0},
1610     {0.35, 0.35, 0.35},
1611     {0.2, 0.2, 0.8},
1612     {0.8, 0.2, 0.2},
1613     {0.7, 0.85, 0.45},
1614     {0.6431372549, 0.6666666667, 0.6784313725},
1615     {0.6, 0.6, 0.6},
1616     {0.6, 0.6, 0.7},
1617     {0.6431372549, 0.6666666667, 0.6784313725},
1618     {0.6901960784, 0.768627451, 0.8705882353},
1619     {0.1, 0.7, 0.3},
1620     {0.95, 0.9, 0.2},
1621     {0.15, 0.5, 0.1},
1622     {0.6431372549, 0.6666666667, 0.6784313725},
1623     {0.5, 0.5, 0.5},
1624     {0.8, 0.8, 0.7},
1625     {0.6431372549, 0.6666666667, 0.6784313725},
1626     {0.6431372549, 0.6666666667, 0.6784313725},
1627     {0.6431372549, 0.6666666667, 0.6784313725},
1628     {0, 0.8, 0},
1629     {0.6431372549, 0.6666666667, 0.6784313725},
1630     {0.5176470588, 0.5764705882, 0.6529411765},
1631     {0.6431372549, 0.6666666667, 0.6784313725},
1632     {0.257254902, 0.2666666667, 0.271372549},
1633     {0.95, 0.7900735294, 0.01385869565},
1634     {0.6431372549, 0.6666666667, 0.6784313725},
1635     {0.9, 0, 1},
1636     {0.6431372549, 0.6666666667, 0.6784313725},
1637     {1, 1, 0.3},
1638     {0.6431372549, 0.6666666667, 0.6784313725},
1639     {0.5, 0.08, 0.12},
1640     {0.6431372549, 0.6666666667, 0.6784313725},
1641     {0.6431372549, 0.6666666667, 0.6784313725},
1642     {0.6431372549, 0.6666666667, 0.6784313725},
1643     {0.6431372549, 0.6666666667, 0.6784313725},
1644     {0.6431372549, 0.6666666667, 0.6784313725},
1645     {0.6431372549, 0.6666666667, 0.6784313725},
1646     {0.6431372549, 0.6666666667, 0.6784313725},
1647     {0.6431372549, 0.6666666667, 0.6784313725},
1648     {0.6431372549, 0.6666666667, 0.6784313725},
1649     {0.6431372549, 0.6666666667, 0.6784313725},
1650     {0.6431372549, 0.6666666667, 0.6784313725},
1651     {0.6431372549, 0.6666666667, 0.6784313725},
1652     {0.6431372549, 0.6666666667, 0.6784313725},
1653     {0.6431372549, 0.6666666667, 0.6784313725},
1654     {0.6431372549, 0.6666666667, 0.6784313725},
1655     {0.6431372549, 0.6666666667, 0.6784313725},
1656     {0.6431372549, 0.6666666667, 0.6784313725},
1657     {0.5, 0.1, 0.5},
1658     {0.6431372549, 0.6666666667, 0.6784313725},
1659     {0.6431372549, 0.6666666667, 0.6784313725},
1660     {0.6431372549, 0.6666666667, 0.6784313725},
1661     {0.6431372549, 0.6666666667, 0.6784313725},
1662     {0.8, 0.8, 0},
1663     {0.6431372549, 0.6666666667, 0.6784313725},
1664     {0.6431372549, 0.6666666667, 0.6784313725},
1665     {0.6431372549, 0.6666666667, 0.6784313725},
1666     {0.6431372549, 0.6666666667, 0.6784313725},
1667     {0.6431372549, 0.6666666667, 0.6784313725},
1668     {1, 0.8431372549, 0},
1669     {0.6431372549, 0.6666666667, 0.6784313725},
1670     {0.6431372549, 0.6666666667, 0.6784313725},
1671     {0.6431372549, 0.6666666667, 0.6784313725},
1672     {0.6431372549, 0.6666666667, 0.6784313725},
1673     {0.6431372549, 0.6666666667, 0.6784313725},
1674     {0.6431372549, 0.6666666667, 0.6784313725},
1675     {0.6431372549, 0.6666666667, 0.6784313725},
1676     {0.6431372549, 0.6666666667, 0.6784313725},
1677     {0.6431372549, 0.6666666667, 0.6784313725},
1678     {0.6431372549, 0.6666666667, 0.6784313725},
1679     {0.6431372549, 0.6666666667, 0.6784313725},
1680     {0.6431372549, 0.6666666667, 0.6784313725},
1681     {0.6431372549, 0.6666666667, 0.6784313725},
1682     {0.6431372549, 0.6666666667, 0.6784313725},
1683     {0.9, 0.8, 0},
1684     {0.6431372549, 0.6666666667, 0.6784313725},
1685     {0.6431372549, 0.6666666667, 0.6784313725},
1686     {0.6431372549, 0.6666666667, 0.6784313725},
1687     {0.6431372549, 0.6666666667, 0.6784313725},
1688     {0.6431372549, 0.6666666667, 0.6784313725},
1689     {0.8, 0.2, 0.2},
1690     {0.6431372549, 0.6666666667, 0.6784313725},
1691     {0.6431372549, 0.6666666667, 0.6784313725},
1692     {0.6431372549, 0.6666666667, 0.6784313725},
1693     {0.6431372549, 0.6666666667, 0.6784313725},
1694     {0.6431372549, 0.6666666667, 0.6784313725},
1695     {0.6431372549, 0.6666666667, 0.6784313725},
1696     {0.6431372549, 0.6666666667, 0.6784313725},
1697     {0.6431372549, 0.6666666667, 0.6784313725},
1698     {0.6431372549, 0.6666666667, 0.6784313725},
1699     {0.6431372549, 0.6666666667, 0.6784313725},
1700     {0.6431372549, 0.6666666667, 0.6784313725},
1701     {0.6431372549, 0.6666666667, 0.6784313725},
1702     {0.1, 0.7, 0.3},
1703     {0.1, 0.3, 0.7},
1704     {0.6431372549, 0.6666666667, 0.6784313725},
1705     {0.6431372549, 0.6666666667, 0.6784313725},
1706     {0.6431372549, 0.6666666667, 0.6784313725},
1707     {0.6431372549, 0.6666666667, 0.6784313725},
1708     {0.6431372549, 0.6666666667, 0.6784313725},
1709     {0.9, 0.8, 0},
1710     {0.6431372549, 0.6666666667, 0.6784313725},
1711     {0.6431372549, 0.6666666667, 0.6784313725},
1712     {0.6431372549, 0.6666666667, 0.6784313725},
1713     {0.6431372549, 0.6666666667, 0.6784313725}
1714   };
1715 
1716   for (int i = 0; i < NELEMENTS; i++)
1717     if (strcmp(element,name[i]) == 0) return rgb[i];
1718   return NULL;
1719 }
1720 
1721 /* ----------------------------------------------------------------------
1722    search the list of element names for the string element
1723    return a pointer to the 3 floating point RGB values
1724    this list is used by AtomEye and is taken from its Mendeleyev.c file
1725 ------------------------------------------------------------------------- */
1726 
element2diam(char * element)1727 double Image::element2diam(char *element)
1728 {
1729   static const char *name[NELEMENTS] = {
1730     "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne",
1731     "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",
1732     "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
1733     "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr",
1734     "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
1735     "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
1736     "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
1737     "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
1738     "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
1739     "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
1740     "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt"
1741   };
1742 
1743   static double diameter[NELEMENTS] = {
1744     0.35, 1.785, 1.45, 1.05, 0.85, 0.72, 0.65, 0.6, 0.5, 1.5662,
1745     1.8, 1.5, 1.4255, 1.07, 1, 1, 1, 1.8597, 2.2, 1.8,
1746     1.6, 1.4, 1.51995, 1.44225, 1.4, 1.43325, 1.35, 1.35, 1.278, 1.35,
1747     1.3, 1.25, 1.15, 1.15, 1.15, 2.0223, 2.35, 2, 1.8, 1.55,
1748     1.6504, 1.3872, 1.35, 1.3, 1.35, 1.4, 1.6, 1.55, 1.55, 1.45,
1749     1.45, 1.4, 1.4, 2.192, 2.6, 2.15, 1.95, 1.85, 1.85, 1.85,
1750     1.85, 1.85, 1.85, 1.8, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75,
1751     1.75, 1.55, 1.6529, 1.5826, 1.35, 1.3, 1.35, 1.35, 1.35, 1.5,
1752     1.9, 1.8, 1.6, 1.9, 1.6, 1.0, 1.0, 2.15, 1.95, 1.8,
1753     1.8, 1.75, 1.75, 1.75, 1.75, 1.0, 1.0, 1.6, 1.6, 1.0,
1754     1.0, 1.0, 1.0, 1.0, 1.6, 1.0, 1.0, 1.0, 1.0
1755   };
1756 
1757   for (int i = 0; i < NELEMENTS; i++)
1758     if (strcmp(element,name[i]) == 0) return diameter[i];
1759   return 0.0;
1760 }
1761 
1762 // ----------------------------------------------------------------------
1763 // ----------------------------------------------------------------------
1764 // ColorMap class
1765 // ----------------------------------------------------------------------
1766 // ----------------------------------------------------------------------
1767 
ColorMap(SPARTA * sparta,Image * caller)1768 ColorMap::ColorMap(SPARTA *sparta, Image *caller) : Pointers(sparta)
1769 {
1770   image = caller;
1771 
1772   // default color map
1773 
1774   dynamic = 1;
1775 
1776   mlo = MINVALUE;
1777   mhi = MAXVALUE;
1778   mstyle = CONTINUOUS;
1779   mrange = FRACTIONAL;
1780 
1781   nentry = 2;
1782   mentry = new MapEntry[nentry];
1783   mentry[0].single = MINVALUE;
1784   mentry[0].color = image->color2rgb("blue");
1785   mentry[1].single = MAXVALUE;
1786   mentry[1].color = image->color2rgb("red");
1787 }
1788 
1789 /* ---------------------------------------------------------------------- */
1790 
~ColorMap()1791 ColorMap::~ColorMap()
1792 {
1793   delete [] mentry;
1794 }
1795 
1796 /* ----------------------------------------------------------------------
1797    redefine color map
1798    args = lo hi style delta N entry1 entry2 ... entryN as defined by caller
1799    return 1 if any error in args, else return 0
1800 ------------------------------------------------------------------------- */
1801 
reset(int narg,char ** arg)1802 int ColorMap::reset(int narg, char **arg)
1803 {
1804   if (!islower(arg[0][0])) {
1805     mlo = NUMERIC;
1806     mlovalue = atof(arg[0]);
1807   } else if (strcmp(arg[0],"min") == 0) mlo = MINVALUE;
1808   else return 1;
1809 
1810   if (!islower(arg[1][0])) {
1811     mhi = NUMERIC;
1812     mhivalue = atof(arg[1]);
1813   } else if (strcmp(arg[1],"max") == 0) mhi = MAXVALUE;
1814   else return 1;
1815 
1816   if (mlo == NUMERIC && mhi == NUMERIC && mlovalue >= mhivalue) return 1;
1817 
1818   if (mlo == MINVALUE || mhi == MAXVALUE) dynamic = 1;
1819   else dynamic = 0;
1820 
1821   if (strlen(arg[2]) != 2) return 1;
1822   if (arg[2][0] == 'c') mstyle = CONTINUOUS;
1823   else if (arg[2][0] == 'd') mstyle = DISCRETE;
1824   else if (arg[2][0] == 's') mstyle = SEQUENTIAL;
1825   else return 1;
1826   if (arg[2][1] == 'a') mrange = ABSOLUTE;
1827   else if (arg[2][1] == 'f') mrange = FRACTIONAL;
1828   else return 1;
1829 
1830   if (mstyle == SEQUENTIAL) {
1831     mbinsize = atof(arg[3]);
1832     if (mbinsize <= 0.0) return 1;
1833     mbinsizeinv = 1.0/mbinsize;
1834   }
1835 
1836   nentry = atoi(arg[4]);
1837   if (nentry < 1) return 1;
1838   delete [] mentry;
1839   mentry = new MapEntry[nentry];
1840 
1841   int n = 5;
1842   for (int i = 0; i < nentry; i++) {
1843     if (mstyle == CONTINUOUS) {
1844       if (n+2 > narg) return 1;
1845       if (!islower(arg[n][0])) {
1846         mentry[i].single = NUMERIC;
1847         mentry[i].svalue = atof(arg[n]);
1848       } else if (strcmp(arg[n],"min") == 0) mentry[i].single = MINVALUE;
1849       else if (strcmp(arg[n],"max") == 0) mentry[i].single = MAXVALUE;
1850       else return 1;
1851       mentry[i].color = image->color2rgb(arg[n+1]);
1852       n += 2;
1853     } else if (mstyle == DISCRETE) {
1854       if (n+3 > narg) return 1;
1855       if (!islower(arg[n][0])) {
1856         mentry[i].lo = NUMERIC;
1857         mentry[i].lvalue = atof(arg[n]);
1858       } else if (strcmp(arg[n],"min") == 0) mentry[i].single = MINVALUE;
1859       else if (strcmp(arg[n],"max") == 0) mentry[i].single = MAXVALUE;
1860       else return 1;
1861       if (!islower(arg[n+1][0])) {
1862         mentry[i].hi = NUMERIC;
1863         mentry[i].hvalue = atof(arg[n+1]);
1864       } else if (strcmp(arg[n+1],"min") == 0) mentry[i].single = MINVALUE;
1865       else if (strcmp(arg[n+1],"max") == 0) mentry[i].single = MAXVALUE;
1866       else return 1;
1867       mentry[i].color = image->color2rgb(arg[n+2]);
1868       n += 3;
1869     } else if (mstyle == SEQUENTIAL) {
1870       if (n+1 > narg) return 1;
1871       mentry[i].color = image->color2rgb(arg[n]);
1872       n += 1;
1873     }
1874     if (mentry[i].color == NULL) return 1;
1875   }
1876 
1877   if (mstyle == CONTINUOUS) {
1878     if (nentry < 2) return 1;
1879     if (mentry[0].single != MINVALUE || mentry[nentry-1].single != MAXVALUE)
1880       return 1;
1881     for (int i = 2; i < nentry-1; i++) {
1882       if (mentry[i].svalue <= mentry[i-1].svalue) return 1;
1883     }
1884   } else if (mstyle == DISCRETE) {
1885     if (nentry < 1) return 1;
1886     if (mentry[nentry-1].lo != MINVALUE || mentry[nentry-1].hi != MAXVALUE)
1887       return 1;
1888   } else if (mstyle == SEQUENTIAL) {
1889     if (nentry < 1) return 1;
1890   }
1891 
1892   // one-time call to minmax if color map is static
1893 
1894   if (!dynamic) return minmax(mlovalue,mhivalue);
1895 
1896   return 0;
1897 }
1898 
1899 /* ----------------------------------------------------------------------
1900    set explicit values for all min/max settings in color map
1901      from min/max dynamic values
1902    set lo/hi current and lvalue/hvalue entries that are MIN/MAX VALUE
1903    called only once if mlo/mhi != MIN/MAX VALUE, else called repeatedly
1904    return 1 = error if any values now overlap incorrectly with dynamic bounds
1905    else return 0
1906 ------------------------------------------------------------------------- */
1907 
minmax(double mindynamic,double maxdynamic)1908 int ColorMap::minmax(double mindynamic, double maxdynamic)
1909 {
1910   if (mlo == MINVALUE) locurrent = mindynamic;
1911   else locurrent = mlovalue;
1912   if (mhi == MAXVALUE) hicurrent = maxdynamic;
1913   else hicurrent = mhivalue;
1914   if (locurrent > hicurrent) return 1;
1915 
1916   if (mstyle == CONTINUOUS) {
1917     if (mrange == ABSOLUTE) mentry[0].svalue = locurrent;
1918     else mentry[0].svalue = 0.0;
1919     if (mrange == ABSOLUTE) mentry[nentry-1].svalue = hicurrent;
1920     else mentry[nentry-1].svalue = 1.0;
1921 
1922     // error in ABSOLUTE mode if new lo/hi current cause
1923     // first/last entry to become lo > hi with adjacent entry
1924 
1925     if (mrange == ABSOLUTE) {
1926       if (mentry[0].svalue > mentry[1].svalue) return 1;
1927       if (mentry[nentry-2].svalue > mentry[nentry-1].svalue) return 1;
1928     }
1929 
1930   // OK if new lo/hi current cause an entry to have lo > hi,
1931   // since last entry will always be a match
1932 
1933   } else if (mstyle == DISCRETE) {
1934     for (int i = 0; i < nentry; i++) {
1935       if (mentry[i].lo == MINVALUE) {
1936         if (mrange == ABSOLUTE) mentry[i].lvalue = locurrent;
1937         else mentry[i].lvalue = 0.0;
1938       }
1939       if (mentry[i].hi == MAXVALUE) {
1940         if (mrange == ABSOLUTE) mentry[i].hvalue = hicurrent;
1941         else mentry[i].hvalue = 1.0;
1942       }
1943     }
1944   }
1945 
1946   return 0;
1947 }
1948 
1949 /* ----------------------------------------------------------------------
1950    convert value into an RGB color via color map
1951    return pointer to 3-vector
1952 ------------------------------------------------------------------------- */
1953 
value2color(double value)1954 double *ColorMap::value2color(double value)
1955 {
1956   double lo;
1957 
1958   value = MAX(value,locurrent);
1959   value = MIN(value,hicurrent);
1960 
1961   if (mrange == FRACTIONAL) {
1962     if (locurrent == hicurrent) value = 0.0;
1963     else value = (value-locurrent) / (hicurrent-locurrent);
1964     lo = 0.0;
1965   } else {
1966     lo = locurrent;
1967   }
1968 
1969   if (mstyle == CONTINUOUS) {
1970     for (int i = 0; i < nentry-1; i++)
1971       if (value >= mentry[i].svalue && value <= mentry[i+1].svalue) {
1972         double fraction = (value-mentry[i].svalue) /
1973           (mentry[i+1].svalue-mentry[i].svalue);
1974         interpolate[0] = mentry[i].color[0] +
1975           fraction*(mentry[i+1].color[0]-mentry[i].color[0]);
1976         interpolate[1] = mentry[i].color[1] +
1977           fraction*(mentry[i+1].color[1]-mentry[i].color[1]);
1978         interpolate[2] = mentry[i].color[2] +
1979           fraction*(mentry[i+1].color[2]-mentry[i].color[2]);
1980         return interpolate;
1981       }
1982   } else if (mstyle == DISCRETE) {
1983     for (int i = 0; i < nentry; i++)
1984       if (value >= mentry[i].lvalue && value <= mentry[i].hvalue)
1985         return mentry[i].color;
1986   } else {
1987     int ibin = static_cast<int> ((value-lo) * mbinsizeinv);
1988     return mentry[ibin%nentry].color;
1989   }
1990 
1991   return NULL;
1992 }
1993