diff --git a/EBSDAnalysis/@EBSD/calcGrains.m b/EBSDAnalysis/@EBSD/calcGrains.m
index 5463bb80e..cd74a5b47 100644
--- a/EBSDAnalysis/@EBSD/calcGrains.m
+++ b/EBSDAnalysis/@EBSD/calcGrains.m
@@ -62,6 +62,7 @@
 % See also
 % GrainReconstruction GrainReconstructionAdvanced
 
+% (1) 
 % subdivide the domain into cells according to the measurement locations,
 % i.e. by Voronoi tessellation or unit cell
 if isa(ebsd,'EBSDsquare') || isa(ebsd,'EBSDhex')
@@ -75,6 +76,7 @@
 % D - cell array of cells
 % I_FD - incidence matrix faces to vertices
 
+% (2)
 % determine which cells to connect
 [A_Db,I_DG] = doSegmentation(I_FD,ebsd,varargin{:});
 % A_db - neighboring cells with (inner) grain boundary
@@ -94,15 +96,58 @@
   spdiags(ebsd.phaseId,0,length(ebsd),length(ebsd)),[],2));
 phaseId(phaseId==0) = 1; % why this is needed?
 
+% (3)
 % compute boundary this gives
 % I_FDext - faces x cells external grain boundaries
 % I_FDint - faces x cells internal grain boundaries
-[I_FDext, I_FDint, Fext, Fint] = calcBoundary;
+% [I_FDext, I_FDint, Fext, Fint] = calcBoundary;
 
+% TODO:
+% should return
+[isExt, isInt] = calcBoundary(A_Db, I_DG, I_FD);
+
+isIntOrExt = isInt | isExt;
+
+% reduce F
+F = F(isIntOrExt,:);
+I_FD = I_FD(isIntOrExt,:);
+isInt = isInt(isIntOrExt);
+isExt = isExt(isIntOrExt);
+
+% reduce V
+[inUse,~,F] = unique(F(:));
+V = V(inUse,:);
+F = reshape(F,[],2);
+
+%(4)
+% detect and sort segments
+% sP  - list of segment points
+% s   - starting point of segements in sP
+% sfF - indices to F of the new face order
+[s, sP, sf, sfF, sfD] = TriplepointSegments(F,V);
+
+
+% * reorder F
+% * split F into Fext and Fint
+% * reorder and reduce I_FD and then split into I_FDext and I_FD_int
+
+F = F(sfF,:);
+isExt = isExt(sfF);
+isInt = isInt(sfF);
+F(sfD == 0,[1,2]) = F(sfD == 0,[2,1]);% flip direction
+I_FD = I_FD(sfF,:);
+
+Fext = F(isExt,:); % external boundary segments
+Fint = F(isInt,:); % internal boundary segments
+I_FDext = I_FD(isExt,:);
+I_FDint = I_FD(isInt,:);
+
+% does this still work? TODO
 if check_option(varargin,'removeQuadruplePoints')
   qAdded = removeQuadruplePoints; 
 end
 
+
 % setup grains
 try
   poly = calcPolygonsC(I_FDext * I_DG,Fext,V);
@@ -258,48 +303,6 @@
 
   end
 
-  function [I_FDext, I_FDint, Fext, Fint] = calcBoundary
-    % distinguish between interior and exterior grain boundaries
-    
-    % cells that have a subgrain boundary, i.e. a boundary with a cell
-    % belonging to the same grain
-    sub = ((A_Db * I_DG) & I_DG)';                 % grains x cell
-    [i,j] = find( diag(any(sub,1))*double(A_Db) ); % all adjacent to those
-    sub = any(sub(:,i) & sub(:,j),1);              % pairs in a grain
-    
-    % split grain boundaries A_Db into interior and exterior
-    A_Db_int = sparse(i(sub),j(sub),1,size(I_DG,1),size(I_DG,1));
-    A_Db_ext = A_Db - A_Db_int;                    % adjacent over grain boundary
-    
-    % create incidence graphs
-    I_FDbg = diag( sum(I_FD,2)==1 ) * I_FD;
-    D_Fbg  = diag(any(I_FDbg,2));
-    
-    [ix,iy] = find(A_Db_ext);
-    D_Fext  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
-    
-    I_FDext = (D_Fext| D_Fbg)*I_FD;
-    
-    [ix,iy] = find(A_Db_int);
-    D_Fsub  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
-    I_FDint = D_Fsub*I_FD;
-    
-    % remove empty lines from I_FD, F, and V
-    isExt = full(any(I_FDext,2));
-    I_FDext = I_FDext.'; I_FDext = I_FDext(:,isExt).';
-
-    isInt = full(any(I_FDint,2));
-    I_FDint = I_FDint.'; I_FDint = I_FDint(:,isInt).';
-      
-    % remove vertices that are not needed anymore
-    [inUse,~,F] = unique(F(isExt | isInt,:));
-    V = V(inUse,:);
-    F = reshape(F,[],2);
-    Fext = F(isExt(isExt | isInt),:); % external boundary segments
-    Fint = F(isInt(isExt | isInt),:); % internal boundary segments
-
-  end
-
   function gB = makeBoundary(F,I_FD)
 
     % compute ebsdInd
