/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
   Contributing author: Nathan Fabian (Sandia)
------------------------------------------------------------------------- */

#include <mpi.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#include "image.h"
#include "math_extra.h"
#include "random_mars.h"
#include "math_const.h"
#include "error.h"
#include "force.h"
#include "memory.h"

#ifdef LAMMPS_JPEG
#include "jpeglib.h"
#endif

#ifdef LAMMPS_PNG
#include <png.h>
#include <zlib.h>
#include <setjmp.h>
#include "version.h"
#endif

using namespace LAMMPS_NS;
using namespace MathConst;

#define NCOLORS 140
#define NELEMENTS 109
#define EPSILON 1.0e-6

enum{NUMERIC,MINVALUE,MAXVALUE};
enum{CONTINUOUS,DISCRETE,SEQUENTIAL};
enum{ABSOLUTE,FRACTIONAL};
enum{NO,YES};

/* ---------------------------------------------------------------------- */

Image::Image(LAMMPS *lmp, int nmap_caller) : Pointers(lmp)
{
  MPI_Comm_rank(world,&me);
  MPI_Comm_size(world,&nprocs);

  // defaults for 3d viz

  width = height = 512;
  theta = 60.0 * MY_PI/180.0;
  phi = 30.0 * MY_PI/180.0;
  zoom = 1.0;
  persp = 0.0;
  shiny = 1.0;
  ssao = NO;

  up[0] = 0.0;
  up[1] = 0.0;
  up[2] = 1.0;

  // colors

  ncolors = 0;
  username = NULL;
  userrgb = NULL;

  boxcolor = color2rgb("yellow");
  background[0] = background[1] = background[2] = 0;

  // define nmap colormaps, all with default settings

  nmap = nmap_caller;
  maps = new ColorMap*[nmap];
  for (int i = 0; i < nmap; i++)
    maps[i] = new ColorMap(lmp,this);

  // static parameters

  FOV = MY_PI/6.0;              // 30 degrees
  ambientColor[0] = 0.0;
  ambientColor[1] = 0.0;
  ambientColor[2] = 0.0;

  keyLightPhi = -MY_PI4;        // -45 degrees
  keyLightTheta = MY_PI/6.0;    // 30 degrees
  keyLightColor[0] = 0.9;
  keyLightColor[1] = 0.9;
  keyLightColor[2] = 0.9;

  fillLightPhi = MY_PI/6.0;     // 30 degrees
  fillLightTheta = 0;
  fillLightColor[0] = 0.45;
  fillLightColor[1] = 0.45;
  fillLightColor[2] = 0.45;

  backLightPhi = MY_PI;         // 180 degrees
  backLightTheta = MY_PI/12.0;  // 15 degrees
  backLightColor[0] = 0.9;
  backLightColor[1] = 0.9;
  backLightColor[2] = 0.9;

  random = NULL;
}

/* ---------------------------------------------------------------------- */

Image::~Image()
{
  for (int i = 0; i < nmap; i++) delete maps[i];
  delete [] maps;

  for (int i = 0; i < ncolors; i++) delete [] username[i];
  memory->sfree(username);
  memory->destroy(userrgb);

  memory->destroy(depthBuffer);
  memory->destroy(surfaceBuffer);
  memory->destroy(imageBuffer);
  memory->destroy(depthcopy);
  memory->destroy(surfacecopy);
  memory->destroy(rgbcopy);

  if (random) delete random;
}

/* ----------------------------------------------------------------------
   allocate image and depth buffers
   called after image size is set
------------------------------------------------------------------------- */

void Image::buffers()
{
  npixels = width * height;
  memory->create(depthBuffer,npixels,"image:depthBuffer");
  memory->create(surfaceBuffer,2*npixels,"image:surfaceBuffer");
  memory->create(imageBuffer,3*npixels,"image:imageBuffer");
  memory->create(depthcopy,npixels,"image:depthcopy");
  memory->create(surfacecopy,2*npixels,"image:surfacecopy");
  memory->create(rgbcopy,3*npixels,"image:rgbcopy");
}

/* ----------------------------------------------------------------------
   reset view parameters
   called once if view is STATIC
   called before every render if view is DYNAMIC
------------------------------------------------------------------------- */

void Image::view_params(double boxxlo, double boxxhi, double boxylo,
                        double boxyhi, double boxzlo, double boxzhi)
{
  // camDir points at the camera, view direction = -camDir

  camDir[0] = sin(theta)*cos(phi);
  camDir[1] = sin(theta)*sin(phi);
  camDir[2] = cos(theta);

  // normalize up vector

  if (up[0] == 0.0 && up[1] == 0.0 && up[2] == 0.0)
    error->all(FLERR,"Invalid image up vector");
  MathExtra::norm3(up);

  // adjust camDir by epsilon if camDir and up are parallel
  // do this by tweaking view direction, not up direction
  // try to insure continuous images as changing view passes thru up
  // sufficient to handle common cases where theta = 0 or 180 is degenerate?

  double dot = MathExtra::dot3(up,camDir);
  if (fabs(dot) > 1.0-EPSILON) {
    if (theta == 0.0) {
      camDir[0] = sin(EPSILON)*cos(phi);
      camDir[1] = sin(EPSILON)*sin(phi);
      camDir[2] = cos(EPSILON);
    } else if (theta == MY_PI) {
      camDir[0] = sin(theta-EPSILON)*cos(phi);
      camDir[1] = sin(theta-EPSILON)*sin(phi);
      camDir[2] = cos(theta-EPSILON);
    } else {
      camDir[0] = sin(theta+EPSILON)*cos(phi);
      camDir[1] = sin(theta+EPSILON)*sin(phi);
      camDir[2] = cos(theta+EPSILON);
    }
  }

  // camUp = camDir x (Up x camDir)

  MathExtra::cross3(up,camDir,camRight);
  MathExtra::norm3(camRight);
  MathExtra::cross3(camDir,camRight,camUp);
  if (camUp[0] == 0.0 && camUp[1] == 0.0 && camUp[2] == 0.0)
    error->all(FLERR,"Invalid image up vector");
  MathExtra::norm3(camUp);

  // zdist = camera distance = function of zoom & bounding box
  // camPos = camera position = function of camDir and zdist

  double delx = 2.0*(boxxhi-boxxlo);
  double dely = 2.0*(boxyhi-boxylo);
  double delz = 2.0*(boxzhi-boxzlo);
  double maxdel = MAX(delx,dely);
  maxdel = MAX(maxdel,delz);

  zdist = maxdel;
  zdist /= tan(FOV);
  zdist += 0.5 * (delx*camDir[0] + dely*camDir[1] + delz*camDir[2]);
  zdist /= zoom;

  camPos[0] = camDir[0] * zdist;
  camPos[1] = camDir[1] * zdist;
  camPos[2] = camDir[2] * zdist;

  // light directions in terms of -camDir = z

  keyLightDir[0] = cos(keyLightTheta) * sin(keyLightPhi);
  keyLightDir[1] = sin(keyLightTheta);
  keyLightDir[2] = cos(keyLightTheta) * cos(keyLightPhi);

  fillLightDir[0] = cos(fillLightTheta) * sin(fillLightPhi);
  fillLightDir[1] = sin(fillLightTheta);
  fillLightDir[2] = cos(fillLightTheta) * cos(fillLightPhi);

  backLightDir[0] = cos(backLightTheta) * sin(backLightPhi);
  backLightDir[1] = sin(backLightTheta);
  backLightDir[2] = cos(backLightTheta) * cos(backLightPhi);

  keyHalfDir[0] = 0 + keyLightDir[0];
  keyHalfDir[1] = 0 + keyLightDir[1];
  keyHalfDir[2] = 1 + keyLightDir[2];
  MathExtra::norm3(keyHalfDir);

  // adjust shinyness of the reflection

  specularHardness = 16.0 * shiny;
  specularIntensity = shiny;

  // adjust strength of the SSAO

  if (ssao) {
    if (!random) random = new RanMars(lmp,seed+me);
    SSAORadius = maxdel * 0.05 * ssaoint;
    SSAOSamples = static_cast<int> (8.0 + 32.0*ssaoint);
    SSAOJitter = MY_PI / 12;
    ambientColor[0] = 0.5;
    ambientColor[1] = 0.5;
    ambientColor[2] = 0.5;
  }

  // param for rasterizing spheres

  tanPerPixel = -(maxdel / (double) height);
}

/* ----------------------------------------------------------------------
   initialize image to background color and depth buffer
   no need to init surfaceBuffer, since will be based on depth
------------------------------------------------------------------------- */

