diff --git a/common/Domain.cpp b/common/Domain.cpp
index 0fa8545a..a88b47da 100644
--- a/common/Domain.cpp
+++ b/common/Domain.cpp
@@ -329,7 +329,6 @@ void Domain::Decomp(std::string Filename)
 	char *SegData = NULL;
 	
 	if (RANK==0){
-
 		printf("Input media: %s\n",Filename.c_str());
 		printf("Relabeling %lu values\n",ReadValues.size());
 		for (int idx=0; idx<ReadValues.size(); idx++){
@@ -485,6 +484,7 @@ void Domain::Decomp(std::string Filename)
 						else{
 							// solid checkers
 							SegData[k*global_Nx*global_Ny+j*global_Nx+i] = 0;
+							
 						}
 					}
 				}
diff --git a/common/Domain.h b/common/Domain.h
index 455a8012..5d84d1d6 100755
--- a/common/Domain.h
+++ b/common/Domain.h
@@ -16,7 +16,6 @@
 #include "common/Communication.h"
 #include "common/Database.h"
 
-
 class Domain;
 template<class TYPE> class PatchData;
 
diff --git a/common/ReadMicroCT.cpp b/common/ReadMicroCT.cpp
index f02163dd..94c728ef 100644
--- a/common/ReadMicroCT.cpp
+++ b/common/ReadMicroCT.cpp
@@ -70,6 +70,8 @@ Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
     auto n = domain.getVector<int>( "n" );
     int rank = comm_rank(MPI_COMM_WORLD);
     auto nproc = domain.getVector<int>( "nproc" );
+	auto ReadValues = domain.getVector<int>( "ReadValues" );
+	auto WriteValues = domain.getVector<int>( "WriteValues" );
     RankInfoStruct rankInfo( rank, nproc[0], nproc[1], nproc[2] );
     
     // Determine the largest file number to get
@@ -81,18 +83,37 @@ Array<uint8_t> readMicroCT( const Database& domain, MPI_Comm comm )
     Array<uint8_t> data;
     RankInfoStruct srcRankInfo( rank, Nfx, Nfy, Nfz );
     if ( srcRankInfo.ix >= 0 ) {
-        auto filename = domain.getScalar<std::string>( "Filename" );
-        char tmp[100];
-        if ( filename.find( "0x_0y_0z.gbd.gz" ) != std::string::npos ) {
-            sprintf( tmp, "%ix_%iy_%iz.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
-            filename = filename.replace( filename.find( "0x_0y_0z.gbd.gz" ), 15, std::string( tmp ) );
-        } else if ( filename.find( "x0_y0_z0.gbd.gz" ) != std::string::npos ) {
-            sprintf( tmp, "x%i_y%i_z%i.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
-            filename = filename.replace( filename.find( "x0_y0_z0.gbd.gz" ), 15, std::string( tmp ) );
-        } else {
-            ERROR( "Invalid name for first file" );
-        }
-        data = readMicroCT( filename );
+    	auto filename = domain.getScalar<std::string>( "Filename" );
+    	char tmp[100];
+    	if ( filename.find( "0x_0y_0z.gbd.gz" ) != std::string::npos ) {
+    		sprintf( tmp, "%ix_%iy_%iz.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
+    		filename = filename.replace( filename.find( "0x_0y_0z.gbd.gz" ), 15, std::string( tmp ) );
+    	} else if ( filename.find( "x0_y0_z0.gbd.gz" ) != std::string::npos ) {
+    		sprintf( tmp, "x%i_y%i_z%i.gbd.gz", srcRankInfo.ix, srcRankInfo.jy, srcRankInfo.kz );
+    		filename = filename.replace( filename.find( "x0_y0_z0.gbd.gz" ), 15, std::string( tmp ) );
+    	} else {
+    		ERROR( "Invalid name for first file" );
+    	}
+    	data = readMicroCT( filename );
+
+    	// Relabel the data
+    	for (int k = 0; k<1024; k++){
+    		for (int j = 0; j<1024; j++){
+    			for (int i = 0; i<1024; i++){
+    				//n = k*Nfx*Nfy + j*Nfx + i;
+    				//char locval = loc_id[n];
+    				char locval = data(i,j,k);
+    				for (int idx=0; idx<ReadValues.size(); idx++){
+    					signed char oldvalue=ReadValues[idx];
+    					signed char newvalue=WriteValues[idx];
+    					if (locval == oldvalue){
+    						data(i,j,k) = newvalue;
+    						idx = ReadValues.size();
+    					}
+    				}
+    			}
+    		}
+    	}
     }
 
     // Redistribute the data
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index c6692312..69b5f485 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -4,6 +4,7 @@ color lattice boltzmann model
 #include "models/ColorModel.h"
 #include "analysis/distance.h"
 #include "analysis/morphology.h"
+#include "common/ReadMicroCT.h"
 #include <stdlib.h>
 #include <time.h>
 
@@ -189,6 +190,10 @@ void ScaLBL_ColorModel::ReadInput(){
 		Mask->Decomp(first_image);
 		IMAGE_INDEX++;
 	}
+	else if (domain_db->keyExists( "GridFile" )){
+	    auto input_id = readMicroCT( *domain_db, MPI_COMM_WORLD );
+	    for (int i=0; i<Nx*Ny*Nz; i++) Mask->id[i] = input_id(i);
+	}
 	else if (domain_db->keyExists( "Filename" )){
 		auto Filename = domain_db->getScalar<std::string>( "Filename" );
 		Mask->Decomp(Filename);