main.pyx

%%cython -a
# cython: cdivision=True
# cython: boundscheck=False
# cython: nonecheck=False
# cython: wraparound=False
# distutils: language = c++

from libcpp.vector cimport vector
from libcpp.pair cimport pair
import skimage.io as io
from skimage.feature import draw_multiblock_lbp
from skimage.feature._texture cimport _multiblock_lbp
from matplotlib import pyplot as plt
import os
import numpy as np

cdef struct MBLBP:

    Py_ssize_t r
    Py_ssize_t c
    Py_ssize_t width
    Py_ssize_t height

cdef struct TreeBranch:
    
    float numerator
    float denominator
    bint is_updated

cdef struct Tree:
    
    TreeBranch[256] branches

cdef vector[float] initialize_weights(int train_data_size):
    
    cdef vector[float] weights
    cdef float initial_weight = 1.0 / train_data_size 
    
    weights.resize(train_data_size, initial_weight)
    
    return weights

cdef vector[float] create_train_exmpl_index(int pos_exmpl_number, int neg_exmpl_number):
    
    cdef:
        vector[float] exmpl_indx
        int overall_number = pos_exmpl_number + neg_exmpl_number
    
    
    exmpl_indx.resize(pos_exmpl_number, 1)
    exmpl_indx.resize(overall_number, -1)
    
    return exmpl_indx
    

cdef pair[float[:, :, ::1], vector[float]] load_train_dataset(pos_exmpls, neg_exmpls):
    
    cdef:
        pair[float[:, :, ::1], vector[float]] result_container
        float[:, :, ::1] integral_img_cube_view
        vector[float] index_arr
        int pos_exmpl_number
        int neg_exmpl_number
    
    img_collection_pos = io.imread_collection(pos_exmpls)
    img_collection_neg = io.imread_collection(neg_exmpls)

    img_pos = io.concatenate_images(img_collection_pos)
    img_neg = io.concatenate_images(img_collection_neg)
    
    pos_exmpl_number = img_pos.shape[0]
    neg_exmpl_number = img_neg.shape[0]
    index_arr = create_train_exmpl_index(pos_exmpl_number, neg_exmpl_number)
    
    img_cube = np.vstack((img_pos, img_neg))
    integral_img_cube = img_cube.cumsum(1).cumsum(2)
    integral_img_cube_view = np.ascontiguousarray(integral_img_cube, dtype=np.float32)
    
    result_container.first = integral_img_cube_view
    result_container.second = index_arr
    
    return result_container

cdef vector[float[:, ::1]] create_integral_images_slices(float[:, :, ::1] int_img_cube):
    
    cdef:
        Py_ssize_t img_amount = int_img_cube.shape[0]
        Py_ssize_t img_number
        vector[float[:, ::1]] slices_container
    
    slices_container.resize(img_amount)
    
    for img_number in range(img_amount):
        slices_container[img_number] = int_img_cube[img_number, :, :]
    
    return slices_container

cdef Py_ssize_t mblbp_features_overall_amount(Py_ssize_t window_width, Py_ssize_t window_height):
    
    cdef:
        Py_ssize_t current_width
        Py_ssize_t current_height
        Py_ssize_t count = 0
    
    for current_width in range(1, window_width / 3 + 1):
        for current_height in range(1, window_height / 3 + 1):
            count += (window_width - current_width*3 + 1) * (window_height - current_height*3 + 1)
    
    return count

cdef vector[MBLBP] precompute_features(Py_ssize_t window_width, Py_ssize_t window_height):
    
    cdef:
        Py_ssize_t current_height
        Py_ssize_t current_width
        Py_ssize_t current_row
        Py_ssize_t current_col
        Py_ssize_t count = 0
        Py_ssize_t overall_amount = mblbp_features_overall_amount(window_width, window_height) 
        vector[MBLBP] features
        MBLBP new_feature

    features.resize(overall_amount)

    for current_width in range(1, window_width / 3 + 1):
        for current_height in range(1, window_height / 3 + 1):
            for current_row in range(window_height - current_height * 3 + 1):
                for current_col in range(window_width - current_width * 3 + 1):

                    new_feature = MBLBP(current_row, current_col, current_width, current_height)
                    features[count] = new_feature
                    count += 1
    
    return features


cdef vector[vector[int]] compute_feature_values(vector[float[:, ::1]] &int_img_slices, vector[MBLBP] &features):
    
    cdef:
        Py_ssize_t img_amount = int_img_slices.size()
        Py_ssize_t features_amount = features.size()
        Py_ssize_t img_count
        Py_ssize_t feature_count
        vector[vector[int]] feature_values
        float[:, ::1] current_slice
        MBLBP current_feature
    
    feature_values.resize(img_amount)

    for img_count in range(img_amount):
        feature_values[img_count].resize(features_amount)
    
    for img_count in range(img_amount):
        for feature_count in range(features_amount):
            
            current_slice = int_img_slices[img_count]
            current_feature = features[feature_count]
            
            feature_values[img_count][feature_count] = _multiblock_lbp(current_slice,
                                                                       current_feature.r,
                                                                       current_feature.c,
                                                                       current_feature.width,
                                                                       current_feature.height)
    
    return feature_values

