diff --git a/drizzle/resample.py b/drizzle/resample.py
index 18f4258..66d461f 100644
--- a/drizzle/resample.py
+++ b/drizzle/resample.py
@@ -105,6 +105,8 @@ def __init__(self, n_images=1, kernel="square",
                 "'out_shape' cannot be None when all output arrays are None."
             )
 
+        # shapes will collect user specified 'out_shape' and shapes of
+        # out_* arrays (if provided) in order to check all shapes are the same.
         shapes = set()
 
         if out_sci is not None:
@@ -128,7 +130,7 @@ def __init__(self, n_images=1, kernel="square",
                     "'out_shape' cannot be None when all output arrays are "
                     "None."
                 )
-            out_shape = shapes[0]
+            out_shape = shapes.pop()
         else:
             shapes.add(out_shape)
 
diff --git a/drizzle/utils.py b/drizzle/utils.py
index 4748689..85174c8 100644
--- a/drizzle/utils.py
+++ b/drizzle/utils.py
@@ -24,6 +24,15 @@ def calc_pixmap(wcs_from, wcs_to, estimate_pixel_scale_ratio=False,
     estimate_pixel_scale_ratio : bool, optional
         Estimate the ratio of "to" to "from" WCS pixel scales.
 
+        .. note::
+            Pixel scale is estimated as the square root of pixel's area
+            (i.e., pixels are assumed to have a square shape) at the reference
+            pixel position which is taken as the center of the bounding box
+            if ``wcs_*`` has a bounding box defined, or as the center of the box
+            defined by the ``pixel_shape`` attribute of the input WCS if
+            ``pixel_shape`` is defined (not `None`), or at pixel coordinates
+            ``(0, 0)``.
+
     refpix_from : numpy.ndarray, tuple, list
         Image coordinates of the reference pixel near which pixel scale should
         be computed in the "from" image. In FITS WCS this could be, for example,
@@ -72,7 +81,10 @@ def _estimate_pixel_scale(wcs, refpix):
     # from https://trs.jpl.nasa.gov/handle/2014/40409
     if refpix is None:
         if wcs.bounding_box is None:
-            refpix = np.ones(wcs.pixel_n_dim)
+            if wcs.pixel_shape:
+                refpix = np.array([(i - 1) // 2 for i in wcs.pixel_shape])
+            else:
+                refpix = np.zeros(wcs.pixel_n_dim)
         else:
             refpix = np.mean(wcs.bounding_box, axis=-1)
     else: