Row detection in crops

2 次查看（过去 30 天）
ssisscu Garcia
ssisscu Garcia 2018-3-23
评论： ssisscu Garcia 2018-4-9
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!
  1 个评论
Image Analyst
Image Analyst 2018-3-28
Somehow the image is no longer attached. But this is not the 1 answer - just a comment - so see John's Answer below.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John BG
John BG 2018-3-26
编辑：John BG 2018-4-2
UNDER CONSTRUCTION - REACHED LIMIT AMOUNT UPLOADS
Hi Siscu
this is John BG <mailto:jgb2012@sky.com jgb2012@sky.com>
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),)
  2 个评论
John BG
John BG 2018-4-8
Siscu
this is John BG <mailto:jgb2012@sky.com jgb2012@sky.com>
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:jgb2012@sky.com jgb2012@sky.com>

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Estimation 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by