diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index d8cc4dfd705..b2ed441ff99 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -34,6 +34,7 @@
     "MorphologyContainer",
     "BaseHillasParametersContainer",
     "CameraHillasParametersContainer",
+    "CameraImageFitParametersContainer",
     "ImageFitParametersContainer",
     "CameraTimingParametersContainer",
     "ParticleClassificationContainer",
diff --git a/ctapipe/image/ellipsoid.py b/ctapipe/image/ellipsoid.py
index 4ec88357552..2e048f411b7 100644
--- a/ctapipe/image/ellipsoid.py
+++ b/ctapipe/image/ellipsoid.py
@@ -41,7 +41,6 @@
 
 class PDFType(Enum):
     gaussian = "gaussian"
-    laplace = "laplace"
     skewed = "skewed"
 
 
@@ -54,9 +53,8 @@ def create_initial_guess(geometry, hillas, size):
     ----------
     geometry: ctapipe.instrument.CameraGeometry
         Camera geometry, the cleaning mask should be applied to improve performance
-    image: ndarray
-        Charge in each pixel, the cleaning mask should already be applied to
-        improve performance.
+    hillas: HillasParametersContainer
+        Hillas parameters
     size: float/int
         Total charge after cleaning and dilation
 
@@ -101,8 +99,8 @@ def boundaries(geometry, image, dilated_mask, hillas, pdf):
         Charge in each pixel, no cleaning mask should be applied
     dilated_mask : ndarray
         mask (array of booleans) after image cleaning and dilation
-    x0 : dict
-       seeds of the fit
+    hillas : HillasParametersContainer
+       Hillas parameters
     pdf: PDFType instance
         e.g. PDFType("gaussian")
 
@@ -157,6 +155,11 @@ def boundaries(geometry, image, dilated_mask, hillas, pdf):
         -0.99, min(0.99, hillas.skewness + 0.3)
     )  # Guess from Hillas unit tests
 
+    if pdf == PDFType.gaussian:
+        amplitude = np.sum(row_image) / (2 * np.pi * width_min * length_min)
+    else:
+        amplitude = np.sum(row_image) / scale / (2 * np.pi * width_min)
+
     bounds = [
         (cogx_min, cogx_max),
         (cogy_min, cogy_max),
@@ -164,18 +167,9 @@ def boundaries(geometry, image, dilated_mask, hillas, pdf):
         (length_min, length_max),
         (width_min, width_max),
         (skew_min, skew_max),
+        (0, amplitude),
     ]
 
-    if pdf == PDFType.gaussian:
-        amplitude = np.sum(row_image) / (2 * np.pi * width_min * length_min)
-    elif pdf == PDFType.skewed:
-        amplitude = np.sum(row_image) / scale / (2 * np.pi * width_min)
-    else:
-        amplitude = (
-            np.sum(row_image) / scale / (np.sqrt(2 * np.pi) * np.sqrt(2) * width_min)
-        )
-
-    bounds.append((0, amplitude))
     return bounds
 
 
@@ -333,8 +327,8 @@ def likelihood(cog_x, cog_y, psi, length, width, skewness, amplitude):
     fit_phi = np.arctan2(pars[1], pars[0])
 
     b = pars[1] ** 2 + pars[0] ** 2
-    A = (-pars[1] / (b)) ** 2
-    B = (pars[0] / (b)) ** 2
+    A = (-pars[1] / b) ** 2
+    B = (pars[0] / b) ** 2
     fit_phi_err = np.sqrt(A * errors[0] ** 2 + B * errors[1] ** 2)
     fit_rcog_err = np.sqrt(
         pars[0] ** 2 / b * errors[0] ** 2 + pars[1] ** 2 / b * errors[1] ** 2
diff --git a/ctapipe/image/tests/test_ellipsoid.py b/ctapipe/image/tests/test_ellipsoid.py
index ae6eebaa5b8..0b187632d2c 100644
--- a/ctapipe/image/tests/test_ellipsoid.py
+++ b/ctapipe/image/tests/test_ellipsoid.py
@@ -213,30 +213,20 @@ def test_truncated(prod5_lst):
 def test_percentage(prod5_lst):
     geom = prod5_lst.camera.geometry
 
-    widths = u.Quantity([0.01, 0.02, 0.03, 0.07], u.m)
-    lengths = u.Quantity([0.1, 0.2, 0.3, 0.4], u.m)
-    intensities = np.array([2000, 5000])
-
-    xs = u.Quantity([0.1, 0.2, -0.1, -0.2], u.m)
-    ys = u.Quantity([-0.2, -0.1, 0.2, 0.1], u.m)
-    psis = Angle([-60, -45, 10, 45, 60], unit="deg")
-
-    for x, y, width, length, intensity in zip(xs, ys, widths, lengths, intensities):
-        for psi in psis:
-            # Gaussian
-            image, clean_mask = create_sample_image(psi="0d", geometry=geom)
+    # Gaussian
+    image, clean_mask = create_sample_image(psi="0d", geometry=geom)
 
-            fit = image_fit_parameters(
-                geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("gaussian")
-            )
+    fit = image_fit_parameters(
+        geom, image, n_row=2, cleaned_mask=clean_mask, pdf=PDFType("gaussian")
+    )
 
-            cleaned_image = image.copy()
-            cleaned_image[~clean_mask] = 0.0
-            conc = concentration_parameters(geom, image, fit)
-            signal_inside_ellipse = conc.core
+    cleaned_image = image.copy()
+    cleaned_image[~clean_mask] = 0.0
+    conc = concentration_parameters(geom, image, fit)
+    signal_inside_ellipse = conc.core
 
-            if fit.is_valid and fit.is_accurate:
-                assert signal_inside_ellipse > 0.3
+    if fit.is_valid and fit.is_accurate:
+        assert signal_inside_ellipse > 0.3
 
 
 def test_with_toy_mst_tel(prod5_mst_flashcam):
@@ -499,6 +489,7 @@ def test_reconstruction_in_telescope_frame(prod5_lst):
     xs = u.Quantity([0.5, 0.5, -0.5, -0.5], u.m)
     ys = u.Quantity([0.5, -0.5, 0.5, -0.5], u.m)
     psis = Angle([-90, -45, 0, 45, 90], unit="deg")
+    skew = 0.5
 
     def distance(coord):
         return np.sqrt(np.diff(coord.x) ** 2 + np.diff(coord.y) ** 2) / 2
@@ -530,7 +521,9 @@ def get_transformed_width(telescope_hillas, telescope_frame, camera_frame):
     for x, y in zip(xs, ys):
         for psi in psis:
             # generate a toy image
-            model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi)
+            model = toymodel.SkewedGaussian(
+                x=x, y=y, width=width, length=length, psi=psi, skewness=skew
+            )
             image, signal, noise = model.generate_image(
                 geom, intensity=intensity, nsb_level_pe=5
             )