/* 

for documentation see display.h 


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.
*/


#include "display.h"

int imagecounter = 0;          // incremental counter to save temporary images

display::display(tensorfield &t_field)
{ 
  t_field.getDelta(_deltax, _deltay, _deltaz);
  t_field.getDimension(_dimx, _dimy, _dimz);
#ifdef debug
  cout << "Debug: Dimensions are: " <<_dimx   <<", "<<_dimy   <<", "<<_dimz<< endl;
  cout << "Debug: Voxelsizes are: " <<_deltax <<", "<<_deltay <<", "<<_deltaz<< endl;
#endif
}

display::display(eigenfield &e_field)
{ 
  e_field.getDelta(_deltax, _deltay, _deltaz);
  e_field.getDimension(_dimx, _dimy, _dimz);
#ifdef debug
  cout << "Debug: Dimensions are: " <<_dimx   <<", "<<_dimy   <<", "<<_dimz<< endl;
  cout << "Debug: Voxelsizes are: " <<_deltax <<", "<<_deltay <<", "<<_deltaz<< endl;
#endif
}

display::display(const int dimx=256, const int dimy=256, const int dimz=1, const float deltax=1.0, const float deltay=1.0, const float deltaz=1.0)
{
  _dimx = dimx;
  _dimy = dimy;
  _dimz = dimz;
  _deltax = deltax;
  _deltay = deltay;
  _deltaz = deltaz;
}

display::~display()
{
  system("rm tmpi_*");
}

int
display::color(short *image, const int dim1=256, const int dim2=256) const
{
  float max=0;
  float min=0;
  int tmp;
  char buffer[10];
  FILE *out_file;

  char tmppath[13];
  char command[18];
  command[0] = 'x';
  command[1] = 'v';
  command[2] = ' ';
  command[3] = tmppath[0] = 't';
  command[4] = tmppath[1] = 'm';
  command[5] = tmppath[2] = 'p';
  command[6] = tmppath[3] = 'i';
  command[7] = tmppath[4] = '_';
  command[8] = tmppath[5] = char( imagecounter/100    + 48);
  command[9] = tmppath[6] = char( imagecounter/10%10  + 48);
  command[10]= tmppath[7] = char( imagecounter%100%10 + 48);
  command[11]= tmppath[8] = '.';
  command[12]= tmppath[9] = 'p';
  command[13]= tmppath[10]= 'p';
  command[14]= tmppath[11]= 'm';
  command[15]= ' ';
  command[16]= '&';
  command[17]= tmppath[12]='\0';
  imagecounter++;

  out_file = fopen(tmppath, "w");
  if(out_file == NULL){
    cout << "\n\t Cannot open file "<<tmppath<< " for writing \tReturning without writing\n";
    return -1;
  }
  if(dim1>999 || dim2>999){
    cout << "Image dimension too big, edit function color() in display.cc\n";
    return -1;
  }
  for(int i=0; i<dim1*dim2*3; i++){
    if(float(image[i]) < min) min = float(image[i]);     // find minimum value
    if(float(image[i]) > max) max = float(image[i]);     // find max value
  }
  if(-min > max) max = -min;
  if(max == 0){
    cout << "WARNING: Maximum value is 0\nNothing to display, return"<<endl;
    fclose(out_file);
    return -1;
  }

  fwrite("P3\n",  sizeof(char), 3, out_file);
  buffer[0] = char( dim1/100    + 48);
  buffer[1] = char( dim1/10%10  + 48);
  buffer[2] = char( dim1%100%10 + 48);
  buffer[3] = ' ';
  buffer[4] = char( dim2/100    + 48);
  buffer[5] = char( dim2/10%10  + 48);
  buffer[6] = char( dim2%100%10 + 48); 
  buffer[7] = '\n';
  buffer[8] = '\0';
  fwrite(buffer,  sizeof(char), 8, out_file);
  fwrite("255\n", sizeof(char), 4, out_file);

  buffer[5] = '\0';

  for(int i=0; i<dim1*dim2*3; i++){
    tmp = int(abs(float(image[i])/max*255));        // tmp should be in range 0...255 now!
#ifdef debug
    if (tmp>255 || tmp<0){
      cout << "error in color() in display.cc, value out of range: "<<tmp<<endl;
    }
#endif
    buffer[0] = char( tmp/100    + 48);
    buffer[1] = char( tmp/10%10  + 48);
    buffer[2] = char( tmp%100%10 + 48);
    fwrite(buffer, sizeof(char), 4, out_file);
    if(i % 15 == 0){
      fwrite("\n", sizeof(char), 1, out_file);
    }
  }
  fclose(out_file);
  system(command);
  return 1;
}

