/*
  for ducumentation see background.h

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 "background.h"

background::background(const char path[], const int dimx=256, const int dimy=256, const int dimz=1)
{
  int length = strlen(path);
  if(path[length-3] =='x' && path[length-2] =='m' && path[length-1] =='l'){ // reading xml-header
    cout << "Initializing background from xml-file\t\t"<< flush;
    this->read_XMLHeader(path);
    _AAI = new unsigned short[_dimx*_dimy*_dimz];
    _ADC = new unsigned short[_dimx*_dimy*_dimz];
    _DWI = new unsigned short[_dimx*_dimy*_dimz];
    _RAI = new unsigned short[_dimx*_dimy*_dimz];
    _T2W = new unsigned short[_dimx*_dimy*_dimz];
    cout << "done\n" << flush;

    length = strlen(_datafile);
    char *datapath = new char[length+5];
    strcpy(datapath, _datafile);
    datapath[length] = '.';
    datapath[length+1] = 'b';
    datapath[length+2] = 'g';
    datapath[length+3] = '\0';
    //cout <<datapath<<endl;
    this->load(datapath);
    return;
  }
  cout << "Initializing background \t\t\t" << flush;                      // reading raw data
  _dimx = dimx;
  _dimy = dimy;
  _dimz = dimz;
  _deltax = _deltay = _deltaz = 1;
  _AAI = new unsigned short[_dimx*_dimy*_dimz];
  _ADC = new unsigned short[_dimx*_dimy*_dimz];
  _DWI = new unsigned short[_dimx*_dimy*_dimz];
  _RAI = new unsigned short[_dimx*_dimy*_dimz];
  _T2W = new unsigned short[_dimx*_dimy*_dimz];
  cout << "done\n" << flush;

  char tmp[length + 4];
  char number[4];
  _sourcefile = new char[length];
  strcpy(_sourcefile, path);
  
  for(int i=1; i<=_dimz; i++){
    if(i<100){
      number[0] = char(i/10+48);
      number[1] = char(i%10+48);
      number[2] = '\0';
    }
    else if(i<1000){
      number[0] = char(i/100+48);
      number[1] = char((i/10)%10+48);
      number[2] = char(i%100%10+48);
      number[3] = '\0';
    }
    else {
      cerr << "Too many slices, cannot process in background::background(char, int, int, int)\nEXIT\n";
      exit(0);
    }
    strcat(strcpy(tmp, path), number);
    read_back(tmp, (i-1));
  }
}

background::background(const int dimx=256, const int dimy=256, const int dimz=1)
{
  _AAI = new unsigned short[_dimx*_dimy*_dimz];
  _ADC = new unsigned short[_dimx*_dimy*_dimz];
  _DWI = new unsigned short[_dimx*_dimy*_dimz];
  _RAI = new unsigned short[_dimx*_dimy*_dimz];
  _T2W = new unsigned short[_dimx*_dimy*_dimz];
  
  for(int i=0; i<_dimx*_dimy*_dimz; i++){
    _AAI[i] = 0;
    _ADC[i] = 0;
    _DWI[i] = 0;
    _RAI[i] = 0;
    _T2W[i] = 0;
  } 
}

background::~background()
{
  delete[] _AAI;
  delete[] _ADC;
  delete[] _DWI;
  delete[] _RAI;
  delete[] _T2W;
  _AAI = NULL;
  _ADC = NULL;
  _DWI = NULL;
  _RAI = NULL;
  _T2W = NULL;
}

void
background::read_back(const char path[], const int slice_no = 0)
{
  unsigned short *aai = new unsigned short[_dimx*_dimy];
  unsigned short *adc = new unsigned short[_dimx*_dimy];
  unsigned short *dwi = new unsigned short[_dimx*_dimy];
  unsigned short *rai = new unsigned short[_dimx*_dimy];
  unsigned short *t2w = new unsigned short[_dimx*_dimy];
  int pos = 0;
  int pos2= 0;

  cout <<"Reading background for slice no. " << slice_no << "\t\t"<<flush;
  char tmp[strlen(path)+8];
  
  strcat(strcpy(tmp, path), "AAI.pic");
  aai = read_pic(aai, tmp);
  strcat(strcpy(tmp, path), "ADC.pic");
  adc = read_pic(adc, tmp);
  strcat(strcpy(tmp, path), "DWI.pic");
  dwi = read_pic(dwi, tmp);
  strcat(strcpy(tmp, path), "RAI.pic");
  rai = read_pic(rai, tmp);
  strcat(strcpy(tmp, path), "T2W.pic");
  t2w = read_pic(t2w, tmp);

  for(int j=0; j < _dimx; j++){
    for(int i=0; i < _dimy; i++){ 
      pos = i+j*_dimx;
      pos2= slice_no*_dimx*_dimy + pos;
      _AAI[pos2] = aai[pos];
      _ADC[pos2] = adc[pos];
      _DWI[pos2] = dwi[pos];
      _RAI[pos2] = rai[pos];
      _T2W[pos2] = t2w[pos];
    }
  } 

  cout <<"done\n" << flush;
  delete[] aai; 
  delete[] adc;
  delete[] dwi; 
  delete[] rai;
  delete[] t2w;

}


void
background::load(const char path[])
{
  cout << "Reading background from disc"<<flush;
  FILE *in_file;
  in_file = fopen(path, "r");
  if(in_file == NULL) {
    cerr << "\nERROR: Cannot open " << path << "\tExit\n";
    exit(0);
  } else {
    fseek(in_file, 0, 0);
    fread(_AAI, sizeof(unsigned short), _dimx*_dimy*_dimz, in_file);
    fread(_ADC, sizeof(unsigned short), _dimx*_dimy*_dimz, in_file);
    fread(_DWI, sizeof(unsigned short), _dimx*_dimy*_dimz, in_file);
    fread(_RAI, sizeof(unsigned short), _dimx*_dimy*_dimz, in_file);
    fread(_T2W, sizeof(unsigned short), _dimx*_dimy*_dimz, in_file);
  }
  fclose(in_file);  
  cout << "\t\t\tdone\n"<<flush;
}

void
background::save(const char path[], int offset[3]=NULL, int windowsize[3]=NULL)
{
  cout << "Writing background to disc"<<flush;
  int length = strlen(path);
  char * filename = new char[length+5];
  char * xmlpath  = new char[length+5];
  strcpy(filename, path);
  strcpy(xmlpath,  path);
  strcpy(_datafile,path);
  if(path[length-3] =='x' && path[length-2] =='m' && path[length-1] =='l'){ //  xmlfile
    filename[length-3] = 'b';
    filename[length-2] = 'g';
    filename[length-1] = '\0';
  }
  else if(path[length-3] =='.' && path[length-2] =='b' && path[length-1] =='g'){// background file
    xmlpath[length-2]  = 'x';
    xmlpath[length-1]  = 'm';
    xmlpath[length  ]  = 'l';
    xmlpath[length+1]  = '\0';
    _datafile[length-3]= '\0';    
  }
  else { // no extensions => add extensions for background file (not needed for xmlfile)
    filename[length]   = '.';
    filename[length+1] = 'b';
    filename[length+2] = 'g';
    filename[length+3] = '\0';
  }
  unsigned short *tmp_ptr;
  FILE *out_file;
  out_file = fopen(filename, "w");
  if(out_file == NULL){
    cout << "\n\t Cannot open file "<<filename<< " for writing \tReturning without writing\n";
    return;
  }    

  if(offset!=NULL || windowsize!=NULL){
    check_boundaries(offset, windowsize);
    for(int z=0; z<windowsize[2]; z++){
      for(int y=0; y<windowsize[1]; y++){
	tmp_ptr = &_AAI[offset[0]+(y+offset[1])*_dimx+(z+offset[2])*_dimx*_dimy];
	fwrite(tmp_ptr, sizeof(unsigned short), windowsize[0], out_file);
      }
    }
    for(int z=0; z<windowsize[2]; z++){
      for(int y=0; y<windowsize[1]; y++){
	tmp_ptr = &_ADC[offset[0]+(y+offset[1])*_dimx+(z+offset[2])*_dimx*_dimy];
	fwrite(tmp_ptr, sizeof(unsigned short), windowsize[0], out_file);
      }
    }
    for(int z=0; z<windowsize[2]; z++){
      for(int y=0; y<windowsize[1]; y++){
	tmp_ptr = &_DWI[offset[0]+(y+offset[1])*_dimx+(z+offset[2])*_dimx*_dimy];
	fwrite(tmp_ptr, sizeof(unsigned short), windowsize[0], out_file);
      }
    }
    for(int z=0; z<windowsize[2]; z++){
      for(int y=0; y<windowsize[1]; y++){
	tmp_ptr = &_RAI[offset[0]+(y+offset[1])*_dimx+(z+offset[2])*_dimx*_dimy];
	fwrite(tmp_ptr, sizeof(unsigned short), windowsize[0], out_file);
      }
    }
    for(int z=0; z<windowsize[2]; z++){
      for(int y=0; y<windowsize[1]; y++){
	tmp_ptr = &_T2W[offset[0]+(y+offset[1])*_dimx+(z+offset[2])*_dimx*_dimy];
	fwrite(tmp_ptr, sizeof(unsigned short), windowsize[0], out_file);
      }
    }
  }

  else {
    fwrite(_AAI, sizeof(unsigned short), _dimx*_dimy*_dimz, out_file);
    fwrite(_ADC, sizeof(unsigned short), _dimx*_dimy*_dimz, out_file);
    fwrite(_DWI, sizeof(unsigned short), _dimx*_dimy*_dimz, out_file);
    fwrite(_RAI, sizeof(unsigned short), _dimx*_dimy*_dimz, out_file);
    fwrite(_T2W, sizeof(unsigned short), _dimx*_dimy*_dimz, out_file);
  }
  fclose(out_file);
  delete[] filename;
  delete[] xmlpath;
  cout <<"\t\t\tdone\n"<<flush;
}

unsigned short 
background::getAAI(const int fpos)
{
  return _AAI[fpos];
}

unsigned short 
background::getADC(const int fpos)
{
  return _ADC[fpos];
}

unsigned short 
background::getDWI(const int fpos)
{
  return _DWI[fpos];
}
  
unsigned short 
background::getRAI(const int fpos)
{
  return _RAI[fpos];
}

unsigned short
background::getT2W(const int fpos)
{
  return _T2W[fpos];
}

unsigned short * 
background::getAAI()
{
  return _AAI;
}
  
unsigned short * 
background::getADC()
{
  return _ADC;
}
  
unsigned short * 
background::getDWI()
{
  return _DWI;
}
  
unsigned short * 
background::getRAI()
{
  return _RAI;
}
  
unsigned short * 
background::getT2W()
{
  return _T2W;
}

void 
background::setAAI(const unsigned short* AAI)
{
  for(int i=0; i<_dimx*_dimy*_dimz; i++){
    _AAI[i] = AAI[i];
  }
}

void 
background::setADC(const unsigned short* ADC)
{
  for(int i=0; i<_dimx*_dimy*_dimz; i++){
    _ADC[i] = ADC[i];
  }
}
  
void 
background::setDWI(const unsigned short* DWI)
{
  for(int i=0; i<_dimx*_dimy*_dimz; i++){
    _DWI[i] = DWI[i];
  }
}
  
void 
background::setRAI(const unsigned short* RAI)
{
  for(int i=0; i<_dimx*_dimy*_dimz; i++){
    _RAI[i] = RAI[i];
  }
}
  
void 
background::setT2W(const unsigned short* T2W)
{
  for(int i=0; i<_dimx*_dimy*_dimz; i++){
    _T2W[i] = T2W[i];
  }
}

void 
background::mask(const float thresh, const int selection, const int slice)
{
  int position=0;
  cout << "Masking all backgrounds; threshold "<<thresh<<"\t\t"<<flush;
  for(int i=0; i < _dimx*_dimy; i++){
    position = i + slice*_dimx*_dimy;
    switch(selection){
    case 0: if(_AAI[position]<thresh) this->setZero(position);    break;
    case 1: if(_ADC[position]<thresh) this->setZero(position);    break;
    case 2: if(_DWI[position]<thresh) this->setZero(position);    break;
    case 3: if(_RAI[position]<thresh) this->setZero(position);    break;
    case 4: if(_T2W[position]<thresh) this->setZero(position);    break;
    default: cout << "WARNING: Wrong selection in mask() in background.h\nReturning without masking\n"; break;
    }
  }
  this->reParameter(slice);
  cout << "done\n";
}

void 
background::masksingle(const float thresh, const int selection, const int slice)
{
  int position=0;
  cout << "Masking background "<<selection<<"; threshold "<<thresh<<"\t\t"<<flush;
  for(int i=0; i < _dimx*_dimy; i++){
    position = i + slice*_dimx*_dimy;
    switch(selection){
    case 0: if(_AAI[position]<thresh) _AAI[position] = 0;    break;
    case 1: if(_ADC[position]<thresh) _ADC[position] = 0;    break;
    case 2: if(_DWI[position]<thresh) _DWI[position] = 0;    break;
    case 3: if(_RAI[position]<thresh) _RAI[position] = 0;    break;
    case 4: if(_T2W[position]<thresh) _T2W[position] = 0;    break;
    default: cout << "WARNING: Wrong selection in mask() in background.h\nReturning without masking\n"; break;
    }
  }
  this->reParameter(slice);
  cout << "done\n";
}

void
background::remask(float * segmentation)
{
  for(int i=0; i < _dimx*_dimy*_dimz; i++){
    if(segmentation[i]==0){
      this->setZero(i);
    }
  }
  this->reParameter();
}

void 
background::setZero(const int fpos)
{
  _AAI[fpos] = 0;
  _ADC[fpos] = 0;
  _DWI[fpos] = 0;
  _RAI[fpos] = 0;
  _T2W[fpos] = 0;
}
 
void 
background::reParameter(const int slice=-1)
{
  unsigned short MAXSHORT = 65535;
  unsigned short aaimin = MAXSHORT;
  unsigned short adcmin = MAXSHORT;
  unsigned short dwimin = MAXSHORT;
  unsigned short raimin = MAXSHORT;
  unsigned short t2wmin = MAXSHORT;

  if(slice==-1){
    for(int i=0; i<_dimx*_dimy*_dimz; i++){
      if(_AAI[i]<aaimin && _AAI[i]!=0 ) aaimin = _AAI[i];
      if(_ADC[i]<adcmin && _ADC[i]!=0 ) adcmin = _ADC[i];
      if(_DWI[i]<dwimin && _DWI[i]!=0 ) dwimin = _DWI[i];
      if(_RAI[i]<raimin && _RAI[i]!=0 ) raimin = _RAI[i];
      if(_T2W[i]<t2wmin && _T2W[i]!=0 ) t2wmin = _T2W[i];
    }
    
    for(int i=0; i<_dimx*_dimy*_dimz; i++){ 
      if( _AAI[i]!=0 ) _AAI[i] = _AAI[i] - aaimin+1;
      if( _ADC[i]!=0 ) _ADC[i] = _ADC[i] - adcmin+1;
      if( _DWI[i]!=0 ) _DWI[i] = _DWI[i] - dwimin+1;
      if( _RAI[i]!=0 ) _RAI[i] = _RAI[i] - raimin+1;
      if( _T2W[i]!=0 ) _T2W[i] = _T2W[i] - t2wmin+1; 
    }
  } 
  else {
    for(int i=slice*_dimx*_dimz; i<(slice+1)*_dimx*_dimy; i++){
      if(_AAI[i]<aaimin && _AAI[i]!=0 ) aaimin = _AAI[i];
      if(_ADC[i]<adcmin && _ADC[i]!=0 ) adcmin = _ADC[i];
      if(_DWI[i]<dwimin && _DWI[i]!=0 ) dwimin = _DWI[i];
      if(_RAI[i]<raimin && _RAI[i]!=0 ) raimin = _RAI[i];
      if(_T2W[i]<t2wmin && _T2W[i]!=0 ) t2wmin = _T2W[i];
    }
    
    for(int i=slice*_dimx*_dimz; i<(slice+1)*_dimx*_dimy; i++){ 
      if( _AAI[i]!=0 ) _AAI[i] = _AAI[i] - aaimin+1;
      if( _ADC[i]!=0 ) _ADC[i] = _ADC[i] - adcmin+1;
      if( _DWI[i]!=0 ) _DWI[i] = _DWI[i] - dwimin+1;
      if( _RAI[i]!=0 ) _RAI[i] = _RAI[i] - raimin+1;
      if( _T2W[i]!=0 ) _T2W[i] = _T2W[i] - t2wmin+1; 
    }
  }
}
 
bool 
background::isZero(const int fpos, const int selection=-1)
{
  if(selection ==-1){
    if(_AAI[fpos] == 0 && 
       _ADC[fpos] == 0 &&
       _DWI[fpos] == 0 &&
       _RAI[fpos] == 0 &&
       _T2W[fpos] == 0) return true;
    return false;
  }
  else if(selection==0){  if(_AAI[fpos] == 0 ) return true;}
  else if(selection==1){  if(_ADC[fpos] == 0 ) return true;}
  else if(selection==2){  if(_DWI[fpos] == 0 ) return true;}
  else if(selection==3){  if(_RAI[fpos] == 0 ) return true;}
  else if(selection==4){  if(_T2W[fpos] == 0 ) return true;}
  return false;
}

void
background::saveAsMRI(const char *path, const int selection)
{
  char * location = new char[300];
  char * buffer   = new char[10];
  float *image    = new float[_dimx*_dimy];

  int jtmp = 0;
  for (int j=0; j<_dimz; j++){
    jtmp = j*_dimx*_dimy;
    for(int i=0; i<_dimx*_dimy; i++){
      switch(selection){
      case 0: image[i] = _AAI[i+jtmp];    break;
      case 1: image[i] = _ADC[i+jtmp];    break;
      case 2: image[i] = _DWI[i+jtmp];    break;
      case 3: image[i] = _RAI[i+jtmp];    break;
      case 4: image[i] = _T2W[i+jtmp];    break;
      default: cout << "WARNING: Wrong selection in saveAsMRI() in background.h\nReturning without saving\n"; return;
      }
    }
      
    switch(selection){
    case 0: strcpy(buffer, "/AAI.");    break;
    case 1: strcpy(buffer, "/ADC.");    break;
    case 2: strcpy(buffer, "/DWI.");    break;
    case 3: strcpy(buffer, "/RAI.");    break;
    case 4: strcpy(buffer, "/T2W.");    break;
    }

    strcpy(location, path);
    strcat(location, buffer);
    sprintf(buffer, "%03d", j+1);
    strcat(location, buffer);
    cout << "location to save MRI data is: "<< location<<endl;
    saveFloatAsMRI(location, image, _dimx*_dimy);
  } 

  delete[] image;
  delete[] location;
}