diff --git a/detectRoof.m b/detectRoof.m
new file mode 100644
index 0000000..919d037
--- /dev/null
+++ b/detectRoof.m
@@ -0,0 +1,231 @@
+function [ClustersOut] = detectRoof(ptCloud,h)
+% DETECTROOF
+%
+% Funtion to detect roof vertices from an aerial point cloud. 
+%
+% Inputs: 
+% - ptCloud: point cloud data 
+% - h: mean curvature of each point. Compute using Beksi (2014). If this
+% value is not provided, the function will compute it (may take a while)
+%
+% Outputs:
+% - ClustersOut: a struct containing the planar surfaces of the roof.
+%
+% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
+
+format long g
+
+normals=ptCloud.Normal;
+
+%% COMPUTE CURVATURE IF NOT PROVIDED
+if nargin<2
+    tree = KDTreeSearcher(ptCloud.Location);
+    radius = 1;
+    f=waitbar(0,'Computing curvatures...');
+    c1=zeros(1,ptCloud.Count);
+    c2=zeros(1,ptCloud.Count);
+    h=zeros(1,ptCloud.Count);
+    for i=1:ptCloud.Count
+        query = [ptCloud.Location(i,1) ptCloud.Location(i,2) ...
+            ptCloud.Location(i,3)];
+        [c1(i), c2(i)] = estimateCurvatures(normals, tree, query, radius);
+        h(i) = (c1(i)+c2(i))/2; %Mean curvature; if = 0 means point is plane 
+        waitbar(i/ptCloud.Count,f)
+    end
+    close(f)
+end
+% if needed, save the results
+%savefile = strcat('curvature_detectRoof.mat');
+%save(savefile, 'h');
+%% NORMALS-BASED REGION GROWING
+% detect the roof patches and store it in a struct
+% see also: regiongrowingnormalsOct.m
+[RegionsOct]=regiongrowingnormalsOct(ptCloud,normals,h);
+
+%% DISTANCE-BASED REGION GROWING 
+% this part is used to clean up the resulting patches from noises
+
+% extract patch names
+nbRegions=numel(fieldnames(RegionsOct));
+f=fieldnames(RegionsOct);
+
+% extract number of patches
+RegionCount=nbRegions;
+for i=1:nbRegions
+    RegionName = f{i};
+    labels=pcsegdist(RegionsOct.(RegionName),0.2);
+    uniqueLabels=unique(labels);
+    [nbUniqueLabels,~]=size(uniqueLabels);
+    for j=1:nbUniqueLabels
+        [~, boolInd]=ismember(labels,uniqueLabels(j,1));
+        indexInd=find(boolInd);
+        newRegion=select(RegionsOct.(RegionName),indexInd);
+        
+        % if the resulting patch consists of less than 1000 points, delete
+        % it and do not use for further processing
+        if newRegion.Count>1000
+            RegionCount=RegionCount+1;
+            newRegionName=strcat('Region',string(RegionCount));
+            RegionsOct.(newRegionName)=newRegion;
+        else
+            continue
+        end
+    end
+    RegionsOct=rmfield(RegionsOct,RegionName);
+end
+
+%% EXTRACT THE PLANE MESHES
+f=fieldnames(RegionsOct);
+nbRegions=numel(fieldnames(RegionsOct));
+objectList=strings(nbRegions,1);
+for i=1:nbRegions
+    RegionName = f{i};
+    thisPtCloud=RegionsOct.(RegionName);
+
+    inPtCloud = thisPtCloud;
+    
+    % perform RANSAC to extract the plane
+    [model,inlierIndices,~] = pcfitplane(inPtCloud,0.2);
+    plane = select(inPtCloud,inlierIndices);
+
+    %convert the Matlab surface model from pcfitplane to geom3d format
+    a=model.Parameters(1,1);
+    b=model.Parameters(1,2);
+    c=model.Parameters(1,3);
+    d=model.Parameters(1,4);
+
+    P1 = [plane.Location(1,1) plane.Location(1,2) 0];
+
+    P1(1,3) = (-d-(a*P1(1,1))-(b*P1(1,2)))/c;
+
+    planeGeom3d = createPlane(P1, model.Normal); %plane in geom3d format
+    
+    %Perform 3D Delaunay Triangulation on the point cloud
+    TR = delaunayTriangulation(double(plane.Location(:,1)),...
+        double(plane.Location(:,2)),double(plane.Location(:,3)));
+    
+    %extract the surface of the Delaunay mesh TR (otherwise stored in
+    %tetrahedral form)
+    [~,xf] = freeBoundary(TR);
+    [size_xf,~] = size(xf);
+    
+    %project the points to the mathematical surface model to get 1 surface
+    %in the form of the RANSAC plane
+    prjtdPts = zeros(size_xf,3);
+    
+    for j=1:size_xf
+        %create a 3D ray for each foint in the mesh surface
+        point = [xf(j,1), xf(j,2), xf(j,3)];
+        
+        %intersect the 3D ray with the mathematical surface
+        prjtdPt = projPointOnPlane(point, planeGeom3d);
+        
+        %store the coordinates of the intersection
+        prjtdPts(j,1) = prjtdPt(1,1);
+        prjtdPts(j,2) = prjtdPt(1,2);
+        prjtdPts(j,3) = prjtdPt(1,3);
+    end
+    
+    %perform 3D Delaunay Triangulation on the projected points
+    TR2 = delaunayTriangulation(double(prjtdPts(:,1)),...
+        double(prjtdPts(:,2)),double(prjtdPts(:,3)));
+    
+    %extract the surface of the Delaunay mesh of the projected points
+    [~,xf2] = freeBoundary(TR2);
+    
+    %further smoothing here
+    [xf3,F2] = meshParallelSimp(xf2, 0.1);
+    
+    %store the mesh properties in a struct
+    Mesh.Vertices = xf3;
+    Mesh.Faces = F2;
+    
+    %create a Matlab alphashape
+    shp = alphaShape(xf2(:,1),xf2(:,2),xf2(:,3),inf);
+
+    %store all of this stuff in the another struct called Object
+    objectList(i,1) = strcat('Object',num2str(i));
+    Object.PtCloud = plane;
+    Object.AlphaShp = shp;
+    Object.Mesh = Mesh;
+    Object.PlaneGeom3D = planeGeom3d;
+    Object.PlaneMatlab = model;
+    
+    Clusters.(objectList{i}) = Object;
+end
+
+%% PERFORM 'SNAPPING'
+SnappedCluster = meshSnap(Clusters,5);
+%% REORDER THE VERTICE LIST TO COMPLY WITH CITYGML
+% CityGML vertices must be sorted clockwise from barycentre
+
+figure('name','Detected Roof Vertices - Patch/Polygon')
+for i=1:nbRegions
+    
+    % name of current object
+    ObjectName = objectList{i};
+    
+    % retrieve the object
+    thisObject=SnappedCluster.(ObjectName);
+    
+    % retrieve the vertices list
+    thisVertexList=thisObject.Mesh.Vertices;
+    
+    % initiate a list of bearing angles
+    [nbVertices,~]=size(thisVertexList);
+    listBearings=zeros(nbVertices,2);
+    listBearings2=zeros(nbVertices,2);
+    
+    % compute the Principal Component Analysis
+    coeffs = pca(thisVertexList);
+
+    % transform the point cloud to PC-aligned system i.e. planar surface
+    TransformPCA = thisVertexList*coeffs(:,1:3);
+    
+    % compute the barycenter of the vertices
+    barycenter=[mean(TransformPCA(:,1)),mean(TransformPCA(:,2)), ...
+        mean(TransformPCA(:,3))];
+    
+    % take only the XY (assume the points are coplanar)
+    Xa=barycenter(1);
+    Ya=barycenter(2);
+    
+    for j=1:nbVertices
+        %indexing
+        listBearings(j,1)=j;
+        
+        % retrieve the vertex's XY 
+        Xb=TransformPCA(j,1);
+        Yb=TransformPCA(j,2);
+        
+        % compute the bearing from the barycenter to each vertex
+        [alpha, ~] = bearing_surv(Xa,Ya,Xb,Yb);
+        
+        listBearings(j,2)=alpha;
+    end
+    
+    % sort the vertice list according to this bearing
+    [sorted,index]=sort(listBearings(:,2),'ascend');
+    listBearings2(:,1)=index;
+    listBearings2(:,2)=sorted;
+    
+    % retrieve the original points, but sorted
+    sortedVertices=zeros(nbVertices,3);
+    for j=1:nbVertices
+        sortedVertices(j,:)=thisVertexList(listBearings2(j,1),:);
+    end
+    
+    % update the result with the new sorted values
+     SnappedCluster.(ObjectName).Mesh.Vertices=sortedVertices;
+    
+    patch(sortedVertices(:,1),sortedVertices(:,2),sortedVertices(:,3),'g')
+    view(3)
+    hold on
+    pcshow(SnappedCluster.(ObjectName).PtCloud)
+    hold on
+end
+
+% copy the results to the output struct
+ClustersOut=SnappedCluster;
+end
+
diff --git a/intersectMultPlanes.m b/intersectMultPlanes.m
new file mode 100644
index 0000000..eb30793
--- /dev/null
+++ b/intersectMultPlanes.m
@@ -0,0 +1,27 @@
+function [vectorXYZ,residuals] = intersectMultPlanes(Planes)
+% INTERSECTMULTPLANES
+%
+% Use classical non-weighted least squares to compute the most probable
+% intersection point of multiple 3D planes.
+%
+% Inputs: 
+% - Planes: an Nx4 matrix containing the each plane's parameters (a,b,c,d).
+% N is the number of available planes
+%
+% Outputs:
+% - vectorXYZ: a vector containing the most likely 3D coordinates of the
+% intersection
+% - residuals: an Nx1 matrix containing the residual distances of vectorXYZ
+% to the respective planes
+%
+% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
+
+A=Planes(:,1:3);
+F=Planes(:,4)*-1;
+
+X_LS=((A'*A)^-1)*(A'*F);
+vectorXYZ=X_LS';
+
+residuals=(A*X_LS)-F;
+end
+
diff --git a/meshParallelSimp.m b/meshParallelSimp.m
new file mode 100644
index 0000000..ac071f9
--- /dev/null
+++ b/meshParallelSimp.m
@@ -0,0 +1,144 @@
+function [vertices2,faces2] = meshParallelSimp(vertices1, tol)
+% MESHPARALLELSIMP
+%
+% Function to simplify a surface mesh. The function will detect parallels 
+% in the shape, and take only the vertices of each of the longest vector in
+% the parallel (i.e. the parallel's extremities) as output
+%
+% Inputs: 
+% - vertices1: list of nx3 coordinates of mesh free boundary vertices. The
+% mesh should already be a planar surface (cf. roof detection)
+% - tol: tolerance for parallelism. The closer the value to 0, the better 
+% the original mesh shape is preserved (but the lesser the size reduction)
+%
+% Outputs:
+% - vertices2: list of nx3 coordinates of the simplified mesh
+% - faces2: triangular faces of vertices2 resulting froma a Delaunay
+% triangulation
+%
+% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
+
+% input points
+pts=vertices1;
+[numPts,~]=size(pts);
+indexList=zeros(numPts,1);
+
+for i=1:numPts 
+    indexList(i,1)=i;
+end
+
+% list all possible vector combinations
+VectorList=nchoosek(indexList,2);
+[nbCombinations,~]=size(VectorList);
+
+% create a list: columns 1 and 2 are the vector points; colums 3 to 5 the
+% xyz non-normalised vector
+for i=1:nbCombinations
+    index1=VectorList(i,1);
+    index2=VectorList(i,2);
+    VectorList(i,3)=pts(index1,1)-pts(index2,1);
+    VectorList(i,4)=pts(index1,2)-pts(index2,2);
+    VectorList(i,5)=pts(index1,3)-pts(index2,3);
+end
+%% Check for parallel vectors
+i=1;
+VectorListRef=VectorList;
+
+% keep looping until the VectorList is empty of 'seeds'
+while ~isempty(VectorList)
+    parallelName = strcat('Parallel',string(i));
+    
+    % 'seed' vector
+    v1 = [VectorList(1,3) VectorList(1,4) VectorList(1,5)];
+    [~, index]=ismember(VectorListRef(:,3:5),v1,'rows');
+    indexv1=find(index,1);
+    
+    % normalise the seed vector (v1)
+    v1 = normalizeVector3d(v1);
+    
+    % number of pairs to check the parallelism
+    [pairsToCheck,~]=size(VectorList);
+    
+    % add this seed to the parallel's list of vectors
+    listParallel = indexv1;
+    ListToDelete=[];
+    
+    for j=2:pairsToCheck 
+        % second vector (to be checked against v1)
+        v2 = [VectorList(j,3) VectorList(j,4) VectorList(j,5)];
+        [~, index2]=ismember(VectorListRef(:,3:5),v2,'rows');
+        indexv2=find(index2,1);
+        
+        % normalise v2
+        v2 = normalizeVector3d(v2);
+        
+        % check the disparity between the two vectors.
+        disparity=abs(v1-v2);
+        
+        % magDisparity should be 0 if they are parallel
+        magDisparity=sqrt(disparity(1)^2+disparity(2)^2+disparity(3)^2);
+        
+        % set parallelism tolerance here
+        if magDisparity < tol
+           listParallel=vertcat(listParallel,indexv2);
+           ListToDelete = vertcat(ListToDelete,j);
+        end
+    end
+    
+    % sanity check. don't delete anything if there is nothing to delete
+    if ~isempty(ListToDelete)
+        VectorList(ListToDelete,:) = [];
+    end
+    VectorList(1,:)=[];
+    
+    % store the vector indices belonging to the same parallel to a struct
+    Struct.(parallelName).VectorIndices=listParallel;
+    
+    % moving on...
+    i=i+1;
+end
+%% compute parallel vector distances
+nbFields=length(fieldnames(Struct));
+verticesSimp = [];
+for i=1:nbFields
+    parallelName = strcat('Parallel',string(i));
+    [nbVectors,~]=size(Struct.(parallelName).VectorIndices);
+    for j=1:nbVectors
+        thisVectorIndex=Struct.(parallelName).VectorIndices(j,1);
+        v=[VectorListRef(thisVectorIndex,3) VectorListRef(thisVectorIndex,4) VectorListRef(thisVectorIndex,5)];
+        distance = sqrt(v(1)^2+v(2)^2+v(3)^2);
+       Struct.(parallelName).Distances(j,1)=distance;
+    end
+    % determine the longest vector in a parallel
+    [M,I] = max(Struct.(parallelName).Distances);
+    Struct.(parallelName).LongestVecMag = M;
+    Struct.(parallelName).LongestVecInd = Struct.(parallelName).VectorIndices(I,1);
+    VectorIndex = Struct.(parallelName).VectorIndices(I,1);
+    
+    % extract the corresponding vertices
+    point1Index = VectorListRef(VectorIndex,1);
+    point2Index = VectorListRef(VectorIndex,2);
+    
+    coordPoint1 = pts(point1Index,:);
+    coordPoint2 = pts(point2Index,:);
+    
+    Struct.(parallelName).ExtremePts = vertcat(coordPoint1,coordPoint2);
+    dummy = Struct.(parallelName).ExtremePts;
+    
+    verticesSimp = vertcat(verticesSimp,dummy);
+end
+
+% delete doubles
+verticesSimp = unique (verticesSimp,'rows');
+
+%% delaunay triangulation
+TR = delaunayTriangulation(verticesSimp(:,1),verticesSimp(:,2),verticesSimp(:,3));
+
+[faces2,vertices2] = freeBoundary(TR);
+[numPts2,~] = size(vertices2);
+
+gain=(abs(numPts2-numPts)/numPts)*100;
+
+%disp(strcat('Heyo! Vertices reduced by:',num2str(gain),'%'));
+end
+
diff --git a/meshSnap.m b/meshSnap.m
new file mode 100644
index 0000000..c6d92c9
--- /dev/null
+++ b/meshSnap.m
@@ -0,0 +1,269 @@
+function [NewClusterStruct] = meshSnap(inputCluster,tol)
+% MESHSNAP
+%
+% Function to merge nearby points within a set of meshes located within a
+% distance of TOL.
+%
+% Inputs: 
+% - inputCluster: a struct containing Objects, containing the point cloud
+% and mesh information. The Objects are separate entities; the function
+% will however snap the mesh vertices by taking into account the
+% global/merged state of all Objects within the inputCluster
+% - tol: distance threshold between a point with its neighbours to be
+% considered as the same point
+%
+% Outputs:
+% - vertices2: list of nx3 coordinates of the simplified mesh
+% - faces2: triangular faces of vertices2 resulting froma a Delaunay
+% triangulation
+%
+% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
+
+NewClusterStruct=inputCluster;
+
+%% MERGE THE MESHES
+nbRegions=numel(fieldnames(inputCluster));
+objectList=strings(nbRegions,1);
+
+% create a list containing Object names (useful later)
+for j=1:nbRegions
+    objectList(j,1)=strcat('Object',num2str(j));
+    [nbVertices,~]=size(inputCluster.(objectList(j,1)).Mesh.Vertices);
+    
+    % create a new field with index to the original object
+    v = zeros(nbVertices, 1);
+    v(:) = j;
+    inputCluster.(objectList(j,1)).Mesh.ObjectIndex = v;
+end
+
+% initialise variables
+MergedVertices=inputCluster.Object1.Mesh.Vertices;
+MergedFaces=inputCluster.Object1.Mesh.Faces;
+MergedIndex=inputCluster.Object1.Mesh.ObjectIndex;
+i=1;
+
+% start the loop
+while i<nbRegions
+    objectName1=objectList(i);
+    objectName2=objectList(i+1);
+    
+    % stack the vertex lists
+    MergedVertices = vertcat(MergedVertices, inputCluster.(objectName2).Mesh.Vertices);
+
+    % update the face list
+    faceUpdate = inputCluster.(objectName2).Mesh.Faces+max(max(MergedFaces));
+    MergedFaces = vertcat(MergedFaces,faceUpdate);
+    
+    % update the object index list
+    MergedIndex = vertcat(MergedIndex, inputCluster.(objectName2).Mesh.ObjectIndex);
+    
+    i=i+1;
+end
+
+%% COMPUTE "SNAPPING"
+%initialise the variables
+ClustersCopy=inputCluster; % "dynamic" copy of the original Cluster struct
+VertexList=MergedVertices; % "dynamic" list of vertices
+MergedIndex2=MergedIndex; %"dynamic" list of indices
+
+% start the loop, continue as long as the list of vertices is not empty
+while ~isempty(VertexList)
+    
+    % take the first seed
+    vertex1 = VertexList(1,:);
+    
+    % change VertexList to Matlab point cloud struct, necessary for the
+    % next function
+    VertexPointCloud=pointCloud(VertexList);
+    
+    % find the seed's 10 nearest neighbours (this value can be changed, but
+    % I doubt it would be necessary)
+    [indices,dists] = findNearestNeighbors(VertexPointCloud,vertex1,10);
+    
+    % for now, delete the seed from the list of nearest neighbours
+    indices(1,:)=[];
+    dists(1,:)=[];
+    
+    % number of neighbours inside the tolerance radius
+    nb = length(indices);
+    
+    % initialise a dynamic list of neighbours and their pointer towards the
+    % original Object
+    NeighborList = []; % empty list of neighbours inside the radius
+    OriginList = []; % empty list of neighbours' origins
+  
+    % start a loop to check each neighbour
+    for i=1:nb
+        % if the neighbour point's distance is less than the tolerance,
+        % consider it as the same point i.e. "snapped"
+        if dists(i,1)<tol
+            a=MergedIndex2;
+            NeighborList(i,1)=indices(i,1);
+            OriginList(i,1)= a(indices(i,1));
+        else
+            continue
+        end
+    end
+    
+    % if the seed has neighbours...
+    if ~isempty(NeighborList)
+        % return the seed to both lists, necessary to compute the
+        % coordinates of the new "snapped" point
+        NeighborList=[NeighborList;1];
+        OriginList=[OriginList;MergedIndex2(1)];
+        
+        % number of items in the list
+        szList=length(NeighborList);
+
+        % how many unique Objects?
+        UniqueOrigins=unique(OriginList);
+        szUnique=length(UniqueOrigins);
+
+        % create a recapitulative table
+        tableToDelete=zeros(szList,2);
+        for x=1:szList
+            originNm=strcat('Object',string(OriginList(x)));
+    
+            % retrieve the vector to be deleted from reference list of 
+            % vertices
+            VectToDelete =  VertexList(NeighborList(x),:);
+    
+            % find its index in the original Object vertices
+            [~, boolDel]=ismember(ClustersCopy.(originNm).Mesh.Vertices,...
+                VectToDelete);
+            indexDel=find(boolDel,1);
+    
+            % update the recapitulative table
+            tableToDelete(x,1) = OriginList(x);
+            tableToDelete(x,2) = indexDel;
+        end
+        
+        % delete the points whose indices are listed in NeighborList in 
+        % respective original Objects
+        for x=1:szUnique
+            originNm=strcat('Object',string(UniqueOrigins(x)));
+    
+            % find the indices of This Object to be deleted from the 
+            % recapitulative table
+            ind1 = tableToDelete(:,1) == UniqueOrigins(x);
+            thisOriginDel = tableToDelete(ind1,2);
+    
+            % delete the said point in the dynamic/copied cluster
+            ClustersCopy.(originNm).Mesh.Vertices(thisOriginDel,:)=[];
+    
+            % keep only unique points and delete (possible) redundants
+            ClustersCopy.(originNm).Mesh.Vertices=...
+            unique(ClustersCopy.(originNm).Mesh.Vertices,'rows');
+        end
+        
+        % initialise XYZ vectors of all the neighbours and the seed 
+        x_vect=zeros(szList,1);
+        y_vect=zeros(szList,1);
+        z_vect=zeros(szList,1);
+        
+        % retrieve the neighbours and seed's XYZ
+        for j=1:szList
+            x_vect(j)=VertexList(NeighborList(j),1);
+            y_vect(j)=VertexList(NeighborList(j),2);
+            z_vect(j)=VertexList(NeighborList(j),3);          
+        end
+    
+        % compute the median to generate a new point
+        % why median? more robust than mean/average
+        x_new = median(x_vect);
+        y_new = median(y_vect);
+        z_new = median(z_vect);
+    
+        new = [x_new y_new z_new];
+        
+        % add planarity constraint
+        listPlanes=zeros(szUnique,9);
+        listPlanesParameters=zeros(szUnique,4);
+        
+        % if the neighbours belong to more than 2 planes...
+        if szUnique>2
+            for k=1:szUnique
+                originNm = strcat('Object',string(UniqueOrigins(k)));
+                listPlanesParameters(k,:)=...
+                    inputCluster.(originNm).PlaneMatlab.Parameters;
+            end
+            
+            % use Geom3D's three planes intersection function
+            %  plane1 = listPlanes(1,:);
+            %  plane2 = listPlanes(2,:);
+            %  plane3 = listPlanes(3,:);
+            %  new2 = intersectThreePlanes(plane1, plane2, plane3);
+
+            % use my least squares multiple planes function
+            [new2,~]=intersectMultPlanes(listPlanesParameters);
+        
+        % if the neighbours belong to two planes...
+        elseif szUnique==2
+            for k=1:2
+                originNm = strcat('Object',string(UniqueOrigins(k)));
+                listPlanes(k,:)=inputCluster.(originNm).PlaneGeom3D;
+            end
+            plane1 = listPlanes(1,:);
+            plane2 = listPlanes(2,:);
+            
+            % intersect the two planes to get a 3D line
+            line = intersectPlanes(plane1, plane2);
+            
+            % project the computed median point into this 3D line
+            new2 = projPointOnLine3d(new, line);
+        else
+            new2=new;
+        end
+        
+        % stact the new point to the un-merged objects
+        if szUnique>1
+            for x=1:szUnique
+                originNm = strcat('Object',string(UniqueOrigins(x)));
+                ClustersCopy.(originNm).Mesh.Vertices= ...
+                    [ClustersCopy.(originNm).Mesh.Vertices;new2];
+            end
+        else
+            originNm = strcat('Object',string(UniqueOrigins(1)));
+            ClustersCopy.(originNm).Mesh.Vertices=...
+                [ClustersCopy.(originNm).Mesh.Vertices;new2];
+        end
+        
+        % delete the snapped points from the two lists
+        indNeighbors=NeighborList';
+        VertexList(indNeighbors,:)=[];
+        MergedIndex2(indNeighbors,:)=[];
+    else        
+        % delete the seed from the two lists if they don't have neighbours   
+        VertexList(1,:)=[];
+        MergedIndex2(1,:)=[];
+    end
+end
+
+%% VISUALISE THE 'SNAPPED' RESULT AND STORE IN NEW STRUCT
+figure('name','Detected Roof Vertices - Delaunay Mesh')
+for x=1:nbRegions
+    originNm=strcat('Object',string(x));
+    
+    % Delaunay triangulation
+    TR = delaunayTriangulation...
+        (double(ClustersCopy.(originNm).Mesh.Vertices(:,1)),...
+        double(ClustersCopy.(originNm).Mesh.Vertices(:,2)),...
+        double(ClustersCopy.(originNm).Mesh.Vertices(:,3)));
+    
+    % retrieve the free boundary mesh
+    [F,xf] = freeBoundary(TR);
+    [V2, F2] = ensureManifoldMesh(xf, F);
+    
+    % store the new vertices and faces in a new struct
+    NewClusterStruct.(originNm).Mesh.Vertices=V2;
+    NewClusterStruct.(originNm).Mesh.Faces=F2;
+    
+    % plot the mesh
+    pcshow(ClustersCopy.(originNm).PtCloud)
+    hold on
+    drawMesh(V2,F2,'b');
+    hold on
+end
+
+end
+
diff --git a/ptCloudCenter.m b/ptCloudCenter.m
new file mode 100644
index 0000000..4413e20
--- /dev/null
+++ b/ptCloudCenter.m
@@ -0,0 +1,41 @@
+function [centeredPtCloud,translation] = ptCloudCenter(ptCloud,X0,Y0,Z0)
+% PTCLOUDCENTER
+%
+% Funtion to center point clouds with large coordinates (typically coor-
+% dinates from projected systems)
+%
+% Inputs: 
+% - ptCloud: point cloud data 
+% - X0,Y0,Z0: translation values. If these values are not available, the
+% function will use the barycentre instead to translate the point cloud.
+%
+% Outputs:
+% - recenteredPtCloud: point cloud with recentred coordinates.
+% - translation: XYZ vector with the translation values.
+%
+% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
+
+sensorCenter=[mean(ptCloud.Location(:,1)),mean(ptCloud.Location(:,2)), ...
+     mean(ptCloud.Location(:,3))];
+ptRGB=ptCloud.Color;
+ptNormals=ptCloud.Normal;
+ptIntensity=ptCloud.Intensity;
+
+if nargin<2
+    translation(1)=sensorCenter(1);
+    translation(2)=sensorCenter(2);
+    translation(3)=sensorCenter(3);
+else
+    translation(1)=X0;
+    translation(2)=Y0;
+    translation(3)=Z0;
+end
+
+x=ptCloud.Location(:,1)-translation(1);
+y=ptCloud.Location(:,2)-translation(2);
+z=ptCloud.Location(:,3)-translation(3);
+
+centeredPtCloud=pointCloud([x y z],'Color',ptRGB,'Normal',ptNormals,...
+    'Intensity',ptIntensity);
+end
+
diff --git a/ptCloudRecenter.m b/ptCloudRecenter.m
new file mode 100644
index 0000000..24230f0
--- /dev/null
+++ b/ptCloudRecenter.m
@@ -0,0 +1,28 @@
+function [recenteredPtCloud] = ptCloudRecenter(ptCloud,translation)
+% PTCLOUDRECENTER
+%
+% Funtion to recentre point clouds previously centered using PTCLOUDCENTER.
+% This function returns the centred point cloud to its original coordinate 
+% system.
+%
+% Inputs: 
+% - ptCloud: point cloud data 
+% - translation: translation vector from PTCLOUDCENTER
+%
+% Outputs:
+% - recenteredPtCloud: point cloud with recentred coordinates.
+%
+% (c) Arnadi Murtiyoso (INSA Strasbourg - ICube-TRIO UMR 7357)
+
+ptRGB=ptCloud.Color;
+ptNormals=ptCloud.Normal;
+ptIntensity=ptCloud.Intensity;
+
+x=ptCloud.Location(:,1)+translation(1);
+y=ptCloud.Location(:,2)+translation(2);
+z=ptCloud.Location(:,3)+translation(3);
+
+recenteredPtCloud=pointCloud([x y z],'Color',ptRGB,'Normal',ptNormals,...
+    'Intensity',ptIntensity);
+end
+