Repository: scikit-image/scikit-image

Issue:
RANSAC returns an invalid model
## RANSAC returns an invalid model

In some cases that I cannot exactly pinpoint, RANSAC returns a model that does not pass the `is_model_valid` test.

## Way to reproduce
```python
import numpy as np
from skimage.measure import LineModelND, ransac

def is_model_valid(model, *random_data) -> bool:
    """Only allow models with a maximum of 10 degree tilt from the vertical"""
    tilt = abs(np.arccos(np.dot(model.params[1], [0, 0, 1])))
    return tilt <= (10 / 180 * np.pi)

np.random.seed(0)
model, inliers = ransac(
    np.linspace([0,0,0], [0.3,0,1], 1000) + np.random.rand(1000, 3) - 0.5,
    LineModelND,
    min_samples=2,
    residual_threshold=0.3,
    max_trials=50,
    random_state=0,
    is_model_valid=is_model_valid
)
print(is_model_valid(model) if model else "No model")
```


## Version information
```python
3.8.8 | packaged by conda-forge | (default, Feb 20 2021, 16:22:27) 
[GCC 9.3.0]
Linux-4.14.198-152.320.amzn2.x86_64-x86_64-with-glibc2.10
scikit-image version: 0.17.2
numpy version: 1.20.3
```

```python
print(is_model_valid(model) if model else "No model")
# False
```


Discussion:

Patch:
diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py
index 369e119cb..5dd3be62b 100644
--- a/skimage/data/__init__.py
+++ b/skimage/data/__init__.py
@@ -18,7 +18,6 @@ from .. import __version__
 
 import os.path as osp
 import os
