algorithm - 如何生成具有函数模式的矩阵？

Question

背景信息(可选阅读):

我正在模拟声波在边界内的反射。空间中点的介质条件使用矩阵设置。假设空间的维度是 N 乘 N 网格，我关心两种声速，c0 和 c1.

现在我正在使用如下代码生成障碍模式

medium.sound_speed = c0*ones(N,N);   % set the speed of sound to be c0 everywhere
medium.sound_speed(:, N/2:N) = c1;   % set the right half of the grid to a different speed
medium.sound_speed(50:70, 50:70) = c1; % set a box to have a different speed

或者

% set all speeds to c0 except set the diagonal to c1
medium.sound_speed = c0*ones(N,N)-(c0*eye(N,N))+c1*eye(N,N); 

但是，我无法生成具有不同曲率的更复杂的边界。

问题

我想以编程方式创建具有反射(reflect)函数模式的矩阵。例如，我想输入 f(x)=2 并为此创建一个看起来像这样的矩阵，假设 N=6。

[ 0 0 0 0 0 0
  0 0 0 0 0 0
  0 0 0 0 0 0 
  1 1 1 1 1 1 
  0 0 0 0 0 0 
  0 0 0 0 0 0 ]

或者f(x)=0.5*x+1

[ 0 0 0 0 0 0
  0 0 0 0 0 0
  0 0 0 0 1 1
  0 0 1 1 0 0 
  1 1 0 0 0 0 
  0 0 0 0 0 0]

我还可以生成类似 f(x)=1/x 的曲线图案，这似乎需要某种形式的 Midpoint circle algorithm , 用于绘制带像素的曲率。

[ 1 0 0 0 0 0
  1 0 0 0 0 0
  0 1 0 0 0 0
  0 0 1 1 0 0
  0 0 0 0 1 1  
  0 0 0 0 0 0 ]

实际上，N 至少为 128，因此手动创建具有一定复杂度的形状矩阵是不切实际的，我认为这是一个有趣的问题。

有没有人知道一些方法来做到这一点，或者对替代方法的建议？

提前谢谢你。

编辑:我修改了this实现 Bresenham 算法以提供具有给定起点和终点的所需直线的矩阵。

function M=bresenham_line(point)

if (abs(point(4)-point(2)) > abs(point(3)-point(1)))       % If the line is steep                                
    x0 = point(2);y0 = point(1); x1 = point(4);y1=point(3);% then it would be converted to 
    token =1;                                              % non steep by changing coordinate
else
    x0 = point(1);y0 = point(2); x1 = point(3);y1=point(4);
    token = 0; 
end
if(x0 >x1)
    temp1 = x0; x0 = x1; x1 = temp1;
    temp2 = y0; y0 = y1; y1 = temp2;
end
dx = abs(x1 - x0) ;                              % Distance to travel in x-direction
dy = abs(y1 - y0);                               % Distance to travel in y-direction
sx = sign(x1 - x0);                              % sx indicates direction of travel in X-dir
sy = sign(y1 - y0);                              % Ensures positive slope line

x = x0; y = y0;                                  % Initialization of line
param = 2*dy - dx ;                              % Initialization of error parameter
for i = 0:dx-1                                   % FOR loop to travel along X
    x_coord(i+1) = x;                            % Saving in matrix form for plot
    y_coord(i+1) = y;

    param = param + 2*dy;                        % parameter value is modified
    if (param >0)                                % if parameter value is exceeded
        y = y +1*sy;                             % then y coordinate is increased
        param = param - 2*(dx );                 % and parameter value is decreased

    end
    x = x + 1*sx;                                % X-coordinate is increased for next point
end

M = zeros(size(x_coord,2), size(y_coord,2));

for i=1:1:size(x_coord,2)
    x = x_coord(i);
    y = y_coord(i);

    M(x,y) = 1;
end

M

这样实现:

c1 = 0;
M = bresenham_line([1 1 Nx/2+1 Ny+1]);
medium.sound_speed = c0*ones(Nx,Ny) - (c0*M) + c1*M; 

在曲线函数形状方面还没有进展。

Answer 1

获得一些类似结果的方法:

f = @(x)0.5*x; %create the function (x should be written even if the function doesn't depend on x: @(x) 0*x + 2)
N = 6; %choose the size of the atrix
M = zeros(N,N); %create an empty matrix
x = (1:N); 
y = round(f(x-1)); %discretization
x(y>N-1|y<0) = [];
y(y>N-1|y<0) = [];
M(sub2ind(size(M),y+1,x)) = 1;
M = flipud(M)

因此您可以选择您的函数，然后矩阵中的结果将看起来像正常 plot 的离散化。