void Image::clear()
{
  int red = background[0];
  int green = background[1];
  int blue = background[2];

  int ix,iy;
  for (iy = 0; iy < height; iy ++)
    for (ix = 0; ix < width; ix ++) {
      imageBuffer[iy * width * 3 + ix * 3 + 0] = red;
      imageBuffer[iy * width * 3 + ix * 3 + 1] = green;
      imageBuffer[iy * width * 3 + ix * 3 + 2] = blue;
      depthBuffer[iy * width + ix] = -1;
    }
}

/* ----------------------------------------------------------------------
   merge image from each processor into one composite image
   done pixel by pixel, respecting depth buffer
   hi procs send to lo procs, cascading down logarithmically
------------------------------------------------------------------------- */

void Image::merge()
{
  MPI_Request requests[3];

  int nhalf = 1;
  while (nhalf < nprocs) nhalf *= 2;
  nhalf /= 2;

  while (nhalf) {
    if (me < nhalf && me+nhalf < nprocs) {
      MPI_Irecv(rgbcopy,npixels*3,MPI_BYTE,me+nhalf,0,world,&requests[0]);
      MPI_Irecv(depthcopy,npixels,MPI_DOUBLE,me+nhalf,0,world,&requests[1]);
      if (ssao)
        MPI_Irecv(surfacecopy,npixels*2,MPI_DOUBLE,
                  me+nhalf,0,world,&requests[2]);
      if (ssao) MPI_Waitall(3,requests,MPI_STATUS_IGNORE);
      else MPI_Waitall(2,requests,MPI_STATUS_IGNORE);

      for (int i = 0; i < npixels; i++) {
        if (depthBuffer[i] < 0 || (depthcopy[i] >= 0 &&
                                   depthcopy[i] < depthBuffer[i])) {
          depthBuffer[i] = depthcopy[i];
          imageBuffer[i*3+0] = rgbcopy[i*3+0];
          imageBuffer[i*3+1] = rgbcopy[i*3+1];
          imageBuffer[i*3+2] = rgbcopy[i*3+2];
          if (ssao) {
            surfaceBuffer[i*2+0] = surfacecopy[i*2+0];
            surfaceBuffer[i*2+1] = surfacecopy[i*2+1];
          }
        }
      }

    } else if (me >= nhalf && me < 2*nhalf) {
      MPI_Send(imageBuffer,npixels*3,MPI_BYTE,me-nhalf,0,world);
      MPI_Send(depthBuffer,npixels,MPI_DOUBLE,me-nhalf,0,world);
      if (ssao) MPI_Send(surfaceBuffer,npixels*2,MPI_DOUBLE,me-nhalf,0,world);
    }

    nhalf /= 2;
  }

  // extra SSAO enhancement
  // bcast full image to all procs
  // each works on subset of pixels
  // gather result back to proc 0

  if (ssao) {
    MPI_Bcast(imageBuffer,npixels*3,MPI_BYTE,0,world);
    MPI_Bcast(surfaceBuffer,npixels*2,MPI_DOUBLE,0,world);
    MPI_Bcast(depthBuffer,npixels,MPI_DOUBLE,0,world);
    compute_SSAO();
    int pixelPart = height/nprocs * width*3;
    MPI_Gather(imageBuffer+me*pixelPart,pixelPart,MPI_BYTE,
               rgbcopy,pixelPart,MPI_BYTE,0,world);
    writeBuffer = rgbcopy;
  } else {
    writeBuffer = imageBuffer;
  }
}

/* ----------------------------------------------------------------------
   draw simulation bounding box as 12 cylinders
------------------------------------------------------------------------- */

void Image::draw_box(double (*corners)[3], double diameter)
{
  draw_cylinder(corners[0],corners[1],boxcolor,diameter,3);
  draw_cylinder(corners[2],corners[3],boxcolor,diameter,3);
  draw_cylinder(corners[0],corners[2],boxcolor,diameter,3);
  draw_cylinder(corners[1],corners[3],boxcolor,diameter,3);
  draw_cylinder(corners[0],corners[4],boxcolor,diameter,3);
  draw_cylinder(corners[1],corners[5],boxcolor,diameter,3);
  draw_cylinder(corners[2],corners[6],boxcolor,diameter,3);
  draw_cylinder(corners[3],corners[7],boxcolor,diameter,3);
  draw_cylinder(corners[4],corners[5],boxcolor,diameter,3);
  draw_cylinder(corners[6],corners[7],boxcolor,diameter,3);
  draw_cylinder(corners[4],corners[6],boxcolor,diameter,3);
  draw_cylinder(corners[5],corners[7],boxcolor,diameter,3);
}

/* ----------------------------------------------------------------------
   draw XYZ axes in red/green/blue
   axes = 4 end points
------------------------------------------------------------------------- */

void Image::draw_axes(double (*axes)[3], double diameter)
{
  draw_cylinder(axes[0],axes[1],color2rgb("red"),diameter,3);
  draw_cylinder(axes[0],axes[2],color2rgb("green"),diameter,3);
  draw_cylinder(axes[0],axes[3],color2rgb("blue"),diameter,3);
}

/* ----------------------------------------------------------------------
   draw sphere at x with surfaceColor and diameter
   render pixel by pixel onto image plane with depth buffering
------------------------------------------------------------------------- */

void Image::draw_sphere(double *x, double *surfaceColor, double diameter)
{
  int ix,iy;
  double projRad;
  double xlocal[3],surface[3];
  double depth;

  xlocal[0] = x[0] - xctr;
  xlocal[1] = x[1] - yctr;
  xlocal[2] = x[2] - zctr;

  double xmap = MathExtra::dot3(camRight,xlocal);
  double ymap = MathExtra::dot3(camUp,xlocal);
  double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(xlocal,camDir);

  double radius = 0.5*diameter;
  double radsq = radius*radius;
  double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
    -tanPerPixel / zoom;
  double pixelRadiusFull = radius / pixelWidth;
  int pixelRadius = static_cast<int> (pixelRadiusFull + 0.5) + 1;

  double xf = xmap / pixelWidth;
  double yf = ymap / pixelWidth;
  int xc = static_cast<int> (xf);
  int yc = static_cast<int> (yf);
  double width_error = xf - xc;
  double height_error = yf - yc;

  // shift 0,0 to screen center (vs lower left)

  xc += width / 2;
  yc += height / 2;

  for (iy = yc - pixelRadius; iy <= yc + pixelRadius; iy++) {
    for (ix = xc - pixelRadius; ix <= xc + pixelRadius; ix++) {
      if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;

      surface[1] = ((iy - yc) - height_error) * pixelWidth;
      surface[0] = ((ix - xc) - width_error) * pixelWidth;
      projRad = surface[0]*surface[0] + surface[1]*surface[1];

      // outside the sphere in the projected image

      if (projRad > radsq) continue;
      surface[2] = sqrt(radsq - projRad);
      depth = dist - surface[2];

      surface[0] /= radius;
      surface[1] /= radius;
      surface[2] /= radius;

      draw_pixel (ix, iy, depth, surface, surfaceColor);
    }
  }
}

/* ----------------------------------------------------------------------
   draw axis oriented cube at x with surfaceColor and diameter in size
   render pixel by pixel onto image plane with depth buffering
------------------------------------------------------------------------- */

