#include <stdio.h>
 
#include "TH2F.h"
 
#include "TCanvas.h"
 
#include "TMath.h"
 
#include "TStyle.h"
 
#include "TColor.h"
 
 
 
void mypalette(){
 
 
 
   Int_t MyPalette[1000];
 
   double stops[] = {0.00, 0.34, 0.61, 0.84, 1.00};
 
   double red[]   = {1.00, 0.84, 0.61, 0.34, 0.00};
 
   double green[] = {1.00, 0.84, 0.61, 0.34, 0.00};
 
   double blue[]  = {1.00, 0.84, 0.61, 0.34, 0.00};
 
 
 
   
 
   Int_t gs = TColor::CreateGradientColorTable(5, stops, red, green, blue, 999);
 
   //for (int i=0;i<999;i++) MyPalette[i] = gs+998-i;
 
 
 
   for (int i=0;i<999;i++) MyPalette[i] = gs+i;
 
   gStyle->SetPalette(999, MyPalette);
 
}
 
 
 
int imageread(char *fname,int nchx=28, int nchy=56, int nchz=0, int format=0, int size=0, char *opt=""){
 
FILE *fp= fopen(fname,"read");
 
const int maxbuf=360000;
 
float *fbuf= new float[maxbuf];
 
unsigned short*ibuf=((unsigned short*) &fbuf[0]);
 
if (size==0) {
 
 if (format) size=sizeof(float);
 
 else size=sizeof(unsigned short);
 
}
 
int ncount=0;
 
 
 
if(fp){
 
 
 
gStyle->SetPalette(100);
 
mypalette();
 
TCanvas *c=new TCanvas("c","c",1500,1500);
 
int sizx=sqrt(nchz);
 
c->Divide(sizx,TMath::Ceil(nchz/double(sizx)) );
 
char hname[128];
 
 
 
while (!feof(fp)){
 
 
 
 
 
int nb=fread(fbuf,size,nchx*nchy,fp);
 
printf("nb=%d\n",nb);
 
if (nb>0){
 
  sprintf(hname,"slice %d",ncount);
 
  TH2F *h=new TH2F("h",hname,nchx,-0.5,nchx-0.5,nchy,-0.5,nchy-0.5);
 
  float val=0;
 
  for (int i=0 ; i<nb;i++){
 
    //h->SetBinContent(i%nchx+1, i/nchx+1, fbuf[i]);
 
    if (format) val=fbuf[i]; else val=float(ibuf[i]);
 
    h->SetBinContent(i/nchy+1,i%nchy+1, val);
 
  } 
 
  c->cd(ncount+1);
 
  h->SetMinimum(0);
 
  h->Draw(opt);
 
}
 
 
 
ncount++;
 
}
 
fclose(fp);
 
c->Modified();
 
c->Update();
 
}
 
 
 
delete fbuf;
 
return 0;
 
}
 
 
 
int interfileread(char *fname, char *opt=""){
 
  FILE *fp=fopen(fname,"r");
 
  if (!fp) {
 
    printf("Cannot open %s\n",fname);
 
    return 0;
 
  }  
 
  int j=0;   
 
  int ndim=400;
 
  char *line= new char[ndim];                
 
  int indx[4]={-1,-1,-1,-1};  
 
  int dim[5]; 
 
  char *cmd,*val,*token;
 
  char fimage[100];                                                                                      
 
  while (fgets(line,ndim,fp)!=NULL) {   
 
     char *v= strstr(line,":=");
 
     if (v!=NULL) {
 
       cmd=line;
 
       v[0]=0;
 
       val = &v[2];
 
       val[strlen(val)-1]=0;
 
       printf("%d #%s#  %s\n",j++, val, cmd);
 
       int axis=-1;
 
       if (strstr(cmd,"matrix axis label")!= NULL){
 
         
 
         sscanf(cmd,"matrix axis label [%d]", &axis);
 
         if(axis>0 && axis<5){
 
           if (strstr(val," x")!= NULL)  indx[axis-1]=1;
 
           if (strstr(val," y")!= NULL)  indx[axis-1]=0;
 
           if (strstr(val," z")!= NULL)  indx[axis-1]=3;
 
 
 
           if (strstr(val,"view")      != NULL)  indx[axis-1]=0;           
 
           if (strstr(val,"tangential")!= NULL)  indx[axis-1]=1;
 
           if (strstr(val,"segment")   != NULL)  indx[axis-1]=2;
 
           if (strstr(val,"axial")     != NULL)  indx[axis-1]=3;
 
           printf("-->> %s axis=%d indx=%d\n",val, axis,indx[axis-1]);
 
           
 
         }
 
       } 
 
       if (strstr(cmd,"!matrix size")!= NULL){
 
           sscanf(cmd,"!matrix size [%d]", &axis);
 
           if(axis>0 && axis<5){
 
               if (indx[axis-1]<3){ dim[indx[axis-1]]=atoi(val);} 
 
               else {
 
                int sum=0; 
 
                token = strtok(val,"{,");
 
                sum+=atoi(token); 
 
                while ( (token=strtok(NULL,"{,")) ) sum+=atoi(token);
 
                printf("total dim=%d\n",sum); 
 
                dim[3]=sum; 
 
               };
 
           }
 
           printf("axis=%d indx=%d\n",axis,indx[axis-1]);          
 
       }
 
       if (strstr(cmd,"name of data file")!= NULL){
 
          sscanf(val,"%s",fimage);
 
          printf("image file=%s\n",fimage);
 
       }
 
       if (strstr(cmd,"!number of bytes per pixel")!= NULL){
 
          dim[4]=atoi(val);
 
          printf("bytes per pixel=%d\n",dim[4]);
 
       } 
 
       if (strstr(cmd,"!number format")!= NULL){
 
          if (strstr(val,"float")!=NULL){ dim[5]=1; };
 
          if (strstr(val,"integer")!=NULL){ dim[5]=0; };
 
          printf("number format=%d\n",dim[5]);
 
       }
 
          
 
     }
 
  }
 
  for (int i=0;i<6;i++) printf ("dim[%d]=%d\n", i,dim[i]);
 
  fclose(fp);
 
  imageread(fimage,dim[0], dim[1], dim[3], dim[5], dim[4], opt);
 
  return 0;
 
}
 
 
 
 
 
int sinoread(char *fname, char *opt=""){
 
  return interfileread(fname,opt);
 
}