/***************************************************************************//**
**  \mainpage Monte Carlo eXtreme - GPU accelerated Monte Carlo Photon Migration \
**      -- OpenCL edition
**  \author Qianqian Fang <q.fang at neu.edu>
**  \copyright Qianqian Fang, 2009-2018
**
**  \section sref Reference:
**  \li \c (\b Yu2018) Leiming Yu, Fanny Nina-Paravecino, David Kaeli, and Qianqian Fang,
**         "Scalable and massively parallel Monte Carlo photon transport simulations 
**         for heterogeneous computing platforms," J. Biomed. Optics, 23(1), 010504 (2018)
**
**  \section slicense License
**          GPL v3, see LICENSE.txt for details
*******************************************************************************/

/***************************************************************************//**
\file    mcx_shapes.c
\brief   3D shape parsing and rasterization unit

In this unit, we load and parse a JSON-formatted shape file,
rasterize the 3D objects and subsequently add to a voxelated volume.
*******************************************************************************/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#include "cjson/cJSON.h"
#include "mcx_shapes.h"

#define MIN(a,b)           ((a)<(b)?(a):(b))
#define MAX(a,b)           ((a)>(b)?(a):(b))

const char *ShapeTags[]={"Name","Origin","Grid","Subgrid","Sphere","Box","XSlabs",
                         "YSlabs","ZSlabs","XLayers","YLayers","ZLayers",
			 "Cylinder","UpperSpace",NULL};
int (*Rasterizers[])(cJSON *obj, Grid3D *g)={NULL,mcx_raster_origin,mcx_raster_grid,mcx_raster_subgrid,
                          mcx_raster_sphere,mcx_raster_box,mcx_raster_slabs,mcx_raster_slabs,
			  mcx_raster_slabs, mcx_raster_layers,mcx_raster_layers,
			  mcx_raster_layers,mcx_raster_cylinder,mcx_raster_upperspace,NULL};
char ErrorMsg[MAX_SHAPE_ERR]={'\0'};


/*******************************************************************************/
/*! \fn int mcx_load_jsonshapes(Grid3D *g, char *fname)

    \brief Load a JSON-formatted shape file and process
    \param g A structure pointing to the volume and dimension data
    \param fname The file name string to the JSON shape file
*/

