本文介绍了：插针式土壤粗糙度仪的数字化步骤，以及具体的Matlab实现代码。这个小程序挺有意思的，可以实现插针式土壤粗糙度仪的交互式数字化。

下图是试验中拍摄一张土壤廓线图，有兴趣的可以拿去练习下。

一、插针式土壤粗糙度仪的数字化步骤

1. 当前目录设置为图片的文件夹目录，运行脚本。首先，在白板上输入图片名称，如下图所示，输入名称：FIELD-2-POINT-3-ALONG-W。然后，将OK转到步骤2。

2. 输入两个点以获得要放大的范围。请注意，范围中应包含四个红点。

3. 输入四个红点的中心位置；请注意，输入的点将标记为绿色，按“Enter”键将结束输入。然后，程序将使用输入的位置消除图像变形。

4. 在重新投影的图片上，再次输入四个红点的中心位置。请注意，输入的点将标记为绿色，按“Enter”键将结束输入。然后，程序将删除由四个点定义的范围之外的图片边框。

5. 逐个描画插针的高度。请注意，输入的点将标记为红色，在完成所有点后，按“Enter”键结束输入。然后，程序将以厘米为单位计算每个点的坐标。结果的X间距将插值为1 cm，然后绘制在新图中。

6. 最后，可以对数字化后的结果进行核验。在弹出的问题对话框中，您可以选择“继续”按钮来数字化下一张图片，或者按“结束任务”结束你的工作。

二、土壤粗糙度拍摄图片的数字化脚本（Matlab）

%% 1.获取要数字化的图片clear, clcclose all
 % 图片所在的文件夹picturedir = 'Z:\Works_UK\Pinboard_photos_20180824\roughness_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 polygonsp(1) = min(floor(p(1)), floor(p(2))); %xminsp(2) = min(floor(p(3)), floor(p(4))); %yminsp(3) = max(ceil(p(1)), ceil(p(2))); %xmaxsp(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];endhold off
 % square sum of the coordinates, and sort the matrixcorSUM = 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-downx_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 pointsx_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
  
 endendhold 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 widthpixels_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
 endend
 % 计算对应的长度和高度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 responseswitch answer
 case 'Continue'
 run roughness_extraction_V2.m;
 case 'End Task'
 disp('Task End');end
 disp(filename)disp(pointName)