/*
************************************************************************
*
*   spec_interface.c -
*
*   Copyright (c) 1995
*
*   ETH Zuerich
*   Institut fuer Molekularbiologie und Biophysik
*   ETH-Hoenggerberg
*   CH-8093 Zuerich
*
*   All Rights Reserved
*
*   Date of last modification : 95/09/15
*   Pathname of SCCS file     : /export/home3/cb/garant-1.0/src/SCCS/s.spec_interface.c
*   SCCS identification       : 1.2
*
************************************************************************
*/
/******************************************************************************
 * interface to access spectral data                                          *
 *                                                                            *
 * Open_Spec_RO .......... Open spectral file for reading and reads parameters*
 * Open_Spec ............. Open spectral file for reading and writing         *
 * Close_Spec ............ Closes spectral file and frees memory              *
 * Read_Spec_Point ....... Reads the intensity at a given coordinate          *
 * Read_Spec_Region ...... Reads the integer intensities in a given region    *
 * Read_Scale_Region ..... Reads integer intensities and scales the size of   *
 *                         the region down                                    *
 * Read_Spec_16bit ....... Reads the 16bit format for a given region          *
 * Read_Spec_Colors ...... Reads selected 8bit points from given region       *
 * Create_Spec ........... Creates a new spectrum file                        *
 * Set_Spec_Calib ........ Sets the Calibration of the spectrum               *
 * Add_Submatrices ....... Adds new submatrices to the spectrum               *
 * Overwrite_Spec_Region . Overwrites the intensities in a given region       *
 * Get_Spec_Param ........ Gets a spectral parameter                          *
 * Calc_8bit_Int ......... Calculate intensity from 8 bit data                *
 * Calc_16bit_Int ........ Calculate intensity from 16 bit data               *
 *                                                                            *
 *                                                                            *
 * programed version 0.1 ............................... 22.8.91 C.Bartels    *
 *                                                                            *
 ******************************************************************************/

/******************************************************************************
 * local routines to access spectral data                                     *
 *                                                                            *
 * Add_Submatrix ...... Adds the information for a  submatrix to the parm file*
 * Comp_Left_Offset ..  Compares left offsets of two subspectra               *
 * Coord_Within_Sub ... Checks if coordinates lie within the subspectrum      *
 * Comp_Overlap_Left_nD Checks if two regions overlap (uses nD data structure *
 * Open_Spec_3D_8 ..... Open spectrum for reading                             *
 * Open_Spec_3D_16 .... Open spectrum for reading                             *
 * Set_Spec_Pos ....... Set position for reading                              *
 * Get_3D_8 ........... Reads byte (char) from 8 bit                          *
 * Get_3D_16 .......... Reads byte from 16 bit and converts into integer      *
 * Reset_Spec_Pos_3D_8 ........... Sets file pointer to position              *
 * Reset_Spec_Pos_3D_16 .......... Sets file pointer to position              *
 * Load_Calibration_Parameters ... Load spectral parameters from disk         *
 * Update_Parameters .. Write spectral parameters onto disk                   *
 * Clip_Copy_Data ..... Clips an submatrix and array of intensities           *
 * Load_Submatrix_Data  Loads a submatrix into an array of intensities        *
 * Load_Data_Vector ... Loads a one dimensional vector of a submatrix         *
 * Set_View_Orientation Permutes parameter of spectrum to match a projection  *
 * Write_Param ........ Creates spectral parameter file with given parameters *
 * Write_Data_Header... Creates spectral data file                            *
 * Write_Sub_Pos ...... Write information of a submatrix into the param. file *
 * logr2 .............. Map integer onto the color scale                      *
 *                                                                            *
 *                                                                            *
 *                                                                            *
 * programed version 0.1 ............................... 22.8.91 C.Bartels    *
 *                                                                            *
 ******************************************************************************/


#include <stdio.h>
#include <math.h>
#include <sys/file.h>
#include <fcntl.h>
#include "constants.h"
#include "data_structures.h"

/* Import from tree_nD.c ***************************************************/
extern struct node_nD *insert_tree();
extern struct node_nD *erase_tree();
extern struct node_nD *search_tree();
extern void apply_tree();
extern void list_tree();

#define READ_INT    1
#define WRITE_INT   2
#define READ_FLOAT  3
#define READ_16     4
#define READ_COLOR  5
#define READ_SCALE  6

#define READ_ONLY   1
#define READ_WRITE  2

float xexp10(x)
float x;
{ return (exp(x*log(10.0))); }
 
int inint(x)
float x;
{ if(x > 0.0) return ( (int)(x+0.5)); else return ( (int)(x-0.5)); }
 

/* Lokale strukturen *******************************************************/

struct rect_nD
{
   int c[DIMENSION];		/* Koordinaten            */
   int s[DIMENSION];		/* Groesse                */
};

struct clip_data_nD
{
   struct rect_nD src;		/* Source cube                       */
   struct rect_nD dst;		/* Destination cube                  */
   int dst_perm[DIMENSION];     /* Permutation of destination        */
   char *flags[DIMENSION];	/* flags whether to read spec. point */
   char *data;			/* Intensity data to read in         */
};

unsigned char *col;
int colsize = 0;

/* Lokale variablen ********************************************************/

static int expon[100];    /* Look up table for 8 bit mode (log to normal)*/
static int reverse[96]; 
static Spec_Param_Data spec_param;
                          /* Spectral parameters to prevent parameter    */
                          /* passing into tree_nD.c subroutines          */
static int rw_mode;       /* Whether to read or write data, used to prevent */
                          /* parameter passing into tree_nD.c subroutines */

/* Lokale routinen *********************************************************/

int Add_Submatrix();
int Clip_Copy_Data();
int Clip_Copy_Colot();
int Comp_Left_Offset();
int Comp_Right_Offset();
int Comp_Overlap_Left_nD();
int Coord_Within_Sub();
int Get_3D_16();
char Get_3D_8();
int Load_Calibration_Parameters();
void Load_Data_Vector();
int Load_Parameters();
void Load_Submatrix_Data();
int Open_Spec_3D_8();
int Open_Spec_3D_16();
void Permute_Subs_Coords();
int Reset_Spec_Pos_3D_16();
int Reset_Spec_Pos_3D_8();
void Set_Spec_Pos();
void Set_View_Orientation();
int Update_Parameters();
void Write_Sub_Pos();
int Write_Param();
int Write_Data_Header();
void Write_Sub_Pos();
int logr2();

/***************************************************************************/
/* global routines                                                         */
/***************************************************************************/

/***************************************************************************/
/* Reads parameters and opens spectrum for raed & write access             */
/* Input: filename                                                         */
/* Output: handle for spectrum                                             */
/* Side effects: opens data file, accesses param file                      */
/***************************************************************************/

Spec_Param_Data *Open_Spec(file)

char file[];               /* File to be loaded                             */

{
   Spec_Param_Data *param; /* Data structure containing spectral parameters */
   char str[MAX_STR];      /* General purpose string            */ 
   int i;

   sprintf(str,"Load spectrum : '%s'",file);

   param = (Spec_Param_Data *)malloc(sizeof(Spec_Param_Data));

   strcpy(param->file,file);
 
   if(Load_Calibration_Parameters(param) == -1) return(NULL);

   for(i=0;i<DIMENSION;i++) {
      param->flags[i] = (char *)malloc(param->size[i] * sizeof(char));
   }

/* Calculate exponentials for 8 bit data */

   for(i = 1; i <= 64; i++)
      expon[i] = (int)(xexp10((float)i/6.6438562)+0.5);

   for(i = 1; i <= 47; i++)
      reverse[i] = (int)(xexp10((float)(i-1)/6.6438562)+0.5);
   for(i = 48; i <= 95; i++)
      reverse[i] = -(int)(xexp10((float)(95 - i)/6.6438562)+0.5);
   reverse[0] = 1;


/* Open file *****************************/

   if (param->file_type == FILE3D_16 &&
       Open_Spec_3D_16(param,READ_WRITE) != FILE3D_16)
      return(NULL);
   else if (param->file_type == FILE3D_8 &&
            Open_Spec_3D_8(param,READ_WRITE) != FILE3D_8)
      return(NULL);
   else return(param);
}

