问向量化4个嵌套的for循环

EN

function RunningFunction
    %% generate regions

    visual_field_width = 10;
    center_r = 2;
    center_c = 3;

    deg_per_pixel = 17;
    N_theta = 2;
    N_e = 5;
    e0_in_deg = 35;
    e_max = 17;

    [raw_r, raw_c, nt, ne] = ndgrid(1:visual_field_width, 1:visual_field_width, 1:N_theta, 1:N_e);
    ind = sub2ind(size(raw_r), raw_r, raw_c, nt, ne);

    r = raw_r - center_r;
    c = raw_c - center_c;

    % convert (r,c) to polar: (eccentricity, angle)
    e = sqrt(r.^2+c.^2)*deg_per_pixel;
    a = mod(atan2(r,c),2*pi);

    regions(ind) = ...
        FS_1(nt-1,a,N_theta) .* ...
        FS_2(ne-1,e,N_e,e0_in_deg, e_max);

    regions = reshape(regions, size(raw_r));

end

function h = FS_1(n,theta,N,t)
    if nargin==2
        N = 9;
        t=1/2;
    elseif nargin==3
        t=1/2;
    end

    w = (2*pi)./N;

    theta = theta + w/4;

    theta(n==0 & theta>(3/2)*pi) = theta(n==0 & theta>(3/2)*pi) - 2*pi;

    h = FS_f((theta - (w*n + 0.5*w*(1-t)))/w);
end

function g = FS_2(n,e,N,e0, e_max)

if nargin==2
    N = 10;
    e0 = .5;
elseif nargin==3
    e0 = .5;
end


w = (log(e_max) - log(e0))/N;
g = FS_f((log(e)-log(e0)-w*(n+1))/w);
end

function f = FS_f(x, t)
    if nargin<2
        t = 0.5;
    end

    f = zeros(size(x));

    % case 1
    idx = x>-(1+t)/2 & x<=(t-1)/2;
    f(idx) = (cos(0.5*pi*((x(idx)-(t-1)/2)/t))).^2;

    % case 2
    idx = x>(t-1)/2 & x<=(1-t)/2;
    f(idx) = 1;

    % case 3
    idx = x>(1-t)/2 & x<=(1+t)/2;
    f(idx) = -(cos(0.5*pi*((x(idx)-(1+t)/2)/t))).^2+1;
end