diff --git a/bob/bio/face/algorithm/GaborJet.py b/bob/bio/face/algorithm/GaborJet.py
index 0e790b36d1d9fec06ca68a548a5a0bf6b46d57bf..54879038e63d6171a80bf53887958191d87f1463 100644
--- a/bob/bio/face/algorithm/GaborJet.py
+++ b/bob/bio/face/algorithm/GaborJet.py
@@ -9,9 +9,11 @@ import numpy
 import math
 
 from bob.bio.base.algorithm import Algorithm
+from bob.bio.face.extractor import GridGraph
 
-class GaborJet (Algorithm):
-  """Computes a comparison of lists of Gabor jets using a similarity function of :py:class:`bob.ip.gabor.Similarity`.
+
+class GaborJet(Algorithm):
+    """Computes a comparison of lists of Gabor jets using a similarity function of :py:class:`bob.ip.gabor.Similarity`.
 
   The model enrollment simply stores all extracted Gabor jets for all enrollment features.
   By default (i.e., ``multiple_feature_scoring = 'max_jet'``), the scoring uses an advanced local strategy.
@@ -40,84 +42,100 @@ class GaborJet (Algorithm):
     Please assure that this class and the :py:class:`bob.bio.face.extractor.GridGraph` class get the same configuration, otherwise unexpected things might happen.
   """
 
-  def __init__(
-      self,
-      # parameters for the tool
-      gabor_jet_similarity_type,
-      multiple_feature_scoring = 'max_jet',
-      # some similarity functions might need a GaborWaveletTransform class, so we have to provide the parameters here as well...
-      gabor_directions = 8,
-      gabor_scales = 5,
-      gabor_sigma = 2. * math.pi,
-      gabor_maximum_frequency = math.pi / 2.,
-      gabor_frequency_step = math.sqrt(.5),
-      gabor_power_of_k = 0,
-      gabor_dc_free = True
-  ):
-
-    # call base class constructor
-    Algorithm.__init__(
+    def __init__(
         self,
-
-        gabor_jet_similarity_type = gabor_jet_similarity_type,
-        multiple_feature_scoring = multiple_feature_scoring,
-        gabor_directions = gabor_directions,
-        gabor_scales = gabor_scales,
-        gabor_sigma = gabor_sigma,
-        gabor_maximum_frequency = gabor_maximum_frequency,
-        gabor_frequency_step = gabor_frequency_step,
-        gabor_power_of_k = gabor_power_of_k,
-        gabor_dc_free = gabor_dc_free,
-
-        multiple_model_scoring = None,
-        multiple_probe_scoring = None
-    )
-
-    # the Gabor wavelet transform; used by (some of) the Gabor jet similarities
-    gwt = bob.ip.gabor.Transform(
-        number_of_scales = gabor_scales,
-        number_of_directions = gabor_directions,
-        sigma = gabor_sigma,
-        k_max = gabor_maximum_frequency,
-        k_fac = gabor_frequency_step,
-        power_of_k = gabor_power_of_k,
-        dc_free = gabor_dc_free
-    )
-
-    # jet comparison function
-    self.similarity_function = bob.ip.gabor.Similarity(gabor_jet_similarity_type, gwt)
-
-    # how to proceed with multiple features per model
-    self.jet_scoring = {
-        'average_model' : None, # compute an average model
-        'average' : numpy.average, # compute the average similarity
-        'min_jet' : min, # for each jet location, compute the minimum similarity
-        'max_jet' : max, # for each jet location, compute the maximum similarity
-        'med_jet' : numpy.median, # for each jet location, compute the median similarity
-        'min_graph' : numpy.average, # for each model graph, compute the minimum average similarity
-        'max_graph' : numpy.average, # for each model graph, compute the maximum average similarity
-        'med_graph' : numpy.average, # for each model graph, compute the median average similarity
-    }[multiple_feature_scoring]
-
-    self.graph_scoring = {
-        'average_model' : None, # compute an average model
-        'average' : numpy.average, # compute the average similarity
-        'min_jet' : numpy.average, # for each jet location, compute the minimum similarity
-        'max_jet' : numpy.average, # for each jet location, compute the maximum similarity
-        'med_jet' : numpy.average, # for each jet location, compute the median similarity
-        'min_graph' : min, # for each model graph, compute the minimum average similarity
-        'max_graph' : max, # for each model graph, compute the maximum average similarity
-        'med_graph' : numpy.median, # for each model graph, compute the median average similarity
-    }[multiple_feature_scoring]
-
-
-  def _check_feature(self, feature):
-    assert isinstance(feature, list)
-    assert len(feature)
-    assert all(isinstance(f, bob.ip.gabor.Jet) for f in feature)
-
-  def enroll(self, enroll_features):
-    """enroll(enroll_features) -> model
+        # parameters for the tool
+        gabor_jet_similarity_type,
+        multiple_feature_scoring="max_jet",
+        # some similarity functions might need a GaborWaveletTransform class, so we have to provide the parameters here as well...
+        gabor_directions=8,
+        gabor_scales=5,
+        gabor_sigma=2.0 * math.pi,
+        gabor_maximum_frequency=math.pi / 2.0,
+        gabor_frequency_step=math.sqrt(0.5),
+        gabor_power_of_k=0,
+        gabor_dc_free=True,
+    ):
+
+        # call base class constructor
+        Algorithm.__init__(
+            self,
+            gabor_jet_similarity_type=gabor_jet_similarity_type,
+            multiple_feature_scoring=multiple_feature_scoring,
+            gabor_directions=gabor_directions,
+            gabor_scales=gabor_scales,
+            gabor_sigma=gabor_sigma,
+            gabor_maximum_frequency=gabor_maximum_frequency,
+            gabor_frequency_step=gabor_frequency_step,
+            gabor_power_of_k=gabor_power_of_k,
+            gabor_dc_free=gabor_dc_free,
+            multiple_model_scoring=None,
+            multiple_probe_scoring=None,
+        )
+
+        self.gabor_jet_similarity_type = gabor_jet_similarity_type
+        self.multiple_feature_scoring = multiple_feature_scoring
+        self.gabor_directions = gabor_directions
+        self.gabor_scales = gabor_scales
+        self.gabor_sigma = gabor_sigma
+        self.gabor_maximum_frequency = gabor_maximum_frequency
+        self.gabor_frequency_step = gabor_frequency_step
+        self.gabor_power_of_k = gabor_power_of_k
+        self.gabor_dc_free = gabor_dc_free
+
+        self.gabor_jet_similarity_type = gabor_jet_similarity_type
+
+        self._init_non_pickables()
+
+    def _init_non_pickables(self):
+        # the Gabor wavelet transform; used by (some of) the Gabor jet similarities
+        self.gwt = bob.ip.gabor.Transform(
+            number_of_scales=self.gabor_scales,
+            number_of_directions=self.gabor_directions,
+            sigma=self.gabor_sigma,
+            k_max=self.gabor_maximum_frequency,
+            k_fac=self.gabor_frequency_step,
+            power_of_k=self.gabor_power_of_k,
+            dc_free=self.gabor_dc_free,
+        )
+
+        # jet comparison function
+        self.similarity_function = bob.ip.gabor.Similarity(
+            self.gabor_jet_similarity_type, self.gwt
+        )
+
+        # how to proceed with multiple features per model
+        self.jet_scoring = {
+            "average_model": None,  # compute an average model
+            "average": numpy.average,  # compute the average similarity
+            "min_jet": min,  # for each jet location, compute the minimum similarity
+            "max_jet": max,  # for each jet location, compute the maximum similarity
+            "med_jet": numpy.median,  # for each jet location, compute the median similarity
+            "min_graph": numpy.average,  # for each model graph, compute the minimum average similarity
+            "max_graph": numpy.average,  # for each model graph, compute the maximum average similarity
+            "med_graph": numpy.average,  # for each model graph, compute the median average similarity
+        }[self.multiple_feature_scoring]
+
+        self.graph_scoring = {
+            "average_model": None,  # compute an average model
+            "average": numpy.average,  # compute the average similarity
+            "min_jet": numpy.average,  # for each jet location, compute the minimum similarity
+            "max_jet": numpy.average,  # for each jet location, compute the maximum similarity
+            "med_jet": numpy.average,  # for each jet location, compute the median similarity
+            "min_graph": min,  # for each model graph, compute the minimum average similarity
+            "max_graph": max,  # for each model graph, compute the maximum average similarity
+            "med_graph": numpy.median,  # for each model graph, compute the median average similarity
+        }[self.multiple_feature_scoring]
+
+    def _check_feature(self, feature):
+
+        assert isinstance(feature, list) or isinstance(feature, numpy.ndarray)
+        assert len(feature)
+        feature = GridGraph.serialize_jets(feature)
+        assert all(isinstance(f, bob.ip.gabor.Jet) for f in feature)
+
+    def enroll(self, enroll_features):
+        """enroll(enroll_features) -> model
 
     Enrolls the model using one of several strategies.
     Commonly, the bunch graph strategy [WFK97]_ is applied, by storing several Gabor jets for each node.
@@ -138,22 +156,26 @@ class GaborJet (Algorithm):
       Each sub-list contains a list of jets, which correspond to the same node.
       When ``multiple_feature_scoring = 'average_model'`` each sub-list contains a single :py:class:`bob.ip.gabor.Jet`.
     """
