<commit_before>/*
 *  mcdespot_main.cpp
 *
 *  Created by <NAME> on 14/02/2012.
 *  Copyright 2012 <NAME>. All rights reserved.
 *
 */

#include <time.h>
#include <getopt.h>
#include <signal.h>
#include <iostream>
#include <atomic>
#include <Eigen/Dense>

#include "NiftiImage.h"
#include "DESPOT.h"
#include "DESPOT_Functors.h"
#include "RegionContraction.h"

#define USE_PROCPAR
#ifdef USE_PROCPAR
	#include "procpar.h"
	using namespace Recon;
#endif

using namespace std;
using namespace Eigen;

//******************************************************************************
// Arguments / Usage
//******************************************************************************
const string credit {
"mcdespot - written by <EMAIL>. \n\
Acknowledgements greatfully received, grant discussions welcome."
};

const string usage {
"Usage is: mcdespot [options]\n\
\n\
The input must consist of at least one line such as:\n\
SPGR path_to_spgr_file (path_to_B1_file) \n\
SSFP path_to_ssfp_file (path_to_B1) (path_to_B0) \n\
\n\
If -d or --drop is specified further input will be prompted for from stdin.\n\
\n\
Phase-cycling is specified in degrees. If no B0/B1 correction is specified,\n\
default values of B0 = 0 Hz and B1 = 1 will be used. B0 value is only used if\n\
--B0 option specified, otherwise it is a free parameter.\n\
\n\
Options:\n\
	--help, -h		  : Print this message.\n\
	--verbose, -v     : Print extra information.\n\
	--no-prompt, -p   : Don't print prompts for input.\n\
	--mask, -m file   : Mask input with specified file.\n\
	--out, -o path    : Add a prefix to the output filenames.\n\
	--1, --2, --3     : Use 1, 2 or 3 component model (default 2).\n\
	--M0 file         : M0 Map file.\n\
	--B0              : Use the default or specified B0 value (don't fit).\n\
	--start_slice n   : Only start processing at slice n.\n\
	--end_slice n     : Finish at slice n-1.\n\
	--drop, -d        : Drop certain flip-angles (Read from stdin).\n\
	--normalise, -n   : Normalise signals to maximum (Ignore M0).\n\
	--samples, -s n   : Use n samples for region contraction (Default 5000).\n\
	--retain, -r  n   : Retain n samples for new boundary (Default 50).\n\
	--contract, -c n  : Contract a maximum of n times (Default 10).\n\
	--expand, -e n    : Re-expand boundary by percentage n (Default 0).\n\
	--tesla, -t 3     : Boundaries suitable for 3T\n\
	            7     : Boundaries suitable for 7T (default)\n\
	            u     : User specified boundaries from stdin.\n"
};

static int verbose = false, prompt = true,
           normalise = false, drop = false, fitB0 = true,
           start_slice = -1, end_slice = -1, slice = 0,
		   samples = 5000, retain = 50, contract = 10,
		   components = 2, tesla = -1;
static double expand = 0.;
static string outPrefix;
static struct option long_options[] =
{
	{"M0", required_argument, 0, 'M'},
	{"mask", required_argument, 0, 'm'},
	{"out", required_argument, 0, 'o'},
	{"verbose", no_argument, 0, 'v'},
	{"no-prompt", no_argument, 0, 'p'},
	{"start_slice", required_argument, 0, 'S'},
	{"end_slice", required_argument, 0, 'E'},
	{"normalise", no_argument, &normalise, true},
	{"tesla", required_argument, 0, 't'},
	{"samples", required_argument, 0, 's'},
	{"retain", required_argument, 0, 'r'},
	{"contract", required_argument, 0, 'c'},
	{"expand", required_argument, 0, 'e'},
	{"drop", no_argument, &drop, true},
	{"help", no_argument, 0, 'h'},
	{"1", no_argument, &components, 1},
	{"2", no_argument, &components, 2},
	{"3", no_argument, &components, 3},
	{"B0", no_argument, &fitB0, false},
	{0, 0, 0, 0}
};
//******************************************************************************
#pragma mark SIGTERM interrupt handler - for ensuring data is saved on a ctrl-c
//******************************************************************************
NiftiImage *paramsHdrs, residualHdr;
vector<double *> paramsData;
double *residualData;

