Matlab中的LBP算法源码和python中scikit-image库中的源码_vision.internal.inputvalidation.checkimagesize

参考文献：《Rotation Invariant Texture Classification using feature distributions》

《Multiresolution gray scale and rotation ivariant texture classification with local binary patterns》

matlab中的源码：

function features = extractLBPFeatures(I,varargin)
%extractLBPFeatures Extract LBP features.
%  features = extractLBPFeatures(I) extracts uniform local binary patterns
%  (LBP) from a grayscale image I and returns the features in
%  a 1-by-N vector. LBP features encode local texture information and can
%  be used for many tasks including classification, detection, and
%  recognition.
%
%  The LBP feature length, N, is based on the image size and the parameter
%  values listed below. See the <a href="matlab:helpview(fullfile(docroot,'toolbox','vision','vision.map'),'lbpFeatureLength')" >documentation</a> for more information.
%
%  [...] = extractLBPFeatures(..., Name, Value) specifies additional
%  name-value pairs described below. LBP algorithm parameters control how
%  local binary patterns are computed for each pixel in I. LBP histogram
%  parameters determine how the distribution of binary patterns is
%  aggregated over I to produce the output features.
%
%  LBP algorithm parameters
%  ------------------------
%
%  'NumNeighbors'  The number of neighbors used to compute the local binary
%                  pattern for each pixel in I. The set of neighbors is
%                  selected from a circularly symmetric pattern around each
%                  pixel. Increase the number of neighbors to encode
%                  greater detail around each pixel. Typical values are
%                  between 4 and 24.
%
%                  Default: 8
%
%  'Radius'        The radius, in pixels, of the circular pattern used to
%                  select neighbors for each pixel in I. Increase the
%                  radius to capture detail over a larger spatial scale.
%                  Typical values range from 1 to 5.
%
%                  Default: 0
%
%  'Upright'       A logical scalar. When set to true, the LBP features do
%                  not encode rotation information. Set 'Upright' to false
%                  when rotationally invariant features are required.
%
%                  Default: true
%
%  'Interpolation' Specify the interpolation method used to compute pixel
%                  neighbors as 'Nearest' or 'Linear'. Use 'Nearest' for
%                  faster computation at the cost of accuracy.
%
%                  Default: 'Linear'
%
%  LBP histogram parameters
%  ------------------------
%
%  'CellSize'      A 2-element vector that partitions I into
%                  floor(size(I)./CellSize) non-overlapping cells.
%                  Select larger cell sizes to collect information over
%                  larger regions at the cost of loosing local detail.
%
%                  Default: size(I)
%
%  'Normalization' Specify the type of normalization applied to the LBP
%                  histograms as 'L2' or 'None'. Select 'None' to apply a
%                  custom normalization method as a post-processing step.
%
%                  Default: 'L2'
%
% Class Support
% -------------
% The input image I can be uint8, uint16, int16, double, single, or
% logical, and it must be real and non-sparse.
%
% Notes
% -----
% This function extracts uniform local binary patterns. Uniform patterns
% have at most two 1-to-0 or 0-to-1 bit transitions.
%
% Example - Differentiate images by texture using LBP features.
% ---------------------------------------------------------------
% % Read images that contain different textures.
% brickWall = imread('bricks.jpg');
% rotatedBrickWall = imread('bricksRotated.jpg');
% carpet  = imread('carpet.jpg');
%
% figure
% imshow(brickWall)
% title('Bricks')
%
% figure
% imshow(rotatedBrickWall)
% title('Rotated bricks')
%
% figure
% imshow(carpet)
% title('Carpet')
%
% % Extract LBP features to encode image texture information.
% lbpBricks1 = extractLBPFeatures(brickWall,'Upright',false);
% lbpBricks2 = extractLBPFeatures(rotatedBrickWall,'Upright',false);
% lbpCarpet = extractLBPFeatures(carpet,'Upright',false);
%
% % Compute the squared error between the LBP features. This helps gauge
% % the similarity between the LBP features.
% brickVsBrick = (lbpBricks1 - lbpBricks2).^2;
% brickVsCarpet = (lbpBricks1 - lbpCarpet).^2;
%
% % Visualize the squared error to compare bricks vs. bricks and bricks vs.
% % carpet. The squared error is smaller when images have similar texture.
% figure
% bar([brickVsBrick; brickVsCarpet]', 'grouped')
% title('Squared error of LBP Histograms')
% xlabel('LBP Histogram Bins')
% legend('Bricks vs Rotated Bricks', 'Bricks vs Carpet')
%
% See also extractHOGFeatures, extractFeatures, detectHarrisFeatures,
% detectFASTFeatures, detectMinEigenFeatures, detectSURFFeatures,
% detectMSERFeatures, detectBRISKFeatures
 
% Copyright 2015 The MathWorks, Inc.
%
% References
% ----------
% Ojala, Timo, Matti Pietikainen, and Topi Maenpaa. "Multiresolution
% gray-scale and rotation invariant texture classification with local
% binary patterns." Pattern Analysis and Machine Intelligence, IEEE
% Transactions on 24.7 (2002): 971-987.
 
%#codegen
 
if isempty(coder.target)
    
    params = parseInputs(I,varargin{:});
    
    lbpImpl  = vision.internal.LBPImpl.getImpl(params);
    
    features = lbpImpl.extractLBPFeatures(I);
    
else
    
    [numNeighbors, radius, interpolation, uniform, upright, cellSize,...
        normalization] = codegenParseInputs(I,varargin{:});
    
    features = vision.internal.LBPImpl.codegenExtractLBPFeatures(... 
        I, numNeighbors, radius, interpolation, ...
        uniform, upright, cellSize, normalization);
end
 
% -------------------------------------------------------------------------
function params = parseInputs(I, varargin)
 
vision.internal.inputValidation.validateImage(I, 'I', 'grayscale');
 
szI = size(I);
 
parser = getInputParser();
parser.parse(varargin{:});
 
userInput = parser.Results;
 
usingDefaultCellSize = ismember('CellSize', parser.UsingDefaults);
    
if usingDefaultCellSize
    userInput.CellSize = szI;  % cell size default is size(I)
end
 
[validInterpolation, validNormalization] = validate(...
    userInput.NumNeighbors, userInput.Radius, userInput.CellSize, ...
    userInput.Upright, userInput.Interpolation, userInput.Normalization);
 
params = setParams(userInput, validInterpolation, validNormalization);
 
crossCheckParams(szI, params.CellSize, params.Radius)
 
% -------------------------------------------------------------------------
function [numNeighbors, radius, interpolation, uniform, upright, ...
    cellSize, normalization] = codegenParseInputs(I, varargin)
 
vision.internal.inputValidation.validateImage(I, 'I', 'grayscale');
eml_invariant(eml_is_const(ismatrix(I)), eml_message('vision:dims:imageNot2D'));
 
szI = size(I);
 
pvPairs = struct( ...
    'NumNeighbors',  uint32(0), ...
    'Radius',        uint32(0), ...    
    'CellSize',      uint32(0), ...   
    'Upright',       uint32(0),...
    'Interpolation', uint32(0),...
    'Normalization', uint32(0));
 
popt = struct( ...
    'CaseSensitivity', false, ...
    'StructExpand'   , true, ...
    'PartialMatching', true);
 
defaults = getParamDefaults();
 
optarg = eml_parse_parameter_inputs(pvPairs, popt, varargin{:});
 
usingDefaultCellSize = ~optarg.CellSize;
 
numNeighbors = eml_get_parameter_value(optarg.NumNeighbors, ...
    defaults.NumNeighbors, varargin{:});
 
radius = eml_get_parameter_value(optarg.Radius, ...
    defaults.Radius, varargin{:});
 
cellSize = eml_get_parameter_value(optarg.CellSize, ...
    defaults.CellSize, varargin{:});
 
upright = eml_get_parameter_value(optarg.Upright, ...
    defaults.Upright, varargin{:});
 
userInterpolation = eml_get_parameter_value(optarg.Interpolation, ...
    defaults.Interpolation, varargin{:});
 
userNormalization = eml_get_parameter_value(optarg.Normalization, ...
    defaults.Normalization, varargin{:});
 
% check const-ness before assigning to struct
vision.internal.errorIfNotConst(numNeighbors,      'NumNeighbors');
vision.internal.errorIfNotConst(radius,            'Radius');
vision.internal.errorIfNotConst(userInterpolation, 'Interpolation');
vision.internal.errorIfNotConst(userNormalization, 'Normalization');
vision.internal.errorIfNotConst(upright,           'Upright');
 
% check const-ness of size
vision.internal.errorIfNotFixedSize(numNeighbors, 'NumNeighbors');
vision.internal.errorIfNotFixedSize(radius,       'Radius');
vision.internal.errorIfNotFixedSize(cellSize,     'CellSize');
  
if usingDefaultCellSize
    cellSize = szI;  % cell size default is size(I)
end
 
[interpolation, normalization] = validate(numNeighbors, radius, cellSize, ...
    upright, userInterpolation, userNormalization);
 
numNeighbors  = single(numNeighbors);
radius        = single(radius);
cellSize      = single(cellSize);
upright       = logical(upright);
uniform       = true;
 
crossCheckParams(szI, cellSize, radius)
 
% -------------------------------------------------------------------------
function params = setParams(userInput, interpMethod, normMethod)
params.NumNeighbors  = single(userInput.NumNeighbors);
params.Radius        = single(userInput.Radius);
params.CellSize      = single(userInput.CellSize);
params.Upright       = logical(userInput.Upright);
params.Interpolation = interpMethod;
params.Normalization = normMethod;
params.Uniform       = true;
params.UseLUT        = false; % reset later based on other params
 
% -------------------------------------------------------------------------
function crossCheckParams(szI, cellSize, radius)
crossCheckImageSizeAndCellSize(szI, cellSize);
    
crossCheckImageSizeAndRadius(szI, radius);
    
crossCheckCellSizeAndRadius(cellSize, radius);
 
% -------------------------------------------------------------------------
function [validInterpolation, validNormalization] = validate(numNeighbors, ...
    radius, cellSize, upright, interpolation, normalization)
 
checkNumNeighbors(numNeighbors);
 
checkRadius(radius);
 
vision.internal.inputValidation.validateLogical(upright, 'Upright');
     
checkCellSize(cellSize);
 
validInterpolation = checkInterpolation(interpolation);
 
validNormalization = checkNormalization(normalization);
 
% -------------------------------------------------------------------------
function checkNumNeighbors(n)
 
vision.internal.errorIfNotFixedSize(n, 'NumNeighbors');
 
validateattributes(n, {'numeric'}, ...
    {'integer', 'real', 'nonsparse', 'scalar', '>=', 2' '<=' 32'},...
    mfilename, 'NumNeighbors'); %#ok<*EMCA>
 
% -------------------------------------------------------------------------
function checkRadius(r)
 
validateattributes(r, {'numeric'}, ...
    {'integer', 'real', 'nonsparse', 'scalar', '>=', 1},...
    mfilename, 'Radius');
 
% -------------------------------------------------------------------------
function str = checkInterpolation(method)
 
str = validatestring(method, {'Nearest', 'Linear'},...
    mfilename, 'Interpolation');
 
% -------------------------------------------------------------------------
function str = checkNormalization(method)
 
str = validatestring(method, {'L2', 'None'},...
    mfilename, 'Normalization');
 
% -------------------------------------------------------------------------
function checkCellSize(sz)
 
validateattributes(sz, {'numeric'}, ...
    {'vector', 'numel', 2, 'positive', 'real', 'integer', 'nonsparse'},...
    mfilename, 'CellSize');
 
% -------------------------------------------------------------------------
function crossCheckCellSizeAndRadius(sz, r)
coder.internal.errorIf(any(sz < (2*r + 1)),...
    'vision:extractLBPFeatures:cellSizeLTRadius');
 
% -------------------------------------------------------------------------
function crossCheckImageSizeAndCellSize(imgSize, cellSize)
 
coder.internal.errorIf(any(cellSize(:) > imgSize(:)), ...
    'vision:extractLBPFeatures:imgSizeLTCellSize');
 
% -------------------------------------------------------------------------
function crossCheckImageSizeAndRadius(sz, r)
coder.internal.errorIf(any(sz < (2*r + 1)),...
    'vision:extractLBPFeatures:imgSizeLTRadius');
 
% -------------------------------------------------------------------------
function parser = getInputParser()
persistent p; % cache parser for speed
 
if isempty(p)
    
    defaults = getParamDefaults();
    
    p = inputParser();
    addParameter(p, 'NumNeighbors',  defaults.NumNeighbors);
    addParameter(p, 'Radius',        defaults.Radius);
    addParameter(p, 'Upright',       defaults.Upright);
    addParameter(p, 'CellSize',      defaults.CellSize);
    addParameter(p, 'Normalization', defaults.Normalization);
    addParameter(p, 'Interpolation', defaults.Interpolation);
end
parser = p;
 
% -------------------------------------------------------------------------
function defaults = getParamDefaults()
 
defaults.NumNeighbors  = single(8);
defaults.Radius        = single(1);
defaults.Upright       = true;
defaults.CellSize      = [3 3];  % default is size(I), but give values here to define dims/type
defaults.Normalization = 'L2';
defaults.Interpolation = 'Linear';

python scikit-image库中的LBP算法源码：

#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as cnp
from libc.math cimport sin, cos, abs
from .._shared.interpolation cimport bilinear_interpolation, round
from .._shared.transform cimport integrate
 
cdef extern from "numpy/npy_math.h":
    double NAN "NPY_NAN"
 
from .._shared.fused_numerics cimport np_anyint as any_int
from .._shared.fused_numerics cimport np_real_numeric
 
 
def _glcm_loop(any_int[:, ::1] image, double[:] distances,
               double[:] angles, Py_ssize_t levels,
               cnp.uint32_t[:, :, :, ::1] out):
    """Perform co-occurrence matrix accumulation.
    Parameters
    ----------
    image : ndarray
        Integer typed input image. Only positive valued images are supported.
        If type is other than uint8, the argument `levels` needs to be set.
    distances : ndarray
        List of pixel pair distance offsets.
    angles : ndarray
        List of pixel pair angles in radians.
    levels : int
        The input image should contain integers in [0, `levels`-1],
        where levels indicate the number of gray-levels counted
        (typically 256 for an 8-bit image).
    out : ndarray
        On input a 4D array of zeros, and on output it contains
        the results of the GLCM computation.
    """
 
    cdef:
        Py_ssize_t a_idx, d_idx, r, c, rows, cols, row, col, start_row,\
                   end_row, start_col, end_col, offset_row, offset_col
        any_int i, j
        cnp.float64_t angle, distance
 
    with nogil:
        rows = image.shape[0]
        cols = image.shape[1]
 
        for a_idx in range(angles.shape[0]):
            angle = angles[a_idx]
            for d_idx in range(distances.shape[0]):
                distance = distances[d_idx]
                offset_row = round(sin(angle) * distance)
                offset_col = round(cos(angle) * distance)
                start_row = max(0, -offset_row)
                end_row = min(rows, rows - offset_row)
                start_col = max(0, -offset_col)
                end_col = min(cols, cols - offset_col)
                for r in range(start_row, end_row):
                    for c in range(start_col, end_col):
                        i = image[r, c]
                        # compute the location of the offset pixel
                        row = r + offset_row
                        col = c + offset_col
                        j = image[row, col]
                        if 0 <= i < levels and 0 <= j < levels:
                            out[i, j, d_idx, a_idx] += 1
 
 
cdef inline int _bit_rotate_right(int value, int length) nogil:
    """Cyclic bit shift to the right.
    Parameters
    ----------
    value : int
        integer value to shift
    length : int
        number of bits of integer
    """
    return (value >> 1) | ((value & 1) << (length - 1))
 
def _local_binary_pattern(double[:, ::1] image,
                          int P, float R, char method=b'D'):
    """Gray scale and rotation invariant LBP (Local Binary Patterns).
    LBP is an invariant descriptor that can be used for texture classification.
    Parameters
    ----------
    image : (N, M) double array
        Graylevel image.
    P : int
        Number of circularly symmetric neighbour set points (quantization of
        the angular space).
    R : float
        Radius of circle (spatial resolution of the operator).
    method : {'D', 'R', 'U', 'N', 'V'}
        Method to determine the pattern.
        * 'D': 'default'
        * 'R': 'ror'
        * 'U': 'uniform'
        * 'N': 'nri_uniform'
        * 'V': 'var'
    Returns
    -------
    output : (N, M) array
        LBP image.
    """
 
    # texture weights
    cdef int[::1] weights = 2 ** np.arange(P, dtype=np.int32)
    # local position of texture elements
    rr = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.double) / P)
    cc = R * np.cos(2 * np.pi * np.arange(P, dtype=np.double) / P)
    cdef double[::1] rp = np.round(rr, 5)
    cdef double[::1] cp = np.round(cc, 5)
 
    # pre-allocate arrays for computation
    cdef double[::1] texture = np.zeros(P, dtype=np.double)
    cdef signed char[::1] signed_texture = np.zeros(P, dtype=np.int8)
    cdef int[::1] rotation_chain = np.zeros(P, dtype=np.int32)
 
    output_shape = (image.shape[0], image.shape[1])
    cdef double[:, ::1] output = np.zeros(output_shape, dtype=np.double)
 
    cdef Py_ssize_t rows = image.shape[0]
    cdef Py_ssize_t cols = image.shape[1]
 
    cdef double lbp
    cdef Py_ssize_t r, c, changes, i
    cdef Py_ssize_t rot_index, n_ones
    cdef cnp.int8_t first_zero, first_one
 
    # To compute the variance features
    cdef double sum_, var_, texture_i
 
    with nogil:
        for r in range(image.shape[0]):
            for c in range(image.shape[1]):
                for i in range(P):
                    bilinear_interpolation[cnp.float64_t, double, double](
                            &image[0, 0], rows, cols, r + rp[i], c + cp[i],
                            b'C', 0, &texture[i])
                # signed / thresholded texture
                for i in range(P):
                    if texture[i] - image[r, c] >= 0:
                        signed_texture[i] = 1
                    else:
                        signed_texture[i] = 0
 
                lbp = 0
 
                # if method == b'var':
                if method == b'V':
                    # Compute the variance without passing from numpy.
                    # Following the LBP paper, we're taking a biased estimate
                    # of the variance (ddof=0)
                    sum_ = 0.0
                    var_ = 0.0
                    for i in range(P):
                        texture_i = texture[i]
                        sum_ += texture_i
                        var_ += texture_i * texture_i
                    var_ = (var_ - (sum_ * sum_) / P) / P
                    if var_ != 0:
                        lbp = var_
                    else:
                        lbp = NAN
                # if method == b'uniform':
                elif method == b'U' or method == b'N':
                    # determine number of 0 - 1 changes
                    changes = 0
                    for i in range(P - 1):
                        changes += (signed_texture[i]
                                    - signed_texture[i + 1]) != 0
                    if method == b'N':
                        # Uniform local binary patterns are defined as patterns
                        # with at most 2 value changes (from 0 to 1 or from 1 to
                        # 0). Uniform patterns can be characterized by their
                        # number `n_ones` of 1.  The possible values for
                        # `n_ones` range from 0 to P.
                        #
                        # Here is an example for P = 4:
                        # n_ones=0: 0000
                        # n_ones=1: 0001, 1000, 0100, 0010
                        # n_ones=2: 0011, 1001, 1100, 0110
                        # n_ones=3: 0111, 1011, 1101, 1110
                        # n_ones=4: 1111
                        #
                        # For a pattern of size P there are 2 constant patterns
                        # corresponding to n_ones=0 and n_ones=P. For each other
                        # value of `n_ones` , i.e n_ones=[1..P-1], there are P
                        # possible patterns which are related to each other
                        # through circular permutations. The total number of
                        # uniform patterns is thus (2 + P * (P - 1)).
                        # Given any pattern (uniform or not) we must be able to
                        # associate a unique code:
                        #
                        # 1. Constant patterns patterns (with n_ones=0 and
                        # n_ones=P) and non uniform patterns are given fixed
                        # code values.
                        #
                        # 2. Other uniform patterns are indexed considering the
                        # value of n_ones, and an index called 'rot_index'
                        # reprenting the number of circular right shifts
                        # required to obtain the pattern starting from a
                        # reference position (corresponding to all zeros stacked
                        # on the right). This number of rotations (or circular
                        # right shifts) 'rot_index' is efficiently computed by
                        # considering the positions of the first 1 and the first
                        # 0 found in the pattern.
                        if changes <= 2:
                            # We have a uniform pattern
                            n_ones = 0  # determines the number of ones
                            first_one = -1  # position was the first one
                            first_zero = -1  # position of the first zero
                            for i in range(P):
                                if signed_texture[i]:
                                    n_ones += 1
                                    if first_one == -1:
                                        first_one = i
                                else:
                                    if first_zero == -1:
                                        first_zero = i
                            if n_ones == 0:
                                lbp = 0
                            elif n_ones == P:
                                lbp = P * (P - 1) + 1
                            else:
                                if first_one == 0:
                                    rot_index = n_ones - first_zero
                                else:
                                    rot_index = P - first_one
                                lbp = 1 + (n_ones - 1) * P + rot_index
                        else:  # changes > 2
                            lbp = P * (P - 1) + 2
                    else:  # method != 'N'
                        if changes <= 2:
                            for i in range(P):
                                lbp += signed_texture[i]
                        else:
                            lbp = P + 1
                else:
                    # method == b'default'
                    for i in range(P):
                        lbp += signed_texture[i] * weights[i]
                    # method == b'ror'
                    if method == b'R':
                        # shift LBP P times to the right and get minimum value
                        rotation_chain[0] = <int>lbp
                        for i in range(1, P):
                            rotation_chain[i] = \
                                _bit_rotate_right(rotation_chain[i - 1], P)
                        lbp = rotation_chain[0]
                        for i in range(1, P):
                            lbp = min(lbp, rotation_chain[i])
                output[r, c] = lbp
    return np.asarray(output)
