<commit_before>#include <irtkGraphCutSegmentation_4D.h>
#include <irtkImage.h>
#include <string>
#include <irtkHistogram_1D.h>
#include <irtkDilation.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <irtkSubcorticalSegmentation_4D.h>
#include <irtkGradientImage.h>

irtkSubcorticalSegmentation_4D::irtkSubcorticalSegmentation_4D(){
	useMask = false;
}

irtkSubcorticalSegmentation_4D::irtkSubcorticalSegmentation_4D(string name, char * structureFile){
	readStructureNames(structureFile);
	imageName = name;
 	useMask = false;
}

void irtkSubcorticalSegmentation_4D::readParameters(char * filename){
	string line;
  ifstream from(filename);

  if (!from) {
    cerr << "Can't open file " << filename
    << endl;
    exit(1);
  }

  if (from.is_open())
  {
    while (! from.eof() )
    {
      getline (from,line);
      this->Read(line);

    }
    from.close();
  }
}

void irtkSubcorticalSegmentation_4D::readStructureNames(char * filename){
  string line;
  ifstream from(filename);
  if (!from) {
    cerr << "Can't open file " << filename
    << endl;
    exit(1);
  }
  if (from.is_open())
  {
    while (! from.eof() )
    {
    	getline (from,line);
       numTissues++;
    }
    from.close();
  }
  numTissues /= 2;
  names = new string*[numTissues];
  structValues = new int*[numTissues];
  ifstream 	from2(filename);  
  cout << "Read " << numTissues << " Structures:" << endl;
  int counter1 = 0;
  int counter2 = 0;
  int pos = 0;
  if (from2.is_open())
  {
    while (! from2.eof() )
    {
    	structValues[counter2] = new int;
    	names[counter1] = new string;
        getline (from2,line);
       
   		if(pos % 2 == 0){
        	*names[counter1++] = line;
        	pos ++;
        	 cout << *names[counter1-1] << ", ";
   		}
        else{
        	*structValues[counter2++] = atoi(line.c_str());
        	pos++;
        	 cout << *structValues[counter2-1] << " " << 	endl;
        }
    }
    from2.close();
  }

}

