ObjectMap.c

/* 
A* -------------------------------------------------------------------
B* This file contains source code for the PyMOL computer program
C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific. 
D* -------------------------------------------------------------------
E* It is unlawful to modify or remove this copyright notice.
F* -------------------------------------------------------------------
G* Please see the accompanying LICENSE file for further information. 
H* -------------------------------------------------------------------
I* Additional authors of this source file include:
-* 
-* 
-*
Z* -------------------------------------------------------------------
*/

#include"os_predef.h"
#include"os_std.h"
#include"os_gl.h"

#include"OOMac.h"
#include"ObjectMap.h"
#include"Base.h"
#include"MemoryDebug.h"
#include"Map.h"
#include"Parse.h"
#include"Isosurf.h"
#include"Vector.h"
#include"Color.h"
#include"main.h"
#include"Scene.h"
#include"PConv.h"
#include"Word.h"
#include"Vector.h"

#define cMapSourceUndefined 0
#define cMapSourceXPLOR 1
#define cMapSourceCCP4 2
#define cMapSourcePHI 3
#define cMapSourceDesc 4
#define cMapSourceFLD 5
#define cMapSourceBRIX 6
#define cMapSourceGRD 7

#ifdef _PYMOL_NUMPY
typedef struct {
  PyObject_HEAD
  char *data;
  int nd;
  int *dimensions, *strides;
  PyObject *base;
  void *descr;
  int flags;
} MyArrayObject;
#endif

int ObjectMapStateGetExcludedStats(ObjectMapState *ms,float *vert_vla, float beyond,float within, float *level)
{
  double sum=0.0,sumsq=0.0;
  float mean,stdev;  
  int cnt = 0;
  int list_size = VLAGetSize(vert_vla)/3;
  float cutoff = beyond;
  MapType *voxelmap = NULL;

  if(cutoff<within)
    cutoff = within;

  if(list_size) 
    voxelmap=MapNew(-cutoff,vert_vla,list_size,NULL);
  if(voxelmap||(!list_size)) {
    int a,b,c;
    int h,k,l,i,j;
    int *fdim = ms->FDim;
    float *v,f_val;
    int within_flag, within_default=false;
    int beyond_flag;
    
    Isofield *field = ms->Field;

    if(list_size)
      MapSetupExpress(voxelmap);  

    within_flag=true;
    beyond_flag=true;

    if(within<R_SMALL4)
      within_default = true;
    for(c=0;c<fdim[2];c++) {
      for(b=0;b<fdim[1];b++) {
        for(a=0;a<fdim[0];a++) {
          if(list_size) {
            within_flag = within_default;
            beyond_flag = true;

            v = F4Ptr(field->points,a,b,c,0);
            
            MapLocus(voxelmap,v,&h,&k,&l);
            i=*(MapEStart(voxelmap,h,k,l));
            if(i) {
              j=voxelmap->EList[i++];
              while(j>=0) {
                if(!within_flag) {
                  if(within3f(vert_vla+3*j,v,within)) {                  
                    within_flag=true;
                  }
                }
                if(within3f(vert_vla+3*j,v,beyond)) {
                  beyond_flag=false;
                  break;
                }
                j=voxelmap->EList[i++];
              }
            }
          }

          if(within_flag&&beyond_flag) { /* point isn't too close to any vertex */
            f_val = F3(field->data,a,b,c);
            sum+=f_val;
            sumsq+=(f_val*f_val);
            cnt++;
          }
        }
      }
    }
    if(voxelmap) 
      MapFree(voxelmap);
  }
  if(cnt) {
    mean = (float)(sum/cnt);
    stdev = (float)sqrt1d((sumsq - (sum*sum/cnt))/(cnt));
    level[1] = mean;
    level[0] = mean-stdev;
    level[2] = mean+stdev;
  }
  return cnt;
}

int ObjectMapInterpolate(ObjectMap *I,int state,float *array,float *result,int n)
{
  int ok=false;
  if(state<0) state=0;
  if(state<I->NState)
    if(I->State[state].Active)
      ok = ObjectMapStateInterpolate(&I->State[state],array,result,n);
  return(ok);
}

static int ObjectMapStateDouble(ObjectMapState *ms)
{
  int div[3];
  int min[3];
  int max[3];
  int fdim[4];
  int a,b,c;
  float v[3],vr[3];
  float *vt;
  float x,y,z;
  float grid[3];

  Isofield *field;

  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    for(a=0;a<3;a++) {
      div[a]=ms->Div[a]*2;
      min[a]=ms->Min[a]*2;
      max[a]=ms->Max[a]*2;
      fdim[a]=ms->FDim[a]*2-1;
    }
    fdim[3]=3;

    field=IsosurfFieldAlloc(fdim);

    for(c=0;c<fdim[2];c++) {
      v[2]=(c+min[2])/((float)div[2]);
      z = (c&0x1) ? 0.5F : 0.0F;
      for(b=0;b<fdim[1];b++) {
        v[1]=(b+min[1])/((float)div[1]);        
        y = (b&0x1) ? 0.5F : 0.0F;
        for(a=0;a<fdim[0];a++) {
          v[0]=(a+min[0])/((float)div[0]);
          transform33f3f(ms->Crystal->FracToReal,v,vr);
          x = (a&0x1) ? 0.5F : 0.0F;
          vt = F4Ptr(field->points,a,b,c,0);
          copy3f(vr,vt);
          if((a&0x1)||(b&0x1)||(c&0x1)) {
            F3(field->data,a,b,c) = FieldInterpolatef(ms->Field->data,
                                    a/2,
                                    b/2,
                                    c/2,x,y,z);
          } else {
            F3(field->data,a,b,c) = F3(ms->Field->data,a/2,b/2,c/2);
          }
        }
      }
    }
    IsosurfFieldFree(ms->Field);
    for(a=0;a<3;a++) {
      ms->Min[a]=min[a];
      ms->Max[a]=max[a];
      ms->FDim[a]=fdim[a];
      ms->Div[a]=div[a];
    }
    ms->Field = field;
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:

    for(a=0;a<3;a++) {
      div[a]=ms->Div[a]*2;
      grid[a]=ms->Grid[a]/2.0F;
      min[a]=ms->Min[a]*2;
      max[a]=ms->Max[a]*2;
      fdim[a]=ms->FDim[a]*2-1;
    }
    fdim[3]=3;

    field=IsosurfFieldAlloc(fdim);

    for(c=0;c<fdim[2];c++) {
      v[2]=ms->Origin[2]+grid[2]*(c+min[2]);
      z = (c&0x1) ? 0.5F : 0.0F;
      for(b=0;b<fdim[1];b++) {
        v[1]=ms->Origin[1]+grid[1]*(b+min[1]);
        y = (b&0x1) ? 0.5F : 0.0F;
        for(a=0;a<fdim[0];a++) {
          v[0]=ms->Origin[0]+grid[0]*(a+min[0]);
          x = (a&0x1) ? 0.5F : 0.0F;
          vt = F4Ptr(field->points,a,b,c,0);
          copy3f(v,vt);
          if((a&0x1)||(b&0x1)||(c&0x1)) {
            F3(field->data,a,b,c) = FieldInterpolatef(ms->Field->data,
                                    a/2,
                                    b/2,
                                    c/2,x,y,z);
          } else {
            F3(field->data,a,b,c) = F3(ms->Field->data,a/2,b/2,c/2);
          }
        }
      }
    }
    IsosurfFieldFree(ms->Field);
    for(a=0;a<3;a++) {
      ms->Min[a]=min[a];
      ms->Max[a]=max[a];
      ms->FDim[a]=fdim[a];
      if(ms->Dim) 
        ms->Dim[a]=fdim[a];
      ms->Div[a]=div[a];
      if(ms->Grid)
        ms->Grid[a]=grid[a];
    }
    ms->Field = field;

    break;
  }
  return 1;
}

int ObjectMapDouble(ObjectMap *I,int state)
{
  int a;
  int result=true;
  if(state<0) {
    for(a=0;a<I->NState;a++) {
      if(I->State[a].Active)
        result = result && ObjectMapStateDouble(&I->State[a]);
    }
  } else if((state>=0)&&(state<I->NState)&&(I->State[state].Active)) {
    ObjectMapStateDouble(&I->State[state]);
  } else {
    PRINTFB(FB_ObjectMap,FB_Errors)
      " ObjectMap-Error: invalidate state.\n"
      ENDFB;
    result=false;
  }
  return(result);
}


int ObjectMapStateInterpolate(ObjectMapState *ms,float *array,float *result,int n)
{
  int ok=true;
  float *inp;
  float *out;
  int a,b,c;
  float x,y,z;
  inp = array;
  out = result;
  
  switch(ms->MapSource) {
  case cMapSourcePHI:
  case cMapSourceFLD:
  case cMapSourceDesc:
    while(n--) {

      x = (inp[0] - ms->Origin[0])/ms->Grid[0];
      y = (inp[1] - ms->Origin[1])/ms->Grid[1];
      z = (inp[2] - ms->Origin[2])/ms->Grid[2];
      inp+=3;
      
      a=(int)floor(x);
      b=(int)floor(y);
      c=(int)floor(z);
      x-=a;
      y-=b;
      z-=c;

      if(a<ms->Min[0]) {
        x=0.0F;
        a=ms->Min[0];
        ok=false;
      } else if(a>=ms->Max[0]) {
        x=1.0F;
        a=ms->Max[0]-1;
        ok=false;
      }

      if(b<ms->Min[1]) {
        y=0.0F;
        b=ms->Min[1];
        ok=false;
      } else if(b>=ms->Max[1]) {
        y=1.0F;
        b=ms->Min[1];
        ok=false;
      }

      if(c<ms->Min[2]) {
        z=0.0F;
        c=ms->Min[2];
        ok=false;
      } else if(c>=ms->Max[2]) {
        z=1.0F;
        c=ms->Max[2]-1;
        ok=false;
      }
      /*      printf("%d %d %d %8.3f %8.3f %8.3f\n",a,b,c,x,y,z);*/
      *(result++)=FieldInterpolatef(ms->Field->data,
                                    a-ms->Min[0],
                                    b-ms->Min[1],
                                    c-ms->Min[2],x,y,z);
    }
    break;
  default:
    ok=false;
    break;
  }
  return(ok);
}

int ObjectMapNumPyArrayToMapState(ObjectMapState *I,PyObject *ary);

static void ObjectMapStateRegeneratePoints(ObjectMapState *ms)
{
  int a,b,c,e;
  float v[3],vr[3];
  switch(ms->MapSource) {
  case cMapSourceXPLOR:
  case cMapSourceCCP4:
  case cMapSourceBRIX:
  case cMapSourceGRD:
    for(c=0;c<ms->FDim[2];c++)
      {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b++) {
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a++) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            for(e=0;e<3;e++) 
              F4(ms->Field->points,a,b,c,e) = vr[e];
          }
        }
      }
    break;
  case cMapSourcePHI:
  case cMapSourceFLD:
    for(c=0;c<ms->FDim[2];c++) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      for(b=0;b<ms->FDim[1];b++) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        for(a=0;a<ms->FDim[0];a++) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          for(e=0;e<3;e++) {
            F4(ms->Field->points,a,b,c,e) = v[e];
          }
        }
      }
  }
  default: 
    break;
  }
}

static PyObject *ObjectMapStateAsPyList(ObjectMapState *I)
{
  PyObject *result = NULL;

  result = PyList_New(15);
  PyList_SetItem(result,0,PyInt_FromLong(I->Active));
  if(I->Crystal) {
    PyList_SetItem(result,1,CrystalAsPyList(I->Crystal));
  } else {
    PyList_SetItem(result,1,PConvAutoNone(Py_None));
  }
  if(I->Origin) {
    PyList_SetItem(result,2,PConvFloatArrayToPyList(I->Origin,3));
  } else {
    PyList_SetItem(result,2,PConvAutoNone(Py_None));
  }
  if(I->Range) {
    PyList_SetItem(result,3,PConvFloatArrayToPyList(I->Range,3));
  } else {
    PyList_SetItem(result,3,PConvAutoNone(Py_None));
  }
  if(I->Dim) {
    PyList_SetItem(result,4,PConvIntArrayToPyList(I->Dim,3));
  } else {
    PyList_SetItem(result,4,PConvAutoNone(Py_None));
  }
  if(I->Grid) {
    PyList_SetItem(result,5,PConvFloatArrayToPyList(I->Grid,3));
  } else {
    PyList_SetItem(result,5,PConvAutoNone(Py_None));
  }
  PyList_SetItem(result,6,PConvFloatArrayToPyList(&I->Corner[0][0],24));
  PyList_SetItem(result,7,PConvFloatArrayToPyList(I->ExtentMin,3));
  PyList_SetItem(result,8,PConvFloatArrayToPyList(I->ExtentMax,3));
  PyList_SetItem(result,9,PyInt_FromLong(I->MapSource));

  PyList_SetItem(result,10,PConvIntArrayToPyList(I->Div,3));
  PyList_SetItem(result,11,PConvIntArrayToPyList(I->Min,3));
  PyList_SetItem(result,12,PConvIntArrayToPyList(I->Max,3));
  PyList_SetItem(result,13,PConvIntArrayToPyList(I->FDim,4));
  
  PyList_SetItem(result,14,IsosurfAsPyList(I->Field));
#if 0
  int Active;
  CCrystal *Crystal;
  int Div[3],Min[3],Max[3],FDim[4];
  Isofield *Field;
  float Corner[8][3];
  int *Dim;
  float *Origin;
  float *Range;
  float *Grid;
  float ExtentMin[3],ExtentMax[3];
#endif

  return(PConvAutoNone(result));  
}