int
display::greylevel(short* image, const int dim1=256, const int dim2=256) const
{
  float max=0;
  float min=0;
  int tmp;
  char buffer[10];
  FILE *out_file;
  char tmppath[13];
  char command[18];

  command[0] = 'x';
  command[1] = 'v';
  command[2] = ' ';
  command[3] = tmppath[0] = 't';
  command[4] = tmppath[1] = 'm';
  command[5] = tmppath[2] = 'p';
  command[6] = tmppath[3] = 'i';
  command[7] = tmppath[4] = '_';
  command[8] = tmppath[5] = char( imagecounter/100    + 48);
  command[9] = tmppath[6] = char( imagecounter/10%10  + 48);
  command[10]= tmppath[7] = char( imagecounter%100%10 + 48);
  command[11]= tmppath[8] = '.';
  command[12]= tmppath[9] = 'p';
  command[13]= tmppath[10]= 'p';
  command[14]= tmppath[11]= 'm';
  command[15]= ' ';
  command[16]= '&';
  command[17]= tmppath[12]='\0';
  imagecounter++;

  out_file = fopen(tmppath, "w");
  if(out_file == NULL){
    cout << "\n\t Cannot open file "<<tmppath<< " for writing \tReturning without writing\n";
    return -1;
  }
  if(dim1>999 || dim2>999){
    cerr << "Image dimension too big, edit function greylevel(short* ) in display.cc\n";
    return -1;
  }
  for(int i=0; i<dim1*dim2; i++){
    if(float(image[i]) < min) min = float(image[i]);     // find minimum value
    if(float(image[i]) > max) max = float(image[i]);     // find maximum value
  }
  if(-min > max) max = -min;
  if(max == 0){
    cout << "WARNING: Maximum value is 0\nNothing to display, return"<<endl;
    fclose(out_file);
    return -1;
  }

  fwrite("P2\n", sizeof(char), 3, out_file);             // write format
  buffer[0] = char( dim1/100    + 48);
  buffer[1] = char( dim1/10%10  + 48);
  buffer[2] = char( dim1%100%10 + 48);
  buffer[3] = ' ';
  buffer[4] = char( dim2/100    + 48);
  buffer[5] = char( dim2/10%10  + 48);
  buffer[6] = char( dim2%100%10 + 48); 
  buffer[7] = '\n';
  buffer[8] = '\0';
  fwrite(buffer,  sizeof(char), 8, out_file);            // write dimensions
  fwrite("255\n", sizeof(char), 4, out_file);            // write colorrange

  buffer[5] = '\0';
  for(int i=0; i<dim1*dim2; i++){
    tmp = int(abs(float(image[i])/max*255.0));           // tmp should be in range 0...255 now!
#ifdef debug
    if (tmp>255 || tmp<0){
      cout << "error in greylevel(short *) in display.cc, value out of range: "<<tmp<<endl;
    }
#endif
    buffer[0] = char( tmp/100    + 48);
    buffer[1] = char( tmp/10%10  + 48);
    buffer[2] = char( tmp%100%10 + 48);
    fwrite(buffer, sizeof(char), 4, out_file);
    if(i % 15 == 0){
      fwrite("\n", sizeof(char), 1, out_file);
    }
  }
  fclose(out_file);
  system(command);
  return 1;
}