/***************************************************************************/
/* Reads parameters and opens spectrum for reading only                    */
/* Input: filename                                                         */
/* Output: handle for spectrum                                             */
/* Side effects: opens data file, accesses param file                      */
/***************************************************************************/

Spec_Param_Data *Open_Spec_RO(file)

char file[];               /* File to be loaded                             */

{
   Spec_Param_Data *param; /* Data structure containing spectral parameters */
   char str[MAX_STR];      /* General purpose string            */ 
   int i;

   sprintf(str,"Load spectrum : '%s'",file);

   param = (Spec_Param_Data *)malloc(sizeof(Spec_Param_Data));

   strcpy(param->file,file);
 
   if(Load_Calibration_Parameters(param) == -1) return(NULL);

   for(i=0;i<DIMENSION;i++) {
      param->flags[i] = (char *)malloc(param->size[i] * sizeof(char));
   }

/* Calculate exponentials for 8 bit data */

   for(i = 1; i <= 64; i++)
      expon[i] = (int)(xexp10((float)i/6.6438562)+0.5);

   for(i = 1; i <= 47; i++)
      reverse[i] = (int)(xexp10((float)(i - 1)/6.6438562)+0.5);
   for(i = 48; i <= 95; i++)
      reverse[i] = -(int)(xexp10((float)(95 - i)/6.6438562)+0.5);
   reverse[0] = 1;


/* Open file *****************************/

   if (param->file_type == FILE3D_16 &&
       Open_Spec_3D_16(param,READ_ONLY) != FILE3D_16)
      return(NULL);
   else if (param->file_type == FILE3D_8 &&
            Open_Spec_3D_8(param,READ_ONLY) != FILE3D_8)
      return(NULL);
   else return(param);
}

/***************************************************************************/
/* Closes spectral file and frees memory                                   */
/* Input: -                                                                */
/* Output: -                                                               */
/* Side effects: closes data file, frees memory                            */
/***************************************************************************/

void
Close_Spec(param)
Spec_Param_Data *param;   /* Data structure containing spectral parameters */

{
   int i;
   erase_tree(param->subs);
   if(param->fd != -1) {
      close(param->fd);
      /* fclose(param->fp); */
   }
   for(i=0;i<DIMENSION;i++) free(param->flags[i]);
   free(param);
}

/***************************************************************************/
/* Reads spectral intensity at a given position in the spectrum            */
/* Input: handle for spectrum                                              */
/*        position within spectrum                                         */
/* Output: intensity at position (integer value)                           */
/* Side effects: resets file pointer and current_sub                       */
/***************************************************************************/

int
Read_Spec_Point(param,pos)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int pos[DIMENSION];       /* Position within spectrum                      */

{
   spec_param = *param;

   if(param->file_type == FILE3D_16)
   {
      if(Reset_Spec_Pos_3D_16(pos) == 1)
         return((int)Get_3D_16());
      else return(0);
   }
   else
   {
      if(Reset_Spec_Pos_3D_8(pos) == 1)
         return(reverse[Get_3D_8()]);
      else return(0);
   }
}

/***************************************************************************/
/* Reads a spectral region of intensity of arbitrary size, position and    */
/* orientation into a integer array                                        */
/* Input: handle for spectrum                                              */
/*        position within spectrum                                         */
/*        size of region to read                                           */
/*        orientation of sepctrum to read                                  */
/*        array for storing intensities                                    */
/* Output: array of integer values                                         */
/* Side effects: resets file pointer and current_sub                       */
/***************************************************************************/

void
Read_Spec_Region(param,pos,size,perm,array)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int pos[DIMENSION];       /* Position within spectrum                      */
int size[DIMENSION];     /* Groesse der Region	                           */
int perm[DIMENSION];	 /* Permutation Region                             */
int array[];		 /* Zu ladender Vektor                             */

{
   struct clip_data_nD clip;
   int i;
   int data_size;
   int *data_pntr;	/* Pointer into Array */
   int newcolsize;

   if(param->file_type == FILE3D_8) {
      newcolsize = param->sub_size[param->permut[0]];
   } else {
      newcolsize = param->sub_size[param->permut[0]] * 2;
   }
   if(newcolsize > colsize) {
      if(colsize>0) free(col);
      colsize = newcolsize;
      col = (unsigned char *)malloc(sizeof(char) * colsize);
   }

/* Set clipping data ********************************************/
   for(i=param->dim;i<DIMENSION;i++) 
   {
      size[i] = 1;
      pos[i] = 0;
      perm[i] = i;
   }

   for(i=0;i<DIMENSION;i++)
   {
      clip.src.c[i] = pos[i];
      clip.src.s[i] = size[i];
      clip.dst.c[i] = 0;
      clip.dst.s[i] = size[i];
      clip.dst_perm[i] = perm[i];
   }
   clip.data = (char *)&array[0];


/* Clear data for the case that not all data is available ******/
   data_size = 1;
   for(i=0;i<DIMENSION;i++) data_size *= clip.src.s[i];
   for(data_pntr=array;data_pntr<array + data_size; data_pntr++)
		*data_pntr = 0;

/* Find and load matching submatrices **************************/

   spec_param = *param;
   rw_mode = READ_INT;
   apply_tree(spec_param.subs,Clip_Copy_Data,
              Comp_Overlap_Left_nD,&clip,0);


}

/***************************************************************************/
/* Reads a spectral region of 16bit data of arbitrary size, position and   */
/* orientation into a character array                                      */
/* Input: handle for spectrum                                              */
/*        position within spectrum                                         */
/*        size of region to read                                           */
/*        orientation of sepctrum to read                                  */
/*        array for storing intensities                                    */
/* Output: array of integer values                                         */
/* Side effects: resets file pointer and current_sub                       */
/***************************************************************************/

void
Read_Spec_16bit(param,pos,size,perm,array)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int pos[DIMENSION];       /* Position within spectrum                      */
int size[DIMENSION];     /* Groesse der Region	                           */
int perm[DIMENSION];	 /* Permutation Region                             */
char array[];		 /* Zu ladender Vektor                             */

{
   struct clip_data_nD clip;
   int i;
   int data_size;
   char *data_pntr;	/* Pointer into Array */
   int newcolsize;

   if(param->file_type == FILE3D_8) {
      newcolsize = param->sub_size[param->permut[0]];
   } else {
      newcolsize = param->sub_size[param->permut[0]] * 2;
   }
   if(newcolsize > colsize) {
      if(colsize>0) free(col);
      colsize = newcolsize;
      col = (unsigned char *)malloc(sizeof(char) * colsize);
   }

/* Set clipping data ********************************************/
   for(i=param->dim;i<DIMENSION;i++) 
   {
      size[i] = 1;
      pos[i] = 0;
      perm[i] = i;
   }

   for(i=0;i<DIMENSION;i++)
   {
      clip.src.c[i] = pos[i];
      clip.src.s[i] = size[i];
      clip.dst.c[i] = 0;
      clip.dst.s[i] = size[i];
      clip.dst_perm[i] = perm[i];
   }
   clip.data = &array[0];


/* Clear data for the case that not all data is available ******/
   data_size = 2;
   for(i=0;i<DIMENSION;i++) data_size *= clip.src.s[i];
   for(data_pntr=array;data_pntr<array + data_size; data_pntr++)
		*data_pntr = 0;

/* Find and load matching submatrices **************************/

   spec_param = *param;
   rw_mode = READ_16;
   apply_tree(spec_param.subs,Clip_Copy_Data,
              Comp_Overlap_Left_nD,&clip,0);


}