-    [self._check_feature(feature) for feature in enroll_features]
-    assert len(enroll_features)
-    assert all(len(feature) == len(enroll_features[0]) for feature in enroll_features)
-
-    # re-organize the jets to have a collection of jets per node
-    jets_per_node = [[enroll_features[g][n] for g in range(len(enroll_features))] for n in range(len(enroll_features[0]))]
+        [self._check_feature(feature) for feature in enroll_features]
+        assert len(enroll_features)
+        assert all(
+            len(feature) == len(enroll_features[0]) for feature in enroll_features
+        )
 
-    if self.jet_scoring is not None:
-      return jets_per_node
+        # re-organize the jets to have a collection of jets per node
+        jets_per_node = [
+            [enroll_features[g][n] for g in range(len(enroll_features))]
+            for n in range(len(enroll_features[0]))
+        ]
 
-    # compute average model, and keep a list with a single jet per node
-    return [[bob.ip.gabor.Jet(jets_per_node[n])] for n in range(len(jets_per_node))]
+        if self.jet_scoring is not None:
+            return jets_per_node
 
+        # compute average model, and keep a list with a single jet per node
+        return [[bob.ip.gabor.Jet(jets_per_node[n])] for n in range(len(jets_per_node))]
 
-  def write_model(self, model, model_file):
-    """Writes the model enrolled by the :py:meth:`enroll` function to the given file.
+    def write_model(self, model, model_file):
+        """Writes the model enrolled by the :py:meth:`enroll` function to the given file.
 
     **Parameters:**
 
