В этом примере показано, как указать изображение к наземной системе координат и создать новое изображение, на которое "геоссылаются". Это требует Image Processing Toolbox™ в дополнение к Mapping Toolbox™.
В этом примере все геосправочные данные находятся в той же наземной системе координат, состояние Массачусетса Плоская система (использующий североамериканскую Данную величину 1 983 в модулях метров). Это задает наши "координаты карты". Геосправочные данные включают слой ортобазового адреса образа и векторный дорожный слой.
Набор данных, на который геосошлются, является цифровой воздушной фотографией, покрывающей части деревни Вест-Конкорд, Массачусетс, собранный в начале пружины, 1997.
Слой ортобазового адреса образа структурирован в 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' 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'
Отобразите изображения в их правильных пространственных положениях.
mapshow(baseImage1,cmap1,R1) ax1 = gca; mapshow(ax1,baseImage2,cmap2,R2) title('Map View, Massachusetts State Plane Coordinates'); xlabel('Easting (meters)'); ylabel('Northing (meters)');
Этот шаг использует функции 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
.
Извлеките пары контрольной точки из структур контрольной точки.
[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
Создайте объект привязки 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'
Примените геометрическое преобразование с помощью 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'));
Отобразите зарегистрированное изображение на том же самом (координата карты) оси как мозаики ортобазового адреса образа. Зарегистрированное изображение не полностью заполняет свою ограничительную рамку, таким образом, оно включает заполненные пустым указателем треугольники. Создайте альфа-маску данных, чтобы сделать эти закрашенные области рендерингом как прозрачные.
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;
Можно оценить регистрацию путем рассмотрения функций, которые охватывают и зарегистрированное изображение и ортофото изображения.
Используйте 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
MapCellsReference
| cpstruct2pairs
| fitgeotrans
| im2uint8
| imread
| imref2d
| intrinsicToWorld
| maprasterref
| transformPointsForward
| worldfileread