void Image::draw_cube(double *x, double *surfaceColor, double diameter)
{
  double xlocal[3],surface[3],normal[3];
  double t,tdir[3];
  double depth;

  xlocal[0] = x[0] - xctr;
  xlocal[1] = x[1] - yctr;
  xlocal[2] = x[2] - zctr;

  double xmap = MathExtra::dot3(camRight,xlocal);
  double ymap = MathExtra::dot3(camUp,xlocal);
  double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(xlocal,camDir);

  double radius = 0.5*diameter;
  double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
    -tanPerPixel / zoom;

  double halfWidth = diameter;
  double pixelHalfWidthFull = halfWidth / pixelWidth;
  int pixelHalfWidth = static_cast<int> (pixelHalfWidthFull + 0.5);

  double xf = xmap / pixelWidth;
  double yf = ymap / pixelWidth;
  int xc = static_cast<int> (xf);
  int yc = static_cast<int> (yf);
  double width_error = xf - xc;
  double height_error = yf - yc;

  // shift 0,0 to screen center (vs lower left)

  xc += width / 2;
  yc += height / 2;

  for (int iy = yc - pixelHalfWidth; iy <= yc + pixelHalfWidth; iy ++) {
    for (int ix = xc - pixelHalfWidth; ix <= xc + pixelHalfWidth; ix ++) {
      if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;

      double sy = ((iy - yc) - height_error) * pixelWidth;
      double sx = ((ix - xc) - width_error) * pixelWidth;
      surface[0] = camRight[0] * sx + camUp[0] * sy;
      surface[1] = camRight[1] * sx + camUp[1] * sy;
      surface[2] = camRight[2] * sx + camUp[2] * sy;

      // iterate through each of the 6 axis-oriented planes of the box
      // only render up to 3 which are facing the camera
      // these checks short circuit a dot product, testing for > 0

      for (int dim = 0; dim < 3; dim ++) {
        if (camDir[dim] > 0) {          // positive faces camera
          t = (radius - surface[dim]) / camDir[dim];
          normal[0] = camRight[dim];
          normal[1] = camUp[dim];
          normal[2] = camDir[dim];
        } else if (camDir[dim] < 0) {   // negative faces camera
          t = -(radius + surface[dim]) / camDir[dim];
          normal[0] = -camRight[dim];
          normal[1] = -camUp[dim];
          normal[2] = -camDir[dim];
        }
        if (camDir[dim] != 0) {
          tdir[0] = camDir[0] * t;
          tdir[1] = camDir[1] * t;
          tdir[2] = camDir[2] * t;

          bool xin = ((surface[0]+tdir[0]) >= -radius) &&
            ((surface[0]+tdir[0]) <= radius);
          bool yin = ((surface[1]+tdir[1]) >= -radius) &&
            ((surface[1]+tdir[1]) <= radius);
          bool zin = ((surface[2]+tdir[2]) >= -radius) &&
            ((surface[2]+tdir[2]) <= radius);

          switch (dim) {
          case 0:
            if (yin & zin) {
              depth = dist - t;
              draw_pixel (ix, iy, depth, normal, surfaceColor);
            }
            break;
          case 1:
            if (xin & zin) {
              depth = dist - t;
              draw_pixel (ix, iy, depth, normal, surfaceColor);
            }
            break;
          case 2:
            if (xin & yin) {
              depth = dist - t;
              draw_pixel (ix, iy, depth, normal, surfaceColor);
            }
            break;
          }
        }
      }
    }
  }
}

/* ----------------------------------------------------------------------
   draw cylinder from x to y with surfaceColor and diameter
   render pixel by pixel onto image plane with depth buffering
   if sflag = 0, draw no end spheres
   if sflag = 1, draw 1st end sphere
   if sflag = 2, draw 2nd end sphere
   if sflag = 3, draw both end spheres
------------------------------------------------------------------------- */

void Image::draw_cylinder(double *x, double *y,
                          double *surfaceColor, double diameter, int sflag)
{
  double surface[3], normal[3];
  double mid[3],xaxis[3],yaxis[3],zaxis[3];
  double camLDir[3], camLRight[3], camLUp[3];
  double zmin, zmax;

  if (sflag % 2) draw_sphere(x,surfaceColor,diameter);
  if (sflag/2) draw_sphere(y,surfaceColor,diameter);

  double radius = 0.5*diameter;
  double radsq = radius*radius;

  zaxis[0] = y[0] - x[0];
  zaxis[1] = y[1] - x[1];
  zaxis[2] = y[2] - x[2];

  double rasterWidth = fabs(MathExtra::dot3(zaxis, camRight)) + diameter;
  double rasterHeight = fabs(MathExtra::dot3(zaxis, camUp)) + diameter;

  mid[0] = (y[0] + x[0]) * 0.5 - xctr;
  mid[1] = (y[1] + x[1]) * 0.5 - yctr;
  mid[2] = (y[2] + x[2]) * 0.5 - zctr;

  double len = MathExtra::len3(zaxis);
  MathExtra::scale3(1.0/len,zaxis);
  len *= 0.5;
  zmax = len;
  zmin = -len;

  double xmap = MathExtra::dot3(camRight,mid);
  double ymap = MathExtra::dot3(camUp,mid);
  double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(mid,camDir);

  double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
    -tanPerPixel / zoom;

  double xf = xmap / pixelWidth;
  double yf = ymap / pixelWidth;
  int xc = static_cast<int> (xf);
  int yc = static_cast<int> (yf);
  double width_error = xf - xc;
  double height_error = yf - yc;

  // shift 0,0 to screen center (vs lower left)

  xc += width / 2;
  yc += height / 2;

  double pixelHalfWidthFull = (rasterWidth * 0.5) / pixelWidth;
  double pixelHalfHeightFull = (rasterHeight * 0.5) / pixelWidth;
  int pixelHalfWidth = static_cast<int> (pixelHalfWidthFull + 0.5);
  int pixelHalfHeight = static_cast<int> (pixelHalfHeightFull + 0.5);

  if (zaxis[0] == camDir[0] && zaxis[1] == camDir[1] && zaxis[2] == camDir[2])
    return;
  if (zaxis[0] == -camDir[0] && zaxis[1] == -camDir[1] &&
      zaxis[2] == -camDir[2]) return;

  MathExtra::cross3(zaxis,camDir,yaxis);
  MathExtra::norm3(yaxis);
  MathExtra::cross3(yaxis,zaxis,xaxis);
  MathExtra::norm3(xaxis);

  camLDir[0] = MathExtra::dot3(camDir,xaxis);
  camLDir[1] = 0.0;
  camLDir[2] = MathExtra::dot3(camDir,zaxis);

  camLRight[0] = MathExtra::dot3(camRight,xaxis);
  camLRight[1] = MathExtra::dot3(camRight,yaxis);
  camLRight[2] = MathExtra::dot3(camRight,zaxis);
  MathExtra::norm3(camLRight);

  camLUp[0] = MathExtra::dot3(camUp,xaxis);
  camLUp[1] = MathExtra::dot3(camUp,yaxis);
  camLUp[2] = MathExtra::dot3(camUp,zaxis);
  MathExtra::norm3(camLUp);

  double a = camLDir[0] * camLDir[0];

  for (int iy = yc - pixelHalfHeight; iy <= yc + pixelHalfHeight; iy ++) {
    for (int ix = xc - pixelHalfWidth; ix <= xc + pixelHalfWidth; ix ++) {
      if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;

      double sy = ((iy - yc) - height_error) * pixelWidth;
      double sx = ((ix - xc) - width_error) * pixelWidth;
      surface[0] = camLRight[0] * sx + camLUp[0] * sy;
      surface[1] = camLRight[1] * sx + camLUp[1] * sy;
      surface[2] = camLRight[2] * sx + camLUp[2] * sy;

      double b = 2 * camLDir[0] * surface[0];
      double c = surface[0] * surface[0] + surface[1] * surface[1] - radsq;

      double partial = b*b - 4*a*c;
      if (partial < 0) continue;
      partial = sqrt (partial);

      double t = (-b + partial) / (2*a);
      double t2 = (-b - partial) / (2*a);
      if (t2 > t) { t = t2; }

      surface[0] += t * camLDir[0];
      surface[1] += t * camLDir[1];
      surface[2] += t * camLDir[2];

      if (surface[2] > zmax || surface[2] < zmin) continue;

      // convert surface into the surface normal

      normal[0] = surface[0] / radius;
      normal[1] = surface[1] / radius;
      normal[2] = 0.0;

      // in camera space

      surface[0] = MathExtra::dot3 (normal, camLRight);
      surface[1] = MathExtra::dot3 (normal, camLUp);
      surface[2] = MathExtra::dot3 (normal, camLDir);

      double depth = dist - t;
      draw_pixel (ix, iy, depth, surface, surfaceColor);
    }
  }
}

/* ----------------------------------------------------------------------
   draw triangle with 3 corner points x,y,z and surfaceColor
------------------------------------------------------------------------- */

