diff --git a/opendm/concurrency.py b/opendm/concurrency.py
index 66fb26e8..dfcf58f4 100644
--- a/opendm/concurrency.py
+++ b/opendm/concurrency.py
@@ -25,6 +25,9 @@ def get_max_memory_mb(minimum = 100, use_at_most = 0.5):
     """
     return max(minimum, (virtual_memory().available / 1024 / 1024) * use_at_most)
 
+def get_total_memory():
+    return virtual_memory().total
+
 def parallel_map(func, items, max_workers=1, single_thread_fallback=True):
     """
     Our own implementation for parallel processing
diff --git a/opendm/dem/commands.py b/opendm/dem/commands.py
index 1631202d..8a1afbcd 100755
--- a/opendm/dem/commands.py
+++ b/opendm/dem/commands.py
@@ -11,7 +11,7 @@ from opendm.system import run
 from opendm import point_cloud
 from opendm import io
 from opendm import system
-from opendm.concurrency import get_max_memory, parallel_map
+from opendm.concurrency import get_max_memory, parallel_map, get_total_memory
 from scipy import ndimage
 from datetime import datetime
 from opendm.vendor.gdal_fillnodata import main as gdal_fillnodata
@@ -69,7 +69,7 @@ error = None
 def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56'], gapfill=True,
                 outdir='', resolution=0.1, max_workers=1, max_tile_size=4096,
                 decimation=None, keep_unfilled_copy=False,
-                apply_smoothing=True):
+                apply_smoothing=True, max_tiles=None):
     """ Create DEM from multiple radii, and optionally gapfill """
     
     global error
@@ -149,6 +149,11 @@ def create_dem(input_point_cloud, dem_type, output_type='max', radiuses=['0.56']
 
                 miny = maxy
             minx = maxx
+    
+    # Safety check
+    if max_tiles is not None:
+        if len(tiles) > max_tiles and final_dem_pixels > get_total_memory() * 10:
+            raise system.ExitException("Max tiles limit exceeded (%s). This is a strong indicator that the reconstruction failed and we would probably run out of memory trying to process this" % max_tiles)
 
     # Sort tiles by increasing radius
     tiles.sort(key=lambda t: float(t['radius']), reverse=True)
diff --git a/opendm/mesh.py b/opendm/mesh.py
index da393b79..58e94378 100644
--- a/opendm/mesh.py
+++ b/opendm/mesh.py
@@ -9,7 +9,7 @@ from opendm import point_cloud
 from scipy import signal
 import numpy as np
 
-def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, available_cores=None, method='gridded', smooth_dsm=True):
+def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=0.05, depth=8, samples=1, maxVertexCount=100000, available_cores=None, method='gridded', smooth_dsm=True, max_tiles=None):
     # Create DSM from point cloud
 
     # Create temporary directory
@@ -31,7 +31,8 @@ def create_25dmesh(inPointCloud, outMesh, radius_steps=["0.05"], dsm_resolution=
             outdir=tmp_directory,
             resolution=dsm_resolution,
             max_workers=available_cores,
-            apply_smoothing=smooth_dsm
+            apply_smoothing=smooth_dsm,
+            max_tiles=max_tiles
         )
 
     if method == 'gridded':
diff --git a/stages/odm_dem.py b/stages/odm_dem.py
index 0d2ea26a..d9029eea 100755
--- a/stages/odm_dem.py
+++ b/stages/odm_dem.py
@@ -100,7 +100,8 @@ class ODMDEMStage(types.ODM_Stage):
                             resolution=resolution / 100.0,
                             decimation=args.dem_decimation,
                             max_workers=args.max_concurrency,
-                            keep_unfilled_copy=args.dem_euclidean_map
+                            keep_unfilled_copy=args.dem_euclidean_map,
+                            max_tiles=math.ceil(len(reconstruction.photos) / 2)
                         )
 
                     dem_geotiff_path = os.path.join(odm_dem_root, "{}.tif".format(product))
diff --git a/stages/odm_meshing.py b/stages/odm_meshing.py
index bfef268a..a113a989 100644
--- a/stages/odm_meshing.py
+++ b/stages/odm_meshing.py
@@ -59,7 +59,8 @@ class ODMeshingStage(types.ODM_Stage):
                     samples=self.params.get('samples'),
                     available_cores=args.max_concurrency,
                     method='poisson' if args.fast_orthophoto else 'gridded',
-                    smooth_dsm=True)
+                    smooth_dsm=True,
+                    max_tiles=math.ceil(len(reconstruction.photos) / 2))
           else:
               log.ODM_WARNING('Found a valid ODM 2.5D Mesh file in: %s' %
                               tree.odm_25dmesh)