如何用函数模式生成矩阵？

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); 

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

问题

我想通过程序生成带有图案的反射功能矩阵。例如，我想输入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实现布氏算法，以提供给定的起点和终点所需的行的矩阵。

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离散化。

Answer 2

这是一个稍微“肮脏”的方式来得到这样的东西，但我最好打赌可能Bresenham的算法。

N = 128; 
[X,Y] = meshgrid(1:N,1:N); 
bound1 = Y<2*X; 
bound2 = Y<2*X+1; 
M = xor(bound1,bound2); 

bound1你可以定义任何功能y=f(x)，并根据它标出的区域。与bound2你选择和稍高一点的区域（向上移动）。一旦你拿到和xor这两个区域，你只需要标记y=f(x)。我认为，为了获得合理的结果，这种转变可能会因更复杂的功能而有所不同。

为了说明我用imagesc（该flipud仅仅是让(0,0)在左下角，而不是左上角）：

imagesc(flipud(M)); 

---

编辑

事实上，对于一些功能这可能不是最好的。例如对于y=x^2，您必须增加班次，并且看起来不太好。

bound1 = Y<X.^2; 
bound2 = Y<X.^2+15; 
M = xor(bound1,bound2);