static PyObject *ObjectMapAllStatesAsPyList(ObjectMap *I)
{
  PyObject *result=NULL;
  int a;
  result = PyList_New(I->NState);
  for(a=0;a<I->NState;a++) {
    if(I->State[a].Active) {
      PyList_SetItem(result,a,ObjectMapStateAsPyList(I->State+a));
    } else {
      PyList_SetItem(result,a,PConvAutoNone(NULL));
    }
  }
  return(PConvAutoNone(result));  
}

static int ObjectMapStateFromPyList(ObjectMapState *I,PyObject *list)
{
  int ok=true;
  int ll;
  PyObject *tmp;
  if(ok) ok=(list!=NULL);
  if(ok) {
    if(!PyList_Check(list))
      I->Active=false;
    else {
      if(ok) ok=PyList_Check(list);
      if(ok) ll = PyList_Size(list);
      /* TO SUPPORT BACKWARDS COMPATIBILITY...
         Always check ll when adding new PyList_GetItem's */

      if(ok) ok = PConvPyIntToInt(PyList_GetItem(list,0),&I->Active);
      if(ok) {
        tmp = PyList_GetItem(list,1);
        if(tmp == Py_None)
          I->Crystal = NULL;
        else 
          ok = ((I->Crystal=CrystalNewFromPyList(tmp))!=NULL);
      }
      if(ok) {
        tmp = PyList_GetItem(list,2);
        if(tmp == Py_None)
          I->Origin = NULL;
        else 
          ok = PConvPyListToFloatArray(tmp,&I->Origin);
      }
      if(ok) {
        tmp = PyList_GetItem(list,3);
        if(tmp == Py_None)
          I->Range = NULL;
        else 
          ok = PConvPyListToFloatArray(tmp,&I->Range);
      }
      if(ok) {
        tmp = PyList_GetItem(list,4);
        if(tmp == Py_None)
          I->Dim = NULL;
        else 
          ok = PConvPyListToIntArray(tmp,&I->Dim);
      }
      if(ok) {
        tmp = PyList_GetItem(list,5);
        if(tmp == Py_None)
          I->Grid = NULL;
        else 
          ok = PConvPyListToFloatArray(tmp,&I->Grid);
      }
      if(ok) ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list,6),&I->Corner[0][0],24);
      if(ok) ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list,7),I->ExtentMin,3);
      if(ok) ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list,8),I->ExtentMax,3);
      if(ok) ok = PConvPyIntToInt(PyList_GetItem(list,9),&I->MapSource);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,10),I->Div,3);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,11),I->Min,3);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,12),I->Max,3);
      if(ok) ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list,13),I->FDim,4);
      if(ok) ok = ((I->Field=IsosurfNewFromPyList(PyList_GetItem(list,14)))!=NULL);
      if(ok) ObjectMapStateRegeneratePoints(I);
    }
  }
  return(ok);
}

static int ObjectMapAllStatesFromPyList(ObjectMap *I,PyObject *list)
{
  int ok=true;
  int a;
  VLACheck(I->State,ObjectMapState,I->NState);
  if(ok) ok=PyList_Check(list);
  if(ok) {
    for(a=0;a<I->NState;a++) {
      ok = ObjectMapStateFromPyList(I->State+a,PyList_GetItem(list,a));
      if(!ok) break;
    }
  }
  return(ok);
}


PyObject *ObjectMapAsPyList(ObjectMap *I)
{
  PyObject *result = NULL;

  result = PyList_New(3);
  PyList_SetItem(result,0,ObjectAsPyList(&I->Obj));
  PyList_SetItem(result,1,PyInt_FromLong(I->NState));
  PyList_SetItem(result,2,ObjectMapAllStatesAsPyList(I));

  return(PConvAutoNone(result));  
}

int ObjectMapNewFromPyList(PyObject *list,ObjectMap **result)
{
  int ok = true;
  int ll;
  ObjectMap *I=NULL;
  (*result) = NULL;
  
  if(ok) ok=(list!=NULL);
  if(ok) ok=PyList_Check(list);
  if(ok) ll = PyList_Size(list);
  /* TO SUPPORT BACKWARDS COMPATIBILITY...
   Always check ll when adding new PyList_GetItem's */
  I=ObjectMapNew();
  if(ok) ok = (I!=NULL);

  if(ok) ok = ObjectFromPyList(PyList_GetItem(list,0),&I->Obj);
  if(ok) ok = PConvPyIntToInt(PyList_GetItem(list,1),&I->NState);
  if(ok) ok = ObjectMapAllStatesFromPyList(I,PyList_GetItem(list,2));
  if(ok) {
    (*result) = I;
    ObjectMapUpdateExtents(I);
  } else {
    /* cleanup? */
  }

  return(ok);
}


ObjectMapState *ObjectMapGetState(ObjectMap *I,int state)
{
  ObjectMapState *result = NULL;

  if(state<I->NState)
    result = &I->State[state];
  return(result);
}

ObjectMapState *ObjectMapStatePrime(ObjectMap *I,int state)
{
  ObjectMapState *ms = NULL;
  if(state<0)
    state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);
  return(ms);
}

ObjectMapState *ObjectMapStateGetActive(ObjectMap *I,int state)
{
  ObjectMapState *ms = NULL;
  if(state>=0) {
    if(state<I->NState) {
      ms=&I->State[state];
      if(!ms->Active)
        ms = NULL;
    }
  }
  return(ms);
}


void ObjectMapUpdateExtents(ObjectMap *I)
{
  int a;
  I->Obj.ExtentFlag=false;
  for(a=0;a<I->NState;a++) {
    if(I->State[a].Active)
      {
        if(!I->Obj.ExtentFlag) {
          copy3f(I->State[a].ExtentMin,I->Obj.ExtentMin);
          copy3f(I->State[a].ExtentMax,I->Obj.ExtentMax);
          I->Obj.ExtentFlag=true;
        } else {
          min3f(I->State[a].ExtentMin,I->Obj.ExtentMin,I->Obj.ExtentMin);
          max3f(I->State[a].ExtentMax,I->Obj.ExtentMax,I->Obj.ExtentMax);
        }
      }
  }
  PRINTFD(FB_ObjectMap)
    " ObjectMapUpdateExtents-DEBUG: ExtentFlag %d\n",I->Obj.ExtentFlag
    ENDFD;
}

int ObjectMapStateSetBorder(ObjectMapState *I,float level)
{
  int result = false;
  int a,b,c;

  c=I->FDim[2]-1;
  for(a=0;a<I->FDim[0];a++) 
    for(b=0;b<I->FDim[1];b++)
      {
        F3(I->Field->data,a,b,0) = level;
        F3(I->Field->data,a,b,c) = level;
      }

  a=I->FDim[0]-1;
  for(b=0;b<I->FDim[1];b++) 
    for(c=0;c<I->FDim[2];c++)
      {
        F3(I->Field->data,0,b,c) = level;
        F3(I->Field->data,a,b,c) = level;
      }

  b=I->FDim[1]-1;
  for(a=0;a<I->FDim[0];a++) 
    for(c=0;c<I->FDim[2];c++)
      {
        F3(I->Field->data,a,0,c) = level;
        F3(I->Field->data,a,b,c) = level;
      }
  return(result);
}

void ObjectMapStatePurge(ObjectMapState *I)
{
  if(I->Field) {
    IsosurfFieldFree(I->Field);
    I->Field=NULL;
  }
  FreeP(I->Origin);
  FreeP(I->Dim);
  FreeP(I->Range);
  FreeP(I->Grid);
  OOFreeP(I->Crystal);
  I->Active=false;
}

static void ObjectMapFree(ObjectMap *I) {

  int a;
  for(a=0;a<I->NState;a++) {
    if(I->State[a].Active)
      ObjectMapStatePurge(I->State+a);
  }
  VLAFreeP(I->State);
  ObjectPurge(&I->Obj);
  OOFreeP(I);
}

static void ObjectMapUpdate(ObjectMap *I) {
  SceneDirty();
}

static void ObjectMapRender(ObjectMap *I,int state,CRay *ray,Pickable **pick,int pass)
{
  ObjectMapState *ms = NULL;
  if(!pass) {
    if(state<I->NState)
      if(I->State[state].Active)
        ms=&I->State[state];
    
    if(ms) {
      ObjectPrepareContext(&I->Obj,ray);

      if(I->Obj.RepVis[cRepExtent]) {
        if(ray) {
          float *vc;
          vc = ColorGet(I->Obj.Color);
          ray->fColor3fv(ray,vc);
          ray->fSausage3fv(ray,ms->Corner[0],ms->Corner[1],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[0],ms->Corner[2],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[2],ms->Corner[3],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[1],ms->Corner[3],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[0],ms->Corner[4],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[1],ms->Corner[5],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[2],ms->Corner[6],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[3],ms->Corner[7],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[4],ms->Corner[5],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[4],ms->Corner[6],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[6],ms->Corner[7],0.20F,vc,vc);
          ray->fSausage3fv(ray,ms->Corner[5],ms->Corner[7],0.20F,vc,vc);
        } else if(pick&&PMGUI) {
        } else if(PMGUI) {
          ObjectUseColor(&I->Obj);
          glDisable(GL_LIGHTING); 
          glBegin(GL_LINES);
          glVertex3fv(ms->Corner[0]);
          glVertex3fv(ms->Corner[1]);
          
          glVertex3fv(ms->Corner[0]);
          glVertex3fv(ms->Corner[2]);
          
          glVertex3fv(ms->Corner[2]);
          glVertex3fv(ms->Corner[3]);
          
          glVertex3fv(ms->Corner[1]);
          glVertex3fv(ms->Corner[3]);
          
          glVertex3fv(ms->Corner[0]);
          glVertex3fv(ms->Corner[4]);
          
          glVertex3fv(ms->Corner[1]);
          glVertex3fv(ms->Corner[5]);
          
          glVertex3fv(ms->Corner[2]);
          glVertex3fv(ms->Corner[6]);
          
          glVertex3fv(ms->Corner[3]);
          glVertex3fv(ms->Corner[7]);
          
          glVertex3fv(ms->Corner[4]);
          glVertex3fv(ms->Corner[5]);
          
          glVertex3fv(ms->Corner[4]);
          glVertex3fv(ms->Corner[6]);
          
          glVertex3fv(ms->Corner[6]);
          glVertex3fv(ms->Corner[7]);
          
          glVertex3fv(ms->Corner[5]);
          glVertex3fv(ms->Corner[7]);
          
          glEnd();
          glEnable(GL_LIGHTING);
        }
      }
    }
  }
}

