Index: trunk/psLib/src/math/psMixtureModels.c
===================================================================
--- trunk/psLib/src/math/psMixtureModels.c	(revision 35469)
+++ trunk/psLib/src/math/psMixtureModels.c	(revision 35785)
@@ -112,4 +112,38 @@
 }
 
+psImage *psMMCensorData(psImage *in,
+			long *Ncensored,
+			psVector *offsets) {
+  // Allocate offset vector
+  offsets = psVectorRecycle(offsets,in->numRows,PS_TYPE_S32);
+
+  *Ncensored = 0;
+
+  for (int i = 0; i < in->numRows; i++) {
+    int isCensored = 0;
+    for (int j = 0; j < in->numCols; j++) {
+      if (!isfinite(in->data.F32[i][j])) {
+	isCensored = 1;
+      }
+    }
+    if (isCensored) {
+      *Ncensored = *Ncensored + 1;
+    }
+    offsets->data.S32[i] = *Ncensored;
+  }
+  psImage *out;
+  if (*Ncensored == 0) {
+    out = psMemIncrRefCounter(in);
+  }
+  else {
+    out = psImageAlloc(in->numCols,in->numRows - *Ncensored,PS_TYPE_F32);
+    for (int i = 0; i < out->numRows; i++) {
+      for (int j = 0; j < out->numCols; j++) {
+	out->data.F32[i][j] = in->data.F32[i + offsets->data.S32[i]][j];
+      }
+    }
+  }
+  return(out);
+}
 
 /////////////////////////////////////////////////////////////////////////////////
@@ -117,17 +151,62 @@
 /////////////////////////////////////////////////////////////////////////////////
 
-
 #define MAX_ITERATIONS 100
 #define THRESHOLD      1e-3
 