int
display::greylevel(float* image, const int dim1=256, const int dim2=256) const
{
  float max=0;
  float min=0;
  int tmp;
  char buffer[10];
  FILE *out_file;
  char tmppath[13];
  char command[18];

  command[0] = 'x';
  command[1] = 'v';
  command[2] = ' ';
  command[3] = tmppath[0] = 't';
  command[4] = tmppath[1] = 'm';
  command[5] = tmppath[2] = 'p';
  command[6] = tmppath[3] = 'i';
  command[7] = tmppath[4] = '_';
  command[8] = tmppath[5] = char( imagecounter/100    + 48);
  command[9] = tmppath[6] = char( imagecounter/10%10  + 48);
  command[10]= tmppath[7] = char( imagecounter%100%10 + 48);
  command[11]= tmppath[8] = '.';
  command[12]= tmppath[9] = 'p';
  command[13]= tmppath[10]= 'p';
  command[14]= tmppath[11]= 'm';
  command[15]= ' ';
  command[16]= '&';
  command[17]= tmppath[12]='\0';
  imagecounter++;

  out_file = fopen(tmppath, "w");
  if(out_file == NULL){
    cout << "\n\t Cannot open file "<<tmppath<< " for writing \tReturning without writing\n";
    return -1;
  }
  if(dim1>999 || dim2>999){
    cout << "Image dimension too big, edit function greylevel(float* ) in display.cc\n";
    return -1;
  }
  for(int i=0; i<dim1*dim2; i++){
    if(image[i] < min) min      =  image[i];             // find minimum value
    if(image[i] > max) max      =  image[i];             // find maximum value
  }
  if(-min > max) max = -min;
  if(max == 0){
    cout << "WARNING: Maximum value is 0\nNothing to display, return"<<endl;
    fclose(out_file);
    return -1;
  }

  fwrite("P2\n", sizeof(char), 3, out_file);             // write format
  buffer[0] = char( dim1/100    + 48);
  buffer[1] = char( dim1/10%10  + 48);
  buffer[2] = char( dim1%100%10 + 48);
  buffer[3] = ' ';
  buffer[4] = char( dim2/100    + 48);
  buffer[5] = char( dim2/10%10  + 48);
  buffer[6] = char( dim2%100%10 + 48); 
  buffer[7] = '\n';
  buffer[8] = '\0';
  fwrite(buffer,  sizeof(char), 8, out_file);            // write dimensions
  fwrite("255\n", sizeof(char), 4, out_file);            // write colorrange

  buffer[5] = '\0';
  for(int i=0; i<dim1*dim2; i++){
    tmp = abs(int(image[i]/max*255.0));                  // tmp should be in range 0...255 now!
#ifdef debug
    if (tmp>255 || tmp<0){ 
      cout << "error in greylevel(float* ) in display.cc, value out of range: "<<tmp<<endl;
    }
#endif
    buffer[0] = char( tmp/100    + 48);
    buffer[1] = char( tmp/10%10  + 48);
    buffer[2] = char( tmp%100%10 + 48);
    fwrite(buffer, sizeof(char), 4, out_file);
    if(i % 15 == 0){
      fwrite("\n", sizeof(char), 1, out_file);
    }
  }
  fclose(out_file);
  system(command);
  return 1;
  /*
  char * number = new char[10];
  strcpy(number, "tmp.");
  sprintf(buffer, "%03d", imagecounter);
  strcat(number, buffer);
  //cout << "File "<<number<<" will be written\n";
  for(int i=0; i<dim1*dim2; i++){
    image[i] = abs(int(image[i]/max*255.0));
  }
  saveFloatAsMRI(number, image, dim1*dim2);
  delete[] number;
  */
}

float
display::backgd(background &bac, int offset[3], int window[3], const int selection, float thresh=0.0)
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  int max =0;
  float lastthresh = 0;
  bac.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  short * image = new short[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	if     (selection==0){ image[ipos]   = bac.getAAI(fpos);}  // greylevel-value
	else if(selection==1){ image[ipos]   = bac.getADC(fpos);}
	else if(selection==2){ image[ipos]   = bac.getDWI(fpos);}
	else if(selection==3){ image[ipos]   = bac.getRAI(fpos);}
	else if(selection==4){ image[ipos]   = bac.getT2W(fpos);}
	else {
	  cout << "Invalid selection of background, returning without displaying\n";
	  delete[] image;
	  return 0;
	}
	if (image[ipos] > max) max = image[ipos];
      }
    }
    lastthresh = thresh*max;
    for( int i=0; i<plane[0]*plane[1]; i++){
      if(image[i] < lastthresh) image[i] = 0;
    }
    max = 0;
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;
  return lastthresh;
}

void 
display::lambda1(eigenfield  &e_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*e_field.get(fpos)).value1();                    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;
}

void 
display::lambda2(     eigenfield  &e_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*e_field.get(fpos)).value2();                    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;
}
  
void 
display::lambda3(     eigenfield  &e_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*e_field.get(fpos)).value3();                    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;
}
 
void 
display::vector1(     eigenfield  &e_field, int offset[3], int window[3], const int element=0) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*e_field.get(fpos)).vector1(element);                    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
}
  
