Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
C
CanopyPorosity
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Model registry
Operate
Environments
Monitor
Incidents
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
GitLab community forum
Contribute to GitLab
Provide feedback
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
asalgado
CanopyPorosity
Commits
6d3d9619
Commit
6d3d9619
authored
5 years ago
by
asalgado
Browse files
Options
Downloads
Patches
Plain Diff
main script including processes for canopy porosity calculation based on gap distribution
parent
fde48fa5
No related branches found
No related tags found
No related merge requests found
Pipeline
#838
canceled
5 years ago
Stage: build
Stage: test
Changes
1
Pipelines
1
Show whitespace changes
Inline
Side-by-side
Showing
1 changed file
canopyPorosityBatch.m
+500
-0
500 additions, 0 deletions
canopyPorosityBatch.m
with
500 additions
and
0 deletions
canopyPorosityBatch.m
0 → 100644
+
500
−
0
View file @
6d3d9619
function
[
results_cp
,
results_fF
,
results_fC
,
results_clump
,
results_lai
]
=
canopyPorosityBatch
(
im
,
rws1
,
cols1
,
thresh1
,
filePathIn
,
imNameIn
)
% % Painting the tree canopy area
% ASA Salgadoe, 21/06/2017
% the image thresholding setting: Suitable RGB shade side Only
% With the help of edge detected, Inverse image and graythresh value
% the canopy area is set as ONEs (3),canopy (1) and Sky Zeros (0)
% clear the camand window
% clear the variables
% target is to count the edges as poros
% two methods involved: 1) edges in a given canopy area 2)edges in a set
% polygon area
% porosity = (sum of edges/sum of canopy area)*100
% clc;
% clear all;
% (1)load the RGB image
% im=input('Type name of the RGB image : ');
% read the image to the im variable
img
=
imread
(
im
);
% For High resolution images*****************
% img=imresize(img,.2);
% *************************************
% (2)1st Cropping
% Crop image from boundaries
% imCrop=cropBorder(img);
% % imCrop=imcrop(img);
rows
=
length
(
img
);
cols
=
length
(
img
(
1
,:,
1
));
minCol
=
cols
*
0.1
;
width
=
cols
*
0.8
;
minRow
=
rows
*
0.1
;
height
=
rows
*
0.8
;
imCrop
=
imcrop
(
img
,[
minCol
minRow
width
height
]);
% Manuall crop;;;;;
% imCrop=img;
% (3)Create the inverse image of the RGB*******************
if
size
(
imCrop
,
3
)
==
3
inv
=
(
255
-
(
imCrop
(:,:,
1
)
+
imCrop
(:,:,
2
)
+
imCrop
(:,:,
3
))/
3
);
else
disp
(
"Sorry Not a RGB image, please load a RGB !!!....."
);
end
% (4)Excute edge detection algorythum, 'Canny'********************
% imEdge=edge(inv,'Canny',[.3 0.65]+-0.13); for mobile
% imEdge=edge(inv,'Canny',[.3 0.7]+-0.35); for analystical
imEdge
=
edge
(
inv
,
'Canny'
,[
.
3
0.65
]
+-
0.13
);
% (5)Get threshold value
% excute Otsu's method of histogram thresholding to get value to segment image
% of equalized histogram
% thresh=graythresh(histeq(inv))*255;
thresh
=
graythresh
(
inv
)
*
255
;
% (4)Image segmentation
% pixels higher thatn threshold get 1 and less 0s
% can decide propotion of the threshold amout
set
=
255
-
thresh
;
% invT=histeq(inv)>(thresh+(set*0));
invT
=
inv
>
thresh
;
% (5)overlay two indexed images
% finally overlay the edge and invC
% rgb2gray results canopy 255 colour code and edges 226 color code
invC
=
imoverlay
(
invT
,
imEdge
);
invGray
=
rgb2gray
(
invC
);
% (6) Re-coding overlayed image
% re-coding edge=3 , canopy =2 and sky=0
for
col
=
1
:
length
(
invGray
(
1
,:,
1
))
% iterating within rows
for
row
=
1
:
length
(
invGray
(:,
1
,
1
))
if
invGray
(
row
,
col
)
==
255
% the canopy area and sky included in threshould get 2 similar to invC recode by 2
invGray
(
row
,
col
)
=
2
;
elseif
invGray
(
row
,
col
)
==
226
% the border edge codes as 3 for future detection
invGray
(
row
,
col
)
=
3
;
end
% complete noise area 0
end
end
% (8)Painting the Only canopy area
% Edge will be same 3, canopy will be 1, and sky will be 0
% iterating through top down in cols
% the graythreshold image will give sky as 1, so have to use edge filtered
% image to itereate across image
for
col
=
1
:
length
(
invGray
(
1
,:,
1
))
% iterating within rows
for
row
=
1
:
length
(
invGray
(:,
1
,
1
))
% Always the border edge(3) is detected when coming from top to down in col.
% or a already recoded (1) area
if
invGray
(
row
,
col
)
==
1
||
invGray
(
row
,
col
)
==
3
if
row
==
length
(
invGray
(:,
1
,
1
))
break
% if 1 its the recoded area, if 2 it is the canopy area
elseif
invGray
(
row
+
1
,
col
)
==
1
||
invGray
(
row
+
1
,
col
)
==
2
invGray
(
row
+
1
,
col
)
=
1
;
elseif
invGray
(
row
+
1
,
col
)
==
0
invGray
(
row
+
1
,
col
)
=
0
;
end
% allowing to iterate manny rows
% checking the next available number or rows before iterate
% with a detected 1, goes up to next seven rows
% looping the following number of rows ahead everytime
% set up to 50 rows ahead always from the point of processing
% after that prossesing point will chose point one step down and loop the
% next following 50 rows
for
i
=
1
:
50
if
row
==
length
(
invGray
(:,
1
,
1
))
-
i
break
elseif
invGray
(
row
+
i
,
col
)
==
1
||
invGray
(
row
+
i
,
col
)
==
2
invGray
(
row
+
i
,
col
)
=
1
;
end
end
else
invGray
(
row
,
col
)
=
0
;
end
end
end
% ----------------------------------------------------
% Edge (3), folia (1), and sky (0)
% ------------------------------------------------------
%
map
=
[
0
,
0
,
0
;
0
,
1
,
0
;
0
,
0
,
1
];
% figure(),
% subplot(1,6,1);imshow(im);title('Original');
% subplot(1,6,2);imshow(invT);title('Colour Segmented');
% subplot(1,6,3);imshow(imEdge);title('Edge Detected');
% subplot(1,6,4);imshow(invGray,map);title('Canopy Painted');
% subplot(1,4,4);imshow(invGray,map);title('Canopy Painted');
% impixelinfo;
% (10) Get Max Point of Canopy
% After edge detection get the max point from the image
% the X and Y values are returned and stored in imX and imY
[
imX
,
imY
]
=
marching_through_image
(
invGray
);
% (9)making all the Region of interest set 1
for
col
=
1
:
length
(
invGray
(
1
,:,
1
))
% iterating within rows
for
row
=
1
:
length
(
invGray
(:,
1
,
1
))
if
invGray
(
row
,
col
)
==
3
invGray
(
row
,
col
)
=
0
;
end
end
end
% ----------------------------------------------------------------------------------------
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%
% height=(length(invGray)-imY);
% % height=200;
% imCrop_final1=reCrop(invGray,0,imY,length(invGray(1,:,1)),height);
% %
%
%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% figure();
% subplot(1,5,1),imshow(im), title("Original");
% subplot(1,5,2),imshow(invT), title("Inverse");
% subplot(1,5,3),imshow(imEdge),title("Edge Detected");
% % subplot(1,5,4),imshow(invGray,map),title("Painted");
% % subplot(1,5,5),imshow(imCrop_final1,map), title("M2: Porosity = "+porosity1+" %");
%
% subplot(1,5,4),imshow(imCrop_final1,map)
% >>>>>>>>>>> BIG GAP segmentation <<<<<<<<<<<<<<<<<<<<<<<<<<
% % >>>>>>>>>>>>>>>>>>>>>>>>>>> subdivision of image
% % Splitting the image by mat2col function
% user determining the devisor of col and row
imCrop_final1
=
invGray
;
% ***************************************
% imwrite(invGray,map,strcat(filePathIn,imNameIn,'.tiff'));
% *****************************************
% set a user defined cols and rows to be subdivided of image
rws
=
rws1
;
cols
=
cols1
;
setThresh
=
thresh1
;
% Best performances
% rws=50;
% cols=40;
numRW
=
0
;
numCL
=
0
;
% get the integer part after devision
numRW
=
fix
(
length
(
imCrop_final1
(:,
1
))/
rws
);
% get the integer part after devision
numCL
=
fix
(
length
(
imCrop_final1
(
1
,:))/
cols
);
% >>>>>>>>>>>>>>>>>>>>>> getting cropping dimensions
cropHeight
=
0
;
cropWidth
=
0
;
setCol
=
1
;
setRow
=
1
;
% Resulting number of rows and cols
resultRW
=
numRW
*
rws
;
resultCL
=
numCL
*
cols
;
% for rows can directly crop from bottom
setRow
=
1
;
cropHeight
=
resultRW
-
1
;
% for cols try to crop evenly from both sides if even remainder if odd from one side
% check if the remainging rows is a even or not even number
% Fix(x/2)->Left crop, (x-left)=right crop ; x=total cols to remove
% decide the cols removing from left edge then from there take desired
% width as total width
totalColRemove
=
(
length
(
imCrop_final1
(
1
,:,
1
))
-
resultCL
);
setCol
=
fix
(
totalColRemove
/
2
);
if
setCol
==
0
setCol
=
1
;
end
cropWidth
=
resultCL
-
1
;
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Execute cropping
imCrop_final2
=
imcrop
(
imCrop_final1
,[
setCol
setRow
cropWidth
cropHeight
]);
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Generate array for subdivision
arrayRW
=
0
;
arrayCL
=
0
;
for
a
=
1
:
rws
arrayRW
(
1
,
a
)
=
numRW
;
end
for
b
=
1
:
cols
arrayCL
(
1
,
b
)
=
numCL
;
end
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> execute image subdivision into cells
cells
=
mat2cell
(
imCrop_final2
,
arrayRW
,
arrayCL
);
% Edge (3), folia (1), and sky (0)
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>> iterate through the cells storing values
% get a duplicate of cells
cellsD
=
cells
;
skycount
=
0
;
% 0 and 3 included
foliacount
=
0
;
% only 1 included
ratio
=
0
;
for
arrayCol
=
1
:
length
(
cells
(
1
,:))
for
arrayRw
=
1
:
length
(
cells
)
% subCell=cells{arrayRw,arrayCol};
skycount
=
0
;
% 0 and 3 included
foliacount
=
0
;
% only 1 included
% iteratin in single cell
for
subCol
=
1
:
length
(
cells
{
arrayRw
,
arrayCol
}(
1
,:))
for
subRow
=
1
:
1
:
length
(
cellsD
{
arrayRw
,
arrayCol
}(:,
1
))
if
cells
{
arrayRw
,
arrayCol
}(
subRow
,
subCol
)
==
0
||
cells
{
arrayRw
,
arrayCol
}(
subRow
,
subCol
)
==
3
skycount
=
skycount
+
1
;
end
if
cells
{
arrayRw
,
arrayCol
}(
subRow
,
subCol
)
==
1
foliacount
=
foliacount
+
1
;
end
end
end
% Edge (3), folia (1), and sky (0)
if
skycount
==
0
skycount
=
1
;
end
if
foliacount
==
0
;
foliacount
=
1
;
end
sky
(
arrayRw
,
arrayCol
)
=
skycount
;
folia
(
arrayRw
,
arrayCol
)
=
foliacount
;
ratio
(
arrayRw
,
arrayCol
)
=
skycount
/(
foliacount
+
skycount
);
end
end
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Re-code GAP area
% get the mean value of total ratio.
% set a threshould for the ratio
% if the ratio is higher than that thresh in a cell set the sky pixels zero
% if the ratio is lower than that thresh in a cell count that sky set as (2)
% GAP is a sky are
% set Edge (3), and sky (0) to small gaps(2)
% Edge and sky included in smallGap is (2),
% Folia (1)
% Edge and sky included big gaps(0)
meanRatio
=
mean
(
mean
(
ratio
));
% setThresh=20;
for
ratioCol
=
1
:
length
(
ratio
(
1
,:))
for
ratioRw
=
1
:
length
(
ratio
(:,
1
))
% best ratios 1.5
if
ratio
(
ratioRw
,
ratioCol
)
<
setThresh
&&
ratio
(
ratioRw
,
ratioCol
)
~=
0.000000
% sometimes the ratio becomes zero for 10x10 subdivisions
% take all the sky and edge pixels and recode if less than thresh (2) if
% higher than thresh (0)
for
subCol
=
1
:
length
(
cellsD
{
ratioRw
,
ratioCol
}(
1
,:))
for
subRow
=
1
:
length
(
cellsD
{
ratioRw
,
ratioCol
}(:,
1
))
if
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
==
0
||
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
==
3
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
=
2
;
end
end
end
elseif
ratio
(
ratioRw
,
ratioCol
)
>
setThresh
&&
ratio
(
ratioRw
,
ratioCol
)
~=
0.000000
for
subCol
=
1
:
length
(
cellsD
{
ratioRw
,
ratioCol
}(
1
,:))
for
subRow
=
1
:
length
(
cellsD
{
ratioRw
,
ratioCol
}(:,
1
))
if
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
==
0
||
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
==
3
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
=
0
;
end
end
end
end
end
end
% Edge(3) and sky(0) included in smallGap is (2),
% Folia (1)
% Edge and sky included big gaps(0)
% >>>>>>>>>>>>>>>>>>>>>>> Counting sky area in the BIG GAPS
skybigGap
=
0
;
skyPxl
=
0
;
for
ratioCol
=
1
:
length
(
ratio
(
1
,:))
for
ratioRw
=
1
:
length
(
ratio
(:,
1
))
skyPxl
=
0
;
if
ratio
(
ratioRw
,
ratioCol
)
>
setThresh
&&
ratio
(
ratioRw
,
ratioCol
)
~=
0.000000
for
subCol
=
1
:
length
(
cellsD
{
ratioRw
,
ratioCol
}(
1
,:))
for
subRow
=
1
:
length
(
cellsD
{
ratioRw
,
ratioCol
}(:,
1
))
% Edge(3) and sky(0)
if
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
==
0
||
cellsD
{
ratioRw
,
ratioCol
}(
subRow
,
subCol
)
==
3
skyPxl
=
skyPxl
+
1
;
end
end
end
end
skybigGap
(
ratioRw
,
ratioCol
)
=
skyPxl
;
end
end
% >>>>>>>>>> Calculations <<<<<<<<<<<<<<<<<<<<<<
tGap
=
sum
(
sum
(
sky
));
tPixl
=
numel
(
imCrop_final2
);
tLGap
=
sum
(
sum
(
skybigGap
));
% Fraction of foliage cover fF
% Fraction of Crown cover fC
% crown porosity cP
% Clumping index clumpIndex
fF
=
0
;
fC
=
0
;
cP
=
0
;
clumpIndex
=
0
;
fF
=
1
-
(
tGap
/
tPixl
);
fC
=
1
-
(
tLGap
/
tPixl
);
% cP=1-(fF/fC);
% *******************************
tFolia
=
sum
(
sum
(
folia
));
tSGap
=
tGap
-
tLGap
;
tCover
=
tSGap
+
tFolia
;
if
tSGap
==
0
cP
=
(
tFolia
/
tCover
)
*
100
;
else
cP
=
(
tSGap
/
tCover
)
*
100
;
end
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% clumping index
% clumpInd
% paper 13
clumpIndex
=
(
1
-
cP
)
*
log
(
1
-
fF
)/
log
(
cP
)
*
fF
;
% LAI (Matlab) paper 13
% extinction-coefficient k (0.5) for euclypt tree
k
=
0.5
;
LAI
=-
fC
*
log
(
cP
)
*
1
/
k
;
% Effective LAI
LAIe
=
LAI
*
clumpIndex
;
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Recreate Full Image
imCreate
=
cell2mat
(
cellsD
);
% figure,imshow(imCreate,map);
% >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% extract= extractAfter(im,'IMG_');
% finalName=strcat(num2str(thresh1),'_',num2str(rws1),'_',extract);
% % imwrite(invT,strcat('C:\Users\asalgado\UNE_Cloud_sync\PhD work\Pilot Test 2_Bundy\Analysis\Matlab\Porosity Work\Batch_images\MobileCam\RGBsmall_selected\Batch1,2&3\',thresh+char(extract),'inv','.tiff'));
% *****************************
% imwrite(imCreate,map,strcat(filePathIn,imNameIn,'_Gap','.tiff'));
% ***********************************
results_cp
=
cP
;
results_fF
=
fF
;
results_fC
=
fC
;
results_clump
=
clumpIndex
;
results_lai
=
LAIe
;
end
% results.clumpIn =clumpIndex;
% results.LAIe =LAIe;
%--------------------------------------------------------------------------
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment