/*

Description: Tensorfields in 2D and 3D
When writing new functions: 
	     1) z-dimension should always default to 1 in a constructor and to 0 in a function, so that
	        tensorfield(256 256, 256); is the constructor for a 256*256*256 field of tensors
		get(0, 0, 0);              returns the first element in this example
		get(255, 255, 255);        returns the last element in this example
		
	     Other functions
	     -	tensorfield(path, 256, 256, 9);      
		reads 9 slices from a serie located in path. 
		Path ends with -S as files created by xpath(?) do
             
Author:      Raimundo Sierra
             www.rsierra.com

Copyright:   Copyright (c) 2000 Raimundo Sierra. All rights reserved.
LICENSE:     This program is free software; you can redistribute it and/or modify
             it under the terms of the GNU General Public License as published by
             the Free Software Foundation; either version 2 of the License, or
             (at your option) any later version.

             This program is distributed in the hope that it will be useful,
             but WITHOUT ANY WARRANTY; without even the implied warranty of
             MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
             GNU General Public License for more details.

             You should have received a copy of the GNU General Public License
             along with this program; if not, write to the Free Software
             Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Institution: Surgical Planning Laboratory
             Department of Radiology
             Brigham and Women's Hospital
             Harvard Medical School
             MA 02115
             USA

Date:        11.2000 - 3.2001

Modifications:
<Date>   <Init>    <Version>    <Description>

*/
#ifndef _tensorfield
#define _tensorfield 1
#include <stdlib.h>
#include <iostream.h>
#include <string.h>
#include <stdio.h>
#include "basefield.h"
#include "background.h"
#include "tensor.h"

#include "vtkTensors.h"
#include "vtkScalars.h"

#ifdef debug
  extern int fieldcounter;                 // counts calls to constructors
  extern int fielddestructorcounter;       // counts calls to destructors
#endif

class tensorfield: public basefield {
 private:
  tensor *_field;                          // pointer to the array of tensors, 3D representation:
                                           // [x + y*_dimx + z*_dimx*_dimz]

  tensorfield(const tensorfield &field);   // do not want a copy constructor to 
                                           // avoid creating unwanted copies of a tensorfield.

                                           // read single image from output of /d/lsdi/bin/LSDIrecon 
                                           // to compute tensorvalues, called from read_tensorslice
  //short* read_pic(short* image, char path_file[], int offset=0); 

 public:
  // constructors
                                           // constructor which reads raw data from disk
                                           // and computes the tensors
                                           // dimz corresponds to the number of slices 
                                           // if path is an xml-file then extract information
                                           // from the xml-file and ignore dimensions
  tensorfield(const char path[], const int dimx=256, const int dimy=256, const int dimz=1);

                                           // constructor that creates empty tensorfield
  tensorfield(const int dimx=256, const int dimy=256, const int dimz=1);

  // destructor
  ~tensorfield();
  
  // basic functions: storage, set and get values

                                           // get tensor at position xyz
  tensor* get(const int x, const int y, const int z=0);
  tensor* get(const int position);

                                           // each array contains on component of the tensorfield in the same order as the field
  void separate(float * array1, float * array2, float * array3, float * array4, float * array5, float * array6);
                                           // opposite function: merge arrays of size of field into tensorfield
  void merge   (float * array1, float * array2, float * array3, float * array4, float * array5, float * array6);
  void merge   (float * array1, float * array2, float * array3, float * array4, float * array5, float * array6, float * array7, float * array8, float * array9);

                                           // set tensor at position xyz to values given by t
  void set(const tensor t, const int x, const int y, const int z=0);

                                           // read images le1x to le3z for a single slice to 
                                           // derive tensors
  void read_tensorslice(const char path[], const int slice_no = 0);

  void load(const char path[]);             // load tensorfield from disc
                                            // save tensorfield to disc
  void save(const char path[], int offset[3]=NULL, int windowsize[3]=NULL); 

                                            // output tensor at position xyz; not really needed since 
                                            // tensor can be printed throu 
                                            // cout << myfield.get(int, int, int);
  void print(const int x, const int y, const int z=0) const;

                                           // transform this classes tensorfield into a vtk-Tensorfield
                                           // no boundarie-check on offset and window
  void tensorField2vtkField(const int offset[3], const int window[3], vtkTensors * newTensors) const;
                                           // create a scalar vtk field
                                           // no boundarie-check on offset and window
  void det2vtkScalars(const int offset[3], const int window[3], vtkScalars *newScalars) const;


  // functions which process the actual tensordata
                                            // remove the strips outside the object (inside=0)
  void remove_strips(const int inside=0);   // inside = 1 also removes the strips inside the object
  void smooth(const int filterlength=3);

                                            // downsample this field into smallfield
                                            // Boundaries are checked
                                            // Post: the tensorfield is smoothed!!
  void downsampleSlices(tensorfield &smallfield, const int factor=2, const int smoothlength=3);

  void roundness(float * pointarray);       // pointarray = det()/trace() at each point
                                            // pointarray = det()^(1/dim)/(trace()+sigma) at each point
  void roundnessCF(float * pointarray, const float thresh=0.01);    

  void trace(float * array) const;          // array = trace() at each point
  void det  (float * array) const;          // array = det()   at each point

  void mask(background &bg);
};

#endif