void 
display::vector2(     eigenfield  &e_field, int offset[3], int window[3], const int element=0) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*e_field.get(fpos)).vector2(element);                    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;
}
  
void 
display::vector3(     eigenfield  &e_field, int offset[3], int window[3], const int element=0) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*e_field.get(fpos)).vector3(element);                    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;
}

void
display::max_eigenval(eigenfield  &e_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0; 
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image
  short * image = new short[3*plane[0]*plane[1]];
  for (int k=0; k<plane[2]; k++){                       // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                           // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                           // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                           // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = 3*(i+j*plane[0]);
	image[ipos]   = short( (*e_field.get(fpos)).value1() * (*e_field.get(fpos)).vector1(0) );  // r-value
	image[ipos+1] = short( (*e_field.get(fpos)).value1() * (*e_field.get(fpos)).vector1(1) );  // g-value
	image[ipos+2] = short( (*e_field.get(fpos)).value1() * (*e_field.get(fpos)).vector1(2) );  // b-value
      }
    }
    // display the new image
    color(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL; 
}

void 
display::greyarrayDifference(float *array1, float *array2, const int dim1=256, const int dim2=256, const int dim3=1)
{
  float * image = new float[dim1*dim2];
  for(int i=0; i<dim3; i++){
    for(int j=0; j<dim1*dim2; j++){
      image[j] = array1[j+i*dim1*dim2]-array2[j+i*dim1*dim2];
    }
    greylevel(image, dim1, dim2);
  }
  delete[] image;
  image = NULL;
}

void 
display::greyarray(float *myarray, const int dim1=256, const int dim2=256, const int dim3=1)
{
  float * image = new float[dim1*dim2];
  for(int i=0; i<dim3; i++){
    for(int j=0; j<dim1*dim2; j++){
      image[j] = myarray[j+i*dim1*dim2];
    }
    greylevel(image, dim1, dim2);
  }
  delete[] image;
  image = NULL;
}


void 
display::greyarray(short *myarray, const int dim1=256, const int dim2=256, const int dim3=1)
{
  short * image = new short[dim1*dim2];
  for(int i=0; i<dim3; i++){
    for(int j=0; j<dim1*dim2; j++){
      image[j] = myarray[j+i*dim1*dim2];
    }
    greylevel(image, dim1, dim2);
  }
  delete[] image;
  image = NULL;
}

void 
display::greyarray(unsigned short *myarray, const int dim1=256, const int dim2=256, const int dim3=1)
{
  float * image = new float[dim1*dim2];
  for(int i=0; i<dim3; i++){
    for(int j=0; j<dim1*dim2; j++){
      image[j] = myarray[j+i*dim1*dim2];
    }
    greylevel(image, dim1, dim2);
  }
  delete[] image;
  image = NULL;
}

void 
display::closeLine(eigenfield  &e_field, int offset[3], int window[3], float thresh=0.8, float tol=0.1) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  float max = 0;
  float min = 1;

  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   =  (*e_field.get(fpos)).close_line();                    // greylevel-value
	if (image[ipos]>max) max = image[ipos];
	if (image[ipos]<min) min = image[ipos];	
      }
    }
    // display the new image
    thresh = thresh*max;
    for (int j=0; j<plane[1]*plane[0]; j++){
      if(image[j]!=0 && image[j] < thresh) image[j] = min;
    }
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;  
}

void 
display::closePlane(eigenfield  &e_field, int offset[3], int window[3], float thresh=0.8, float tol=0.1) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  float max = 0;
  float min = 1;

  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   =  (*e_field.get(fpos)).close_plane();                    // greylevel-value
	if (image[ipos]>max) max = image[ipos];
	if (image[ipos]<min) min = image[ipos];	
      }
    }
    // display the new image
    thresh = thresh*max;
    for (int j=0; j<plane[1]*plane[0]; j++){
      if(image[j]!=0 && image[j] < thresh) image[j] = min;
    }
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;
  image = NULL;  
}

void 
display::closeSphere(eigenfield  &e_field, int offset[3], int window[3], float thresh=0.8, float tol=0.1) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  float max = 0;
  float min = 1;

  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   =  (*e_field.get(fpos)).close_sphere();                    // greylevel-value
	if (image[ipos]>max) max = image[ipos];
	if (image[ipos]<min) min = image[ipos];	
      }
    }
    // display the new image
    thresh = thresh*max;
    for (int j=0; j<plane[1]*plane[0]; j++){
      if(image[j]!=0 && image[j] < thresh) image[j] = min;
    }
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image; 
  image = NULL; 
}