void irtkSubcorticalSegmentation_4D::generateGaussians(double minPerc){

    int _padding = -1;
	vector<double> mi_num(numTissues);
	vector<double> denom(numTissues);
	vector<double> sigma_num(numTissues);
	vector<double> mi(numTissues);
	vector<double> structureValues[numTissues];
	irtkRealPixel *ptr = input.GetPointerToVoxels();
	int min[numTissues];
	int max[numTissues];
	irtkRealPixel maxPerc[numTissues];
	for(int k = 0; k < numTissues; k++){
		min[k]=32767;
	  	max[k]=-32768;
	  	irtkRealPixel tmp;
	  	atlasI[k]->GetMinMax(&tmp,&maxPerc[k]);
	  	cout << "Structure " << k << ": Max atlas value: " << maxPerc[k] << "." << endl;
	}
	irtkRealPixel ** tissuePtr = new irtkRealPixel*[3];
	for(int i = 0; i < 3; i++){
		tissuePtr[i] = new irtkRealPixel;
		tissuePtr[i] = tissuePriors[i]->GetPointerToVoxels();
	}
	for(int i = 0; i < 3 ; i++){
		counter[i]=tissueMy[i]=tissueS[i]=0;
	}
	int atlasVoxelNum = atlasI[0]->GetNumberOfVoxels();
    int inputVoxelNum =  input.GetNumberOfVoxels();

	int whiteVoxelNum = tissuePriors[0]->GetNumberOfVoxels();
	if(whiteVoxelNum != atlasVoxelNum || atlasVoxelNum != inputVoxelNum   || inputVoxelNum != atlasVoxelNum){
		cerr << "Different Voxel Numbers!" << endl;
		exit(1);
 
	} 
	irtkRealPixel **atlasPtr;
	atlasPtr = new irtkRealPixel*[numTissues];
	for(int i = 0; i < numTissues; i++){
		atlasPtr[i] = new irtkRealPixel;
		atlasPtr[i] = atlasI[i]->GetPointerToVoxels();
	}
	for (int i = 0; i < input.GetNumberOfVoxels(); i++){

			for(int ii = 0; ii < 3; ii++){
				counter[ii] += (double)*tissuePtr[ii];
				tissueMy[ii] += *ptr * (double)*tissuePtr[ii];
				tissuePtr[ii]++;
			}
	  		for (int k = 0; k < numTissues; k++){
				if(*atlasPtr[k] >= maxPerc[k]*minPerc && *ptr > _padding){
					structureValues[k].push_back(*ptr);
					if(*ptr>max[k])
						max[k] = *ptr;
					if(*ptr<min[k])
						min[k] = *ptr;
				}
				atlasPtr[k]++;
	  		}
	    ptr++;
	}

	for(int i = 0; i < 3; i++){
		tissueMy[i] /= counter[i];
	}

	ptr = input.GetPointerToVoxels();
	for(int ii = 0; ii < 3; ii++){
		tissuePtr[ii] = tissuePriors[ii]->GetPointerToVoxels();
	}
	for (int i = 0; i < input.GetNumberOfVoxels(); i++){
		for(int ii = 0; ii < 3; ii++){
			if(*tissuePtr[ii] > 0){
				tissueS[ii] += (tissueMy[ii]-*ptr)*(tissueMy[ii]-*ptr);
			}
			tissuePtr[ii]++;
		}
	ptr++;
	}
	cout << endl << "Tissue Classes: " << endl; 
	for(int i = 0; i < 3 ; i++){
		cout   << i+1 << ": my: " << tissueMy[i] << " sigma: " << tissueS[i] /counter[i]<< endl;
	}



	for (int k = 0; k < numTissues; k++){
		for(unsigned int i = 0; i < structureValues[k].size(); i++){
			double val = structureValues[k][i];
				mi_num[k] += val;
				denom[k]  += 1;
		}
	    if (denom[k] != 0){
	      	mi[k] = mi_num[k] / denom[k];
    }
    else{
      	cerr << "Division by zero while computing tissue mean!" << endl;
      	exit(1);
    	}
  	}
	    for (int k = 0; k <numTissues; k++){
	    	for(unsigned int i = 0; i < structureValues[k].size(); i++){
	    		double val = structureValues[k][i];
	    			sigma_num[k] += ((val - mi[k]) * (val - mi[k]));
	    	}

	   }
 
    cout << endl << "Structures: " << endl ;
	for (int k = 0; k <numTissues; k++){
		_my[k] = mi[k];
		_sigma[k] = sigma_num[k] / denom[k];
		cout << *names[k] << ": my: " << k << " " << _my[k] << " sigma: " << _sigma[k] << "...no of voxels considered: " << denom[k]<< endl;
		cout << *names[k] << ": my: " << k << " " << _my[k] << " sigma: " << _sigma[k] << "...no of voxels considered: " << denom[k]<< endl;
  	}
  	cout << endl;
}