/***************************************************************************/
/* Reads points specified by  the flags from spectral region of arbitrary  */
/* size, position and orientation into a character array                   */
/* Input: handle for spectrum                                              */
/*        position within spectrum                                         */
/*        size of region to read                                           */
/*        orientation of sepctrum to read                                  */
/*	  sizes of the destination array (must be smaller or equal than the*/
/*	     size of the region to read                                    */
/*        array for storing intensities                                    */
/* Output: array of integer values                                         */
/* Side effects: resets file pointer and current_sub                       */
/***************************************************************************/

void
Read_Spec_Colors(param,pos,size,perm,dstsize,array)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int pos[DIMENSION];       /* Position within spectrum                      */
int size[DIMENSION];     /* Groesse der Region	                           */
int perm[DIMENSION];	 /* Permutation Region                             */
int dstsize[DIMENSION];  /* destination size                               */
char array[];		 /* Zu ladender Vektor                             */

{
   struct clip_data_nD clip;
   int i,j;
   int data_size;
   char *data_pntr;	/* Pointer into Array */
   float factor,offset;
   int newcolsize;
   int von,bis,set;

   if(param->file_type == FILE3D_8) {
      newcolsize = param->sub_size[param->permut[0]];
   } else {
      newcolsize = param->sub_size[param->permut[0]] * 2;
   }
   if(newcolsize > colsize) {
      if(colsize>0) free(col);
      colsize = newcolsize;
      col = (unsigned char *)malloc(sizeof(char) * colsize);
   }

/* Set clipping data ********************************************/

   for(i=param->dim;i<DIMENSION;i++) 
   {
      size[i] = 1;
      dstsize[i] = 1;
      pos[i] = 0;
      perm[i] = i;
   }

   for(i=0;i<DIMENSION;i++)
   {
      clip.src.c[i] = pos[i];
      clip.src.s[i] = size[i];
      clip.dst.c[i] = 0;
      clip.dst.s[i] = dstsize[i];
      clip.dst_perm[i] = perm[i];
      clip.flags[i] = param->flags[i];
      von = pos[i]; if(von<0) von=0;
      bis = pos[i]+size[i]; if(bis>param->size[i]) bis = param->size[i];
      for(j=von;j<bis;j++) clip.flags[i][j] = 0;
      factor = (float)size[i] / (float)dstsize[i];
      offset = pos[i];
      for(j=0;j<dstsize[i];j++) {
         set = (int)((j + 0.4999) * factor + offset);
         if(set>=0 && set<bis) 
            clip.flags[i][set] = 1;
      }
   }
   clip.data = &array[0];


/* Clear data for the case that not all data is available ******/
   data_size = 1;
   for(i=0;i<DIMENSION;i++) data_size *= clip.dst.s[i];
   for(data_pntr=array;data_pntr<array + data_size; data_pntr++)
		*data_pntr = 0;

/* Find and load matching submatrices **************************/

   spec_param = *param;
   rw_mode = READ_COLOR;
   apply_tree(spec_param.subs,Clip_Copy_Data,
              Comp_Overlap_Left_nD,&clip,0);


}

/***************************************************************************/
/* Reads points specified by  destination size from spectral region of     */
/* arbitrary size, position and orientation into an integer array          */
/* Input: handle for spectrum                                              */
/*        position within spectrum                                         */
/*        size of region to read                                           */
/*        orientation of sepctrum to read                                  */
/*	  sizes of the destination array (must be smaller or equal than the*/
/*	     size of the region to read                                    */
/*        array for storing intensities                                    */
/* Output: array of integer values                                         */
/* Side effects: resets file pointer and current_sub                       */
/***************************************************************************/

void
Read_Scale_Region(param,pos,size,perm,dstsize,array)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int pos[DIMENSION];       /* Position within spectrum                      */
int size[DIMENSION];     /* Groesse der Region	                           */
int perm[DIMENSION];	 /* Permutation Region                             */
int dstsize[DIMENSION];  /* destination size                               */
int array[];		 /* Zu ladender Vektor                             */

{
   struct clip_data_nD clip;
   int i,j;
   int data_size;
   int *data_pntr;	/* Pointer into Array */
   float factor,offset;
   int newcolsize;
   int von,bis,set;

   if(param->file_type == FILE3D_8) {
      newcolsize = param->sub_size[param->permut[0]];
   } else {
      newcolsize = param->sub_size[param->permut[0]] * 2;
   }
   if(newcolsize > colsize) {
      if(colsize>0) free(col);
      colsize = newcolsize;
      col = (unsigned char *)malloc(sizeof(char) * colsize);
   }

/* Set clipping data ********************************************/

   for(i=param->dim;i<DIMENSION;i++) 
   {
      size[i] = 1;
      dstsize[i] = 1;
      pos[i] = 0;
      perm[i] = i;
   }

   for(i=0;i<DIMENSION;i++)
   {
      clip.src.c[i] = pos[i];
      clip.src.s[i] = size[i];
      clip.dst.c[i] = 0;
      clip.dst.s[i] = dstsize[i];
      clip.dst_perm[i] = perm[i];
      clip.flags[i] = param->flags[i];
      von = pos[i]; if(von<0) von=0;
      bis = pos[i]+size[i]; if(bis>param->size[i]) bis = param->size[i];
      for(j=von;j<bis;j++) clip.flags[i][j] = 0;
      factor = (float)size[i] / (float)dstsize[i];
      offset = pos[i];
      for(j=0;j<dstsize[i];j++) {
         set = (int)((j + 0.4999) * factor + offset);
         if(set>=0 && set<bis) 
            clip.flags[i][set] = 1;
      }
   }
   clip.data = (char *)&array[0];

/* Clear data for the case that not all data is available ******/
   data_size = 1;
   for(i=0;i<DIMENSION;i++) data_size *= clip.dst.s[i];
   for(data_pntr=array;data_pntr<array + data_size; data_pntr++)
		*data_pntr = 0;

/* Find and load matching submatrices **************************/

   spec_param = *param;
   rw_mode = READ_SCALE;
   apply_tree(spec_param.subs,Clip_Copy_Data,
              Comp_Overlap_Left_nD,&clip,0);


}

/***************************************************************************/
/* Overwrites an allready existing  spectral region of intensity of        */
/* arbitrary size, position and orientation                                */
/* Input: handle for spectrum                                              */
/*        position within spectrum                                         */
/*        size of region to read                                           */
/*        orientation of the array of intensities                          */
/*        array with intensities                                           */
/* Output: -                                                               */
/* Side effects: changes data file, resets file pointer and current_sub    */
/***************************************************************************/

void
Overwrite_Spec_Region(param,pos,size,perm,array)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int pos[DIMENSION];       /* Position within spectrum                      */
int size[DIMENSION];     /* Groesse der Region	                           */
int perm[DIMENSION];	 /* Permutation Region                             */
int array[];		 /* Zu schreibende Daten                           */

{
   struct clip_data_nD clip;
   int i;
   int newcolsize;

   if(param->file_type == FILE3D_8) {
      newcolsize = param->sub_size[param->permut[0]];
   } else {
      newcolsize = param->sub_size[param->permut[0]] * 2;
   }
   if(newcolsize > colsize) {
      if(colsize>0) free(col);
      colsize = newcolsize;
      col = (unsigned char *)malloc(sizeof(char) * colsize);
   }

   for(i=param->dim;i<DIMENSION;i++) 
   {
      size[i] = 1;
      pos[i] = 0;
      perm[i] = i;
   }

/* Set clipping data ********************************************/
   for(i=0;i<DIMENSION;i++)
   {
      clip.src.c[i] = pos[i];
      clip.src.s[i] = size[i];
      clip.dst.c[i] = 0;
      clip.dst.s[i] = size[i];
      clip.dst_perm[i] = perm[i];
   }
   clip.data = (char *)&array[0];

/* Find and load matching submatrices **************************/

   spec_param = *param;
   rw_mode = WRITE_INT;
   apply_tree(spec_param.subs,Clip_Copy_Data,
              Comp_Overlap_Left_nD,&clip,0);
}