void 
display::anisotropy(eigenfield  &e_field, int offset[3], int window[3], float thresh=0.8, float tol=0.1) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  float max = 0;
  float min = 1;

  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   =  (*e_field.get(fpos)).anisotropy();                    // greylevel-value
	if (image[ipos]>max) max = image[ipos];
	if (image[ipos]<min) min = image[ipos];	
      }
    }
    // display the new image
    thresh = thresh*max;
    for (int j=0; j<plane[1]*plane[0]; j++){
      if(image[j]!=0 && image[j] < thresh) image[j] = min;
    }
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image; 
  image = NULL; 
}

void 
display::element(tensorfield &t_field, int offset[3], int window[3], const int elementx, const int elementy) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  t_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   =  (*t_field.get(fpos)).get(elementx, elementy);    // greylevel-value
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image;  
  image = NULL;
}

void 
display::det(tensorfield &t_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  t_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*t_field.get(fpos)).det();                    // greylevel-value
      }
    }
    // display the new image
    if(greylevel(image, plane[0], plane[1])==-1){
      cout << "In display.det()\n";
    }
  }
  delete[] image;  
  image = NULL;
} 


void 
display::trace(tensorfield &t_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  t_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = (*t_field.get(fpos)).trace();                    // greylevel-value
      }
    }
    // display the new image
    if(greylevel(image, plane[0], plane[1])==-1){
      cout << "In display.trace()\n";
    }
  }
  delete[] image;  
  image = NULL;
} 

void 
display::euclid(tensorfield &t_field, int offset[3], int window[3]) const
{
  int fpos=0;
  int ipos=0;
  int plane[4];
  t_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  short * image = new short[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	image[ipos]   = short( (*t_field.get(fpos)).euclid());                    // greylevel-value
      }
    }
    // display the new image
    if(greylevel(image, plane[0], plane[1])==-1){
      cout << "In display.euclid()\n";
    }
  }
  delete[] image; 
  image = NULL; 
}

