Matlab - 创建两个考虑两个对象所有可能值的数组
Matlab - Creating two arrays that consider all possible values for two objects
我有两个物理"objects."我用两个不同的数组表示它们的位置。
• 对象 1 仅在 xy 平面内移动
• 对象 2 在所有三个物理维度上移动
Objective: 向量化四个 for
循环而不扭曲数据。此外,目的是对要与对象 2 进行比较的对象 1 的所有可能值执行此操作。
这是代码片段:
Npos = 21;
Nsam = 200;
% dummy initialisation
AX = rand(1, Npos);
AY = zeros(1, Npos);
AZ = rand(1, Npos);
Bx = rand(Nsam);
By = rand(Nsam);
Bz = rand(Nsam);
for qx = 1 : Npos
for yx = 1 : Npos
for zx = 1 : Nsam
for cx = 1 : Nsam
Tx2Array( qx, yx, zx, cx ) = sqrt( ( AX( qx ) - Bx( zx, cx ) ).^2 + ( AY( yx ) - By( zx, cx ) ).^2 + ( AZ( yx ) - Bz( zx, cx ) ).^2 );
end
end
end
end
% Result is a 21 x 21 x 200 x 200 matrix filled with all real numbers
图例
AX、AY、AZ 是 1 x 21 数组,代表对象 1 的 (x,y=0,z)
AY 全部为零,但仍包括可读性(因此没有第五个循环!)
Bx、By、Bz 都是 200 x 200 数组,表示对象 2 的 (x,y,z)
Npos = 21; Nsam = 200;
上面使用的公式是:
sqrt( (a1-b1)^2 + (a2-b2)^2 + (a3-b3)^2 )
您可以通过将 zx
和 cx
替换为 :
来避免内部循环,如下所示:
Tx2Array = zeros(Npos, Npos, Nsam, Nsam); % preallocate memory
for qx = 1 : Npos
for yx = 1 : Npos
Tx2Array( qx, yx, :, : ) = sqrt( ( AX( qx ) - Bx( :, : ) ).^2 + ( AY( yx ) - By( :, : ) ).^2 + ( AZ( yx ) - Bz( :, : ) ).^2 );
end
end
这样,最大的维度就被矢量化了。所以,最大的改进已经完成了。
通过将您的 B*
转换为 4D 并为您的 A*
矩阵生成网格,您甚至可以删除所有 for 循环,如下所示:
[AX_, AZ_] = meshgrid(AX, AZ);
AX_ = AX_';
AZ_ = AZ_';
AY_ = zeros(Npos);
Bx_(1, 1, :, :) = Bx;
By_(1, 1, :, :) = By;
Bz_(1, 1, :, :) = Bz;
Tx2Array2 = sqrt( ( AX_ - Bx_ ).^2 + ( AY_ - By_ ).^2 + ( AZ_ - Bz_ ).^2 );
您可以使用以下检查来检查结果是否相同:
max(max(max(max(abs(Tx2Array - Tx2Array2))))) < eps
如果您有可用的统计工具箱,您可以使用pdist2计算对象 1 的每个坐标与对象 2 的每个坐标之间的距离:
[X1, Z1] = ndgrid(AX(:), AZ(:)); % X1 and Z1 will be 21x21
D = pdist2([X1(:), zeros(size(X1(:))), Z1(:)], [Bx(:), By(:), Bz(:)]);
这种情况下的输出将是一个 441 x 40,000
数组,其中 D(i, j)
给出对象 1 的点 i
和对象 2 的点 j
之间的距离,都使用线性索引。
如果数组正确初始化,您的任务将非常简单:
用正确的维度初始化数组
AX = rand( Npos,1);
AY = zeros(1, Npos);
AZ = rand(1, Npos);
Bx = rand(1,1,Nsam,Nsam);
By = rand(1,1,Nsam,Nsam);
Bz = rand(1,1,Nsam,Nsam);
然后在 MATLAB r2016b / Octave 中你可以简单地写:
Tx2Array = sqrt( ( AX - Bx ).^2 + ( AY - By ).^2 + ( AZ - Bz ).^2 );
在 r2016b 之前你可以使用 bsxfun:
Tx2Array = sqrt(bsxfun(@plus,bsxfun(@plus,bsxfun(@minus,AX , Bx).^2,bsxfun(@minus,AY , By).^2),bsxfun(@minus,AZ , Bz).^2));
我有两个物理"objects."我用两个不同的数组表示它们的位置。
• 对象 1 仅在 xy 平面内移动
• 对象 2 在所有三个物理维度上移动
Objective: 向量化四个 for
循环而不扭曲数据。此外,目的是对要与对象 2 进行比较的对象 1 的所有可能值执行此操作。
这是代码片段:
Npos = 21;
Nsam = 200;
% dummy initialisation
AX = rand(1, Npos);
AY = zeros(1, Npos);
AZ = rand(1, Npos);
Bx = rand(Nsam);
By = rand(Nsam);
Bz = rand(Nsam);
for qx = 1 : Npos
for yx = 1 : Npos
for zx = 1 : Nsam
for cx = 1 : Nsam
Tx2Array( qx, yx, zx, cx ) = sqrt( ( AX( qx ) - Bx( zx, cx ) ).^2 + ( AY( yx ) - By( zx, cx ) ).^2 + ( AZ( yx ) - Bz( zx, cx ) ).^2 );
end
end
end
end
% Result is a 21 x 21 x 200 x 200 matrix filled with all real numbers
图例
AX、AY、AZ 是 1 x 21 数组,代表对象 1 的 (x,y=0,z)
AY 全部为零,但仍包括可读性(因此没有第五个循环!)
Bx、By、Bz 都是 200 x 200 数组,表示对象 2 的 (x,y,z)
Npos = 21; Nsam = 200;
上面使用的公式是:
sqrt( (a1-b1)^2 + (a2-b2)^2 + (a3-b3)^2 )
您可以通过将 zx
和 cx
替换为 :
来避免内部循环,如下所示:
Tx2Array = zeros(Npos, Npos, Nsam, Nsam); % preallocate memory
for qx = 1 : Npos
for yx = 1 : Npos
Tx2Array( qx, yx, :, : ) = sqrt( ( AX( qx ) - Bx( :, : ) ).^2 + ( AY( yx ) - By( :, : ) ).^2 + ( AZ( yx ) - Bz( :, : ) ).^2 );
end
end
这样,最大的维度就被矢量化了。所以,最大的改进已经完成了。
通过将您的 B*
转换为 4D 并为您的 A*
矩阵生成网格,您甚至可以删除所有 for 循环,如下所示:
[AX_, AZ_] = meshgrid(AX, AZ);
AX_ = AX_';
AZ_ = AZ_';
AY_ = zeros(Npos);
Bx_(1, 1, :, :) = Bx;
By_(1, 1, :, :) = By;
Bz_(1, 1, :, :) = Bz;
Tx2Array2 = sqrt( ( AX_ - Bx_ ).^2 + ( AY_ - By_ ).^2 + ( AZ_ - Bz_ ).^2 );
您可以使用以下检查来检查结果是否相同:
max(max(max(max(abs(Tx2Array - Tx2Array2))))) < eps
如果您有可用的统计工具箱,您可以使用pdist2计算对象 1 的每个坐标与对象 2 的每个坐标之间的距离:
[X1, Z1] = ndgrid(AX(:), AZ(:)); % X1 and Z1 will be 21x21
D = pdist2([X1(:), zeros(size(X1(:))), Z1(:)], [Bx(:), By(:), Bz(:)]);
这种情况下的输出将是一个 441 x 40,000
数组,其中 D(i, j)
给出对象 1 的点 i
和对象 2 的点 j
之间的距离,都使用线性索引。
如果数组正确初始化,您的任务将非常简单:
用正确的维度初始化数组
AX = rand( Npos,1);
AY = zeros(1, Npos);
AZ = rand(1, Npos);
Bx = rand(1,1,Nsam,Nsam);
By = rand(1,1,Nsam,Nsam);
Bz = rand(1,1,Nsam,Nsam);
然后在 MATLAB r2016b / Octave 中你可以简单地写:
Tx2Array = sqrt( ( AX - Bx ).^2 + ( AY - By ).^2 + ( AZ - Bz ).^2 );
在 r2016b 之前你可以使用 bsxfun:
Tx2Array = sqrt(bsxfun(@plus,bsxfun(@plus,bsxfun(@minus,AX , Bx).^2,bsxfun(@minus,AY , By).^2),bsxfun(@minus,AZ , Bz).^2));