cdef void normalize_tree(Tree &tree):
    
    cdef:
        Py_ssize_t current_branch_number
        float current_numerator
    
    for current_branch_number in range(256):
        if tree.branches[current_branch_number].is_updated:
            current_numerator = tree.branches[current_branch_number].numerator
            current_numerator /= tree.branches[current_branch_number].denominator
            tree.branches[current_branch_number].numerator = current_numerator
            tree.branches[current_branch_number].denominator = 1
            

cdef void normalize_trees(vector[Tree] &trees):
    
    cdef:
        Py_ssize_t trees_number = trees.size()
        Py_ssize_t current_tree_number 
    
    for current_tree_number in range(trees_number):
        normalize_tree(trees[current_tree_number])

cdef vector[Tree] fit_trees(vector[float] &weights, vector[vector[int]] &feature_values, vector[float] &pos_neg_data_indx):
    
    cdef:
        vector[Tree] trees
        Py_ssize_t img_number
        Py_ssize_t feature_number
        int current_branch_number
        float current_numerator
        float current_denominator
        Py_ssize_t data_size = weights.size()
        Py_ssize_t features_amount = 0
    
    if data_size > 0:
        features_amount = feature_values[0].size()
    
    trees.resize(features_amount)
    
    for img_number in range(data_size):
        for feature_number in range(features_amount):
        
            current_branch_number = feature_values[img_number][feature_number]
            current_numerator = weights[img_number] * pos_neg_data_indx[img_number]
            current_denominator = weights[img_number]
            trees[feature_number].branches[current_branch_number].numerator += current_numerator
            trees[feature_number].branches[current_branch_number].denominator += current_denominator
            trees[feature_number].branches[current_branch_number].is_updated = True
    
    normalize_trees(trees)
    
    return trees

cdef vector[float] weighted_squared_error(vector[float] &pos_neg_data_indx, vector[float] &weights, 
                                          vector[Tree] &trees, vector[vector[int]] &feature_values):
    
    cdef:
        Py_ssize_t features_amount = trees.size()
        Py_ssize_t data_size = pos_neg_data_indx.size()
        Py_ssize_t feature_number
        Py_ssize_t img_number
        int current_feature_value
        float current_branch_value
        vector[float] error
    
    error.resize(features_amount)
    
    for feature_number in range(features_amount):
        for img_number in range(data_size):
            
            current_feature_value = feature_values[img_number][feature_number]
            current_branch_value = trees[feature_number].branches[current_feature_value].numerator
            error[feature_number] += weights[img_number]*(pos_neg_data_indx[img_number] - current_branch_value)**2
            
    return error

cdef Py_ssize_t find_min_error_feature_index(vector[float] &error):
    
    cdef:
        Py_ssize_t features_amount = error.size()
        Py_ssize_t current_feature
        Py_ssize_t min_error_feature_index
        float min_error
        float current_error
    
    if features_amount:
        min_error_feature_index = 0
        min_error = error[0]
    
    for current_feature in range(features_amount):
        
        current_error = error[current_feature]
        
        if current_error < min_error:
            
            min_error = current_error
            min_error_feature_index = current_feature
    
    return min_error_feature_index
            
        
cdef:
    pair[float[:, :, ::1], vector[float]] output = load_train_dataset('24_24_faces_web/*.jpg', '24_24_nonfaces_aflw/*.jpg')
    float[:, :, ::1] int_img_cube = output.first
    vector[float] pos_neg_data_indx = output.second
    
    int data_size = int_img_cube.shape[0]
    vector[float] weights = initialize_weights(data_size)
    
    Py_ssize_t height = int_img_cube.shape[1]
    Py_ssize_t width = int_img_cube.shape[2]
    vector[MBLBP] features = precompute_features(width, height)
    
    Py_ssize_t features_amount = features.size()
    vector[float[:, ::1]] int_img_slices = create_integral_images_slices(int_img_cube)
    vector[vector[int]] feature_values = compute_feature_values(int_img_slices, features)
    vector[Tree] trees = fit_trees(weights, feature_values, pos_neg_data_indx)
    vector[float] errors = weighted_squared_error(pos_neg_data_indx, weights, trees, feature_values)
    Py_ssize_t min_error_feature_index = find_min_error_feature_index(errors)
    
    
print min_error_feature_index
print features[min_error_feature_index]