void
display::classes_cl_cp_cs(eigenfield  &e_field, int offset[3], int window[3])
{
  int fpos=0;
  int ipos=0;
  int plane[4];

  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  // create image  
  float cl = 0;
  float cp = 0;
  float cs = 0;

  float * image = new float[plane[0]*plane[1]];
  for( int k=0; k<plane[2]; k++){                      // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                          // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                          // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                          // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	cl   =  (*e_field.get(fpos)).close_line();                    // greylevel-value
	cp   =  (*e_field.get(fpos)).close_plane();
	cs   =  (*e_field.get(fpos)).close_sphere();
	if((*e_field.get(fpos)).isZero()) image[ipos]   = 0;
	else if(cl >= cp && cl >= cs) image[ipos]   = 240;
	else if(cp >= cl && cp >= cs) image[ipos]   = 110;
	else if(cs >= cl && cs >= cp) image[ipos]   = 20;
	else image[ipos]   = 0;
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
  delete[] image; 
  image = NULL; 
}

void 
display::labels(eigenfield  &e_field, int offset[3], int window[3], float thresh=0.8, float tol=0.1, const int method=0, const int similarity=0, const int colors=0) const
{ 
  int plane[4];
  e_field.check_boundaries(offset, window, plane);

  if(plane[2]>20){
    cout << "WARNING: do not display more than 20 images;\nClipping to 20 images\n";
    plane[2]=20;
  }
  short * image = NULL;
  int * label   = NULL;

  switch(method){
  case 0:  label = e_field.structure(label, tol, thresh, similarity); break;
  case 1:  label = e_field.labelAll( label, tol, thresh, similarity); break;
  default: label = e_field.growLabel(label, tol, thresh, similarity); break;    
  }

  if(colors == 1) {
    image = new short[3*plane[0]*plane[1]];  // create image
    this->colorlabels(image, label, offset, plane);
  }
  else if(colors == 0) {
    image = new short[plane[0]*plane[1]];  // create image
    this->greylabels( image, label, offset, plane);
  }

  delete[] label;
  label = NULL;
  delete[] image;
  image = NULL;
}

void 
display::colorlabels(short * image, int * label, int * offset, int * plane) const
{
  int fpos=0;
  int ipos=0;
  for (int k=0; k<plane[2]; k++){                       // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                           // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                           // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                           // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = 3*(i+j*plane[0]);
	if(label[fpos]==0){
	  image[ipos]   = 0;    // r-value
	  image[ipos+1] = 0;    // g-value
	  image[ipos+2] = 0;    // b-value
	}
	else if(label[fpos]==1){
	  image[ipos]   = 100;    // r-value
	  image[ipos+1] = 100;    // g-value
	  image[ipos+2] = 100;    // b-value
	}
	else if(label[fpos]<=100){
	  image[ipos]   = 0;    // r-value
	  image[ipos+1] = 255;  // g-value
	  image[ipos+2] = 0;    // b-value
	}
	else if(label[fpos]<=200){
	  image[ipos]   = 255;    // r-value
	  image[ipos+1] = 0;      // g-value
	  image[ipos+2] = 0;      // b-value
	}
	else if(label[fpos]<=300){
	  image[ipos]   = 0;   
	  image[ipos+1] = 0; 
	  image[ipos+2] = 255;  
	} 
	else if(label[fpos]<=1000){
	  image[ipos]   = 255;   
	  image[ipos+1] = 0; 
	  image[ipos+2] = 255;  
	} 
	else {
	  image[ipos]   = 100;    // r-value
	  image[ipos+1] = 100;  // g-value
	  image[ipos+2] = 100;    // b-value
	} 
      }
    }
    // display the new image
    color(image, plane[0], plane[1]);
  }
}
void 
display::greylabels(short * image, int * label, int * offset, int * plane) const
{
  int fpos=0;
  int ipos=0;
  for (int k=0; k<plane[2]; k++){                       // for all images
    for (int j=0; j<plane[1]; j++){ 
      for (int i=0; i<plane[0]; i++){
	if     (plane[3]==2){                           // xy plane
	  fpos = (i+offset[0]) + (j+offset[1])*_dimx + (k+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==1){                           // xz plane
	  fpos = (i+offset[0]) + (k+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	} 
	else if(plane[3]==0){                           // yz plane
	  fpos = (k+offset[0]) + (i+offset[1])*_dimx + (j+offset[2])*_dimx*_dimy;
	}
	ipos = i+j*plane[0];
	if(label[fpos]<=32767){  // limited to short numbers ... display the largest 255 groups
	  image[ipos]   = label[fpos];    
	}
	else {
	  image[ipos]   = 1;
	}
      }
    }
    // display the new image
    greylevel(image, plane[0], plane[1]);
  }
}

void
display::myVtkTensor(tensorfield &t_field, int offset[3], int window[3], vtkScalars *scalars=NULL, const int resolution=6, const int quality=0, const int background=0)
{
  t_field.check_boundaries(offset, window);

  float origx , origy , origz;
  float spacingx, spacingy, spacingz;

  if(_deltax>=0)  spacingx = 10*_deltax; else spacingx = 10;
  if(_deltay>=0)  spacingy = 10*_deltay; else spacingy = 10;
  if(_deltaz>=0)  spacingz = 10*_deltaz; else spacingz = 10;

  origx = - ((window[0]-1)*spacingx)/2;
  origy = - ((window[1]-1)*spacingy)/2;
  origz = - ((window[2]-1)*spacingz)/2;
  #ifdef debug
  cout << "Spacing is: ("<< spacingx <<", "<<spacingy<<", "<<spacingz<<")\n";
  cout << "Origin is:  ("<< origx <<", "<<origy<<", "<<origz<<")\n";
  cout << "Window is:  ("<<window[0]<<", "<<window[1]<<", "<<window[2]<<")\n";
  cout << "Offset is:  ("<<offset[0]<<", "<<offset[1]<<", "<<offset[2]<<")\n";
  cout << flush;
  #endif
  // create a window, renderer and interactor
  vtkRenderWindow *renWin           = vtkRenderWindow::New();
  vtkRenderer *ren1                 = vtkRenderer::New();
  renWin->AddRenderer(ren1);
  renWin->SetSize(800, 800);

  vtkRenderWindowInteractor *iren   = vtkRenderWindowInteractor::New();
  iren->SetRenderWindow(renWin);

  // create sphere geometry
  vtkSphereSource *sphere           = vtkSphereSource::New();
  sphere->SetThetaResolution(resolution);
  sphere->SetPhiResolution(resolution);
  
  // create tensorstructure
  vtkTensors *newTensors            = vtkTensors::New();
  vtkStructuredPoints *volume       = vtkStructuredPoints::New();

  volume->SetDimensions(window[0], window[1], window[2]);
  volume->SetOrigin(origx, origy, origz);
  volume->SetSpacing(spacingx, spacingy, spacingz);

  // transform
  t_field.tensorField2vtkField( offset, window, newTensors);

  volume->GetPointData()->SetTensors(newTensors);
  volume->GetPointData()->SetScalars(scalars);

  vtkTensorGlyph *ellipsoids        = vtkTensorGlyph::New();
  ellipsoids->SetInput(volume);
  ellipsoids->SetSource(sphere->GetOutput());
  ellipsoids->SetScaleFactor(0.007);
  ellipsoids->ClampScalingOn();
  
  vtkPolyDataNormals *ellipNormals  = vtkPolyDataNormals::New();
  ellipNormals->SetInput(ellipsoids->GetOutput());

  // Map contour
  vtkLookupTable *lut               = vtkLookupTable::New();
  if(scalars!=NULL){
    lut->SetNumberOfColors(3);
    lut->Build();
    if(background==0){                                 //black
      lut->SetTableValue(0, 0.0, 0.0, 1.0, 0.2);
      lut->SetTableValue(1, 0.0, 0.0, 1.0, 1);
      lut->SetTableValue(2, 1.0, 0.8, 0.0, 1);
    }
    else {                                             //white
      lut->SetTableValue(0, 1.0, 0.8, 0.0, 0.2);
      lut->SetTableValue(1, 1.0, 0.8, 0.0, 1);
      lut->SetTableValue(2, 0.0, 0.0, 1.0, 1);
    }
  }
  
  vtkPolyDataMapper *ellipMapper    = vtkPolyDataMapper::New();
  ellipMapper->SetInput(ellipNormals->GetOutput()); // if you change ellipNormals to ellipsoids 
                                                    // it gets faster but with errors: some tensors are black
  //ellipMapper->ImmediateModeRenderingOn();
  ellipMapper->SetScalarRange(volume->GetScalarRange());
  if(scalars==NULL){    ellipMapper->ScalarVisibilityOff(); }
  else             {    ellipMapper->SetLookupTable(lut);   }

  vtkActor *ellipActor              = vtkActor::New();
  vtkLight *light                   = vtkLight::New();
  ellipActor->SetMapper(ellipMapper);
  if(background==0){
    ellipActor->GetProperty()->SetColor(1,0.8,0);    //black
  }
  else {
    ellipActor->GetProperty()->SetColor(1.0,1.0,0.0);  //white
  }
  if(quality == 1){
    if(background==0){
      ellipActor->GetProperty()->SetDiffuse(0.7);    //black
      ellipActor->GetProperty()->SetSpecular(1.0);
      ellipActor->GetProperty()->SetSpecularPower(20);
    }
    else {
      ellipActor->GetProperty()->SetDiffuse(0.5);    //white
      ellipActor->GetProperty()->SetSpecular(1.0);
      ellipActor->GetProperty()->SetSpecularPower(30);
    }
    //light->SetPosition(0, 0, 0);
    //light->SetFocalPoint(1.0, 0, 0);
    //ren1->AddLight(light);
  }


  ren1->AddActor(ellipActor);
  // Background color white
  if(background==0){
    ren1->SetBackground(0,0,0);    //black
  } 
  else ren1->SetBackground(1,1,1); //white
  // set camera and lightings and so on ...
  ren1->GetActiveCamera()->Zoom(1.5);
  ren1->GetActiveCamera()->Azimuth(180);
  ren1->GetActiveCamera()->Roll(180);
  //ren1->GetActiveCamera()->Elevation(180);
  // Begin mouse interaction
  iren->Start();

  // clean up
  light->Delete();
  renWin->Delete();
  ren1->Delete();
  iren->Delete();
  sphere->Delete();
  newTensors->Delete();
  volume->Delete();
  ellipsoids->Delete();
  ellipNormals->Delete();
  lut->Delete();
  ellipMapper->Delete();
  ellipActor->Delete();
  //scalars->Delete();
}

/*
void
display::simpleVtkWithBack()
{

}
*/

void
display::myVtkFibers(tensorfield &t_field, int offset[3], int window[3], const int maxNumberOfFibers=1)
{
  t_field.check_boundaries(offset, window);

  float origx , origy , origz;
  float spacingx, spacingy, spacingz;

  if(_deltax>=0)  spacingx = 10*_deltax; else spacingx = 10;
  if(_deltay>=0)  spacingy = 10*_deltay; else spacingy = 10;
  if(_deltaz>=0)  spacingz = 10*_deltaz; else spacingz = 10;

  origx = - ((window[0]-1)*spacingx)/2;
  origy = - ((window[1]-1)*spacingy)/2;
  origz = - ((window[2]-1)*spacingz)/2;
#ifdef debug
  cout << "Spacing is: ("<< spacingx <<", "<<spacingy<<", "<<spacingz<<")\n";
  cout << "Origin is:  ("<< origx <<", "<<origy<<", "<<origz<<")\n";
  cout << "Window is:  ("<<window[0]<<", "<<window[1]<<", "<<window[2]<<")\n";
  cout << "Offset is:  ("<<offset[0]<<", "<<offset[1]<<", "<<offset[2]<<")\n";
#endif
  // create a window, renderer and interactor
  vtkRenderWindow *renWin           = vtkRenderWindow::New();
  vtkRenderer *ren1                 = vtkRenderer::New();
  renWin->AddRenderer(ren1);
  renWin->SetSize(800, 800);

  vtkRenderWindowInteractor *iren   = vtkRenderWindowInteractor::New();
  iren->SetRenderWindow(renWin);

  // create sphere geometry
  vtkSphereSource *sphere           = vtkSphereSource::New();
  sphere->SetThetaResolution(8);
  sphere->SetPhiResolution(8);
  
  // create tensorstructure
  vtkTensors *newTensors            = vtkTensors::New();
  vtkStructuredPoints *volume       = vtkStructuredPoints::New();

  volume->SetDimensions(window[0], window[1], window[2]);
  volume->SetOrigin(origx, origy, origz);
  volume->SetSpacing(spacingx, spacingy, spacingz);

  
  // transform
  t_field.tensorField2vtkField( offset, window, newTensors);

  volume->GetPointData()->SetTensors(newTensors);

  vtkPolyDataCollection *appendAll  = vtkPolyDataCollection::New();
  vtkPolyDataMapper *s1mapper       = vtkPolyDataMapper::New();
  vtkHyperStreamline *s1            = vtkHyperStreamline::New();
  vtkActor *s1Actor                 = vtkActor::New();

  for( int i=0; i<2; i++){
  s1->SetInput(volume);
  s1->SetStartPosition(0, 2*i, origz);  
  s1->IntegrateMinorEigenvector();
  s1->SetMaximumPropagationDistance(150.0);
  s1->SetIntegrationStepLength(0.1);
  s1->SetStepLength(0.01);
  s1->SetRadius(0.25);
  s1->SetNumberOfSides(18);
  s1->SetIntegrationDirection(2);
  //s1->SetTerminalEigenvalue();
  s1->Update();

  //appendAll->AddItem(s1->GetOutput());
  }

  s1mapper->SetInput(s1->GetOutput());
  s1mapper->ScalarVisibilityOff();
  s1Actor->SetMapper(s1mapper);  
  s1Actor->GetProperty()->SetColor(1,0.8,0);

  ren1->AddActor(s1Actor);
  // Create outline around data
  vtkOutlineFilter *outline         = vtkOutlineFilter::New();
  outline->SetInput(volume);

  vtkPolyDataMapper *outlineMapper  = vtkPolyDataMapper::New();
  outlineMapper->SetInput(outline->GetOutput());

  vtkActor *outlineActor = vtkActor::New();
  outlineActor->SetMapper(outlineMapper);
  outlineActor->GetProperty()->SetColor(0, 0, 0);
  ren1->AddActor(outlineActor);
  
  // Background color white
  ren1->SetBackground(1,1,1);
  // set camera and lightings and so on ...
  ren1->GetActiveCamera()->Zoom(1.5);
  ren1->GetActiveCamera()->Roll(180);
  ren1->GetActiveCamera()->Azimuth(180);
  // Begin mouse interaction
  iren->Start();

  // clean up
  renWin->Delete();
  ren1->Delete();
  iren->Delete();
  sphere->Delete(); 
  newTensors->Delete();
  volume->Delete();
  s1mapper->Delete();
  s1Actor->Delete();
  outline->Delete();
  outlineMapper->Delete();
  outlineActor->Delete();
  appendAll->Delete();
}