@@ -163,20 +185,19 @@ class GaborJet (Algorithm):
     model_file : str or :py:class:`bob.io.base.HDF5File`
       The name of the file or the file opened for writing.
     """
-    f = bob.io.base.HDF5File(model_file, 'w')
-    # several model graphs
-    f.set("NumberOfNodes", len(model))
-    for g in range(len(model)):
-      name = "Node-" + str(g+1)
-      f.create_group(name)
-      f.cd(name)
-      bob.ip.gabor.save_jets(model[g], f)
-      f.cd("..")
-    f.close()
-
-
-  def read_model(self, model_file):
-    """read_model(model_file) -> model
+        f = bob.io.base.HDF5File(model_file, "w")
+        # several model graphs
+        f.set("NumberOfNodes", len(model))
+        for g in range(len(model)):
+            name = "Node-" + str(g + 1)
+            f.create_group(name)
+            f.cd(name)
+            bob.ip.gabor.save_jets(model[g], f)
+            f.cd("..")
+        f.close()
+
+    def read_model(self, model_file):
+        """read_model(model_file) -> model
 
     Reads the model written by the :py:meth:`write_model` function from the given file.
 
@@ -190,19 +211,18 @@ class GaborJet (Algorithm):
     model : [[:py:class:`bob.ip.gabor.Jet`]]
       The list of Gabor jets read from file.
     """
