MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Математика\Mapping Toolbox

В оглавление \ К предыдущему разделу \ К следующему разделу

Географическая привязка на изображениях

В этом примере рассмотрено один из путей географической привязки на изображениях путем регистрации изображений в земной системе координат и создании нового “привязанного” изображения. Для просмотра данного демонстрационного примера необходимо наличие Image Processing Toolbox.

Все географические данные в этом демонстрационном примере представлены в земной системе координат. Данные представляют собой цифровые аэрофотографические изображения части West Concord, (Massachusetts), полученные ранней весной 1997 года.

Содержание:

Шаг 1: Нанесение на основную часть ортоизображения картографических координат

Структурно ортоизображение состоит из 4000х4000 пиксельных частей, каждая из которых покрывает один квадратный метр в картографических координатах. Из таких частей, которые представлены в формате TIFF, состоит карта мира. Аэрофотография West Concord состоит из двух перекрывающихся частей. ( Для удобства в этом демонстрационном примере используются две части изображения размером 2000х2000, которые получены из исходного оригинала с размером 4000х4000 пикселей. Эти части имеют одинаковый пиксельный размер.)

Считаем данные для двух частей

currentFormat = get(0,'format');
format bank
R1 = worldfileread('concord_ortho_w.tfw')
R2 = worldfileread('concord_ortho_e.tfw')
format(currentFormat)

R1 =

          0         -1.00
       1.00             0
  206999.50     913000.50

R2 =

          0         -1.00
       1.00             0
  208999.50     913000.50

Считаем две основные части ортофото, которые необходимы для покрытия аэрофотографического изображения.

[baseImage1,cmap1] = imread('concord_ortho_w.tif');
[baseImage2,cmap2] = imread('concord_ortho_e.tif');

Отобразим изображения в соответствующих им пространственных координатах.

mapshow(baseImage1,cmap1,R1)
ax1 = gca;
mapshow(ax1,baseImage2,cmap2,R2)
title('Map View, Massachusetts State Plane Coordinates');
xlabel('Easting (meters)');
ylabel('Northing (meters)');

Отображение изображения в соответствующих им пространственных координатах

Шаг 2: Регистрация аэрофотоснимков в картографических координатах

На этом этапе используются функции cpselect, cpstruct2pairs, cp2tform, tformfwd и imtransform из приложения Image Processing Toolbox. Они используются в сочетании с функцией pix2map. Вместе они дают возможность реализации географической привязки на основе пар контрольных точек, которые совмещают аэрофотосъемку и ортоизображения.

Считаем аэрофотоизображения

inputImage = imread('concord_aerial_sw.jpg');
figure, imshow(inputImage)
title('Unregistered Aerial Photograph')
Warning: Image is too big to fit on screen; displaying at 67%

Считаем аэрофотоизображения

Оба ортоизображения представляют собой полноцветные изображения, которые далее преобразовываются в полутоновые.

baseGray1 = im2uint8(ind2gray(baseImage1,cmap1));
baseGray2 = im2uint8(ind2gray(baseImage2,cmap2));

Обрабатываются изображения с применением пар контрольных точек.

n = 2; % коэффициент разделения n
load mapexreg.mat % считывание выбранных контрольных точек

Выбранные контрольные точки можно редактировать с помощью функции cpselect. Отметим, что для работы необходимо использовать отдельные контрольные точки для каждой части изображения.

cpselect(inputImage(1:n:end,1:n:end,1),...
     baseGray1(1:n:end,1:n:end),cpstruct1);
cpselect(inputImage(1:n:end,1:n:end,1),...
     baseGray2(1:n:end,1:n:end),cpstruct2);

Legend

Control point selection tool

Control point selection tool 2

Это приложение помогает работать с контрольными точками для совмещения и привязки изображений. Для запоминания контрольных точек необходимо выбрать меню File, далее опцию Save Points to Workspace.

Шаг 3: Определение и применение геометрических преобразований

Извлечем пары контрольных точек из структуры контрольных точек.

[input1,base1] = cpstruct2pairs(cpstruct1);
[input2,base2] = cpstruct2pairs(cpstruct2);

Определим субдискретизацию для коэффициента n.

input1 = n*input1 - 1;
base1  = n*base1  - 1;
input2 = n*input2 - 1;
base2  = n*base2  - 1;

Преобразуем координаты изображения в координаты карты.

spatial1 = pix2map(R1,fliplr(base1));
spatial2 = pix2map(R2,fliplr(base2));

Объединим два набора контрольных точек и проанализируем проективные преобразования.

input   = [input1;     input2];
spatial = [spatial1; spatial2];

tform = cp2tform(input,spatial,'projective');

Вычислим и отобразим результат.

residuals = tformfwd(input,tform) - spatial;
figure
plot(residuals(:,1),residuals(:,2),'.')
title('Control Point Residuals');
xlabel('Easting offset (meters)');
ylabel('Northing offset (meters)');

Отображение результата

