用2023a版matlab写一段将四副256*256的灰度图像I1、I2、I3、I4进行图像归一化处理，图像归一化的原理是：首先计算并确定原始图像对各种仿射变换具有较强不变性的图像矩，然后根据得到的图像矩确定变换函数的具体参数，最后利用变换函数将原始图像转换为具有标准形式的图像，其归一化过程主要包括4个步骤：平移归一化、X-shearing归一化、Y-shearing归一化和缩放归一化，原始图像在经过以上四个步骤后，便能够转化为标准图像，从而能够降低仿射变换对图像的影响

以下是用 Matlab 编写的图像归一化处理代码：

% 读取四幅图像
I1 = imread('image1.jpg');
I2 = imread('image2.jpg');
I3 = imread('image3.jpg');
I4 = imread('image4.jpg');

% 计算图像矩
m10_1 = moment(I1, 1, 0);
m01_1 = moment(I1, 0, 1);
m20_1 = moment(I1, 2, 0);
m02_1 = moment(I1, 0, 2);
m11_1 = moment(I1, 1, 1);

m10_2 = moment(I2, 1, 0);
m01_2 = moment(I2, 0, 1);
m20_2 = moment(I2, 2, 0);
m02_2 = moment(I2, 0, 2);
m11_2 = moment(I2, 1, 1);

m10_3 = moment(I3, 1, 0);
m01_3 = moment(I3, 0, 1);
m20_3 = moment(I3, 2, 0);
m02_3 = moment(I3, 0, 2);
m11_3 = moment(I3, 1, 1);

m10_4 = moment(I4, 1, 0);
m01_4 = moment(I4, 0, 1);
m20_4 = moment(I4, 2, 0);
m02_4 = moment(I4, 0, 2);
m11_4 = moment(I4, 1, 1);

% 平移归一化
tx_1 = 0.5 - m10_1 / m00_1;
ty_1 = 0.5 - m01_1 / m00_1;
I1 = imtranslate(I1, [tx_1, ty_1]);

tx_2 = 0.5 - m10_2 / m00_2;
ty_2 = 0.5 - m01_2 / m00_2;
I2 = imtranslate(I2, [tx_2, ty_2]);

tx_3 = 0.5 - m10_3 / m00_3;
ty_3 = 0.5 - m01_3 / m00_3;
I3 = imtranslate(I3, [tx_3, ty_3]);

tx_4 = 0.5 - m10_4 / m00_4;
ty_4 = 0.5 - m01_4 / m00_4;
I4 = imtranslate(I4, [tx_4, ty_4]);

% X-shearing归一化
theta_x_1 = atan2(2 * (m11_1 / m00_1 - ty_1 * tx_1), m20_1 / m00_1 - tx_1^2);
I1 = imrotate(I1, -theta_x_1 * 180 / pi, 'crop');

theta_x_2 = atan2(2 * (m11_2 / m00_2 - ty_2 * tx_2), m20_2 / m00_2 - tx_2^2);
I2 = imrotate(I2, -theta_x_2 * 180 / pi, 'crop');

theta_x_3 = atan2(2 * (m11_3 / m00_3 - ty_3 * tx_3), m20_3 / m00_3 - tx_3^2);
I3 = imrotate(I3, -theta_x_3 * 180 / pi, 'crop');

theta_x_4 = atan2(2 * (m11_4 / m00_4 - ty_4 * tx_4), m20_4 / m00_4 - tx_4^2);
I4 = imrotate(I4, -theta_x_4 * 180 / pi, 'crop');

% Y-shearing归一化
theta_y_1 = atan2(2 * (m11_1 / m00_1 - ty_1 * tx_1), m02_1 / m00_1 - ty_1^2);
I1 = imrotate(I1, -theta_y_1 * 180 / pi, 'crop');

theta_y_2 = atan2(2 * (m11_2 / m00_2 - ty_2 * tx_2), m02_2 / m00_2 - ty_2^2);
I2 = imrotate(I2, -theta_y_2 * 180 / pi, 'crop');

theta_y_3 = atan2(2 * (m11_3 / m00_3 - ty_3 * tx_3), m02_3 / m00_3 - ty_3^2);
I3 = imrotate(I3, -theta_y_3 * 180 / pi, 'crop');

theta_y_4 = atan2(2 * (m11_4 / m00_4 - ty_4 * tx_4), m02_4 / m00_4 - ty_4^2);
I4 = imrotate(I4, -theta_y_4 * 180 / pi, 'crop');

% 缩放归一化
s_1 = sqrt(m00_1);
I1 = imresize(I1, [256, 256] * s_1);

s_2 = sqrt(m00_2);
I2 = imresize(I2, [256, 256] * s_2);

s_3 = sqrt(m00_3);
I3 = imresize(I3, [256, 256] * s_3);

s_4 = sqrt(m00_4);
I4 = imresize(I4, [256, 256] * s_4);

% 显示归一化后的图像
imshow(I1);
imshow(I2);
imshow(I3);
imshow(I4);

需要注意的是，上述代码中的 `moment` 函数需要自己实现，可以参考下面的代码：

function m = moment(I, p, q)
    [h, w] = size(I);
    [X, Y] = meshgrid(1:w, 1:h);
    m = sum(sum(double(I) .* X.^p .* Y.^q));
end