-    f = bob.io.base.HDF5File(model_file)
-    count = f.get("NumberOfNodes")
-    model = []
-    for g in range(count):
-      name = "Node-" + str(g+1)
-      f.cd(name)
-      model.append(bob.ip.gabor.load_jets(f))
-      f.cd("..")
-    return model
-
-
-  def score(self, model, probe):
-    """score(model, probe) -> score
+        f = bob.io.base.HDF5File(model_file)
+        count = f.get("NumberOfNodes")
+        model = []
+        for g in range(count):
+            name = "Node-" + str(g + 1)
+            f.cd(name)
+            model.append(GridGraph.serialize_jets(bob.ip.gabor.load_jets(f)))
+            f.cd("..")
+        return model
+
+    def score(self, model, probe):
+        """score(model, probe) -> score
 
     Computes the score of the probe and the model using the desired Gabor jet similarity function and the desired score fusion strategy.
 
@@ -219,19 +239,23 @@ class GaborJet (Algorithm):
     score : float
       The fused similarity score.
     """
-    self._check_feature(probe)
-    [self._check_feature(m) for m in model]
-    assert len(model) == len(probe)
-
-    # select jet score averaging function
-    jet_scoring = numpy.average if self.jet_scoring is None else self.jet_scoring
-    graph_scoring = numpy.average if self.graph_scoring is None else self.graph_scoring
-    local_scores = [jet_scoring([self.similarity_function(m, pro) for m in mod]) for mod, pro in zip(model, probe)]
-    return graph_scoring(local_scores)
-
-
-  def score_for_multiple_probes(self, model, probes):
-    """score(model, probes) -> score
+        self._check_feature(probe)
+        [self._check_feature(m) for m in model]
+        assert len(model) == len(probe)
+
+        # select jet score averaging function
+        jet_scoring = numpy.average if self.jet_scoring is None else self.jet_scoring
+        graph_scoring = (
+            numpy.average if self.graph_scoring is None else self.graph_scoring
+        )
+        local_scores = [
+            jet_scoring([self.similarity_function(m, pro) for m in mod])
+            for mod, pro in zip(model, probe)
+        ]
+        return graph_scoring(local_scores)
+
+    def score_for_multiple_probes(self, model, probes):
+        """score(model, probes) -> score
 
     This function computes the score between the given model graph(s) and several given probe graphs.
     The same local scoring strategy as for several model jets is applied, but this time the local scoring strategy is applied between all graphs from the model and probes.
@@ -251,22 +275,83 @@ class GaborJet (Algorithm):
     score : float
       The fused similarity score.
     """
