Шаг 1: Визуализация основы ортоизображение в координатах карты

Базовый слой ортоизображения структурирован в плитки 4000 на 4 000 пикселей с каждым пикселем, покрывающим точно один квадратный метр в координатах карты. Каждая плитка сохранена как изображение TIFF с файлом привязки. Воздушная фотография West Concord перекрывает две плитки в базовом слое ортоизображения. (Для удобства этот пример фактически работает с двумя подплитками 2000 на 2000, извлеченными из больших оригиналов 4000 на 4000. Они имеют одинаковый размер пикселя, но покрывают только интересующую область.)

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

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

Ознакомьтесь с файлами привязки для двух плиток

currentFormat = get(0,'format');
format short g
R1 = worldfileread('concord_ortho_w.tfw','planar',size(baseImage1))
R2 = worldfileread('concord_ortho_e.tfw','planar',size(baseImage2))

R1 = 

  MapCellsReference with properties:

            XWorldLimits: [207000 209000]
            YWorldLimits: [911000 913000]
              RasterSize: [2000 2000]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 1
      CellExtentInWorldY: 1
    RasterExtentInWorldX: 2000
    RasterExtentInWorldY: 2000
        XIntrinsicLimits: [0.5 2000.5]
        YIntrinsicLimits: [0.5 2000.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'
            ProjectedCRS: []



R2 = 

  MapCellsReference with properties:

            XWorldLimits: [209000 211000]
            YWorldLimits: [911000 913000]
              RasterSize: [2000 2000]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 1
      CellExtentInWorldY: 1
    RasterExtentInWorldX: 2000
    RasterExtentInWorldY: 2000
        XIntrinsicLimits: [0.5 2000.5]
        YIntrinsicLimits: [0.5 2000.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'
            ProjectedCRS: []

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

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, fitgeotrans, и imwarp, и метод projective2d/transformPointsForward, из Image Processing Toolbox вместе с объектами растровой привязки карты из Mapping Toolbox. Вместе они позволяют регистрирование на основе пар точек управления, которые связывают аэрофотоснимок с базовым слоем ортоизображения.

Прочитайте фото с воздуха.

inputImage = imread('concord_aerial_sw.jpg');
figure
imshow(inputImage)
title('Unregistered Aerial Photograph')

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

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

Понижайте значения изображений в 2 раза, затем выбирайте два отдельных набора пар контрольных точек: один для точек в воздушном изображении, которые появляются в первой плитке, и другой для точек, которые появляются во второй плитке. Сохраните пары управляющих точек в базовом рабочем пространстве как структуры управляющих точек с именами cpstruct1 и cpstruct2.

n = 2; % downsample by a factor n
load mapexreg.mat % load some points that were already picked

При необходимости отредактируйте или добавьте к предварительно выбранным точкам с помощью 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);

Этот инструмент помогает вам выбрать пары соответствующих контрольных точек. Контрольные точки являются ориентирами, которые можно найти на обоих изображениях, таких как перекресток дороги или естественная функция. Для каждой плитки ортофото уже выбрано три пары контрольных точек. Если вы хотите продолжить работу с этими точками, перейдите к шагу 3: Вывод и примените геометрическое преобразование. Если вы хотите добавить несколько дополнительных пар точек, сделайте это, определив ориентиры и щелкнув по изображениям. Сохраните управляющие точки, выбрав меню Файл (File), затем опцию Сохранить точки в рабочей области (Save Points to Workspace). Сохраните точки, перезаписав переменные cpstruct1 и cpstruct2.

Шаг 3: Вывод и применение геометрического преобразования

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

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

Учитывайте понижающую дискретизацию по фактору n.

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

Преобразуйте координаты базового изображения в координаты карты (State Plane).

[x1, y1] = intrinsicToWorld(R1, base1(:,1), base1(:,2));
[x2, y2] = intrinsicToWorld(R2, base2(:,1), base2(:,2));

Объедините два набора управляющих точек и выведите проективное преобразование. (Проективное преобразование должно быть разумным выбором, поскольку воздушное изображение с кадровой камеры, и местность в этой области довольно пология.)

input = [input1; input2];
spatial = [x1 y1; x2 y2];

tform = fitgeotrans(input,spatial,'projective')

tform = 

  projective2d with properties:

                 T: [3x3 double]
    Dimensionality: 2

Вычислите и постройте график ( 2D) невязок.

residuals = transformPointsForward(tform, input) - spatial;
figure
plot(residuals(:,1),residuals(:,2),'+')
title('Control Point Residuals');
xlabel('Easting offset (meters)');
ylabel('Northing offset (meters)');
xlim([-4 4])
ylim([-4 4])
axis equal

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

mInput = size(inputImage,1);
nInput = size(inputImage,2);

inputCorners = 0.5 ...
    + [0        0;
       0        mInput;
       nInput   mInput;
       nInput   0;
       0        0];

outputCornersSpatial = transformPointsForward(tform, inputCorners);

outputCornersX = outputCornersSpatial(:,1);
outputCornersY = outputCornersSpatial(:,2);

figure(ax1.Parent)
line(outputCornersX,outputCornersY,'Color','w')

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

pixelSize = [hypot( ...
    outputCornersX(2) - outputCornersX(1), ...
    outputCornersY(2) - outputCornersY(1)) / mInput, ...
 hypot( ...
    outputCornersX(4) - outputCornersX(5), ...
    outputCornersY(4) - outputCornersY(5)) / nInput]

pixelSize =

      0.90963      0.89054

Переменные pixelSize задает начальную точку, из которой можно выбрать ширину и высоту для пикселей в нашем георганизованном выходном изображении. В принципе мы могли выбрать любой размер для наших выходных пикселей. Однако, если мы сделаем их слишком маленькими, мы теряем память и время расчета, заканчивая большим (много строк и столбцов) размытым изображением. Если мы сделаем их слишком большими, мы рискуем сглаживать изображение, а также ненужно отбрасывать информацию в оригинальное изображение. В целом мы также хотим, чтобы наши пиксели были квадратными. Разумная эвристика состоит в том, чтобы выбрать размер пикселя, который немного больше max(pixelSize) и представляет собой «разумное» число (т.е. 0,95 или 1,0, а не 0,9096). Здесь мы выбрали 1, что означает, что каждый пиксель на нашем георегистрированном изображении будет покрывать один квадратный метр на земле. «Приятно» иметь георгизированные изображения, которые выравниваются по целочисленным координатам карты для отображения и анализа.

outputPixelSize = 1;

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

xWorldLimits = outputPixelSize ...
    * [floor(min(outputCornersX) / outputPixelSize), ...
        ceil(max(outputCornersX) / outputPixelSize)]

yWorldLimits = outputPixelSize ...
    * [floor(min(outputCornersY) / outputPixelSize), ...
        ceil(max(outputCornersY) / outputPixelSize)]

xWorldLimits =

      208154      209693


yWorldLimits =

      911435      912583

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

line(xWorldLimits([1 1 2 2 1]),yWorldLimits([2 1 1 2 2]),'Color','red')

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

mOutput = diff(yWorldLimits) / outputPixelSize
nOutput = diff(xWorldLimits) / outputPixelSize

mOutput =

        1148


nOutput =

        1539

Создайте объект привязки Image Processing Toolbox для зарегистрированного изображения.

R = imref2d([mOutput nOutput],xWorldLimits,yWorldLimits)

R = 

  imref2d with properties:

           XWorldLimits: [208154 209693]
           YWorldLimits: [911435 912583]
              ImageSize: [1148 1539]
    PixelExtentInWorldX: 1
    PixelExtentInWorldY: 1
    ImageExtentInWorldX: 1539
    ImageExtentInWorldY: 1148
       XIntrinsicLimits: [0.5 1539.5]
       YIntrinsicLimits: [0.5 1148.5]

Создайте объект растровой привязки, который является аналогом Mapping Toolbox с объектом привязки на Image Processing Toolbox.

Rmap = maprasterref('RasterSize',R.ImageSize, ...
    'XWorldLimits',R.XWorldLimits,'YWorldLimits',R.YWorldLimits, ...
    'ColumnsStartFrom','north')

Rmap = 

  MapCellsReference with properties:

            XWorldLimits: [208154 209693]
            YWorldLimits: [911435 912583]
              RasterSize: [1148 1539]
    RasterInterpretation: 'cells'
        ColumnsStartFrom: 'north'
           RowsStartFrom: 'west'
      CellExtentInWorldX: 1
      CellExtentInWorldY: 1
    RasterExtentInWorldX: 1539
    RasterExtentInWorldY: 1148
        XIntrinsicLimits: [0.5 1539.5]
        YIntrinsicLimits: [0.5 1148.5]
      TransformationType: 'rectilinear'
    CoordinateSystemType: 'planar'
            ProjectedCRS: []

Применить геометрическое преобразование можно используя imwarp. Разверните его выход так, чтобы столбцы шли с севера на юг.

registered = flipud(imwarp(inputImage, tform, 'OutputView', R));
figure
imshow(registered)

format(currentFormat)

Вы можете записать зарегистрированное изображение как TIFF с файлом привязки, тем самым привязывая его к нашим координатам карты.

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

Шаг 4: отображение зарегистрированного изображения в координатах карты

Отображение зарегистрированного изображения на тех же осях (координаты карты), что и на базовых плитках ортоизображения. Зарегистрированное изображение не полностью заполняет свой ограничивающий прямоугольник, поэтому включает ядро заполненных треугольников. Создайте альфа-маску данных, чтобы эти области заливки стали прозрачными.

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

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

mInput = mapshow(ax2,registered,Rmap);
mInput.AlphaData = alphaData;

Вы можете оценить регистрацию, посмотрев на функции, которые охватывают как зарегистрированное изображение, так и ортофото изображения.

Шаг 5: Наложение векторного дорожного слоя (из Shapefile)

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

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

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

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

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

Заметьте, что roads Линия довольно хорошо с дорогами на изображениях. Двумя очевидными линейными функциями на изображениях являются не дороги, а железные дороги. Линейной функцией, которая имеет тенденцию примерно к востоку-западу и пересекает обе плитки основания, является линия пригородной железной дороги Фитчбурга транспортного агентства Массачусетского залива. Линейной функцией, которая имеет тенденцию примерно к северо-западу-юго-востоку, является бывшая вторичная линия Фрамингем-Лоуэлл.

Кредиты

concord_orthow_w.tif, concord_ortho_e.tif, concord_roads.shp:

  Office of Geographic and Environmental Information (MassGIS),
  Commonwealth of Massachusetts  Executive Office of Environmental Affairs
  http://www.state.ma.us/mgis

  For more information, run:

  >> type concord_ortho.txt
  >> type concord_roads.txt

concord_aerial_sw.jpg

  Visible color aerial photograph courtesy of mPower3/Emerge.

  For more information, run:

  >> type concord_aerial_sw.txt