Шаг 1: представьте мозаики ортобазового адреса образа в координатах карты

Слой ортобазового адреса образа структурирован в 4000 4000 пиксельные мозаики с каждым пикселем, покрывающим точно один квадратный метр в координатах карты. Каждая мозаика хранится как изображение TIFF с файлом привязки. Воздушная фотография Западного Согласия перекрывает две мозаики в базовом слое ортоизображений. (Для удобства этот пример на самом деле работает с два 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));

Downsample изображения на коэффициент 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, затем Точки сохранения к опции Рабочей области. Сохраните точки, перезаписав переменные 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;

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

[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) и "разумный" номер (i.e., 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: слой Оверлей Вектор-Роуд (из файла форм)

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

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

Создайте изображение условными знаками на основе Атрибута класса (тип дороги). Обратите внимание на то, что очень незначительные дороги имеют значения КЛАССА, равные 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