How to get line count from an image?

3 次查看（过去 30 天）
Shariful Islam
Shariful Islam 2016-5-13
评论： Image Analyst 2018-1-19
I have UAV tree image as attached herewith. Now I want to count the number of rows of trees from that image, How can i do this?
It is an ORCHARD image (part). Here are two types of rows. Rows are the Longest lines that adds trees. At right part (of image) rows have about 60 degrees slope(about 67 rows) and at left part rows have about 120 degrees slope(about 14 rows). You can include isolated trees into the lines at which they can intersect- just think it is part of an orchard. I have thousand of images like this so it needs automation.
I am interested in this row count since i want to know how many rows of trees degraded after some time.
Counting tree is another process of check to measure this on which i am working.
Original image-
Some sample rows (RED lines)-

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

回答（5 个）

John BG
John BG 2016-5-19
编辑：John BG 2016-5-26
ANSWER LINES COUNT - PART II
In this part II different edge detection functions are tested, in an effort to prepare the image for tree lines counting, counting lines always along a line of choice.
BW=imread('lines_count00.jpg')
WHAT WORKS:
1.- Image Processing Toolbox edge with canny,
BW2=edge(BW(:,:,1),'canny'); imshow(BW2);
this edge function is different than the not the edge function in the tatistics and Machine learning Toolbox
2.- edge with log
BW2 = edge(BW(:,:,1), 'log');
imshow(BW2);
3.- canny.m from Mathworks file exchange, by Gagan (2011) http://uk.mathworks.com/matlabcentral/fileexchange/30621-canny-edge-detection/content/canny.m
BW2=canny('lines_count00.jpg',4,.5,.25);
imshow(BW2)
4.- canny2d.m from Mathworks file exchange, by Deshan Yang (2012) http://uk.mathworks.com/matlabcentral/fileexchange/38924-canny-edge-detection-2d/content/canny2d.m
BW2=canny2d('lines_count00.jpg',4,.5,.5);
imshow(BW2)
BW2=canny2d('lines_count00.jpg',4,.5,.25);
imshow(BW2)
5.- canny_v0.m from Mathworks file exchange, by Rachmawan (2014) http://uk.mathworks.com/matlabcentral/fileexchange/46859-canny-edge-detection/content/Canny/Canny_v0.m
WHAT MAY WORK
imgradientxy on red layer
[Gx, Gy] = imgradientxy(BW(:,:,1));
[Gmag, Gdir] = imgradient(Gx, Gy);
figure; imshowpair(Gmag, Gdir, 'montage'); axis off;
title('Gradient Magnitude, Gmag (left), and Gradient Direction, Gdir (right), using Sobel method')
figure; imshowpair(Gx, Gy, 'montage'); axis off;
title('Directional Gradients, Gx and Gy, using Sobel method')
WHAT DOES NOT WORK:
1.-
BW2 = edge(BW(:,:,1), 'roberts');
imshow(BW2);
2.-
BW3=edge(BW(:,:,1),'sobel');
imshow(BW3);
3.- edge on the sum of RGB layers
BW3=edge(BW(:,:,1)+BW(:,:,2)+BW(:,:,3),'sobel');
imshow(BW3);
4.- edge on the CCIR pondered sum of RGB layers
BW3=edge(299*BW(:,:,1)+.587*BW(:,:,2)+.114*BW(:,:,3),'sobel'); imshow(BW3);
5.- bwboundaries, with all the tree lines and only an outer frame added
[B,L,N,A]=bwboundaries(BW(:,:,1))
figure(3), imshow(B)
figure; imshow(BW); hold on;
% Loop through object boundaries
for k = 1:N
% Boundary k is the parent of a hole if the k-th column
% of the adjacency matrix A contains a non-zero element
if (nnz(A(:,k)) > 0)
boundary = B{k};
plot(boundary(:,2),boundary(:,1),'r','LineWidth',2);
% Loop through the children of boundary k
for l = find(A(:,k))'
boundary = B{l};
plot(boundary(:,2),boundary(:,1),'g','LineWidth',2);
end
end
end
this picasso doesn't need a frame
6.- bwboundaries with noholes, another modernist pointless painting
[B,L,N,A]=bwboundaries(BW(:,:,1),'noholes');
imshow(label2rgb(L, @jet, [.5 .5 .5]))
hold on
for k = 1:length(B)
boundary = B{k};
plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)
end
7.- this does not work either:
BW3 = bwboundaries(BW(:,:,1));imshow(BW);hold on;visboundaries(B);
8.- imgradient on sum of RGB layerrs
[Gmag, Gdir] = imgradient(BW(:,:,1)+BW(:,:,2)+BW(:,:,3),'prewitt');
figure; imshowpair(Gmag, Gdir, 'montage');
title('Gradient Magnitude, Gmag (left), and Gradient Direction, Gdir (right), using Prewitt method');axis off;
doesn't work as well as with the coins example, imgradient without option takes sobel by default, same result as the 'prewitt' gradient, similar results when using options 'central' and 'intermediate'
9.- imgradient on single layer
[Gx, Gy] = imgradientxy(BW(:,:,1),'prewitt');
figure; imshowpair(Gx, Gy, 'montage'); axis off;
10.-
[Gx, Gy] = imgradientxy(BW(:,:,1),'intermediate');
figure; imshowpair(Gx, Gy, 'montage');axis off;
11.- canny.m (from Mathworks file exchange) with too high threshold http://uk.mathworks.com/matlabcentral/fileexchange/30621-canny-edge-detection/content/canny.m
BW2=canny('lines_count00.jpg',4,.5,.5);
imshow(BW2)
BW2=canny('lines_count00.jpg',4,.5,.75);
imshow(BW2)
12.- canny2d.m (from Mathworks file exchange) with too high threshold http://uk.mathworks.com/matlabcentral/fileexchange/38924-canny-edge-detection-2d/content/canny2d.m
BW2=canny2d('lines_count00.jpg',4,.5,.75);
imshow(BW2)
13.- the standard starter example to detect cell contours here doesn't work because after inflating, and even without the inflating step, the filling renders the image useless
% 1.- load image
BW = imread('lines_count00.jpg');
figure(1); imshow(BW); % title('original image');
% 2.- detect entire cell
BW2 = edge(BW(:,:,1), 'canny');
imshow(BW2)
threshold=.5;
fudgeFactor = .5;
BWs = edge(BW(:,:,1),'canny', threshold * fudgeFactor);
figure, imshow(BWs), title('binary gradient mask');
% 3.- dilate image
se90 = strel('line', 3, 90);
se0 = strel('line', 3, 0);
BWsdil = imdilate(BWs, [se90 se0]);
figure, imshow(BWsdil), title('dilated gradient mask');
% 4.- fill interior gaps
BWdfill = imfill(BWsdil, 'holes');
figure, imshow(BWdfill);title('binary image with filled holes');
% 5.- remove connected objects on border
% BWnobord = imclearborder(BWdfill, 4);
% figure, imshow(BWnobord), title('cleared border image');
BWnobord = imclearborder(BWsdil, 4);
figure, imshow(BWnobord), title('cleared border image');
% 6.- smoothen image
seD = strel('diamond',1);
BWfinal = imerode(BWnobord,seD);
BWfinal = imerode(BWfinal,seD);
figure, imshow(BWfinal), title('segmented image');
seD = strel('diamond',1);
BWfinal = imerode(BWs,seD);
BWfinal = imerode(BWfinal,seD);
figure, imshow(BWfinal), title('segmented image');
% 7.- optional, place outline around segmented cell
BWoutline = bwperim(BWfinal);
Segout = I;
Segout(BWoutline) = 255;
figure, imshow(Segout), title('outlined original image');
14.- Laplacian with del2 on B3
A=imread('lines_count00.jpg');
A2=(A(:,:,1)+A(:,:,2)+A(:,:,3));
resamp = makeresampler({'nearest','cubic'},'fill');
stretch = maketform('affine',[1.8 0; 0 1.8; 0 0]);
B3 = imtransform(A2,stretch,resamp);
del2 on single RGB layer
del2 on average (A(:,:,1)+A(:,:,2)+A(:,:,3))/3 of all RGB layers without interpolation, going South
insert
  3 个评论