/***************************************************************************/
/* Creates a spectrum                                                      */
/* Input: filename                                                         */
/*        number of dimensions                                             */
/*        size of spectrum                                                 */
/*        size of the submatrices                                          */
/*        permutation of the submatrices                                   */
/*        filetype (i.e. 8 or 16 bit file)                                 */
/* Output: returns handle for spectrum                                     */
/* Side effects: creates data and parameterfile                            */
/***************************************************************************/

int
Create_Spec(file,type,dim,size,sub_size,sub_perm)

char file[];
int type;
int  dim;
int  size[DIMENSION],sub_size[DIMENSION],sub_perm[DIMENSION];

{
   Spec_Param_Data param;
   int i;
   char fn_param[MAX_STR];
   char fn[MAX_STR];
   int status;

/* Define and open output files *********************************************/
   strcpy(fn_param,file);
   strcat(fn_param,".3D.param");
   param.fp_param = fopen(fn_param,"w");
   if (param.fp_param == 0)
   {
      printf("\nProblems opening parameter file");
      return(0);
   }

   if(type == FILE3D_8)
   {
      strcpy(fn,file);
      strcat(fn,".3D.8");
      param.fp = fopen(fn,"w");
      if(param.fp == 0)
      {
         printf("\nProblems opening '8' bit file");
         return(0);
      }
   }

   if(type == FILE3D_16)
   {
      strcpy(fn,file);
      strcat(fn,".3D.16");
      param.fp = fopen(fn,"w");

      if(param.fp == 0)
      {
         printf("Problems opening 16 bit file\n");
         return(0);
      }
   }

   strcpy(param.file,file);
   param.file_type = type;
   param.dim = dim;
   param.subs = NULL;

   for(i=0;i<dim;i++) param.size[i] = size[i];
   for(i=0;i<dim;i++) param.sub_size[i] = sub_size[i];
   for(i=0;i<dim;i++) param.sub_perm[i] = sub_perm[i];

   for(i=0;i<dim;i++) param.sw[i] = 0.0;
   for(i=0;i<dim;i++) param.max_ppm[i] = 0.0;
   for(i=0;i<dim;i++) param.w0[i] = 0.0;
   for(i=0;i<dim;i++) param.folding[i] = RSH;

   /* write version 1 parameter file */
   status = (Write_Param(&param,1) && Write_Data_Header(&param,1));

   fclose(param.fp_param);
   fclose(param.fp);

   return(status);
}


/***************************************************************************/
/* Sets the calibration of a spectrum                                      */
/* Input: spectrum handle                                                  */
/*        maximal ppm in each dimension                                    */
/*        sweep width in each dimension                                    */
/* Output: flag: 1=success, 0=failure                                      */
/* Side effects: modifies data and parameterfile                            */
/***************************************************************************/

int
Set_Spec_Calib(file,sf,max_ppm,sw,type)

char file[];
float sf[DIMENSION],max_ppm[DIMENSION],sw[DIMENSION];
char type;

{
   char filename_param[MAX_STR];  /* Full filename for spectrum */
   Spec_Param_Data param;
   int i, status;
   int version;

   strcpy(param.file,file);

   version = Load_Calibration_Parameters(&param);
   if(version == -1) return(0);

   for(i=0;i<param.dim;i++) param.sw[i] = sw[i];
   for(i=0;i<param.dim;i++) param.max_ppm[i] = max_ppm[i];
   for(i=0;i<param.dim;i++) param.w0[i] = sf[i];
   param.type = type;   

   sprintf(filename_param,"%s.3D.param",param.file);
   if((param.fp_param = fopen(filename_param,"w")) == NULL)
   {
      fclose(param.fp_param);
      return(-1);
   }

   status = Write_Param(&param,version);

   fclose(param.fp_param);

   if(status) return(1); else return(0);
}


/***************************************************************************/
/* Gets spectral parameter                                                 */
/* Input: handle of spectrum                                               */
/*        type of parameter                                                */
/*        dimension of parameter (optional)                                */
/* Output: integer value of desired parameter                              */
/***************************************************************************/

void
Get_Spec_Param(param,what,value)

Spec_Param_Data *param;
char *value;
int what;

{
   int i;
   int *int_ptr;
   float *flt_ptr;

   int_ptr = (int *)value;
   flt_ptr = (float *)value;

   if(what == SP_DIM) *int_ptr = param->dim;
   if(what == SP_SIZE) for(i=0;i<param->dim;i++) int_ptr[i] = param->size[i];
   if(what == SUB_SIZE)
      for(i=0;i<param->dim;i++) int_ptr[i] = param->sub_size[i];
   if(what == SUB_PERM)
      for(i=0;i<param->dim;i++) int_ptr[i] = param->sub_perm[i];
   if(what == SP_SW) for(i=0;i<param->dim;i++) flt_ptr[i] = param->sw[i];
   if(what == SP_MAXPPM) for(i=0;i<param->dim;i++)
                               flt_ptr[i] = param->max_ppm[i];
   if(what == SP_SFO1) for(i=0;i<param->dim;i++) flt_ptr[i] = param->w0[i];
   if(what == SP_TYPE) *value = param->type;
   if(what == SP_FTYPE) *int_ptr = param->file_type;
   if(what == SP_FOLDT) for(i=0;i<param->dim;i++)
                            int_ptr[i] = param->folding[i];

}


/*****************************************************************************/
/* Adds submatrices of a given total size at a given position to the         */
/* spectrum. The calling routine  is responsible, that there is              */
/* no overlap between different submatrices                                  */
/* Input: handle for spectrum                                                */
/*        offset of submatrices within spectrum                              */
/*        total size of region to add submatrices                            */
/* Output: -                                                                 */
/* Side effects: An empty submatrix is appended to the .param file and the   */
/*               size is adjusted                                            */
/*****************************************************************************/

int
Add_Submatrices(file,sub_off,size)

char file[];
int sub_off[DIMENSION];
int size[DIMENSION];

{
   Spec_Param_Data param;
   char filename_param[MAX_STR];  /* Full filename for spectrum */
   int i;
   int status,incdim,inc;
   int cur_off[DIMENSION];
   int version;
   long tot_size;

   strcpy(param.file,file);
   version = Load_Calibration_Parameters(&param);
   if(version == -1) return(0);

   if(param.file_type == FILE3D_8)
   {
      sprintf(filename_param,"%s.3D.8",param.file);
      if((param.fp = fopen(filename_param,"a")) == NULL)
      {
         fclose(param.fp);
         return(0);
      }
   }

   if(param.file_type == FILE3D_16)
   {
      sprintf(filename_param,"%s.3D.16",param.file);
      if((param.fp = fopen(filename_param,"a")) == NULL)
      {
         fclose(param.fp);
         return(0);
      }
   }

   status = 0;

   if(version == 0) {
      sprintf(filename_param,"%s.3D.param",param.file);
      if((param.fp_param = fopen(filename_param,"a")) == NULL)
      {
         fclose(param.fp_param);
         return(0);
      }

      for(i=0;i<param.dim;i++) cur_off[i]=sub_off[i];

      do
      {
         status |= Add_Submatrix(&param,cur_off);

         /* Increment position of submatrix */
         incdim=0;
         inc=FALSE;
         do
         {
            cur_off[incdim] += param.sub_size[incdim];
            if(cur_off[incdim] >= sub_off[incdim] + size[incdim])
            {
               cur_off[incdim] = sub_off[incdim];
               incdim++;
            }
            else inc = TRUE;
         }
         while(!inc && incdim < param.dim);
      }
      while(inc);

      fclose(param.fp_param);

   } else if (version == 1) {
      if(param.file_type == FILE3D_8) tot_size = 1;
      else tot_size = 2;
      for(i=0;i<param.dim;i++) tot_size *= param.size[i];

      status |= fseek(param.fp,tot_size,1);
   }

   /* Write a character to the end of the file to reserve space */
   status |= fseek(param.fp,-1,1);
   status |= (putc(0,param.fp) != 0);

   fclose(param.fp);

   if(status) return(0);
   else return(1);
}


