*** empty log message ***

slicetimer.cc

@@ -55,7 +55,7 @@ Option<string> tcustom(string("--tcustom"), string(""),
 			  string("filename of single-column custom interleave timing file"),
 			  false, requires_argument);
 Option<string> ocustom(string("--ocustom"), string(""),
-			  string("filename of single-column custom interleave order file"),
+			  string("filename of single-column custom interleave order file (first slice is referred to as 1 not 0)"),
 			  false, requires_argument);
 Option<int> direction(string("-d,--direction"), 3,
 		      string("direction of slice acquisition (x=1,y=2,z=3) - default is z"),
@@ -68,11 +68,10 @@ int do_slice_correction()
 {
   volume4D<float> timeseries;
   volumeinfo in_info;
-  ColumnVector voxeltimeseries, interpseries;
   Matrix timings;
   
   int no_volumes, no_slices;
-  float repeat_time, offset, slice_spacing, recenter;
+  float repeat_time, offset=0, slice_spacing;
 
   if (inputname.set()) {
     if (verbose.value()) { cout << "Reading input volume" << endl; }
@@ -106,16 +105,13 @@ int do_slice_correction()
     timings = read_ascii_matrix(ocustom.value(), timeseries.zsize(), 1);
   }
 
-  offset = 0.0;
-
   ColumnVector userkernel = sinckernel1D("hanning", 7, 1201);
-  
-  //  if (down.value())
-  //  recenter = (((float) no_slices)/2 + 0.5)/ no_slices;
-  //else
-  recenter = (((float) no_slices)/2 - 0.5)/ no_slices;
+  // for(int i=1; i<=1201; i++) cout << i << " " << userkernel(i) << endl;
+
+  float recenter = (((float) no_slices)/2 - 0.5)/ no_slices;
   
   for (int slice=1; slice<=no_slices; slice++) {
+
     if (tcustom.set()) {
       offset = -timings(slice, 1);
     } else if (ocustom.set()) { 
@@ -140,15 +136,17 @@ int do_slice_correction()
     
     for (int x_pos = 0; x_pos < timeseries. xsize(); x_pos++)
       for (int y_pos = 0; y_pos < timeseries. ysize(); y_pos++){
-	voxeltimeseries = timeseries.voxelts(x_pos,y_pos,slice-1);
-	interpseries = voxeltimeseries;
-	for (int time_step=1; time_step <= no_volumes; time_step++)
-	  interpseries(time_step) = interpolate_1d(voxeltimeseries, time_step + offset + recenter);
+	ColumnVector voxeltimeseries = timeseries.voxelts(x_pos,y_pos,slice-1);
+	ColumnVector interpseries = voxeltimeseries;
+	for (int time_step=1; time_step <= no_volumes; time_step++){
+	  // interpseries(time_step) = interpolate_1d(voxeltimeseries, time_step + offset + recenter);
+	  interpseries(time_step) = kernelinterpolation_1d(voxeltimeseries, time_step + offset + recenter, userkernel, 3);
+	}
 	timeseries.setvoxelts(interpseries,x_pos,y_pos,slice-1);
       }
     
     if (verbose.value())
-      cerr << "Slice / Offset (s) " << slice << " / " << offset *repeat_time * -1 << endl;
+      cerr << "Slice " << slice << " offset " << offset + recenter << endl;
   }
   
   if (direction.value() == 1) timeseries. swapdimensions(3,2,1); // reverse Flip z and x