function [domains boundaries interiorBoundaries] = geomAnalyze(geom)

    objects = geom.objectNames();

    objCoords = zeros(objects.length,2);
    for i=1:objects.length
        obj = geom.object(objects(i));
        coords = obj.getVertexCoord;
        
        firstCoord = coords(:,1);
        for j=2:length(coords)
            if coords(1,j) < firstCoord(1) || (coords(1,j) == firstCoord(1) && coords(2,j) < firstCoord(2))
                firstCoord = coords(:,j);
            end
        end
        objCoords(i,1)= firstCoord(1);
        objCoords(i,2)= firstCoord(2);
        
    end
    
    taken = zeros(1,objects.length);
    clear domains; 
    for k=1:objects.length
        
        selectedCoord = 0;
        selectedIndex = 0;
        for i=1:length(objCoords)
            if taken(i)==0
                if selectedIndex == 0 || objCoords(i,1) < selectedCoord(1) || (objCoords(i,1) == selectedCoord(1) && objCoords(i,2) < selectedCoord(2))
                   selectedCoord = objCoords(i,:);
                   selectedIndex = i;
                end
            end
        end
        
        taken(selectedIndex) = 1;
        domains(k) = objects(selectedIndex);
    end
    
    clear boundaryCoords boundaries; k=1;
    for i=1:domains.length
        obj = geom.object(domains(i));
        edges = obj.getEdge; edges=edges(1:2,:);
        vertices = obj.getVertex;

        for j = 1:length(edges)
                boundaryCoords(1,k) = vertices(1,edges(1,j));
                boundaryCoords(2,k) = vertices(2,edges(1,j));
                boundaryCoords(3,k) = vertices(1,edges(2,j));
                boundaryCoords(4,k) = vertices(2,edges(2,j));
                boundaryCoords(5,k) = i;                
                k=k+1;
        end
    end
    
    a=geom.feature.create('all','Union');
    a.selection('input').set( geom.objectNames) ; 
    a.set('keep','on');
    geom.run; 
    
    finalEdges = geom.object('all').getEdge;
    finalVertices = geom.object('all').getVertex;
    interiorBoundaries = zeros(1,length(finalEdges));    
    for i=1:length(finalEdges)
        
        vertA = finalVertices(1:2, finalEdges(1,i) );
        vertB = finalVertices(1:2, finalEdges(2,i) );
        vertA = round(vertA.*10^10)./10^10;
        vertB = round(vertB.*10^10)./10^10;
        
        interior = 0;
        bestCandidate = 0;
        bestPerc = 0;
        for j=1:length(boundaryCoords)
            boundaryCoords = round(boundaryCoords.*10^10)./10^10;
            ol = overlapLength([boundaryCoords(1,j),boundaryCoords(3,j)],[boundaryCoords(2,j),boundaryCoords(4,j)],[vertA(1),vertB(1)],[vertA(2),vertB(2)]);
            if ol > 0
                interior = (bestPerc > 0);
                if (ol > bestPerc)
                    bestPerc=ol;
                    bestCandidate=j;
                end
            end
            
        end
            
        if (bestCandidate == 0)
            fprintf('couldnt find boundary for %d\n',i);
        else
            boundaries(i) =domains( boundaryCoords(5,bestCandidate) );
            if (interior ~= 0)
                interiorBoundaries(i) = 1;
            end
        end
    end
    geom.feature().remove('all');
    geom.run;
end