/***************************************************************************/
/* local routines                                                          */
/***************************************************************************/

/*****************************************************************************/
/* Adds a submatrix at a given position to the spectrum                      */
/* Input: handle for spectrum                                                */
/*        offset of the submatrix within spectrum                            */
/* Output: -                                                                 */
/* Side effects: An empty submatrix is appended to the .param file and the   */
/*****************************************************************************/

int
Add_Submatrix(param,sub_off)

Spec_Param_Data *param;
int sub_off[DIMENSION];

{
   long tot_size;
   int i;

   for(i=0;i<param->dim;i++)
   {
      putw(sub_off[i],param->fp_param);
   }
   putw((int)ftell(param->fp),param->fp_param);

   if(param->file_type == FILE3D_8) tot_size = 1;
   else tot_size = 2;
   for(i=0;i<param->dim;i++) tot_size *= param->sub_size[i];

   if(fseek(param->fp,tot_size,1) != 0) return(1); 
   else return(0); 
}


/*****************************************************************************/
/* Compares left offsets of two subspectra                                   */
/*****************************************************************************/

int
Comp_Left_Offset(a,b,dim)

Subspectrum *a,*b;
int dim;

{
   return( (*b).offset[dim] - (*a).offset[dim]);
}



/*****************************************************************************
 * Checks if coordinates lie within the subspectrum                          *
 * result is offset difference of left offsets                               *
 *****************************************************************************/

int
Coord_Within_Sub(sub,coordinate,dim)

Subspectrum *sub;
int *coordinate;
int dim;

{
   int dif;

   dif = coordinate[dim] - sub->offset[dim];
   if ((dif >= 0) && (dif < spec_param.sub_size[dim])) dif = 0;
   return(dif);
}


/*****************************************************************************
 * Checks if left offset of subspectrum lies within the clipping region      *
 * result is difference of left offsets                                      *
 *****************************************************************************/

int
Comp_Overlap_Left_nD(sub,clip,dim)

Subspectrum *sub;
struct clip_data_nD *clip;
int dim;

{
   int dif;

   dif = clip->src.c[dim] - sub->offset[dim];
   if((dif > (-clip->src.s[dim])) && (dif < spec_param.sub_size[dim])) 
		dif = 0;

   return(dif);
}
      

/****************************************************************************/
/* Open spectrum for reading and overwriting                                */
/****************************************************************************/

int
Open_Spec_3D_8(param,mode)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int mode;		  /* open for reading or reading & writing ??      */

{
   int i,j;
   long multiplier;
   char filename[MAX_STR]; /* Full name for data file */

   sprintf(filename,"%s.3D.8",param->file);
   if(mode == READ_WRITE) {
      if(!(param->fd = open(filename,O_RDWR,0)))
      {
         return(-1);
      } 
   } else {
      if(!(param->fd = open(filename,O_RDONLY,0)))
      {
         return(-1);
      } 
   }

   multiplier = 1;
   for(i=0;i<DIMENSION;i++) {
      for(j=0;j<DIMENSION;j++) if(param->sub_perm[j] == i) {
         param->increment[j] = multiplier;
         multiplier *= param->sub_size[j];
         param->permut[i] = j;
      }
   }

   param->current_sub = NULL;
   return(FILE3D_8);
}

/****************************************************************************/
/* Open spectrum for reading and overwriting                                */
/****************************************************************************/

int
Open_Spec_3D_16(param,mode)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */
int mode;		  /* open for reading or reading & writing ??      */

{
   int i,j;
   long multiplier;
   char filename[MAX_STR]; /* Full name for data file */

   sprintf(filename,"%s.3D.16",param->file);
   if(mode == READ_WRITE) {
      if(!(param->fd = open(filename,O_RDWR,0)))
      {
         return(-1);
      } 
   } else {
      if(!(param->fd = open(filename,O_RDONLY,0)))
      {
         return(-1);
      } 
   }

   multiplier = 2;
   for(i=0;i<DIMENSION;i++) {
      for(j=0;j<DIMENSION;j++) if(param->sub_perm[j] == i) {
         param->increment[j] = multiplier;
         multiplier *= param->sub_size[j];
         param->permut[i] = j;
      }
   }

   param->current_sub = NULL;
   return(FILE3D_16);
}

/***************************************************************************/
/* Read byte from 8 bit file 3D Version and convert into integer           */
/***************************************************************************/

char
Get_3D_8()

{
   char ch;

   read(spec_param.fd,&ch,1);
   return ch;
}

/***************************************************************************/
/* Read byte from 16 bit file 3D Version and convert into integer          */
/***************************************************************************/

int
Get_3D_16()


{
   unsigned char a;     /* Mantisse of 16 bit format                        */
   unsigned char n;     /* Exponent of 16 bit format (base sqrt(2))         */ 

   read(spec_param.fd,&a,1);
   read(spec_param.fd,&n,1);

   return((int)((a+615)*reverse[n]/721));
}


/***************************************************************************/
/* sets file pointer to position:                                          */
/*       find adaquate subspectrum                                         */
/*       calculate offset for filepointer                                  */
/*       set filepointer                                                   */
/***************************************************************************/
 
int
Reset_Spec_Pos_3D_8(pos)

int pos[DIMENSION];       /* Position within spectrum                      */

{
   long file_offset;
   int i,res;
   struct node_nD *n;

/******* find adaquate subspectrum                                         */
   if (n = search_tree(spec_param.subs,Coord_Within_Sub,pos,0))
   {
/******* calculate offset for filepointer                                  */
      spec_param.current_sub = (Subspectrum *)n->element;
      file_offset = spec_param.current_sub->pos;
      for (i=0; i<DIMENSION; i++)
      {
         file_offset += spec_param.increment[i]
                        * (pos[i] - spec_param.current_sub->offset[i]);
      }

/****** set absolute filepointer                                          */
      res = (lseek(spec_param.fd, file_offset, 0) != -1);

/****** check status and return                                           */
      if (!res) return(0);
      else return(1);
   }

   else return(0);
}



/***************************************************************************/
/* sets file pointer to position:                                          */
/*       find adaquate subspectrum                                         */
/*       calculate offset for filepointer                                  */
/*       set filepointer                                                   */
/***************************************************************************/
 
int
Reset_Spec_Pos_3D_16(pos)

int pos[DIMENSION];       /* Position within spectrum                      */

{
   long file_offset;
   int i,res;
   struct node_nD *n;

/******* find adaquate subspectrum                                         */
   if (n = search_tree(spec_param.subs,Coord_Within_Sub,pos,0))
   {
/******* calculate offset for filepointer                                  */
      spec_param.current_sub = (Subspectrum *)n->element;
      file_offset = spec_param.current_sub->pos;
      for (i=0; i<DIMENSION; i++)
      {
         file_offset += spec_param.increment[i]
                        * (pos[i] - spec_param.current_sub->offset[i]);
      }

/****** set absolute filepointer                                          */
      res = (lseek(spec_param.fd, file_offset, 0) != -1);

/****** check status and return                                           */
      if (!res) return(0);
      else return(1);
   }
   else return(0);
}


/*****************************************************************************
*                     Load spectral parameters from disk                     *
*****************************************************************************/

int
Load_Calibration_Parameters(param)

Spec_Param_Data *param;   /* Data structure containing spectral parameters */