void Image::draw_triangle(double *x, double *y, double *z, double *surfaceColor)
{
  double d1[3], d1len, d2[3], d2len, normal[3], invndotd;
  double xlocal[3], ylocal[3], zlocal[3];
  double surface[3];
  double depth;

  xlocal[0] = x[0] - xctr;
  xlocal[1] = x[1] - yctr;
  xlocal[2] = x[2] - zctr;
  ylocal[0] = y[0] - xctr;
  ylocal[1] = y[1] - yctr;
  ylocal[2] = y[2] - zctr;
  zlocal[0] = z[0] - xctr;
  zlocal[1] = z[1] - yctr;
  zlocal[2] = z[2] - zctr;

  MathExtra::sub3 (xlocal, ylocal, d1);
  d1len = MathExtra::len3 (d1);
  MathExtra::scale3 (1.0 / d1len, d1);
  MathExtra::sub3 (zlocal, ylocal, d2);
  d2len = MathExtra::len3 (d2);
  MathExtra::scale3 (1.0 / d2len, d2);

  MathExtra::cross3 (d1, d2, normal);
  MathExtra::norm3 (normal);
  invndotd = 1.0 / MathExtra::dot3(normal, camDir);

  // invalid triangle (parallel)

  if (invndotd == 0) return;

  double r[3],u[3];

  r[0] = MathExtra::dot3(camRight,xlocal);
  r[1] = MathExtra::dot3(camRight,ylocal);
  r[2] = MathExtra::dot3(camRight,zlocal);

  u[0] = MathExtra::dot3(camUp,xlocal);
  u[1] = MathExtra::dot3(camUp,ylocal);
  u[2] = MathExtra::dot3(camUp,zlocal);

  double rasterLeft = r[0] - MIN(r[0],MIN(r[1],r[2]));
  double rasterRight = MAX(r[0],MAX(r[1],r[2])) - r[0];
  double rasterDown = u[0] - MIN(u[0],MIN(u[1],u[2]));
  double rasterUp = MAX(u[0],MAX(u[1],u[2])) - u[0];

  double xmap = MathExtra::dot3(camRight,xlocal);
  double ymap = MathExtra::dot3(camUp,xlocal);
  double dist = MathExtra::dot3(camPos,camDir) - MathExtra::dot3(xlocal,camDir);

  double pixelWidth = (tanPerPixel > 0) ? tanPerPixel * dist :
    -tanPerPixel / zoom;

  double xf = xmap / pixelWidth;
  double yf = ymap / pixelWidth;
  int xc = static_cast<int> (xf);
  int yc = static_cast<int> (yf);
  double width_error = xf - xc;
  double height_error = yf - yc;

  // shift 0,0 to screen center (vs lower left)

  xc += width / 2;
  yc += height / 2;

  double pixelLeftFull = rasterLeft / pixelWidth;
  double pixelRightFull = rasterRight / pixelWidth;
  double pixelDownFull = rasterDown / pixelWidth;
  double pixelUpFull = rasterUp / pixelWidth;
  int pixelLeft = static_cast<int> (pixelLeftFull + 0.5);
  int pixelRight = static_cast<int> (pixelRightFull + 0.5);
  int pixelDown = static_cast<int> (pixelDownFull + 0.5);
  int pixelUp = static_cast<int> (pixelUpFull + 0.5);

  for (int iy = yc - pixelDown; iy <= yc + pixelUp; iy ++) {
    for (int ix = xc - pixelLeft; ix <= xc + pixelRight; ix ++) {
      if (iy < 0 || iy >= height || ix < 0 || ix >= width) continue;

      double sy = ((iy - yc) - height_error) * pixelWidth;
      double sx = ((ix - xc) - width_error) * pixelWidth;
      surface[0] = camRight[0] * sx + camUp[0] * sy;
      surface[1] = camRight[1] * sx + camUp[1] * sy;
      surface[2] = camRight[2] * sx + camUp[2] * sy;

      double t = -MathExtra::dot3(normal,surface) * invndotd;

      // test inside the triangle

      double p[3];
      p[0] = xlocal[0] + surface[0] + camDir[0] * t;
      p[1] = xlocal[1] + surface[1] + camDir[1] * t;
      p[2] = xlocal[2] + surface[2] + camDir[2] * t;

      double s1[3], s2[3], s3[3];
      double c1[3], c2[3];

      MathExtra::sub3 (zlocal, xlocal, s1);
      MathExtra::sub3 (ylocal, xlocal, s2);
      MathExtra::sub3 (p, xlocal, s3);
      MathExtra::cross3 (s1, s2, c1);
      MathExtra::cross3 (s1, s3, c2);
      if (MathExtra::dot3 (c1, c2) <= 0) continue;

      MathExtra::sub3 (xlocal, ylocal, s1);
      MathExtra::sub3 (zlocal, ylocal, s2);
      MathExtra::sub3 (p, ylocal, s3);
      MathExtra::cross3 (s1, s2, c1);
      MathExtra::cross3 (s1, s3, c2);
      if (MathExtra::dot3 (c1, c2) <= 0) continue;

      MathExtra::sub3 (ylocal, zlocal, s1);
      MathExtra::sub3 (xlocal, zlocal, s2);
      MathExtra::sub3 (p, zlocal, s3);
      MathExtra::cross3 (s1, s2, c1);
      MathExtra::cross3 (s1, s3, c2);
      if (MathExtra::dot3 (c1, c2) <= 0) continue;

      double cNormal[3];
      cNormal[0] = MathExtra::dot3(camRight, normal);
      cNormal[1] = MathExtra::dot3(camUp, normal);
      cNormal[2] = MathExtra::dot3(camDir, normal);

      depth = dist - t;
      draw_pixel(ix,iy,depth,cNormal,surfaceColor);
    }
  }
}

/* ---------------------------------------------------------------------- */

void Image::draw_pixel(int ix, int iy, double depth,
                           double *surface, double *surfaceColor)
{
  double diffuseKey,diffuseFill,diffuseBack,specularKey;
  if (depth < 0 || (depthBuffer[ix + iy*width] >= 0 &&
                    depth >= depthBuffer[ix + iy*width])) return;
  depthBuffer[ix + iy*width] = depth;

  // store only the tangent relative to the camera normal (0,0,-1)

  surfaceBuffer[0 + ix * 2 + iy*width * 2] = surface[1];
  surfaceBuffer[1 + ix * 2 + iy*width * 2] = -surface[0];

  diffuseKey = saturate(MathExtra::dot3(surface, keyLightDir));
  diffuseFill = saturate(MathExtra::dot3(surface, fillLightDir));
  diffuseBack = saturate(MathExtra::dot3(surface, backLightDir));
  specularKey = pow(saturate(MathExtra::dot3(surface, keyHalfDir)),
                    specularHardness) * specularIntensity;

  double c[3];
  c[0] = surfaceColor[0] * ambientColor[0];
  c[1] = surfaceColor[1] * ambientColor[1];
  c[2] = surfaceColor[2] * ambientColor[2];

  c[0] += surfaceColor[0] * keyLightColor[0] * diffuseKey;
  c[1] += surfaceColor[1] * keyLightColor[1] * diffuseKey;
  c[2] += surfaceColor[2] * keyLightColor[2] * diffuseKey;

  c[0] += keyLightColor[0] * specularKey;
  c[1] += keyLightColor[1] * specularKey;
  c[2] += keyLightColor[2] * specularKey;

  c[0] += surfaceColor[0] * fillLightColor[0] * diffuseFill;
  c[1] += surfaceColor[1] * fillLightColor[1] * diffuseFill;
  c[2] += surfaceColor[2] * fillLightColor[2] * diffuseFill;

  c[0] += surfaceColor[0] * backLightColor[0] * diffuseBack;
  c[1] += surfaceColor[1] * backLightColor[1] * diffuseBack;
  c[2] += surfaceColor[2] * backLightColor[2] * diffuseBack;

  c[0] = saturate(c[0]);
  c[1] = saturate(c[1]);
  c[2] = saturate(c[2]);

  imageBuffer[0 + ix*3 + iy*width*3] = static_cast<int>(c[0] * 255.0);
  imageBuffer[1 + ix*3 + iy*width*3] = static_cast<int>(c[1] * 255.0);
  imageBuffer[2 + ix*3 + iy*width*3] = static_cast<int>(c[2] * 255.0);
}

/* ---------------------------------------------------------------------- */

