diff --git a/src/libcode/vx_shapedata/engine.cc b/src/libcode/vx_shapedata/engine.cc
index 9fefb4c738..766aa7e324 100644
--- a/src/libcode/vx_shapedata/engine.cc
+++ b/src/libcode/vx_shapedata/engine.cc
@@ -6,7 +6,6 @@
 // ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
 // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 
-
 ///////////////////////////////////////////////////////////////////////
 
 using namespace std;
@@ -17,6 +16,8 @@ using namespace std;
 #include <stdlib.h>
 #include <string.h>
 #include <cmath>
+#include <set>
+#include <map>
 
 #include "engine.h"
 #include "mode_columns.h"
@@ -591,12 +592,13 @@ void ModeFuzzyEngine::do_fcst_convolution() {
    // Apply a circular convolution to the raw field
    //
 
-   mlog << Debug(4) << " Before fcst convolution: " << fcst_conv->sdebug_examine() << "\n";
-      
+   mlog << Debug(4) << " Before fcst convolution: "
+        << fcst_conv->sdebug_examine() << "\n";
 
    if(r > 0) fcst_conv->conv_filter_circ(2*r + 1, conf_info.Fcst->vld_thresh);
 
-   mlog << Debug(4) << " After fcst convolution: " << fcst_conv->sdebug_examine() << "\n";
+   mlog << Debug(4) << " After fcst convolution: "
+        << fcst_conv->sdebug_examine() << "\n";
 
    need_fcst_conv       = false;
    need_fcst_thresh     = true;
@@ -625,14 +627,16 @@ void ModeFuzzyEngine::do_obs_convolution() {
    mlog << Debug(3) << "Applying circular convolution of radius "
         << r << " to the observation field.\n";
 
-   mlog << Debug(4) << " Before obs convolution: " << obs_conv->sdebug_examine() << "\n";
+   mlog << Debug(4) << " Before obs convolution: "
+        << obs_conv->sdebug_examine() << "\n";
 
    //
    // Apply a circular convolution to the raw field
    //
    if(r > 0) obs_conv->conv_filter_circ(2*r + 1, conf_info.Obs->vld_thresh);
 
-   mlog << Debug(4) << " After obs convolution: " << obs_conv->sdebug_examine() << "\n";
+   mlog << Debug(4) << " After obs convolution: "
+        << obs_conv->sdebug_examine() << "\n";
 
    need_obs_conv       = false;
    need_obs_thresh     = true;
@@ -657,13 +661,17 @@ void ModeFuzzyEngine::do_fcst_thresholding() {
 
    *fcst_thresh = *fcst_raw;
    fcst_thresh->threshold(conf_info.Fcst->conv_thresh);
-   mlog << Debug(4) << " After thresholding raw fcst field: " << fcst_thresh->sdebug_examine() << "\n";
+
+   mlog << Debug(4) << " After thresholding raw fcst field: "
+        << fcst_thresh->sdebug_examine() << "\n";
 
    //
    // Threshold the convolved field
    //
    fcst_mask->threshold(conf_info.Fcst->conv_thresh);
-   mlog << Debug(4) << " After thresholding convolved fcst field: " << fcst_mask->sdebug_examine() << "\n";
+
+   mlog << Debug(4) << " After thresholding convolved fcst field: "
+        << fcst_mask->sdebug_examine() << "\n";
 
    mlog << Debug(3) << "Applying convolution threshold "
         << conf_info.Fcst->conv_thresh.get_str()
@@ -692,13 +700,17 @@ void ModeFuzzyEngine::do_obs_thresholding() {
 
    *obs_thresh = *obs_raw;
    obs_thresh->threshold(conf_info.Obs->conv_thresh);
-   mlog << Debug(4) << " After thresholding raw obs field: " << obs_thresh->sdebug_examine() << "\n";
+
+   mlog << Debug(4) << " After thresholding raw obs field: "
+        << obs_thresh->sdebug_examine() << "\n";
 
    //
    // Threshold the convolved field
    //
    obs_mask->threshold(conf_info.Obs->conv_thresh);
-   mlog << Debug(4) << " After thresholding convolved obs field: " << obs_mask->sdebug_examine() << "\n";
+
+   mlog << Debug(4) << " After thresholding convolved obs field: "
+        << obs_mask->sdebug_examine() << "\n";
 
    mlog << Debug(3) << "Applying convolution threshold "
         << conf_info.Obs->conv_thresh.get_str()
@@ -733,7 +745,9 @@ void ModeFuzzyEngine::do_fcst_filtering() {
                                 fcst_raw, conf_info.Fcst->conv_thresh, grid,
                                 conf_info.Fcst->var_info->is_precipitation());
 
-      mlog << Debug(4) << " After attribute filtering of fcst field: " << fcst_mask->sdebug_examine() << "\n";
+      mlog << Debug(4) << " After attribute filtering of fcst field: "
+           << fcst_mask->sdebug_examine() << "\n";
+
       mlog << Debug(3) << "Applying object attribute filtering"
            << " resulted in " <<  fcst_mask->n_objects()
            << " simple forecast objects.\n";
@@ -769,7 +783,9 @@ void ModeFuzzyEngine::do_obs_filtering() {
                                obs_raw, conf_info.Obs->conv_thresh, grid,
                                conf_info.Obs->var_info->is_precipitation());
 
-      mlog << Debug(4) << " After attribute filtering of obs field: " << obs_mask->sdebug_examine() << "\n";
+      mlog << Debug(4) << " After attribute filtering of obs field: "
+           << obs_mask->sdebug_examine() << "\n";
+
       mlog << Debug(3) << "Applying object attribute filtering"
            << " resulted in " <<  obs_mask->n_objects()
            << " simple observation objects.\n";
@@ -797,7 +813,9 @@ void ModeFuzzyEngine::do_fcst_splitting() {
    if(!need_fcst_split) return;
 
    *fcst_split = split(*fcst_mask, n_fcst);
-   mlog << Debug(4) << " After splitting of fcst field: " << fcst_split->sdebug_examine() << "\n";
+
+   mlog << Debug(4) << " After splitting of fcst field: "
+        << fcst_split->sdebug_examine() << "\n";
 
    need_fcst_split      = false;
    need_fcst_merge      = true;
@@ -816,7 +834,9 @@ void ModeFuzzyEngine::do_obs_splitting() {
    if(!need_obs_split) return;
 
    *obs_split = split(*obs_mask, n_obs);
-   mlog << Debug(4) << " After splitting of obs field: " << obs_split->sdebug_examine() << "\n";
+
+   mlog << Debug(4) << " After splitting of obs field: "
+        << obs_split->sdebug_examine() << "\n";
 
    need_obs_split      = false;
    need_obs_merge      = true;
@@ -1268,15 +1288,11 @@ void ModeFuzzyEngine::do_match_merge() {
 ///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_fcst_merge_thresh() {
-   int j, k, x, y;
-   int n_fcst_merge, intersection;
-   int count, first_k;
+   const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_fcst_merge;
    ShapeData fcst_merge_mask, fcst_merge_split;
 
-   // TODO: reimplement this function to avoid allocating these ShapeData objects
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * fcst_merge_shape = (ShapeData *) 0;
-
    do_fcst_splitting();
 
    //
@@ -1296,90 +1312,60 @@ void ModeFuzzyEngine::do_fcst_merge_thresh() {
    fcst_merge_split = split(fcst_merge_mask, n_fcst_merge);
 
    //
-   // Allocate space for all of the simple forecast shapes and
-   // forecast merge shapes
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   fcst_shape       = new ShapeData [n_fcst];
-   fcst_merge_shape = new ShapeData [n_fcst_merge];
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-   if(!fcst_shape || !fcst_merge_shape) {
+   // Add map entries for each object
+   for(j=1; j<=n_fcst;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_fcst_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-      mlog << Error << "\nModeFuzzyEngine::do_fcst_merge_thresh() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
+   // Loop over all the data points
+   for(j=0; j<fcst_split->data.nxy(); j++) {
+
+      oid = fcst_split->data.buf()[j];
+      mid = fcst_merge_split.data.buf()[j];
+
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
+
+   } // end for j
 
    //
-   // Select all of the simple forecast shapes and
-   // forecast merge shapes
+   // Each simple object should be fully contained in a merge object
    //
-   for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-   }
-
-   for(j=0; j<n_fcst_merge; j++) {
-      fcst_merge_shape[j] = select(fcst_merge_split, j+1);
+   for(oid=1; oid<=n_fcst; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No forecast threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object! "
+              << "Recommend modifying the \"fcst.merge_thresh\" setting.\n\n";
+         return;
+      }
    }
 
    //
    // Calculate the composite object sets
    //
-   for(j=0; j<n_fcst_merge; j++) {
-
-      count = first_k = 0;
-
-      for(k=0; k<n_fcst; k++) {
-
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
-
-         for(x=0; x<fcst_merge_mask.data.nx(); x++) {
-            for(y=0; y<fcst_merge_mask.data.ny(); y++) {
+   for(mid=1; mid<=n_fcst_merge; mid++) {
 
-               if(fcst_merge_shape[j].s_is_on(x, y) &&
-                  fcst_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
-
-         //
-         // Add to set collection only if the fcst object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= fcst area
-         //
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-         if(intersection >= fcst_shape[k].area()) {
-
-            collection.make_room();
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(*it, -1);
+      }
 
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
+      // Increment set counter
+      collection.n_sets++;
 
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(first_k, -1);
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
-
-            count++;
-         }
-      }  // end for k
-
-      if(count > 0) collection.n_sets++;
-   }  // end for j
-
-   //
-   // Done
-   //
-
-   delete [] fcst_shape;       fcst_shape = (ShapeData *) 0;
-   delete [] fcst_merge_shape; fcst_merge_shape = (ShapeData *) 0;
+   } // end for mid
 
    return;
 }
@@ -1393,15 +1379,11 @@ void ModeFuzzyEngine::do_fcst_merge_thresh() {
 ///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_obs_merge_thresh() {
-   int j, k, x, y;
-   int n_obs_merge, intersection;
-   int count, first_k;
+   const char *method_name = "ModeFuzzyEngine::do_obs_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_obs_merge;
    ShapeData obs_merge_mask, obs_merge_split;
 
-   // TODO: reimplement this function to avoid allocating these ShapeData objects
-   ShapeData * obs_shape = (ShapeData *) 0;
-   ShapeData * obs_merge_shape = (ShapeData *) 0;
-
    do_obs_splitting();
 
    //
@@ -1421,89 +1403,60 @@ void ModeFuzzyEngine::do_obs_merge_thresh() {
    obs_merge_split = split(obs_merge_mask, n_obs_merge);
 
    //
-   // Allocate space for all of the simple observation shapes and
-   // observation merge shapes
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   obs_shape       = new ShapeData [n_obs];
-   obs_merge_shape = new ShapeData [n_obs_merge];
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-   if(!obs_shape || !obs_merge_shape) {
+   // Add map entries for each object
+   for(j=1; j<=n_obs;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_obs_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-      mlog << Error << "\nModeFuzzyEngine::do_obs_merge_thresh() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
+   // Loop over all the data points
+   for(j=0; j<obs_split->data.nxy(); j++) {
+
+      oid = obs_split->data.buf()[j];
+      mid = obs_merge_split.data.buf()[j];
+
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
+
+   } // end for j
 
    //
-   // Select all of the simple observation shapes and
-   // simple merge shapes
+   // Each simple object should be fully contained in a merge object
    //
-   for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-   }
-
-   for(j=0; j<n_obs_merge; j++) {
-      obs_merge_shape[j] = select(obs_merge_split, j+1);
+   for(oid=1; oid<=n_obs; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No observation threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object! "
+              << "Recommend modifying the \"obs.merge_thresh\" setting.\n\n";
+         return;
+      }
    }
 
    //
    // Calculate the composite object sets
    //
-   for(j=0; j<n_obs_merge; j++) {
+   for(mid=1; mid<=n_obs_merge; mid++) {
 
-      count = first_k = 0;
-
-      for(k=0; k<n_obs; k++) {
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
-
-         for(x=0; x<obs_merge_mask.data.nx(); x++) {
-            for(y=0; y<obs_merge_mask.data.ny(); y++) {
-
-               if(obs_merge_shape[j].s_is_on(x, y) &&
-                  obs_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
-
-         //
-         // Add to set collection only if the obs object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= obs area
-         //
-         if(intersection >= obs_shape[k].area()) {
-
-            collection.make_room();
-
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
-
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(-1, first_k);
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-
-            count++;
-         }
-      }  // end for k
-
-      if(count > 0) collection.n_sets++;
-   }  // end for j
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(-1, *it);
+      }
 
-   //
-   // Done
-   //
+      // Increment set counter
+      collection.n_sets++;
 
-   delete [] obs_shape;       obs_shape = (ShapeData *) 0;
-   delete [] obs_merge_shape; obs_merge_shape = (ShapeData *) 0;
+   } // end for mid
 
    return;
 }
@@ -1511,124 +1464,88 @@ void ModeFuzzyEngine::do_obs_merge_thresh() {
 ///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_fcst_merge_thresh(const ShapeData &merge_data) {
-   int j, k, x, y;
-   int n_fcst_merge, intersection;
-   int count, first_k;
-   //ShapeData fcst_merge_mask;
+   const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_fcst_merge;
    ShapeData fcst_merge_split;
 
-   // TODO: reimplement this function to avoid allocating these ShapeData objects
-   ShapeData * fcst_shape = (ShapeData *) 0;
-   ShapeData * fcst_merge_shape = (ShapeData *) 0;
-
    do_fcst_splitting();
 
    //
-   // Define the forecast merge field by applying the specified threshold
-   // to the convolved field
-   //
-   fcst_merge_split = merge_data;
-   n_fcst_merge = fcst_merge_split.n_objects();
-      
-   // fcst_merge_mask = *fcst_conv;
-
-   // //
-   // // Threshold the forecast merge field
-   // //
-   // fcst_merge_mask.threshold(conf_info.Fcst->merge_thresh);
-
-   // //
-   // // Split up the forecast merge field
-   // //
-   // fcst_merge_split = split(fcst_merge_mask, n_fcst_merge);
-
-   //
-   // Allocate space for all of the simple forecast shapes and
-   // forecast merge shapes
+   // Check dimensions
    //
-   fcst_shape       = new ShapeData [n_fcst];
-   fcst_merge_shape = new ShapeData [n_fcst_merge];
-
-   if(!fcst_shape || !fcst_merge_shape) {
-
-      mlog << Error << "\nModeFuzzyEngine::do_fcst_merge_thresh() -> "
-           << "memory allocation error\n\n";
+   if(fcst_split->data.nx() != merge_data.data.nx() ||
+      fcst_split->data.ny() != merge_data.data.ny()) {
+      mlog << Error << "\n" << method_name
+           << "grid dimensions do not match (" << fcst_split->data.nx() << ", "
+           << fcst_split->data.ny() << ") != (" << merge_data.data.nx() << ", "
+           << merge_data.data.ny() << ")!\n\n";
       exit(1);
    }
 
    //
-   // Select all of the simple forecast shapes and
-   // forecast merge shapes
+   // Assume that the input merge data has already been split
    //
-   for(j=0; j<n_fcst; j++) {
-      fcst_shape[j] = select(*fcst_split, j+1);
-   }
-
-   for(j=0; j<n_fcst_merge; j++) {
-      fcst_merge_shape[j] = select(fcst_merge_split, j+1);
-   }
+   fcst_merge_split = merge_data;
+   n_fcst_merge = fcst_merge_split.n_objects();
 
    //
-   // Calculate the composite object sets
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   for(j=0; j<n_fcst_merge; j++) {
-
-      count = first_k = 0;
-
-      for(k=0; k<n_fcst; k++) {
-
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-         for(x=0; x<fcst_merge_split.data.nx(); x++) {
-            for(y=0; y<fcst_merge_split.data.ny(); y++) {
+   // Add map entries for each object
+   for(j=1; j<=n_fcst;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_fcst_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-               if(fcst_merge_shape[j].s_is_on(x, y) &&
-                  fcst_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
+   // Loop over all the data points
+   for(j=0; j<fcst_split->data.nxy(); j++) {
 
-         //
-         // Add to set collection only if the fcst object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= fcst area
-         //
+      oid = fcst_split->data.buf()[j];
+      mid = fcst_merge_split.data.buf()[j];
 
-         if(intersection >= fcst_shape[k].area()) {
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
 
-            collection.make_room();
+   } // end for j
 
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
+   //
+   // Each simple object should be fully contained in a merge object
+   //
+   for(oid=1; oid<=n_fcst; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No forecast threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object!\n\n";
+         return;
+      }
+   }
 
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(first_k, -1);
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(k+1, -1);
-            }
+   //
+   // Calculate the composite object sets
+   //
+   for(mid=1; mid<=n_fcst_merge; mid++) {
 
-            count++;
-         }
-      }  // end for k
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-      if(count > 0) collection.n_sets++;
-   }  // end for j
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(*it, -1);
+      }
 
-   //
-   // Done
-   //
+      // Increment set counter
+      collection.n_sets++;
 
-   delete [] fcst_shape;       fcst_shape = (ShapeData *) 0;
-   delete [] fcst_merge_shape; fcst_merge_shape = (ShapeData *) 0;
+   } // end for mid
 
-   mlog << Debug(4) << " After merging of fcst field: " << fcst_split->sdebug_examine() << "\n";
+   mlog << Debug(4) << " After merging of fcst field: "
+        << fcst_split->sdebug_examine() << "\n";
 
    return;
 }
@@ -1636,119 +1553,88 @@ void ModeFuzzyEngine::do_fcst_merge_thresh(const ShapeData &merge_data) {
 ///////////////////////////////////////////////////////////////////////
 
 void ModeFuzzyEngine::do_obs_merge_thresh(const ShapeData &merge_data) {
-   int j, k, x, y;
-   int n_obs_merge, intersection;
-   int count, first_k;
-   //ShapeData obs_merge_mask;
+   const char *method_name = "ModeFuzzyEngine::do_fcst_merge_thresh() -> ";
+   int j, mid, oid;
+   int n_obs_merge;
    ShapeData obs_merge_split;
-   ShapeData * obs_shape = (ShapeData *) 0;
-   ShapeData * obs_merge_shape = (ShapeData *) 0;
 
    do_obs_splitting();
 
    //
-   // Define the forecast merge field by applying the specified threshold
-   // to the convolved field
+   // Check dimensions
    //
-   // obs_merge_mask = *obs_conv;
-
-   //
-   // Threshold the forecast merge field
-   //
-   // obs_merge_mask.threshold(conf_info.Obs->merge_thresh);
+   if(obs_split->data.nx() != merge_data.data.nx() ||
+      obs_split->data.ny() != merge_data.data.ny()) {
+      mlog << Error << "\n" << method_name
+           << "grid dimensions do not match (" << obs_split->data.nx() << ", "
+           << obs_split->data.ny() << ") != (" << merge_data.data.nx() << ", "
+           << merge_data.data.ny() << ")!\n\n";
+      exit(1);
+   }
 
    //
-   // Split up the forecast merge field
+   // Assume that the input merge data has already been split
    //
-   obs_merge_split = merge_data; //split(obs_merge_mask, n_obs_merge);
+   obs_merge_split = merge_data;
    n_obs_merge = obs_merge_split.n_objects();
 
    //
-   // Allocate space for all of the simple observation shapes and
-   // observation merge shapes
+   // Mapping of simple object id to merge object id values and vice-versa
    //
-   obs_shape       = new ShapeData [n_obs];
-   obs_merge_shape = new ShapeData [n_obs_merge];
+   map<int, std::set<int> > obj_to_mrg_map;
+   map<int, std::set<int> > mrg_to_obj_map;
+   std::set<int> empty_set;
 
-   if(!obs_shape || !obs_merge_shape) {
+   // Add map entries for each object
+   for(j=1; j<=n_obs;       j++) obj_to_mrg_map[j] = empty_set;
+   for(j=1; j<=n_obs_merge; j++) mrg_to_obj_map[j] = empty_set;
 
-      mlog << Error << "\nModeFuzzyEngine::do_obs_merge_thresh() -> "
-           << "memory allocation error\n\n";
-      exit(1);
-   }
+   // Loop over all the data points
+   for(j=0; j<obs_split->data.nxy(); j++) {
+
+      oid = obs_split->data.buf()[j];
+      mid = obs_merge_split.data.buf()[j];
+
+      // Update object id maps
+      if(oid > 0)            obj_to_mrg_map[oid].insert(mid);
+      if(oid > 0 && mid > 0) mrg_to_obj_map[mid].insert(oid);
+
+   } // end for j
 
    //
-   // Select all of the simple observation shapes and
-   // simple merge shapes
+   // Each simple object should be fully contained in a merge object
    //
-   for(j=0; j<n_obs; j++) {
-      obs_shape[j] = select(*obs_split, j+1);
-   }
-
-   for(j=0; j<n_obs_merge; j++) {
-      obs_merge_shape[j] = select(obs_merge_split, j+1);
+   for(oid=1; oid<=n_obs; oid++) {
+      if(obj_to_mrg_map[oid].size() != 1) {
+         mlog << Warning << "\n" << method_name
+              << "No observation threshold merging done because simple object "
+              << oid << " is not fully contained within one merge object!\n\n";
+         return;
+      }
    }
 
    //
    // Calculate the composite object sets
    //
-   for(j=0; j<n_obs_merge; j++) {
+   for(mid=1; mid<=n_obs_merge; mid++) {
 
-      count = first_k = 0;
-
-      for(k=0; k<n_obs; k++) {
+      // Ignore merge objects containing a single simple object
+      if(mrg_to_obj_map[mid].size() <= 1) continue;
 
-         //
-         // Calculate intersection area
-         //
-         intersection = 0;
-
-         for(x=0; x<obs_merge_split.data.nx(); x++) {
-            for(y=0; y<obs_merge_split.data.ny(); y++) {
-
-               if(obs_merge_shape[j].s_is_on(x, y) &&
-                  obs_shape[k].s_is_on(x, y)) intersection++;
-            }
-         }
-
-         //
-         // Add to set collection only if the obs object is
-         // completely contained in the merge object.  Meaning,
-         // intersection area >= obs area
-         //
-         if(intersection >= obs_shape[k].area()) {
-
-            collection.make_room();
-
-            //
-            // Keep track of the first embedded shape.  You only want to
-            // create a composite if there are more than 1 shapes in it.
-            //
-            if(count == 0) first_k = k+1;
-
-            else if(count == 1) {
-               collection.set[collection.n_sets].add_pair(-1, first_k);
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-            else {
-               collection.set[collection.n_sets].add_pair(-1, k+1);
-            }
-
-            count++;
-         }
-      }  // end for k
-
-      if(count > 0) collection.n_sets++;
-   }  // end for j
+      // Add all simple objects into one new set
+      collection.make_room();
+      for(auto it  = mrg_to_obj_map[mid].begin();
+               it != mrg_to_obj_map[mid].end(); it++) {
+         collection.set[collection.n_sets].add_pair(-1, *it);
+      }
 
-   //
-   // Done
-   //
+      // Increment set counter
+      collection.n_sets++;
 
-   delete [] obs_shape;       obs_shape = (ShapeData *) 0;
-   delete [] obs_merge_shape; obs_merge_shape = (ShapeData *) 0;
+   } // end for mid
 
-   mlog << Debug(4) << " After merging of obs field: " << obs_split->sdebug_examine() << "\n";
+   mlog << Debug(4) << " After merging of obs field: "
+        << obs_split->sdebug_examine() << "\n";
 
    return;
 }