/*

Description: Tensor computation as eigensystem
             Class representing tensorfield as eigenvalues and eigenvectors;
             
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 _eigenfield
#define _eigenfield 1
#include <iostream.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include "basefield.h"
#include "eigen.h"
#include "tensorfield.h"
#include "background.h"
#include "vtkScalars.h"

#include <math.h>
#include <assert.h>

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

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

 public:
                                           // constructor for empty eigenfield  
  eigenfield(const int dimx=256, const int dimy=256, const int dimz = 1);
                                           // constructor that derives field from tensorfield  
  eigenfield(tensorfield &t, int offset[3], int window[3]);
                                           // constructor that loads field from disc
  eigenfield(char path[], const int dimx=256, const int dimy=256, const int dimz=1);

  ~eigenfield();                           // destructor

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

                                           // set eigen at position xyz to values given by e
  void set(const eigen e, const int x, const int y, const int z=0);
  void set(const eigen e, const int position);
  void setVectors(const float *vec1, const float *vec2, const float *vec3, const int position);

  void tensor2eigen(tensorfield &t);       // converts a tensorfield into an eigensystem 
  void eigen2tensor(tensorfield &t) const; // converts an eigensystem into a tensorfield
  
  void read_eigenslice(char path[], const int slice_no = 0);

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

  // functions
                                           // set values depending on property of eigensystem
  void measure2vtkScalars(int offset[3], int window[3], vtkScalars *scalars, tensorfield &t);

  // LABELING

                                           // returns a short array of same dimension as _field
                                           // which is a labelmap of _field
  int* structure(int* labels, const float tolerance = 0.1, const float threshold = 0.8, const int similarity = 0); 
  int* labelAll (int* labels, const float tolerance = 0.1, const float threshold = 0.8, const int similarity = 0);
  int* growLabel(int* labels, const float tolerance = 0.1, const float threshold = 0.8, const int similarity = 0);

  void setLabel (int* labels, const int &position, const float &tolerance, int &lastlabel);
  void setLabel1(int* labels, const int &position, const float &threshold, const float &tolerance, int &lastlabel, const int &A, const int &similarity);
  void setLabel2(int* labels, const int &position, const float &threshold, const float &tolerance, int &lastlabel, const int &A, const int &B, const int &similarity);
  void setLabel3(int* labels, const int &position, const float &threshold, const float &tolerance, int &lastlabel, const int &A, const int &B, const int &C, const int &similarity);
  void setLabel4(int* labels, const int &position, const float &threshold, const float &tolerance, int &lastlabel, const int &A, const int &B, const int &C, const int &D, const int &similarity);

                                           // reorder labels so that
                                           // labels[0] is background
                                           // labels[1] is rest group
                                           // labels[2] .. labels[lastlabel] are groups ordered by groupsize 
  int* reorder(int* labels, const int lastlabel);


                                           // merge to labels and decrease number of labels
  void merge2(int *labelmap, int labelA, int labelB, const int position, int &lastlabel);
  void merge3(int *labelmap, int labelA, int labelB, int labelC, const int position, int &lastlabel);
  void merge4(int *labelmap, int labelA, int labelB, int labelC, int labelD,const int position, int &lastlabel);

  void grow(const int position, int *labels, const float &tolerance, const int &similarity);

  int  check_neighbours(const int position, const int* labels, int& neighbour) const;
  int  neighbours(const int *labels, const int position) const;

  // FILTERS
  void smooth();                           // smooth the field (does not apply on eigenvectors)
  void median(const int filterlength);     // apply medianfilter on eigenvalues (does not apply on eigevectors)
  void remove_strips(const int inside=0);  // simpler than median: remove locations where data is 0 value 0

  void setZero();                          // set eigensystem to zero if !(all lambdas !=0)

  //eigenfield* downsample(eigenfield &e_field, const int ratio);  // set values of e_field so it represents a downsampled version of this field

  void correctVectors();                   // correct the eigenvectors; post: eigensystem is symmetric 
                                           // and all vectors are orthogonal
  void check();


  void mask(background &bg);                // set to zero where background is zero

  float *neighbours_direction(float *results); // values in array represent anisotropy of largest eigenvector in a neighbourhood
                                               // values range: 0 .. 16 (-8 .. 8 is offset by 8 so all values are >0 )

                                               // rotate all eigensystems by angles pitch, roll, yaw in grad
  void rotate(const float pitch, const float roll, const float yaw);
                                               // scale all eigensystems by factors scalex, scaley, scalez
  void scale( const float scalex, const float scaley, const float scalez);
};

#endif
