Пересечение луча и поверхности

Я новичок в Matlab, и, возможно, эта проблема несложная, но у меня есть проблемы. Мой код выглядит так

x=-25:0.1:25; y=-25:0.1:25; [X,Y] = meshgrid(x,y);
R=sqrt(X.^2+Y.^2);
Z=sin(R);

faces=delaunay(X,Y);
alpha(0.3);
trisurf(faces,X,Y,Z,'Edgecolor','none');
axis ([-20 20 -20 20 -1 51]);
axis equal;
xlabel('X');
ylabel('Y');
alpha(1);
hold on
syms t;
x_source=20;  
x_dir=-2;
y_source=0;   
y_dir=3;
z_source=50;  
z_dir=-10;

height_above_plane = @(t) z_source + t * z_dir - interp2(X, Y, Z, ...
    x_source + t*x_dir, y_source + t*y_dir);
t_intercept = fzero(height_above_plane, 0)

x_ray = x_source + t_intercept * x_dir;
y_ray = y_source + t_intercept * y_dir;
z_ray = z_source + t_intercept * z_dir;
plot(x_ray,y_ray,'r.','MarkerSize',10);

Как будто луч начинается в позиции 20;0;50. То, что я пытаюсь сделать, это процедура виртуального лазерного сканирования. fzero дает мне NaN, когда x_dir превышает небольшое число, и я не могу понять, почему... Будем очень признательны за любую помощь!

Плохой цикл после решения проблемы:

t=0;
x_source=20;
y_source=0;
z_source=50;
z_dir=-10;
for alfa=26:0.5:6
 x_dir=tand(alfa)*z_dir;
    for beta=16:0.5:-16
    y_dir=tand(beta)*z_dir;    
  y_dir=3;
    height_above_plane = @(t) z_source + t * z_dir - interp2(X, Y, Z, ...
    x_source + t*x_dir, y_source + t*y_dir,'linear',0);
    t_intercept = fmincon(@(t) abs(height_above_plane(t)),0,[],[],[],[],4.5,6)

    x_ray = x_source + t_intercept * x_dir;
    y_ray = y_source + t_intercept * y_dir;
    z_ray = z_source + t_intercept * z_dir;
    plot(x_ray,y_ray,'k.','MarkerSize', 20);
    end;
end;

person Paulius Razma    schedule 07.12.2015    source источник


Ответы (1)


Полная ошибка

Exiting fzero: aborting search for an interval containing a sign change
    because NaN or Inf function value encountered during search.
(Function value at -2.56 is NaN.)
Check function or try again with a different starting value.

и ваша функция, оцененная как -2,56, действительно возвращает NaN. Это потому, что вы делаете interp2(X,Y,Z,x_source + t*x_dir, y_source + t*y_dir), но x_source-2.56*x_dir=25.12. Поскольку ваши данные относятся только к x между -25 и 25, вам нужно экстраполировать, чтобы получить здесь интерполянт. По умолчанию Matlab возвращает NaN для экстраполяции, но вы можете изменить значение с помощью необязательного пятого аргумента на interp2, EXTRAPVAL (добавьте ,'linear',0 к вызову interp2, проверьте документацию (0 тоже может быть не лучшим выбором)).

Учитывая, что вам нужно решение корня с ограничениями, рассмотрите возможность использования вместо этого fmincon, например так:

t_intercept = fmincon(@(t) abs(height_above_plane(t)),0,[],[],[],[],-2.5,25/3)

который медленнее, но более надежен, если у вас есть способ определить, какими должны быть ограничения для t.

person David    schedule 07.12.2015
comment
Большое спасибо, это сработало! Однако моя попытка поместить это в цикл for не удалась. Однако это не дает мне никаких ошибок, поэтому я снова не знаю, что делать. Я отредактировал исходный пост, чтобы показать цикл. - person Paulius Razma; 08.12.2015
comment
Вам придется изменить 4,5 и 6 для каждой итерации, которую я представляю. - person David; 08.12.2015
comment
О, я думал, что это интервал для моей переменной t. решение находится в диапазоне lb‹= x ‹=ub это то, что говорит справка Matlab - person Paulius Razma; 08.12.2015
comment
Да, это диапазон для t, но диапазон для t устанавливается путем обеспечения -25 < x_source + t*x_dir < 25, а также для y_source, чтобы можно было выполнить интерполяцию. Может быть, попробовать max((x_source-25)/x_dir,(y_source-25)/y_dir) для нижней границы и что-то подобное (но с min) для верхней границы. - person David; 08.12.2015
comment
Поскольку z_dir является константой -10, а минимальное расстояние от источника равно 49,5, мы знаем, что t никогда не будет меньше 4,5, но забавно то, что даже такой цикл, как | for blob = 45:0.5:5 disp(blob) end | у меня в отдельном скрипте сейчас работает на матлабе, наверное пора спать... - person Paulius Razma; 08.12.2015
comment
Спасибо за помощь! - person Paulius Razma; 08.12.2015