Row detection in crops

I am trying to detect rows (as a line) in georeferenced ortophotos. I am able to do it by computing hough transform, but only in images without coordinates and what I need is to extract a georeferenced vector file consisting on the different rows from the images.

Any ideas?

I attach a sample of the image, but if you want the whole ortophoto send me a private message and I can send it to you.

Thanks!

[Image Analyst]

Somehow the image is no longer attached. But this is not the 1 answer - just a comment - so see John's Answer below.

Answer 1

[John BG]

UNDER CONSTRUCTION - REACHED LIMIT AMOUNT UPLOADS

Hi Siscu

this is John BG <mailto:[e-mail] [e-mail]>

Allow me to answer your crop measurements question in 2 parts.

1st part the count of amount of rows between 2 points with enough accuracy is defined. Different approaches are shown.

In the second part a method is suggested to determine the positions and lengths of all those gaps along each line that is wider than the required threshold of 1.3m.

It's important to bear in mind that higher accuracy of the pictures would allow to actually detect each plant / bush / tree, and therefore it would be far easier to actually measure everything.

PART I : Counting rows on single cross section line.

1.-

Acquiring image

I have used the .jpg copy of your sample.

A=imread('001-1.jpg');
figure(1);imshow(A)

%P001

If the image had more resolution one would correlate 2D the entire image against small images of trees and bushes to find, but it's not the case, is it?

2.-

One way to count crop lines is through the variance