void Image::compute_SSAO()
{
  // used for rasterizing the spheres

  double delTheta = 2.0*MY_PI / SSAOSamples;

  // typical neighborhood value for shading

  double pixelWidth = (tanPerPixel > 0) ? tanPerPixel :
        -tanPerPixel / zoom;
  int pixelRadius = (int) trunc (SSAORadius / pixelWidth + 0.5);

  int x,y,s;
  int hPart = height / nprocs;
  int index = me * hPart * width;
  for (y = me * hPart; y < (me + 1) * hPart; y ++) {
    for (x = 0; x < width; x ++, index ++) {
      double cdepth = depthBuffer[index];
      if (cdepth < 0) { continue; }

      double sx = surfaceBuffer[index * 2 + 0];
      double sy = surfaceBuffer[index * 2 + 1];
      double sin_t = -sqrt(sx*sx + sy*sy);

      double mytheta = random->uniform() * SSAOJitter;
      double ao = 0.0;

      for (s = 0; s < SSAOSamples; s ++) {
        double hx = cos(mytheta);
        double hy = sin(mytheta);
        mytheta += delTheta;

        // multiply by z cross surface tangent
        // so that dot (aka cos) works here

        double scaled_sin_t = sin_t * (hx*sy + hy*sx);

        // Bresenham's line algorithm to march over depthBuffer

        int dx = static_cast<int> (hx * pixelRadius);
        int dy = static_cast<int> (hy * pixelRadius);
        int ex = x + dx;
        if (ex < 0) { ex = 0; } if (ex >= width) { ex = width - 1; }
        int ey = y + dy;
        if (ey < 0) { ey = 0; } if (ey >= height) { ey = height - 1; }
        double delta;
        int small, large;
        double lenIncr;
        if (fabs(hx) > fabs(hy)) {
          small = (hx > 0) ? 1 : -1;
          large = (hy > 0) ? width : -width;
          delta = fabs(hy / hx);
        } else {
          small = (hy > 0) ? width : -width;
          large = (hx > 0) ? 1 : -1;
          delta = fabs(hx / hy);
        }
        lenIncr = sqrt (1 + delta * delta) * pixelWidth;

        // initialize with one step
        // because the center point doesn't need testing

        int end = ex + ey * width;
        int ind = index + small;
        double len = lenIncr;
        double err = delta;
        if (err >= 1.0) {
          ind += large;
          err -= 1.0;
        }

        double minPeak = -1;
        double peakLen = 0.0;
        int stepsTaken = 1;
        while ((small > 0 && ind <= end) || (small < 0 && ind >= end)) {
          if (ind < 0 || ind >= (width*height)) {
            break;
          }

          // cdepth - depthBuffer B/C we want it in the negative z direction

          if (minPeak < 0 || (depthBuffer[ind] >= 0 &&
                              depthBuffer[ind] < minPeak)) {
            minPeak = depthBuffer[ind];
            peakLen = len;
          }
          ind += small;
          len += lenIncr;
          err += delta;
          if (err >= 1.0) {
            ind += large;
            err -= 1.0;
          }
          stepsTaken ++;
        }

        if (peakLen > 0) {
          double h = atan ((cdepth - minPeak) / peakLen);
          ao += saturate(sin (h) - scaled_sin_t);
        } else {
          ao += saturate(-scaled_sin_t);
        }
      }
      ao /= (double)SSAOSamples;

      double c[3];
      c[0] = (double) (*(unsigned char *) &imageBuffer[index * 3 + 0]);
      c[1] = (double) (*(unsigned char *) &imageBuffer[index * 3 + 1]);
      c[2] = (double) (*(unsigned char *) &imageBuffer[index * 3 + 2]);
      c[0] *= (1.0 - ao);
      c[1] *= (1.0 - ao);
      c[2] *= (1.0 - ao);
      imageBuffer[index * 3 + 0] = (int) c[0];
      imageBuffer[index * 3 + 1] = (int) c[1];
      imageBuffer[index * 3 + 2] = (int) c[2];
    }
  }
}

/* ---------------------------------------------------------------------- */

void Image::write_JPG(FILE *fp)
{
#ifdef LAMMPS_JPEG
  struct jpeg_compress_struct cinfo;
  struct jpeg_error_mgr jerr;
  JSAMPROW row_pointer;

  cinfo.err = jpeg_std_error(&jerr);
  jpeg_create_compress(&cinfo);
  jpeg_stdio_dest(&cinfo,fp);
  cinfo.image_width = width;
  cinfo.image_height = height;
  cinfo.input_components = 3;
  cinfo.in_color_space = JCS_RGB;

  jpeg_set_defaults(&cinfo);
  jpeg_set_quality(&cinfo,85,TRUE);
  jpeg_start_compress(&cinfo,TRUE);

  while (cinfo.next_scanline < cinfo.image_height) {
    row_pointer = (JSAMPROW)
      &writeBuffer[(cinfo.image_height - 1 - cinfo.next_scanline) * 3 * width];
    jpeg_write_scanlines(&cinfo,&row_pointer,1);
  }

  jpeg_finish_compress(&cinfo);
  jpeg_destroy_compress(&cinfo);
#endif
}

/* ---------------------------------------------------------------------- */

void Image::write_PNG(FILE *fp)
{
#ifdef LAMMPS_PNG
  png_structp png_ptr;
  png_infop info_ptr;

  png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);

  if (!png_ptr) return;

  info_ptr = png_create_info_struct(png_ptr);
  if (!info_ptr) {
    png_destroy_write_struct(&png_ptr, NULL);
    return;
  }

  if (setjmp(png_jmpbuf(png_ptr))) {
    png_destroy_write_struct(&png_ptr, &info_ptr);
    return;
  }

  png_init_io(png_ptr, fp);
  png_set_compression_level(png_ptr,Z_BEST_COMPRESSION);
  png_set_IHDR(png_ptr,info_ptr,width,height,8,PNG_COLOR_TYPE_RGB,
    PNG_INTERLACE_NONE,PNG_COMPRESSION_TYPE_DEFAULT,PNG_FILTER_TYPE_DEFAULT);

  png_text text_ptr[2];
  memset(text_ptr,0,2*sizeof(png_text));

  char key0[]  = "Software";
  char text0[] = "LAMMPS " LAMMPS_VERSION;
  char key1[]  = "Description";
  char text1[] = "Dump image snapshot";
  text_ptr[0].key = key0;
  text_ptr[0].text = text0;
  text_ptr[1].key = key1;
  text_ptr[1].text = text1;
  text_ptr[0].compression = PNG_TEXT_COMPRESSION_NONE;
  text_ptr[1].compression = PNG_TEXT_COMPRESSION_NONE;

  png_set_text(png_ptr,info_ptr,text_ptr,1);
  png_write_info(png_ptr,info_ptr);

  png_bytep row_pointers[height];
  for (int i=0; i < height; ++i)
    row_pointers[i] = (png_bytep) &writeBuffer[(height-i-1)*3*width];

  png_write_image(png_ptr, row_pointers);
  png_write_end(png_ptr, info_ptr);

  png_destroy_write_struct(&png_ptr, &info_ptr);
#endif
}

/* ---------------------------------------------------------------------- */

void Image::write_PPM(FILE *fp)
{
  fprintf(fp,"P6\n%d %d\n255\n",width,height);

  int y;
  for (y = height-1; y >= 0; y--)
    fwrite(&writeBuffer[y*width*3],3,width,fp);
}

/* ----------------------------------------------------------------------
   return static/dynamic status of color map index
------------------------------------------------------------------------- */

int Image::map_dynamic(int index)
{
  return maps[index]->dynamic;
}

/* ----------------------------------------------------------------------
   redefine properties of the color map index
   return 1 if any error in args, else return 0
------------------------------------------------------------------------- */

int Image::map_reset(int index, int narg, char **arg)
{
  return maps[index]->reset(narg,arg);
}

/* ----------------------------------------------------------------------
   set min/max bounds of dynamic color map index
------------------------------------------------------------------------- */

int Image::map_minmax(int index, double mindynamic, double maxdynamic)
{
  return maps[index]->minmax(mindynamic,maxdynamic);
}

/* ----------------------------------------------------------------------
   return 3-vector color corresponding to value from color map index
------------------------------------------------------------------------- */

double *Image::map_value2color(int index, double value)
{
  return maps[index]->value2color(value);
}

/* ----------------------------------------------------------------------
   add a new color to username and userrgb
   redefine RGB values in userrgb if name already exists
   return 1 if RGB values are invalid, else return 0
------------------------------------------------------------------------- */

int Image::addcolor(char *name, double r, double g, double b)
{
  int icolor;
  for (icolor = 0; icolor < ncolors; icolor++)
    if (strcmp(name,username[icolor]) == 0) break;

  if (icolor == ncolors) {
    username = (char **)
      memory->srealloc(username,(ncolors+1)*sizeof(char *),"image:username");
    memory->grow(userrgb,ncolors+1,3,"image:userrgb");
    ncolors++;
  }

  int n = strlen(name) + 1;
  username[icolor] = new char[n];
  strcpy(username[icolor],name);

  if (r < 0.0 || r > 1.0 || g < 0.0 || g > 1.0 || b < 0.0 || b > 1.0)
    return 1;

  userrgb[icolor][0] = r;
  userrgb[icolor][1] = g;
  userrgb[icolor][2] = b;

  return 0;
}

