--- a
+++ b/combinedDeepLearningActiveContour/minFunc/mcholC.c
@@ -0,0 +1,199 @@
+#include <math.h>
+#include "mex.h"
+
+double mymax(double x, double y)
+{
+    if (x > y)
+        return x;
+    else
+        return y;
+}
+
+double absolute(double x)
+{
+    if (x >= -x)
+        return x;
+    else
+        return -x;
+}
+
+void permuteInt(int *x, int p, int q)
+{
+    int temp;
+    temp = x[p];
+    x[p] = x[q];
+    x[q] = temp;
+}
+
+void permute(double *x, int p, int q)
+{
+    double temp;
+    temp = x[p];
+    x[p] = x[q];
+    x[q] = temp;
+}
+
+void permuteRows(double *x, int p, int q,int n)
+{
+    int i;
+    double temp;
+    for(i = 0; i < n; i++)
+    {
+        temp = x[p+i*n];
+        x[p+i*n] = x[q+i*n];
+        x[q+i*n] = temp;
+    }
+}
+
+void permuteCols(double *x, int p, int q,int n)
+{
+    int i;
+    double temp;
+    for(i = 0; i < n; i++)
+    {
+        temp = x[i+p*n];
+        x[i+p*n] = x[i+q*n];
+        x[i+q*n] = temp;
+    }
+}
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{
+    int n,sizL[2],sizD[2],i,j,q,s,
+    *P;
+    
+    double mu,gamma,xi,delta,beta,maxVal,theta,
+    *c,    *H, *L, *D, *A;
+    
+    /* Input */
+    H = mxGetPr(prhs[0]);
+    if (nrhs == 1)
+    {
+        mu = 1e-12;
+    }
+    else
+    {
+        mu = mxGetScalar(prhs[1]);
+    }
+    
+    /* Compute Sizes */
+    n = mxGetDimensions(prhs[0])[0];
+    
+    /* Form Output */
+    sizL[0] = n;
+    sizL[1] = n;
+    plhs[0] = mxCreateNumericArray(2,sizL,mxDOUBLE_CLASS,mxREAL);
+    L = mxGetPr(plhs[0]);
+    sizD[0] = n;
+    sizD[1] = 1;
+    plhs[1] = mxCreateNumericArray(2,sizD,mxDOUBLE_CLASS,mxREAL);
+    D = mxGetPr(plhs[1]);
+    plhs[2] = mxCreateNumericArray(2,sizD,mxINT32_CLASS,mxREAL);
+    P = (int*)mxGetData(plhs[2]);
+    
+    /* Initialize */
+    c = mxCalloc(n*n,sizeof(double));
+    A = mxCalloc(n*n,sizeof(double));
+    
+    for (i = 0; i < n; i++)
+    {
+        P[i] = i;
+        for (j = 0;j < n; j++)
+        {
+            A[i+n*j] = H[i+n*j];
+        }
+    }
+    
+    gamma = 0;
+    for (i = 0; i < n; i++)
+    {
+        L[i+n*i] = 1;
+        c[i+n*i] = A[i+n*i];
+    }
+    
+    /* Compute modification parameters */
+    gamma = -1;
+    xi = -1;
+    for (i = 0; i < n; i++)
+    {
+        gamma = mymax(gamma,absolute(A[i+n*i]));
+        for (j = 0;j < n; j++)
+        {
+            //printf("A(%d,%d) = %f, %f\n",i,j,A[i+n*j],absolute(A[i+n*j]));
+            if (i != j)
+                xi = mymax(xi,absolute(A[i+n*j]));
+        }
+    }
+    delta = mu*mymax(gamma+xi,1);
+    
+    if (n > 1)
+    {
+        beta = sqrt(mymax(gamma,mymax(mu,xi/sqrt(n*n-1))));
+    }
+    else
+    {
+        beta = sqrt(mymax(gamma,mu));
+    }
+    
+    for (j = 0; j < n; j++)
+    {
+        
+    /* Find q that results in Best Permutation with j */
+        maxVal = -1;
+        q = 0;
+        for(i = j; i < n; i++)
+        {
+            if (absolute(c[i+n*i]) > maxVal)
+            {
+                maxVal = mymax(maxVal,absolute(c[i+n*i]));
+                q = i;
+            }
+        }
+        
+        /* Permute D,c,L,A,P */
+        permute(D,j,q);
+        permuteInt(P,j,q);
+        permuteRows(c,j,q,n);
+        permuteCols(c,j,q,n);
+        permuteRows(L,j,q,n);
+        permuteCols(L,j,q,n);
+        permuteRows(A,j,q,n);
+        permuteCols(A,j,q,n);
+        
+        for(s = 0; s <= j-1; s++)
+            L[j+n*s] = c[j+n*s]/D[s];
+        
+        for(i = j+1; i < n; i++)
+        {
+            c[i+j*n] = A[i+j*n];
+            for(s = 0; s <= j-1; s++)
+            {
+                c[i+j*n] -= L[j+n*s]*c[i+n*s];
+            }
+        }
+        
+        theta = 0;
+        if (j < n-1)
+        {
+            for(i = j+1;i < n; i++)
+                theta = mymax(theta,absolute(c[i+n*j]));
+        }
+        
+        D[j] = mymax(absolute(c[j+n*j]),mymax(delta,theta*theta/(beta*beta)));
+        
+        if (j < n-1)
+        {
+            for(i = j+1; i < n; i++)
+            {
+                c[i+n*i] = c[i+n*i] - c[i+n*j]*c[i+n*j]/D[j];
+            }
+        }
+        
+    }
+    
+    for(i = 0; i < n; i++)
+        P[i]++;
+    
+    mxFree(c);
+    mxFree(A);
+}
\ No newline at end of file