QA0002
ZIP with notes and MATLAB script:
Q: How to estimate amount of crop rows out of bird eye single photograph.
A: In industrial agriculture there are many modern application to control the amount of trees in plantations.
I once read an application of a drone flying close enough to orange trees to count tree trunks.
Here a single photo is supplied and the amount of rows is needed.
An answer to a similar question is also supplied here:
https://es.mathworks.com/matlabcentral/answers/390241-row-detection-in-crops
In the ZIP pack download from top right Dropbox link there's an application note by R. Marino W. Davis as example of the more efficient foliage penetrating laser application to detect vehicles.
That application shows of the importance to 'peal' the scenario under test, and the surface layer is by no means anyway enough data to obtain sought information. Human eye (visible light band) does not return solid objects underneath tree foliage. An adequate laser is needed to 'comb' a volume at different power levels and that would be a way to accurately locate and count trees trunks, or the ground surface itself. From there one can infer the amount of rows and lines, that often split, merge, or may simply end up far from field edges.
% counting crop rows.
% spotting gaps along crop rows that are wider than given threshold.
clear all;close all;clc
A=imread('001.jpg');
figure(1);imshow(A)
%P001
%% PART I : counting rows on single cross section line
[im_1 im_2 im_3]=size(A)
varA=var(double(A),0,3);
stdA=std(double(A),0,3);
figure(2);h1=surf(varA,'EdgeColor','None');
figure(3);h2=surf(stdA,'EdgeColor','None');
% P002
% remove mean off each colour channel
A10=A(:,:,1)-min(A(:,:,1));
A20=A(:,:,2)-min(A(:,:,2));
A30=A(:,:,3)-min(A(:,:,3));
% picking single colour channels the result of narrowing contrast is quite the same
% narrow contrast on Red-mean(Red)
figure(4);imshow(A10);imcontrast
% manual imcontrast can be automated with commands imadjust or strechlim
A2=A(:,:,2);
A_2=imadjust(A2,[.49;.6],[ 1 ;0]);
figure(6);imshow(A_2)
% A_2 is the start image now. We start to appreciate the gaps along rows that you want to measure
%% let's count bush rows between the 2 thick trenches G-mean(G) along a single cross section
close all;
hf1=figure(1);imshow(A_2)
ax=hf1.CurrentAxes
% NOW choose 2 points between the 2 thick trenches or roads at each side of figure(1):
p3=ginput(2);p3=uint64(floor(p3)) % mouse input 2 points
hold(ax,'on');plot(ax,p3(1,1),p3(1,2),'ro');plot(ax,p3(2,1),p3(2,2),'ro')
plot(ax,p3(:,1),p3(:,2),'r')
%P007-2
py_ref=min(p3(2,2),p3(1,2))
im_line3_1=A(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))],1) % Red
im_line3_2=A(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))],2) % Green
im_line3_3=A(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))],3) % Blue
v3_ref=[1:1:(abs(double(p3(1,1))-double(p3(2,1))))+1]
figure(2);plot(v3_ref,im_line3_1,'r',v3_ref,im_line3_2,'g',v3_ref,im_line3_3,'b');grid on
im_line3_1=255-im_line3_1
im_line3_2=255-im_line3_2
im_line3_3=255-im_line3_3
figure(3);plot(v3_ref,im_line3_1,'r',v3_ref,im_line3_2,'g',v3_ref,im_line3_3,'b')
% P008
% again note that there's need for higher resolution.
% TRY: build mix line for each pix if B>125 get B if B<125 get R
% Green component does not help that much, the rows are green
% but the space between rows is Cyan and easier to count
% so carry on with Red and Blue only, (R+B) and amplify a bit
L1=3*(double(im_line3_1)+double(im_line3_3)) % amplify a bit
nL1=[1:1:length(L1)]
figure(4);plot(L1);grid on
% .- remove DC
L1=L1-mean(L1)
hold on;figure(4);plot(L1);grid on
% 9.- 1st assault: counting peaks, trying different findpeaks threshold, values
% returns different values, not a robust approach
numel(findpeaks(L1,'minpeakheight',80)) = 7
numel(findpeaks(L1,'minpeakheight',100)) = 7
numel(findpeaks(L1,'minpeakheight',140)) = 7
numel(findpeaks(L1,'minpeakheight',170)) = 7
numel(findpeaks(L1,'minpeakheight',190)) = 7
% sweep findpeaks parameter 'minpeakheight'
peaks_prob_1=0
for k=80:1:250
peaks_prob_1=[peaks_prob_1 numel(findpeaks(L1,'minpeakheight',k))];
end
figure(5);plot(peaks_prob_1);grid on
% no need for this figure, fundpeaks sweep is flat, all 7
% ignore the rest of the script, adapting answer to another question, because DTFT
% and FFT need further processing, do not use FFT for such low amount of input data
% Yet the following filter works straight away
% define the filter
b=[0.028 0.053 0.071 0.053 0.028]
a=[1.000 -2.026 2.148 -1.159 0.279]
L1_BPfilt=filter(b,a,L1)
nL1_BPfilt=[1:length(L1_BPfilt)]
hold all;figure;plot(nL1,L1,nL1_BPfilt,L1_BPfilt);grid on
L1_BPfilt_detr=detrend(L1_BPfilt,'linear',[95 289])
% with
abs_fft1=abs(fft(L1_BPfilt_detr))
plot(abs_fft1);grid on
r_L1=xcorr(L1_BPfilt_detr,L1_BPfilt_detr)
r_L1=r_L1([floor(.5*numel(r_L1)):end]) % self correlations are symmetric, we only one one side, taking right figure(2);plot(r_L1);grid on
[pks,locs]=findpeaks(r_L1,'minpeakheight',0)
numel(pks) = 7
To address bifurcations, line crossings, bent lines, cut lines and even varying width lines, one has to
repeat the above detailed cross sectioning for enough cross sections and then build
determine a line centre point (for each cross section)
determine a line width (for each cross section)
and then sweep these points to build the skeleton that allows to control varying field of crop lines.
