diff --git a/merge_parts_gpu.cc b/merge_parts_gpu.cc
new file mode 100644
index 0000000000000000000000000000000000000000..35eaa002300508121f52768bac37b59537565891
--- /dev/null
+++ b/merge_parts_gpu.cc
@@ -0,0 +1,141 @@
+/*  merge_parts_gpu.cc
+
+    Tim Behrens, Saad Jbabdi, Stam Sotiropoulos, Moises Hernandez  - FMRIB Image Analysis Group
+
+    Copyright (C) 2005 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#include "xfibresoptions.h"
+#include "newmat.h"
+#include "newimage/newimageall.h"
+
+using namespace Xfibres;
+
+
+void join_Parts(NEWIMAGE::volume<float> mask, string name_in, string name_out, string subjdir, int nvox, int nsamples, int nParts, float max, float min){
+		int size_parts = nvox/nParts;
+		int last_part = nvox - ((nParts-1)*size_parts);
+
+		//if(mean) nsamples=1;
+
+		Matrix result(nsamples,0);
+		Matrix part;
+
+		for(int i=0;i<(nParts-1);i++){
+			part.ReSize(nsamples,size_parts);
+			std::ostringstream num;
+			num << i;
+			std::string part_number;
+			part_number.assign(num.str());
+			std::string aux;
+			while(part_number.size()<4){
+				aux = "0" + part_number;
+				part_number.assign(aux);
+			}
+			std::string file_name;
+			file_name.assign(subjdir);
+			file_name += ".bedpostX/diff_parts/data_part_";
+			file_name += part_number; 
+			file_name += "/"; 
+			file_name += name_in; 
+			file_name += "J"; 
+
+			ifstream in;
+			in.open (file_name.data(), ios::in | ios::binary);
+			in.read((char*)&part(1,1), size_parts*nsamples*sizeof(Real));
+			in.close();
+
+			result = result | part;
+		}
+		part.ReSize(nsamples,last_part);
+		std::ostringstream num;
+		num << nParts-1;;
+		std::string part_number;
+		part_number.assign(num.str());
+		std::string aux;
+		while(part_number.size()<4){
+			aux = "0" + part_number;
+			part_number.assign(aux);
+		}
+		std::string file_name;
+		file_name.assign(subjdir);
+		file_name += ".bedpostX/diff_parts/data_part_";
+		file_name += part_number; 
+		file_name += "/"; 
+		file_name += name_in;
+		file_name += "J";  
+
+		ifstream in;
+		in.open (file_name.data(), ios::in | ios::binary);
+		in.read((char*)&part(1,1), last_part*nsamples*sizeof(Real));
+		in.close();
+		result = result | part;
+
+		NEWIMAGE::volume4D<float> tmp;
+      		tmp.setmatrix(result,mask);
+		if(max==-10) max=tmp.max();
+		if(min==-10) min=tmp.min(); 
+		tmp.setDisplayMaximumMinimum(max,min);
+     		save_volume4D(tmp,subjdir+".bedpostX/"+name_out);
+}
+
+//////////////////////////////////////////////////////////
+//       MERGE THE OUTPUTS FILES OF BEDPOSTX
+//////////////////////////////////////////////////////////
+//parameters:
+// argc - 3 nvox
+// argc - 2 nParts
+// argc - 1 subjdir
+
+int main(int argc, char *argv[])
+{
+	// Setup logging:
+    	Log& logger = LogSingleton::getInstance();
+    	xfibresOptions& opts = xfibresOptions::getInstance();
+	opts.parse_command_line(argc-3,argv,logger);
+
+    	NEWIMAGE::volume<float> mask;
+    	read_volume(mask,opts.maskfile.value());
+
+	cout << opts.maskfile.value() << endl;
+    	
+	int nvox = atoi(argv[argc-3]);
+	int nParts = atoi(argv[argc-2]);
+	string subjdir = argv[argc-1];
+
+	int nsamples = opts.njumps.value()/opts.sampleevery.value();
+	
+	//////////////////////////////////////////////////////////////
+	////////// JOIN Results of the Parts //////////////////////
+	//////////////////////////////////////////////////////////////
+
+	if(opts.modelnum.value()==1){
+		join_Parts(mask,"mean_dsamples","mean_dsamples",subjdir, nvox, 1, nParts, -10, 0);
+	}else if(opts.modelnum.value()==2){
+		join_Parts(mask,"mean_dsamples","mean_dsamples",subjdir, nvox, 1, nParts, -10, 0);
+		join_Parts(mask,"mean_d_stdsamples","mean_d_stdsamples",subjdir, nvox, 1, nParts, -10, 0);
+		//join_Parts(mask,"dsamples","dsamples",subjdir, nvox, nsamples, nParts, -10, 0);
+		//join_Parts(mask,"d_stdsamples","d_stdsamples",subjdir, nvox, nsamples, nParts, -10, 0);
+	}	
+	if (opts.f0.value()){
+      		join_Parts(mask,"mean_f0samples","mean_f0samples",subjdir, nvox, 1, nParts, 1, 0);
+		//join_Parts(mask,"f0samples","f0samples",subjdir, nvox, nsamples, nParts, 1, 0);
+    	}
+    	if (opts.rician.value()){
+		join_Parts(mask,"mean_tausamples","mean_tausamples",subjdir, nvox, 1, nParts, -10, 0);
+    	}
+	join_Parts(mask,"mean_S0samples","mean_S0samples",subjdir, nvox, 1, nParts, -10, 0);
+
+	for(int f=0;f<opts.nfibres.value();f++){
+		join_Parts(mask,"th"+num2str(f+1)+"samples","merged_th"+num2str(f+1)+"samples",subjdir, nvox, nsamples, nParts, -10, -10);
+		join_Parts(mask,"ph"+num2str(f+1)+"samples","merged_ph"+num2str(f+1)+"samples",subjdir, nvox, nsamples, nParts, -10, -10);
+		join_Parts(mask,"f"+num2str(f+1)+"samples","merged_f"+num2str(f+1)+"samples",subjdir, nvox, nsamples, nParts, 1, 0);
+
+		//join_Parts(mask,"mean_f"+num2str(f+1)+"samples",subjdir, nvox, 1, nParts, 1, 0);
+		//join_Parts(mask,"dyads"+num2str(f+1),subjdir, nvox, nsamples, nParts, 1, -1);
+	}	
+		
+  	return 0;
+}
+