-    [self._check_feature(probe) for probe in probes]
-    [self._check_feature(m) for m in model]
-    assert all(len(model) == len(probe) for probe in probes)
-
-    jet_scoring = numpy.average if self.jet_scoring is None else self.jet_scoring
-    graph_scoring = numpy.average if self.graph_scoring is None else self.graph_scoring
-    local_scores = [jet_scoring([self.similarity_function(m, probe[n]) for m in model[n] for probe in probes]) for n in range(len(model))]
-    return graph_scoring(local_scores)
-
-
-  # overwrite functions to avoid them being documented.
-  def train_projector(*args, **kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
-  def load_projector(*args, **kwargs) : pass
-  def project(*args, **kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
-  def write_feature(*args, **kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
-  def read_feature(*args, **kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
-  def train_enroller(*args, **kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
-  def load_enroller(*args, **kwargs) : pass
-  def score_for_multiple_models(*args, **kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
+        [self._check_feature(probe) for probe in probes]
+        [self._check_feature(m) for m in model]
+        assert all(len(model) == len(probe) for probe in probes)
+
+        jet_scoring = numpy.average if self.jet_scoring is None else self.jet_scoring
+        graph_scoring = (
+            numpy.average if self.graph_scoring is None else self.graph_scoring
+        )
+        local_scores = [
+            jet_scoring(
+                [
+                    self.similarity_function(m, probe[n])
+                    for m in model[n]
+                    for probe in probes
+                ]
+            )
+            for n in range(len(model))
+        ]
+        return graph_scoring(local_scores)
+
+    def score_for_multiple_models(self, models, probe):
+
+        self._check_feature(probe)
+        [self._check_feature(m) for model in models for m in model]
+
+        jet_scoring = numpy.average if self.jet_scoring is None else self.jet_scoring
+        graph_scoring = (
+            numpy.average if self.graph_scoring is None else self.graph_scoring
+        )
+        scores = []
+        for model in models:
+            local_scores = local_scores = [
+                jet_scoring([self.similarity_function(m, pro) for m in mod])
+                for mod, pro in zip(model, probe)
+            ]
+            scores.append(graph_scoring(local_scores))
+
+        return scores
+
+    # overwrite functions to avoid them being documented.
+    def train_projector(*args, **kwargs):
+        raise NotImplementedError(
+            "This function is not implemented and should not be called."
+        )
+
+    def load_projector(*args, **kwargs):
+        pass
+
+    def project(*args, **kwargs):
+        raise NotImplementedError(
+            "This function is not implemented and should not be called."
+        )
+
+    def write_feature(*args, **kwargs):
+        raise NotImplementedError(
+            "This function is not implemented and should not be called."
+        )
+
+    def read_feature(*args, **kwargs):
+        raise NotImplementedError(
+            "This function is not implemented and should not be called."
+        )
+
+    def train_enroller(*args, **kwargs):
+        raise NotImplementedError(
+            "This function is not implemented and should not be called."
+        )
+
+    def load_enroller(*args, **kwargs):
+        pass
+
+    def __getstate__(self):
+        d = dict(self.__dict__)
+        d.pop("gwt")
+        d.pop("similarity_function")
+        return d
+
+    def __setstate__(self, d):
+        self.__dict__ = d
+        self._init_non_pickables()
diff --git a/bob/bio/face/extractor/GridGraph.py b/bob/bio/face/extractor/GridGraph.py
index 387d86bec66b3b5e4f38d09240d50233c7b6508b..a3b6e24807f60a9cfed7ed9e847a925f7950be31 100644
--- a/bob/bio/face/extractor/GridGraph.py
+++ b/bob/bio/face/extractor/GridGraph.py
@@ -9,8 +9,9 @@ import numpy
 import math
 from bob.bio.base.extractor import Extractor
 
-class GridGraph (Extractor):
-  """Extracts Gabor jets in a grid structure [GHW12]_ using functionalities from :ref:`bob.ip.gabor <bob.ip.gabor>`.
+
+class GridGraph(Extractor):
+    """Extracts Gabor jets in a grid structure [GHW12]_ using functionalities from :ref:`bob.ip.gabor <bob.ip.gabor>`.
 
   The grid can be either aligned to the eye locations (in which case the grid might be rotated), or a fixed grid graph can be extracted.
 
@@ -47,129 +48,155 @@ class GridGraph (Extractor):
     If ``None``, it is calculated automatically to equally cover the whole image.
   """
 
-  def __init__(
-      self,
-      # Gabor parameters
-      gabor_directions = 8,
-      gabor_scales = 5,
-      gabor_sigma = 2. * math.pi,
-      gabor_maximum_frequency = math.pi / 2.,
-      gabor_frequency_step = math.sqrt(.5),
-      gabor_power_of_k = 0,
-      gabor_dc_free = True,
-
-      # what kind of information to extract
-      normalize_gabor_jets = True,
-
-      # setup of the aligned grid
-      eyes = None, # if set, the grid setup will be aligned to the eye positions {'leye' : LEFT_EYE_POS, 'reye' : RIGHT_EYE_POS},
-      nodes_between_eyes = 4,
-      nodes_along_eyes = 2,
-      nodes_above_eyes = 3,
-      nodes_below_eyes = 7,
-
-      # setup of static grid
-      node_distance = None,    # one or two integral values
-      first_node = None,       # one or two integral values, or None -> automatically determined
-  ):
-
-    # call base class constructor
-    Extractor.__init__(
+    def __init__(
         self,
-
-        gabor_directions = gabor_directions,
-        gabor_scales = gabor_scales,
-        gabor_sigma = gabor_sigma,
-        gabor_maximum_frequency = gabor_maximum_frequency,
-        gabor_frequency_step = gabor_frequency_step,
-        gabor_power_of_k = gabor_power_of_k,
-        gabor_dc_free = gabor_dc_free,
-        normalize_gabor_jets = normalize_gabor_jets,
-        eyes = eyes,
-        nodes_between_eyes = nodes_between_eyes,
-        nodes_along_eyes = nodes_along_eyes,
-        nodes_above_eyes = nodes_above_eyes,
-        nodes_below_eyes = nodes_below_eyes,
-        node_distance = node_distance,
-        first_node = first_node
-    )
-
-    # create Gabor wavelet transform class
-    self.gwt = bob.ip.gabor.Transform(
-        number_of_scales = gabor_scales,
-        number_of_directions = gabor_directions,
-        sigma = gabor_sigma,
-        k_max = gabor_maximum_frequency,
-        k_fac = gabor_frequency_step,
-        power_of_k = gabor_power_of_k,
-        dc_free = gabor_dc_free
-    )
-
-    # create graph extractor
-    if eyes is not None:
-      self._aligned_graph = bob.ip.gabor.Graph(
-          righteye = [int(e) for e in eyes['reye']],
-          lefteye = [int(e) for e in eyes['leye']],
-          between = int(nodes_between_eyes),
-          along = int(nodes_along_eyes),
-          above = int(nodes_above_eyes),
-          below = int(nodes_below_eyes)
-      )
-    else:
-      if node_distance is None:
-        raise ValueError("Please specify either 'eyes' or the grid parameters 'node_distance' (and 'first_node')!")
-      self._aligned_graph = None
-      self._last_image_resolution = None
-      self.first_node = first_node
-      self.node_distance = node_distance
-      if isinstance(self.node_distance, (int, float)):
-         self.node_distance = (int(self.node_distance), int(self.node_distance))
-
-    self.normalize_jets = normalize_gabor_jets
-    self.trafo_image = None
-
-  def _extractor(self, image):
-    """Creates an extractor based on the given image.
+        # Gabor parameters
+        gabor_directions=8,
+        gabor_scales=5,
+        gabor_sigma=2.0 * math.pi,
+        gabor_maximum_frequency=math.pi / 2.0,
+        gabor_frequency_step=math.sqrt(0.5),
+        gabor_power_of_k=0,
+        gabor_dc_free=True,
+        # what kind of information to extract
+        normalize_gabor_jets=True,
+        # setup of the aligned grid
+        eyes=None,  # if set, the grid setup will be aligned to the eye positions {'leye' : LEFT_EYE_POS, 'reye' : RIGHT_EYE_POS},
+        nodes_between_eyes=4,
+        nodes_along_eyes=2,
+        nodes_above_eyes=3,
+        nodes_below_eyes=7,
+        # setup of static grid
+        node_distance=None,  # one or two integral values
+        first_node=None,  # one or two integral values, or None -> automatically determined
+    ):
+
+        # call base class constructor
+        Extractor.__init__(
+            self,
+            gabor_directions=gabor_directions,
+            gabor_scales=gabor_scales,
+            gabor_sigma=gabor_sigma,
+            gabor_maximum_frequency=gabor_maximum_frequency,
+            gabor_frequency_step=gabor_frequency_step,
+            gabor_power_of_k=gabor_power_of_k,
+            gabor_dc_free=gabor_dc_free,
+            normalize_gabor_jets=normalize_gabor_jets,
+            eyes=eyes,
+            nodes_between_eyes=nodes_between_eyes,
+            nodes_along_eyes=nodes_along_eyes,
+            nodes_above_eyes=nodes_above_eyes,
+            nodes_below_eyes=nodes_below_eyes,
+            node_distance=node_distance,
+            first_node=first_node,
+        )
+
+        self.gabor_directions = gabor_directions
+        self.gabor_scales = gabor_scales
+        self.gabor_sigma = gabor_sigma
+        self.gabor_maximum_frequency = gabor_maximum_frequency
+        self.gabor_frequency_step = gabor_frequency_step
+        self.gabor_power_of_k = gabor_power_of_k
+        self.gabor_dc_free = gabor_dc_free
+        self.normalize_gabor_jets = normalize_gabor_jets
+        self.eyes = eyes
+        self.nodes_between_eyes = nodes_between_eyes
+        self.nodes_along_eyes = nodes_along_eyes
+        self.nodes_above_eyes = nodes_above_eyes
+        self.nodes_below_eyes = nodes_below_eyes
+        self.node_distance = node_distance
+        self.first_node = first_node
+
+        self.normalize_jets = normalize_gabor_jets
+        self.trafo_image = None
+        self._init_non_pickables()
+
+    def _init_non_pickables(self):
+
+        # create Gabor wavelet transform class
+        self.gwt = bob.ip.gabor.Transform(
+            number_of_scales=self.gabor_scales,
+            number_of_directions=self.gabor_directions,
+            sigma=self.gabor_sigma,
+            k_max=self.gabor_maximum_frequency,
+            k_fac=self.gabor_frequency_step,
+            power_of_k=self.gabor_power_of_k,
+            dc_free=self.gabor_dc_free,
+        )
+
+        # create graph extractor
+        if self.eyes is not None:
+            self._aligned_graph = bob.ip.gabor.Graph(
+                righteye=[int(e) for e in self.eyes["reye"]],
+                lefteye=[int(e) for e in self.eyes["leye"]],
+                between=int(self.nodes_between_eyes),
+                along=int(self.nodes_along_eyes),
+                above=int(self.nodes_above_eyes),
+                below=int(self.nodes_below_eyes),
+            )
+        else:
+            if self.node_distance is None:
+                raise ValueError(
+                    "Please specify either 'eyes' or the grid parameters 'node_distance' (and 'first_node')!"
+                )
+            self._aligned_graph = None
+            self._last_image_resolution = None
+            if isinstance(self.node_distance, (int, float)):
+                self.node_distance = (int(self.node_distance), int(self.node_distance))
+
+    def _extractor(self, image):
+        """Creates an extractor based on the given image.
     If an aligned graph was specified in the constructor, it is simply returned.
     Otherwise the resolution of the given image is used to create a graph extractor.
     If the ``first_node`` was not specified, it is calculated automatically.
     """
 
-    if self.trafo_image is None or self.trafo_image.shape[1:3] != image.shape:
-      # create trafo image
-      self.trafo_image = numpy.ndarray((self.gwt.number_of_wavelets, image.shape[0], image.shape[1]), numpy.complex128)
-
-    if self._aligned_graph is not None:
-      return self._aligned_graph
-
-    # check if a new extractor needs to be created
-    if self._last_image_resolution != image.shape:
-      self._last_image_resolution = image.shape
-      if self.first_node is None:
-        # automatically compute the first node
-        first_node = [0,0]
-        for i in (0,1):
-          offset = int((image.shape[i] - int(image.shape[i]/self.node_distance[i])*self.node_distance[i]) / 2)
-          if offset < self.node_distance[i]//2: # This is not tested, but should ALWAYS be the case.
-            offset += self.node_distance[i]//2
-          first_node[i] = offset
-      else:
-        first_node = self.first_node
-      # .. and the last node
-      last_node = tuple([int(image.shape[i] - max(first_node[i],1)) for i in (0,1)])
-
-      # take the specified nodes
-      self._graph = bob.ip.gabor.Graph(
-          first = first_node,
-          last = last_node,
-          step = self.node_distance
-      )
-
-    return self._graph
-
-
-  def __call__(self, image):
-    """__call__(image) -> feature
+        if self.trafo_image is None or self.trafo_image.shape[1:3] != image.shape:
+            # create trafo image
+            self.trafo_image = numpy.ndarray(
+                (self.gwt.number_of_wavelets, image.shape[0], image.shape[1]),
+                numpy.complex128,
+            )
+
+        if self._aligned_graph is not None:
+            return self._aligned_graph
+
+        # check if a new extractor needs to be created
+        if self._last_image_resolution != image.shape:
+            self._last_image_resolution = image.shape
+            if self.first_node is None:
+                # automatically compute the first node
+                first_node = [0, 0]
+                for i in (0, 1):
+                    offset = int(
+                        (
+                            image.shape[i]
+                            - int(image.shape[i] / self.node_distance[i])
+                            * self.node_distance[i]
+                        )
+                        / 2
+                    )
+                    if (
+                        offset < self.node_distance[i] // 2
+                    ):  # This is not tested, but should ALWAYS be the case.
+                        offset += self.node_distance[i] // 2
+                    first_node[i] = offset
+            else:
+                first_node = self.first_node
+            # .. and the last node
+            last_node = tuple(
+                [int(image.shape[i] - max(first_node[i], 1)) for i in (0, 1)]
+            )
+
+            # take the specified nodes
+            self._graph = bob.ip.gabor.Graph(
+                first=first_node, last=last_node, step=self.node_distance
+            )
+
+        return self._graph
+
+    def __call__(self, image):
+        """__call__(image) -> feature
 
     Returns a list of Gabor jets extracted from the given image.
 
@@ -184,27 +211,28 @@ class GridGraph (Extractor):
       The list of Gabor jets extracted from the image.
       The 2D location of the jet's nodes is not returned.
     """
-    assert image.ndim == 2
-    assert isinstance(image, numpy.ndarray)
-    assert image.dtype == numpy.float64
 
-    extractor = self._extractor(image)
+        assert image.ndim == 2
+        assert isinstance(image, numpy.ndarray)
+        image = image.astype(numpy.float64)
+        assert image.dtype == numpy.float64
 
-    # perform Gabor wavelet transform
-    self.gwt.transform(image, self.trafo_image)
-    # extract face graph
-    jets = extractor.extract(self.trafo_image)
+        extractor = self._extractor(image)
 
-    # normalize the Gabor jets of the graph only
-    if self.normalize_jets:
-      [j.normalize() for j in jets]
+        # perform Gabor wavelet transform
+        self.gwt.transform(image, self.trafo_image)
+        # extract face graph
+        jets = extractor.extract(self.trafo_image)
 
-    # return the extracted face graph
-    return jets
+        # normalize the Gabor jets of the graph only
+        if self.normalize_jets:
+            [j.normalize() for j in jets]
 
+        # return the extracted face graph
+        return self.__class__.serialize_jets(jets)
 
-  def write_feature(self, feature, feature_file):
-    """Writes the feature extracted by the `__call__` function to the given file.
+    def write_feature(self, feature, feature_file):
+        """Writes the feature extracted by the `__call__` function to the given file.
 
     **Parameters:**
 
@@ -214,12 +242,15 @@ class GridGraph (Extractor):
     feature_file : str or :py:class:`bob.io.base.HDF5File`
       The name of the file or the file opened for writing.
     """
-    feature_file = feature_file if isinstance(feature_file, bob.io.base.HDF5File) else bob.io.base.HDF5File(feature_file, 'w')
-    bob.ip.gabor.save_jets(feature, feature_file)
+        feature_file = (
+            feature_file
+            if isinstance(feature_file, bob.io.base.HDF5File)
+            else bob.io.base.HDF5File(feature_file, "w")
+        )
+        bob.ip.gabor.save_jets(feature, feature_file)
 
-
-  def read_feature(self, feature_file):
-    """read_feature(feature_file) -> feature
+    def read_feature(self, feature_file):
+        """read_feature(feature_file) -> feature
 
     Reads the feature written by the :py:meth:`write_feature` function from the given file.
 
@@ -233,8 +264,36 @@ class GridGraph (Extractor):
     feature : [:py:class:`bob.ip.gabor.Jet`]
       The list of Gabor jets read from file.
     """
-    return bob.ip.gabor.load_jets(bob.io.base.HDF5File(feature_file))
-
-  # re-define the train function to get it non-documented
-  def train(*args,**kwargs) : raise NotImplementedError("This function is not implemented and should not be called.")
-  def load(*args,**kwargs) : pass
+        return self.__class__.serialize_jets(
+            bob.ip.gabor.load_jets(bob.io.base.HDF5File(feature_file))
+        )
+
+    # re-define the train function to get it non-documented
+    def train(*args, **kwargs):
+        raise NotImplementedError(
+            "This function is not implemented and should not be called."
+        )
+
+    def load(*args, **kwargs):
+        pass
+
+    def __getstate__(self):
+        d = dict(self.__dict__)
+        d.pop("gwt")
+        d.pop("_aligned_graph")
+        if "_graph" in d:
+            d.pop("_graph")
+        return d
+
+    def __setstate__(self, d):
+        self.__dict__ = d
+        self._init_non_pickables()
+
+    @staticmethod
+    def serialize_jets(jets):
+        serialize_jets = []
+        for jet in jets:
+            sj = bob.ip.gabor.Jet(jet.length)
+            sj.jet = jet.jet
+            serialize_jets.append(sj)
+        return serialize_jets