diff --git a/pysaliency/__init__.py b/pysaliency/__init__.py
index e8c9f54..1d4f639 100755
--- a/pysaliency/__init__.py
+++ b/pysaliency/__init__.py
@@ -5,6 +5,7 @@
 from . import models
 from . import external_models
 from . import external_datasets
+from . import utils
 
 from .datasets import (
     Fixations,
diff --git a/pysaliency/external_models/__init__.py b/pysaliency/external_models/__init__.py
index 69803eb..e976563 100644
--- a/pysaliency/external_models/__init__.py
+++ b/pysaliency/external_models/__init__.py
@@ -11,4 +11,9 @@
     CovSal,
 )
 
-from .utils import ExternalModelMixin
\ No newline at end of file
+from .deepgaze import (
+    DeepGazeI,
+    DeepGazeIIE,
+)
+
+from .utils import ExternalModelMixin
diff --git a/pysaliency/models.py b/pysaliency/models.py
index 0366151..5ad5cca 100755
--- a/pysaliency/models.py
+++ b/pysaliency/models.py
@@ -416,7 +416,9 @@ def log_likelihoods(self, stimuli, fixations, verbose=False):
         for stimulus_index in stimulus_indices:
             stimulus_shapes[stimulus_index] = stimuli.stimulus_objects[stimulus_index].size
 
-        stimulus_log_likelihoods = -np.log(stimulus_shapes).sum(axis=1)
+        with np.errstate(divide='ignore'):  # ignore log(0) warnings, we won't use them anyway
+            stimulus_log_likelihoods = -np.log(stimulus_shapes).sum(axis=1)
+
         return stimulus_log_likelihoods[fixations.n]
 
 
diff --git a/pysaliency/saliency_map_models.py b/pysaliency/saliency_map_models.py
index b02c0ac..ff7d005 100644
--- a/pysaliency/saliency_map_models.py
+++ b/pysaliency/saliency_map_models.py
@@ -410,9 +410,14 @@ def AUC_per_image(self, stimuli, fixations, nonfixations='uniform', thresholds='
             nonfixations = FullShuffledNonfixationProvider(stimuli, fixations)
 
         for n in tqdm(range(len(stimuli)), disable=not verbose):
+            inds = fixations.n == n
+            if not inds.sum():
+                rocs_per_image.append(np.nan)
+                continue
+
             out = self.saliency_map(stimuli.stimulus_objects[n])
             check_prediction_shape(out, stimuli[n])
-            inds = fixations.n == n
+
             positives = np.asarray(out[fixations.y_int[inds], fixations.x_int[inds]])
             if nonfixations == 'uniform':
                 negatives = out.flatten()
@@ -482,6 +487,14 @@ def AUC(self, stimuli, fixations, nonfixations='uniform', average='fixation', th
 
             return np.average(aucs, weights=weights)
         elif average == 'image':
+            stimulus_indices = set(fixations.n)
+            nan_value_indices = np.nonzero(np.isnan(aucs))[0]
+
+            if stimulus_indices.intersection(nan_value_indices):
+                raise ValueError("Some images with fixations returned AUC of nan, which should not happen")
+
+            aucs = aucs[~np.isnan(aucs)]
+
             return np.mean(aucs)
         else:
             raise ValueError(average)
diff --git a/tests/test_saliency_map_models.py b/tests/test_saliency_map_models.py
index 20f4e98..3131949 100644
--- a/tests/test_saliency_map_models.py
+++ b/tests/test_saliency_map_models.py
@@ -184,6 +184,31 @@ def test_auc_gauss(stimuli, fixation_trains):
     np.testing.assert_allclose(aucs_single, aucs_combined)
 
 
+def test_auc_per_image(stimuli, fixation_trains):
+    gsmm = GaussianSaliencyMapModel()
+
+    aucs = gsmm.AUC_per_image(stimuli, fixation_trains)
+    np.testing.assert_allclose(aucs,
+                               [0.196625, 0.313125],
+                               rtol=1e-6)
+
+
+def test_auc_per_image_images_without_fixations(stimuli, fixation_trains):
+    gsmm = GaussianSaliencyMapModel()
+
+    aucs = gsmm.AUC_per_image(stimuli, fixation_trains[:5],)
+    np.testing.assert_allclose(aucs,
+                               [0.196625, np.nan],
+                               rtol=1e-6)
+
+
+def test_auc_image_average_with_images_without_fixations(stimuli, fixation_trains):
+    gsmm = GaussianSaliencyMapModel()
+
+    auc = gsmm.AUC(stimuli, fixation_trains[:5], average='image')
+    np.testing.assert_allclose(auc, 0.196625, rtol=1e-6)
+
+
 def test_nss_gauss(stimuli, fixation_trains):
     gsmm = GaussianSaliencyMapModel()