本文介紹了：插針式土壤粗糙度儀的數字化步驟，以及具體的Matlab實現代碼。這個小程序寫的挺有意思的，可以實現插針式土壤粗糙度儀的互動式數字化。

下圖是試驗中拍攝一張土壤廓線圖，有興趣的可以拿去練習下。

一、插針式土壤粗糙度儀的數字化步驟

1. 當前目錄設置為圖片的文件夾目錄，運行腳本。首先，在白板上輸入圖片名稱，如下圖所示，輸入名稱：FIELD-2-POINT-3-ALONG-W。然後，將OK轉到步驟2。

2. 輸入兩個點以獲得要放大的範圍。請注意，範圍中應包含四個紅點。

3. 輸入四個紅點的中心位置；請注意，輸入的點將標記為綠色，按「Enter」鍵將結束輸入。然後，程序將使用輸入的位置消除圖像變形。

4. 在重新投影的圖片上，再次輸入四個紅點的中心位置。請注意，輸入的點將標記為綠色，按「Enter」鍵將結束輸入。然後，程序將刪除由四個點定義的範圍之外的圖片邊框。

5. 逐個描畫插針的高度。請注意，輸入的點將標記為紅色，在完成所有點後，按「Enter」鍵結束輸入。然後，程序將以釐米為單位計算每個點的坐標。結果的X間距將插值為1 cm，然後繪製在新圖中。

6. 最後，可以對數字化後的結果進行核驗。在彈出的問題對話框中，您可以選擇「繼續」按鈕來數字化下一張圖片，或者按「結束任務」結束你的工作。

二、土壤粗糙度拍攝圖片的數字化腳本（Matlab）

%% 1.獲取要數字化的圖片
clear, clc
close all

% 圖片所在的文件夾
picturedir = Z:Works_UKPinboard_photos_20180824
oughness_extraction;

% 選擇要數字化的圖片
filename = uigetfile(*.JPG,Pick a file,picturedir);
roughness_picture = importdata([picturedir filename]);

%% 2.獲取有效範圍
screensize = get(0,ScreenSize); % get screen resolution
screensize2 = screensize;
screensize2(1) = screensize(1)+screensize(3);

% 最大化圖片
figure(Name,filename,position,screensize)
imagesc(roughness_picture);

% 輸入採樣點信息
pointName = cell2mat(inputdlg(Please input the point name:,Input,[1 50]));
pointName = matlab.lang.makeValidName(pointName);

% 放大至有效作圖範圍
uiwait(msgbox(Please input two points, then enlarge to that extend));

p = ginput(2);
% the corners of the polygon
sp(1) = min(floor(p(1)), floor(p(2))); %xmin
sp(2) = min(floor(p(3)), floor(p(4))); %ymin
sp(3) = max(ceil(p(1)), ceil(p(2))); %xmax
sp(4) = max(ceil(p(3)), ceil(p(4))); %ymax

% 顯示結果
MM = roughness_picture(sp(2):sp(4), sp(1): sp(3),:);

imagesc(MM);

%% 3.糾正圖片變形
uiwait(msgbox(Please input the central locations of the four red point to wrap the picture!));
% [x_pixs,y_pixs] = ginput;
x_pix4 = [];
y_pix4 = [];
finish=false;
set(gcf,CurrentCharacter,@);

while ~finish
% get the coordinates of the four points
% check for keys
key = get(gcf,CurrentCharacter);
children = get(gca, children);

if key~=@ % has it changed from the dummy character?
set(gcf,CurrentCharacter,@); % reset the character

if strcmp (key, char(13))
% finish=true;
break;
% if input delete, then delete the current point
elseif strcmp(key, char(127)) && length(children)>2

delete(children(1));
delete(children(2));
x_pix4 = x_pix4(1:end-2);
y_pix4 = y_pix4(1:end-2);
continue;
end
end

[x_pix1,y_pix1] = ginput(1);
hold on
plot(x_pix1,y_pix1,*g)

x_pix4 = [x_pix4 x_pix1];
y_pix4 = [y_pix4 y_pix1];
end
hold off

% square sum of the coordinates, and sort the matrix
corSUM = x_pix4.^2+y_pix4.^2;
[corSUM_sort,I] = sort(corSUM);

% exchange (swap) the second and third elements
I([2 3]) = I([3 2]);

