diff --git a/point_cloud_colorize/las_colorize.py b/point_cloud_colorize/las_colorize.py
index be43a55..c545e65 100644
--- a/point_cloud_colorize/las_colorize.py
+++ b/point_cloud_colorize/las_colorize.py
@@ -2,15 +2,17 @@
 """
 Python3
 
-@author: Chris Lucas
+@author: Chris Lucas, Arno Timmer
 """
 
 import argparse
 import json
 from pathlib import Path
-
+import shutil
 import pdal
-
+import datetime
+from joblib import Parallel, delayed
+import os
 
 PDAL_PIPELINE = """{{
   "pipeline":[
@@ -73,26 +75,114 @@ def run_pdal(input_path, output_path, las_srs, wms_url,
                 'wms_pixel_size': wms_pixel_size,
                 'wms_max_image_size': wms_max_image_size,
                 'las_srs': las_srs}
+
     pdalargs_str = json.dumps(pdalargs).replace('"', '\\"')
 
-    path = Path(__file__)
+    path = Path(os.getcwd())
+
     pipeline_json = PDAL_PIPELINE.format(input_file=input_path.as_posix(),
                                          output_file=output_path.as_posix(),
                                          srs=las_srs,
                                          pdalargs=pdalargs_str,
-                                         directory=path.parent.as_posix())
+                                         directory=path.as_posix())
     pipeline = pdal.Pipeline(pipeline_json)
     pipeline.validate()
     pipeline.execute()
 
 
+def parallel_coloring(f, i, verbose, tmp_col_path, pdalargs, tmp_col='tmp_col_{}.laz'):
+    tmp_col = Path(tmp_col_path.joinpath(tmp_col.format(i)))
+    run_pdal(Path(f), tmp_col,
+             pdalargs['las_srs'], pdalargs['wms_url'],
+             pdalargs['wms_layer'], pdalargs['wms_srs'],
+             pdalargs['wms_version'], pdalargs['wms_format'],
+             pdalargs['wms_pixel_size'], pdalargs['wms_max_image_size'])
+
+    if verbose:
+        print(f'colored {i} parts of the las at {datetime.datetime.now()}')
+
+
+def process_files_parallel(input, output, las_srs,
+                           wms_url, wms_layer, wms_srs,
+                           wms_version, wms_format,
+                           wms_pixel_size, wms_max_image_size,
+                           verbose):
+    """
+    :param input_path:
+    :param output_path:
+    :param las_srs:
+    :param verbose:
+    :return:
+    """
+    if verbose:
+        print(f'started colorizing parts at {datetime.datetime.now()}')
+    output_path = Path(output)
+    if not output_path.is_dir():
+        raise ValueError('Output should be a directory')
+
+    input_path = Path(input)
+    tmp_div_path = Path(output_path.parent.joinpath('tmp_div'))
+
+    if tmp_div_path.exists():
+        if verbose:
+            print('Temporary path exists, deleting.')
+        shutil.rmtree(tmp_div_path)
+    tmp_div_path.mkdir(parents=True, exist_ok=True)
+
+    divide_pipeline = """{{
+                          "pipeline":[
+                            {{
+                              "type": "readers.las",
+                              "filename": "{input_file}"
+                            }},
+                            {{
+                              "type": "filters.divider",
+                              "count": "6"
+                            }},
+                            {{
+                              "type": "writers.las",
+                              "a_srs": "{srs}",
+                              "filename": "{output_path}/tmp_#.laz"
+                            }}
+                          ]}}"""
+
+    # create pipeline for dividing the las
+    div_pipeline_json = divide_pipeline.format(input_file=input_path.as_posix(),
+                                               srs=las_srs,
+                                               output_path=tmp_div_path.as_posix())
+    pipeline = pdal.Pipeline(div_pipeline_json)
+    pipeline.validate()
+    pipeline.execute()
+
+    if verbose:
+        print(f'las is divided at {datetime.datetime.now()}')
+
+    pdalargs = {'wms_url': wms_url,
+                'wms_layer': wms_layer,
+                'wms_srs': wms_srs,
+                'wms_version': wms_version,
+                'wms_format': wms_format,
+                'wms_pixel_size': wms_pixel_size,
+                'wms_max_image_size': wms_max_image_size,
+                'las_srs': las_srs}
+
+    # for each of the created las-parts
+    if verbose:
+        print(f'start parallel processing at {datetime.datetime.now()}')
+
+    Parallel(n_jobs=6)(delayed(parallel_coloring)(f, i, verbose, output_path, pdalargs) for i, f in
+                       enumerate(Path(tmp_div_path).iterdir(), 1))
+
+    if verbose:
+        print(f'colorizing in parts finished at {datetime.datetime.now()}')
+
+
 def process_files(input_path, output_path, las_srs,
                   wms_url, wms_layer, wms_srs,
                   wms_version, wms_format, wms_pixel_size,
                   wms_max_image_size, verbose=False):
     """
     Run the pdal pipeline for the input files.
-
     Parameters
     ----------
     input_path : str
@@ -224,7 +314,14 @@ def argument_parser():
                               '(int, default: 1000)'),
                         required=False,
                         default=1000)
