import numpy as np
from scipy.ndimage import convolve, maximum_filter


def derivative_filters():
    """ Create derivative filters for x and y direction

    Returns:
        fx: derivative filter in x direction
        fy: derivative filter in y direction
    """
    fx = np.array([[0.5, 0, -0.5]])
    fy = fx.transpose()
    return fx, fy


def gauss_2d(sigma, fsize):
    """ Create a 2D Gaussian filter

    Args:
        sigma: width of the Gaussian filter
        fsize: (h, w) dimensions of the filter
    Returns:
        *normalized* Gaussian filter as (h, w) np.array
    """
    m, n = fsize
    x = np.arange(-m / 2 + 0.5, m / 2)
    y = np.arange(-n / 2 + 0.5, n / 2)
    xx, yy = np.meshgrid(x, y, sparse=True)
    g = np.exp(-(xx ** 2 + yy ** 2) / (2 * sigma ** 2))
    return g / np.sum(g)


def compute_hessian(img, gauss, fx, fy):
    """ Compute elements of the Hessian matrix

    Args:
        img: numpy array with the image
        gauss: Gaussian filter
        fx: derivative filter in x direction
        fy: derivative filter in y direction

    Returns:
        I_xx: (h, w) np.array of 2nd derivatives in x direction
        I_yy: (h, w) np.array of 2nd derivatives in y direction
        I_xy: (h, w) np.array of 2nd derivatives in x-y direction
    """
    # Apply filters to the image
    img_smooth = convolve(img, gauss, mode='mirror')
    img_smooth_x = convolve(img_smooth, fx, mode='mirror')
    img_smooth_y = convolve(img_smooth, fy, mode='mirror')

    # Compute second derivatives
    I_xx = convolve(img_smooth_x, fx, mode='mirror')
    I_yy = convolve(img_smooth_y, fy, mode='mirror')
    I_xy = convolve(img_smooth_x, fy, mode='mirror')

    return I_xx, I_yy, I_xy


def compute_criterion(I_xx, I_yy, I_xy, sigma):
    """ Compute criterion function

    Args:
        I_xx: (h, w) np.array of 2nd derivatives in x direction
        I_yy: (h, w) np.array of 2nd derivatives in y direction
        I_xy: (h, w) np.array of 2nd derivatives in x-y direction
        sigma: scaling factor

    Returns:
        criterion: (h, w) np.array of scaled determinant of Hessian matrix
    """
    # Compute determinant of the Hessian matrix
    det_Hessian = (I_xx * I_yy) - (I_xy ** 2)

    # Create a 2D Gaussian filter
    size = int(6 * sigma)
    if size % 2 == 0:
        size += 1
    x = np.arange(-size // 2 + 1., size // 2 + 1.)
    gauss = np.exp(-0.5 * (x / sigma) ** 2)
    gauss /= gauss.sum()

    # Apply Gaussian filter to the determinant for smoothing
    criterion = convolve(det_Hessian, np.outer(gauss, gauss), mode='mirror')

    # Compute the scaled determinant
    criterion /= (sigma ** 4)

    return criterion


def non_max_suppression(criterion, threshold):
    """ Apply non-maximum suppression to criterion values
        and return Hessian interest points

        Args:
            criterion: (h, w) np.array of criterion function values
            threshold: criterion threshold
        Returns:
            rows: (n,) np.array with y-positions of interest points
            cols: (n,) np.array with x-positions of interest points
    """
    # Apply non-maximum suppression using convolution and maximum_filter
    local_maxima = (criterion == maximum_filter(criterion, footprint=np.ones((3, 3)), mode='constant'))

    # Apply threshold
    local_maxima[criterion < threshold] = False

    # Extract rows and columns of interest points
    rows, cols = np.where(local_maxima)

    return rows, cols
correct the code to get the intrest points