Этот пример показывает, как указать изображение к наземной системе координат и создать новое изображение, на которое "геоссылаются". Это требует 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');
Считайте worldfiles для этих двух мозаик
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));
Субдискретизируйте изображения фактором 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