/* ----------------------------------------------------------------------
   if index > 0, return ptr to index-1 color from rgb
   if index < 0, return ptr to -index-1 color from userrgb
   if index = 0, search the 2 lists of color names for the string color
   search user-defined color names first, then the list of NCOLORS names
   return a pointer to the 3 floating point RGB values or NULL if didn't find
------------------------------------------------------------------------- */

double *Image::color2rgb(const char *color, int index)
{
  static const char *name[NCOLORS] = {
    "aliceblue",
    "antiquewhite",
    "aqua",
    "aquamarine",
    "azure",
    "beige",
    "bisque",
    "black",
    "blanchedalmond",
    "blue",
    "blueviolet",
    "brown",
    "burlywood",
    "cadetblue",
    "chartreuse",
    "chocolate",
    "coral",
    "cornflowerblue",
    "cornsilk",
    "crimson",
    "cyan",
    "darkblue",
    "darkcyan",
    "darkgoldenrod",
    "darkgray",
    "darkgreen",
    "darkkhaki",
    "darkmagenta",
    "darkolivegreen",
    "darkorange",
    "darkorchid",
    "darkred",
    "darksalmon",
    "darkseagreen",
    "darkslateblue",
    "darkslategray",
    "darkturquoise",
    "darkviolet",
    "deeppink",
    "deepskyblue",
    "dimgray",
    "dodgerblue",
    "firebrick",
    "floralwhite",
    "forestgreen",
    "fuchsia",
    "gainsboro",
    "ghostwhite",
    "gold",
    "goldenrod",
    "gray",
    "green",
    "greenyellow",
    "honeydew",
    "hotpink",
    "indianred",
    "indigo",
    "ivory",
    "khaki",
    "lavender",
    "lavenderblush",
    "lawngreen",
    "lemonchiffon",
    "lightblue",
    "lightcoral",
    "lightcyan",
    "lightgoldenrodyellow",
    "lightgreen",
    "lightgrey",
    "lightpink",
    "lightsalmon",
    "lightseagreen",
    "lightskyblue",
    "lightslategray",
    "lightsteelblue",
    "lightyellow",
    "lime",
    "limegreen",
    "linen",
    "magenta",
    "maroon",
    "mediumaquamarine",
    "mediumblue",
    "mediumorchid",
    "mediumpurple",
    "mediumseagreen",
    "mediumslateblue",
    "mediumspringgreen",
    "mediumturquoise",
    "mediumvioletred",
    "midnightblue",
    "mintcream",
    "mistyrose",
    "moccasin",
    "navajowhite",
    "navy",
    "oldlace",
    "olive",
    "olivedrab",
    "orange",
    "orangered",
    "orchid",
    "palegoldenrod",
    "palegreen",
    "paleturquoise",
    "palevioletred",
    "papayawhip",
    "peachpuff",
    "peru",
    "pink",
    "plum",
    "powderblue",
    "purple",
    "red",
    "rosybrown",
    "royalblue",
    "saddlebrown",
    "salmon",
    "sandybrown",
    "seagreen",
    "seashell",
    "sienna",
    "silver",
    "skyblue",
    "slateblue",
    "slategray",
    "snow",
    "springgreen",
    "steelblue",
    "tan",
    "teal",
    "thistle",
    "tomato",
    "turquoise",
    "violet",
    "wheat",
    "white",
    "whitesmoke",
    "yellow",
    "yellowgreen"
  };

  static double rgb[NCOLORS][3] = {
    {240/255.0, 248/255.0, 255/255.0},
    {250/255.0, 235/255.0, 215/255.0},
    {0/255.0, 255/255.0, 255/255.0},
    {127/255.0, 255/255.0, 212/255.0},
    {240/255.0, 255/255.0, 255/255.0},
    {245/255.0, 245/255.0, 220/255.0},
    {255/255.0, 228/255.0, 196/255.0},
    {0/255.0, 0/255.0, 0/255.0},
    {255/255.0, 255/255.0, 205/255.0},
    {0/255.0, 0/255.0, 255/255.0},
    {138/255.0, 43/255.0, 226/255.0},
    {165/255.0, 42/255.0, 42/255.0},
    {222/255.0, 184/255.0, 135/255.0},
    {95/255.0, 158/255.0, 160/255.0},
    {127/255.0, 255/255.0, 0/255.0},
    {210/255.0, 105/255.0, 30/255.0},
    {255/255.0, 127/255.0, 80/255.0},
    {100/255.0, 149/255.0, 237/255.0},
    {255/255.0, 248/255.0, 220/255.0},
    {220/255.0, 20/255.0, 60/255.0},
    {0/255.0, 255/255.0, 255/255.0},
    {0/255.0, 0/255.0, 139/255.0},
    {0/255.0, 139/255.0, 139/255.0},
    {184/255.0, 134/255.0, 11/255.0},
    {169/255.0, 169/255.0, 169/255.0},
    {0/255.0, 100/255.0, 0/255.0},
    {189/255.0, 183/255.0, 107/255.0},
    {139/255.0, 0/255.0, 139/255.0},
    {85/255.0, 107/255.0, 47/255.0},
    {255/255.0, 140/255.0, 0/255.0},
    {153/255.0, 50/255.0, 204/255.0},
    {139/255.0, 0/255.0, 0/255.0},
    {233/255.0, 150/255.0, 122/255.0},
    {143/255.0, 188/255.0, 143/255.0},
    {72/255.0, 61/255.0, 139/255.0},
    {47/255.0, 79/255.0, 79/255.0},
    {0/255.0, 206/255.0, 209/255.0},
    {148/255.0, 0/255.0, 211/255.0},
    {255/255.0, 20/255.0, 147/255.0},
    {0/255.0, 191/255.0, 255/255.0},
    {105/255.0, 105/255.0, 105/255.0},
    {30/255.0, 144/255.0, 255/255.0},
    {178/255.0, 34/255.0, 34/255.0},
    {255/255.0, 250/255.0, 240/255.0},
    {34/255.0, 139/255.0, 34/255.0},
    {255/255.0, 0/255.0, 255/255.0},
    {220/255.0, 220/255.0, 220/255.0},
    {248/255.0, 248/255.0, 255/255.0},
    {255/255.0, 215/255.0, 0/255.0},
    {218/255.0, 165/255.0, 32/255.0},
    {128/255.0, 128/255.0, 128/255.0},
    {0/255.0, 128/255.0, 0/255.0},
    {173/255.0, 255/255.0, 47/255.0},
    {240/255.0, 255/255.0, 240/255.0},
    {255/255.0, 105/255.0, 180/255.0},
    {205/255.0, 92/255.0, 92/255.0},
    {75/255.0, 0/255.0, 130/255.0},
    {255/255.0, 240/255.0, 240/255.0},
    {240/255.0, 230/255.0, 140/255.0},
    {230/255.0, 230/255.0, 250/255.0},
    {255/255.0, 240/255.0, 245/255.0},
    {124/255.0, 252/255.0, 0/255.0},
    {255/255.0, 250/255.0, 205/255.0},
    {173/255.0, 216/255.0, 230/255.0},
    {240/255.0, 128/255.0, 128/255.0},
    {224/255.0, 255/255.0, 255/255.0},
    {250/255.0, 250/255.0, 210/255.0},
    {144/255.0, 238/255.0, 144/255.0},
    {211/255.0, 211/255.0, 211/255.0},
    {255/255.0, 182/255.0, 193/255.0},
    {255/255.0, 160/255.0, 122/255.0},
    {32/255.0, 178/255.0, 170/255.0},
    {135/255.0, 206/255.0, 250/255.0},
    {119/255.0, 136/255.0, 153/255.0},
    {176/255.0, 196/255.0, 222/255.0},
    {255/255.0, 255/255.0, 224/255.0},
    {0/255.0, 255/255.0, 0/255.0},
    {50/255.0, 205/255.0, 50/255.0},
    {250/255.0, 240/255.0, 230/255.0},
    {255/255.0, 0/255.0, 255/255.0},
    {128/255.0, 0/255.0, 0/255.0},
    {102/255.0, 205/255.0, 170/255.0},
    {0/255.0, 0/255.0, 205/255.0},
    {186/255.0, 85/255.0, 211/255.0},
    {147/255.0, 112/255.0, 219/255.0},
    {60/255.0, 179/255.0, 113/255.0},
    {123/255.0, 104/255.0, 238/255.0},
    {0/255.0, 250/255.0, 154/255.0},
    {72/255.0, 209/255.0, 204/255.0},
    {199/255.0, 21/255.0, 133/255.0},
    {25/255.0, 25/255.0, 112/255.0},
    {245/255.0, 255/255.0, 250/255.0},
    {255/255.0, 228/255.0, 225/255.0},
    {255/255.0, 228/255.0, 181/255.0},
    {255/255.0, 222/255.0, 173/255.0},
    {0/255.0, 0/255.0, 128/255.0},
    {253/255.0, 245/255.0, 230/255.0},
    {128/255.0, 128/255.0, 0/255.0},
    {107/255.0, 142/255.0, 35/255.0},
    {255/255.0, 165/255.0, 0/255.0},
    {255/255.0, 69/255.0, 0/255.0},
    {218/255.0, 112/255.0, 214/255.0},
    {238/255.0, 232/255.0, 170/255.0},
    {152/255.0, 251/255.0, 152/255.0},
    {175/255.0, 238/255.0, 238/255.0},
    {219/255.0, 112/255.0, 147/255.0},
    {255/255.0, 239/255.0, 213/255.0},
    {255/255.0, 239/255.0, 213/255.0},
    {205/255.0, 133/255.0, 63/255.0},
    {255/255.0, 192/255.0, 203/255.0},
    {221/255.0, 160/255.0, 221/255.0},
    {176/255.0, 224/255.0, 230/255.0},
    {128/255.0, 0/255.0, 128/255.0},
    {255/255.0, 0/255.0, 0/255.0},
    {188/255.0, 143/255.0, 143/255.0},
    {65/255.0, 105/255.0, 225/255.0},
    {139/255.0, 69/255.0, 19/255.0},
    {250/255.0, 128/255.0, 114/255.0},
    {244/255.0, 164/255.0, 96/255.0},
    {46/255.0, 139/255.0, 87/255.0},
    {255/255.0, 245/255.0, 238/255.0},
    {160/255.0, 82/255.0, 45/255.0},
    {192/255.0, 192/255.0, 192/255.0},
    {135/255.0, 206/255.0, 235/255.0},
    {106/255.0, 90/255.0, 205/255.0},
    {112/255.0, 128/255.0, 144/255.0},
    {255/255.0, 250/255.0, 250/255.0},
    {0/255.0, 255/255.0, 127/255.0},
    {70/255.0, 130/255.0, 180/255.0},
    {210/255.0, 180/255.0, 140/255.0},
    {0/255.0, 128/255.0, 128/255.0},
    {216/255.0, 191/255.0, 216/255.0},
    {253/255.0, 99/255.0, 71/255.0},
    {64/255.0, 224/255.0, 208/255.0},
    {238/255.0, 130/255.0, 238/255.0},
    {245/255.0, 222/255.0, 179/255.0},
    {255/255.0, 255/255.0, 255/255.0},
    {245/255.0, 245/255.0, 245/255.0},
    {255/255.0, 255/255.0, 0/255.0},
    {154/255.0, 205/255.0, 50/255.0}
  };

  if (index > 0) {
    if (index > NCOLORS) return NULL;
    return rgb[index-1];
  }
  if (index < 0) {
    if (-index > ncolors) return NULL;
    return userrgb[-index-1];
  }

  for (int i = 0; i < ncolors; i++)
    if (strcmp(color,username[i]) == 0) return userrgb[i];
  for (int i = 0; i < NCOLORS; i++)
    if (strcmp(color,name[i]) == 0) return rgb[i];
  return NULL;
}