void ObjectMapStateInit(ObjectMapState *I) 
{
  ObjectMapStatePurge(I);
  I->Crystal = CrystalNew();
  I->Field = NULL;
  I->Origin = NULL;
  I->Dim = NULL;
  I->Range = NULL;
  I->Grid = NULL;
}
int ObjectMapGetNStates(ObjectMap *I)     
{
  return(I->NState);
}
/*========================================================================*/
ObjectMap *ObjectMapNew(void)
{
  OOAlloc(ObjectMap);

  ObjectInit((CObject*)I);
  I->Obj.type = cObjectMap;

  
  I->NState = 0;
  I->State=VLAMalloc(1,sizeof(ObjectMapState),5,true); /* autozero important */

  I->Obj.RepVis[cRepExtent]=true; 
  I->Obj.fFree = (void (*)(struct CObject *))ObjectMapFree;
  I->Obj.fUpdate =  (void (*)(struct CObject *)) ObjectMapUpdate;
  I->Obj.fRender =(void (*)(struct CObject *, int, CRay *, Pickable **,int))ObjectMapRender;

  I->Obj.fGetNFrame = (int (*)(struct CObject *)) ObjectMapGetNStates;

  return(I);
}
/*========================================================================*/
ObjectMapState *ObjectMapNewStateFromDesc(ObjectMap *I,ObjectMapDesc *md,int state)
{
  int ok=true;
  float v[3];
  int a,b,c,d;
  float *fp;
  ObjectMapState *ms = NULL;

  ms=ObjectMapStatePrime(I,state);
  
  if(I) {
    ms->Origin=Alloc(float,3);
    ms->Range=Alloc(float,3);
    ms->Dim=Alloc(int,3);
    ms->Grid=Alloc(float,3);
    ms->MapSource=cMapSourceDesc;
  }
  switch(md->mode) {
  case cObjectMap_OrthoMinMaxGrid: /* Orthorhombic: min, max, spacing, centered over range  */

    subtract3f(md->MaxCorner,md->MinCorner,v);
    for(a=0;a<3;a++) { if(v[a]<0.0) swap1f(md->MaxCorner+a,md->MinCorner+a); };
    subtract3f(md->MaxCorner,md->MinCorner,v);
    for(a=0;a<3;a++) {
      md->Dim[a] = (int)(v[a]/md->Grid[a]);
      if(md->Dim[a]<1) md->Dim[a]=1;
      if((md->Dim[a]*md->Grid[a])<v[a]) md->Dim[a]++;
    }

    PRINTFB(FB_ObjectMap,FB_Blather)
      " ObjectMap: Dim %d %d %d\n",md->Dim[0],md->Dim[1],md->Dim[2]
      ENDFB;

    average3f(md->MaxCorner,md->MinCorner,v);
    for(a=0;a<3;a++) { md->MinCorner[a] = v[a]-0.5F*md->Dim[a]*md->Grid[a]; }

    if(Feedback(FB_ObjectMap,FB_Blather)) {
      dump3f(md->MinCorner," ObjectMap: MinCorner:");
      dump3f(md->MaxCorner," ObjectMap: MaxCorner:");
      dump3f(md->Grid," ObjectMap: Grid:");
    }
    
    /* now populate the map data structure */

    copy3f(md->MinCorner,ms->Origin);
    copy3f(md->Grid,ms->Grid);
    for(a=0;a<3;a++) ms->Range[a] = md->Grid[a] * (md->Dim[a]-1);

    /* these maps start at zero */
    for(a=0;a<3;a++) {
      ms->Min[a]=0; 
      ms->Max[a]=md->Dim[a]-1;
    }
    
    /* define corners */

    for(a=0;a<8;a++) copy3f(ms->Origin,ms->Corner[a]);

    d = 0;
    for(c=0;c<2;c++) {
      {
        v[2] = (c ? ms->Range[2] : 0.0F);
        for(b=0;b<2;b++) {
          v[1]= (b ? ms->Range[1] : 0.0F);
          for(a=0;a<2;a++) {
            v[0]= (a ? ms->Range[0] : 0.0F);
            add3f(v,ms->Corner[d],ms->Corner[d]);
            d++;
          }
        }
      }
    }
    for(a=0;a<3;a++) ms->FDim[a] = ms->Max[a]-ms->Min[a]+1;
    ms->FDim[3] = 3; 

    ms->Field=IsosurfFieldAlloc(ms->FDim);
    if(!ms->Field) 
      ok=false;
    else {
      for(a=0;a<md->Dim[0];a++) {
        v[0] = md->MinCorner[0] + a * md->Grid[0];
        for(b=0;b<md->Dim[1];b++) {
          v[1] = md->MinCorner[1] + b * md->Grid[1];
          for(c=0;c<md->Dim[2];c++) {
            v[2] = md->MinCorner[2] + c * md->Grid[2];
            fp = F4Ptr(ms->Field->points,a,b,c,0);
            copy3f(v,fp);
          }
        }
      }
    }
    break;
  default:
    ok = false;
  }
  if(ok) {
    switch(md->init_mode) {
    case 0:
      for(a=0;a<md->Dim[0];a++) {      
        for(b=0;b<md->Dim[1];b++) {
          for(c=0;c<md->Dim[2];c++) {
            F3(ms->Field->data,a,b,c)=0.0F;
          }
        }
      }
      break;
    case 1:
      for(a=0;a<md->Dim[0];a++) {      
        for(b=0;b<md->Dim[1];b++) {
          for(c=0;c<md->Dim[2];c++) {
            F3(ms->Field->data,a,b,c)=1.0F;
          }
        }
      }
      break;
    case -2: /* for testing... */
      for(a=0;a<md->Dim[0];a++) {      
        for(b=0;b<md->Dim[1];b++) {
          for(c=0;c<md->Dim[2];c++) {
            F3(ms->Field->data,a,b,c)=(float)sqrt1d(a*a+b*b+c*c);
          }
        }
      }
      break;
    }
  }
  
  if(ok) {
    copy3f(ms->Origin,ms->ExtentMin);
    copy3f(ms->Origin,ms->ExtentMax);
    add3f(ms->Range,ms->ExtentMax,ms->ExtentMax);
    ObjectMapUpdateExtents(I);
  }
  if(!ok) {
    ErrMessage("ObjectMap","Unable to create map");
    ObjectMapFree(I);
    I=NULL;
  } else {
    PRINTFB(FB_ObjectMap,FB_Actions) 
      " ObjectMap: Map created.\n"
      ENDFB;
  }
  
  return(ms);
}
/*========================================================================*/
int ObjectMapCCP4StrToMap(ObjectMap *I,char *CCP4Str,int bytes,int state) {
  
  char *p;
  int *i;
  unsigned int *u;
  unsigned char *uc,c0,c1,c2,c3;
  float *f;
  float dens;
  int a,b,c,d,e;
  float v[3],vr[3],maxd,mind;
  int ok = true;
  int little_endian = 1,map_endian;
  /* CCP4 named from their docs */
  int nc,nr,ns;
  int map_mode;
  int ncstart,nrstart,nsstart;
  int nx,ny,nz;
  float xlen,ylen,zlen,alpha,beta,gamma;
  int n_skip;
  int mapc,mapr,maps;
  int cc[3],xref[3];
  int n_pts;
  double sum,sumsq;
  float mean,stdev;
  int normalize;
  ObjectMapState *ms;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);

  normalize=(int)SettingGet(cSetting_normalize_ccp4_maps);
  maxd = FLT_MIN;
  mind = FLT_MAX;
  p=CCP4Str;
  little_endian = *((char*)&little_endian);
  map_endian = (*p||*(p+1));

  if(bytes<256*sizeof(int)) {
    PRINTFB(FB_ObjectMap,FB_Errors)
      " ObjectMapCCP4: Map appears to be truncated -- aborting."
      ENDFB;
    return(0);
  }
  if(little_endian!=map_endian) {
    PRINTFB(FB_ObjectMap,FB_Blather)
      " ObjectMapCCP4: Map appears to be reverse endian, swapping...\n"
      ENDFB;
    c = bytes;
    u = (unsigned int*)p;
    uc = (unsigned char *)u;
    while(c>3) {
      c0 = *(uc++);
      c1 = *(uc++);
      c2 = *(uc++);
      c3 = *(uc++);
      uc = (unsigned char *)u;
      *(uc++) = c3;
      *(uc++) = c2;
      *(uc++) = c1;
      *(uc++) = c0;
      u++;
      c-=4;
    }
  }
  i = (int*)p;
  nc=*(i++); /* columns */
  nr=*(i++); /* rows */
  ns=*(i++); /* sections */
  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: NC %d   NR %d   NS %d\n",
    nc,nr,ns
    ENDFB;
  map_mode = *(i++); /* mode */

  
  if(map_mode!=2) {
    PRINTFB(FB_ObjectMap,FB_Errors)
      "ObjectMapCCP4-ERR: Only map mode 2 currently supported (this map is mode %d)",map_mode
      ENDFB;
    return(0);
  }

  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: Map is mode %d.\n",map_mode
    ENDFB;

  ncstart = *(i++);
  nrstart = *(i++);
  nsstart = *(i++);

  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: NCSTART %d   NRSTART %d   NSSTART  %d\n",
    ncstart,nrstart,nsstart
    ENDFB;

  nx = *(i++);
  ny = *(i++);
  nz = *(i++);

  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: NX %d   NY %d   NZ  %d \n",
    nx,ny,nz
    ENDFB;

  xlen = *(float*)(i++);
  ylen = *(float*)(i++);
  zlen = *(float*)(i++);

  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: X %8.3f   Y %8.3f  Z  %8.3f \n",
    xlen,ylen,zlen
    ENDFB;

  alpha = *(float*)(i++);
  beta = *(float*)(i++);
  gamma = *(float*)(i++);

  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: alpha %8.3f   beta %8.3f  gamma %8.3f \n",
    alpha,beta,gamma
    ENDFB;

  mapc = *(i++);
  mapr = *(i++);
  maps = *(i++);

  PRINTFB(FB_ObjectMap,FB_Blather)
    " ObjectMapCCP4: MAPC %d   MAPR %d  MAPS  %d \n",
    mapc,mapr,maps
    ENDFB;

  i+=4;
  n_skip = *(i++);
  
  if(*(i++)) {
    PRINTFB(FB_ObjectMap,FB_Errors)
      "ObjectMapCCP4-ERR: PyMOL doesn't know how to handle skewed maps. Sorry!\n"
      ENDFB;
    return(0);
  }

  n_pts = nc*ns*nr;
  if((unsigned)bytes<(n_skip + sizeof(int)*(256+n_pts))) {
    PRINTFB(FB_ObjectMap,FB_Errors)
      " ObjectMapCCP4: Map appears to be truncated -- aborting.\n"
      ENDFB;
    return(0);
  }


  if(n_pts>1) {
    f = (float*)(p+(sizeof(int)*256)+n_skip);
    c = n_pts;
    sum = 0.0;
    sumsq = 0.0;
    while(c--) {
      sumsq+=(*f)*(*f);
      sum+=*f++;
    }
    mean = (float)(sum/n_pts);
    stdev = (float)sqrt1d((sumsq - (sum*sum/n_pts))/(n_pts-1));

    if(normalize) {
      PRINTFB(FB_ObjectMap,FB_Details)
        " ObjectMapCCP4: Normalizing with mean = %8.6f and stdev = %8.6f.\n",
        mean,stdev
        ENDFB;
    } else {
      PRINTFB(FB_ObjectMap,FB_Details)
        " ObjectMapCCP4: Map will not be normalized.\n ObjectMapCCP4: Current mean = %8.6f and stdev = %8.6f.\n",
        mean,stdev
        ENDFB;
    }

    if(stdev<0.000001)
      stdev = 1.0;
    
  } else {
    mean = 1.0;
    stdev = 1.0;
  }
  
  f = (float*)(p+(sizeof(int)*256)+n_skip);
  mapc--; /* convert to C indexing... */
  mapr--;
  maps--;

  xref[maps]=0;
  xref[mapr]=1;
  xref[maps]=2;

  xref[maps]=0;
  xref[mapr]=1;
  xref[mapc]=2;

  ms->Div[0] = nx;
  ms->Div[1] = ny;
  ms->Div[2] = nz;

  ms->FDim[mapc] = nc;
  ms->Min[mapc] = ncstart;
  ms->Max[mapc] = nc+ncstart-1;

  ms->FDim[mapr] = nr;
  ms->Min[mapr] = nrstart;
  ms->Max[mapr] = nr+nrstart-1;

  ms->FDim[maps] = ns;
  ms->Min[maps] = nsstart;
  ms->Max[maps] = ns+nsstart-1;

  if(Feedback(FB_ObjectMap,FB_Blather)) {
    dump3i(ms->Div,"ms->Div");
    dump3i(ms->Min,"ms->Min");
    dump3i(ms->Max,"ms->Max");
    dump3i(ms->FDim,"ms->FDim");
  }

  ms->Crystal->Dim[0] = xlen;
  ms->Crystal->Dim[1] = ylen;
  ms->Crystal->Dim[2] = zlen;

  ms->Crystal->Angle[0] = alpha;
  ms->Crystal->Angle[1] = beta;
  ms->Crystal->Angle[2] = gamma;


  ms->FDim[3]=3;
  if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
    ok=false;
  else {
    CrystalUpdate(ms->Crystal);
    CrystalDump(ms->Crystal);
    fflush(stdout);
    ms->Field=IsosurfFieldAlloc(ms->FDim);
    ms->MapSource = cMapSourceCCP4;
    ms->Field->save_points=false;

    for(cc[maps]=0;cc[maps]<ms->FDim[maps];cc[maps]++)
      {
        v[maps]=(cc[maps]+ms->Min[maps])/((float)ms->Div[maps]);

        for(cc[mapr]=0;cc[mapr]<ms->FDim[mapr];cc[mapr]++) {
          v[mapr]=(cc[mapr]+ms->Min[mapr])/((float)ms->Div[mapr]);

          for(cc[mapc]=0;cc[mapc]<ms->FDim[mapc];cc[mapc]++) {
            v[mapc]=(cc[mapc]+ms->Min[mapc])/((float)ms->Div[mapc]);

            if(normalize) 
              dens = (*f-mean)/stdev;
            else 
              dens = *f;
            F3(ms->Field->data,cc[0],cc[1],cc[2]) = dens;
            if(maxd<*f) maxd = dens;
            if(mind>*f) mind = dens;
            f++;
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            for(e=0;e<3;e++) 
              F4(ms->Field->points,cc[0],cc[1],cc[2],e) = vr[e];
          }
        }
      }
  }
  if(ok) {
    d = 0;
    for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
      {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            copy3f(vr,ms->Corner[d]);
            d++;
          }
        }
      }
  }

  PRINTFB(FB_Details,FB_ObjectMap) 
    " ObjectMapCCP4: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
    ENDFB;
  
  if(ok) {
    v[2]=(ms->Min[2])/((float)ms->Div[2]);
    v[1]=(ms->Min[1])/((float)ms->Div[1]);
    v[0]=(ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
    
    v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
    v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
    v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
  }
#ifdef _UNDEFINED
  printf("%d %d %d %d %d %d %d %d %d\n",
         ms->Div[0],
         ms->Min[0],
         ms->Max[0],
         ms->Div[1],
         ms->Min[1],
         ms->Max[1],
         ms->Div[2],
         ms->Min[2],
         ms->Max[2]);
  printf("Okay? %d\n",ok);
  fflush(stdout);
#endif
  if(!ok) {
    ErrMessage("ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
  }
  return(ok);
}
/*========================================================================*/
static int ObjectMapPHIStrToMap(ObjectMap *I,char *PHIStr,int bytes,int state) {
  
  char *p;
  float dens,dens_rev;
  int a,b,c,d,e;
  float v[3],maxd,mind;
  int ok = true;
  int little_endian = 1;
  /* PHI named from their docs */
  int map_endian = 0;
  int map_dim;
  int map_bytes;

  ObjectMapState *ms;

  char cc[MAXLINELEN];
  char *rev;

  little_endian = *((char*)&little_endian);


  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);

  maxd = FLT_MIN;
  mind = FLT_MAX;
  p=PHIStr;

  if(*p)  /* use FORMATTED IO record to determine map endiness */
    map_endian = 1;
  else
    map_endian = 0;

  p+=4;

  ParseNCopy(cc,p,20);
  PRINTFB(FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB;
  p+=20;
  p+=4;

  p+=4;
  ParseNCopy(cc,p,10);
  PRINTFB(FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB;
  p+=10;
  ParseNCopy(cc,p,60);
  PRINTFB(FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB;
  p+=60;
  p+=4;

  rev = (char*)&dens_rev;

  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];
  }

  map_bytes = *((int*)rev);

  map_dim = (int)(pow((map_bytes/4.0),1/3.0)+0.5);

  if((4*map_dim*map_dim*map_dim)!=map_bytes) /* consistency check */
    map_dim = 65;

  PRINTFB(FB_Details,FB_ObjectMap) 
      " PHIStrToMap: Map Size %d x %d x %d\n",map_dim,map_dim,map_dim
      ENDFB;
  p+=4;

  ms->FDim[0] = map_dim;
  ms->FDim[1] = map_dim;
  ms->FDim[2] = map_dim;
  ms->FDim[3] = 3;

  ms->Div[0] = (map_dim-1)/2;
  ms->Div[1] = (map_dim-1)/2;
  ms->Div[2] = (map_dim-1)/2;
  ms->Min[0] = -ms->Div[0];
  ms->Min[1] = -ms->Div[1];
  ms->Min[2] = -ms->Div[2];
  ms->Max[0] = ms->Div[0];
  ms->Max[1] = ms->Div[1];
  ms->Max[2] = ms->Div[2];

  ms->Field=IsosurfFieldAlloc(ms->FDim);
  ms->MapSource = cMapSourcePHI;
  ms->Field->save_points=false;


  for(c=0;c<ms->FDim[2];c++) { /* z y x ordering into c b a  so that x = a, etc. */
    for(b=0;b<ms->FDim[1];b++) {
      for(a=0;a<ms->FDim[0];a++) {

        if(little_endian!=map_endian) {
          rev[0]=p[3];
          rev[1]=p[2];
          rev[2]=p[1];
          rev[3]=p[0];
        } else {
          rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
          rev[1]=p[1];
          rev[2]=p[2];
          rev[3]=p[3];
        }
        dens = *((float*)rev);
        F3(ms->Field->data,a,b,c) = dens;
        if(maxd<dens) maxd = dens;
        if(mind>dens) mind = dens;
        p+=4;
      }
    }
  }
  p+=4;


  p+=4;
  ParseNCopy(cc,p,16);
  PRINTFB(FB_ObjectMap,FB_Details)
    " PHIStrToMap: %s\n",cc
    ENDFB;
  p+=16;
  p+=4;

  ms->Grid = Alloc(float,3);
  p+=4;
  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];
  }
  ms->Grid[0] = 1.0F/(*((float*)rev));
  ms->Grid[1] = ms->Grid[0];
  ms->Grid[2] = ms->Grid[0];
  p+=4;

  ms->Origin = Alloc(float,3);
  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];;
  }
  ms->Origin[0] = *((float*)rev);
  p+=4;

  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];;
  }
  ms->Origin[1] = *((float*)rev);
  p+=4;
  if(little_endian!=map_endian) {
    rev[0]=p[3];
    rev[1]=p[2];
    rev[2]=p[1];
    rev[3]=p[0];
  } else {
    rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
    rev[1]=p[1];
    rev[2]=p[2];
    rev[3]=p[3];;
  }
  ms->Origin[2] = *((float*)rev);
  p+=4;

  p+=4;

  if(ok) {
    for(e=0;e<3;e++) {
      ms->ExtentMin[e] = ms->Origin[e]+ms->Grid[e]*ms->Min[e];
      ms->ExtentMax[e] = ms->Origin[e]+ms->Grid[e]*ms->Max[e];
    }
  }

  for(c=0;c<ms->FDim[2];c++) {
    v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
    for(b=0;b<ms->FDim[1];b++) {
      v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
      for(a=0;a<ms->FDim[0];a++) {
        v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
        for(e=0;e<3;e++) {
          F4(ms->Field->points,a,b,c,e) = v[e];
        }
      }
    }
  }

  d=0;
  for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
    v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);

    for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
      v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);

      for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
        v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
        copy3f(v,ms->Corner[d]);
        d++;
      }
    }
  }

  /* interpolation test code 
  { 
    float test[3];
    float result;
    float cmp;

    for(c=0;c<ms->FDim[0]-1;c++) {
      for(b=0;b<ms->FDim[1]-1;b++) {
        for(a=0;a<ms->FDim[2]-1;a++) {
          for(e=0;e<3;e++) {
            test[e] = (F4(ms->Field->points,a,b,c,e)+
                       F4(ms->Field->points,a,b,c,e))/2.0;
          }
          ObjectMapStateInterpolate(ms,test,&result,1);
          cmp = (F3(ms->Field->data,a,b,c)+
                 F3(ms->Field->data,a,b,c))/2.0;
          if(fabs(cmp-result)>0.001) {
            printf("%d %d %d\n",a,b,c);
            printf("%8.5f %8.5f\n",
                   cmp,
                   result);
          }
        }
      }
    }
  }
  */

  if(!ok) {
    ErrMessage("ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);
  }
  return(ok);
}
/*========================================================================*/
int ObjectMapXPLORStrToMap(ObjectMap *I,char *XPLORStr,int state) {
  
  char *p;
  int a,b,c,d,e;
  float v[3],vr[3],dens,maxd,mind;
  char cc[MAXLINELEN];
  int n;
  int ok = true;
  ObjectMapState *ms;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }

  ms=&I->State[state];
  ObjectMapStateInit(ms);

  maxd = FLT_MIN;
  mind = FLT_MAX;
  p=XPLORStr;

  while(*p) {
    p = ParseNCopy(cc,p,8);
    if(!*cc) 
      p = ParseNextLine(p);
    else if(sscanf(cc,"%i",&n)==1) {
      p=ParseWordCopy(cc,p,MAXLINELEN);
      if(strstr(cc,"!NTITLE")||(!*cc)) {
        p=ParseNextLine(p);
        while(n--) {
          p=ParseNextLine(p);          
        } 
      } else if(strstr(cc,"REMARKS")) {
        p=ParseNextLine(p);          
      } else {
        break;
      }
    }
  }
  if(*p) { /* n contains first dimension */
    ms->Div[0]=n;
    if(sscanf(cc,"%i",&ms->Min[0])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Max[0])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Div[1])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Min[1])!=1) ok=false;    
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Max[1])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Div[2])!=1) ok=false;
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Min[2])!=1) ok=false;    
    p = ParseNCopy(cc,p,8); if(sscanf(cc,"%i",&ms->Max[2])!=1) ok=false;
    p=ParseNextLine(p);
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Dim[0])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Dim[1])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Dim[2])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Angle[0])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Angle[1])!=1) ok=false;
    p = ParseNCopy(cc,p,12); if(sscanf(cc,"%f",&ms->Crystal->Angle[2])!=1) ok=false;
    p=ParseNextLine(p);
    p = ParseNCopy(cc,p,3);
    if(strcmp(cc,"ZYX")) ok=false;
    p=ParseNextLine(p);
    
  } else {
    ok=false;
  }
  if(ok) {
    
    ms->FDim[0]=ms->Max[0]-ms->Min[0]+1;
    ms->FDim[1]=ms->Max[1]-ms->Min[1]+1;
    ms->FDim[2]=ms->Max[2]-ms->Min[2]+1;
    ms->FDim[3]=3;
    if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
      ok=false;
    else {
      CrystalUpdate(ms->Crystal);
      ms->Field=IsosurfFieldAlloc(ms->FDim);
      ms->MapSource = cMapSourceXPLOR;
      ms->Field->save_points=false;
      for(c=0;c<ms->FDim[2];c++)
        {
          v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
          p=ParseNextLine(p);
          for(b=0;b<ms->FDim[1];b++) {
            v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
            for(a=0;a<ms->FDim[0];a++) {
              v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
              p=ParseNCopy(cc,p,12);
              if(!cc[0]) {
                p=ParseNextLine(p);
                p=ParseNCopy(cc,p,12);                
              }
              if(sscanf(cc,"%f",&dens)!=1) {
                ok=false;
              } else {
                F3(ms->Field->data,a,b,c) = dens;
                if(maxd<dens) maxd = dens;
                if(mind>dens) mind = dens;
              }
              transform33f3f(ms->Crystal->FracToReal,v,vr);
              for(e=0;e<3;e++) {
                F4(ms->Field->points,a,b,c,e) = vr[e];
              }
            }
          }
          p=ParseNextLine(p);
        }
      if(ok) {
        d = 0;
        for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
          {
            v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
            for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
              v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
              for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                transform33f3f(ms->Crystal->FracToReal,v,vr);
                copy3f(vr,ms->Corner[d]);
                d++;
              }
            }
          }
      }
    }
  }

  if(ok) {
    v[2]=(ms->Min[2])/((float)ms->Div[2]);
    v[1]=(ms->Min[1])/((float)ms->Div[1]);
    v[0]=(ms->Min[0])/((float)ms->Div[0]);

    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
    
    v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
    v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
    v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);

    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);

  }