-import stat
 
 __all__ = ['data_dir',
            'download_all',
diff --git a/skimage/feature/_cascade.pyx b/skimage/feature/_cascade.pyx
index 5fd1a03ba..255658473 100644
--- a/skimage/feature/_cascade.pyx
+++ b/skimage/feature/_cascade.pyx
@@ -424,13 +424,37 @@ cdef class Cascade:
     window_height : Py_ssize_t
         The height of a detection window.
     stages : Stage*
-        A link to the c array that stores stages information using
+        A pointer to the C array that stores stages information using a
         Stage struct.
     features : MBLBP*
-        Link to the c array that stores MBLBP features using MBLBP struct.
+        A pointer to the C array that stores MBLBP features using an MBLBP
+        struct.
     LUTs : cnp.uint32_t*
-        The ling to the array with look-up tables that are used by trained
+        A pointer to the C array with look-up tables that are used by trained
         MBLBP features (MBLBPStumps) to evaluate a particular region.
+
+    Notes
+    -----
+    The cascade approach was first described by Viola and Jones [1]_, [2]_,
+    although these initial publications used a set of Haar-like features. This
+    implementation instead uses multi-scale block local binary pattern (MB-LBP)
+    features [3]_.
+
+    References
+    ----------
+    .. [1] Viola, P. and Jones, M. "Rapid object detection using a boosted
+           cascade of simple features," In: Proceedings of the 2001 IEEE
+           Computer Society Conference on Computer Vision and Pattern
+           Recognition. CVPR 2001, pp. I-I.
+           :DOI:`10.1109/CVPR.2001.990517`
+    .. [2] Viola, P. and Jones, M.J, "Robust Real-Time Face Detection",
+           International Journal of Computer Vision 57, 137–154 (2004).
+           :DOI:`10.1023/B:VISI.0000013087.49260.fb`
+    .. [3] Liao, S. et al. Learning Multi-scale Block Local Binary Patterns for
+           Face Recognition. International Conference on Biometrics (ICB),
+           2007, pp. 828-837. In: Lecture Notes in Computer Science, vol 4642.
+           Springer, Berlin, Heidelberg.
+           :DOI:`10.1007/978-3-540-74549-5_87`
     """
 
     cdef:
diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx
index e9703e4af..9fa8157e5 100644
--- a/skimage/feature/_texture.pyx
+++ b/skimage/feature/_texture.pyx
@@ -310,10 +310,11 @@ cpdef int _multiblock_lbp(np_floats[:, ::1] int_image,
 
     References
     ----------
-    .. [1] Face Detection Based on Multi-Block LBP
-           Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao,
-           Stan Z. Li
+    .. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
+           on Multi-Block LBP Representation", In Proceedings: Advances in
+           Biometrics, International Conference, ICB 2007, Seoul, Korea.
            http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
+           :DOI:`10.1007/978-3-540-74549-5_2`
     """
 
     cdef:
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 08d3c2c61..c7324ae33 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -224,7 +224,7 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=0.5,
 
     Parameters
     ----------
-    image : 2D or 3D ndarray
+    image : ndarray
         Input grayscale image, blobs are assumed to be light on dark
         background (white on black).
     min_sigma : scalar or sequence of scalars, optional
@@ -268,7 +268,7 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=0.5,
     -------
     A : (n, image.ndim + sigma) ndarray
         A 2d array with each row representing 2 coordinate values for a 2D
-        image, and 3 coordinate values for a 3D image, plus the sigma(s) used.
+        image, or 3 coordinate values for a 3D image, plus the sigma(s) used.
         When a single sigma is passed, outputs are:
         ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or
         ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard
@@ -376,7 +376,7 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=0.5,
 
     # Catch no peaks
     if local_maxima.size == 0:
-        return np.empty((0, 3))
+        return np.empty((0, image.ndim + (1 if scalar_sigma else image.ndim)))
 
     # Convert local_maxima to float64
     lm = local_maxima.astype(float_dtype)
@@ -408,7 +408,7 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
 
     Parameters
     ----------
-    image : 2D or 3D ndarray
+    image : ndarray
         Input grayscale image, blobs are assumed to be light on dark
         background (white on black).
     min_sigma : scalar or sequence of scalars, optional
@@ -456,7 +456,7 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
     -------
     A : (n, image.ndim + sigma) ndarray
         A 2d array with each row representing 2 coordinate values for a 2D
-        image, and 3 coordinate values for a 3D image, plus the sigma(s) used.
+        image, or 3 coordinate values for a 3D image, plus the sigma(s) used.
         When a single sigma is passed, outputs are:
         ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or
         ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard
@@ -542,7 +542,7 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
 
     # Catch no peaks
     if local_maxima.size == 0:
-        return np.empty((0, 3))
+        return np.empty((0, image.ndim + (1 if scalar_sigma else image.ndim)))
 
     # Convert local_maxima to float64
     lm = local_maxima.astype(float_dtype)
diff --git a/skimage/feature/texture.py b/skimage/feature/texture.py
index c5495d7ca..a19602119 100644
--- a/skimage/feature/texture.py
+++ b/skimage/feature/texture.py
@@ -3,9 +3,10 @@ Methods to characterize image textures.
 """
 
 import numpy as np
+
 from .._shared.utils import check_nD
-from ..util import img_as_float
 from ..color import gray2rgb
+from ..util import img_as_float
 from ._texture import (_glcm_loop,
                        _local_binary_pattern,
                        _multiblock_lbp)
@@ -58,14 +59,15 @@ def graycomatrix(image, distances, angles, levels=None, symmetric=False,
 
     References
     ----------
-    .. [1] The GLCM Tutorial Home Page,
-           http://www.fp.ucalgary.ca/mhallbey/tutorial.htm
-    .. [2] Haralick, RM.; Shanmugam, K.,
-           "Textural features for image classification"
-           IEEE Transactions on systems, man, and cybernetics 6 (1973): 610-621.
+    .. [1] M. Hall-Beyer, 2007. GLCM Texture: A Tutorial
+           https://prism.ucalgary.ca/handle/1880/51900
+           DOI:`10.11575/PRISM/33280`
+    .. [2] R.M. Haralick, K. Shanmugam, and I. Dinstein, "Textural features for
+           image classification", IEEE Transactions on Systems, Man, and
+           Cybernetics, vol. SMC-3, no. 6, pp. 610-621, Nov. 1973.
            :DOI:`10.1109/TSMC.1973.4309314`
-    .. [3] Pattern Recognition Engineering, Morton Nadler & Eric P.
-           Smith
+    .. [3] M. Nadler and E.P. Smith, Pattern Recognition Engineering,
+           Wiley-Interscience, 1993.
     .. [4] Wikipedia, https://en.wikipedia.org/wiki/Co-occurrence_matrix
 
 
@@ -147,7 +149,7 @@ def graycomatrix(image, distances, angles, levels=None, symmetric=False,
     # normalize each GLCM
     if normed:
         P = P.astype(np.float64)
-        glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1))
+        glcm_sums = np.sum(P, axis=(0, 1), keepdims=True)
         glcm_sums[glcm_sums == 0] = 1
         P /= glcm_sums
 
@@ -170,8 +172,8 @@ def graycoprops(P, prop='contrast'):
         .. math:: \\sum_{i,j=0}^{levels-1} P_{i,j}\\left[\\frac{(i-\\mu_i) \\
                   (j-\\mu_j)}{\\sqrt{(\\sigma_i^2)(\\sigma_j^2)}}\\right]
 
-    Each GLCM is normalized to have a sum of 1 before the computation of texture
-    properties.
+    Each GLCM is normalized to have a sum of 1 before the computation of
+    texture properties.
 
     Parameters
     ----------
@@ -193,8 +195,10 @@ def graycoprops(P, prop='contrast'):
 
     References
     ----------
-    .. [1] The GLCM Tutorial Home Page,
-           http://www.fp.ucalgary.ca/mhallbey/tutorial.htm
+    .. [1] M. Hall-Beyer, 2007. GLCM Texture: A Tutorial v. 1.0 through 3.0.
+           The GLCM Tutorial Home Page,
+           https://prism.ucalgary.ca/handle/1880/51900
+           DOI:`10.11575/PRISM/33280`
 
     Examples
     --------
@@ -225,7 +229,7 @@ def graycoprops(P, prop='contrast'):
 
     # normalize each GLCM
     P = P.astype(np.float64)
-    glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1))
+    glcm_sums = np.sum(P, axis=(0, 1), keepdims=True)
     glcm_sums[glcm_sums == 0] = 1
     P /= glcm_sums
 
@@ -244,23 +248,20 @@ def graycoprops(P, prop='contrast'):
 
     # compute property for each GLCM
     if prop == 'energy':
-        asm = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
+        asm = np.sum(P ** 2, axis=(0, 1))
         results = np.sqrt(asm)
     elif prop == 'ASM':
-        results = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
+        results = np.sum(P ** 2, axis=(0, 1))
     elif prop == 'correlation':
         results = np.zeros((num_dist, num_angle), dtype=np.float64)
         I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
         J = np.array(range(num_level)).reshape((1, num_level, 1, 1))
-        diff_i = I - np.apply_over_axes(np.sum, (I * P), axes=(0, 1))[0, 0]
-        diff_j = J - np.apply_over_axes(np.sum, (J * P), axes=(0, 1))[0, 0]
+        diff_i = I - np.sum(I * P, axis=(0, 1))
+        diff_j = J - np.sum(J * P, axis=(0, 1))
 
-        std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i) ** 2),
-                                           axes=(0, 1))[0, 0])
-        std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j) ** 2),
-                                           axes=(0, 1))[0, 0])
-        cov = np.apply_over_axes(np.sum, (P * (diff_i * diff_j)),
-                                 axes=(0, 1))[0, 0]
+        std_i = np.sqrt(np.sum(P * (diff_i) ** 2, axis=(0, 1)))
+        std_j = np.sqrt(np.sum(P * (diff_j) ** 2, axis=(0, 1)))
+        cov = np.sum(P * (diff_i * diff_j), axis=(0, 1))
 
         # handle the special case of standard deviations near zero
         mask_0 = std_i < 1e-15
@@ -268,11 +269,11 @@ def graycoprops(P, prop='contrast'):
         results[mask_0] = 1
 
         # handle the standard case
-        mask_1 = mask_0 == False
+        mask_1 = ~mask_0
         results[mask_1] = cov[mask_1] / (std_i[mask_1] * std_j[mask_1])
     elif prop in ['contrast', 'dissimilarity', 'homogeneity']:
         weights = weights.reshape((num_level, num_level, 1, 1))
-        results = np.apply_over_axes(np.sum, (P * weights), axes=(0, 1))[0, 0]
+        results = np.sum(P * weights, axis=(0, 1))
 
     return results
 
@@ -302,7 +303,7 @@ def local_binary_pattern(image, P, R, method='default'):
             finer quantization of the angular space which is gray scale and
             rotation invariant.
         * 'nri_uniform': non rotation-invariant uniform patterns variant
-            which is only gray scale invariant [2]_.
+            which is only gray scale invariant [2]_, [3]_.
         * 'var': rotation invariant variance measures of the contrast of local
             image texture which is rotation but not gray scale invariant.
 
@@ -313,14 +314,21 @@ def local_binary_pattern(image, P, R, method='default'):
 
     References
     ----------
-    .. [1] Multiresolution Gray-Scale and Rotation Invariant Texture
-           Classification with Local Binary Patterns.
-           Timo Ojala, Matti Pietikainen, Topi Maenpaa.
-           http://www.ee.oulu.fi/research/mvmp/mvg/files/pdf/pdf_94.pdf, 2002.
-    .. [2] Face recognition with local binary patterns.
-           Timo Ahonen, Abdenour Hadid, Matti Pietikainen,
-           http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.214.6851,
-           2004.
+    .. [1] T. Ojala, M. Pietikainen, T. Maenpaa, "Multiresolution gray-scale
+           and rotation invariant texture classification with local binary
+           patterns", IEEE Transactions on Pattern Analysis and Machine
+           Intelligence, vol. 24, no. 7, pp. 971-987, July 2002
+           :DOI:`10.1109/TPAMI.2002.1017623`
+    .. [2] T. Ahonen, A. Hadid and M. Pietikainen. "Face recognition with
+           local binary patterns", in Proc. Eighth European Conf. Computer
+           Vision, Prague, Czech Republic, May 11-14, 2004, pp. 469-481, 2004.
+           http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.214.6851
+           :DOI:`10.1007/978-3-540-24670-1_36`
+    .. [3] T. Ahonen, A. Hadid and M. Pietikainen, "Face Description with
+           Local Binary Patterns: Application to Face Recognition",
+           IEEE Transactions on Pattern Analysis and Machine Intelligence,
+           vol. 28, no. 12, pp. 2037-2041, Dec. 2006
+           :DOI:`10.1109/TPAMI.2006.244.`
     """
     check_nD(image, 2)
 
@@ -371,10 +379,11 @@ def multiblock_lbp(int_image, r, c, width, height):
 
     References
     ----------
-    .. [1] Face Detection Based on Multi-Block LBP
-           Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao,
-           Stan Z. Li
+    .. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
+           on Multi-Block LBP Representation", In Proceedings: Advances in
+           Biometrics, International Conference, ICB 2007, Seoul, Korea.
            http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
+           :DOI:`10.1007/978-3-540-74549-5_2`
     """
 
     int_image = np.ascontiguousarray(int_image, dtype=np.float32)
@@ -431,10 +440,11 @@ def draw_multiblock_lbp(image, r, c, width, height,
 
     References
     ----------
-    .. [1] Face Detection Based on Multi-Block LBP
-           Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao,
-           Stan Z. Li
+    .. [1] L. Zhang, R. Chu, S. Xiang, S. Liao, S.Z. Li. "Face Detection Based
+           on Multi-Block LBP Representation", In Proceedings: Advances in
+           Biometrics, International Conference, ICB 2007, Seoul, Korea.
            http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
+           :DOI:`10.1007/978-3-540-74549-5_2`
     """
 
     # Default colors for regions.
diff --git a/skimage/filters/_gabor.py b/skimage/filters/_gabor.py
index 803144c29..9c0949cd1 100644
--- a/skimage/filters/_gabor.py
+++ b/skimage/filters/_gabor.py
@@ -1,3 +1,5 @@
+import math
+
 import numpy as np
 from scipy import ndimage as ndi
 
@@ -9,7 +11,7 @@ __all__ = ['gabor_kernel', 'gabor']
 def _sigma_prefactor(bandwidth):
     b = bandwidth
     # See http://www.cs.rug.nl/~imaging/simplecell.html
-    return 1.0 / np.pi * np.sqrt(np.log(2) / 2.0) * \
+    return 1.0 / np.pi * math.sqrt(math.log(2) / 2.0) * \
         (2.0 ** b + 1) / (2.0 ** b - 1)
 
 
@@ -83,19 +85,26 @@ def gabor_kernel(frequency, theta=0, bandwidth=1, sigma_x=None, sigma_y=None,
     if np.dtype(dtype).kind != 'c':
         raise ValueError("dtype must be complex")
 
-    x0 = np.ceil(max(np.abs(n_stds * sigma_x * np.cos(theta)),
-                     np.abs(n_stds * sigma_y * np.sin(theta)), 1))
-    y0 = np.ceil(max(np.abs(n_stds * sigma_y * np.cos(theta)),
-                     np.abs(n_stds * sigma_x * np.sin(theta)), 1))
-    y, x = np.mgrid[-y0:y0 + 1, -x0:x0 + 1]
-
-    rotx = x * np.cos(theta) + y * np.sin(theta)
-    roty = -x * np.sin(theta) + y * np.cos(theta)
-
-    g = np.zeros(y.shape, dtype=dtype)
-    g[:] = np.exp(-0.5 * (rotx ** 2 / sigma_x ** 2 + roty ** 2 / sigma_y ** 2))
-    g /= 2 * np.pi * sigma_x * sigma_y
-    g *= np.exp(1j * (2 * np.pi * frequency * rotx + offset))
+    ct = math.cos(theta)
+    st = math.sin(theta)
+    x0 = math.ceil(
+        max(abs(n_stds * sigma_x * ct), abs(n_stds * sigma_y * st), 1)
+    )
+    y0 = math.ceil(
+        max(abs(n_stds * sigma_y * ct), abs(n_stds * sigma_x * st), 1)
+    )
+    y, x = np.meshgrid(np.arange(-y0, y0 + 1),
+                       np.arange(-x0, x0 + 1),
+                       indexing='ij',
+                       sparse=True)
+    rotx = x * ct + y * st
+    roty = -x * st + y * ct
+
+    g = np.empty(roty.shape, dtype=dtype)
+    np.exp(-0.5 * (rotx ** 2 / sigma_x ** 2 + roty ** 2 / sigma_y ** 2)
+           + 1j * (2 * np.pi * frequency * rotx + offset),
+           out=g)
+    g *= 1 / (2 * np.pi * sigma_x * sigma_y)
 
     return g
 
diff --git a/skimage/measure/fit.py b/skimage/measure/fit.py
index 406dfeccb..4d8a53937 100644
--- a/skimage/measure/fit.py
+++ b/skimage/measure/fit.py
@@ -16,11 +16,6 @@ def _check_data_atleast_2D(data):
         raise ValueError('Input data must be at least 2D.')
 
 
-def _norm_along_axis(x, axis):
-    """NumPy < 1.8 does not support the `axis` argument for `np.linalg.norm`."""
-    return np.sqrt(np.einsum('ij,ij->i', x, x))
-
-
 class BaseModel(object):
 
     def __init__(self):
@@ -129,7 +124,7 @@ class LineModelND(BaseModel):
         origin, direction = params
         res = (data - origin) - \
               ((data - origin) @ direction)[..., np.newaxis] * direction
-        return _norm_along_axis(res, axis=1)
+        return np.linalg.norm(res, axis=1)
 
     def predict(self, x, axis=0, params=None):
         """Predict intersection of the estimated line model with a hyperplane
@@ -618,21 +613,14 @@ def _dynamic_max_trials(n_inliers, n_samples, min_samples, probability):
     if n_inliers == 0:
         return np.inf
 
-    nom = 1 - probability
-    if nom == 0:
+    if probability == 1:
         return np.inf
 
-    inlier_ratio = n_inliers / float(n_samples)
-    denom = 1 - inlier_ratio ** min_samples
-    if denom == 0:
+    if n_inliers == n_samples:
         return 1
-    elif denom == 1:
-        return np.inf
 
-    nom = np.log(nom)
-    denom = np.log(denom)
-    if denom == 0:
-        return 0
+    nom = math.log(1 - probability)
+    denom = math.log(1 - (n_inliers / n_samples) ** min_samples)
 
     return int(np.ceil(nom / denom))
 
@@ -779,8 +767,9 @@ def ransac(data, model_class, min_samples, residual_threshold,
     >>> sum(inliers) > 40
     True
 
-    RANSAC can be used to robustly estimate a geometric transformation. In this section,
-    we also show how to use a proportion of the total samples, rather than an absolute number.
+    RANSAC can be used to robustly estimate a geometric
+    transformation. In this section, we also show how to use a
+    proportion of the total samples, rather than an absolute number.
 
     >>> from skimage.transform import SimilarityTransform
     >>> rng = np.random.default_rng()
@@ -806,10 +795,11 @@ def ransac(data, model_class, min_samples, residual_threshold,
 
     """
 
-    best_model = None
     best_inlier_num = 0
     best_inlier_residuals_sum = np.inf
-    best_inliers = None
+    best_inliers = []
+    validate_model = is_model_valid is not None
+    validate_data = is_data_valid is not None
 
     random_state = np.random.default_rng(random_state)
 
@@ -819,7 +809,7 @@ def ransac(data, model_class, min_samples, residual_threshold,
     num_samples = len(data[0])
 
     if not (0 < min_samples < num_samples):
-        raise ValueError("`min_samples` must be in range (0, <number-of-samples>)")
+        raise ValueError(f"`min_samples` must be in range (0, {num_samples})")
 
     if residual_threshold < 0:
         raise ValueError("`residual_threshold` must be greater than zero")
@@ -831,74 +821,78 @@ def ransac(data, model_class, min_samples, residual_threshold,
         raise ValueError("`stop_probability` must be in range [0, 1]")
 
     if initial_inliers is not None and len(initial_inliers) != num_samples:
-        raise ValueError("RANSAC received a vector of initial inliers (length %i)"
-                         " that didn't match the number of samples (%i)."
-                         " The vector of initial inliers should have the same length"
-                         " as the number of samples and contain only True (this sample"
-                         " is an initial inlier) and False (this one isn't) values."
-                         % (len(initial_inliers), num_samples))
+        raise ValueError(
+            f"RANSAC received a vector of initial inliers (length "
+            f"{len(initial_inliers)}) that didn't match the number of "
+            f"samples ({num_samples}). The vector of initial inliers should "
+            f"have the same length as the number of samples and contain only "
+            f"True (this sample is an initial inlier) and False (this one "
+            f"isn't) values.")
 
     # for the first run use initial guess of inliers
     spl_idxs = (initial_inliers if initial_inliers is not None
-                else random_state.choice(num_samples, min_samples, replace=False))
+                else random_state.choice(num_samples, min_samples,
+                                         replace=False))
+
+    # estimate model for current random sample set
+    model = model_class()
 
     for num_trials in range(max_trials):
         # do sample selection according data pairs
         samples = [d[spl_idxs] for d in data]
-        # for next iteration choose random sample set and be sure that no samples repeat
+
+        # for next iteration choose random sample set and be sure that
+        # no samples repeat
         spl_idxs = random_state.choice(num_samples, min_samples, replace=False)
 
         # optional check if random sample set is valid
-        if is_data_valid is not None and not is_data_valid(*samples):
+        if validate_data and not is_data_valid(*samples):
             continue
 
-        # estimate model for current random sample set
-        sample_model = model_class()
-
-        success = sample_model.estimate(*samples)
+        success = model.estimate(*samples)
         # backwards compatibility
         if success is not None and not success:
             continue
 
         # optional check if estimated model is valid
-        if is_model_valid is not None and not is_model_valid(sample_model, *samples):
+        if validate_model and not is_model_valid(model, *samples):
             continue
 
-        sample_model_residuals = np.abs(sample_model.residuals(*data))
+        residuals = np.abs(model.residuals(*data))
         # consensus set / inliers
-        sample_model_inliers = sample_model_residuals < residual_threshold
-        sample_model_residuals_sum = np.sum(sample_model_residuals ** 2)
+        inliers = residuals < residual_threshold
+        residuals_sum = residuals.dot(residuals)
 
         # choose as new best model if number of inliers is maximal
-        sample_inlier_num = np.sum(sample_model_inliers)
+        inliers_count = np.count_nonzero(inliers)
         if (
             # more inliers
-            sample_inlier_num > best_inlier_num
+            inliers_count > best_inlier_num
             # same number of inliers but less "error" in terms of residuals
-            or (sample_inlier_num == best_inlier_num
-                and sample_model_residuals_sum < best_inlier_residuals_sum)
-        ):
-            best_model = sample_model
-            best_inlier_num = sample_inlier_num
-            best_inlier_residuals_sum = sample_model_residuals_sum
-            best_inliers = sample_model_inliers
+            or (inliers_count == best_inlier_num
+                and residuals_sum < best_inlier_residuals_sum)):
+            best_inlier_num = inliers_count
+            best_inlier_residuals_sum = residuals_sum
+            best_inliers = inliers
             dynamic_max_trials = _dynamic_max_trials(best_inlier_num,
                                                      num_samples,
                                                      min_samples,
                                                      stop_probability)
             if (best_inlier_num >= stop_sample_num
-                or best_inlier_residuals_sum <= stop_residuals_sum
-                or num_trials >= dynamic_max_trials):
+                    or best_inlier_residuals_sum <= stop_residuals_sum
+                    or num_trials >= dynamic_max_trials):
                 break
 
     # estimate final model using all inliers
-    if best_inliers is not None and any(best_inliers):
+    if any(best_inliers):
         # select inliers for each data array
         data_inliers = [d[best_inliers] for d in data]
-        best_model.estimate(*data_inliers)
+        model.estimate(*data_inliers)
+        if validate_model and not is_model_valid(model, *data_inliers):
+            warn("Estimated model is not valid. Try increasing max_trials.")
     else:
-        best_model = None
+        model = None
         best_inliers = None
         warn("No inliers found. Model not fitted")
 
-    return best_model, best_inliers
+    return model, best_inliers


Test patch:
diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
index f1dbba756..291207ab6 100644
--- a/skimage/feature/tests/test_blob.py
+++ b/skimage/feature/tests/test_blob.py
@@ -4,6 +4,7 @@ import numpy as np
 import pytest
 from numpy.testing import assert_almost_equal
 
+from skimage import feature
 from skimage.draw import disk
 from skimage.draw.draw3d import ellipsoid
 from skimage.feature import blob_dog, blob_doh, blob_log
@@ -169,6 +170,23 @@ def test_blob_dog_excl_border():
     assert blobs.shape[0] == 0, msg
 
 
+@pytest.mark.parametrize('anisotropic', [False, True])
+@pytest.mark.parametrize('ndim', [1, 2, 3, 4])
+@pytest.mark.parametrize('function_name', ['blob_dog', 'blob_log'])
+def test_nd_blob_no_peaks_shape(function_name, ndim, anisotropic):
+    # uniform image so no blobs will be found
+    z = np.zeros((16,) * ndim, dtype=np.float32)
+    if anisotropic:
+        max_sigma = 8 + np.arange(ndim)
+    else:
+        max_sigma = 8
+    blob_func = getattr(feature, function_name)
+    blobs = blob_func(z, max_sigma=max_sigma)
+    # z.ndim +  (z.ndim sigmas if anisotropic, only one sigma otherwise)
+    expected_shape = 2 * z.ndim if anisotropic else z.ndim + 1
+    assert blobs.shape == (0, expected_shape)
+
+
 @pytest.mark.parametrize(
     'dtype', [np.uint8, np.float16, np.float32, np.float64]
 )
diff --git a/skimage/measure/tests/test_fit.py b/skimage/measure/tests/test_fit.py
index eddecd354..e1149b1ba 100644
--- a/skimage/measure/tests/test_fit.py
+++ b/skimage/measure/tests/test_fit.py
@@ -426,3 +426,20 @@ def test_ransac_with_no_final_inliers():
                                 residual_threshold=0, random_state=1523427)
     assert inliers is None
     assert model is None
+
+
+def test_ransac_non_valid_best_model():
+    """Example from GH issue #5572"""
+    def is_model_valid(model, *random_data) -> bool:
+        """Allow models with a maximum of 10 degree tilt from the vertical
+
+        """
+        tilt = abs(np.arccos(np.dot(model.params[1], [0, 0, 1])))
+        return tilt <= (10 / 180 * np.pi)
+
+    rnd = np.random.RandomState(1)
+    data = np.linspace([0, 0, 0], [0.3, 0, 1], 1000) + rnd.rand(1000, 3) - 0.5
+    with expected_warnings(["Estimated model is not valid"]):
+        ransac(data, LineModelND, min_samples=2,
+               residual_threshold=0.3, max_trials=50, random_state=0,
+               is_model_valid=is_model_valid)