/* ----------------------------------------------------------------------
   return number of default colors
------------------------------------------------------------------------- */

int Image::default_colors()
{
  return NCOLORS;
}

/* ----------------------------------------------------------------------
   search the list of element names for the string element
   return a pointer to the 3 floating point RGB values
   this list is used by AtomEye and is taken from its Mendeleyev.c file
------------------------------------------------------------------------- */

double *Image::element2color(char *element)
{
  static const char *name[NELEMENTS] = {
    "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne",
    "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",
    "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
    "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr",
    "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
    "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
    "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
    "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
    "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
    "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
    "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt"
  };

  static double rgb[NELEMENTS][3] = {
    {0.8, 0.8, 0.8},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.7, 0.7, 0.7},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.9, 0.4, 0},
    {0.35, 0.35, 0.35},
    {0.2, 0.2, 0.8},
    {0.8, 0.2, 0.2},
    {0.7, 0.85, 0.45},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6, 0.6, 0.6},
    {0.6, 0.6, 0.7},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6901960784, 0.768627451, 0.8705882353},
    {0.1, 0.7, 0.3},
    {0.95, 0.9, 0.2},
    {0.15, 0.5, 0.1},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.5, 0.5, 0.5},
    {0.8, 0.8, 0.7},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0, 0.8, 0},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.5176470588, 0.5764705882, 0.6529411765},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.257254902, 0.2666666667, 0.271372549},
    {0.95, 0.7900735294, 0.01385869565},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.9, 0, 1},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {1, 1, 0.3},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.5, 0.08, 0.12},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.5, 0.1, 0.5},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.8, 0.8, 0},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {1, 0.8431372549, 0},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.9, 0.8, 0},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.8, 0.2, 0.2},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.1, 0.7, 0.3},
    {0.1, 0.3, 0.7},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.9, 0.8, 0},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725},
    {0.6431372549, 0.6666666667, 0.6784313725}
  };

  for (int i = 0; i < NELEMENTS; i++)
    if (strcmp(element,name[i]) == 0) return rgb[i];
  return NULL;
}

/* ----------------------------------------------------------------------
   search the list of element names for the string element
   return a pointer to the 3 floating point RGB values
   this list is used by AtomEye and is taken from its Mendeleyev.c file
------------------------------------------------------------------------- */

double Image::element2diam(char *element)
{
  static const char *name[NELEMENTS] = {
    "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne",
    "Na", "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca",
    "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
    "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr",
    "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
    "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
    "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
    "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
    "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
    "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
    "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt"
  };

  static double diameter[NELEMENTS] = {
    0.35, 1.785, 1.45, 1.05, 0.85, 0.72, 0.65, 0.6, 0.5, 1.5662,
    1.8, 1.5, 1.4255, 1.07, 1, 1, 1, 1.8597, 2.2, 1.8,
    1.6, 1.4, 1.51995, 1.44225, 1.4, 1.43325, 1.35, 1.35, 1.278, 1.35,
    1.3, 1.25, 1.15, 1.15, 1.15, 2.0223, 2.35, 2, 1.8, 1.55,
    1.6504, 1.3872, 1.35, 1.3, 1.35, 1.4, 1.6, 1.55, 1.55, 1.45,
    1.45, 1.4, 1.4, 2.192, 2.6, 2.15, 1.95, 1.85, 1.85, 1.85,
    1.85, 1.85, 1.85, 1.8, 1.75, 1.75, 1.75, 1.75, 1.75, 1.75,
    1.75, 1.55, 1.6529, 1.5826, 1.35, 1.3, 1.35, 1.35, 1.35, 1.5,
    1.9, 1.8, 1.6, 1.9, 1.6, 1.0, 1.0, 2.15, 1.95, 1.8,
    1.8, 1.75, 1.75, 1.75, 1.75, 1.0, 1.0, 1.6, 1.6, 1.0,
    1.0, 1.0, 1.0, 1.0, 1.6, 1.0, 1.0, 1.0, 1.0
  };

  for (int i = 0; i < NELEMENTS; i++)
    if (strcmp(element,name[i]) == 0) return diameter[i];
  return 0.0;
}

// ----------------------------------------------------------------------
// ----------------------------------------------------------------------
// ColorMap class
// ----------------------------------------------------------------------
// ----------------------------------------------------------------------

ColorMap::ColorMap(LAMMPS *lmp, Image *caller) : Pointers(lmp)
{
  image = caller;

  // default color map

  dynamic = 1;

  mlo = MINVALUE;
  mhi = MAXVALUE;
  mstyle = CONTINUOUS;
  mrange = FRACTIONAL;

  nentry = 2;
  mentry = new MapEntry[nentry];
  mentry[0].single = MINVALUE;
  mentry[0].color = image->color2rgb("blue");
  mentry[1].single = MAXVALUE;
  mentry[1].color = image->color2rgb("red");
}

/* ---------------------------------------------------------------------- */

ColorMap::~ColorMap()
{
  delete [] mentry;
}

/* ----------------------------------------------------------------------
   redefine color map
   args = lo hi style delta N entry1 entry2 ... entryN as defined by caller
   return 1 if any error in args, else return 0
------------------------------------------------------------------------- */