{
   Spec_Param_Data temp_spec;
   FILE *fp;                      /* parameter file             */
   char filename_param[MAX_STR];  /* Full filename for spectrum */
   int status;			  /* Whether file is ok         */
   int temp;
   char string[5];
   int i,dummy;
   Subspectrum *subspec;

   int pos[DIMENSION];		  /* position in spectrum       */
   int incdim, inc, incperm;
   int nr_sub, tot_size;
   int version;


/* Open parameter file ****************************/

   sprintf(filename_param,"%s.3D.param",param->file);
   if((fp = fopen(filename_param,"r")) == NULL)
   {
      fclose(fp);
      return(-1);
   }

   temp_spec = *param;

/* Load parameter file and store parameters *******/

   status = (fscanf(fp,"Version ....................... %d\n",&version) == 1);
   if(version == 0) {
      status &= (fscanf(fp,"Number of dimensionen ......... %d\n",
		&temp_spec.dim) == 1);
      status &= (fscanf(fp,"16 or 8 bit file type ......... %s\n",string) == 1);
      if (string[0] == '1') temp_spec.file_type = FILE3D_16;
      else temp_spec.file_type = FILE3D_8;
      for (i=0;i<temp_spec.dim;i++)
         status &= (fscanf(fp,"Spectrometer frequency in w%d .. %f\n",&dummy,
	   &temp_spec.w0[i]) == 2);
      for (i=0;i<temp_spec.dim;i++)
        status &= (fscanf(fp,"Spectral sweep width in w%d .... %f\n",&dummy,
               &temp_spec.sw[i]) == 2);
      for (i=0;i<temp_spec.dim;i++)
        status &= (fscanf(fp,"Maximum chemical shift in w%d .. %f\n",&dummy,
                &temp_spec.max_ppm[i]) == 2);
      for (i=0;i<temp_spec.dim;i++)
        status &= (fscanf(fp,"Size of spectrum in w%d ........ %d\n",&dummy,
                &temp_spec.size[i]) == 2);
      for ( ;i<DIMENSION; i++) temp_spec.size[i] = 1;
      status &= (fscanf(fp,"Type of spectrum .............. %c\n",
                &temp_spec.type) == 1);

/* Default setting of folding to RSH ****************************************/

      for(i=0;i<DIMENSION;i++) temp_spec.folding[i] = RSH;

/* Read subspectra information ******************************************/

      status &= (fscanf(fp,"\n\n\nInformation for subspectra: \n") == 0);

      for (i=0;i<temp_spec.dim;i++)
      {
         temp_spec.sub_size[i] = getw(fp);
      }
      for ( ;i<DIMENSION;i++)
         temp_spec.sub_size[i] = 1;

      for (i=0;i<temp_spec.dim;i++)
         temp_spec.sub_perm[i] = getw(fp);
      for ( ;i<DIMENSION;i++)
         temp_spec.sub_perm[i] = i;

      temp_spec.subs = NULL;		/* init tree */


      temp = getw(fp);	/* read to look if end of file has been reached */

      while (!feof(fp))
      {
         subspec = (Subspectrum *)malloc(sizeof(Subspectrum));
						/* allocate memory */
         (*subspec).offset[0] = temp;
         for (i=1; i<temp_spec.dim; i++)
            (*subspec).offset[i] = getw(fp);
         for ( ;i<DIMENSION;i++)
            (*subspec).offset[i] = 0;


         (*subspec).pos = getw(fp);   /* position of the file pointer */

         if (ferror(fp))
         {
            fclose(fp);
            free(subspec);
            temp_spec.subs = erase_tree(temp_spec.subs);
            status &= FALSE;
         }

         dummy = FALSE;
         temp_spec.subs = insert_tree(subspec,temp_spec.subs,&dummy,
                                           Comp_Left_Offset,0);

         temp = getw(fp);   /* read to see if end of file has been reached */
      }
   } else if(version == 1) {
      status &= (fscanf(fp,"Number of dimensions ..........%d\n",
		&temp_spec.dim) == 1);
      status &= (fscanf(fp,"16 or 8 bit file type .........%s\n",string) == 1);
      if (string[0] == '1') temp_spec.file_type = FILE3D_16;
      else temp_spec.file_type = FILE3D_8;
      for (i=0;i<temp_spec.dim;i++)
         status &= (fscanf(fp,"Spectrometer frequency in w%d ..%f\n",&dummy,
	   &temp_spec.w0[i]) == 2);
      for (i=0;i<temp_spec.dim;i++)
         status &= (fscanf(fp,"Spectral sweep width in w%d ....%f\n",&dummy,
               &temp_spec.sw[i]) == 2);
      for (i=0;i<temp_spec.dim;i++)
         status &= (fscanf(fp,"Maximum chemical shift in w%d ..%f\n",&dummy,
                &temp_spec.max_ppm[i]) == 2);
      for (i=0;i<temp_spec.dim;i++)
         status &= (fscanf(fp,"Size of spectrum in w%d ........%d\n",&dummy,
                &temp_spec.size[i]) == 2);
      for ( ;i<DIMENSION; i++) temp_spec.size[i] = 1;
      for (i=0;i<temp_spec.dim;i++)
         status &= (fscanf(fp,"Submatrix size in w%d ..........%d\n",&dummy,
                &temp_spec.sub_size[i]) == 2);
      for ( ;i<DIMENSION; i++) temp_spec.sub_size[i] = 1;
      for(i=0;i<DIMENSION;i++) temp_spec.sub_perm[i] = i;
      for (i=0;i<temp_spec.dim;i++) {
         status &= (fscanf(fp,"Permutation for w%d ............%d\n",&dummy,
                &temp_spec.sub_perm[i]) == 2);
         temp_spec.sub_perm[i] -= 1;
      }
      for (i=0;i<temp_spec.dim;i++) {
         status &= (fscanf(fp,"Folding in w%d .................%10s\n",&dummy,
                           string) == 2);
         if(strcmp(string,"TPPI") == 0) temp_spec.folding[i] = TPPI;
         else temp_spec.folding[i] = RSH;
      }
      for ( ;i<DIMENSION; i++) temp_spec.folding[i] = RSH;
      status &= (fscanf(fp,"Type of spectrum ..............%c\n",
                &temp_spec.type) == 1);

/* Init subspectra information **********************************************/

      if(temp_spec.file_type == FILE3D_8) tot_size = 1;
      else tot_size = 2;
      for(i=0;i<DIMENSION;i++) tot_size *= temp_spec.sub_size[i];

      nr_sub = 0;

      temp_spec.subs = NULL;		/* init tree */

/* loop over spectrum */
      for(i=0;i<DIMENSION;i++) pos[i] = 0;
      do
      {
         subspec = (Subspectrum *)malloc(sizeof(Subspectrum));
         for (i=0; i<temp_spec.dim; i++)
            (*subspec).offset[i] = pos[i];
         for ( ;i<DIMENSION;i++)
            (*subspec).offset[i] = 0;

         (*subspec).pos = nr_sub * tot_size;

         dummy = FALSE;
         temp_spec.subs = insert_tree(subspec,temp_spec.subs,&dummy,
                                           Comp_Left_Offset,0);

         /* Increment position of submatrix */
         nr_sub++;
         incperm = 0;
         for(i=0;i<temp_spec.dim;i++) if(temp_spec.sub_perm[i] == incperm)
            incdim = i;
         inc=FALSE;
         do
         {
            pos[incdim] += temp_spec.sub_size[incdim];
            if(pos[incdim] >= temp_spec.size[incdim])
            {
               pos[incdim] = 0;
               incperm++;
               for(i=0;i<temp_spec.dim;i++)
                  if(temp_spec.sub_perm[i] == incperm)
                     incdim = i;
            }
            else inc = TRUE;
         }
         while(!inc && incperm < temp_spec.dim);
      }
      while(inc);

   } else {
      status = 0;
   }

   fclose(fp);

/* Set permutation for global view *****************************************/

   for(i=0;i<DIMENSION;i++) temp_spec.view_perm[i] = i;

   if (status)
   {
      *param = temp_spec;
       return(version);
   }

/* Bad parameter file *************************/

   else return(-1);
}