[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

If you walk across the crop lines, through the gaps the Blue is always at least 20 points below those local values of Read and Green. But the variance and standard deviation look quite noisy, I prefer the FFT of the walk perpendicular across the crop lines, let me explain:

3.-

Choose the best colour channel or combination of them, it depends upon the input picture.

% 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
% P003
% P004
% applying a narrow contrast one can get the following
%P005

%  narrow contrast on Green-mean(G)
% imshow(A20);imcontrast
%  narrow contrast on Blue-mean(B)
% imshow(A30);imcontrast
% colour channels without removing respective mean seem to be a bit more solid.
% A1=A(:,:,1);
% A2=A(:,:,2);
% A3=A(:,:,3);
%  narrow contrast on Red
% imshow(A1);imcontrast
% quite the same as Red
%  narrow contrast on Green
% imshow(A2);imcontrast
%  narrow contrast on Blue
% imshow(A3);imcontrast
% R+G=Y
% R+B=M
% B+G=C
% Y=A;Y(:,:,3)=0;
% M=A;Y(:,:,2)=0;
% C=A;Y(:,:,1)=0;
%  narrow contrast on Red-mean(Red)
% Y2=rgb2gray(Y);imshow(Y2);imcontrast
% quite the same
%  narrow contrast on Green-mean(G)
% M2=rgb2gray(M);imshow(M2);imcontrast
%  narrow contrast on Blue-mean(B)
% C2=rgb2gray(C);imshow(C2);imcontrast

4.-

Manual imcontrast can be automated with imadjust

 either with strechlim

A_2=imadjust(A,stretchlim(A),[]);
figure(5);imshow(A_2)
% P006
% A_3=imadjust(A,stretchlim(A_2),[]);figure;imshow(A_3)

or directly telling the contrast limits

%  J=imadjust(I,[low_in;high_in],[low_out;high_out])

in previous imcontrast parameters screen shot the narrow filter has shoulders at .49 and .6

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.

5.-

Let's count crop 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

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 it's important to mention that there's serious need for higher resolution.

With more pixels per object, one would be able to correlate the entire image to single bush images and exactly allocate each shoot and root, making every thins else easier.

6.-

Less is more

Green component does not help that much, the rows are dark green but the spaces between rows are kind of darkish Cyan, easier to count.

So we carry on with Red and Blue only, (R+B). Amplify a bit

L1=3*(double(im_line3_1)+double(im_line3_3)) % amplify a bit
nL1=[1:1:length(L1)]
figure(8);plot(L1);grid on

7 .-

Remove DC

L1=L1-mean(L1)
hold on;figure(9);plot(L1);grid on

8.-

1st assault: counting peaks, trying different findpeaks threshold, values returns different values, not a robust approach

numel(findpeaks(L1,'minpeakheight',80))
% =
%    199

numel(findpeaks(L1,'minpeakheight',100))
%  =
%    185

numel(findpeaks(L1,'minpeakheight',140))
% =
%    114
%    
numel(findpeaks(L1,'minpeakheight',170))

%  =
%     60

numel(findpeaks(L1,'minpeakheight',190))

%  =
%     35

% 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
% P011

I manually counted 141.

There is already a slightly flat section between ~142 and 152, where the sought result 141 lays, counted manually.

So the valid result already shows up but there should be no doubt what value to go for.

Such flat run is too narrow.

v_prob_1=[1:length(peaks_prob_1)]

figure(6);hist(peaks_prob_1,v_prob_1)

% P012

Still misleading: the histogram shows 1st peak on 166, 2nd and 3rd peaks have same weight 160 and 147.

4th and 5th peaks have values 61 and 40 while sharing same weight as % 2nd and 3rd, which is again misleading, there are far many more rows than 40 or 61.

9.-

Let's use DFT:

Spectrum analysis needs as many useful samples as possible, as well as the data having to show some kind of cycle or repetition.

nL1=[1:length(L1)] % just copied lines from another script with DTFT
x=L1;n=[1:1:length(L1)]
x=double(x)
k=0:500;w=(pi/500)*k
X=x*(exp(-j*pi/500)).^(n'*k)
X=x*(exp(-j*pi/500)).^(n'*k)
magX=abs(X);angX=angle(X);realX=real(X);imagX=imag(X)
figure(15);subplot(2,2,1);plot(k/500,magX);grid
xlabel('freq[rad]');title('Magnitude(X)')
figure(15);subplot(2,2,3);plot(k/500,angX/pi);grid
xlabel('freq[rad]');title('Phase(X)')
figure(15);subplot(2,2,2);plot(k/500,realX);grid
xlabel('freq[rad]');title('Real(X)')
figure(15);subplot(2,2,4);plot(k/500,imagX);grid
xlabel('freq[pi]');title('Imag(X)')

% P013

counting cycles:

2 px/cycle is max frequency of DFT

amount_f_max_cycles_between_p3_points=(p3(2,1)-p3(1,1))/2
DFT_fundamental=.112
rows_count =DFT_fundamental*amount_f_max_cycles_between_p3_points
% =  
% 138

X.-

what about with FFT?

NFFT=2^nextpow2(length(L1)) % 512 freq points [0 2pi]
absFFT_L1=abs(fft(L1,NFFT))
figure(5);plot(absFFT_L1);grid on

% P014

nfmax=length(absFFT_L1)/2
nf1=233  % marker on abs(FFT) graph
nf1/nfmax

Again: 2 consecutive pixels of image per cycle constitute the highest frequency possible.

rows_count=amount_f_max_cycles_between_p3_points*nf1/nfmax
% 
% = 140

The expected result, is 141, only for this particular cross section.

11.-

Now, let's have a look at time domain analysis, let's count the peaks of the cross section sample

command FINDPEAKS and parameter 'minpeakprominence'

A single shot with findpeaks may be too narrow of a measurement

[pks,locs]=findpeaks(L1,'minpeakprominence',60);figure(7);plot(nL1,L1,locs,pks+4,'kv');grid on;numel(pks) 
%=31

the length of peaks converges to an constant amount

pks_range={};
for k=10:1:220
    [pks_,locs]=findpeaks(L1,'minpeakprominence',k);
    pks_range=[pks_range pks_];
end

% P015

lengths_findpeaks=0
for k=1:1:size(pks_range,2)
    L_=pks_range{k};
    lengths_findpeaks=[lengths_findpeaks length(L_)];
end
lengths_findpeaks(1)=[];
hf20=figure;hhist=histogram(lengths_findpeaks);grid on
hhist.BinEdges=[10:1:221]-.5;
hf20.CurrentAxes.XLim=[110 220];
hhist_y=hhist.Values;
hhist_x=hhist.Data;
rows_count=hhist_x(find(hhist_y==max(hhist_y)))
% =
% 137

 not as accurate as the FFT, and fast at the same time, winning combination.

% P016

12.-

Now frequency contents with command PLOMB

Fs=1;Fmax=1/2;[Pxx,f]=plomb(double(L1),Fs,Fmax,10,'power');figure(7);plot(Pxx);grid on;

% P017

[Pxx_pks,f0]=findpeaks(Pxx,f,'MinPeakHeight',400)

% Pxx_pks =
%    1.0e+03 *
%    0.596695269876198
%    8.345979694883447
%    1.401890369308839
% f0 =
%    0.056034132466477
%    0.056603006907761
%    0.057212515237708

the peak of interest is:

f0(2)
% =
% f0=0.056603006907761

Now there's no need to halve the length of the spectrum because it's one sided spectrum, no mirror.

To this particular purpose, do not use f0(2), but

f_peak=find(Pxx==max(Pxx))/length(Pxx)
amount_f_max_cycles_between_p3_points=(p3(2,1)-p3(1,1))/2
rows_count=amount_f_max_cycles_between_p3_points*f_peak
% =
% 140

13.-

telling plomb() what Fs and Fmax to use:

[Pxx,f]=plomb(double(L1),nL1);figure(7);plot(Pxx);grid on

% P018

% marker on f=123, started all over with different line, now length(L1)=379 % tried with different windows, Welch for instance broadens peak of interest, and the frequency spike % now it's sharper

f_peak=find(Pxx==max(Pxx))/length(Pxx)
amount_f_max_cycles_between_p3_points=(p3(2,1)-p3(1,1))/2
rows_count=amount_f_max_cycles_between_p3_points*f_peak
% =
% 139, still FFT seems to be best option

14.-

Correlate signal with expected pulse, have a look at a single received pulse try different pulses, gradually make the pulse narrower and taller

pulse1=[0 200 200 200 200 200 200 200 200 200 200 200 0] % poor correlation
pulse2=[0 50 100 100 100 100 100 100 100 100 100 50 0] % a bit better
pulse3=[0 0 50 100 100 100 100 100 100 100 50 0 0] % a bit better, just missing 1 pulse
pulse4=[0 0 0 50 100 100 100 100 100 50 0 0 0]
pulse5=[0 0 0 100 200 200 200 200 200 100 0 0 0]
pulse6=[0 0 0 0 100 200 200 200 100 0 0 0 0] % the count of peaks is now correct, do pulse1=[0 0 0 0 200 400 400 400 200 0 0 0 0]
pulse7=[0 0 0 0 0 200 400 200 0 0 0 0 0] 
pulse8=[0 50 0 -200 0 800 1200 800 0 -200 0 50 0] 

pulse_test=[pulse1;pulse2;pulse3;pulse4;pulse5;pulse6;pulse7;pulse8];
% pulse_test=[pulse1;pulse7;pulse8]

for k=1:1:size(pulse_test,1)
  r_L1_pulse1=conv(L1,pulse_test(k,:))
  hold on;figure(5);plot(r_L1_pulse1);grid on
end
grid on

% P019

% the blue trace is the enhanced one

% if you manually check with findpeaks, at every next pulse there is a bit improvement. % Yet, depending on chosen sample line, there are still missed rows. % the closer to a sin(x) or cos(x) that the pulse signal looks like, which % is what I tend to do along the test pulses, the chances to catch all rows % increases, yet, if a pulse resembling many sin(x) cycles is going to be % use, does it make sense to fo straight to the FFT? after all such % transform is a correlation with a comb of sin(n*x) signals

prob2=[]
for k=1:20
    [pks_xcorr1puls,locs_xcorr1puls]=findpeaks(r_L1_pulse1,'minpeakdistance',k)
    prob2=[prob2 numel(pks_xcorr1puls)]
end
% figure(9);plot(prob2);grid on
[pks_xcorr1puls,locs_xcorr1puls]=findpeaks(r_L1_pulse1,'minpeakdistance',6)
% remove start and end counts
numel(pks_xcorr1puls)-2 % -2, as mentioned, for first and last of previous picture

1st half of r_L1_pulse1 does not have relevant information

r_L1_pulse1=r_L1_pulse1([375:750])
n_r_L1_pulse1=[1:1:length(r_L1_pulse1)]

spectrum of the autocorrelation

[Pxcorr,f]=plomb(r_L1_pulse1,n_r_L1_pulse1);figure(10);plot(Pxcorr);grid on

% P020

f_peak=find(Pxcorr==max(Pxcorr))/length(Pxcorr)
amount_f_max_cycles_between_p3_points=(p3(2,1)-p3(1,1))/2
rows_count=amount_f_max_cycles_between_p3_points*f_peak
% =
% 141, same

PART II : measure positions of gaps > 1.3m along each row

close all;
hf1=figure(1);imshow(A_2)
ax=hf1.CurrentAxes

disp('key 2 points, 1st up 2nd down, that define left side vertical limiting row')
P_left=ginput(2);P_left=uint64(floor(P_left)) % mouse input 2 points 
hold(ax,'on');plot(ax,P_left(1,1),P_left(1,2),'bo','LineWidth',2);plot(ax,P_left(2,1),P_left(2,2),'bo','LineWidth',2)
plot(ax,P_left(:,1),P_left(:,2),'b','LineWidth',2)

disp('key 2 points, 1st up 2nd down, that define right side vertical limiting row')
P_right=ginput(2);P_right=uint64(floor(P_right)) % mouse input 2 points 
hold(ax,'on');plot(ax,P_right(1,1),P_right(1,2),'bo','LineWidth',2);plot(ax,P_right(2,1),P_right(2,2),'bo','LineWidth',2)
plot(ax,P_right(:,1),P_right(:,2),'b','LineWidth',2)

hold(ax,'on');plot(ax,[P_left(1,1) P_right(1,1)],[P_left(1,2) P_right(1,2)],'r-','LineWidth',2)
hold(ax,'on');plot(ax,[P_left(2,1) P_right(2,1)],[P_left(2,2) P_right(2,2)],'r-','LineWidth',2)

%P2_001

% L_1=linspace(P_left(1,1),P_right(1,1),)

[John BG]

Siscu

this is John BG <mailto:[e-mail] [e-mail]>

I was wondering that perhaps given the length of the above 1st part, would it be ok to you if the 2nd part is posted in another question, just for the sake of saving scroll-down time to the readers.

Please let me know, thanks in advance

John BG

<mailto:[e-mail] [e-mail]>

[ssisscu Garcia]

You can find it here: https://es.mathworks.com/matlabcentral/answers/393650-row-detection-in-crops-2