显示 1更早的评论
Shariful Islam
Shariful Islam 2016-6-2
编辑：Shariful Islam 2016-6-2
@ImageAnalyst HaHa- I am following (and upvoted all) all these and finding if these answers my analytical process. BTW it posited many new techniques! But i am thundered a prominent software with larger community experts can not solve this!- this the feeling in my gut as newbie in this field
John BG
John BG 2016-6-2
I am happy for both supplying a useful answer to Shariful, and having the attention of the expert Image Analyst. I was exploring all the possibilities.
It's my opinion that when answering these kind of questions, most of people find it useful to have the reasoning disclosed rather than posting a compact telegraphic answer.
Regards

请先登录，再进行评论。

John BG
John BG 2016-5-16
编辑：John BG 2016-5-23
ANSWER LINES COUNT - PART III
1.- Initially the RGB image has been translated into single layer with
A(:,:,1)+A(:,:,2)+A(:,:,3)
Some may consider that the black & white CCIR correction may be a better option
A42=.299*A(:,:,1)+.587*A(:,:,2)+.114*A(:,:,3); % CCIR
figure(3);imshow(42);
2.- threshold with imcontrast
figure(2);h1=imshow(A(:,:,1)+A(:,:,2)+A(:,:,3));imcontrast(h1)
figure(1);h2=imshow(A(:,:,1));imcontrast(h2)
since h2 the output handle of imcontrast when called as h2=imcontrast(h1) does not point to the resulting image, (it points to the imconstrast display, but what percentage of imcontrast users really need a handle to the GUI compared to the processed image?) the second contrast windowing with [146 147] is equivalent to for instance imcontrast [235 237] on the initial A(:,:,1)+A(:,:,2)+A(:,:,3)
3.- to actually obtain the resulting image either multithresh and imquantize
% A2=A(:,:,1)+A(:,:,2)+A(:,:,3)
level = multithresh(A2);
or
level = multithresh(A2,2);
both show quite same image
seg_A2=imquantize(A2,level);
figure(3);imshow(seg_A2,[])
4.- or simply use imshow that for this particular case it gets to same place as imcontrast, multithresh or quantize. I tried
figure(4);imshow(A2,[135 137])
figure(5);imshow(A42,[125 126])
figure(6);imshow(A42,[135 137])
5.- focus on single RGB layer and narrow image to ribbon: In this case less data is more information. So instead of using all RGB layers, the red layer alone seems to have more contrast than using all the available data. Observe the central patch, where in the original image looks like a Cesna crashed, no trees: with A(:,:,1) such patch is pitch black while using all 3 RGB layers, either with CCIR correction or adding them up, image noise in that area increases.
Let's focuss on left hand side ribbon, just before the Cesna crash like hole starts, and note that we are all the time counting shadows, not the tree lines themselves as their shadows are more consistent:
D=imread('ribbon.jpg')
[sz1D sz2D sz3D]=size(D)
sz1D =
484.00
sz2D =
165.00
sz3D =
3.00
resamp = makeresampler({'nearest','cubic'},'fill');
stretch = maketform('affine',[1.8 0; 0 1.8; 0 0]);
D2 = imtransform(D(:,:,1),stretch,resamp);
[sz1D2 sz2D2 sz3D2]=size(D2)
sz1D2 =
871.00
sz2D2 =
297.00
sz3D2 =
1.00
figure(2);imshow(D2);
6.- Simplify: ignore tree lines tilt angle. A single line is
L1=D2(:,20);
From the plot of the 1st sample line inverted (part I), it seems obvious that if the sample line (the count line running perpendicular to the tree lines), has a strong enough Dynamic range http://info.adimec.com/blogposts/bid/102325/Dynamic-Range-DNR-and-Signal-to-Noise-Ratio-SNR-for-CCD-and-CMOS-image-sensors
the angle between the sample line and the tree lines (as long as 'far enough' from radon rays running parallel ~ theta>45º) doesn't really matter.
7.- radon: Since R2016 release, the Image Processing Toolbox is included in MATLAB base pack, which is great news, because in the past some ocasional MATLAB users seeking image processing solutions, a one-off visit to MATLAB, were put off by not being able to use IP Toolbox functions on the MATLAB package, any one has a chance to test drive the IP Toolbox.
radon basics:
In this case there are 2 ways to use radon • parallel to the tree lines, theta~80º • perpendicular to the tree lines, theta~-20º
the variation (jitter) of spacing between adjacent 'tree lines' is enough to spoil radon on the narrower ribbon, have a look, run the following lines, and press for instance Enter until all angles plotted. Place the plot figure besides MATLAB window so you can see what happens:
theta=[70:.1:80];
[R,xp]=radon(D2,theta)
[sz1R sz2R]=size(R)
for k=1:sz2R
plot(R(:,k));grid on
disp(k);
str=input(' ','s');
end
And this was on the narrow ribbon. No way to get a clear comb with one and only one sharp peak for each tree line.
radon applied to the main image is worse, no matter colour, BW, BW CCIR, single RGB layer: while you may catch a long burst of pulses, there are long enough sections of the sample line that render peak counting useless. Have a look:
8.- radon on ribbon with horizontal tree lines:
We can appreciate that while the top tree lines have slight positive slope, the tree lines on the lower portion of the ribbon have negative slope, the flying platform carrying the camera was probably turning left when the photo was taken.
Really, radon seems an extremely useful function, but in this particular case, at least to me, it make more sense to count peaks along single lines running perpendicular to the tree lines.
9.- radon on even narrower ribbon:
theta=[85:.1:95];
P=imread('ribbon_02_narrower.jpg');
[R,xp]=radon(P(:,:,1),theta);
[sz1R sz2R]=size(R)
for k=1:sz2R
plot(R(:,k));grid on
disp(k);
str=input(' ','s');
end
theta =86.5º
theta = 90º theta close to 90º
Now we have clear sharp peaks, apparently one per tree line, or so we think:
theta=90º is k=50. Clipping
L1=R(:,50)'
L1(L1>2000)=2000 % clipping
[pks,locs]=findpeaks(L1) % has 78 peaks. There are about 15 pixels between adjacent pulse ramps.
plot(L1,'LineWidth',1.5,'Color','b');grid on
hold all;plot(locs,pks,'x','Color','r')
ax=gca
ax.YLim=[0 2100]
numel(pks)
=
78.00
now remove false peaks
pks2=pks(pks>1600)
numel(pks2)
=
48.00
It looks good, but what if the chosen sample line suffers a few more gaps like this one?
On the other hand, DTFT on a sample line with a few missing peaks is more robust, and will still show a frequency peak on the right frequency. This is why DTFT, in this particular case, is a better choice.
10.- So, DTFT line 20 of the ribbon, no need to have horizontal tree lines.
L1=D2(:,20)
length(L1)
=
871
r_L1=xcorr(L1);r_L12=r_L1([floor(length(r_L1)/2):end])
x=r_L12';n=[1:1:length(r_L12)] % just reusing code
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);plot(magX)
since we are after the cycle of a simple tone, expected around 100 on the previous graph, then let's find the exact position within such narrower frequency range
[pks_magX loc_pks_magX]=findpeaks(magX([90:110]),[90:110])
loc_pks_magX(find(pks_magX==max(pks_magX)))
pks_magX =
1.0e+08 *
Columns 1 through 4
0.031378416088323 0.021018552598383 2.416607712099700 0.584907695507022
Column 5
0.598161849635904
loc_pks_magX =
92 94 102 106 108
the first peak above DC
f_peak=loc_pks_magX(find(pks_magX==max(pks_magX))) % peak amplitude is 57
=
102
with this DTFT Fs=
numel(magX)
=
501
the maximum frequency that could occur along a sample line is floor(numel(L1)/2). So the answer is
f_peak/numel(magX)*numel(L1)/2 % 57/501*871/2
ans =
49.547904191616766
which is quite Carling, isn't it?
Following, different attempts to improve accuracy.
On the following graph, the 1st peak above DC, the one with the marker on in previous DTFT graph is not sharp enough, moving the marker around the peak you notice there are at least 2 samples with quite similar value.
11.- Attempting to improve result (see 50.0) just increasing
frequency resolution is not going to work:
k=0:2500;w=(pi/2500)*k
X=x*(exp(-j*pi/2500)).^(n'*k)
X=x*(exp(-j*pi/2500)).^(n'*k)
magX=abs(X)
plot(magX); grid on
[pks_magX loc_pks_magX]=findpeaks(magX([260:310]),[260:310])
f_peak=loc_pks_magX(find(pks_magX==max(pks_magX)))
=
284 % peak value 4.026e+08
f_peak/numel(magX)*numel(L1)/2
=
49.453018792483007
Again, because the part I single line had 1665 samples, the spectrum is more accurate than previously, yet is there really need for such added accuracy?
12.-more wood, further interpolating
resamp = makeresampler({'nearest','cubic'},'fill');
stretch = maketform('affine',[5.1 0; 0 5.1; 0 0]);
D3 = imtransform(D(:,:,1),stretch,resamp);
imshow(D3)
imstransform stops updating the pixel coordinates, displaying a decimated version, and not even updating the size of the picture. But the warning message has a way around:
Warning: IMTRANSFORM has increased the scale of your output pixels in X-Y space to keep
the dimensions of the output image from getting too large. This automatic scale change
will be removed in a future release.
To ensure that output pixel scale matches input pixel scale specify the 'XYScale'
parameter. For example, call IMTRANSFORM like this:
B = IMTRANSFORM(A,T,'XYScale',1)
> In imtransform>compute_xyscale (line 323)
In imtransform>parse_inputs (line 463)
In imtransform (line 265)
Warning: Image is too big to fit on screen; displaying at 67%
> In images.internal.initSize (line 71)
In imshow (line 309)
add: 'XYScale',1 increasing the value of the XYScale actually reduces the size of the output image that is supposed to be larger, not smaller
resamp = makeresampler({'nearest','cubic'},'fill');
stretch = maketform('affine',[5.1 0; 0 5.1; 0 0]);
D3 = imtransform(D(:,:,1),stretch,resamp,'XYScale',1);
L1=D3(:,100)
plot(L1);grid on
r_L1=xcorr(L1)
r_L12=r_L1([floor(length(r_L1)/2):end])
x=r_L12';n=[1:1:length(r_L12)] % just reusing code
r_L1=xcorr(L1)
r_L12=r_L1([floor(length(r_L1)/2):end])
k=0:2500;w=(pi/2500)*k
X=x*(exp(-j*pi/2500)).^(n'*k)
X=x*(exp(-j*pi/2500)).^(n'*k)
magX=abs(X)
plot(magX); grid on
[pks_magX loc_pks_magX]=findpeaks(magX([95:105]),[95:105])
f_peak=loc_pks_magX(find(pks_magX==max(pks_magX)))
f_peak/numel(magX)*numel(L1)/2
=
101
f_peak/numel(magX)*numel(L1)/2
=
49.773090763694526
k=0:2500;w=(pi/2500)*k
X=x*(exp(-j*pi/2500)).^(n'*k)
X=x*(exp(-j*pi/2500)).^(n'*k)
magX=abs(X)
plot(magX); grid on
[pks_magX loc_pks_magX]=findpeaks(magX([95:105]),[95:105])
f_peak=loc_pks_magX(find(pks_magX==max(pks_magX)))
f_peak/numel(magX)*numel(L1)/2
s=15500
[590:650]
=
49.296819560028389
13.- FFT, sharper peak
magX=abs(fft(r_L12))
plot(magX); grid on
50/(2466/2)*2465/2
=
49.979724249797236
14.- PLOMB
[Pxx,F]=plomb(x,n)
197/4932*2465/2
=
49.230028386050279
15.- tried different DTFT values dtft_range=[500:500:2500]
f_peaks_a=zeros(numel(dtft_range,sz1D2))
for s=[1:1:numel(dtft_range)]
for k=1:sz1D2
L1=D2(:,k)
r_L1=xcorr(L1)
r_L12=r_L1([floor(length(r_L1)/2):end])
k=0:(s*500);w=(pi/(s*500))*k
X=x*(exp(-j*pi/(s*500))).^(n'*k)
X=x*(exp(-j*pi/(s*500))).^(n'*k)
magX=abs(X)
plot(magX); grid on
magX=abs(X);
[pks_magX loc_pks_magX]=findpeaks(magX([90:110]))
end
end
So it is reasonable to answer that on vertical sample line 20 there are 49 tree lines.
Certainty statements are beyond the scope of this QA. But, at least in this case, the spectral basic analysis of the autocorrelation of a single line hatt has many missing peaks (tree lines) is more robust than radon.
So this is it: 49.5 tree lines on a single count line is as good as it gets, and if you think it twice, it may be that there is half cycle lost on either head or tail of the vertical sample line.
Other spectrum analysis tools:
16.- spectrogram, can we spot the frequency variation?
spectrogram(r_L12,'yaxis')
17.- periodogram directly on sample line L1
periodogram(double(L1'))
18.- periodogram on r_L1
periodogram(R_L12)
  2 个评论
Image Analyst
Image Analyst 2016-5-16
Wow! He should be paying you for the amount of work you're putting in on his project.
John BG
John BG 2016-5-18
Mr Image Analyst
I see it the other way round, we should be paying you for reading our questions and answers. Having had a quick look to some of your answers it's a privilege for us to receive advice from you.
And my test, it's just a hobby, gym. Now that R2016 comes with the Image Processing toolbox in the base pack, there are some tasks that I want to get used to solve with the IP toolbox.
I should have a readable more compact answer including the filtering by this Friday.
And if you don't mind me asking, but I am missing using del2 and perhaps other contour detection functions from the beginning.

请先登录，再进行评论。

John BG
John BG 2016-5-14
编辑：John BG 2016-5-19
Shariful
PART I
1.- acquire image
A=imread('lines_count00.jpg');
[im_1 im_2 im_3]=size(A);
% imshow(A)
2.- you may want to consider processing the variance of the RGB components
surf(var(double(A),0,3))
but there is simpler way
3.- sample a single line, choose 2 points make 1st point as close to [0 max(Y)] as possible, inside the image
p3=ginput(2); p3=uint64(floor(p3))
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);p=plot(v3_ref,im_line3_1,'r',v3_ref,im_line3_2,'g',v3_ref,im_line3_3,'b');grid on
p(1).LineWidth=1.5;p(2).LineWidth=1.5;p(3).LineWidth=1.5;
4.- inverting
im_line3_1=255-im_line3_1;
im_line3_2=255-im_line3_2;
im_line3_3=255-im_line3_3;
figure(3);p2=plot(v3_ref,im_line3_1,'r',v3_ref,im_line3_2,'g',v3_ref,im_line3_3,'b');grid on
p2(1).LineWidth=2;p2(2).LineWidth=2;p2(3).LineWidth=2;
5.- amplify
figure(3);p3=plot(v3_ref,20*log10(double(im_line3_1)),'r',v3_ref,20*log10(double(im_line3_2)),'g',v3_ref,20*log10(double(im_line3_3)),'b');grid on
p3(1).LineWidth=1.5;p3(2).LineWidth=1.5;p3(3).LineWidth=1.5;
amplifying doesn't solve the lack of resolution: there are about 52 'horizontal' lines, rought visual count, from the bottom left corner of the original picture moving up perpendicular to the 'tree lines'. The reference vector, the vector used to plot such sample line is 145 elements long only.
We need more samples, I would be more comfortable with at least 5 samples per cycle, now only 3 roughly, too close to the bare minimum of 2 samples per cycle of highest frequency components, that failing to hold sampling above such threshold, aliasing occurs.
%%%
Comment: You don't mind me pointing out, do you, that you have neither corrected perspective distortion nor rotation.
When processing airborne acquired imagery it's common practice to correct these 2 points before doing any measurement that may change had the camera and platform had a different attitude when the picture taken.
%%%
6.- switch to black & white
% A=imread('lines_count00.jpg');
% [im_1 im_2 im_3]=size(A);
% figure(1);imshow(A)
A2=(A(:,:,1)+A(:,:,2)+A(:,:,3));
figure(3);imshow(A2);
[sz1_A2 sz2_A2 sz3_A2]=size(A2)
note that the size of A is 484x701x3 and the size of A2 is 484x701, it would be great to have more samples.
% DC block stretching to fill up available range [0 255]
%
% minA=min(min(min(A)));
% maxA=max(max(max(A)));
% A2=A2-minA;
% A2=A2.*maxA;
% A2=uint8(A2)
% figure(2);imshow(A2)
7.- increase amount of samples with cubic interpolation, the matrix in maketform controls the stretching 1.8 times on both X and Y directions
resamp = makeresampler({'nearest','cubic'},'fill');
stretch = maketform('affine',[1.8 0; 0 1.8; 0 0]);
B3 = imtransform(A2,stretch,resamp);
figure(4);imshow(B3)
now we have enough samples
[sz1_B3 sz2_B3]=size(B3)
sz1_B3 =
871
sz2_B3 =
1261
8.- to take advantage of script I already have to solve similar problem where lines are vertical, use the marker and read the pixel position on the right hand side of the sample line to know the exact angle needed turn vertical the lines of interest.
Rotating image 0.350942428742509 rad angle=20.10752 approx. 20º06'
B4=imrotate(B3,-(90-20.10752));
figure(5);imshow(B4);
but
not perfect, can't be, can't tell whether the tree lines are really parallel or it's just the uncorrected perspective, but let's carry on.
Apparently the rotation has not reduced the amount of samples
[sz1_B4 sz2_B4]=size(B4)
sz1_B4 =
1484
sz2_B4 =
1252
9.- since the rotation angle is not multiple of 90º or 180º, there is, a loss of samples caused by the rotation function that compresses the rotated image within the frame of the initial interpolated image.
Roughly, we are back to the initial 3 samples per cycle.
So let's interpolate again:
resamp = makeresampler({'nearest','cubic'},'fill');
stretch = maketform('affine',[1.8 0; 0 1.8; 0 0]);
B5 = imtransform(B4,stretch,resamp);
figure(6);imshow(B5)
[sz1_B5 sz2_B5]=size(B5)
sz1_B5 =
2671
sz2_B5 =
2253
10.- sample 1 line
p3=ginput(2); p3=uint64(floor(p3))
py_ref=min(p3(2,2),p3(1,2))
im_line_1=B5(py_ref,[min(p3(1,1),p3(2,1)):max(p3(1,1),p3(2,1))])
v3_ref=[1:1:(abs(double(p3(1,1))-double(p3(2,1))))+1]
figure(2);plot(v3_ref,im_line_1,'r','LineWidth',1.5);grid on
11.- invert
im_line_1=255-im_line_1;
figure(3);plot(v3_ref,im_line_1,'r','LineWidth',1.5);grid on
numel(findpeaks(im_line_1,'minpeakheight',25))
ans =
61
peaks_prob_1=0
for k=10:1:250
peaks_prob_1=[peaks_prob_1 numel(findpeaks(double(im_line_1),'minpeakheight',k))];
end
figure(7);plot(peaks_prob_1);grid on
v_prob_1=[1:length(peaks_prob_1)];
figure(8);hist(peaks_prob_1,v_prob_1);
both cummulative and histogram
both show the 2 highest reasonable counts 46 and 48, yet visual inspection shows the right amount is 52 or 53
12.-
r_L1=xcorr(im_line_1);figure(3);
plot(r_L1([floor(length(r_L1)/2):end]));grid on
from the previous plot, the average distance in pixels is 32, let's sweep around this expected value with findpeaks and option 'minpeakdistance':
r_L12=r_L1([floor(length(r_L1)/2):end])
prob2=[]
for k=20:1:40
[pks_xcorr1puls,locs_xcorr1puls]=findpeaks(r_L12,'minpeakdistance',k)
prob2=[prob2 numel(pks_xcorr1puls) ]
end
figure(9);plot(prob2);grid on
jiggling the distance between adjacent peaks of the auto correlation we observe that tops there are 55 pixels between 2 adjacent tree lines.
The sample line was taken from B5, not B4, in red
there at least 1665/55~30 tree lines
13.-
%nL1=[1:length(r_L1)] % just copied lines from another script with DTFT
x=r_L1;n=[1:1:length(r_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)')
1665/2 is the amount of pixels per cycle on Fs
the peak that the DTFT detects is 0.064 normalized frequency
1665/2*0.064=53.28 in the sampled line there are 53 cycles of what can only be section across tree lines.
If you only want a sample measurement, perhaps the previous steps are ok to you.
If you want to refine, then you can
1.- band pass filter around the expected value range of let's say [43 63].
2.- process a bunch of lines, average the results
3. make the measurement independent of the points measured by stating that the photo has tree lines density of #lines/cm or #lines/inch of photograph.
I hope you don't mind me repeating that without perspective correction the tree line measurement on 1 line, or averaged along a group of lines, may not be used on any other line of the photo because at this point we don't know whether the tree lines are parallel or not.
If you correct the perspective distortion, then is when you can tell if the the tree lines are parallel (or not), and if parallel then you can use 1 or a few measurements all over the photo, and then approximate trees counting if you want.
If you want I can help you these 3 points closing points (filter, average multiple lines, obtain lines density assuming tree rows are parallel) if the previous measurement on a sample lines is not enough.
I can also help you to draft a simple perspective correction, that needs
  1. altimeter reading at the moment of taking the photo
  2. altitude difference between altimeter and camera (you need to know the aircraft)
  3. camera down angle,
  4. (and camera CCD aspect ratio may also help)
If you find this answer of any help solving your question,
please click on the thumbs-up vote link,
thanks in advance
John
jgb2012@sky.com
  3 个评论
显示 1更早的评论
Manjiree Waikar
Manjiree Waikar 2017-11-3
Sir which technique is better to use over here, Radon Transform, Wavelet Transform or Hough Transform to extract the line parameters?
Image Analyst
Image Analyst 2017-11-3
Over where? You have given us no information about you or your application.

请先登录，再进行评论。

Image Analyst
Image Analyst 2016-5-14
It's not clear to me, even if I manually counted them, much less tried to do it automatically. Plus it's somewhat arbitrary. You just shift the photo a little bit and suddenly you have a dramatically different number of rows. And the number of rows changes a lot if you just rotate the image.
And what makes a row? 100 trees? 10 trees? How about 2 trees? Are the isolated trees in the blue field part of the same row as other trees farther away but on the same extension of the line? Or are they their own separate, different row? What about where the horizontal rows interweave with the vertical rows? Do you count those trees twice because they're in rows of different angles?
No, I think counting the number of rows is not a useful metric. I think something like area fraction would be a more useful, stable, and robust measurement.
  5 个评论
显示 3更早的评论
Shariful Islam
Shariful Islam 2016-5-14
编辑：Shariful Islam 2016-5-14
Thanks, I will measure some indices (e.g NDVI), Leaf chlorophyll content, Leaf area index, tree count etc- to estimate plant vigor, health etc. But it is an simple measure to count rows and i think if i count/identify the rows i can also measure the length of each rows. Again in some image (not this) there are rows that fully degraded.
I am stumbling on devicing a generalized procedure for row count.
Image Analyst
Image Analyst 2016-5-14
I have not done this. My only suggestions remain the ones I already gave. Try radon(). Try hough() or houghlines(). And find a paper in the link I gave you that does the same thing. I don't have any further suggestions because I think one of those should work.

请先登录，再进行评论。

ATHINISWARYRAM RAMADAS
how to count the lines
  1 个评论
Image Analyst
Image Analyst 2018-1-19
This is not an answer to Shariful's question.
Start your own question with your own image after you read this link.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Feature Detection and Extraction 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by