请写一个完整的MATLAB的ORB算法，不使用工具箱

ORB（Oriented FAST and Rotated BRIEF）算法是一种用于图像特征点描述子提取的方法，它在FAST算法和BRIEF算法的基础上进行了改进，可以实现旋转不变性和尺度不变性。下面是一个简单的MATLAB实现，主要分为四个部分：关键点检测、方向分配、描述子计算和匹配。

1. 关键点检测

首先，我们需要使用FAST算法来检测图像中的关键点。FAST算法是一种高效的角点检测算法，它能够快速地检测出图像中的角点，并且它的计算量比较小。

以下是 FAST 算法的 MATLAB 实现：

function [keypoints] = detectFAST(image, threshold)
% 设置 FAST 算法中的阈值
if nargin < 2
    threshold = 10;
end

% 在图像中搜索特征点
corners = detectMinEigenFeatures(image);

% 选择符合条件的角点
scores = corners.Metric;
strong_corners = corners.selectStrongest(length(corners), scores > threshold);

% 将角点的坐标存储到 keypoints 数组中
keypoints = strong_corners.Location;
end

2. 方向分配

在ORB算法中，我们需要为每个关键点计算一个方向，以保证描述子的旋转不变性。我们可以使用图像的梯度信息来计算每个关键点的方向。

以下是方向分配的 MATLAB 实现：

function [keypoints, orientations] = assignOrientation(image, keypoints, patch_size)
% 计算图像的梯度
[dx, dy] = gradient(image);

% 对角点周围的像素计算方向直方图
orientations = zeros(size(keypoints, 1), 36);
for i = 1:size(keypoints, 1)
    x = round(keypoints(i, 1));
    y = round(keypoints(i, 2));

    % 提取关键点周围的图像块
    patch = image(y-patch_size:y+patch_size, x-patch_size:x+patch_size);

    % 计算图像块中每个像素的梯度方向
    patch_dx = dx(y-patch_size:y+patch_size, x-patch_size:x+patch_size);
    patch_dy = dy(y-patch_size:y+patch_size, x-patch_size:x+patch_size);
    angles = atan2(patch_dy, patch_dx);
    magnitudes = sqrt(patch_dx.^2 + patch_dy.^2);

    % 将梯度方向加入方向直方图中
    for j = 1:numel(angles)
        bin = floor((angles(j) + pi) / (2*pi/36)) + 1;
        orientations(i, bin) = orientations(i, bin) + magnitudes(j);
    end
end

% 找到每个关键点的主要方向
[~, max_orientations] = max(orientations, [], 2);
orientations = max_orientations * (2*pi/36) - pi;

% 将关键点的方向存储到 keypoints 数组中
keypoints = [keypoints, orientations];
end

3. 描述子计算

在ORB算法中，我们使用BRIEF算法来计算每个关键点的描述子。BRIEF算法会比较关键点周围的像素对，并生成一个二进制字符串作为描述子。为了保证描述子的旋转不变性，我们需要对关键点的周围像素进行旋转，并使用旋转后的像素对生成描述子。

以下是描述子计算的 MATLAB 实现：

function [descriptors] = computeDescriptors(image, keypoints, patch_size, num_bits)
% 提取关键点周围的像素对，并旋转坐标系
pairs = generatePairs(num_bits);
descriptors = zeros(size(keypoints, 1), num_bits);
for i = 1:size(keypoints, 1)
    x = round(keypoints(i, 1));
    y = round(keypoints(i, 2));
    orientation = keypoints(i, 3);

    % 对关键点周围的像素对进行旋转
    rotated_pairs = rotatePairs(pairs, orientation);

    % 提取关键点周围的图像块
    patch = image(y-patch_size:y+patch_size, x-patch_size:x+patch_size);

    % 对关键点周围的像素对计算描述子
    for j = 1:num_bits
        pair = rotated_pairs(j, :);
        x1 = x + pair(1);
        y1 = y + pair(2);
        x2 = x + pair(3);
        y2 = y + pair(4);
        value1 = patch(y1-patch_size, x1-patch_size);
        value2 = patch(y2-patch_size, x2-patch_size);
        descriptors(i, j) = value1 < value2;
    end
end
end

function [pairs] = generatePairs(num_bits)
% 生成关键点周围的像素对
pairs = randi([-8, 8], num_bits, 4);
end

function [rotated_pairs] = rotatePairs(pairs, orientation)
% 旋转坐标系，使得关键点方向为水平方向
theta = -orientation;
rotated_pairs = pairs * [cos(theta), -sin(theta); sin(theta), cos(theta)];
end

4. 匹配

最后，我们需要使用描述子来进行匹配，找到两幅图像中的相似关键点。我们可以使用简单的欧几里得距离来计算描述子之间的差异。

以下是匹配的 MATLAB 实现：

function [matches] = matchDescriptors(descriptors1, descriptors2, threshold)
% 对两幅图像中的描述子进行匹配
distances = pdist2(descriptors1, descriptors2, 'hamming');
[~, indices] = sort(distances, 2);
matches = [];
for i = 1:size(indices, 1)
    if distances(i, indices(i, 1)) < threshold * distances(i, indices(i, 2))
        matches = [matches; i, indices(i, 1)];
    end
end
end

完整的ORB算法实现如下：

function [matches] = ORB(image1, image2, threshold)
% 设置算法参数
if nargin < 3
    threshold = 0.8;
end
patch_size = 16;
num_bits = 256;

% 检测关键点
keypoints1 = detectFAST(image1);
keypoints2 = detectFAST(image2);

% 分配方向
[keypoints1, orientations1] = assignOrientation(image1, keypoints1, patch_size);
[keypoints2, orientations2] = assignOrientation(image2, keypoints2, patch_size);

% 计算描述子
descriptors1 = computeDescriptors(image1, keypoints1, patch_size, num_bits);
descriptors2 = computeDescriptors(image2, keypoints2, patch_size, num_bits);

% 进行匹配
matches = matchDescriptors(descriptors1, descriptors2, threshold);
end

注意：这只是一个简单的实现，还有很多可以进行优化的地方。例如，可以使用图像金字塔来实现尺度不变性，可以使用Harris角点检测算法来替代FAST算法等等。