#ifdef _UNDEFINED
    printf("%d %d %d %d %d %d %d %d %d\n",
           ms->Div[0],
           ms->Min[0],
           ms->Max[0],
           ms->Div[1],
           ms->Min[1],
           ms->Max[1],
           ms->Div[2],
           ms->Min[2],
           ms->Max[2]);
    printf("Okay? %d\n",ok);
    fflush(stdout);
#endif
  if(!ok) {
    ErrMessage("ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
  }
    
  return(ok);
}
/*========================================================================*/
static int ObjectMapFLDStrToMap(ObjectMap *I,char *PHIStr,int bytes,int state) 
{
  char *p;
  float dens,dens_rev;
  int a,b,c,d,e;
  float v[3],maxd,mind;
  int ok = true;
  int little_endian = 1;
  /* PHI named from their docs */
  int map_endian = 1;
  char cc[MAXLINELEN];
  char *rev;
  int ndim=0;
  int veclen=0;
  int nspace=0;
  ObjectMapState *ms;
    int got_ndim=false;
    int got_dim1=false;
    int got_dim2=false;
    int got_dim3=false;
    int got_data=false;
    int got_veclen=false;
    int got_min_ext=false;
    int got_max_ext=false;
    int got_field=false;
    int got_nspace=false;


  little_endian = *((char*)&little_endian);
  map_endian = little_endian;

  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);

  maxd = FLT_MIN;
  mind = FLT_MAX;

  p = PHIStr;

  while((*p)&&(!(
                 got_ndim&&
                 got_dim1&&
                 got_dim2&&
                 got_dim3&&
                 got_veclen&&
                 got_min_ext&&
                 got_max_ext&&
                 got_field&&
                 got_nspace&&
                 got_data))) {

    if(!got_ndim) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"ndim")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ndim)==1)
          got_ndim=true;
      }
    }
      
    if(!got_dim1) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"dim1")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[0])==1)
          got_dim1=true;
      }
    }
      
    if(!got_dim2) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"dim2")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[1])==1)
          got_dim2=true;
      }
    }
      
    if(!got_dim3) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"dim3")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[2])==1)
          got_dim3=true;
      }
    }
       
    if(!got_nspace) {
      ParseWordCopy(cc,p,6);
      if(!strcmp(cc,"nspace")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&nspace)==1)
          got_nspace=true;
      }
    }

    if(!got_veclen) {
      ParseWordCopy(cc,p,6);
      if(!strcmp(cc,"veclen")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&veclen)==1)
          got_veclen=true;
      }
    }

    if(!got_data) {
      ParseWordCopy(cc,p,4);
      if(!strcmp(cc,"data")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(!strcmp(cc,"float"))
          got_data=true;
      }
    }

    if(!got_min_ext) {
      ParseWordCopy(cc,p,7);
      if(!strcmp(cc,"min_ext")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%f",&ms->ExtentMin[0])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%f",&ms->ExtentMin[1])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->ExtentMin[2])==1) {
              got_min_ext=true;
            }
          }
        }
      }
    }

    if(!got_max_ext) {
      ParseWordCopy(cc,p,7);
      if(!strcmp(cc,"max_ext")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%f",&ms->ExtentMax[0])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%f",&ms->ExtentMax[1])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->ExtentMax[2])==1) {
              got_max_ext=true;
            }
          }
        }
      }
    }

    if(!got_field) {
      ParseWordCopy(cc,p,5);
      if(!strcmp(cc,"field")) {
        p = ParseSkipEquals(p);
        p = ParseWordCopy(cc,p,50);
        if(!strcmp(cc,"uniform"))
          got_field=true;
      }
    }
    p=ParseNextLine(p);
  }
    
  if(got_ndim&&
     got_dim1&&
     got_dim2&&
     got_dim3&&
     got_veclen&&
     got_min_ext&&
     got_max_ext&&
     got_field&&
     got_nspace&&
     got_data) {

    int pass = 0;

    ms->Origin = Alloc(float,3);
    ms->Range = Alloc(float,3);
    ms->Grid = Alloc(float,3);

    copy3f(ms->ExtentMin,ms->Origin);
    subtract3f(ms->ExtentMax,ms->ExtentMin,ms->Range);
    ms->FDim[3] = 3;

    PRINTFB(FB_Details,FB_ObjectMap) 
      " FLDStrToMap: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB;
      
    for(a=0;a<3;a++) {
      ms->Min[a] = 0;
      ms->Max[a] = ms->FDim[a]-1;
      ms->Div[a] = ms->FDim[a];
        
      if(ms->FDim[a]) 
        ms->Grid[a]=ms->Range[a]/(ms->Max[a]);
      else
        ms->Grid[a]=0.0F;
    }

    ms->Field=IsosurfFieldAlloc(ms->FDim);
    ms->MapSource = cMapSourceFLD;
    ms->Field->save_points=false;

    
    while(1) {
      maxd = FLT_MIN;
      mind = FLT_MAX;
      p=PHIStr;
      
      while(*p) { /* ^L^L sentinel */
        if((p[0]==12)&&(p[1]==12)) {
          p+=2;
          break;
        }
        p++;
      }
      
      rev = (char*)&dens_rev;
      for(c=0;c<ms->FDim[2];c++) { /* z y x ordering into c b a  so that x = a, etc. */
        for(b=0;b<ms->FDim[1];b++) {
          for(a=0;a<ms->FDim[0];a++) {
            
            if(little_endian!=map_endian) {
              rev[0]=p[3];
              rev[1]=p[2];
              rev[2]=p[1];
              rev[3]=p[0];
            } else {
              rev[0]=p[0]; /* gotta go char by char because of memory alignment issues ... */
              rev[1]=p[1];
              rev[2]=p[2];
              rev[3]=p[3];
            }
            dens = *((float*)rev);
            F3(ms->Field->data,a,b,c) = dens;
            if(maxd<dens) maxd = dens;
            if(mind>dens) mind = dens;
            p+=4;
          }
        }
      }

      /* There's no way to determine the original handedness of input
          field files.  So instead, we simplymake an educated guess about
       whether we're byte-swapped based on the range of the density
       values obtained. */
      
      if(((maxd/FLT_MAX)>0.1F)&&((mind/(-FLT_MAX))>0.1F)) {
        if(pass==0) {
          map_endian = (!map_endian); /* okay, try again swapped */
        } else if(pass==1) {
          /* didn't help, so resort to original order */
          map_endian = (!map_endian); 
        } else {
          break;
        }
      } else {
        break;
      }
      pass++;
    }
    
    for(c=0;c<ms->FDim[2];c++) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      for(b=0;b<ms->FDim[1];b++) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        for(a=0;a<ms->FDim[0];a++) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          for(e=0;e<3;e++) {
            F4(ms->Field->points,a,b,c,e) = v[e];
          }
        }
      }
    }
    
    d=0;
    for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      if(!c) 
        v[2]=ms->ExtentMin[2];
      else 
        v[2]=ms->ExtentMax[2];
      
      for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);

      if(!b) 
        v[1]=ms->ExtentMin[1];
      else 
        v[1]=ms->ExtentMax[1];
        
        for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);

      if(!a) 
        v[0]=ms->ExtentMin[0];
      else 
        v[0]=ms->ExtentMax[0];

          copy3f(v,ms->Corner[d]);
          d++;
        }
      }
    }
    
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);
      
  } else {
    PRINTFB(FB_Errors,FB_ObjectMap) 
      " Error: unable to read FLD file.\n"
      ENDFB;
            /*  printf(" got_ndim %d\n",got_ndim);
    printf(" got_dim1 %d\n",got_dim1);
    printf(" got_dim2 %d\n",got_dim2);
    printf(" got_dim3 %d\n",got_dim3);
    printf(" got_veclen %d\n",got_veclen);
    printf(" got_min_ext %d\n",got_min_ext);
    printf(" got_max_ext %d\n",got_max_ext);
    printf(" got_field %d\n",got_field);
    printf(" got_nspace %d\n",got_nspace);
    printf(" got_data %d\n",got_data);
            */
  }
    

  return(ok);

}
/*========================================================================*/
static int ObjectMapBRIXStrToMap(ObjectMap *I,char *BRIXStr,int bytes,int state) 
{
  char *p,*pp;
  float dens;
  int a,b,c,d,e;
  float v[3],vr[3],maxd,mind;
  int ok = true;
  char cc[MAXLINELEN];
  ObjectMapState *ms;
  int got_origin=false;
  int got_extent=false;
  int got_grid=false;
  int got_cell=false;
  int got_prod=false;
  int got_plus=false;
  int got_sigma=false;
  int got_smiley=false;
  float prod=1.0F;
  float plus=0.0F;
  float sigma=0.0F;
  float calc_sigma=0.0F;
  float calc_mean=0.0F;
  int normalize;
  int swap_bytes;

  normalize=(int)SettingGet(cSetting_normalize_o_maps);
  swap_bytes=(int)SettingGet(cSetting_swap_dsn6_bytes);
  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);

  maxd = FLT_MIN;
  mind = FLT_MAX;
  
  p = BRIXStr;

  ParseNCopy(cc,p,3);
  got_smiley = (strcmp(cc,":-)")==0);

  if(got_smiley) {

    /* ASCII BRIX format header */

    while((*p)&&(!(
                   got_origin&&
                   got_extent&&
                   got_grid&&
                   got_cell&&
                   got_prod&&
                   got_plus&&
                   got_sigma))) {
      
      if(!got_origin) {
        pp=ParseWordCopy(cc,p,6);
        if(WordMatch("origin",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%d",&ms->Min[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Min[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Min[2])==1) {
                got_origin=true;
              }
            }
          }
        }
      }
      
      if(!got_extent) {
        pp=ParseWordCopy(cc,p,6);
        if(WordMatch("extent",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%d",&ms->Max[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Max[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Max[2])==1) {
                got_extent=true;
                ms->Max[0]+=ms->Min[0]-1;
                ms->Max[1]+=ms->Min[1]-1;
                ms->Max[2]+=ms->Min[2]-1;
              }
            }
          }
        }
      }
      
      
      if(!got_grid) {
        pp=ParseWordCopy(cc,p,4);
        if(WordMatch("grid",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%d",&ms->Div[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Div[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Div[2])==1) {
                got_grid=true;
              }
            }
          }
        }
      }
      
      if(!got_cell) {
        pp = ParseWordCopy(cc,p,4);
        if(WordMatch("cell",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          
          if(sscanf(cc,"%f",&ms->Crystal->Dim[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->Crystal->Dim[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%f",&ms->Crystal->Dim[2])==1) {
                p = ParseWordCopy(cc,p,50);
                if(sscanf(cc,"%f",&ms->Crystal->Angle[0])==1) {
                  p = ParseWordCopy(cc,p,50);
                  if(sscanf(cc,"%f",&ms->Crystal->Angle[1])==1) {
                    p = ParseWordCopy(cc,p,50);
                    if(sscanf(cc,"%f",&ms->Crystal->Angle[2])==1) {
                      got_cell=true;
                    }
                  }
                }
              }
            }
          }
        }
      }
      
      if(!got_plus) {
        pp=ParseWordCopy(cc,p,4);
        if(WordMatch("plus",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%f",&plus)==1) {
            got_plus=true;
          }
        }
      }
      
      if(!got_prod) {
        pp=ParseWordCopy(cc,p,4);
        if(WordMatch("prod",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%f",&prod)==1) {
            got_prod=true;
          }
        }
      }
      
      if(!got_sigma) {
        pp = ParseWordCopy(cc,p,5);
        if(WordMatch("sigma",cc,true)<0) {
          p = ParseWordCopy(cc,pp,50);
          if(sscanf(cc,"%f",&sigma)==1) {
            got_sigma=true;
          }
        }
      }
      
      p++;
    }
  } else {
    /* Binary DSN6 format */

    /* first, determine whether or not we need to byte-swap the header... */

    int passed_endian_check = false;
    short int *shint_ptr;
    float scale_factor;
    char tmp_char;
    int a;
    int start_swap_at = 0;
    int stop_swap_at = bytes;

    shint_ptr = (short int*)(p+18*2);
    if(*shint_ptr==100) {
      passed_endian_check = true;
      start_swap_at = 512; /* don't byte-swap header */
    }
    
    if(!swap_bytes) {
      stop_swap_at = 512;
    }
    /* for DSN6, always swap the bricks */

    p = BRIXStr + start_swap_at;
    for( a=start_swap_at; a<stop_swap_at; a+=2) {
      tmp_char = *p;
      *p = *(p+1);
      *(p+1) = tmp_char;
      p+=2;
    }

    p = BRIXStr;

    if(*shint_ptr==100) {
      passed_endian_check = true;
    }
    
    if(!passed_endian_check) {
      PRINTFB(FB_Errors,FB_ObjectMap) 
        " Error: This looks like a DSN6 map file, but I can't match endianness.\n"
        ENDFB;
    } else {
      shint_ptr = (short int*)p;
      
      ms->Min[0] = *(shint_ptr++);
      ms->Min[1] = *(shint_ptr++);
      ms->Min[2] = *(shint_ptr++);
      got_origin = true;
      
      ms->Max[0] = *(shint_ptr++);
      ms->Max[1] = *(shint_ptr++);
      ms->Max[2] = *(shint_ptr++);
      got_extent=true;
      ms->Max[0]+=ms->Min[0]-1;
      ms->Max[1]+=ms->Min[1]-1;
      ms->Max[2]+=ms->Min[2]-1;
      
      ms->Div[0] = *(shint_ptr++);
      ms->Div[1] = *(shint_ptr++);
      ms->Div[2] = *(shint_ptr++);
      got_grid = true;
      
      ms->Crystal->Dim[0] = (float)(*(shint_ptr++));
      ms->Crystal->Dim[1] = (float)(*(shint_ptr++));
      ms->Crystal->Dim[2] = (float)(*(shint_ptr++));
      ms->Crystal->Angle[0] = (float)(*(shint_ptr++));
      ms->Crystal->Angle[1] = (float)(*(shint_ptr++));
      ms->Crystal->Angle[2] = (float)(*(shint_ptr++));
      got_cell = true;
      
      prod = (float)(*(shint_ptr++)) / 100.0F;
      got_prod = true;
      
      plus = (float)(*(shint_ptr++));
      got_plus = true;
      
      scale_factor = (float)(*(shint_ptr++));
      
      for(a=0;a<3;a++) {
        ms->Crystal->Dim[a]/=scale_factor;
        ms->Crystal->Angle[a]/= scale_factor;
      }
    }
  }
  
  if(got_origin&&
     got_extent&&
     got_grid&&
     got_cell&&
     got_plus&&
     got_prod) {

    
    PRINTFB(FB_Blather,FB_ObjectMap) 
      " BRIXStrToMap: Prod = %8.3f, Plus = %8.3f\n",prod,plus
      ENDFB;

    ms->FDim[0]=ms->Max[0]-ms->Min[0]+1;
    ms->FDim[1]=ms->Max[1]-ms->Min[1]+1;
    ms->FDim[2]=ms->Max[2]-ms->Min[2]+1;
    ms->FDim[3]=3;
    if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
      ok=false;
    else {
      CrystalUpdate(ms->Crystal);
      ms->Field=IsosurfFieldAlloc(ms->FDim);
      ms->MapSource = cMapSourceBRIX;
      ms->Field->save_points=false;
      
      {
        int block_size = 8;
        int xblocks = ((ms->FDim[0]-1)/block_size)+1;
        int yblocks = ((ms->FDim[1]-1)/block_size)+1;
        int zblocks = ((ms->FDim[2]-1)/block_size)+1;
        int cc,bb,aa,ia,ib,ic,xa,xb,xc;
        double sum = 0.0;
        double sumsq = 0.0;
        int n_pts = 0;

        p=BRIXStr + 512;
        for(cc=0;cc<zblocks;cc++) {
          for(bb=0;bb<yblocks;bb++) {
            for(aa=0;aa<xblocks;aa++) {
              
              for(c=0;c<block_size;c++) {
                xc = c + cc*block_size;
                ic = xc + ms->Min[2];
                v[2]=ic/((float)ms->Div[2]);
                
                for(b=0;b<block_size;b++) {
                  xb = b + bb*block_size;
                  ib = xb + ms->Min[1];
                  v[1]=ib/((float)ms->Div[1]);
                  
                  for(a=0;a<block_size;a++) {
                    xa = a + aa*block_size;
                    ia = xa + ms->Min[0];
                    v[0]=ia/((float)ms->Div[0]);
                    
                    dens = (((float)(*((unsigned char*)(p++)))) - plus) / prod;
                    if((ia<=ms->Max[0])&&(ib<=ms->Max[1])&&(ic<=ms->Max[2])) {
                      F3(ms->Field->data,xa,xb,xc) = dens;
                      sumsq+=dens*dens;
                      sum+=dens;
                      n_pts++;
                      if(maxd<dens) maxd = dens;
                      if(mind>dens) mind = dens;
                      transform33f3f(ms->Crystal->FracToReal,v,vr);
                      for(e=0;e<3;e++) {
                        F4(ms->Field->points,xa,xb,xc,e) = vr[e];
                      }
                      
                    }
                  }
                }
              }
            }
          }
        }

        if(n_pts>1) {
          calc_mean = (float)(sum/n_pts);
          calc_sigma = (float)sqrt1d((sumsq - (sum*sum/n_pts))/(n_pts-1));
        }
        
      }
      
      if(ok) {
        d = 0;
        for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
          {
            v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
            for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
              v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
              for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                transform33f3f(ms->Crystal->FracToReal,v,vr);
                copy3f(vr,ms->Corner[d]);
                d++;
              }
            }
          }
      }

      if(ok&&normalize) {

        if(calc_sigma>R_SMALL8) {
          for(c=0;c<ms->FDim[2];c++) {
            for(b=0;b<ms->FDim[1];b++) {
              for(a=0;a<ms->FDim[0];a++) {
                F3(ms->Field->data,a,b,c) =  (F3(ms->Field->data,a,b,c)-calc_mean)/calc_sigma;
              }
            }
          }
        }
      }

      v[2]=(ms->Min[2])/((float)ms->Div[2]);
      v[1]=(ms->Min[1])/((float)ms->Div[1]);
      v[0]=(ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
      
      v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
      v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
      v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
      
      PRINTFB(FB_Details,FB_ObjectMap) 
        " BRIXStrToMap: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
        ENDFB;

      if(got_sigma) {
        PRINTFB(FB_Details,FB_ObjectMap) 
          " BRIXStrToMap: Reported Sigma = %8.3f\n",sigma
          ENDFB;
      }

      PRINTFB(FB_Details,FB_ObjectMap)
        " BRIXStrToMap: Range = %5.6f to %5.6f\n",mind,maxd
        ENDFB;

      PRINTFB(FB_Details,FB_ObjectMap) 
        " BRIXStrToMap: Calculated Mean = %8.3f, Sigma = %8.3f\n",calc_mean,calc_sigma
        ENDFB;

      if(normalize) {
        PRINTFB(FB_Details,FB_ObjectMap) 
          " BRIXStrToMap: Normalizing...\n"
          ENDFB;
      }

      ms->Active=true;
      ObjectMapUpdateExtents(I);

    }
  } else {
    PRINTFB(FB_Errors,FB_ObjectMap) 
      " Error: unable to read BRIX/DSN6 file.\n"
      ENDFB;
  }
  
  
  return(ok);
  
}
/*========================================================================*/
static int ObjectMapGRDStrToMap(ObjectMap *I,char *GRDStr,int bytes,int state) 
{
  char *p;
  float dens;
  int a,b,c,d,e;
  float v[3],vr[3],maxd,mind;
  int ok = true;
  char cc[MAXLINELEN];
  ObjectMapState *ms;
  int got_cell=false;
  int got_brick=false;
  int fast_axis=1;
  int got_ranges=false;
  int normalize = false;
  float mean,stdev;
  double sum = 0.0;
  double sumsq = 0.0;
  int n_pts = 0;
  
  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);
  normalize=(int)SettingGet(cSetting_normalize_grd_maps);
  maxd = FLT_MIN;
  mind = FLT_MAX;
  
  p = GRDStr;

  /* print map title */

  p = ParseNCopy(cc,p,100);

  PRINTFB(FB_ObjectMap,FB_Details)
    " ObjectMap: %s\n",cc
    ENDFD;

  p = ParseNextLine(p);

  /* skip format -- we're reading float regardless... */

  p = ParseNextLine(p);

  /* read unit cell */

  if(ok) {
    p = ParseWordCopy(cc,p,50);
    if(sscanf(cc,"%f",&ms->Crystal->Dim[0])==1) {
      p = ParseWordCopy(cc,p,50);
      if(sscanf(cc,"%f",&ms->Crystal->Dim[1])==1) {
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%f",&ms->Crystal->Dim[2])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%f",&ms->Crystal->Angle[0])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%f",&ms->Crystal->Angle[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%f",&ms->Crystal->Angle[2])==1) {
                got_cell=true;
              }
            }
          }
        }
      }
      ok = got_cell;
    }
  }

  p=ParseNextLine(p);

  /* read brick dimensions */

  if(ok) {
    p = ParseWordCopy(cc,p,50);
    if(sscanf(cc,"%d",&ms->FDim[0])==1) {
      p = ParseWordCopy(cc,p,50);
      if(sscanf(cc,"%d",&ms->FDim[1])==1) {
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->FDim[2])==1) {
          got_brick=true;
          ms->FDim[0]++;
          ms->FDim[1]++;
          ms->FDim[2]++; /* stupid 0-based counters */
        }
      }
    }
    ok = got_brick;
  }

  p=ParseNextLine(p);

  /* read ranges */

  if(ok) {
    p = ParseWordCopy(cc,p,50);
    if(sscanf(cc,"%d",&fast_axis)==1) {
      p = ParseWordCopy(cc,p,50);      
      if(sscanf(cc,"%d",&ms->Min[0])==1) {
        p = ParseWordCopy(cc,p,50);
        if(sscanf(cc,"%d",&ms->Div[0])==1) {
          p = ParseWordCopy(cc,p,50);
          if(sscanf(cc,"%d",&ms->Min[1])==1) {
            p = ParseWordCopy(cc,p,50);
            if(sscanf(cc,"%d",&ms->Div[1])==1) {
              p = ParseWordCopy(cc,p,50);
              if(sscanf(cc,"%d",&ms->Min[2])==1) {
                p = ParseWordCopy(cc,p,50);
                if(sscanf(cc,"%d",&ms->Div[2])==1) {
                  got_ranges=true;
                }
              }
            }
          }
        }
      }
    }
    ok = got_ranges;
  }
  
  if(ok) {

    ms->Div[0]= ms->Div[0] - ms->Min[0];
    ms->Div[1]= ms->Div[1] - ms->Min[1];
    ms->Div[2]= ms->Div[2] - ms->Min[2];

    ms->Max[0]=ms->Min[0]+ms->FDim[0]-1;
    ms->Max[1]=ms->Min[1]+ms->FDim[1]-1;
    ms->Max[2]=ms->Min[2]+ms->FDim[2]-1;

    ms->FDim[3]=3;

    if(Feedback(FB_ObjectMap,FB_Blather)) {
      dump3i(ms->Div,"ms->Div");
      dump3i(ms->Min,"ms->Min");
      dump3i(ms->Max,"ms->Max");
      dump3i(ms->FDim,"ms->FDim");
    }

    CrystalUpdate(ms->Crystal);
    ms->Field=IsosurfFieldAlloc(ms->FDim);
    ms->MapSource = cMapSourceGRD;
    ms->Field->save_points=false;

    switch(fast_axis) {
    case 3: /* Fast Y - UNTESTED! */
    default:
    case 1: /* Fast X */
      for(c=0;c<ms->FDim[2];c++)
        {
          v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
          for(b=0;b<ms->FDim[1];b++) {
            v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
            for(a=0;a<ms->FDim[0];a++) {
              v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
              p=ParseNextLine(p);
              p=ParseNCopy(cc,p,24);                
              if(sscanf(cc,"%f",&dens)!=1) {
                ok=false;
              } else {
                sumsq+=dens*dens;
                sum+=dens;
                n_pts++;
                F3(ms->Field->data,a,b,c) = dens;
                if(maxd<dens) maxd = dens;
                if(mind>dens) mind = dens;
              }
              transform33f3f(ms->Crystal->FracToReal,v,vr);
              for(e=0;e<3;e++) {
                F4(ms->Field->points,a,b,c,e) = vr[e];
              }
            }
          }
        }
      break;
    }
  }   

  if(n_pts>1) {
    mean = (float)(sum/n_pts);
    stdev = (float)sqrt1d((sumsq - (sum*sum/n_pts))/(n_pts-1));

    if(normalize) {
      PRINTFB(FB_ObjectMap,FB_Details)
        " ObjectMapGRDStrToMap: Normalizing: mean = %8.6f and stdev = %8.6f.\n",
        mean,stdev
        ENDFB;
      if(stdev<R_SMALL8)
        stdev = 1.0F;
      for(c=0;c<ms->FDim[2];c++)
        for(b=0;b<ms->FDim[1];b++) {
          for(a=0;a<ms->FDim[0];a++) {
            dens = F3(ms->Field->data,a,b,c);
            F3(ms->Field->data,a,b,c) = (dens-mean)/stdev;
          }
        }
    } else {
      PRINTFB(FB_ObjectMap,FB_Details)
        " ObjectMapGRDStrToMap: Mean = %8.6f and stdev = %8.6f.\n",
        mean,stdev
        ENDFB;
    }
  }

  if(ok) {
    d = 0;
    for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
      {
        v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
        for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
          v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
          for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
            v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
            transform33f3f(ms->Crystal->FracToReal,v,vr);
            copy3f(vr,ms->Corner[d]);
            d++;
          }
        }
      }
    
    v[2]=(ms->Min[2])/((float)ms->Div[2]);
    v[1]=(ms->Min[1])/((float)ms->Div[1]);
    v[0]=(ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
    
    v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
    v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
    v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
    
    transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);
    
    PRINTFB(FB_Details,FB_ObjectMap) 
      " GRDXStrToMap: Map Size %d x %d x %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB;
    
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);

  } else {
    PRINTFB(FB_Errors,FB_ObjectMap) 
      " Error: unable to read GRD file.\n"
      ENDFB;
  }
  
  return(ok);
  
}
/*========================================================================*/
ObjectMap *ObjectMapReadXPLORStr(ObjectMap *I,char *XPLORStr,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapXPLORStrToMap(I,XPLORStr,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}
/*========================================================================*/
ObjectMap *ObjectMapReadCCP4Str(ObjectMap *I,char *XPLORStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapCCP4StrToMap(I,XPLORStr,bytes,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}
/*========================================================================*/
ObjectMap *ObjectMapLoadCCP4File(ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage("ObjectMapLoadCCP4File","Unable to open file!");
  else
       {
            if(Feedback(FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadCCP4File: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadCCP4Str(obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
#ifdef _UNDEFINED
  {int a;
  for(a=0;a<9;a++)
    printf("%10.5f\n",mat[a]);
  }
#endif
  return(I);

}
/*========================================================================*/
static ObjectMap *ObjectMapReadFLDStr(ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapFLDStrToMap(I,MapStr,bytes,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}

/*========================================================================*/
static ObjectMap *ObjectMapReadBRIXStr(ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapBRIXStrToMap(I,MapStr,bytes,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}
/*========================================================================*/
static ObjectMap *ObjectMapReadGRDStr(ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapGRDStrToMap(I,MapStr,bytes,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}

/*========================================================================*/
static ObjectMap *ObjectMapReadPHIStr(ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapPHIStrToMap(I,MapStr,bytes,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}
/*========================================================================*/
ObjectMap *ObjectMapLoadPHIFile(ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage("ObjectMapLoadPHIFile","Unable to open file!");
  else
       {
            if(Feedback(FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadPHIFile: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadPHIStr(obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/
static int ObjectMapDXStrToMap(ObjectMap *I,char *DXStr,int bytes,int state) {

  int n_items = 0;

  char *p;
  float dens;
  int a,b,c,d,e;
  float v[3],maxd,mind;
  int ok = true;
  /* DX named from their docs */

  int stage = 0;

  ObjectMapState *ms;

  char cc[MAXLINELEN];


  if(state<0) state=I->NState;
  if(I->NState<=state) {
    VLACheck(I->State,ObjectMapState,state);
    I->NState=state+1;
  }
  ms=&I->State[state];
  ObjectMapStateInit(ms);

  ms->Origin=Alloc(float,3);
  ms->Grid = Alloc(float,3);

  maxd = FLT_MIN;
  mind = FLT_MAX;
  p=DXStr;

  /* get the dimensions */

  ms->FDim[3] = 3;

  while(ok&&(*p)&&(stage==0)) {
    p = ParseNCopy(cc,p,35);
    if(strcmp(cc,"object 1 class gridpositions counts")==0) {
      p = ParseWordCopy(cc,p,10);      
      if(sscanf(cc,"%d",&ms->FDim[0])==1) {
        p = ParseWordCopy(cc,p,10);      
        if(sscanf(cc,"%d",&ms->FDim[1])==1) {
          p = ParseWordCopy(cc,p,10);      
          if(sscanf(cc,"%d",&ms->FDim[2])==1) {
            stage = 1;
          }
        }
      }
    }
    p = ParseNextLine(p);
  }

  if(ok&&(stage==1)) {
    PRINTFB(FB_ObjectMap,FB_Details)
      " DXStrToMap: Dimensions: %d %d %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]
      ENDFB;
  }

  /* get the origin */

  while(ok&&(*p)&&(stage==1)) {  
    p = ParseNCopy(cc,p,6);
    if(strcmp(cc,"origin")==0) {
      p = ParseWordCopy(cc,p,20);      
      if(sscanf(cc,"%f",&ms->Origin[0])==1) {
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%f",&ms->Origin[1])==1) {
          p = ParseWordCopy(cc,p,20);      
          if(sscanf(cc,"%f",&ms->Origin[2])==1) {
            stage = 2;
          }
        }
      }
    }
    p = ParseNextLine(p);
    
  }

  if(ok&&(stage==2)) {
    PRINTFB(FB_ObjectMap,FB_Details)
      " DXStrToMap: Origin %8.3f %8.3f %8.3f\n",ms->Origin[0],ms->Origin[1],ms->Origin[2]
      ENDFB;
  }


  while(ok&&(*p)&&(stage==2)) {  
    p = ParseNCopy(cc,p,5);
    if(strcmp(cc,"delta")==0) {
      p = ParseWordCopy(cc,p,20);      
      if(sscanf(cc,"%f",&ms->Grid[0])==1) {
        p = ParseNextLine(p);
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%f",&ms->Grid[1])==1) {
        p = ParseNextLine(p);
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          p = ParseWordCopy(cc,p,20);      
          if(sscanf(cc,"%f",&ms->Grid[2])==1) {
            stage = 3;
          }
        }
      }
    }
    p = ParseNextLine(p);
    
  }

  if(ok&&(stage==3)) {
    PRINTFB(FB_ObjectMap,FB_Details)
      " DXStrToMap: Grid %8.3f %8.3f %8.3f\n",ms->Grid[0],ms->Grid[1],ms->Grid[2]
      ENDFB;
  }

  while(ok&&(*p)&&(stage==3)) {  
    p = ParseNCopy(cc,p,6);
    if(strcmp(cc,"object")==0) {
      p = ParseWordCopy(cc,p,20);      
      p = ParseNTrim(cc,p,29);      
      if(strcmp(cc,"class array type double rank")==0) {
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        p = ParseWordCopy(cc,p,20);      
        if(sscanf(cc,"%d",&n_items)==1) {
          if(n_items == ms->FDim[0]*ms->FDim[1]*ms->FDim[2])
            stage = 4;
        }
      }
    }
    p = ParseNextLine(p);
  }

  if(stage == 4) {

    if(ok&&(stage==4)) {
      PRINTFB(FB_ObjectMap,FB_Details)
        " DXStrToMap: %d data points.\n",n_items
        ENDFB;
    }
    
    ms->Field=IsosurfFieldAlloc(ms->FDim);
    ms->MapSource = cMapSourcePHI;
    ms->Field->save_points=false;
    
    for(a=0;a<3;a++) {
      ms->Div[a] = ms->FDim[a] - 1;
      ms->Min[a] = 0;
      ms->Max[a] = ms->FDim[a] - 1;
    }
    
    for(a=0;a<ms->FDim[0];a++) {
      for(b=0;b<ms->FDim[1];b++) {
        for(c=0;c<ms->FDim[2];c++) { 
          
          p = ParseWordCopy(cc,p,20);      
          if(!cc[0]) {
            p = ParseNextLine(p);           
            p = ParseWordCopy(cc,p,20);      
          }
          if(sscanf(cc,"%f",&dens)==1) {
            if(maxd<dens) maxd = dens;
            if(mind>dens) mind = dens;
            F3(ms->Field->data,a,b,c) = dens;
          } else {
            ok=false;
          }
        }
      }

    }
    
    for(e=0;e<3;e++) {
      ms->ExtentMin[e] = ms->Origin[e]+ms->Grid[e]*ms->Min[e];
      ms->ExtentMax[e] = ms->Origin[e]+ms->Grid[e]*ms->Max[e];
    }
    
    for(c=0;c<ms->FDim[2];c++) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      for(b=0;b<ms->FDim[1];b++) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        for(a=0;a<ms->FDim[0];a++) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          for(e=0;e<3;e++) {
            F4(ms->Field->points,a,b,c,e) = v[e];
          }
        }
      }
    }
    
    d=0;
    for(c=0;c<ms->FDim[2];c+=ms->FDim[2]-1) {
      v[2]=ms->Origin[2]+ms->Grid[2]*(c+ms->Min[2]);
      
      for(b=0;b<ms->FDim[1];b+=ms->FDim[1]-1) {
        v[1]=ms->Origin[1]+ms->Grid[1]*(b+ms->Min[1]);
        
        for(a=0;a<ms->FDim[0];a+=ms->FDim[0]-1) {
          v[0]=ms->Origin[0]+ms->Grid[0]*(a+ms->Min[0]);
          copy3f(v,ms->Corner[d]);
          d++;
        }
      }
    }
    if(ok)
      stage = 5;
  }

  if(stage!=5)
    ok=false;

  if(!ok) {
    ErrMessage("ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    ObjectMapUpdateExtents(I);
    printf(" ObjectMap: Map Read.  Range = %5.6f to %5.6f\n",mind,maxd);
  }
  return(ok);
}

/*========================================================================*/
static ObjectMap *ObjectMapReadDXStr(ObjectMap *I,char *MapStr,int bytes,int state)
{
  int ok=true;
  int isNew = true;

  if(!I) 
       isNew=true;
  else 
       isNew=false;
  if(ok) {
       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    ObjectMapDXStrToMap(I,MapStr,bytes,state);
    SceneChanged();
    SceneCountFrames();
  }
  return(I);
}

/*========================================================================*/
ObjectMap *ObjectMapLoadDXFile(ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage("ObjectMapLoadDXFile","Unable to open file!");
  else
       {
            if(Feedback(FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadDXFile: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadDXStr(obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}

/*========================================================================*/
ObjectMap *ObjectMapLoadFLDFile(ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
       ok=ErrMessage("ObjectMapLoadFLDFile","Unable to open file!");
  else
       {
            if(Feedback(FB_ObjectMap,FB_Actions))
              {
                  printf(" ObjectMapLoadFLDFile: Loading from '%s'.\n",fname);
              }
            
            fseek(f,0,SEEK_END);
      size=ftell(f);
            fseek(f,0,SEEK_SET);

            buffer=(char*)mmalloc(size);
            ErrChkPtr(buffer);
            p=buffer;
            fseek(f,0,SEEK_SET);
            fread(p,size,1,f);
            fclose(f);

            I=ObjectMapReadFLDStr(obj,buffer,size,state);

            mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}

/*========================================================================*/
ObjectMap *ObjectMapLoadBRIXFile(ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
    ok=ErrMessage("ObjectMapLoadBRIXFile","Unable to open file!");
  if(ok)
       {
      if(Feedback(FB_ObjectMap,FB_Actions))
        {
          printf(" ObjectMapLoadBRIXFile: Loading from '%s'.\n",fname);
        }
      
      fseek(f,0,SEEK_END);
      size=ftell(f);
      fseek(f,0,SEEK_SET);
      
      buffer=(char*)mmalloc(size+255);
      ErrChkPtr(buffer);
      p=buffer;
      fseek(f,0,SEEK_SET);
      fread(p,size,1,f);
      p[size]=0;
      fclose(f);
      
            I=ObjectMapReadBRIXStr(obj,buffer,size,state);

      mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/
ObjectMap *ObjectMapLoadGRDFile(ObjectMap *obj,char *fname,int state)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  long size;
  char *buffer,*p;
  float mat[9];

  f=fopen(fname,"rb");
  if(!f)
    ok=ErrMessage("ObjectMapLoadGRDFile","Unable to open file!");
  if(ok)
       {
      if(Feedback(FB_ObjectMap,FB_Actions))
        {
          printf(" ObjectMapLoadGRDFile: Loading from '%s'.\n",fname);
        }
      
      fseek(f,0,SEEK_END);
      size=ftell(f);
      fseek(f,0,SEEK_SET);
      
      buffer=(char*)mmalloc(size+255);
      ErrChkPtr(buffer);
      p=buffer;
      fseek(f,0,SEEK_SET);
      fread(p,size,1,f);
      p[size]=0;
      fclose(f);
      
            I=ObjectMapReadGRDStr(obj,buffer,size,state);

      mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/
ObjectMap *ObjectMapLoadXPLORFile(ObjectMap *obj,char *fname,int state,int is_file)
{
  ObjectMap *I = NULL;
  int ok=true;
  FILE *f = NULL;
  long size;
  char *buffer,*p;
  float mat[9];

  if(is_file) {
    f=fopen(fname,"rb");
    if(!f)
      ok=ErrMessage("ObjectMapLoadXPLORFile","Unable to open file!");
  }
  if(ok)
       {
      if(Feedback(FB_ObjectMap,FB_Actions))
        {
          if(is_file) {
            printf(" ObjectMapLoadXPLORFile: Loading from '%s'.\n",fname);
          } else {
            printf(" ObjectMapLoadXPLORFile: Loading...\n");
          }
          
        }
      
      if(is_file) {
        fseek(f,0,SEEK_END);
        size=ftell(f);
        fseek(f,0,SEEK_SET);
        
        buffer=(char*)mmalloc(size+255);
        ErrChkPtr(buffer);
        p=buffer;
        fseek(f,0,SEEK_SET);
        fread(p,size,1,f);
        p[size]=0;
        fclose(f);
      } else {
        buffer = fname;
      }
      
            I=ObjectMapReadXPLORStr(obj,buffer,state);

      if(is_file) 
        mfree(buffer);
      if(state<0)
        state=I->NState-1;
      if(state<I->NState) {
        ObjectMapState *ms;
        ms = &I->State[state];
        if(ms->Active) {
          CrystalDump(ms->Crystal);
          multiply33f33f(ms->Crystal->FracToReal,ms->Crystal->RealToFrac,mat);
        }
      }
    }
  return(I);

}
/*========================================================================*/
int ObjectMapSetBorder(ObjectMap *I,float level)
{
  int a;
  int result=false;
  for(a=0;a<I->NState;a++) {
    if(I->State[a].Active)
      result = result && ObjectMapStateSetBorder(&I->State[a],level);
  }
  return(result);
}
/*========================================================================*/
int ObjectMapNumPyArrayToMapState(ObjectMapState *ms,PyObject *ary) {
  
  int a,b,c,d,e;
  float v[3],dens,maxd,mind;
  int ok = true;
#ifdef _PYMOL_NUMPY
  MyArrayObject *pao;
#endif


#ifdef _PYMOL_NUMPY
  pao = (MyArrayObject*)ary;
#endif
  maxd = FLT_MIN;
  mind = FLT_MAX;
  if(ok) {
    ms->FDim[0]=ms->Dim[0];
    ms->FDim[1]=ms->Dim[1];
    ms->FDim[2]=ms->Dim[2];
    ms->FDim[3]=3;

    if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
      ok=false;
    else {
      ms->Field=IsosurfFieldAlloc(ms->FDim);
      for(c=0;c<ms->FDim[2];c++)
        {
          v[2]=ms->Origin[2]+ms->Grid[2]*c;
          for(b=0;b<ms->FDim[1];b++) {
            v[1]=ms->Origin[1]+ms->Grid[1]*b;
            for(a=0;a<ms->FDim[0];a++) {
              v[0]=ms->Origin[0]+ms->Grid[0]*a;
#ifdef _PYMOL_NUMPY
              dens = (float)(*((double*)
                (pao->data+
                 (pao->strides[0]*a)+
                 (pao->strides[1]*b)+
                 (pao->strides[2]*c))));
#else
              dens = 0.0;
#endif
              F3(ms->Field->data,a,b,c) = dens;
              if(maxd<dens) maxd = dens;
              if(mind>dens) mind = dens;
              for(e=0;e<3;e++) 
                F4(ms->Field->points,a,b,c,e) = v[e];
            }
          }
        }
      d = 0;
      for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
        {
          v[2]=ms->Origin[2]+ms->Grid[2]*c;
          for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
            v[1]=ms->Origin[1]+ms->Grid[1]*b;
            for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
              v[0]=ms->Origin[0]+ms->Grid[0]*a;
              copy3f(v,ms->Corner[d]);
              d++;
            }
          }
        }
    }
  }
  if(ok) {
    copy3f(ms->Origin,ms->ExtentMin);
    copy3f(ms->Origin,ms->ExtentMax);
    add3f(ms->Range,ms->ExtentMax,ms->ExtentMax);
  }
  if(!ok) {
    ErrMessage("ObjectMap","Error reading map");
  } else {
    ms->Active=true;
    if(Feedback(FB_ObjectMap,FB_Actions)) {
      printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
    }
  }
  return(ok);
}
/*========================================================================*/
ObjectMap *ObjectMapLoadChemPyBrick(ObjectMap *I,PyObject *Map,
                                           int state,int discrete)
{
  int ok=true;
  int isNew = true;
  PyObject *tmp;
  ObjectMapState *ms;


  if(!I) 
       isNew=true;
  else 
       isNew=false;

  if(ok) {

       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }

    if(state<0) state=I->NState;
    if(I->NState<=state) {
      VLACheck(I->State,ObjectMapState,state);
      I->NState=state+1;
    }
    ms=&I->State[state];
    ObjectMapStateInit(ms);
    
    if(PyObject_HasAttrString(Map,"origin")&&
       PyObject_HasAttrString(Map,"dim")&&
       PyObject_HasAttrString(Map,"range")&&
       PyObject_HasAttrString(Map,"grid")&&
       PyObject_HasAttrString(Map,"lvl"))
      {
        tmp = PyObject_GetAttrString(Map,"origin");
        if(tmp) {
          PConvPyListToFloatArray(tmp,&ms->Origin);
          Py_DECREF(tmp);
        } else 
          ok=ErrMessage("ObjectMap","missing brick origin.");
        tmp = PyObject_GetAttrString(Map,"dim");
        if(tmp) {
          PConvPyListToIntArray(tmp,&ms->Dim);
          Py_DECREF(tmp);
        } else 
          ok=ErrMessage("ObjectMap","missing brick dimension.");
        tmp = PyObject_GetAttrString(Map,"range");
        if(tmp) {
          PConvPyListToFloatArray(tmp,&ms->Range);
          Py_DECREF(tmp);
        } else 
          ok=ErrMessage("ObjectMap","missing brick range.");
        tmp = PyObject_GetAttrString(Map,"grid");
        if(tmp) {
          PConvPyListToFloatArray(tmp,&ms->Grid);
          Py_DECREF(tmp);
        } else
          ok=ErrMessage("ObjectMap","missing brick grid.");
        tmp = PyObject_GetAttrString(Map,"lvl");
        if(tmp) {

          ObjectMapNumPyArrayToMapState(ms,tmp);       

          Py_DECREF(tmp);
        } else
          ok=ErrMessage("ObjectMap","missing brick density.");

      }
    SceneChanged();
    SceneCountFrames();
    if(ok) {
      ms->Active=true;
      ObjectMapUpdateExtents(I);
          
    }
  }
  return(I);
}

/*========================================================================*/
ObjectMap *ObjectMapLoadChemPyMap(ObjectMap *I,PyObject *Map,
                                  int state,int discrete)
{

  int ok=true;
  int isNew = true;
  float *cobj;
  WordType format;
  float v[3],vr[3],dens,maxd,mind;
  int a,b,c,d,e;
  ObjectMapState *ms;



  /*  
  double test[1000];
  for(a=0;a<1000;a++) {
      test[a]=rand()/(1.0+INT_MAX);
      }
      PyObject_SetAttrString(Map,"c_object",
      PyCObject_FromVoidPtr(test,NULL));
  */
  maxd = FLT_MIN;
  mind = FLT_MAX;

  if(!I) 
       isNew=true;
  else 
       isNew=false;

  if(ok) {

       if(isNew) {
            I=(ObjectMap*)ObjectMapNew();
            isNew = true;
       } else {
            isNew = false;
       }
    
    if(state<0) state=I->NState;
    if(I->NState<=state) {
      VLACheck(I->State,ObjectMapState,state);
      I->NState=state+1;
    }
    ms=&I->State[state];
    ObjectMapStateInit(ms);

    if(!PConvAttrToStrMaxLen(Map,"format",format,sizeof(WordType)-1))
      ok=ErrMessage("LoadChemPyMap","bad 'format' parameter.");
    else if(!PConvAttrToFloatArrayInPlace(Map,"cell_dim",ms->Crystal->Dim,3))
      ok=ErrMessage("LoadChemPyMap","bad 'cell_dim' parameter.");
    else if(!PConvAttrToFloatArrayInPlace(Map,"cell_ang",ms->Crystal->Angle,3))
      ok=ErrMessage("LoadChemPyMap","bad 'cell_ang' parameter.");
    else if(!PConvAttrToIntArrayInPlace(Map,"cell_div",ms->Div,3))
      ok=ErrMessage("LoadChemPyMap","bad 'cell_div' parameter.");
    else if(!PConvAttrToIntArrayInPlace(Map,"first",ms->Min,3))
      ok=ErrMessage("LoadChemPyMap","bad 'first' parameter.");
    else if(!PConvAttrToIntArrayInPlace(Map,"last",ms->Max,3))
      ok=ErrMessage("LoadChemPyMap","bad 'last' parameter.");

    if(ok) {
      if (strcmp(format,"CObjectZYXfloat")==0) {
        ok = PConvAttrToPtr(Map,"c_object",(void**)&cobj);
        if(!ok)
          ErrMessage("LoadChemPyMap","CObject unreadable.");        
      } else {
        ok=ErrMessage("LoadChemPyMap","unsupported format.");        
      }
    }
    /* good to go */

    if(ok) {
      if (strcmp(format,"CObjectZYXfloat")==0) {

        ms->FDim[0]=ms->Max[0]-ms->Min[0]+1;
        ms->FDim[1]=ms->Max[1]-ms->Min[1]+1;
        ms->FDim[2]=ms->Max[2]-ms->Min[2]+1;
        if(Feedback(FB_ObjectMap,FB_Actions)) {
          printf(" LoadChemPyMap: CObjectZYXdouble %dx%dx%d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]);        
        }
        ms->FDim[3]=3;
        if(!(ms->FDim[0]&&ms->FDim[1]&&ms->FDim[2])) 
          ok=false;
        else {
          CrystalUpdate(ms->Crystal);
          ms->Field=IsosurfFieldAlloc(ms->FDim);
          for(c=0;c<ms->FDim[2];c++)
            {
              v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
              for(b=0;b<ms->FDim[1];b++) {
                v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
                for(a=0;a<ms->FDim[0];a++) {
                  v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                  
                  dens = *(cobj++);

                  F3(ms->Field->data,a,b,c) = dens;
                  if(maxd<dens) maxd = dens;
                  if(mind>dens) mind = dens;
                  transform33f3f(ms->Crystal->FracToReal,v,vr);
                  for(e=0;e<3;e++) 
                    F4(ms->Field->points,a,b,c,e) = vr[e];
                }
              }
            }

          if(ok) {
            d = 0;
            for(c=0;c<ms->FDim[2];c+=(ms->FDim[2]-1))
              {
                v[2]=(c+ms->Min[2])/((float)ms->Div[2]);
                for(b=0;b<ms->FDim[1];b+=(ms->FDim[1]-1)) {
                  v[1]=(b+ms->Min[1])/((float)ms->Div[1]);
                  for(a=0;a<ms->FDim[0];a+=(ms->FDim[0]-1)) {
                    v[0]=(a+ms->Min[0])/((float)ms->Div[0]);
                    transform33f3f(ms->Crystal->FracToReal,v,vr);
                    copy3f(vr,ms->Corner[d]);
                    d++;
                  }
                }
              }
          }
        }
      }
    }
    
    if(ok) {
      CrystalDump(ms->Crystal);
      
      v[2]=(ms->Min[2])/((float)ms->Div[2]);
      v[1]=(ms->Min[1])/((float)ms->Div[1]);
      v[0]=(ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMin);
      
      v[2]=((ms->FDim[2]-1)+ms->Min[2])/((float)ms->Div[2]);
      v[1]=((ms->FDim[1]-1)+ms->Min[1])/((float)ms->Div[1]);
      v[0]=((ms->FDim[0]-1)+ms->Min[0])/((float)ms->Div[0]);
      
      transform33f3f(ms->Crystal->FracToReal,v,ms->ExtentMax);

    }

    if(!ok) {
      ErrMessage("ObjectMap","Error reading map");
    } else {
      ms->Active=true;
      ObjectMapUpdateExtents(I);
            if(Feedback(FB_ObjectMap,FB_Actions)) {
        printf(" ObjectMap: Map Read.  Range = %5.3f to %5.3f\n",mind,maxd);
      }
    }

    if(ok) {
      SceneChanged();
      SceneCountFrames();
    }
  }
  return(I);
}