@@ -449,6 +452,48 @@
 
 end
 
+function [isExt, isInt] = calcBoundary(A_Db, I_DG, I_FD)
+    % distinguish between interior and exterior grain boundaries
+    
+    % cells that have a subgrain boundary, i.e. a boundary with a cell
+    % belonging to the same grain
+    sub = ((A_Db * I_DG) & I_DG)';                 % grains x cell
+    [i,j] = find( diag(any(sub,1))*double(A_Db) ); % all adjacent to those
+    sub = any(sub(:,i) & sub(:,j),1);              % pairs in a grain
+    
+    % split grain boundaries A_Db into interior and exterior
+    A_Db_int = sparse(i(sub),j(sub),1,size(I_DG,1),size(I_DG,1));
+    A_Db_ext = A_Db - A_Db_int;                    % adjacent over grain boundary
+    
+    % create incidence graphs
+    I_FDbg = diag( sum(I_FD,2)==1 ) * I_FD;
+    D_Fbg  = diag(any(I_FDbg,2));
+    
+    [ix,iy] = find(A_Db_ext);
+    D_Fext  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
+    
+    I_FDext = (D_Fext| D_Fbg)*I_FD;
+    
+    [ix,iy] = find(A_Db_int);
+    D_Fsub  = diag(sum(abs(I_FD(:,ix)) & abs(I_FD(:,iy)),2)>0);
+    I_FDint = D_Fsub*I_FD;
+    
+    % remove empty lines from I_FD, F, and V
+    isExt = full(any(I_FDext,2));
+    %I_FDext = I_FDext.'; I_FDext = I_FDext(:,isExt).';
+
+    isInt = full(any(I_FDint,2));
+    %I_FDint = I_FDint.'; I_FDint = I_FDint(:,isInt).';
+      
+    % remove vertices that are not needed anymore
+    %[inUse,~,F] = unique(F(isExt | isInt,:));
+    %V = V(inUse,:);
+    %F = reshape(F,[],2);
+    %Fext = F(isExt(isExt | isInt),:); % external boundary segments
+    %Fint = F(isInt(isExt | isInt),:); % internal boundary segments
+
+  end
+
 function poly = calcPolygons(I_FG,F,V)
 %
 % Input:
diff --git a/mex/TriplepointSegments.mexw64 b/mex/TriplepointSegments.mexw64
new file mode 100644
index 000000000..28f9a1bea
Binary files /dev/null and b/mex/TriplepointSegments.mexw64 differ
diff --git a/mex/mex_install.m b/mex/mex_install.m
index 6e24771b8..3896fe50d 100644
--- a/mex/mex_install.m
+++ b/mex/mex_install.m
@@ -13,8 +13,8 @@ function mex_install(varargin)
   'extern/jcvoronoi/',...
   'extern/libDirectional/numerical',...
   'SO3Fun/@SO3FunHarmonic/private/wignerTrafo',...
-  'tools/graph_tools/EulerCyclesC'
-  };
+  'tools/graph_tools/EulerCyclesC',...
+  'tools/graph_tools/TriplepointSegments'};
   
 % TODO: Check for mex-Compiler
 