-
-bool psMMkmeans(psImage *D,
-		int dim,
-		long N,
-		psVector *modes,
-		psImage *means,
-		int m,
-		int *iterations,
-		double *V) {
+bool psMMkmeans(psImage *D,      //< Input data to fit.  					 
+		int dim,	 //< dimension of the data					 
+		long N,		 //< number of data points					 
+		psVector *modes, //< Output vector of which mode each data point belongs to	 
+		psImage *means,	 //< Matrix containing mode mean values for each dimension	 
+		int m,		 //< Number of modes to fit					 
+		int *iterations, //< Number of iterations to perform				 
+		double *V) {	 //< Deviation value.  Pseudo chi^2                           
+  PS_ASSERT_IMAGE_NON_NULL(D,false);
+  PS_ASSERT_IMAGE_TYPE(D,PS_TYPE_F32,false);
+  // We need to scan the input for invalid d
+  long Ncensored;
+  psVector *offsets = psVectorAlloc(N,PS_TYPE_S32);
+  psImage  *Dcensored = psMMCensorData(D,&Ncensored,offsets);
+
+  psVector *tempModes = psVectorAlloc(N - Ncensored,PS_TYPE_F32);
+  bool success = psMMkmeansUncensored(Dcensored,
+				      dim,
+				      N - Ncensored,
+				      tempModes,
+				      means,
+				      m,
+				      iterations,
+				      V);
+  if (Ncensored == 0) { // We can directly return the modes and means calculated
+    if (modes) {
+      psFree(modes);
+    }
+    modes = psMemIncrRefCounter(tempModes);
+  }
+  else {
+    if (!modes) {
+      modes = psVectorAlloc(N,PS_TYPE_F32);
+    }
+
+    psVectorInit(modes,NAN);
+    for (int i = 0; i < N - Ncensored; i++) {
+      modes->data.F32[i - offsets->data.S32[i]] = tempModes->data.F32[i];
+    }
+  }
+  psFree(offsets);
+  psFree(Dcensored);
+  psFree(tempModes);
+  return(success);
+}
+
+
+bool psMMkmeansUncensored(psImage *D,       //< Input data to fit.  
+			  int dim,          //< dimension of the data
+			  long N,           //< number of data points
+			  psVector *modes,  //< Output vector of which mode each data point belongs to
+			  psImage *means,   //< Matrix containing mode mean values for each dimension
+			  int m,            //< Number of modes to fit
+			  int *iterations,  //< Number of iterations to perform
+			  double *V) {      //< Deviation value.  Pseudo chi^2                           
   // Assertions
   PS_ASSERT_IMAGE_NON_NULL(D,false);
@@ -135,7 +214,8 @@
   //  PS_ASSERT_IMAGE_SIZE(D,N,dim,false);
 
+  
   // Allocate things we can allocate here if they're not already allocated
   if (!modes) {
-    modes = psVectorAlloc(m,PS_TYPE_F32);
+    modes = psVectorAlloc(N,PS_TYPE_F32);
   }
   if (!means) {
@@ -277,15 +357,82 @@
 }
 
-bool psMMGMM(psImage *D,
-	     int dim,
-	     long N,
-	     psVector *modes,
-	     psImage *means,
-	     psArray *sigma,
-	     psVector *pi,
-	     psImage *P,
-	     int m,
-	     int *iterations,
-	     double *V) {
+bool psMMGMM(psImage *D,           //< Input data to fit
+	     int dim,              //< Dimension of each data vector
+	     long N,               //< Number of data vectors
+	     psVector *modes,      //< Output vector assigning each data point to a mode
+	     psImage *means,       //< Output matrix containing vectors of the means for each mode
+	     psArray *sigma,       //< Output array containing covariance matrices for each mode
+	     psVector *pi,         //< Output vector contining population fractions for each mode
+	     psImage *P,           //< Output matrix containing the probabilities a data point is in each mode
+	     int m,                //< Number of modes to fit
+	     int *iterations,      //< Number of iterations performed
+	     double *V) {          //< Log-likelihood of final solution.
+
+  PS_ASSERT_IMAGE_NON_NULL(D,false);
+  PS_ASSERT_IMAGE_TYPE(D,PS_TYPE_F32,false);
+  // We need to scan the input for invalid d
+  long Ncensored;
+  psVector *offsets = psVectorAlloc(N,PS_TYPE_F32);
+  psImage  *Dcensored = psMMCensorData(D,&Ncensored,offsets);
+
+  psVector *tempModes = psVectorAlloc(N - Ncensored,PS_TYPE_F32);
+  psImage  *tempP     = psImageAlloc(m,N - Ncensored,PS_TYPE_F32);
+  bool success = psMMGMMUncensored(Dcensored,
+				   dim,
+				   N - Ncensored,
+				   tempModes,
+				   means,
+				   sigma,
+				   pi,
+				   tempP,
+				   m,
+				   iterations,
+				   V);
+  if (Ncensored == 0) { // We can directly return the modes and means calculated
+    if (modes) {
+      psFree(modes);
+    }
+    if (P) {
+      psFree(P);
+    }
+    modes = psMemIncrRefCounter(tempModes);
+    P     = psMemIncrRefCounter(tempP);
+  }
+  else {
+    if (!modes) {
+      modes = psVectorAlloc(N,PS_TYPE_F32);
+    }
+    if (!P) {
+      P = psImageAlloc(N,m,PS_TYPE_F32);
+    }
+    psVectorInit(modes,NAN);
+    psImageInit(P,NAN);
+
+    for (int i = 0; i < N - Ncensored; i++) {
+      modes->data.F32[i - offsets->data.S32[i]] = tempModes->data.F32[i];
+      for (int j = 0; j < dim; j++) {
+	P->data.F32[i - offsets->data.S32[i]][j] = tempP->data.F32[i][j];
+      }
+    }
+  }
+  psFree(offsets);
+  psFree(Dcensored);
+  psFree(tempModes);
+  psFree(tempP);
+  return(success);
+}
+
+
+bool psMMGMMUncensored(psImage *D,           //< Input data to fit
+		       int dim,              //< Dimension of each data vector
+		       long N,               //< Number of data vectors
+		       psVector *modes,      //< Output vector assigning each data point to a mode
+		       psImage *means,       //< Output matrix containing vectors of the means for each mode
+		       psArray *sigma,       //< Output array containing covariance matrices for each mode
+		       psVector *pi,         //< Output vector contining population fractions for each mode
+		       psImage *P,           //< Output matrix containing the probabilities a data point is in each mode
+		       int m,                //< Number of modes to fit
+		       int *iterations,      //< Number of iterations performed
+		       double *V) {          //< Log-likelihood of final solution.
   // Assertions
   PS_ASSERT_IMAGE_NON_NULL(D,false);