diff --git a/zUtil_Cxx/internal_gtAssignGrainToVol.cpp b/zUtil_Cxx/internal_gtAssignGrainToVol.cpp
index 2cc4c532fc5bad79e06d8eebc8cd791537dc04b1..126fa44f82ac7972e9f21a9b86b5bda8ea2afc32 100755
--- a/zUtil_Cxx/internal_gtAssignGrainToVol.cpp
+++ b/zUtil_Cxx/internal_gtAssignGrainToVol.cpp
@@ -1,7 +1,7 @@
 /**
  * Adds a grain to the volume, without taking care of the superposition between
  * grains.
- * If two grains overlap, the region becomes the sum of their indexes
+ * If two grains overlap, the region is overwritten by the new volume
  *
  * Nicola Vigano', 2012, ID11 @ ESRF vigano@esrf.eu
  */
@@ -29,6 +29,36 @@ public:
   { }
 };
 
+template<typename TypeOut, typename TypeIn>
+class GtMatlabVolumeAssignValue : public GtMatlabVolumeLooper<TypeOut, TypeIn> {
+protected:
+  void doWork(TypeOut & dest, const TypeIn & orig)
+  {
+    if (orig) { dest = (TypeOut) orig; }
+  }
+public:
+  GtMatlabVolumeAssignValue(void * _out, void * _vol, const double * limits,
+      const mwSize * dims_out, const mwSize * dims_vol)
+  : GtMatlabVolumeLooper<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol, 0)
+  { }
+};
+
+template<typename TypeOut, typename TypeIn>
+inline void
+assign(void * _out, void * _vol, const double * limits,
+    const mwSize * dims_out, const mwSize * dims_vol, const double & indx)
+{
+  GtMatlabVolumeLooper<TypeOut, TypeIn> * looper = NULL;
+  if (indx) {
+    looper = new GtMatlabVolumeAssign<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol, indx);
+  } else {
+    looper = new GtMatlabVolumeAssignValue<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol);
+  }
+  looper->loopOver3D();
+  delete looper;
+}
+
+
 void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
 {
   /* Incoming image dimensions */
@@ -100,63 +130,43 @@ extDoAssignment(mxArray * output, void * vol, const double * limits,
 {
   switch(mxGetClassID(output)) {
     case mxDOUBLE_CLASS: {
-      GtMatlabVolumeAssign<double, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<double, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxSINGLE_CLASS: {
-      GtMatlabVolumeAssign<float, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<float, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT8_CLASS: {
-      GtMatlabVolumeAssign<char, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<char, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT8_CLASS: {
-      GtMatlabVolumeAssign<unsigned char, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned char, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT16_CLASS: {
-      GtMatlabVolumeAssign<short int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<short int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT16_CLASS: {
-      GtMatlabVolumeAssign<unsigned short int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned short int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT32_CLASS: {
-      GtMatlabVolumeAssign<int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT32_CLASS: {
-      GtMatlabVolumeAssign<unsigned int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT64_CLASS: {
-      GtMatlabVolumeAssign<long int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<long int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT64_CLASS: {
-      GtMatlabVolumeAssign<unsigned long int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned long int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     default:
diff --git a/zUtil_Cxx/internal_gtAssignGrainToVol_interf.cpp b/zUtil_Cxx/internal_gtAssignGrainToVol_interf.cpp
index a435f2d798f660af6659b3bd3cc26ae530f47d1c..69763c78a578ba25de7826b0dfa2533df0bf6f15 100755
--- a/zUtil_Cxx/internal_gtAssignGrainToVol_interf.cpp
+++ b/zUtil_Cxx/internal_gtAssignGrainToVol_interf.cpp
@@ -29,6 +29,35 @@ public:
   { }
 };
 
+template<typename TypeOut, typename TypeIn>
+class GtMatlabVolumeInterferenceValue : public GtMatlabVolumeLooper<TypeOut, TypeIn> {
+protected:
+  void doWork(TypeOut & dest, const TypeIn & orig)
+  {
+    if (orig) { dest = (TypeOut) (-1 + (! dest) * orig); }
+  }
+public:
+  GtMatlabVolumeInterferenceValue(void * _out, void * _vol, const double * limits,
+      const mwSize * dims_out, const mwSize * dims_vol)
+  : GtMatlabVolumeLooper<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol, 0)
+  { }
+};
+
+template<typename TypeOut, typename TypeIn>
+inline void
+assign(void * _out, void * _vol, const double * limits,
+    const mwSize * dims_out, const mwSize * dims_vol, const double & indx)
+{
+  GtMatlabVolumeLooper<TypeOut, TypeIn> * looper = NULL;
+  if (indx) {
+    looper = new GtMatlabVolumeInterference<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol, indx);
+  } else {
+    looper = new GtMatlabVolumeInterferenceValue<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol);
+  }
+  looper->loopOver3D();
+  delete looper;
+}
+
 void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
 {
   /* Incoming image dimensions */
@@ -102,67 +131,46 @@ extDoAssignment(mxArray * output, void * vol, const double * limits,
 {
   switch(mxGetClassID(output)) {
     case mxDOUBLE_CLASS: {
-      GtMatlabVolumeInterference<double, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<double, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxSINGLE_CLASS: {
-      GtMatlabVolumeInterference<float, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<float, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT8_CLASS: {
-      GtMatlabVolumeInterference<char, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<char, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT8_CLASS: {
-      GtMatlabVolumeInterference<unsigned char, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned char, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT16_CLASS: {
-      GtMatlabVolumeInterference<short int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<short int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT16_CLASS: {
-      GtMatlabVolumeInterference<unsigned short int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned short int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT32_CLASS: {
-      GtMatlabVolumeInterference<int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT32_CLASS: {
-      GtMatlabVolumeInterference<unsigned int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT64_CLASS: {
-      GtMatlabVolumeInterference<long int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<long int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT64_CLASS: {
-      GtMatlabVolumeInterference<unsigned long int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned long int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     default:
       break;
   }
 }
-
diff --git a/zUtil_Cxx/internal_gtAssignGrainToVol_sum.cpp b/zUtil_Cxx/internal_gtAssignGrainToVol_sum.cpp
index 37c51ac620e04e0a77c184a3c518bc8a424b512f..83c2cc39fd6369dbefc3a906bebdbf5a81d44b2a 100755
--- a/zUtil_Cxx/internal_gtAssignGrainToVol_sum.cpp
+++ b/zUtil_Cxx/internal_gtAssignGrainToVol_sum.cpp
@@ -29,6 +29,35 @@ public:
   { }
 };
 
+template<typename TypeOut, typename TypeIn>
+class GtMatlabVolumeSumValue : public GtMatlabVolumeLooper<TypeOut, TypeIn> {
+protected:
+  void doWork(TypeOut & dest, const TypeIn & orig)
+  {
+    if (orig) { dest += (TypeOut) orig; }
+  }
+public:
+  GtMatlabVolumeSumValue(void * _out, void * _vol, const double * limits,
+      const mwSize * dims_out, const mwSize * dims_vol)
+  : GtMatlabVolumeLooper<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol, 0)
+  { }
+};
+
+template<typename TypeOut, typename TypeIn>
+inline void
+assign(void * _out, void * _vol, const double * limits,
+    const mwSize * dims_out, const mwSize * dims_vol, const double & indx)
+{
+  GtMatlabVolumeLooper<TypeOut, TypeIn> * looper = NULL;
+  if (indx) {
+    looper = new GtMatlabVolumeSum<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol, indx);
+  } else {
+    looper = new GtMatlabVolumeSumValue<TypeOut, TypeIn>(_out, _vol, limits, dims_out, dims_vol);
+  }
+  looper->loopOver3D();
+  delete looper;
+}
+
 void mexFunction( int nlhs, mxArray * plhs[], int nrhs, const mxArray * prhs[] )
 {
   /* Incoming image dimensions */
@@ -100,63 +129,43 @@ extDoAssignment(mxArray * output, void * vol, const double * limits,
 {
   switch(mxGetClassID(output)) {
     case mxDOUBLE_CLASS: {
-      GtMatlabVolumeSum<double, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<double, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxSINGLE_CLASS: {
-      GtMatlabVolumeSum<float, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<float, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT8_CLASS: {
-      GtMatlabVolumeSum<char, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<char, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT8_CLASS: {
-      GtMatlabVolumeSum<unsigned char, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned char, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT16_CLASS: {
-      GtMatlabVolumeSum<short int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<short int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT16_CLASS: {
-      GtMatlabVolumeSum<unsigned short int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned short int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT32_CLASS: {
-      GtMatlabVolumeSum<int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT32_CLASS: {
-      GtMatlabVolumeSum<unsigned int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxINT64_CLASS: {
-      GtMatlabVolumeSum<long int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<long int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     case mxUINT64_CLASS: {
-      GtMatlabVolumeSum<unsigned long int, TypeIn>
-                looper(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
-      looper.loopOver3D();
+      assign<unsigned long int, TypeIn>(mxGetData(output), vol, limits, dimsOut, dimsVol, index);
       break;
     }
     default:
diff --git a/zUtil_Imaging/gtPlaceSubVolume.m b/zUtil_Imaging/gtPlaceSubVolume.m
index 4144d0fec921ef3762cc792ba852dc287b4d5ab6..c355194734b99fc807117c97ca949d6dfb03c92b 100644
--- a/zUtil_Imaging/gtPlaceSubVolume.m
+++ b/zUtil_Imaging/gtPlaceSubVolume.m
@@ -15,23 +15,18 @@ function output = gtPlaceSubVolume(output, input, shift, index, assign_op, use_c
 %     Note - this used to add on volume to another, rather than place one image
 %     in the other.  I have modified so it no longer adds.  Hope this doesn't
 %     break anything.
+%
+%     Note 2 - If index is 0, the volume will be copied as it is, otherwise all
+%     the voxels will have value of the index
 
     if (~exist('use_c_functions', 'var'))
         use_c_functions = true;
     end
-
     if (~exist('index', 'var'))
-        elem = input(find(input ~= 0, 1));
-        if (all(input(input ~= 0) == elem))
-            index = elem;
-        else
-            warning('PLACE_SUB_VOL:low_performance', ...
-                    'Not using C functions, could be slower');
-            use_c_functions = false;
-        end
+        index = 0;
     end
     if (~exist('assign_op', 'var'))
-        assign_op = 'sum';
+        assign_op = 'assign';
     end
 
     % output volume limits