# Constant values that are used by `_multiblock_lbp` function.
# Values represent offsets of neighbour rectangles relative to central one.
# It has order starting from top left and going clockwise.
cdef:
    Py_ssize_t[::1] mlbp_r_offsets = np.asarray([-1, -1, -1, 0, 1, 1, 1, 0], dtype=np.intp)
    Py_ssize_t[::1] mlbp_c_offsets = np.asarray([-1, 0, 1, 1, 1, 0, -1, -1], dtype=np.intp)
cpdef int _multiblock_lbp(float[:, ::1] int_image,
                          Py_ssize_t r,
                          Py_ssize_t c,
                          Py_ssize_t width,
                          Py_ssize_t height) nogil:
    """Multi-block local binary pattern (MB-LBP) [1]_.
    Parameters
    ----------
    int_image : (N, M) float array
        Integral image.
    r : int
        Row-coordinate of top left corner of a rectangle containing feature.
    c : int
        Column-coordinate of top left corner of a rectangle containing feature.
    width : int
        Width of one of 9 equal rectangles that will be used to compute
        a feature.
    height : int
        Height of one of 9 equal rectangles that will be used to compute
        a feature.
    Returns
    -------
    output : int
        8-bit MB-LBP feature descriptor.
    References
    ----------
    .. [1] Face Detection Based on Multi-Block LBP
           Representation. Lun Zhang, Rufeng Chu, Shiming Xiang, Shengcai Liao,
           Stan Z. Li
           http://www.cbsr.ia.ac.cn/users/scliao/papers/Zhang-ICB07-MBLBP.pdf
    """
    cdef:
        # Top-left coordinates of central rectangle.
        Py_ssize_t central_rect_r = r + height
        Py_ssize_t central_rect_c = c + width
        Py_ssize_t r_shift = height - 1
        Py_ssize_t c_shift = width - 1
        Py_ssize_t current_rect_r, current_rect_c
        Py_ssize_t element_num, i
        double current_rect_val
        int has_greater_value
        int lbp_code = 0
    # Sum of intensity values of central rectangle.
    cdef float central_rect_val = integrate(int_image, central_rect_r, central_rect_c,
                                            central_rect_r + r_shift,
                                            central_rect_c + c_shift)
    for element_num in range(8):
        current_rect_r = central_rect_r + mlbp_r_offsets[element_num]*height
        current_rect_c = central_rect_c + mlbp_c_offsets[element_num]*width
        current_rect_val = integrate(int_image, current_rect_r, current_rect_c,
                                     current_rect_r + r_shift,
                                     current_rect_c + c_shift)
        has_greater_value = current_rect_val >= central_rect_val
        # If current rectangle's intensity value is bigger
        # make corresponding bit to 1.
        lbp_code |= has_greater_value << (7 - element_num)
 
    return lbp_code