Короче говоря, вы получите очень плохую оценку размера, используя ваш метод. Края размытых областей, независимо от того, какое пороговое значение вы выберете, не соответствуют краям инселбергов, которые вы хотите измерить. Вместо этого я предлагаю вам следовать рецепту ниже. Я использую DIPlib в Python для этого (отказ от ответственности: я автор). Привязки Python представляют собой тонкий слой в библиотеке C ++, довольно просто перевести приведенный ниже код Python в C ++ (мне просто проще разрабатывать его интерактивно в Python).
I загрузили исходные данные высоты по предоставленной вами ссылке (а не по склону). DIPlib может напрямую читать TIFF-файлы с плавающей запятой, поэтому нет необходимости в специальном преобразовании. Я обрезал область, аналогичную той, которая используется OP для этой демонстрации, но нет причин не применять метод ко всей плитке.
import PyDIP as dip
height = dip.ImageRead('17S42_ZN.tif')
height.SetPixelSize(0.000278, 'rad') # not really radian, but we don't have degrees
height = height[3049:3684, 2895:3513];
Код также устанавливает размер пикселя в соответствии с данными в файл TIFF (в единицах радиана, потому что DIPlib не делает градусов).
height">
Затем я применяю фильтр top-hat с конкретный c диаметр (25 пикселей). Это изолирует все пики диаметром 25 пикселей или менее. Отрегулируйте этот размер в соответствии с тем, какая максимальная ширина должна быть у Инсельберга.
local_height = dip.Tophat(height, 25)
В результате получается локальная высота, высота над некоторой базовой линией, определяемая размером фильтра.
local_height">
Далее я применяю порог гистерезиса (двойной порог). Это дает бинарное изображение с порогом на 100 м выше базовой линии, где ландшафт поднимается выше 200 м выше этой базовой линии. То есть я решил, что инсельберг должен быть как минимум на 200 м выше базовой линии, но отрезал каждый из них на расстоянии 100 м. На этой высоте мы будем измерять размер (площадь). Опять же, отрегулируйте пороги так, как считаете нужным.
inselbergs = dip.HysteresisThreshold(local_height, 100, 200)
inselbergs">
Теперь все, что осталось, это измерение найденных областей:
labels = dip.Label(inselbergs)
result = dip.MeasurementTool.Measure(labels, features=['Size', 'Center'])
print(result)
Это выводит:
| Size | Center |
-- | ---------- | ----------------------- |
| | dim0 | dim1 |
| (rad²) | (rad) | (rad) |
-- | ---------- | ---------- | ---------- |
1 | 1.863e-05 | 0.1514 | 0.01798 |
2 | 4.220e-05 | 0.1376 | 0.02080 |
3 | 6.214e-05 | 0.09849 | 0.04429 |
4 | 6.492e-06 | 0.1282 | 0.04710 |
5 | 3.022e-05 | 0.1354 | 0.04925 |
6 | 4.274e-05 | 0.1510 | 0.05420 |
7 | 2.218e-05 | 0.1228 | 0.05802 |
8 | 1.932e-05 | 0.1420 | 0.05689 |
9 | 7.690e-05 | 0.1493 | 0.06960 |
10 | 3.285e-05 | 0.1120 | 0.07089 |
11 | 5.248e-05 | 0.1389 | 0.07851 |
12 | 4.637e-05 | 0.1096 | 0.09016 |
13 | 3.787e-05 | 0.07146 | 0.1012 |
14 | 2.133e-05 | 0.09046 | 0.09908 |
15 | 3.895e-05 | 0.08553 | 0.1064 |
16 | 3.308e-05 | 0.09972 | 0.1143 |
17 | 3.277e-05 | 0.05312 | 0.1174 |
18 | 2.581e-05 | 0.07298 | 0.1167 |
19 | 1.955e-05 | 0.04038 | 0.1304 |
20 | 4.846e-05 | 0.03657 | 0.1448 |
(Помните, где говорится, что "rad" - это действительно градусы.) Площадь в квадратных градусах немного странная, но вы можете преобразовать ее в квадратные метры. так как вы знаете местоположение на земном шаре. На самом деле может быть проще перевести размеры пикселей в метры перед вычислениями.
Значения, приведенные здесь для «Центра», относятся к верхнему левому пикселю, если мы не обрезали плитку до Начнем с того, что мы могли бы добавить координаты для плитки (как можно получить из соответствующего тега в файле TIFF): (-42.0, -17.0).
В C ++ должен выглядеть код как это:
#include <diplib/simple_file_io.h>
#include <diplib/morphology.h>
#include <diplib/segmentation.h>
#include <diplib/regions.h>
#include <diplib/measurement.h>
//...
dip::Image height = dip::ImageRead("17S42_ZN.tif");
height.SetPixelSize(0.000278 * dip::Units::Radian());
height = height.At(dip::Range(3049, 3684), dip::Range(2895, 3513));
dip::Image local_height = dip::Tophat(height, 25);
dip::Image inselbergs = dip::HysteresisThreshold(local_height, 100, 200);
dip::Image labels = dip::Label(inselbergs);
dip::MeasurementTool measurementTool;
dip::Measurement result = measurementTool.Measure(labels, {}, {"Size", "Center"});
std::cout << result;