How to get line count from an image?

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)-

Answer 1

[John BG]

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

[Image Analyst]

You've done an extensive amount of work over the past two weeks - perhaps even setting an Answers record. If I were you though, I'd wait to see if Shariful responds before doing anymore. It seems he has totally forgotten about this discussion, not having responded in almost two weeks in spite of repeated, very generous updates from you. It's your choice though.

[Shariful Islam]

@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]

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

Answer 2

[John BG]

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)

[Image Analyst]

Wow! He should be paying you for the amount of work you're putting in on his project.

[John BG]

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.

Answer 3

[John BG]

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

[e-mail]

[Image Analyst]

Are you going to try the radon transform?

[Manjiree Waikar]

Sir which technique is better to use over here, Radon Transform, Wavelet Transform or Hough Transform to extract the line parameters?

[Image Analyst]

Over where? You have given us no information about you or your application.

Answer 4

[Image Analyst]

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.

[Shariful Islam]

Thanks, Ah! rows are not arbitrary- it is an ORCHARD. 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.

[Image Analyst]

I know it's an orchard/plantation. I see I did not get my point across. So....

Convert the image to YCbCr color space. Take the Cb channel and threshold it at around 100. Then try radon() or hough() on it. See the hough example for gantry crane and see how it found lines.

[Image Analyst]

(My last attempt at persuading you.) Don't you realize that counting lines will not tell you how many trees are in the image or how healthy the orchard is or how much fruit it will produce? If you count lines, you'd get the same answer with your posted image, or with an image where the blue field was gone and replaced by healthy tree lines.

You say "i want to know how many rows of trees degraded after some time." So, if there's a drought or insect infestation, perhaps every single tree line is degraded, yet the number of lines remains the same. So measuring tree lines is not a good measure unless the whole entire line is gone. I'm not sure why you don't believe that area fraction and spectral signature is a better approach of tree health and counting.

To see how others do it, see section 22.5.9 here: http://www.visionbib.com/bibliography/contentscartog.html#Cartography,%20Aerial%20Images,%20Remote%20Sensing,%20Buildings,%20Roads,%20Terrain,%20ATR

[Shariful Islam]

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]

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.

Answer 5

[ATHINISWARYRAM RAMADAS]

how to count the lines

[Image Analyst]

This is not an answer to Shariful's question.

Start your own question with your own image after you read this link.