-
Notifications
You must be signed in to change notification settings - Fork 0
/
process_image_ndvi.m
65 lines (44 loc) · 1.6 KB
/
process_image_ndvi.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
function [floresta_km2, floresta_pct] = process_image_ndvi(im, area_km2)
im = im2double(im);
s = size(im);
pixel_count = s(1) * s(2);
pixel_area_km2 = area_km2 / pixel_count;
figure, imshow(im), title("Imagem original");
%% NDVI
nir = im(:,:,1);
red = im(:,:,2);
ndvi = (nir - red) ./ (nir + red);
%% Autocontraste
ndvi = imadjust(ndvi);
figure, imshow(ndvi);
% , title("Normalized Difference Vegetation Index (NDVI)");
%% Filtro gaussiano
h = fspecial("gaussian", 6, 0.8);
ndvi = imfilter(ndvi, h);
% figure, imshow(ndvi), title("Após filtro gaussiano");
%% K-Means
K = 3;
% reorganiza pra passar pro k-means
nr = size(ndvi,1);
nc = size(ndvi,2);
np = nr*nc;
ndvi_cols = reshape(ndvi, np, 1);
% K-means
labeled = kmeans(ndvi_cols, K, 'Distance', 'sqeuclidean');
labeled = reshape(labeled, nr, nc);
labeled_rgb = labeloverlay(im, labeled, "Transparency", 0.8, "Colormap", "jet");
figure, imshow(labeled_rgb), title('Imagem segmentada. Clique em uma área de floresta')
point = drawpoint('Color', "green");
point.Label = "Floresta";
y = floor(point.Position(2));
x = floor(point.Position(1));
delete(point)
label_floresta = labeled(y, x);
%% Calcular area
areas = zeros(K);
for label = 1:K
areas(label) = sum(labeled == label, "all") * pixel_area_km2;
end
floresta_km2 = areas(label_floresta);
floresta_pct = (floresta_km2 / area_km2) * 100;
end