Спрогнозируем размещение углов регистрируемого изображения в координатах карты и свяжем с контурами изображения.

w = size(inputImage,2);
h = size(inputImage,1);
inputCorners = [0  0;
          w  0;
          w  h;
          0  h;
          0  0] + .5;
outputCornersSpatial = tformfwd(inputCorners,tform);
figure(get(ax1,'Parent'))
line(outputCornersSpatial(:,1),outputCornersSpatial(:,2),'Color','w')

Размещение углов регистрируемого изображения в координатах карты

Вычислим средний размер пикселей для исходного изображения (в единицах карты).

d = [outputCornersSpatial(2,:) - outputCornersSpatial(1,:);...
     outputCornersSpatial(5,:) - outputCornersSpatial(4,:)];
[theta,lengths] = cart2pol(d(:,1),d(:,2));
pixelSize = lengths./[w; h]
pixelSize =

    0.8906
    0.9096

Переменная pixelSize представляет исходные точки исходя из выбора высоты и ширины пикселей результирующего изображения. Сначала мы выбираем произвольный размер для всех пикселей. Однако, если поставить исходный размер очень маленьким, то это приведет к нерациональному использованию памяти и время вычислений будет очень большим. Если же исходный размер пикселей сделать слишком большим, то существует опасность ступенчатого представления изображения, которое неизбежно приводит к потере информации. В общем, мы хотим, чтобы наши пиксели были квадратными.

outputPixelSize = 1;

Выбор ограничивающего прямоугольника.

outputBoundingBox = [floor(min(outputCornersSpatial)); ...
                     ceil(max(outputCornersSpatial))];
left   = outputBoundingBox(1);
right  = outputBoundingBox(2);
bottom = outputBoundingBox(3);
top    = outputBoundingBox(4);
outputBoundingBoxClose = [left  top;
                left  bottom;
                right bottom;
                right top;
                left  top];

Отобразим ограничивающий прямоугольник на изображении.

line(outputBoundingBoxClose(:,1),outputBoundingBoxClose(:,2),'Color','r')

Отображение ограничивающего прямоугольника на изображении

Вычислим координаты верхнего левого и нижнего правого углов.

outputCornerCenters = ...
    [left top; right bottom]...
    + [0.5 -0.5; -0.5 0.5] * outputPixelSize
outputCornerCenters =

  1.0e+005 *

    2.0815    9.1258
    2.0969    9.1144

Преобразуем аэрофотоизображение в изображение, которое совпадает с ограничивающим прямоугольником.

registered = imtransform(inputImage,tform,...
           'XData',outputCornerCenters(:,1)',...
           'Ydata',outputCornerCenters(:,2)',...
           'XYScale',outputPixelSize);
figure, imshow(registered)
Warning: Image is too big to fit on screen; displaying at 67%

Преобразование аэрофотоизображение в изображение

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

alphaData = registered(:,:,1);
alphaData(alphaData~=0) = 255;

Создадим маску для регистрируемого изображения.

Rregistered = makerefmat(outputCornerCenters(1,1),...
             outputCornerCenters(1,2),...
             outputPixelSize,...
            -outputPixelSize);

Проверим новую пространственную структуру для сравнения с выделенным прямоугольником.

sizeRegistered = size(registered);
boundingBoxCheck = mapbbox(Rregistered,sizeRegistered(1:2));
isequal(boundingBoxCheck,outputBoundingBox)
ans =

     1

Регистрируемое изображение можно сохранить в формате TIFF.

imwrite(registered,'concord_aerial_sw_reg.tif');
worldfilewrite(Rregistered,getworldfilename('concord_aerial_sw_reg.tif'));

Шаг 4: Визуализация сформированного изображения в картографических координатах

Визуализируем регистрируемое изображение в тех же координатах карты, что и основная часть ортоизображения.

mapshow(baseImage1,cmap1,R1)
ax2 = gca;
mapshow(ax2,baseImage2,cmap2,R2)
title('Map View, Massachusetts State Plane Coordinates');
xlabel('Easting (meters)');
ylabel('Northing (meters)');

h = mapshow(ax2,registered,Rregistered);
set(h,'AlphaData',alphaData)

Визуализация регистрируемого изображения

Шаг 5: Нанесение данных вектора дорог (из файла форм)

Используя функции shapeinfo и shaperead мы можем получить данные из вектора дорог.

roadsfile  = 'concord_roads.shp';
roadinfo   = shapeinfo(roadsfile);
roads = shaperead(roadsfile);

Создадим символьное изображение дорог.

roadspec = makesymbolspec('Line',...
               {'CLASS',6,'Visible','off'});

mapshow(ax2,roads,'SymbolSpec',roadspec,'Color','cyan')

Создание символьного изображения дорог

Отметим, что сформированное символьное обозначение дорог совпадает с их реальным размещением на изображении.


В оглавление \ К предыдущему разделу \ К следующему разделу

 


Поиск по сайту:

Система Orphus

Яндекс.Метрика