int
Clip_Copy_Data(submatrix,clip)

struct clip_data_nD *clip;
Subspectrum *submatrix;

{
   struct clip_data_nD new_clip;
   int i;
   int dc,ds;

   new_clip = *clip;

/* Set submatrix **********************************************************/

   spec_param.current_sub = submatrix;

/* Clip spectral coordinates **********************************************/

   for(i=0;i<spec_param.dim;i++)
   {
      dc = clip->src.c[i] - spec_param.current_sub->offset[i];
      if (dc<0) new_clip.src.c[i] = spec_param.current_sub->offset[i];

      ds = dc + clip->src.s[i] - spec_param.sub_size[i];
      if (ds<0) new_clip.src.s[i] = clip->src.c[i] + clip->src.s[i]
				    - new_clip.src.c[i];
      else new_clip.src.s[i] = spec_param.current_sub->offset[i]
                   + spec_param.sub_size[i] - new_clip.src.c[i];
   }

/* Calculate offsets and size in destination array *************************/
/* the size contains the total size of the destination array ***************/

   for(i=0;i<spec_param.dim;i++)
   {
      if(clip->src.s[i] != clip->dst.s[i])
         new_clip.dst.c[i] = 
            inint(((float)new_clip.src.c[i] - (float)clip->src.c[i]) /
                  (float)clip->src.s[i] * (float)clip->dst.s[i]);
      else 
         new_clip.dst.c[i] = new_clip.src.c[i] - clip->src.c[i];
      new_clip.dst.s[i] = clip->dst.s[i];	/* Total size of dst array ! */
   }

/* Load Submatrix *********************************************************/

   Load_Submatrix_Data(&new_clip);
}

void
Load_Submatrix_Data(clip)

struct clip_data_nD *clip;

{
   char *data_pos;		/* Pointer in the destination data           */
   int file_pos;		/* Position within file                      */
   int data_inc[DIMENSION];     /* Increment to advance one point in the data*/
   int file_inc[DIMENSION];     /* Increment to advance one point in the file*/
   int index;
   int i;

/* Set data and file pointer to start of array resp. submatrix ***************/
   file_pos = spec_param.current_sub->pos;
   data_pos = clip->data;

/* Calculate increments for file and destination array ***********************/

   if(rw_mode == READ_INT) data_inc[clip->dst_perm[0]] = sizeof(int);
   else if(rw_mode == READ_16) data_inc[clip->dst_perm[0]] = 2;
   else if(rw_mode == WRITE_INT) data_inc[clip->dst_perm[0]] = sizeof(int);
   else if(rw_mode == READ_COLOR) data_inc[clip->dst_perm[0]] = 1;
   else if(rw_mode == READ_SCALE) data_inc[clip->dst_perm[0]] = sizeof(int);

   for (i=0; i<spec_param.dim; i++)
   {
      file_inc[i] = spec_param.increment[i];
      file_pos += file_inc[i]
               * (clip->src.c[i] - spec_param.current_sub->offset[i]);

      index = clip->dst_perm[i];
      data_pos += data_inc[index] * clip->dst.c[index];
      if (i<(spec_param.dim - 1))
         data_inc[clip->dst_perm[i+1]] = clip->dst.s[index] * data_inc[index];

   }


/* Recursively load one dimensional vectors using the increments ************/

   Load_Data_Vector(clip,spec_param.dim-1,data_pos,file_pos,data_inc,file_inc);
}


/****************************************************************************/
/* Recursively load one dimensional vectors using the increments            */
/****************************************************************************/

void
Load_Data_Vector(clip,dim,data_pos,file_pos,data_inc,file_inc)

struct clip_data_nD *clip;   /* Clipping information			  */
int dim;	 	     /* No. of dimension to load                  */
char *data_pos;		     /* Pointer in the destination data           */
int file_pos;		     /* Position within file                      */
int data_inc[DIMENSION];     /* Increment to advance one point in the data*/
int file_inc[DIMENSION];     /* Increment to advance one point in the file*/

{
   int i;
   int index;
   unsigned char a;   /* Mantisse of 16 bit format                        */
   unsigned char n;   /* Exponent of 16 bit format (base sqrt(2))         */ 

   index = spec_param.permut[dim];

/* Dimension 0 --> read data (file access is sequential) ******************/
   if(dim == 0)
   {
/****** set absolute filepointer                                          */
      lseek(spec_param.fd, file_pos, 0);

/****** read integer data from file ***************************************/
      if(rw_mode == READ_INT)
      {
         if(spec_param.file_type == FILE3D_8)
         {
            read(spec_param.fd,col,clip->src.s[index]);
            for (i=0;i<clip->src.s[index];i++)
            {
               (*((int *)data_pos)) = reverse[col[i]];
               data_pos += data_inc[index];
            }
         }
         else /* 16 bit file */
         {
            read(spec_param.fd,col,clip->src.s[index] * 2);
            for (i=0;i<2*clip->src.s[index];i += 2)
            {
               (*((int *)data_pos)) =
                   (int)(((int)col[i]+615)*reverse[col[i+1]]/721);
               data_pos += data_inc[index];
            }
         }
      }
/****** read 16 bit data from file ***************************************/
      else if(rw_mode == READ_16)
      {
         if(spec_param.file_type == FILE3D_8)
         {
            read(spec_param.fd,col,clip->src.s[index]);
            for (i=0;i<clip->src.s[index];i++)
            {
               *data_pos = 128;
               *(data_pos+1) = col[i];
               data_pos += data_inc[index];
            }
         }
         else /* 16 bit file */
         {
            read(spec_param.fd,col,clip->src.s[index] * 2);
            for (i=0;i<2*clip->src.s[index];i += 2)
            {
               *data_pos = col[i];
               *(data_pos+1) = col[i+1];
               data_pos += data_inc[index];
            }
         }
      }
/****** read 8 bit color data from file, scale size **********************/
      else if(rw_mode == READ_COLOR)
      {
         if(spec_param.file_type == FILE3D_8)
         {
            read(spec_param.fd,col,clip->src.s[index]);
            for (i=0;i<clip->src.s[index];i++)
            {
               if(clip->flags[index][i+clip->src.c[index]]) {
                  *data_pos = col[i];
                  data_pos += data_inc[index];
               }
            }
         }
         else /* 16 bit file */
         {
            read(spec_param.fd,col,clip->src.s[index] * 2);
            for (i=0;i<2*clip->src.s[index];i += 2)
            {
               if(clip->flags[index][(i>>1)+clip->src.c[index]]) {
                  *data_pos = col[i+1];
                  data_pos += data_inc[index];
               }
            }
         }
      }
/****** read integer data from file, scale size ***************************/
      else if(rw_mode == READ_SCALE)
      {
         if(spec_param.file_type == FILE3D_8)
         {
            read(spec_param.fd,col,clip->src.s[index]);
            for (i=0;i<clip->src.s[index];i++)
            {
               if(clip->flags[index][i+clip->src.c[index]]) {
                  (*((int *)data_pos)) = reverse[col[i]];
                  data_pos += data_inc[index];
               }
            }
         }
         else /* 16 bit file */
         {
            read(spec_param.fd,col,clip->src.s[index] * 2);
            for (i=0;i<2*clip->src.s[index];i += 2)
            {
               if(clip->flags[index][(i>>1)+clip->src.c[index]]) {
                  (*((int *)data_pos)) =
                      (int)(((int)col[i]+615)*reverse[col[i+1]]/721);
                  data_pos += data_inc[index];
               }
            }
         }
      }
/****** write data to file ************************************************/
      else if(rw_mode == WRITE_INT)
      {
         if(spec_param.file_type == FILE3D_8)
         {
            for (i=0;i<clip->src.s[index];i++)
            {
               n = (unsigned char)logr2(*((int *)data_pos),expon);
               write(spec_param.fd,&n,1);

               data_pos += data_inc[index];
            }
         }
         else /* 16 bit file */
         {
            for (i=0;i<clip->src.s[index];i++)
            {
               n = (unsigned char)logr2(*((int *)data_pos),expon);
               if(*((int *)data_pos))
               {
                  if(*((int *)data_pos) == 5 || *((int *)data_pos) == -5) a=128;
                  else a=(unsigned char)inint((float)(*((int *)data_pos) * 721)
                          / reverse[n] - 615);
               }
               else a=0;
               write(spec_param.fd,&a,1);
               write(spec_param.fd,&n,1);

               data_pos += data_inc[index];
            }
         }
      }
   }

/* Dimension > 0 ---> set positions and make recursive call **************/

   else
   {
      for (i=0;i<clip->src.s[index];i++)
      {
         if((rw_mode != READ_COLOR && rw_mode != READ_SCALE) || 
             clip->flags[index][i+clip->src.c[index]]) {
            Load_Data_Vector(clip,dim-1,data_pos,file_pos,data_inc,file_inc);
            data_pos += data_inc[index];
         }
	 file_pos += file_inc[index];
      }
   }
}

