Шаг 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, из панели инструментов обработки изображений вместе с растровыми опорными объектами карты из панели инструментов сопоставления. Вместе они обеспечивают георегистрацию на основе пар контрольных точек, которые связывают аэрофотосъемку с базовым слоем ортоизображений.

Читайте аэрофотоснимок.

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: Создать (Infer) и примените геометрическое преобразование. Чтобы добавить дополнительные пары точек, определите ориентиры и щелкните по изображениям. Сохраните управляющие точки, выбрав пункт меню Файл (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;

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

[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

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

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]

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

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);

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

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

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

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

Кредиты

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