int ColorMap::reset(int narg, char **arg)
{
  if (!islower(arg[0][0])) {
    mlo = NUMERIC;
    mlovalue = force->numeric(FLERR,arg[0]);
  } else if (strcmp(arg[0],"min") == 0) mlo = MINVALUE;
  else return 1;

  if (!islower(arg[1][0])) {
    mhi = NUMERIC;
    mhivalue = force->numeric(FLERR,arg[1]);
  } else if (strcmp(arg[1],"max") == 0) mhi = MAXVALUE;
  else return 1;

  if (mlo == NUMERIC && mhi == NUMERIC && mlovalue >= mhivalue) return 1;

  if (mlo == MINVALUE || mhi == MAXVALUE) dynamic = 1;
  else dynamic = 0;

  if (strlen(arg[2]) != 2) return 1;
  if (arg[2][0] == 'c') mstyle = CONTINUOUS;
  else if (arg[2][0] == 'd') mstyle = DISCRETE;
  else if (arg[2][0] == 's') mstyle = SEQUENTIAL;
  else return 1;
  if (arg[2][1] == 'a') mrange = ABSOLUTE;
  else if (arg[2][1] == 'f') mrange = FRACTIONAL;
  else return 1;

  if (mstyle == SEQUENTIAL) {
    mbinsize = force->numeric(FLERR,arg[3]);
    if (mbinsize <= 0.0) return 1;
    mbinsizeinv = 1.0/mbinsize;
  }

  nentry = force->inumeric(FLERR,arg[4]);
  if (nentry < 1) return 1;
  delete [] mentry;
  mentry = new MapEntry[nentry];

  int expandflag = 0;

  int n = 5;
  for (int i = 0; i < nentry; i++) {
    if (mstyle == CONTINUOUS) {
      if (n+2 > narg) return 1;
      if (!islower(arg[n][0])) {
        mentry[i].single = NUMERIC;
        mentry[i].svalue = force->numeric(FLERR,arg[n]);
      } else if (strcmp(arg[n],"min") == 0) mentry[i].single = MINVALUE;
      else if (strcmp(arg[n],"max") == 0) mentry[i].single = MAXVALUE;
      else return 1;
      mentry[i].color = image->color2rgb(arg[n+1]);
      n += 2;
    } else if (mstyle == DISCRETE) {
      if (n+3 > narg) return 1;
      if (!islower(arg[n][0])) {
        mentry[i].lo = NUMERIC;
        mentry[i].lvalue = force->numeric(FLERR,arg[n]);
      } else if (strcmp(arg[n],"min") == 0) mentry[i].lo = MINVALUE;
      else if (strcmp(arg[n],"max") == 0) mentry[i].lo = MAXVALUE;
      else return 1;
      if (!islower(arg[n+1][0])) {
        mentry[i].hi = NUMERIC;
        mentry[i].hvalue = force->numeric(FLERR,arg[n+1]);
      } else if (strcmp(arg[n+1],"min") == 0) mentry[i].hi = MINVALUE;
      else if (strcmp(arg[n+1],"max") == 0) mentry[i].hi = MAXVALUE;
      else return 1;
      mentry[i].color = image->color2rgb(arg[n+2]);
      n += 3;
    } else if (mstyle == SEQUENTIAL) {
      // NOTE: this is unfinished code, not sure how useful it is
      // idea is to allow a list of colors to be specified with a single arg
      // problem is that sequential colors in ALL are not very different
      // e.g. ALL or USER or ALL5:10 or USER1:10:2
      // current code is just 1st nentry values of ALL or USER
      // need to comment out error check in DumpImage::modify_param()
      //   for amap check on (narg < n) to get it to work
      // need to add extra logic here to check not accessing undefined colors
      if (i == 0) {
        if (n+1 > narg) return 1;
        if (strcmp(arg[n],"ALL") == 0) expandflag = 1;
        if (strcmp(arg[n],"USER") == 0) expandflag = 2;
      }
      if (expandflag == 0) {
        if (n+1 > narg) return 1;
        mentry[i].color = image->color2rgb(arg[n]);
      } else if (expandflag == 1) {
        mentry[i].color = image->color2rgb(NULL,i+1);
      } else if (expandflag == 2) {
        mentry[i].color = image->color2rgb(NULL,-(i+1));
      }
      n += 1;
    }
    if (mentry[i].color == NULL) return 1;
  }

  if (mstyle == CONTINUOUS) {
    if (nentry < 2) return 1;
    if (mentry[0].single != MINVALUE || mentry[nentry-1].single != MAXVALUE)
      return 1;
    for (int i = 2; i < nentry-1; i++)
      if (mentry[i].svalue <= mentry[i-1].svalue) return 1;
  } else if (mstyle == DISCRETE) {
    if (nentry < 1) return 1;
    if (mentry[nentry-1].lo != MINVALUE || mentry[nentry-1].hi != MAXVALUE)
      return 1;
  } else if (mstyle == SEQUENTIAL) {
    if (nentry < 1) return 1;
  }

  // one-time call to minmax if color map is static

  if (!dynamic) return minmax(mlovalue,mhivalue);

  return 0;
}

/* ----------------------------------------------------------------------
   set explicit values for all min/max settings in color map
     from min/max dynamic values
   set lo/hi current and lvalue/hvalue entries that are MIN/MAX VALUE
   called only once if mlo/mhi != MIN/MAX VALUE, else called repeatedly
   return 1 = error if any values now overlap incorrectly with dynamic bounds
   else return 0
------------------------------------------------------------------------- */

int ColorMap::minmax(double mindynamic, double maxdynamic)
{
  if (mlo == MINVALUE) locurrent = mindynamic;
  else locurrent = mlovalue;
  if (mhi == MAXVALUE) hicurrent = maxdynamic;
  else hicurrent = mhivalue;
  if (locurrent > hicurrent) return 1;

  if (mstyle == CONTINUOUS) {
    if (mrange == ABSOLUTE) mentry[0].svalue = locurrent;
    else mentry[0].svalue = 0.0;
    if (mrange == ABSOLUTE) mentry[nentry-1].svalue = hicurrent;
    else mentry[nentry-1].svalue = 1.0;

    // error in ABSOLUTE mode if new lo/hi current cause
    // first/last entry to become lo > hi with adjacent entry

    if (mrange == ABSOLUTE) {
      if (mentry[0].svalue > mentry[1].svalue) return 1;
      if (mentry[nentry-2].svalue > mentry[nentry-1].svalue) return 1;
    }

  // OK if new lo/hi current cause an entry to have lo > hi,
  // since last entry will always be a match

  } else if (mstyle == DISCRETE) {
    for (int i = 0; i < nentry; i++) {
      if (mentry[i].lo == MINVALUE) {
        if (mrange == ABSOLUTE) mentry[i].lvalue = locurrent;
        else mentry[i].lvalue = 0.0;
      }
      if (mentry[i].hi == MAXVALUE) {
        if (mrange == ABSOLUTE) mentry[i].hvalue = hicurrent;
        else mentry[i].hvalue = 1.0;
      }
    }
  }

  return 0;
}

/* ----------------------------------------------------------------------
   convert value into an RGB color via color map
   return pointer to 3-vector
------------------------------------------------------------------------- */

double *ColorMap::value2color(double value)
{
  double lo;//,hi;

  value = MAX(value,locurrent);
  value = MIN(value,hicurrent);

  if (mrange == FRACTIONAL) {
    if (locurrent == hicurrent) value = 0.0;
    else value = (value-locurrent) / (hicurrent-locurrent);
    lo = 0.0;
    //hi = 1.0;
  } else {
    lo = locurrent;
    //hi = hicurrent;
  }

  if (mstyle == CONTINUOUS) {
    for (int i = 0; i < nentry-1; i++)
      if (value >= mentry[i].svalue && value <= mentry[i+1].svalue) {
        double fraction = (value-mentry[i].svalue) /
          (mentry[i+1].svalue-mentry[i].svalue);
        interpolate[0] = mentry[i].color[0] +
          fraction*(mentry[i+1].color[0]-mentry[i].color[0]);
        interpolate[1] = mentry[i].color[1] +
          fraction*(mentry[i+1].color[1]-mentry[i].color[1]);
        interpolate[2] = mentry[i].color[2] +
          fraction*(mentry[i+1].color[2]-mentry[i].color[2]);
        return interpolate;
      }
  } else if (mstyle == DISCRETE) {
    for (int i = 0; i < nentry; i++)
      if (value >= mentry[i].lvalue && value <= mentry[i].hvalue)
        return mentry[i].color;
  } else {
    int ibin = static_cast<int> ((value-lo) * mbinsizeinv);
    return mentry[ibin%nentry].color;
  }

  return NULL;
}