/*****************************************************************************
*      Write submatrix information from a node into the parameter file       *
*****************************************************************************/

void
Write_Sub_Pos(sub)

Subspectrum *sub;

{
   int i;

   for(i=0;i<spec_param.dim;i++)
      putw(sub->offset[i],spec_param.fp_param);
   putw(sub->pos,spec_param.fp_param);
}

/************************************************************************/
/* Write parameter file                                                 */
/************************************************************************/

int
Write_Param(param,version)

Spec_Param_Data *param;
int version;

{
   int i;

/* Write standard header to parameter file ********************************/

   if(version == 0) {
      fprintf(param->fp_param,"Version ....................... 0\n");
      fprintf(param->fp_param,"Number of dimensionen ......... %d\n",param->dim);
      if(param->file_type == FILE3D_8)
         fprintf(param->fp_param,"16 or 8 bit file type ......... 8\n");
      if(param->file_type == FILE3D_16)
         fprintf(param->fp_param,"16 or 8 bit file type ......... 16\n");
      for(i=0;i<param->dim;i++)
         fprintf(param->fp_param,
                 "Spectrometer frequency in w%d .. %f\n",i+1,param->w0[i]);
      for(i=0;i<param->dim;i++)
         fprintf(param->fp_param,
                 "Spectral sweep width in w%d .... %f\n",i+1,param->sw[i]);
      for(i=0;i<param->dim;i++)
         fprintf(param->fp_param,"Maximum chemical shift in w%d .. %f\n",i+1,
                                                            param->max_ppm[i]);
      for(i=0;i<param->dim;i++)
         fprintf(param->fp_param,"Size of spectrum in w%d ........ %d\n",i+1,
                                                            param->size[i]);
      fprintf(param->fp_param,
               "Type of spectrum .............. %c\n",param->type);
      fprintf(param->fp_param,"\n\n\nInformation for subspectra: \n");
      for(i=0;i<param->dim;i++)
         putw(param->sub_size[i],param->fp_param);
      for(i=0;i<param->dim;i++)   /* Order of dimensions  */
         putw(param->sub_perm[i],param->fp_param);

      spec_param = *param;   /* Use global variable to prevent parameter passing */
      spec_param.fp_param = param->fp_param;
      list_tree(spec_param.subs,Write_Sub_Pos);
   } else if (version == 1) {
/* Print version, file type and dimension *****************************/

      fprintf(param->fp_param,"Version ....................... 1\n");
      fprintf(param->fp_param,"Number of dimensions .......... %d\n",
	  	   param->dim);

      if(param->file_type == FILE3D_16)
         fprintf(param->fp_param,"16 or 8 bit file type ......... 16\n");
      else
         fprintf(param->fp_param,"16 or 8 bit file type ......... 8\n");

/* Write calibration data *************************/

      for (i=0;i<param->dim;i++)
         fprintf(param->fp_param,"Spectrometer frequency in w%d .. %f\n",i+1,
	       param->w0[i]);
      for (i=0;i<param->dim;i++)
         fprintf(param->fp_param,"Spectral sweep width in w%d .... %f\n",i+1,
               param->sw[i]);
      for (i=0;i<param->dim;i++)
        fprintf(param->fp_param,"Maximum chemical shift in w%d .. %f\n",i+1,
               param->max_ppm[i]);
      for (i=0;i<param->dim;i++)
        fprintf(param->fp_param,"Size of spectrum in w%d ........ %d\n",i+1,
               param->size[i]);
      for (i=0;i<param->dim;i++)
        fprintf(param->fp_param,"Submatrix size in w%d .......... %d\n",i+1,
               param->sub_size[i]);
      for (i=0;i<param->dim;i++)
        fprintf(param->fp_param,"Permutation for w%d ............ %d\n",i+1,
               param->sub_perm[i] + 1);
      for (i=0;i<param->dim;i++) {
        if(param->folding[i] == RSH)
           fprintf(param->fp_param,
                   "Folding in w%d ................. RSH\n",i+1);
        else if(param->folding[i] == TPPI)
           fprintf(param->fp_param,
                   "Folding in w%d ................. TPPI\n",i+1);
      }
      fprintf(param->fp_param,"Type of spectrum .............. %c\n",
               param->type);
   }

   return(1);
}

/************************************************************************/
/* Write empty data file                                                */
/************************************************************************/

int
Write_Data_Header(param,version)

Spec_Param_Data *param;
int version;

{
   int i;

   if(version == 0) {
      fprintf(param->fp,"Version 0");         /* Versionsnummer      */
      putw(param->dim,param->fp);   	      /* Anzahl Dimensionen  */
      for(i=0;i<param->dim;i++)
         putw(param->size[i],param->fp);      /* Size of spectrum    */
      if(param->file_type == FILE3D_8) putw(8,param->fp);
      if(param->file_type == FILE3D_16)putw(16,param->fp);
      for(i=0;i<param->dim;i++)
         putw(param->sub_size[i],param->fp);   /* Size of submatrix   */
      for(i=0;i<param->dim;i++)
         putw(param->sub_perm[i],param->fp);   /* Order of dimensions */
   }

   return(1);
}


/*****************************************************************************/
/** Calculate the logarithm of each data point and map onto the color scale **/
/*****************************************************************************/
int
logr2(num,expon)

register int num;
register int *expon;

{
	register i;
	register cnt;

	cnt = 16;
	i = 32;

	if(num >= 1) /* Positive number */
	{
   	   while((num < expon[i] || num >= expon[i+1]) && cnt != 0)
	   {
	      if(num > expon[i])
                i += cnt;
	      else
                i -= cnt;
              cnt >>= 1;
	   }

	   if((expon[i+1] - num) <= (num - expon[i]))
	      i++;

	   i++;	
	   if(i > 47) i = 47;
	}
	else if(num <= -1) /* Negative number */
	{
	   num = -num;

	   while((num < expon[i] || num >= expon[i+1]) && cnt != 0)
	   {
	      if(num > expon[i])
                i += cnt;
	      else
                i -= cnt;
              cnt >>= 1;
	   }

	   if((expon[i+1] - num) <= (num - expon[i]))
	      i++;

	   if(i > 47) i = 47;
	   i = 95 - i;
	}
	else if(num == 0)
	   i = 1;

	return(i);
}


int Calc_8bit_Int(a)
int a;
{
   return expon[a];
}
      
int Calc_16bit_Int(a,n)
int a;
int n;
{
   return((int)((a+615)*reverse[n]/721));
}