diff --git a/TimeSeries/s2theia.py b/TimeSeries/s2theia.py
index fa6665ddf9e7000572a87326e363eb341d870c4a..39593e2c02b22f55eae467783853a3f5bbdeecc6 100644
--- a/TimeSeries/s2theia.py
+++ b/TimeSeries/s2theia.py
@@ -793,9 +793,9 @@ class S2TheiaPipeline:
     def set_output_epsg(self, epsg):
         self.output_epsg = epsg
 
-    def extract_feature_set(self, out_fld, feat_list=None, mosaicking=None, store_gapfill=False,
-                            align=False, align_to=None, align_to_band=3, align_using_band=3,
-                            warp_to=None, warp_to_band=1, warp_using_band=3):
+    def extract_feature_set(self, out_fld, feat_list=None, mosaicking=None, store_gapfill=True,
+                            align=False, align_to=None, align_to_band=3, align_using_band=3, output_aligned=None,
+                            warp_to=None, warp_to_band=1, warp_using_band=3, output_warped=None):
         out = []
         stack_name = ''
         if self.output_dates is not None:
@@ -814,9 +814,17 @@ class S2TheiaPipeline:
                             raise ValueError("Provided string is not a valid date nor a valid file.")
                     elif align_to is None:
                         t.rigid_align(match_band=align_using_band-1)
+                    if output_aligned is not None:
+                        if not os.path.exists(output_aligned):
+                            os.makedirs(output_aligned)
+                        t.write_outputs(output_aligned, update_pipe=True, flag_nodata=True)
                 t.reproject(self.output_epsg)
                 if warp_to is not None:
                     t.coregister(warp_to, warp_to_band, warp_using_band, t.temp_fld)
+                    if output_warped is not None:
+                        if not os.path.exists(output_warped):
+                            os.makedirs(output_warped)
+                        t.write_outputs(output_warped, update_pipe=True, flag_nodata=True)
                 t.gapfill(self.output_dates, store_gapfill)
                 stack_name = t.generate_feature_stack(feat_list)
                 out.append(t.write_outputs(out_fld))