int mcx_load_jsonshapes(Grid3D *g, char *fname){
    FILE *fp=fopen(fname,"rb");
    if(fp==NULL){
        sprintf(ErrorMsg,"Can not read the JSON file");
        return -2;
    }
    if(g==NULL){
        sprintf(ErrorMsg,"The background grid is not initialized");
        return -3;
    }
    if(strstr(fname,".json")!=NULL){
    	char *jbuf;
    	int len,err;

    	fseek (fp, 0, SEEK_END);
    	len=ftell(fp)+1;
    	jbuf=(char *)malloc(len);
    	rewind(fp);
    	if(fread(jbuf,len-1,1,fp)!=1){
            sprintf(ErrorMsg,"Failed when reading a JSON file from %s",fname);
    	    return -1;
	}
    	jbuf[len-1]='\0';
        fclose(fp);
    	if((err=mcx_parse_shapestring(g,jbuf))){ /*error msg is generated inside*/
	   free(jbuf);
	   return err;
	}
    	free(jbuf);
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_parse_shapestring(Grid3D *g, char *shapedata)

    \brief Load JSON-formatted shape definitions from a string
    \param g A structure pointing to the volume and dimension data
    \param shapedata A string containg the JSON shape data
*/

int mcx_parse_shapestring(Grid3D *g, char *shapedata){
    if(g && shapedata){
    	cJSON *jroot = cJSON_Parse(shapedata);
    	if(jroot){
	    int err;
	    if(g && g->vol && *(g->vol))
	        *(g->vol)=(unsigned int*)realloc((void *)*(g->vol),g->dim->x*g->dim->y*g->dim->z*sizeof(int));
	    else
	        *(g->vol)=(unsigned int*)calloc(g->dim->x*g->dim->y,g->dim->z*sizeof(int));

    	    if((err=mcx_parse_jsonshapes(jroot,g))) /*error msg is generated inside*/
	       return err;

    	    cJSON_Delete(jroot);
    	}else{
    	    char *ptrold, *ptr=(char*)cJSON_GetErrorPtr();
    	    if(ptr) ptrold=strstr(shapedata,ptr);

    	    if(ptr && ptrold){
    	       char *offs=(ptrold-shapedata>=50) ? ptrold-50 : shapedata;
    	       while(offs<ptrold){
    		  fprintf(stderr,"%c",*offs);
    		  offs++;
    	       }
    	       fprintf(stderr,"<error>%.50s\n",ptrold);
    	    }
            sprintf(ErrorMsg,"Invalid JSON file");
    	    return -2;
    	}
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_parse_jsonshapes(cJSON *root, Grid3D *g)

    \brief Parse a JSON-formatted shape file and rasterize the objects to the volume
    \param root A cJSON pointer points to the root obj of the shape block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_parse_jsonshapes(cJSON *root, Grid3D *g){
    cJSON *shapes;
    int id,objcount=1;
    int (*raster)(cJSON *obj, Grid3D *g)=NULL;

    shapes  = cJSON_GetObjectItem(root,"Shapes");
    if(shapes){
    	shapes=shapes->child;
    	while(shapes && shapes->child){
            id=mcx_find_shapeid(shapes->child->string);
            if(id>=0 && id<sizeof(ShapeTags)/sizeof(char*))
        	raster=Rasterizers[id];
            else{
                sprintf(ErrorMsg,"The #%d element in the Shapes section has an undefined tag %s",
		    objcount,shapes->child->string);
        	return -(objcount+100);
	    }
            if(raster){
               id=raster(shapes->child,g);
               if(id)
        	   return id;
            }
            objcount++;
            shapes=shapes->next;
        }
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_origin(cJSON *obj, Grid3D *g)

    \brief Reset the origin of the domain, default is [0,0,0]
    \param obj A cJSON pointer points to the Origin obj block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_origin(cJSON *obj, Grid3D *g){
    if(obj && cJSON_GetArraySize(obj)==3){
       g->orig.x=obj->child->valuedouble;
       g->orig.y=obj->child->next->valuedouble;
       g->orig.z=obj->child->next->next->valuedouble;
    }else{
       sprintf(ErrorMsg,"An Origin record does not contain a triplet");
       return 1;
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_sphere(cJSON *obj, Grid3D *g)

    \brief Rasterize a 3D sphere and add to the volume
    \param obj A cJSON pointer points to the sphere obj block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_sphere(cJSON *obj, Grid3D *g){
    float O[3],R,R2,dx,dy,dz;
    int i,j,k,dimxy,dimyz,tag=0;

    cJSON *val=cJSON_GetObjectItem(obj,"O");
    if(val && cJSON_GetArraySize(val)==3){
       O[0]=val->child->valuedouble;
       O[1]=val->child->next->valuedouble;
       O[2]=val->child->next->next->valuedouble;
    }else{
       sprintf(ErrorMsg,"A Sphere command misses O field");
       return 1;
    }

    val=cJSON_GetObjectItem(obj,"R");
    if(val){
       R=val->valuedouble;
    }else{
       sprintf(ErrorMsg,"A Sphere command misses R field");
       return 2;
    }

    val=cJSON_GetObjectItem(obj,"Tag");
    if(val){
       tag=val->valueint;
    }

    R2=R*R;
    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    for(k=0;k<g->dim->z;k++){
      dz=(k+0.5f)-O[2];
      for(j=0;j<g->dim->y;j++){
    	dy=(j+0.5f)-O[1];
    	for(i=0;i<g->dim->x;i++){
               dx=(i+0.5f)-O[0];
               if(dx*dx+dy*dy+dz*dz<=R2)
        	   (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
        }
      }
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_subgrid(cJSON *obj, Grid3D *g)

    \brief Rasterize a 3D rectangular region and add to the volume
    \param obj A cJSON pointer points to the rectangular obj block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_subgrid(cJSON *obj, Grid3D *g){
    int O[3]={0},S[3]={0};
    int i,j,k,dimxy,dimyz,tag=0;

    cJSON *val=cJSON_GetObjectItem(obj,"O");
    if(val && cJSON_GetArraySize(val)==3){
       O[0]=val->child->valueint-1-g->orig.x;  /*Subgrid: numbers start from 1*/
       O[1]=val->child->next->valueint-1-g->orig.y;
       O[2]=val->child->next->next->valueint-1-g->orig.z;
    }else{
       sprintf(ErrorMsg,"A Subgrid command misses O field");
       return 1;
    }
    val=cJSON_GetObjectItem(obj,"Size");
    if(val && cJSON_GetArraySize(val)==3){
       S[0]=val->child->valueint;
       S[1]=val->child->next->valueint;
       S[2]=val->child->next->next->valueint;
    }else{
       sprintf(ErrorMsg,"A Box command misses Size field");
       return 2;
    }

    val=cJSON_GetObjectItem(obj,"Tag");
    if(val)
       tag=val->valueint;

    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    for(k=0;k<g->dim->z;k++){
      if(k<O[2] || k>O[2]+S[2]) continue;
      for(j=0;j<g->dim->y;j++){
        if(j<O[1] || j>O[1]+S[1]) continue;
    	for(i=0;i<g->dim->x;i++){
            if(i<O[0] || i>O[0]+S[0]) continue;
            (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
        }
      }
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_box(cJSON *obj, Grid3D *g)

    \brief Rasterize a 3D rectangular region and add to the volume
    \param obj A cJSON pointer points to the rectangular obj block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_box(cJSON *obj, Grid3D *g){
    float O[3]={0.f},S[3]={0.f},dx,dy,dz;
    int i,j,k,dimxy,dimyz,tag=0;

    cJSON *val=cJSON_GetObjectItem(obj,"O");
    if(val && cJSON_GetArraySize(val)==3){
       O[0]=val->child->valuedouble-g->orig.x;
       O[1]=val->child->next->valuedouble-g->orig.y;
       O[2]=val->child->next->next->valuedouble-g->orig.z;
    }else{
       sprintf(ErrorMsg,"A Box command misses O field");
       return 1;
    }
    val=cJSON_GetObjectItem(obj,"Size");
    if(val && cJSON_GetArraySize(val)==3){
       S[0]=val->child->valuedouble;
       S[1]=val->child->next->valuedouble;
       S[2]=val->child->next->next->valuedouble;
    }else{
       sprintf(ErrorMsg,"A Box command misses Size field");
       return 2;
    }

    val=cJSON_GetObjectItem(obj,"Tag");
    if(val)
       tag=val->valueint;

    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    for(k=0;k<g->dim->z;k++){
      dz=(k+0.5f);
      if(dz<O[2] || dz>O[2]+S[2]) continue;
      for(j=0;j<g->dim->y;j++){
    	dy=(j+0.5f);
        if(dy<O[1] || dy>O[1]+S[1]) continue;
    	for(i=0;i<g->dim->x;i++){
            dx=(i+0.5f);
            if(dx<O[0] || dx>O[0]+S[0]) continue;
            (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
        }
      }
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_cylinder(cJSON *obj, Grid3D *g)

    \brief Rasterize a finite 3D cylindrical region and add to the volume
    \param obj A cJSON pointer points to the cylindrical obj block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_cylinder(cJSON *obj, Grid3D *g){
    float C0[3],C1[3],v[3],d0,R,R2,d,dx,dy,dz;
    int i,j,k,dimxy,dimyz,tag=0;

    cJSON *val=cJSON_GetObjectItem(obj,"C0");
    if(val && cJSON_GetArraySize(val)==3){
       C0[0]=val->child->valuedouble-g->orig.x;
       C0[1]=val->child->next->valuedouble-g->orig.y;
       C0[2]=val->child->next->next->valuedouble-g->orig.z;
    }else{
       sprintf(ErrorMsg,"A Cylinder command misses C0 field");
       return 1;
    }

    val=cJSON_GetObjectItem(obj,"C1");
    if(val && cJSON_GetArraySize(val)==3){
       C1[0]=val->child->valuedouble-g->orig.x;
       C1[1]=val->child->next->valuedouble-g->orig.y;
       C1[2]=val->child->next->next->valuedouble-g->orig.z;
    }else{
       sprintf(ErrorMsg,"A Cylinder command misses C1 field");
       return 1;
    }
    v[0]=C1[0]-C0[0];
    v[1]=C1[1]-C0[1];
    v[2]=C1[2]-C0[2];
    d0=sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
    if(d0==0.f){
       sprintf(ErrorMsg,"Coincident end points in the definition of Cylinder command");
       return 1;
    }
    v[0]/=d0;v[1]/=d0;v[2]/=d0;

    val=cJSON_GetObjectItem(obj,"R");
    if(val){
       R=val->valuedouble;
    }else{
       sprintf(ErrorMsg,"A Sphere command misses R field");
       return 2;
    }

    val=cJSON_GetObjectItem(obj,"Tag");
    if(val){
       tag=val->valueint;
    }

    R2=R*R;
    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    for(k=0;k<g->dim->z;k++){
      dz=(k+0.5f)-C0[2];
      for(j=0;j<g->dim->y;j++){
    	dy=(j+0.5f)-C0[1];
    	for(i=0;i<g->dim->x;i++){
            dx=(i+0.5f)-C0[0];
	    d=v[0]*dx+v[1]*dy+v[2]*dz; /* (|PC0|*cos(theta)) */
	    if(d>d0 || d<0.f) continue;
	    d=dx*dx+dy*dy+dz*dz-d*d;   /* (|PC0|*sin(theta))^2 */
            if(d<=R2)
             	(*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
        }
      }
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_slabs(cJSON *obj, Grid3D *g)

    \brief Rasterize a 3D layered-slab structure and add to the volume
    \param obj A cJSON pointer points to the layered-slab structure block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_slabs(cJSON *obj, Grid3D *g){
    float *bd=NULL;
    int i,j,k,dimxy,dimyz,tag=0,num=0,num2,dir=-1,p;
    cJSON *item, *val;

    if(strcmp(obj->string,"XSlabs")==0)      dir=0;
    else if(strcmp(obj->string,"YSlabs")==0) dir=1;
    else if(strcmp(obj->string,"ZSlabs")==0) dir=2;
    else{
        sprintf(ErrorMsg,"Unsupported layer command");
        return 1;
    }
    val=cJSON_GetObjectItem(obj,"Bound");
    if(val){
       num=cJSON_GetArraySize(val);
       if(num==0) return 0;
       if(num==2 && obj->child->type!=cJSON_Array){
          item=obj;
	  bd=(float *)malloc(cJSON_GetArraySize(obj)*sizeof(float));
	  num=1;
       }else{
          item=obj->child;
          bd=(float *)malloc(cJSON_GetArraySize(obj)*2*sizeof(float));
       }
       for(i=0;i<num;i++){
           if(cJSON_GetArraySize(item)!=2){
              sprintf(ErrorMsg,"The Bound field must contain number pairs");
	      return 2;
	   }
	   bd[i<<1]    =MAX(item->child->valuedouble-0.5,0.0); /*inclusive of both ends*/
	   bd[(i<<1)+1]=MIN(item->child->next->valuedouble,(&(g->dim->x))[dir]);
	   if(bd[(i<<1)+1] < bd[i<<1]){
	      float tmp=bd[(i<<1)+1];
	      bd[(i<<1)+1]=bd[i<<1];
	      bd[i<<1]=tmp;
	   }
           item=item->next;
       }
    }else{
       sprintf(ErrorMsg,"A %s command misses Bound field",obj->string);
       return 1;
    }

    if(num==0) return 0;
    num2=num<<1;

    val=cJSON_GetObjectItem(obj,"Tag");
    if(val){
       tag=val->valueint;
    }

    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    if(dir==0){
        for(p=0;p<num2;p+=2)
	  for(k=0;k<g->dim->z;k++)
	    for(j=0;j<g->dim->y;j++)
	       for(i=(int)(bd[p]);i<(int)(bd[p+1]);i++)
                  (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
    }else if(dir==1){
        for(p=0;p<num2;p+=2)
	  for(k=0;k<g->dim->z;k++)
	    for(j=(int)(bd[p]);j<(int)(bd[p+1]);j++)
	       for(i=0;i<g->dim->x;i++)
                  (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
    }else if(dir==2){
        for(p=0;p<num2;p+=2)
	  for(k=(int)(bd[p]);k<(int)(bd[p+1]);k++)
	    for(j=0;j<g->dim->y;j++)
	       for(i=0;i<g->dim->x;i++)
                  (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
    }
    if(bd) free(bd);
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_layers(cJSON *obj, Grid3D *g)

    \brief Rasterize a 3D layer structure and add to the volume
    \param obj A cJSON pointer points to the layer structure block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_layers(cJSON *obj, Grid3D *g){
    int *bd=NULL;
    int i,j,k,dimxy,dimyz,num=0,num3,dir=-1,p;
    cJSON *item;

    if(strcmp(obj->string,"XLayers")==0)      dir=0;
    else if(strcmp(obj->string,"YLayers")==0) dir=1;
    else if(strcmp(obj->string,"ZLayers")==0) dir=2;
    else{
        sprintf(ErrorMsg,"Unsupported command %s",obj->string);
        return 1;
    }
    if(obj && obj->type==cJSON_Array){
       num=cJSON_GetArraySize(obj);
       if(num==0) return 0;
       if(num==3 && obj->child->type!=cJSON_Array){
          item=obj;
	  bd=(int *)malloc(cJSON_GetArraySize(obj)*sizeof(int));
	  num=1;
       }else{
          item=obj->child;
          bd=(int *)malloc(cJSON_GetArraySize(obj)*3*sizeof(int));
       }
       for(i=0;i<num;i++){
           if(cJSON_GetArraySize(item)!=3){
              sprintf(ErrorMsg,"The %s must contain integer triplets",obj->string);
	      return 2;
	   }
	   bd[i*3]  =MAX(item->child->valueint,1)-1; /*inclusive of both ends*/
	   bd[i*3+1]=MIN(item->child->next->valueint,(&(g->dim->x))[dir]);
	   bd[i*3+2]=item->child->next->next->valueint;
	   if(bd[i*3+1] < bd[i*3]){
	      float tmp=bd[i*3+1];
	      bd[i*3+1]=bd[i*3];
	      bd[i*3]=tmp;
	   }
           item=item->next;
       }
    }else{
       sprintf(ErrorMsg,"A %s object must be an array",obj->string);
       return 1;
    }

    if(num==0) return 0;
    num3=num*3;

    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    if(dir==0){
        for(p=0;p<num3;p+=3)
	  for(k=0;k<g->dim->z;k++)
	    for(j=0;j<g->dim->y;j++)
	       for(i=bd[p];i<bd[p+1];i++)
                  (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=bd[p+2];
    }else if(dir==1){
        for(p=0;p<num3;p+=3)
	  for(k=0;k<g->dim->z;k++)
	    for(j=bd[p];j<bd[p+1];j++)
	       for(i=0;i<g->dim->x;i++)
                  (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=bd[p+2];
    }else if(dir==2){
        for(p=0;p<num3;p+=3)
	  for(k=bd[p];k<bd[p+1];k++)
	    for(j=0;j<g->dim->y;j++)
	       for(i=0;i<g->dim->x;i++)
                  (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=bd[p+2];
    }
    if(bd) free(bd);
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_upperspace(cJSON *obj, Grid3D *g)

    \brief Rasterize a 3D semi-space region and add to the volume
    \param obj A cJSON pointer points to the semi-space object block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_upperspace(cJSON *obj, Grid3D *g){
    float C[4],dx,dy,dz;
    int i,j,k,dimxy,dimyz,tag=0;

    cJSON *val=cJSON_GetObjectItem(obj,"Coef");
    if(val && cJSON_GetArraySize(val)==4){
       C[0]=val->child->valuedouble;
       C[1]=val->child->next->valuedouble;
       C[2]=val->child->next->next->valuedouble;
       C[3]=val->child->next->next->next->valuedouble;
    }else{
       sprintf(ErrorMsg,"An UpperSpace command misses Coef field");
       return 1;
    }
    val=cJSON_GetObjectItem(obj,"Tag");
    if(val){
       tag=val->valueint;
    }
    dimxy=g->dim->x*g->dim->y;
    dimyz=g->dim->y*g->dim->z;

    for(k=0;k<g->dim->z;k++){
      dz=(k+0.5f);
      for(j=0;j<g->dim->y;j++){
    	dy=(j+0.5f);
    	for(i=0;i<g->dim->x;i++){
               dx=(i+0.5f);
               if(C[0]*dx+C[1]*dy+C[2]*dz>C[3])
        	   (*(g->vol))[g->rowmajor? i*dimyz+j*g->dim->z+k : k*dimxy+j*g->dim->x+i]=tag;
        }
      }
    } 
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_raster_grid(cJSON *obj, Grid3D *g)

    \brief Recreate the background grid with a different dimension or medium
    \param obj A cJSON pointer points to the background grid block
    \param g  A structure pointing to the volume and dimension data
*/

int mcx_raster_grid(cJSON *obj, Grid3D *g){
    int dimxy,dimxyz;
    unsigned int tag=0;

    cJSON *val=cJSON_GetObjectItem(obj,"Size");
    if(val && cJSON_GetArraySize(val)==3){
       g->dim->x=val->child->valuedouble;
       g->dim->y=val->child->next->valuedouble;
       g->dim->z=val->child->next->next->valuedouble;
       dimxy=g->dim->x*g->dim->y;
       dimxyz=dimxy*g->dim->z;
       if(dimxy*g->dim->z>0){
    	  if(g->vol && *g->vol) free(*(g->vol));
          *(g->vol)=(unsigned int *)calloc(dimxy*sizeof(int),g->dim->z);
       }else
          *(g->vol)=NULL;
    }else{
       sprintf(ErrorMsg,"A Grid command misses Size field");
       return 1;
    }

    val=cJSON_GetObjectItem(obj,"Tag");
    if(val){
       int i;
       tag=val->valueint;
       if(g->vol && *(g->vol))
          for(i=0;i<dimxyz;i++)
              (*(g->vol))[i]=tag;
    }
    return 0;
}

/*******************************************************************************/
/*! \fn int mcx_find_shapeid(char *shapename)

    \brief Look up the JSON object tag and return the index to the processing function
    \param shapename The string of the JSON shape object
*/

int mcx_find_shapeid(char *shapename){
    int i=0;
    while(ShapeTags[i]){
        if(strcmp(shapename,ShapeTags[i])==0){
                return i;
        }
        i++;
    }
    return -1;
}

/*******************************************************************************/
/*! \fn char * mcx_last_shapeerror()

    \brief return the last error message encountered in the processing
*/

char * mcx_last_shapeerror(){
    return ErrorMsg;
}
