diff --git a/Cargo.toml b/Cargo.toml
index 5f60f8a7..ed1384f8 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -9,12 +9,16 @@ license = "Apache-2.0/MIT"
 keywords = ["math", "random"]
 authors = ["The Noise-rs Developers."]
 
+[profile.release]
+debug = true
+
 [lib]
 name = "noise"
 
 [dependencies]
 rand = "0.5"
 image = { version = "0.18", optional = true }
+rayon = "1.2.0"
 
 [features]
 default = ["image"]
diff --git a/src/utils/noise_map_builder.rs b/src/utils/noise_map_builder.rs
index c2feb422..135260e6 100644
--- a/src/utils/noise_map_builder.rs
+++ b/src/utils/noise_map_builder.rs
@@ -1,6 +1,10 @@
+extern crate rayon;
+
 use math::interp;
 use noise_fns::NoiseFn;
 use utils::noise_map::NoiseMap;
+use std::sync::{Arc, RwLock};
+use self::rayon::prelude::*;
 
 pub trait NoiseMapBuilder<'a> {
     fn set_size(self, width: usize, height: usize) -> Self;
@@ -205,34 +209,45 @@ impl<'a> NoiseMapBuilder<'a> for PlaneMapBuilder<'a> {
         let x_step = x_extent / width as f64;
         let y_step = y_extent / height as f64;
 
+        let x_bounds = self.x_bounds;
+        let y_bounds = self.y_bounds;
+        let is_seamless = self.is_seamless;
+
         for y in 0..height {
             let current_y = self.y_bounds.0 + y_step * y as f64;
-
-            for x in 0..width {
-                let current_x = self.x_bounds.0 + x_step * x as f64;
-
-                let final_value = if self.is_seamless {
-                    let sw_value = self.source_module.get([current_x, current_y, 0.0]);
-                    let se_value = self.source_module
-                        .get([current_x + x_extent, current_y, 0.0]);
-                    let nw_value = self.source_module
-                        .get([current_x, current_y + y_extent, 0.0]);
-                    let ne_value =
-                        self.source_module
-                            .get([current_x + x_extent, current_y + y_extent, 0.0]);
-
-                    let x_blend = 1.0 - ((current_x - self.x_bounds.0) / x_extent);
-                    let y_blend = 1.0 - ((current_y - self.y_bounds.0) / y_extent);
-
-                    let y0 = interp::linear(sw_value, se_value, x_blend);
-                    let y1 = interp::linear(nw_value, ne_value, x_blend);
-
-                    interp::linear(y0, y1, y_blend)
+            let locked = RwLock::new(self);
+
+            // split tasks into parallel threads, then collect
+            let final_values: Vec<_> = (0..width).into_par_iter().map(|x| {
+                let current_x = x_bounds.0 + x_step * x as f64;
+                let readable =  locked.read().unwrap();
+
+                let final_value = if is_seamless {
+                  let sw_value = readable.source_module.get([current_x, current_y, 0.0]);
+                  let se_value = readable.source_module
+                    .get([current_x + x_extent, current_y, 0.0]);
+                  let nw_value = readable.source_module
+                    .get([current_x, current_y + y_extent, 0.0]);
+                  let ne_value =
+                    readable.source_module
+                      .get([current_x + x_extent, current_y + y_extent, 0.0]);
+                  let x_blend = 1.0 - ((current_x - x_bounds.0) / x_extent);
+                  let y_blend = 1.0 - ((current_y - y_bounds.0) / y_extent);
+
+                  let y0 = interp::linear(sw_value, se_value, x_blend);
+                  let y1 = interp::linear(nw_value, ne_value, x_blend);
+
+                  interp::linear(y0, y1, y_blend)
                 } else {
-                    self.source_module.get([current_x, current_y, 0.0])
+                  readable.source_module.get([current_x, current_y, 0.0])
                 };
 
-                result_map.set_value(x, y, final_value);
+                final_value
+            }).collect();
+
+            // write results back to data structure
+            for x in 0..width {
+              result_map.set_value(x, y, final_values[x]);
             }
         }