% from left-right and up-down
x_pix4 = x_pix4(I);
y_pix4 = y_pix4(I);

topLeft = [x_pix4(1), y_pix4(1)];
topRight = [x_pix4(2), y_pix4(2)];
botRight = [x_pix4(4), y_pix4(4)];
botLeft = [x_pix4(3), y_pix4(3)];

U = [topLeft; topRight; botRight; botLeft];
width = size(MM,2);
height = size(MM,1);
topLeftNew = [1 1];
topRightNew = [width,1];
bottomLeftNew = [1,height];
bottomRightNew = [width,height];
X = double([topLeftNew; topRightNew; bottomRightNew; bottomLeftNew;]);
tform = fitgeotrans(U, X, projective);
MM_wrap = imwarp(MM, tform);

% 顯示糾偏後的結果
imagesc(MM_wrap)

%% 4.裁剪多餘區域
uiwait(msgbox(Please input the central locations of the four red points to cut the picture!));
% get the coordinates of the four points
x_pix4 = [];
y_pix4 = [];

finish=false;
set(gcf,CurrentCharacter,@);

while ~finish
% get the coordinates of the four points
[x_pix1,y_pix1] = ginput(1);
hold on
plot(x_pix1,y_pix1,*g)

x_pix4 = [x_pix4 x_pix1];
y_pix4 = [y_pix4 y_pix1];

% check for keys
key = get(gcf,CurrentCharacter);
children = get(gca, children);

if key~=@ % has it changed from the dummy character?
set(gcf,CurrentCharacter,@); % reset the character
if strcmp (key, char(13))
break;

elseif strcmp(key, char(127)) && length(children)>2
delete(children(1));
delete(children(2));
x_pix4 = x_pix4(1:end-2);
y_pix4 = y_pix4(1:end-2);
continue;
end

end
end
hold off

yMin = int16(min(y_pix4));
yMax = int16(max(y_pix4));
xMin = int16(min(x_pix4));
xMax = int16(max(x_pix4));

MM_new = MM_wrap(yMin:yMax,xMin:xMax,:);
imagesc(MM_new)

%% 5.逐插針描畫其位置
% 獲取1釐米代表多少個像素點
pixels_height = size(MM_new,1)/26;

% pixels cover 1 cm width
pixels_width = size(MM_new,2)/106;

uiwait(msgbox(describe the profile: from 1st to the last bin))

x_pixs = [];
y_pixs = [];

finish=false;
set(gcf,CurrentCharacter,@);

while ~finish

[x_pix,y_pix] = ginput(1);
hold on
plot(x_pix,y_pix,*r)

x_pixs = [x_pixs x_pix];
y_pixs = [y_pixs y_pix];

% check for keys
key = get(gcf,CurrentCharacter);
children = get(gca, children);

if key~=@ % has it changed from the dummy character?
set(gcf,CurrentCharacter,@); % reset the character
if strcmp (key, char(13)),
break;
elseif strcmp(key, char(127)) && length(children)>2
delete(children(1));
delete(children(2));
x_pixs = x_pixs(1:end-2);
y_pixs = y_pixs(1:end-2);
continue;
end
end
end

% 計算對應的長度和高度
x_widths = (x_pixs-x_pixs(1))/pixels_width;
y_heights = -(y_pixs - max(y_pixs))/pixels_height;

x_intp = transpose(0:100);
y_intp = interp1(x_widths,y_heights,x_intp,spline);

%% 6. 確認數字化後的結果
SW0 = screensize(:,3)*46/1920;
SW = screensize(:,3)*1833/1920;
SH0 = screensize(:,4)*735/1080;
SH = screensize(:,4)*201/1080;

figure(Name,filename,position,[SW0 SH0 SW SH])

plot(x_intp,y_intp)

xlabel(Length (cm))
ylabel(Height (cm))

%% 保存結果

dir_result = [picturedir,2018_greece];

if ~isdir(dir_result)

mkdir(dir_result)

end

save([dir_result pointName,.mat],x_intp,y_intp)

%% continue or end?
answer = questdlg(Would you like continue?, ...
Make a choose, ...
Continue,End Task,Continue);

% Handle response
switch answer
case Continue
run roughness_extraction_V2.m;
case End Task
disp(Task End);
end

disp(filename)
disp(pointName)

推薦閱讀：

相關文章