-    parser.add_argument('-V', '--verbose', default=False, action="store_true",
+    parser.add_argument('-d', '--divide',
+                        default=5,
+                        action="store_true",
+                        help='Divide the point cloud in a given number of '
+                             'smaller areas which are colored seperately')
+    parser.add_argument('-V', '--verbose',
+                        default=False,
+                        action="store_true",
                         help='Set verbose.')
     args = parser.parse_args()
     return args
@@ -235,11 +332,18 @@ def main():
     Run the application.
     """
     args = argument_parser()
-    process_files(args.input, args.output, args.las_srs,
-                  args.wms_url, args.wms_layer, args.wms_srs,
-                  args.wms_version, args.wms_format,
-                  args.wms_pixel_size, args.wms_max_image_size,
-                  args.verbose)
+    if args.divide:
+        process_files_parallel(args.input, args.output, args.las_srs,
+                               args.wms_url, args.wms_layer, args.wms_srs,
+                               args.wms_version, args.wms_format,
+                               args.wms_pixel_size, args.wms_max_image_size,
+                               args.verbose)
+    else:
+        process_files(args.input, args.output, args.las_srs,
+                      args.wms_url, args.wms_layer, args.wms_srs,
+                      args.wms_version, args.wms_format,
+                      args.wms_pixel_size, args.wms_max_image_size,
+                      args.verbose)
 
 
 if __name__ == '__main__':
diff --git a/point_cloud_colorize/tests/__init__.py b/point_cloud_colorize/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/point_cloud_colorize/tests/create_test_data.py b/point_cloud_colorize/tests/create_test_data.py
new file mode 100644
index 0000000..cf2c6ac
--- /dev/null
+++ b/point_cloud_colorize/tests/create_test_data.py
@@ -0,0 +1,26 @@
+import json
+from typing import Tuple
+
+import pdal
+import wget
+
+
+def crop_tile(bounds: Tuple[float, float, float, float] = (139703.2, 556571.8, 139717.1, 556581.8)) -> None:
+    x_min, y_min, x_max, y_max = bounds
+    pdal_pipeline = {
+        "pipeline": [
+            "09HN2.LAZ",
+            {"type": "filters.crop", "bounds": f"([{x_min}, {y_min}], [{x_max}, {y_max}])"},
+            {"type": "writers.las", "filename": "test.laz"}
+        ]
+    }
+    pipeline = pdal.Pipeline(json.dumps(pdal_pipeline))
+    pipeline.validate()
+    pipeline.execute()
+
+
+if __name__ == '__main__':
+    url = 'https://geodata.nationaalgeoregister.nl/ahn3/extract/ahn3_laz/C_09HN2.LAZ'
+    wget.download(url, '09HN2.LAZ')
+
+    crop_tile()
diff --git a/point_cloud_colorize/tests/data/input/README.md b/point_cloud_colorize/tests/data/input/README.md
new file mode 100644
index 0000000..fab9743
--- /dev/null
+++ b/point_cloud_colorize/tests/data/input/README.md
@@ -0,0 +1,2 @@
+# Input data directory
+This is the directory for storing the input data
\ No newline at end of file
diff --git a/point_cloud_colorize/tests/data/input/simple.laz b/point_cloud_colorize/tests/data/input/simple.laz
new file mode 100644
index 0000000..a867f27
Binary files /dev/null and b/point_cloud_colorize/tests/data/input/simple.laz differ
diff --git a/point_cloud_colorize/tests/data/output/README.md b/point_cloud_colorize/tests/data/output/README.md
new file mode 100644
index 0000000..3ceef65
--- /dev/null
+++ b/point_cloud_colorize/tests/data/output/README.md
@@ -0,0 +1,2 @@
+# Output data directory
+This is the directory for storing the output data
\ No newline at end of file
diff --git a/point_cloud_colorize/tests/test.laz b/point_cloud_colorize/tests/test.laz
new file mode 100644
index 0000000..52d57cd
Binary files /dev/null and b/point_cloud_colorize/tests/test.laz differ
diff --git a/point_cloud_colorize/tests/test_colorize.py b/point_cloud_colorize/tests/test_colorize.py
new file mode 100644
index 0000000..7a4079c
--- /dev/null
+++ b/point_cloud_colorize/tests/test_colorize.py
@@ -0,0 +1,42 @@
+import os
+import unittest
+
+from point_cloud_colorize.las_colorize import process_files, process_files_parallel
+
+class TestColorizePointcloud(unittest.TestCase):
+    def test_unparrallel_processing(self) -> None:
+        with self.subTest('It fails if the colorizing fails'):
+            input_laz = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'test.laz')
+            output_laz = 'test_out.laz'
+
+            process_files(input_path=input_laz,
+                          output_path=output_laz,
+                          las_srs='EPSG:28992',
+                          wms_url='https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wms?',
+                          wms_layer='Actueel_ortho25',
+                          wms_srs='EPSG:28992',
+                          wms_version='1.3.0',
+                          wms_format='image/png',
+                          wms_pixel_size=0.25,
+                          wms_max_image_size=1000,
+                          verbose=False)
+
+    def test_parrallel_processing(self) -> None:
+        with self.subTest('It fails if the colorizing fails'):
+
+            input_laz = os.path.join(os.path.dirname(os.path.abspath(__file__)),'test.laz')
+            output_laz = 'test_out_dir'
+
+            process_files_parallel(input=input_laz,
+                                   output=output_laz,
+                                   las_srs='EPSG:28992',
+                                   wms_url='https://geodata.nationaalgeoregister.nl/luchtfoto/rgb/wms?',
+                                   wms_layer='Actueel_ortho25',
+                                   wms_srs='EPSG:28992',
+                                   wms_version='1.3.0',
+                                   wms_format='image/png',
+                                   wms_pixel_size=0.25,
+                                   wms_max_image_size=1000,
+                                   verbose=False)
+
+