diff --git a/tools/graph_tools/TriplepointSegments.c b/tools/graph_tools/TriplepointSegments.c
new file mode 100644
index 000000000..56c261480
--- /dev/null
+++ b/tools/graph_tools/TriplepointSegments.c
@@ -0,0 +1,291 @@
+
+// mex TriplepointSegments.c -R2018a
+
+#include <string.h>
+#include <stdarg.h>
+#include "mex.h"
+
+// maximum number of neighbors of a vertex
+#define MAX_NEIGHBORS ( 8 )
+
+#define ERROR_ID_PREFIX ( "MTEX:TriplepointSegments" )
+
+// ------------ structs and types ------------
+
+struct IntStack {
+    int* buffer;
+    int index;
+    int maxSize;
+};
+typedef struct IntStack IntStack;
+
+struct WorkingData {
+    double *F_doubles;
+    int F_M;// number of faces
+
+    int V_M;// number of vertices
+
+    int *VN;// VN[i*MAX_NEIGHBORS...] = [neigbor1ofPoint_i, neighbor2ofPoint_i, ...]
+    int *VNFindex;// VNFindex[i*MAX_NEIGHBORS...] = [faceToNeighbor1ofPoint_i, faceToNeighbor2ofPoint_i, ...]
+    int *VNSize;// VNSize[i] = number of neighbors of point i
+
+    IntStack startPoints;
+
+    IntStack s;// segments
+    IntStack sP;// segmentPoints (indice)
+    IntStack sf;// segments (faces)
+    IntStack sfF;// segmentFaces (indices)
+    IntStack sfD;// segment face directions
+
+    int* F_used;
+};
+typedef struct WorkingData WorkingData;
+
+
+// ------------ function headers ------------
+
+void calcTriplepointSegments(const mxArray *F, const mxArray *V, mxArray **s, mxArray **sP, mxArray **sf, mxArray **sfF, mxArray **sfD);
+
+void buildNeighborList(WorkingData* w);
+
+void calcSegmentsStartingFrom(WorkingData *w, int startPoint);
+
+void calcSegmentStartingFrom(WorkingData *w, int startPoint, int nIndex);
+
+static inline void push(IntStack *intStack, int value);
+
+void initIntStack(IntStack *intStack, int maxSize);
+
+void error(const char *id, const char *format, ...);
+
+
+
+
+
+// ------------ function implementations ------------
+
+void mexFunction(int nOutputs, mxArray *outputs[], int nInputs, const mxArray *inputs[]) {
+    if (nInputs != 2) {
+        error("nInputs","2 Inputs required! nInputs=%d",nInputs);
+    }
+    if (nOutputs != 5) {
+        error("nOutputs","5 Outputs required! nOutputs=%d",nOutputs);
+    }
+
+    calcTriplepointSegments(inputs[0], inputs[1], &outputs[0], &outputs[1], &outputs[2], &outputs[3], &outputs[4]);
+}
+
+void calcTriplepointSegments(const mxArray *F, const mxArray *V, mxArray **s, mxArray **sP, mxArray **sf, mxArray **sfF, mxArray **sfD) {    
+    WorkingData w;
+    
+    // ------ read input data -----
+
+    w.F_M = mxGetM(F);
+    w.F_doubles = mxGetDoubles(F);
+
+    w.V_M = mxGetM(V);// number of vertices
+
+    // ------ build data structures ------
+
+    buildNeighborList(&w);
+
+    initIntStack(&(w.startPoints), w.V_M);
+    
+    // calculate if segment start point
+    for (int i = 0; i < w.V_M; i++) {
+        if (w.VNSize[i] != 2 && w.VNSize[i] != 0) {
+            // if number of neighbours not 2 or 0: start point
+            push(&(w.startPoints), i);
+            //printf("segmentstart: %d (%d)\n",i,w.VNSize[i]);
+        }
+    }
+    
+    initIntStack(&(w.s), w.F_M + 1);
+    initIntStack(&(w.sP), w.F_M*2);
+    initIntStack(&(w.sf), w.F_M + 1);
+    initIntStack(&(w.sfF), w.F_M);
+    initIntStack(&(w.sfD), w.F_M);
+
+    w.F_used = mxMalloc(sizeof(int)*w.F_M);
+    memset(w.F_used, 0, sizeof(int)*w.F_M);
+
+    // triple points
+    for (int i = 0; i < w.startPoints.index; i++) {
+        calcSegmentsStartingFrom(&w, w.startPoints.buffer[i]);
+    }
+
+    // remaining: circles
+    for (int i = 0; i < w.F_M; i++) {
+        if (w.F_used[i] == 0) {
+            calcSegmentStartingFrom(&w, (int)(w.F_doubles[i] - 1), 0);
+        }
+    }
+
+    // close data structure
+    push(&(w.s), w.sP.index);
+    push(&(w.sf), w.sfF.index);
+    
+    // ------ generate output data ------
+
+    *s = mxCreateDoubleMatrix(w.s.index, 1, mxREAL);
+    double* s_doubles = mxGetDoubles(*s);
+    for (int i = 0; i < w.s.index; i++) {
+        s_doubles[i] = w.s.buffer[i] + 1;
+    }
+
+    *sP = mxCreateDoubleMatrix(w.sP.index, 1, mxREAL);
+    double* sP_doubles = mxGetDoubles(*sP);
+    for (int i = 0; i < w.sP.index; i++) {
+        sP_doubles[i] = w.sP.buffer[i] + 1;
+    }
+
+    *sf = mxCreateDoubleMatrix(w.sf.index, 1, mxREAL);
+    double* sf_doubles = mxGetDoubles(*sf);
+    for (int i = 0; i < w.sf.index; i++) {
+        sf_doubles[i] = w.sf.buffer[i] + 1;
+    }
+
+    *sfF = mxCreateDoubleMatrix(w.sfF.index, 1, mxREAL);
+    double* sfF_doubles = mxGetDoubles(*sfF);
+    for (int i = 0; i < w.sfF.index; i++) {
+        sfF_doubles[i] = w.sfF.buffer[i] + 1;
+    }
+
+    *sfD = mxCreateDoubleMatrix(w.sfD.index, 1, mxREAL);
+    double* sfD_doubles = mxGetDoubles(*sfD);
+    for (int i = 0; i < w.sfD.index; i++) {
+        sfD_doubles[i] = w.sfD.buffer[i];
+    }
+}
+
+void buildNeighborList(WorkingData* w) {
+    // allocate buffer
+    w->VN = (int *) mxMalloc(sizeof(int)*w->V_M*MAX_NEIGHBORS);// neighbor vertices
+    w->VNFindex = (int *) mxMalloc(sizeof(int)*w->V_M*MAX_NEIGHBORS);// neighbor faces
+    w->VNSize = (int *) mxMalloc(sizeof(int)*w->V_M);// number of neighbors
+
+    // fill with zeros
+    memset(w->VN, 0, sizeof(int)*w->V_M*MAX_NEIGHBORS);
+    memset(w->VNFindex, 0, sizeof(int)*w->V_M*MAX_NEIGHBORS);
+    memset(w->VNSize, 0, sizeof(int)*w->V_M);
+    
+    // for each face
+    for (int i = 0; i < w->F_M; i++) {
+        int vStart = (int)w->F_doubles[i] - 1;
+        int vEnd = (int)w->F_doubles[i + w->F_M] - 1;
+
+        // validate face
+        if (vStart < 0 || vStart >= w->V_M) {
+            error("F_V_out_of_bounds", "face contains node index out of bounds: vStart=%d is not between 1,%d", vStart+1, w->F_M);
+        }
+        if (vEnd < 0 || vEnd >= w->V_M) {
+            error("F_V_out_of_bounds", "face contains node index out of bounds: vEnd=%d is not between 1,%d", vEnd+1, w->F_M);
+        }
+
+        // prevent buffer overflow
+        if (w->VNSize[vStart] >= MAX_NEIGHBORS || w->VNSize[vEnd] >= MAX_NEIGHBORS) {
+            error("VN_max","A point has too many neighbors! MAX_NEIGHBORS=%d", MAX_NEIGHBORS);
+        }
+
+        w->VN[vStart*MAX_NEIGHBORS + w->VNSize[vStart]] = vEnd;// add neighbor
+        w->VNFindex[vStart*MAX_NEIGHBORS + w->VNSize[vStart]] = i;// set corresponding neighbor face
+        w->VNSize[vStart] += 1;
+
+        if (vStart != vEnd) {
+            // if no loop
+            w->VN[vEnd*MAX_NEIGHBORS + w->VNSize[vEnd]] = vStart;// add neighbor
+            w->VNFindex[vEnd*MAX_NEIGHBORS + w->VNSize[vEnd]] = i;// set corresponding neighbor face
+            w->VNSize[vEnd] += 1;
+        }
+    }
+}
+
+void calcSegmentsStartingFrom(WorkingData *w, int startPoint) {
+    //printf("calcSegmentsStartingFrom(w, %d)\n",startPoint);
+    for (int i = 0; i < w->VNSize[startPoint]; i++) {
+        int face = w->VNFindex[startPoint*MAX_NEIGHBORS + i];
+        if (w->F_used[face] == 0) {
+            calcSegmentStartingFrom(w, startPoint, i);
+        }
+    }
+}
+
+void calcSegmentStartingFrom(WorkingData *w, int startPoint, int nIndex) {
+    //printf("calcSegmentStartingFrom(w, %d, %d)\n",startPoint,nIndex);
+    push(&(w->s), w->sP.index);
+    push(&(w->sf), w->sfF.index);
+
+    // first Point
+    push(&(w->sP), startPoint);
+    //printf("first Point: %d\n",startPoint);
+
+    // first Face
+    int face = w->VNFindex[startPoint*MAX_NEIGHBORS + nIndex];
+    int faceDirection = (int)(w->F_doubles[face]) == startPoint + 1;
+    push(&(w->sfF), face);
+    push(&(w->sfD), faceDirection);
+
+    w->F_used[face] = 1;// use face
+
+    // next Points:
+    int currentPoint = w->VN[startPoint*MAX_NEIGHBORS + nIndex];
+    while (1) {
+        push(&(w->sP), currentPoint);
+        //printf("currentPoint: %d\n",currentPoint);
+
+        if (w->VNSize[currentPoint] != 2) {
+            // if point has not 2 neigbours: end of segment
+            return;
+        }
+
+        face = w->VNFindex[currentPoint*MAX_NEIGHBORS];// first option
+        nIndex = 0;
+
+        if (w->F_used[face]) {
+            face = w->VNFindex[currentPoint*MAX_NEIGHBORS + 1];// second option
+            nIndex = 1;
+        }
+
+        if (w->F_used[face]) {
+            return;// if both faces already used: stop
+        }
+
+        faceDirection = (int)(w->F_doubles[face]) == currentPoint + 1;
+        push(&(w->sfF), face);
+        push(&(w->sfD), faceDirection);
+
+        w->F_used[face] = 1;// use face
+
+        currentPoint = w->VN[currentPoint*MAX_NEIGHBORS + nIndex];// next point
+    }
+}
+
+static inline void push(IntStack *intStack, int value) {
+    //printf("push(%d) at index=%d maxSize=%d\n",value,index,maxSize);
+    if (intStack->index >= intStack->maxSize) {
+        error("internal:IntStackBufferFull", "internal Error: IntStack buffer full! intStack->maxSize=%d",intStack->maxSize);
+    }
+    intStack->buffer[intStack->index] = value;
+    intStack->index += 1;
+}
+
+void initIntStack(IntStack *intStack, int maxSize) {
+    intStack->maxSize = maxSize;
+    intStack->buffer = mxMalloc(sizeof(int)*maxSize);
+    intStack->index = 0;
+}
+
+void error(const char *id, const char *format, ...) {
+    char errorMessage[1024];
+    char fullId[128];
+
+    snprintf(fullId, 128, "%s:%s", ERROR_ID_PREFIX, id);
+    
+    // ---- snprintf(fullId, 128, format, ...) ----
+    va_list args;
+    va_start(args, format);
+    vsnprintf(errorMessage, 1024, format, args);
+    va_end(args);
+    
+    mexErrMsgIdAndTxt(fullId, errorMessage);
+}
\ No newline at end of file
diff --git a/tools/graph_tools/triplepointSegmentsExampleData.mat b/tools/graph_tools/triplepointSegmentsExampleData.mat
new file mode 100644
index 000000000..a25e4489d
Binary files /dev/null and b/tools/graph_tools/triplepointSegmentsExampleData.mat differ
diff --git a/tools/graph_tools/triplepointSegmentsExampleScript.m b/tools/graph_tools/triplepointSegmentsExampleScript.m
new file mode 100644
index 000000000..0430cecac
--- /dev/null
+++ b/tools/graph_tools/triplepointSegmentsExampleScript.m
@@ -0,0 +1,18 @@
+load triplepointSegmentsExampleData.mat
+
+[s, sP, sf, sfF, sfD] = TriplepointSegments(F,V);
+
+clf;
+hold on;
+for i=1:(size(s,1)-1)
+    current_sP = sP(s(i):(s(i+1)-1));
+    Vsection = V(current_sP,:);
+    plot(Vsection(:,1),Vsection(:,2));
+end
+
+uiwait;
+
+hold on;
+i = 1:(size(s,1)-1);
+plot([V(sP(s(i)),1)';V(sP(s(i+1)-1),1)'],[V(sP(s(i)),2)';V(sP(s(i+1)-1),2)']);
+scatter([V(sP(s(i)),1)';V(sP(s(i+1)-1),1)'],[V(sP(s(i)),2)';V(sP(s(i+1)-1),2)']);