First commit

LICENSE.md

@@ -25,4 +25,4 @@ DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

arcade.m

@@ -0,0 +1,141 @@
+%This is the main MATLAB file used to call other functions for the purpose
+%of crater detection, segmentation, and computation of features: centroids, areas, and
+%perimeters of detected craters.
+%Crater shells are detected using a computer vision cascade object detector.
+%Trained features are stored in detector object that is created using the craterDetector.xml file 
+%Detected images are segmented using the Chan-Vese active contours without edges approach
+
+%ARtillery Crater Analysis and Detection Engine (ARCADE)
+%Developed by Ali M Bukar for Rudiment.info
+%Centre for Visual Computing
+%University of Bradford, UK
+
+clc
+
+clear all
+
+
+    %Read xml file and create detector object
+
+    detector3 = vision.CascadeObjectDetector('craterDetector.xml');
+
+    %read the textfile containing image url, lat1 lon1 and lat2 lon2
+    filename = 'ReadCraters.txt';
+    fileID = fopen(filename);
+    c = textscan(fileID, '%s %f %f %f %f', 'Delimiter', '\t');
+    fclose(fileID);
+
+    datalength = length(c{1});
+
+    filename2 = 'Settings.txt';
+    fileID2 = fopen(filename2);
+    c2 = textscan(fileID2, '%d');
+    fclose(fileID2);
+
+    itr = c2{1};
+
+    %Iteratively perform detection, segmentation, and feature extraction on
+    %each data (row) retrieved from ReadCraters.txt
+
+    for ii = 1:datalength
+
+        imname = c{1}{ii};
+        lata = c{2}(ii);
+        lona = c{3}(ii);
+        latb = c{4}(ii);
+        lonb = c{5}(ii);
+
+        numst = num2str(ii);
+
+            %call detect and segment function
+            [number_craters, detectedIm, segmentedIm, img] = detect_n_segment_craters(imname, detector3, itr);
+
+            %format output name
+            [pathstr, name, ext] = fileparts(imname);
+            outputname = strcat(name,'_',numst,'.txt');
+            outputname2 = strcat(name,'_',numst,'.xlsx');
+            if(number_craters == 0)
+
+                    %myformat = '%f\t%f\t%f\t%f\r\n';
+                    T = table(0,0,0,0);
+                    write(T, outputname, 'Delimiter', '\t');
+                    write(T, outputname2);
+
+            else
+
+                    [pathstr, name, ext] = fileparts(imname);
+
+                    detectedstr = strcat('detected','_', name, ext);
+                    segmentedstr = strcat('segmented','_', name, ext);
+
+                    %save image showing detected craters annotated
+                    %also save the segmented image to working directory
+
+                    imwrite(detectedIm, detectedstr);
+                    imwrite(segmentedIm, segmentedstr);
+
+                    %compute properties of segmented craters
+                    [num_craters, centx, centy, area, perimeter] = compute_crater_properties(segmentedIm);
+
+
+                    %create mapping structure to be used to convert pixels
+                    %coordinates to latitude and longitude
+
+%                     lat1 = lata; lon1 = lonb;
+%                     lat2 = lata; lon2 = lona;
+%                     lat4 = latb; lon4 = lonb;
+
+                    %if img is coloured convert to gray
+                    if(size(img,3) == 3)
+
+                        img = rgb2gray(img);
+
+                    end
+
+                    R = compute_map_struct(img, lata,lona, latb,lonb);
+
+
+                    %convert centroid pixels to latitude and longitude
+                         %preallocation
+                         Clat = zeros(num_craters, 1);
+                         Clon = zeros(num_craters, 1);
+
+                    for iii = 1:num_craters
+
+                        Centx = centx(iii);
+                        Centy = centy(iii);
+                        [Clat(iii), Clon(iii)] = compute_centroid_latlon(R, Centx, Centy);
+
+
+                    end
+
+                    %output a text file with
+                    %centroidlat[tab]centroidlon[tab]area[tab]perimeter
+                    %XX = [Clat, Clon, area, perimeter];
+                    %dlmwrite(outputname,XX, 'delimiter','\t', 'newline', 'pc');%for windows
+                    %dlmwrite(outputname,XX, 'delimiter','\t', 'newline', 'unix');%for unix
+
+                    %write to spread sheet
+                    %xlswrite(outputname2, XX);
+
+                    T = table(Clat,Clon,area,perimeter);
+                    write(T, outputname, 'Delimiter', '\t');
+                    write(T, outputname2);
+
+                    %superimpose centroids on segmented image
+                    imshow(segmentedIm);
+
+                    hold on
+
+                    plot(centx, centy, 'r*');
+
+                    str2 = strcat('segmented','_', name,'_centroids', ext);
+                    options.Format = 'jpeg';
+                    hgexport(gcf, str2, options);
+
+                    hold off
+            end
+
+    end
+
+

compute_centroid_latlon.m

@@ -0,0 +1,11 @@
+%This function converts centroid (x, and y coordintaes) to latitude and longitude
+%It reads in the structuring element R, as well as the x, and y coordinaates
+
+%ARtillery Crater Analysis and Detection Engine (ARCADE)
+%Developed by Ali M Bukar for Rudiment.info
+%Centre for Visual Computing
+%University of Bradford, UK
+
+function [Clat, Clon] = compute_centroid_latlon(R, centx, centy)
+
+	[Clat, Clon] = pix2latlon(R,centx, centy);

