diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp
index 7821a526b0..f9d18e1560 100644
--- a/Libs/Mesh/Mesh.cpp
+++ b/Libs/Mesh/Mesh.cpp
@@ -1393,6 +1393,7 @@ bool Mesh::prepareFFCFields(std::vector<std::vector<Eigen::Vector3d>> boundaries
       auto absvalues = setDistanceToBoundaryValueFieldForFFCs(values, points, boundaryVerts, inout, V, F);
       this->poly_data_->GetPointData()->SetActiveScalars("value");
       std::vector<Eigen::Matrix3d> face_grad = this->setGradientFieldForFFCs(absvalues, V, F);
+
     }
 
   }  // Per boundary for loop end
@@ -1718,4 +1719,17 @@ Mesh& Mesh::operator+=(const Mesh& otherMesh) {
   return *this;
 }
 
+void Mesh::computeFFCGradients(vtkSmartPointer<vtkDoubleArray> inout){
+  // Extract mesh vertices and faces
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  vtkSmartPointer<vtkPoints> points;
+  points = getIGLMesh(V, F);
+  auto values = vtkSmartPointer<vtkDoubleArray>::New();
+  auto absvalues = setDistanceToBoundaryValueFieldForFFCs(values, points, boundaryVerts, inout, V, F);
+  this->poly_data_->GetPointData()->SetActiveScalars("value");
+  std::vector<Eigen::Matrix3d> face_grad = this->setGradientFieldForFFCs(absvalues, V, F);
+}
+
 }  // namespace shapeworks