void int_handler(int sig);
void int_handler(int sig)
{
	fprintf(stdout, "Processing terminated. Writing currently processed data.\n");
	for (int p = 0; p < mcDESPOT::nP(components); p++) {
		paramsHdrs[p].writeSubvolume(0, 0, slice, 0, -1, -1, slice + 1, 1, paramsData[p]);
		paramsHdrs[p].close();
	}
	residualHdr.writeSubvolume(0, 0, slice, 0, -1, -1, slice + 1, 1, residualData);
	residualHdr.close();
	exit(EXIT_FAILURE);
}

//******************************************************************************
#pragma mark Read in all required files and data from cin
//******************************************************************************
NiftiImage *parseInput(vector<mcDESPOT::SignalType> &signalTypes, vector<VectorXd> &angles,
                       vector<DESPOTConstants> &consts,
				       vector<NiftiImage *> &signalFiles,
				       vector<NiftiImage *> &B1_files,
				       vector<NiftiImage *> &B0_files);
NiftiImage *parseInput(vector<mcDESPOT::SignalType> &signalTypes, vector<VectorXd> &angles,
                       vector<DESPOTConstants> &consts,
				       vector<NiftiImage *> &signalFiles,
				       vector<NiftiImage *> &B1_files,
				       vector<NiftiImage *> &B0_files) {
	NiftiImage *inHdr, *savedHeader;
	string type, path;
	if (prompt) cout << "Specify next image type (SPGR/SSFP): " << flush;
	while (getline(cin, type) && !type.empty()) {
		
		if (type == "SPGR")
			signalTypes.push_back(mcDESPOT::SignalSPGR);
		else if (type == "SSFP")
			signalTypes.push_back(mcDESPOT::SignalSSFP);
		else {
			cerr << "Unknown signal type: " << type << endl;
			exit(EXIT_FAILURE);
		}
		if (prompt) cout << "Enter image path: " << flush;
		getline(cin, path);
		inHdr = new NiftiImage(path, NiftiImage::NIFTI_READ);
		if (verbose) cout << "Opened " << type << " header: " << path << endl;
		if (signalFiles.size() == 0) // First time through save the header for comparison
			savedHeader = inHdr;
		inHdr->checkVoxelsCompatible(*savedHeader);
		signalFiles.push_back(inHdr);
		double inTR = 0., inPhase = 0.;
		VectorXd inAngles(inHdr->dim(4));
		#ifdef USE_PROCPAR
		ParameterList pars;
		if (ReadProcpar(inHdr->basename() + ".procpar", pars)) {
			inTR = RealValue(pars, "tr");
			for (int i = 0; i < inAngles.size(); i++)
				inAngles[i] = RealValue(pars, "flip1", i);
			if (signalTypes.back() == mcDESPOT::SignalSSFP)
				inPhase = RealValue(pars, "rfphase");
		} else
		#endif
		{
			if (prompt) cout << "Enter TR: " << flush;
			cin >> inTR;
			if (prompt) cout << "Enter Flip-angles: " << flush;
			for (int i = 0; i < inAngles[i]; i++)
				cin >> inAngles[i];
			getline(cin, path); // Just to eat the newline
			if (signalTypes.back() == mcDESPOT::SignalSSFP) {<commit_msg>Fixed silly mcDESPOT input bug when reading input angles.<commit_after>/*
 *  mcdespot_main.cpp
 *
 *  Created by <NAME> on 14/02/2012.
 *  Copyright 2012 <NAME>. All rights reserved.
 *
 */

#include <time.h>
#include <getopt.h>
#include <signal.h>
#include <iostream>
#include <atomic>
#include <Eigen/Dense>

#include "NiftiImage.h"
#include "DESPOT.h"
#include "DESPOT_Functors.h"
#include "RegionContraction.h"

#define USE_PROCPAR
#ifdef USE_PROCPAR
	#include "procpar.h"
	using namespace Recon;
#endif

using namespace std;
using namespace Eigen;

//******************************************************************************
// Arguments / Usage
//******************************************************************************
const string credit {
"mcdespot - written by <EMAIL>. \n\
Acknowledgements greatfully received, grant discussions welcome."
};

const string usage {
"Usage is: mcdespot [options]\n\
\n\
The input must consist of at least one line such as:\n\
SPGR path_to_spgr_file (path_to_B1_file) \n\
SSFP path_to_ssfp_file (path_to_B1) (path_to_B0) \n\
\n\
If -d or --drop is specified further input will be prompted for from stdin.\n\
\n\
Phase-cycling is specified in degrees. If no B0/B1 correction is specified,\n\
default values of B0 = 0 Hz and B1 = 1 will be used. B0 value is only used if\n\
--B0 option specified, otherwise it is a free parameter.\n\
\n\
Options:\n\
	--help, -h		  : Print this message.\n\
	--verbose, -v     : Print extra information.\n\
	--no-prompt, -p   : Don't print prompts for input.\n\
	--mask, -m file   : Mask input with specified file.\n\
	--out, -o path    : Add a prefix to the output filenames.\n\
	--1, --2, --3     : Use 1, 2 or 3 component model (default 2).\n\
	--M0 file         : M0 Map file.\n\
	--B0              : Use the default or specified B0 value (don't fit).\n\
	--start_slice n   : Only start processing at slice n.\n\
	--end_slice n     : Finish at slice n-1.\n\
	--drop, -d        : Drop certain flip-angles (Read from stdin).\n\
	--normalise, -n   : Normalise signals to maximum (Ignore M0).\n\
	--samples, -s n   : Use n samples for region contraction (Default 5000).\n\
	--retain, -r  n   : Retain n samples for new boundary (Default 50).\n\
	--contract, -c n  : Contract a maximum of n times (Default 10).\n\
	--expand, -e n    : Re-expand boundary by percentage n (Default 0).\n\
	--tesla, -t 3     : Boundaries suitable for 3T\n\
	            7     : Boundaries suitable for 7T (default)\n\
	            u     : User specified boundaries from stdin.\n"
};

static int verbose = false, prompt = true,
           normalise = false, drop = false, fitB0 = true,
           start_slice = -1, end_slice = -1, slice = 0,
		   samples = 5000, retain = 50, contract = 10,
		   components = 2, tesla = -1;
static double expand = 0.;
static string outPrefix;
static struct option long_options[] =
{
	{"M0", required_argument, 0, 'M'},
	{"mask", required_argument, 0, 'm'},
	{"out", required_argument, 0, 'o'},
	{"verbose", no_argument, 0, 'v'},
	{"no-prompt", no_argument, 0, 'p'},
	{"start_slice", required_argument, 0, 'S'},
	{"end_slice", required_argument, 0, 'E'},
	{"normalise", no_argument, &normalise, true},
	{"tesla", required_argument, 0, 't'},
	{"samples", required_argument, 0, 's'},
	{"retain", required_argument, 0, 'r'},
	{"contract", required_argument, 0, 'c'},
	{"expand", required_argument, 0, 'e'},
	{"drop", no_argument, &drop, true},
	{"help", no_argument, 0, 'h'},
	{"1", no_argument, &components, 1},
	{"2", no_argument, &components, 2},
	{"3", no_argument, &components, 3},
	{"B0", no_argument, &fitB0, false},
	{0, 0, 0, 0}
};
//******************************************************************************
#pragma mark SIGTERM interrupt handler - for ensuring data is saved on a ctrl-c
//******************************************************************************
NiftiImage *paramsHdrs, residualHdr;
vector<double *> paramsData;
double *residualData;

void int_handler(int sig);
void int_handler(int sig)
{
	fprintf(stdout, "Processing terminated. Writing currently processed data.\n");
	for (int p = 0; p < mcDESPOT::nP(components); p++) {
		paramsHdrs[p].writeSubvolume(0, 0, slice, 0, -1, -1, slice + 1, 1, paramsData[p]);
		paramsHdrs[p].close();
	}
	residualHdr.writeSubvolume(0, 0, slice, 0, -1, -1, slice + 1, 1, residualData);
	residualHdr.close();
	exit(EXIT_FAILURE);
}

//******************************************************************************
#pragma mark Read in all required files and data from cin
//******************************************************************************
NiftiImage *parseInput(vector<mcDESPOT::SignalType> &signalTypes, vector<VectorXd> &angles,
                       vector<DESPOTConstants> &consts,
				       vector<NiftiImage *> &signalFiles,
				       vector<NiftiImage *> &B1_files,
				       vector<NiftiImage *> &B0_files);
NiftiImage *parseInput(vector<mcDESPOT::SignalType> &signalTypes, vector<VectorXd> &angles,
                       vector<DESPOTConstants> &consts,
				       vector<NiftiImage *> &signalFiles,
				       vector<NiftiImage *> &B1_files,
				       vector<NiftiImage *> &B0_files) {
	NiftiImage *inHdr, *savedHeader;
	string type, path;
	if (prompt) cout << "Specify next image type (SPGR/SSFP): " << flush;
	while (getline(cin, type) && !type.empty()) {
		
		if (type == "SPGR")
			signalTypes.push_back(mcDESPOT::SignalSPGR);
		else if (type == "SSFP")
			signalTypes.push_back(mcDESPOT::SignalSSFP);
		else {
			cerr << "Unknown signal type: " << type << endl;
			exit(EXIT_FAILURE);
		}
		if (prompt) cout << "Enter image path: " << flush;
		getline(cin, path);
		inHdr = new NiftiImage(path, NiftiImage::NIFTI_READ);
		if (verbose) cout << "Opened " << type << " header: " << path << endl;
		if (signalFiles.size() == 0) // First time through save the header for comparison
			savedHeader = inHdr;
		inHdr->checkVoxelsCompatible(*savedHeader);
		signalFiles.push_back(inHdr);
		double inTR = 0., inPhase = 0.;
		VectorXd inAngles(inHdr->dim(4));
		#ifdef USE_PROCPAR
		ParameterList pars;
		if (ReadProcpar(inHdr->basename() + ".procpar", pars)) {
			inTR = RealValue(pars, "tr");
			for (int i = 0; i < inAngles.size(); i++)
				inAngles[i] = RealValue(pars, "flip1", i);
			if (signalTypes.back() == mcDESPOT::SignalSSFP)
				inPhase = RealValue(pars, "rfphase");
		} else
		#endif
		{
			if (prompt) cout << "Enter TR: " << flush;
			cin >> inTR;
			if (prompt) cout << "Enter Flip-angles: " << flush;
			for (int i = 0; i < inAngles.size(); i++)
				cin >> inAngles[i];
			getline(cin, path); // Just to eat the newline
			if (signalTypes.back() == mcDESPOT::SignalSSFP) {
{'fraction_non_alphanumeric': 0.21510228640192539, 'fraction_numerical': 0.017448856799037304, 'mean_word_length': 3.660708026638626, 'pattern_counts': {'":': 0, '<': 117, '<?xml version=': 0, '>': 61, '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': '3466554', 'n_tokens_mistral': 4867, 'n_tokens_neox': 4426, 'n_words': 1336}