void irtkSubcorticalSegmentation_4D::doSegmentation(double l, double g, double c, double sigmaG){
	cout << "Do segmentation with c = " << c << ", lambda = " << l << ", alpha = " << g << ", sigmaG = " << sigmaG << endl;
	cout << "Get Gradient Magnitude " << endl;
	irtkGradientImageX<irtkRealPixel> gradientX;
	irtkGradientImageY<irtkRealPixel> gradientY;
	irtkGradientImageZ<irtkRealPixel> gradientZ;
	int _xDim = input.GetX();
	int _yDim = input.GetY();
	int _zDim = input.GetZ();
	irtkGenericImage<irtkRealPixel> *_xGrad;
	irtkGenericImage<irtkRealPixel> *_yGrad;
	irtkGenericImage<irtkRealPixel> *_zGrad;
	int _tDim = input.GetT();
 	_xGrad = new irtkGenericImage<irtkRealPixel>(_xDim, _yDim, _zDim, _tDim);
	_yGrad = new irtkGenericImage<irtkRealPixel>(_xDim, _yDim, _zDim, _tDim);
	_zGrad = new irtkGenericImage<irtkRealPixel>(_xDim, _yDim, _zDim, _tDim);
	gradientX.SetInput(&input);
	gradientY.SetInput(&input);
	gradientZ.SetInput(&input);
	gradientX.SetOutput(_xGrad);
	gradientY.SetOutput(_yGrad);
	gradientZ.SetOutput(_zGrad);
	gradientX.Run();
	gradientY.Run();
	gradientZ.Run();
	irtkGenericImage<irtkRealPixel> gradMagnitude;
	gradMagnitude.Initialize(input.GetImageAttributes());<commit_msg>Avoiding run-time initialization of arrays<commit_after>#include <irtkGraphCutSegmentation_4D.h>
#include <irtkImage.h>
#include <string>
#include <irtkHistogram_1D.h>
#include <irtkDilation.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <irtkSubcorticalSegmentation_4D.h>
#include <irtkGradientImage.h>

irtkSubcorticalSegmentation_4D::irtkSubcorticalSegmentation_4D(){
	useMask = false;
}

irtkSubcorticalSegmentation_4D::irtkSubcorticalSegmentation_4D(string name, char * structureFile){
	readStructureNames(structureFile);
	imageName = name;
 	useMask = false;
}

void irtkSubcorticalSegmentation_4D::readParameters(char * filename){
	string line;
  ifstream from(filename);

  if (!from) {
    cerr << "Can't open file " << filename
    << endl;
    exit(1);
  }

  if (from.is_open())
  {
    while (! from.eof() )
    {
      getline (from,line);
      this->Read(line);

    }
    from.close();
  }
}

void irtkSubcorticalSegmentation_4D::readStructureNames(char * filename){
  string line;
  ifstream from(filename);
  if (!from) {
    cerr << "Can't open file " << filename
    << endl;
    exit(1);
  }
  if (from.is_open())
  {
    while (! from.eof() )
    {
    	getline (from,line);
       numTissues++;
    }
    from.close();
  }
  numTissues /= 2;
  names = new string*[numTissues];
  structValues = new int*[numTissues];
  ifstream 	from2(filename);  
  cout << "Read " << numTissues << " Structures:" << endl;
  int counter1 = 0;
  int counter2 = 0;
  int pos = 0;
  if (from2.is_open())
  {
    while (! from2.eof() )
    {
    	structValues[counter2] = new int;
    	names[counter1] = new string;
        getline (from2,line);
       
   		if(pos % 2 == 0){
        	*names[counter1++] = line;
        	pos ++;
        	 cout << *names[counter1-1] << ", ";
   		}
        else{
        	*structValues[counter2++] = atoi(line.c_str());
        	pos++;
        	 cout << *structValues[counter2-1] << " " << 	endl;
        }
    }
    from2.close();
  }

}

void irtkSubcorticalSegmentation_4D::generateGaussians(double minPerc){

    int _padding = -1;
	vector<double> mi_num(numTissues);
	vector<double> denom(numTissues);
	vector<double> sigma_num(numTissues);
	vector<double> mi(numTissues);
	vector<double> *structureValues = new vector<double>[numTissues];
	irtkRealPixel *ptr = input.GetPointerToVoxels();
	int min[numTissues];
	int max[numTissues];
	irtkRealPixel *maxPerc = new irtkRealPixel[numTissues];
	for(int k = 0; k < numTissues; k++){
		min[k]=32767;
	  	max[k]=-32768;
	  	irtkRealPixel tmp;
	  	atlasI[k]->GetMinMax(&tmp,&maxPerc[k]);
	  	cout << "Structure " << k << ": Max atlas value: " << maxPerc[k] << "." << endl;
	}
	irtkRealPixel ** tissuePtr = new irtkRealPixel*[3];
	for(int i = 0; i < 3; i++){
		tissuePtr[i] = new irtkRealPixel;
		tissuePtr[i] = tissuePriors[i]->GetPointerToVoxels();
	}
	for(int i = 0; i < 3 ; i++){
		counter[i]=tissueMy[i]=tissueS[i]=0;
	}
	int atlasVoxelNum = atlasI[0]->GetNumberOfVoxels();
    int inputVoxelNum =  input.GetNumberOfVoxels();

	int whiteVoxelNum = tissuePriors[0]->GetNumberOfVoxels();
	if(whiteVoxelNum != atlasVoxelNum || atlasVoxelNum != inputVoxelNum   || inputVoxelNum != atlasVoxelNum){
		cerr << "Different Voxel Numbers!" << endl;
		exit(1);
 
	} 
	irtkRealPixel **atlasPtr;
	atlasPtr = new irtkRealPixel*[numTissues];
	for(int i = 0; i < numTissues; i++){
		atlasPtr[i] = new irtkRealPixel;
		atlasPtr[i] = atlasI[i]->GetPointerToVoxels();
	}
	for (int i = 0; i < input.GetNumberOfVoxels(); i++){

			for(int ii = 0; ii < 3; ii++){
				counter[ii] += (double)*tissuePtr[ii];
				tissueMy[ii] += *ptr * (double)*tissuePtr[ii];
				tissuePtr[ii]++;
			}
	  		for (int k = 0; k < numTissues; k++){
				if(*atlasPtr[k] >= maxPerc[k]*minPerc && *ptr > _padding){
					structureValues[k].push_back(*ptr);
					if(*ptr>max[k])
						max[k] = *ptr;
					if(*ptr<min[k])
						min[k] = *ptr;
				}
				atlasPtr[k]++;
	  		}
	    ptr++;
	}

	for(int i = 0; i < 3; i++){
		tissueMy[i] /= counter[i];
	}

	ptr = input.GetPointerToVoxels();
	for(int ii = 0; ii < 3; ii++){
		tissuePtr[ii] = tissuePriors[ii]->GetPointerToVoxels();
	}
	for (int i = 0; i < input.GetNumberOfVoxels(); i++){
		for(int ii = 0; ii < 3; ii++){
			if(*tissuePtr[ii] > 0){
				tissueS[ii] += (tissueMy[ii]-*ptr)*(tissueMy[ii]-*ptr);
			}
			tissuePtr[ii]++;
		}
	ptr++;
	}
	cout << endl << "Tissue Classes: " << endl; 
	for(int i = 0; i < 3 ; i++){
		cout   << i+1 << ": my: " << tissueMy[i] << " sigma: " << tissueS[i] /counter[i]<< endl;
	}



	for (int k = 0; k < numTissues; k++){
		for(unsigned int i = 0; i < structureValues[k].size(); i++){
			double val = structureValues[k][i];
				mi_num[k] += val;
				denom[k]  += 1;
		}
	    if (denom[k] != 0){
	      	mi[k] = mi_num[k] / denom[k];
    }
    else{
      	cerr << "Division by zero while computing tissue mean!" << endl;
      	exit(1);
    	}
  	}
	    for (int k = 0; k <numTissues; k++){
	    	for(unsigned int i = 0; i < structureValues[k].size(); i++){
	    		double val = structureValues[k][i];
	    			sigma_num[k] += ((val - mi[k]) * (val - mi[k]));
	    	}

	   }
 
    cout << endl << "Structures: " << endl ;
	for (int k = 0; k <numTissues; k++){
		_my[k] = mi[k];
		_sigma[k] = sigma_num[k] / denom[k];
		cout << *names[k] << ": my: " << k << " " << _my[k] << " sigma: " << _sigma[k] << "...no of voxels considered: " << denom[k]<< endl;
		cout << *names[k] << ": my: " << k << " " << _my[k] << " sigma: " << _sigma[k] << "...no of voxels considered: " << denom[k]<< endl;
  	}
  	cout << endl;

  	delete []structureValues;
  	delete []maxPerc;
}

void irtkSubcorticalSegmentation_4D::doSegmentation(double l, double g, double c, double sigmaG){
	cout << "Do segmentation with c = " << c << ", lambda = " << l << ", alpha = " << g << ", sigmaG = " << sigmaG << endl;
	cout << "Get Gradient Magnitude " << endl;
	irtkGradientImageX<irtkRealPixel> gradientX;
	irtkGradientImageY<irtkRealPixel> gradientY;
	irtkGradientImageZ<irtkRealPixel> gradientZ;
	int _xDim = input.GetX();
	int _yDim = input.GetY();
	int _zDim = input.GetZ();
	irtkGenericImage<irtkRealPixel> *_xGrad;
	irtkGenericImage<irtkRealPixel> *_yGrad;
	irtkGenericImage<irtkRealPixel> *_zGrad;
	int _tDim = input.GetT();
 	_xGrad = new irtkGenericImage<irtkRealPixel>(_xDim, _yDim, _zDim, _tDim);
	_yGrad = new irtkGenericImage<irtkRealPixel>(_xDim, _yDim, _zDim, _tDim);
	_zGrad = new irtkGenericImage<irtkRealPixel>(_xDim, _yDim, _zDim, _tDim);
	gradientX.SetInput(&input);
	gradientY.SetInput(&input);
	gradientZ.SetInput(&input);
	gradientX.SetOutput(_xGrad);
	gradientY.SetOutput(_yGrad);
	gradientZ.SetOutput(_zGrad);
	gradientX.Run();
	gradientY.Run();
	gradientZ.Run();
	irtkGenericImage<irtkRealPixel> gradMagnitude;
	gradMagnitude.Initialize(input.GetImageAttributes());
{'fraction_non_alphanumeric': 0.18067819247607034, 'fraction_numerical': 0.011723677376270683, 'mean_word_length': 3.595294919877259, 'pattern_counts': {'":': 8, '<': 366, '<?xml version=': 0, '>': 76, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3225355', 'n_tokens_mistral': 5726, 'n_tokens_neox': 5300, 'n_words': 1041}