compute_crater_properties.m

@@ -0,0 +1,34 @@
+%This function  uses Matlab regionprops to compute crater-blob properties
+
+%ARtillery Crater Analysis and Detection Engine (ARCADE)
+%Developed by Ali M Bukar for Rudiment.info
+%Centre for Visual Computing
+%University of Bradford, UK
+
+function [num, centx, centy, area, perimeter] = compute_crater_properties(segmentedIm)
+
+
+	[labeledim, num] = bwlabel(segmentedIm,8);%label the blobs(s)
+
+if (num == 0)%if there is no blob, return zeros
+
+        centx = 0; centy =0; area =0; perimeter = 0;
+
+else
+
+        s = regionprops(labeledim, 'Centroid', 'Eccentricity','MajorAxisLength', 'MinorAxisLength', 'Orientation', 'Area','Perimeter');
+
+        centroids = cat(1,s.Centroid);%concatenated centroids
+
+        centx = centroids(:,1);
+
+        centy = centroids(:,2);
+
+        for ii= 1:num
+            area(ii) = s(ii).Area;
+            perimeter(ii) = s(ii).Perimeter;
+        end
+
+        area= area';
+        perimeter = perimeter';
+end

compute_map_struct.m

@@ -0,0 +1,17 @@
+%This function computes a mapping structure used for the conversion of
+%pixel coordinates to latitude and longitudde
+
+%ARtillery Crater Analysis and Detection Engine (ARCADE)
+%Developed by Ali M Bukar for Rudiment.info
+%Centre for Visual Computing
+%University of Bradford, UK
+
+function R = compute_map_struct(img, lat1,lon1, lat2,lon2)
+
+	[m, n] = size(img);
+
+	dlon = (lon2 - lon1)/(1-n);
+
+	dlat =(lat1 - lat2)/(1-m);
+
+	R = makerefmat(lon2,lat1, dlon,dlat);

detect_n_segment_craters.m

@@ -0,0 +1,89 @@
+%This function first detects blast crater shells using cascade object detection
+%Then active contours without edges approach is used for segmentation
+%It reads in an image filename and the detector object
+%It returns cnt: number of craters detected
+%           dI: an image showing annotation of detected craters
+%           I: A binary image showing segmented craters as white and
+%           backgroud in black
+%           im: The original image
+
+%ARtillery Crater Analysis and Detection Engine (ARCADE)
+%Developed by Ali M Bukar for Rudiment.info
+%Centre for Visual Computing
+%University of Bradford, UK
+
+
+
+function [cnt, dI, I, im] = detect_n_segment_craters(filename, detector3, itr)
+
+
+    %load trained detector object it should be found in working directory
+
+   %create detector object
+
+
+
+    %read image url
+
+    im = imread(filename,'jpg');
+
+    I = im;
+
+    if(isa (I, 'double'))
+        I = im2uint8(I);
+    end
+
+    %%convert to gray
+    if(size(I,3) == 3)
+
+        I = rgb2gray(I);
+
+    end
+    bbox = step(detector3, I);
+    I3 = I;
+    I = insertObjectAnnotation(I,'rectangle',bbox,'c');
+
+    dI = I; %annotated images
+
+    %count number of craters
+    cnt = size(bbox,1);
+
+    %define number of iterations for active contour
+    if ~exist('itr','var')|| isempty(itr)
+
+        iter = 500;
+    else
+        iter = itr;
+    end
+
+
+    [m,n] = size(I3);
+
+    cc = zeros(m,n);
+
+    k = size(bbox,1);
+
+
+    %Iteratively segment each annotated image
+    for i = 1:k
+
+        x = bbox(i,:);
+        a = x(2); b= x(2)+x(4); c = x(1); d = x(1)+x(3);
+        cc(a:b, c:d) = I3(a:b, c:d);
+
+        Im = cc(a:b, c:d);
+
+        mask = ones(x(4)+1, x(3)+1);%define a mask
+
+        ck = activecontour(Im,mask,iter,'Chan-Vese');
+        ck = get_largestblob(ck);
+        cc(a:b, c:d) = ck;
+    end
+
+    %convert cc to logical/binary
+
+    I = im2bw(cc); %this is the segmented image
+
+
+
+

get_largestblob.m

@@ -0,0 +1,15 @@
+%This is a helper function used to get largest blob 
+
+%ARtillery Crater Analysis and Detection Engine (ARCADE)
+%Developed by Ali M Bukar for Rudiment.info
+%Centre for Visual Computing
+%University of Bradford, UK
+
+function imbin = get_largestblob(imbin)
+    [labeledIm numL] = bwlabel(imbin);
+    blobM = regionprops(labeledIm, 'area');
+    allAreas = [blobM.Area];
+    %sort descend
+    [sortedAreas, sortind] = sort(allAreas, 'descend');
+    biggestB = ismember(labeledIm, sortind(1));
+    imbin = biggestB>0;