Index: trunk/psLib/src/math/psMixtureModels.c
===================================================================
--- trunk/psLib/src/math/psMixtureModels.c	(revision 35767)
+++ 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);
Index: trunk/psLib/src/math/psMixtureModels.h
===================================================================
--- trunk/psLib/src/math/psMixtureModels.h	(revision 35767)
+++ trunk/psLib/src/math/psMixtureModels.h	(revision 35785)
@@ -2,4 +2,8 @@
 #define PS_MIXTUREMODELS_H
 
+// These are the recommended functions to call.  They in turn call code
+// that scans the input data matrix/image and censored data points with
+// non-finite values.  The censored data is then passed to the next set
+// of functions.
 
 bool psMMkmeans(psImage *D,
@@ -24,4 +28,34 @@
 	     double *V);
 
+// These functions assume that the input data matrix/image is cleaned of
+// non-finite values.  They should work identically as the ones above on
+// data that is well defined.
+bool psMMkmeansUncensored(psImage *D,
+			  int dim,
+			  long N,
+			  psVector *modes,
+			  psImage *means,
+			  int m,
+			  int *iterations,
+			  double *V);
+
+bool psMMGMMUncensored(psImage *D,
+		       int dim,
+		       long N,
+		       psVector *modes,
+		       psImage *means,
+		       psArray *sigma,
+		       psVector *pi,
+		       psImage *P,
+		       int m,
+		       int *iterations,
+		       double *V);
+
+// This function does the classification in five stages:
+//  1) pass data to K-means code with m modes for initial values.
+//  2) pass data+K-mean values to GMM code with m modes for best values.
+//  3) repeat K-means pass with m-1 modes.
+//  4) repeat GMM pass with m-1 modes.
+//  5) calculate Pvalue that the model with m modes is better than that with m-1
 bool psMMClass(psImage *D,
 	       int dim,
@@ -34,5 +68,11 @@
 	       int m,
 	       double *Pval);
-  
+
+
+// These functions are one dimensional equivalent wrappers of the above functions.
+// The psVector data is inserted into the appropriate sized psImage structure to
+// avoid having to do so in each program calling this library.  Since this is the
+// only special case (2-dim is a matrix, n-dim is matrix as well), no other wrappers
+// are needed.
 bool psMM1Dkmeans(psVector *D,
 		  long N,