diff --git a/Libs/Mesh/Mesh.h b/Libs/Mesh/Mesh.h
index 84abcd9517..7752f6161d 100644
--- a/Libs/Mesh/Mesh.h
+++ b/Libs/Mesh/Mesh.h
@@ -238,6 +238,8 @@ class Mesh {
   bool prepareFFCFields(std::vector<std::vector<Eigen::Vector3d>> boundaries, Eigen::Vector3d query,
                         bool onlyGenerateInOut = false);
 
+  void computeFFCGradients(vtkSmartPointer<vtkDoubleArray> inout);
+
   /// Gets values for FFCs
   double getFFCValue(Eigen::Vector3d query) const;
 
diff --git a/Libs/Optimize/Constraints.cpp b/Libs/Optimize/Constraints.cpp
index 07524b311f..237c22ef62 100644
--- a/Libs/Optimize/Constraints.cpp
+++ b/Libs/Optimize/Constraints.cpp
@@ -784,8 +784,8 @@ void Constraints::clipMesh(Mesh& mesh) {
 
   if (getFreeformConstraint().isSet()) {
     auto& ffc = getFreeformConstraint();
-    mesh.prepareFFCFields(ffc.boundaries(), ffc.getQueryPoint(), true);
-    mesh = Mesh(mesh.clipByField("inout", 1.0));
+    ffc.applyToPolyData(mesh.getVTKMesh());
+    mesh = Mesh(mesh.clipByField("ffc_paint", 1.0));
   }
 }
 
diff --git a/Libs/Optimize/FreeFormConstraint.cpp b/Libs/Optimize/FreeFormConstraint.cpp
index 5d9f6a60f6..bb8bf0b921 100644
--- a/Libs/Optimize/FreeFormConstraint.cpp
+++ b/Libs/Optimize/FreeFormConstraint.cpp
@@ -3,6 +3,7 @@
 #include <vtkContourFilter.h>
 #include <vtkContourLoopExtraction.h>
 #include <vtkFloatArray.h>
+#include <vtkKdTreePointLocator.h>
 #include <vtkPointData.h>
 
 namespace shapeworks {
@@ -30,7 +31,7 @@ void FreeFormConstraint::setDefinition(vtkSmartPointer<vtkPolyData> polyData) {
 void FreeFormConstraint::applyToPolyData(vtkSmartPointer<vtkPolyData> polyData) {
   Mesh mesh(polyData);
 
-  if (boundaries_.empty()) {
+  if (!inoutPolyData_) {
     auto array = createFFCPaint(polyData);
     if (!painted_) {
       array->FillComponent(0, 1.0);
@@ -38,15 +39,22 @@ void FreeFormConstraint::applyToPolyData(vtkSmartPointer<vtkPolyData> polyData)
     return;
   }
 
-  mesh.prepareFFCFields(boundaries_, queryPoint_, true);
-  auto inout = mesh.getField("inout", Mesh::FieldType::Point);
+  vtkFloatArray* inout = vtkFloatArray::SafeDownCast(inoutPolyData_->GetPointData()->GetArray("ffc_paint"));
 
   auto array = createFFCPaint(polyData);
   array->FillComponent(0, 1.0);
+
+  // create a kd tree point locator
+  vtkSmartPointer<vtkKdTreePointLocator> kdTree = vtkSmartPointer<vtkKdTreePointLocator>::New();
+  kdTree->SetDataSet(inoutPolyData_);
+  kdTree->BuildLocator();
+
+  // iterate over all points in polyData and find the closest point in the kd tree
   for (int i = 0; i < polyData->GetNumberOfPoints(); i++) {
-    double* value = inout->GetTuple(i);
-    float f = value[0];
-    array->SetTuple(i, &f);
+    double p[3];
+    polyData->GetPoint(i, p);
+    vtkIdType id = kdTree->FindClosestPoint(p);
+    array->SetTuple(i, inout->GetTuple(id));
   }
 }
 
@@ -151,6 +159,42 @@ void FreeFormConstraint::createInoutPolyData() {
   inoutPolyData_ = vtkSmartPointer<vtkPolyData>::New();
   inoutPolyData_->DeepCopy(definitionPolyData_);
 }
+
+//-----------------------------------------------------------------------------
+void FreeFormConstraint::computeGradientFields(std::shared_ptr<Mesh> mesh) {
+  applyToPolyData(mesh->getVTKMesh());
+  setDefinition(mesh->getVTKMesh());
+  computeBoundaries();
+
+
+  std::vector<size_t> boundaryVerts;
+  for (auto& boundary : boundaries_) {
+    for (auto& pos : boundary) {
+      boundaryVerts.push_back(pos);
+    }
+  }
+
+  // Extract mesh vertices and faces
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  vtkSmartPointer<vtkPoints> points;
+  points = mesh->getIGLMesh(V, F);
+  auto values = vtkSmartPointer<vtkDoubleArray>::New();
+  auto absvalues = setDistanceToBoundaryValueFieldForFFCs(values, points, boundaryVerts, getInOutScalars(), V, F);
+  mesh->getVTKMesh()->GetPointData()->SetActiveScalars("value");
+  std::vector<Eigen::Matrix3d> face_grad = mesh->setGradientFieldForFFCs(absvalues, V, F);
+}
+
+//-----------------------------------------------------------------------------
+vtkFloatArray* FreeFormConstraint::getInOutScalars() {
+  if (!inoutPolyData_) {
+    return nullptr;
+  }
+  vtkFloatArray* array = vtkFloatArray::SafeDownCast(inoutPolyData_->GetPointData()->GetArray("ffc_paint"));
+  return array;
+}
+
 //-----------------------------------------------------------------------------
 
 }  // namespace shapeworks
diff --git a/Libs/Optimize/FreeFormConstraint.h b/Libs/Optimize/FreeFormConstraint.h
index 08083d3777..ff5b80c225 100644
--- a/Libs/Optimize/FreeFormConstraint.h
+++ b/Libs/Optimize/FreeFormConstraint.h
@@ -69,7 +69,11 @@ class FreeFormConstraint : public Constraint {
   //! Reset to initial state
   void reset();
 
+  void computeGradientFields(std::shared_ptr<Mesh> mesh);
+
  private:
+
+  vtkFloatArray* getInOutScalars();
   vtkFloatArray* createFFCPaint(vtkSmartPointer<vtkPolyData> polyData);
 
   std::shared_ptr<shapeworks::Mesh> mesh_;
diff --git a/Libs/Optimize/OptimizeParameterFile.cpp b/Libs/Optimize/OptimizeParameterFile.cpp
index 4bb6f5ad92..25bf9613b0 100644
--- a/Libs/Optimize/OptimizeParameterFile.cpp
+++ b/Libs/Optimize/OptimizeParameterFile.cpp
@@ -540,8 +540,8 @@ bool OptimizeParameterFile::read_mesh_inputs(TiXmlHandle* docHandle, Optimize* o
             std::cout << "ffcssize " << ffcs.size() << std::endl;
           }
             if (index < ffcs.size()) {
-              mesh.prepareFFCFields(ffcs[index].boundaries, ffcs[index].query);
-              mesh = Mesh(mesh.clipByField("inout", 0.0));
+              ffcs[index].applyToPolyData(mesh.getVTKMesh());
+              mesh = Mesh(mesh.clipByField("ffc_paint", 0.0));
             }
         }
 
@@ -822,7 +822,7 @@ bool OptimizeParameterFile::read_constraints(TiXmlHandle* doc_handle, Optimize*
     return false;
   }
 
-  return this->read_cutting_ffcs(doc_handle, optimize);
+  return true;
 }
 
 //---------------------------------------------------------------------------
@@ -1162,116 +1162,6 @@ bool OptimizeParameterFile::read_cutting_spheres(TiXmlHandle* doc_handle, Optimi
   return true;
 }
 
-bool OptimizeParameterFile::read_cutting_ffcs(TiXmlHandle* docHandle, Optimize* optimize)
-{
-  TiXmlElement* elem;
-  std::istringstream inputsBuffer;
-  int num_inputs = this->get_num_inputs(docHandle);
-
-  int num_ffcs = num_inputs;
-
-  if (this->verbosity_level_ > 1) {
-    std::cout << "Number of free-form constraints " << num_ffcs << std::endl;
-  }
-  // Check that if ffcs are not given, return true. If we expected ffcs(num_ffcs_per_input was given), print a warning to cerr.
-  elem = docHandle->FirstChild("ffcs").Element();
-  if (!elem) {
-    return true;
-  }
-
-  inputsBuffer.str(elem->GetText());
-  auto flags = optimize->GetDomainFlags();
-
-  // load ffc filenames
-  std::vector<std::string> ffcFiles;
-  std::string ffcfilename;
-  while (inputsBuffer >> ffcfilename) {
-    ffcFiles.push_back(ffcfilename);
-  }
-
-  size_t count = 0;
-  std::string buffer;
-  if (ffcFiles.size() != num_ffcs) {
-    std::cerr << "ERROR: Error reading free-form constraint(ffc) data! Expected " << num_ffcs
-              << " but instead got " << ffcFiles.size()
-              << " instead. Can't use ffcs, please fix the input xml." << std::endl;
-    throw 1;
-  }
-  else {
-    for (size_t shapeCount = 0; shapeCount < num_inputs; shapeCount++) {
-      // Reading in ffc file for shape shapeCount # ffcCount
-      std::string fn = ffcFiles[count];
-      //std::cout << "Shape " << shapeCount << " ffc num " << ffcCount << " filename " <<  fn << std::endl;
-      bool query_read = true;
-
-      std::vector<std::vector<Eigen::Vector3d>> boundaries;
-      int boundary_count = -1;
-      Eigen::Vector3d query;
-
-      fstream newfile;
-      newfile.open(fn, ios::in);
-
-      if (newfile.is_open()) {   //checking whether the file is open
-        std::string tp;
-        while (getline(newfile, tp)) { //read data from file object and put it into string.
-          //cout << tp << "\n"; //print the data of the string
-          if (tp == "query") {
-            query_read = true;
-          }
-          else if (tp == "boundary_pts") {
-            query_read = false;
-            std::vector<Eigen::Vector3d> boundary;
-            boundaries.push_back(boundary);
-            boundary_count++;
-          }
-          else {
-            try {
-              if (query_read) {
-                std::istringstream iss(tp);
-                iss >> buffer;
-                query(0) = std::stod(buffer);
-                iss >> buffer;
-                query(1) = std::stod(buffer);
-                iss >> buffer;
-                query(2) = std::stod(buffer);
-                //cout << "Added query " << query.transpose() << "\n";
-              }
-              else {
-                Eigen::Vector3d bpt;
-                std::istringstream iss(tp);
-                iss >> buffer;
-                bpt(0) = std::stod(buffer);
-                iss >> buffer;
-                bpt(1) = std::stod(buffer);
-                iss >> buffer;
-                bpt(2) = std::stod(buffer);
-                boundaries[boundary_count].push_back(bpt);
-                //cout << "Added boundary " << bpt.transpose() << "\n";
-              }
-            }
-            catch (std::exception& e) {
-              cout << e.what() << '\n';
-              std::cerr << "ERROR: File " << fn
-                        << " threw an exception. Please check the correctness and formating of the file."
-                        << std::endl;
-              throw 1;
-            }
-          }
-        }
-        newfile.close(); //close the file object.
-      }
-      else {
-        std::cerr << "ERROR: File " << fn
-                  << " could not be opened. Please check that the file is available." << std::endl;
-        throw 1;
-      }
-      optimize->GetSampler()->AddFreeFormConstraint(shapeCount, boundaries, query);
-      count++;
-    }
-  }
-
-  return true;
-}
 
 //---------------------------------------------------------------------------
 bool OptimizeParameterFile::read_explanatory_variables(TiXmlHandle* doc_handle, Optimize* optimize)
diff --git a/Libs/Optimize/OptimizeParameterFile.h b/Libs/Optimize/OptimizeParameterFile.h
index 4c7e180e8c..1affd0a52a 100644
--- a/Libs/Optimize/OptimizeParameterFile.h
+++ b/Libs/Optimize/OptimizeParameterFile.h
@@ -59,8 +59,6 @@ class OptimizeParameterFile {
 
   bool read_cutting_spheres(TiXmlHandle* doc_handle, Optimize* optimize);
 
-  bool read_cutting_ffcs(TiXmlHandle* docHandle, Optimize* optimize);
-
   bool read_explanatory_variables(TiXmlHandle* doc_handle, Optimize* optimize);
 
   bool read_flag_particles(TiXmlHandle* doc_handle, Optimize* optimize);
diff --git a/Libs/Optimize/OptimizeParameters.cpp b/Libs/Optimize/OptimizeParameters.cpp
index 1c712ce4b8..d3f01fde8a 100644
--- a/Libs/Optimize/OptimizeParameters.cpp
+++ b/Libs/Optimize/OptimizeParameters.cpp
@@ -391,7 +391,7 @@ bool OptimizeParameters::set_up_optimize(Optimize *optimize) {
 
         auto& ffc = constraint.getFreeformConstraint();
         if (ffc.isSet()) {
-          optimize->GetSampler()->AddFreeFormConstraint(domain_count, ffc.boundaries(), ffc.getQueryPoint());
+          optimize->GetSampler()->AddFreeFormConstraint(domain_count, ffc);
         }
       }
 
diff --git a/Libs/Optimize/Sampler.cpp b/Libs/Optimize/Sampler.cpp
index 2f0e93e122..7fdde44b19 100644
--- a/Libs/Optimize/Sampler.cpp
+++ b/Libs/Optimize/Sampler.cpp
@@ -383,14 +383,11 @@ void Sampler::AddSphere(unsigned int i, vnl_vector_fixed<double, Dimension>& c,
   }
 }
 
-void Sampler::AddFreeFormConstraint(unsigned int i, const std::vector<std::vector<Eigen::Vector3d> > boundaries,
-                                    const Eigen::Vector3d query) {
-  if (m_FFCs.size() < i + 1) {
-    m_FFCs.resize(i + 1);
+void Sampler::AddFreeFormConstraint(int domain, const FreeFormConstraint &ffc) {
+  if (m_FFCs.size() < domain + 1) {
+    m_FFCs.resize(domain + 1);
   }
-
-  m_FFCs[i].boundaries = boundaries;
-  m_FFCs[i].query = query;
+  m_FFCs[domain] = ffc;
 }
 
 void Sampler::AddImage(ImageType::Pointer image, double narrow_band, std::string name) {
@@ -419,15 +416,14 @@ bool Sampler::initialize_ffcs(size_t dom) {
   auto mesh = std::make_shared<Mesh>(m_meshes[dom]);
   if (m_verbosity >= 1) std::cout << "dom " << dom << " point count " << mesh->numPoints() << " faces " << mesh->numFaces() << std::endl;
 
-  if (m_FFCs[dom].boundaries.size() > 0) {
-    if (m_verbosity >= 1) {
-      std::cout << "Splitting mesh FFC for domain " << dom << " shape with query point "
-                << m_FFCs[dom].query.transpose() << std::endl;
-    }
-    mesh->prepareFFCFields(m_FFCs[dom].boundaries, m_FFCs[dom].query);
+  if (m_FFCs[dom].isSet()) {
     this->m_DomainList[dom]->GetConstraints()->addFreeFormConstraint(mesh);
+    m_FFCs[dom].
   }
 
+  // todo, initialize ffc gradient fields
+
+
 
 #if defined(VIZFFC)
   MeshUtils mutil;
diff --git a/Libs/Optimize/Sampler.h b/Libs/Optimize/Sampler.h
index 07f0872cad..f0c293c6e0 100644
--- a/Libs/Optimize/Sampler.h
+++ b/Libs/Optimize/Sampler.h
@@ -68,12 +68,6 @@ class Sampler {
     double radius;
   };
 
-  /** Convenient typedef for storing free form constraint information */
-  struct FFCType {
-     std::vector< std::vector< Eigen::Vector3d > > boundaries;
-     Eigen::Vector3d query;
-  };
-
   /** Returns the particle system used in the surface sampling. */
   itkGetObjectMacro(ParticleSystem, itk::ParticleSystem);
 
@@ -177,9 +171,7 @@ class Sampler {
                        const vnl_vector_fixed<double, Dimension>& va,
                        const vnl_vector_fixed<double, Dimension>& vb,
                        const vnl_vector_fixed<double, Dimension>& vc);
-  void AddFreeFormConstraint(unsigned int i,
-                             const std::vector< std::vector< Eigen::Vector3d > > boundaries,
-                             const Eigen::Vector3d query);
+  void AddFreeFormConstraint(int domain, const FreeFormConstraint &ffc);
 
   /** Transform a cutting plane based on procrustes transformation */
   void TransformCuttingPlanes(unsigned int i);
@@ -452,7 +444,7 @@ class Sampler {
 
   }
 
-  std::vector<FFCType> GetFFCs() { return m_FFCs; }
+  std::vector<FreeFormConstraint> GetFFCs() { return m_FFCs; }
 
   void SetMeshFFCMode(bool mesh_ffc_mode) {m_meshFFCMode = mesh_ffc_mode;}
 
@@ -548,7 +540,7 @@ class Sampler {
   std::string m_PrefixTransformFile;
   std::vector<std::vector<CuttingPlaneType>> m_CuttingPlanes;
   std::vector<std::vector<SphereType>> m_Spheres;
-  std::vector<FFCType> m_FFCs;
+  std::vector<FreeFormConstraint> m_FFCs;
   std::vector<vtkSmartPointer<vtkPolyData>> m_meshes;
   bool m_meshFFCMode = false;