Расширение изображения с преобразованием наклонных дальностей в земные координаты

clear, clc
 
A = imread('2013.08.01_13.54_.bmp'); % считываем фото
[h w] = size(A); % размеры изображения
Z = double(A); % переводим в числовой формат
Z = medfilt2( Z/max(Z(:)), [10 10]); % нормируем и сглаживаем шумы
% Z(Z<0.05) = 0; %обнуляем все что меньше 0.05
H = 600; % высота
Rmin = 933; % мин. растояние
Rmax = 2366; % макс. расстояние
Rx = linspace(Rmin, Rmax, w); % масив точек по ширине
dr = Rx(2)-Rx(1); % шаг длинны
Ry = 933 + (0:h-1)*dr; % масив точек по высоте
 
% расчет дальностей по поверхности
Dx = sqrt(Rx.^2 - H^2);
Dy = sqrt(Ry.^2 - H^2);
 
subplot(211)
surf(Rx,Ry,-Z)
shading flat
colormap gray
view(0,90)
xlim([500 2500])
ylim([600 1800])
axis equal
title('Z(R)')
xlabel('Rx')
ylabel('Ry')
 
subplot(212)
surf(Dx,Dy,-Z)
shading flat
view(0,90)
xlim([500 2500])
ylim([600 1800])
axis equal
title('Z(D)')
xlabel('Dx')
ylabel('Dy')

clc; clear all;
format short g
fname='2013.08.01_13.54_.bmp'; 
I = imread(fname);
[h w] = size(I); % размеры изображения
%Z = double(I); % переводим в числовой формат
%Z = medfilt2( Z/max(Z(:)), [10 10]); % нормируем и сглаживаем шумы
R = (1:1024).*1.4 + 715;
D = sqrt(R.^2 - 600^2);
X = linspace(min(D), max(D),10*size(I,2)); % новый расширенный массив 
d = 1;
A = zeros(size(I,1), size(I,2));
for i = 2:size(I,2)
    [u d(i)] = min( abs(X-D(i)) ); 
    % d - номер столбца с единицами, номер элемента в массиве m
    j = d(i-1):d(i); 
    A(:,j) = I(:,(i-1)*ones(size(j))); % записываем столбец
end
%imagesc(I)
%colormap(gray)
%figure;
imagesc( A(:,1:10:end));
colormap(gray)
axis off
%imwrite( uint8(I) ,'D:\Шмель\15654.png')