Fix lookuptable

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "extern/argh"]
+	path = extern/argh
+	url = https://github.com/adishavit/argh.git

CMakeLists.txt

@@ -1,15 +1,22 @@
+cmake_minimum_required(VERSION 3.20)
+
 project(
   contour-creator
   LANGUAGES CXX)    
   
-find_package(GDAL CONFIG REQUIRED)
-
 add_executable(${PROJECT_NAME}
 	src/HeightMap.cpp src/CellMap.cpp src/main.cpp
 )
+
+# Argh is a simple argrument parser
+add_subdirectory(extern/argh)
+target_link_libraries(${PROJECT_NAME} PRIVATE argh)
+
+# Gdal is used for geodata IO
+find_package(GDAL CONFIG REQUIRED)
+target_link_libraries(${PROJECT_NAME} PRIVATE GDAL::GDAL)
+
 find_package(OpenMP)
 if(OpenMP_CXX_FOUND)
     target_link_libraries(${PROJECT_NAME} PUBLIC OpenMP::OpenMP_CXX)
 endif()
-
-target_link_libraries(${PROJECT_NAME} GDAL::GDAL)

Doxyfile

@@ -42,7 +42,7 @@ DOXYFILE_ENCODING      = UTF-8
 # title of most generated pages and in a few other places.
 # The default value is: My Project.
 
-PROJECT_NAME           = "Marching Squares"
+PROJECT_NAME           = "Contour Creator"
 
 # The PROJECT_NUMBER tag can be used to enter a project or revision number. This
 # could be handy for archiving the generated documentation or if some version

README.md

@@ -6,8 +6,19 @@ This is our project for the INF205: Resource-efficient programming course
 - [x] Run the marching squares algorithm and produce a "cellmap"
 - [ ] Use a lookuptable to produce a vector file from the "cellmap"
 ## Dependencies
-- GDAL
+- GDAL >= 3.5
 - OpenMP
+- CMake
+
+If you want to clone the repo with the example files you need [git-lfs](https://git-lfs.com/) installed and activated with ´git lfs install´
+To install the packages on Fedora run  
+```
+dnf install g++ git-lfs cmake gdal-devel
+```
+on Debian:
+```
+apt install g++ git-lfs cmake libgdal-dev
+```
 ## How to build:
 ```
 mkdir build
@@ -16,3 +27,8 @@ cmake -DCMAKE_BUILD_TYPE=Debug ..
 cmake --build . --parallel
 ```
 Then you can run `./contour-creator PATH/TO/HEIGTHMAP.TIF`  
+Run an example file: `./contour-creator "../example_files/Follo 2014-dtm.tif"`  
+Currently the program outputs to ´out.shp´ in the current directory.
+
+## Docs
+[Doxygen output](https://trygve.gitlab.io/contour-creator/)

build&run.sh

@@ -3,4 +3,4 @@ rm -rf build
 mkdir build
 cmake -DCMAKE_BUILD_TYPE=Debug -B build
 cmake --build build --parallel 
-build/contour-creator 'example_files/crop.tif'
+build/contour-creator --interval=1 --blur example_files/crop.tif

documentation/ER_diagram.gaphor

@@ -113,7 +113,7 @@
 </Property>
 <Class id="4a8d98db-e76f-11ee-9ded-8bb162b1502a">
 <name>
-<val>Cells</val>
+<val>Cell</val>
 </name>
 <ownedAttribute>
 <reflist>

documentation/ER_diagram.svg

@@ -35,9 +35,6 @@
 <g id="glyph-0-10">
 <path d="M 1.171875 -10.640625 L 3.625 -10.640625 L 3.625 0 L 1.171875 0 Z M 1.171875 -10.640625 "/>
 </g>
-<g id="glyph-0-11">
-<path d="M 7.15625 -7.421875 L 7.15625 -5.5625 C 6.632812 -5.78125 6.128906 -5.941406 5.640625 -6.046875 C 5.148438 -6.160156 4.691406 -6.21875 4.265625 -6.21875 C 3.796875 -6.21875 3.445312 -6.15625 3.21875 -6.03125 C 3 -5.914062 2.890625 -5.738281 2.890625 -5.5 C 2.890625 -5.300781 2.972656 -5.148438 3.140625 -5.046875 C 3.304688 -4.941406 3.613281 -4.863281 4.0625 -4.8125 L 4.484375 -4.765625 C 5.742188 -4.597656 6.585938 -4.332031 7.015625 -3.96875 C 7.453125 -3.601562 7.671875 -3.03125 7.671875 -2.25 C 7.671875 -1.4375 7.367188 -0.820312 6.765625 -0.40625 C 6.160156 0 5.265625 0.203125 4.078125 0.203125 C 3.566406 0.203125 3.039062 0.160156 2.5 0.078125 C 1.96875 -0.00390625 1.414062 -0.125 0.84375 -0.28125 L 0.84375 -2.140625 C 1.332031 -1.898438 1.832031 -1.71875 2.34375 -1.59375 C 2.851562 -1.476562 3.375 -1.421875 3.90625 -1.421875 C 4.382812 -1.421875 4.742188 -1.488281 4.984375 -1.625 C 5.222656 -1.757812 5.34375 -1.957031 5.34375 -2.21875 C 5.34375 -2.4375 5.257812 -2.597656 5.09375 -2.703125 C 4.925781 -2.804688 4.597656 -2.890625 4.109375 -2.953125 L 3.671875 -3.015625 C 2.578125 -3.148438 1.8125 -3.398438 1.375 -3.765625 C 0.9375 -4.140625 0.71875 -4.703125 0.71875 -5.453125 C 0.71875 -6.265625 0.992188 -6.863281 1.546875 -7.25 C 2.109375 -7.644531 2.960938 -7.84375 4.109375 -7.84375 C 4.566406 -7.84375 5.039062 -7.804688 5.53125 -7.734375 C 6.03125 -7.671875 6.570312 -7.566406 7.15625 -7.421875 Z M 7.15625 -7.421875 "/>
-</g>
 <g id="glyph-1-0">
 <path d="M 6.4375 -8.78125 L 6.4375 -4.96875 L 10.25 -4.96875 L 10.25 -3.8125 L 6.4375 -3.8125 L 6.4375 0 L 5.296875 0 L 5.296875 -3.8125 L 1.484375 -3.8125 L 1.484375 -4.96875 L 5.296875 -4.96875 L 5.296875 -8.78125 Z M 6.4375 -8.78125 "/>
 </g>
@@ -251,11 +248,10 @@
 </g>
 <path fill="none" stroke-width="2" stroke-linecap="butt" stroke-linejoin="round" stroke="rgb(0%, 0%, 0%)" stroke-opacity="1" stroke-miterlimit="10" d="M 0.00103399 0 L 311.001034 0 L 311.001034 117 L 0.00103399 117 Z M 0.00103399 0 Z M 0.00103399 0 " transform="matrix(1, 0, 0, 1, 363.924747, 26.042969)"/>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">
-<use xlink:href="#glyph-0-9" x="500.925781" y="57.042969"/>
+<use xlink:href="#glyph-0-9" x="504.925781" y="57.042969"/>
-<use xlink:href="#glyph-0-1" x="510.925781" y="57.042969"/>
+<use xlink:href="#glyph-0-1" x="514.925781" y="57.042969"/>
-<use xlink:href="#glyph-0-10" x="519.925781" y="57.042969"/>
+<use xlink:href="#glyph-0-10" x="523.925781" y="57.042969"/>
-<use xlink:href="#glyph-0-10" x="524.925781" y="57.042969"/>
+<use xlink:href="#glyph-0-10" x="528.925781" y="57.042969"/>
-<use xlink:href="#glyph-0-11" x="529.925781" y="57.042969"/>
 </g>
 <path fill="none" stroke-width="2" stroke-linecap="butt" stroke-linejoin="round" stroke="rgb(0%, 0%, 0%)" stroke-opacity="1" stroke-miterlimit="10" d="M 0.00103399 41 L 311.001034 41 " transform="matrix(1, 0, 0, 1, 363.924747, 26.042969)"/>
 <g fill="rgb(0%, 0%, 0%)" fill-opacity="1">

documentation/lab5.md

@@ -1,18 +1,59 @@
-# 29. Glossary:
+% INF205 Lab 5
+% Esther and Trygve
+% April 23 2024
+
+# 29. Glossary contributions:
 ## dynamic library
-Dynamic libraries are libraries are pre compiled and often provided by the operating system, as opposed to compiling it as part of your project. This saves space as each program dont need its own "copy".
+Dynamic libraries are libraries that are pre compiled and often provided by the operating system, as opposed to compiling them as part of your project. This saves space as each program dont need its own "copy".
+
+## command-line arguments
+Command-line arguments are extra commands that we enter after the program's executable name so that the functionality of the program changes. 
+
 
 # 30. Draft slides:
 See slides.pdf
 
 # 31. Basic functionality and validation per
-...
+
+We ran a simple benchmark against gdal_contour as the reference implementation using `time` in the shell. My laptop has a 6 core AMD Ryzen 5 PRO 5675U running linux 6.8. The average of three tests was:
+```
+Ours:      33,53 s
+Reference:  4,46 s
+```
+This is despite the reference gdal implementation beinge single threaded and our being paritally paralell. We think we will gaing a significant speedup from optimizing the part of he program drawing the contours into a file. Right now its a placeholder that just creates points, not lines.  
+This innital gap in performance leaves us a lot of room to improve and it will be interesting to see what the final performance will be.  
+The plan is to create a more sophisticated benchmark program in c++, but for now we just use the time command
+
+Memory usage is also very big because of our datastructure for CellMap. It is storing essentially a line in a 2d array, so most of the items are 0. We will try to fix this by making a cell type that stores its coordinates or by getting rid of the cellmap entirely.
+
+The source code is available at [https://gitlab.com/Trygve/contour-creator](https://gitlab.com/Trygve/contour-creator) with instructions in the readme for building and running.
+If you want to view the .shp file you can use this simple webapp: [https://mapshaper.org/](https://mapshaper.org/)
 
 # 32 Special interest functionality and responsibilities
-Per now we dont have time to add any extra functionality.
+It does not seem like we will have time to add any extra functionality.
 
-33. Programming project: Progress on data structure implementation
+# 33. Programming project: Progress on data structure implementation
 
-...
+![ER diagram](ER_diagram.svg)
 
-34. Parallelization
+We need to create destructors to call free on data in HeightMap and cells in Cells as those are arrays allocated with malloc. The default copy behaviur is fine for our program.
+
+\newpage
+
+# 34. Parallelization
+We are using openMP because we want multiple threads using the same data and the files are small enough to fit into memory (~500mb). To scale it to multiple machines we would use mpi to distrubute seperate iamge tiles to each node. 
+We haven't had time to test the performance and we have only parallelized half of what is possible. But we can see that it is using all the cores.  
+The following code was used for parallelization using openMP, each thread produces one "layer" of the cellmap:
+```cpp
+#pragma omp parallel
+    {
+        std::vector<CellMap> vec_private;
+        #pragma omp for
+        for (int i = 1; i <= num_contours; i++) 
+        {
+            vec_private.push_back(produce_cellmap(heightmap, heightmap->min + interval*i));
+        }
+        #pragma omp critical
+        vector_contours.insert(vector_contours.end(), vec_private.begin(), vec_private.end());
+    }
+```

documentation/slides.md

@@ -0,0 +1,35 @@
+---
+title:
+- "INF205: creating contours using the marching squares algorithm"
+author:
+- Trygve og Esther
+
+---
+
+# Output of our program:
+![Contour map of Ås](kurver_ås.png)
+
+# Marching squares
+![By Nicoguaro; Adroitwhiz - File:Marching_squares_algorithm.svg, CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=143418710](ms_wikipedia.svg)
+
+# The logical flow of our program
+```cpp
+int main(int argc, const char* argv[])
+{
+    const char* filepath = argv[1];
+    HeightMap map(filepath); 
+    map.blur(0.8)
+    auto lines = create_lines(&map, 5);
+    write_output_file(cellmaps, "out.geojson", &map);
+}
+```
+
+# Performance
+|                    |       Time         |
+|--------------------|--------------------|
+|  12 threads        |     41s            |
+|   1 thread         |     116s           |
+
+
+# Problems we encountered
+

src/CellMap.cpp

@@ -4,13 +4,13 @@
 #include "contour_creator.hh"
 
 
-CellMap::CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system)
+CellMap::CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system, double* geotransform)
 {
     this->width = width;
     this->height = height;
     this->cells = cells;
     this->reference_system = reference_system;
-
+    this->geotransform = geotransform;
 }
 int CellMap::get_cell(int x, int y)
 {

src/HeightMap.cpp

@@ -1,6 +1,7 @@
+#include <gdal/cpl_conv.h>
 #include <iostream>
 #include <stdexcept>
-
+#include <omp.h>
 #include <gdal/gdal.h>
 #include "gdal/gdal_priv.h"
 #include <gdal/gdal_frmts.h>
@@ -28,8 +29,12 @@ HeightMap::HeightMap(const char* filepath)
     this->height = band->GetYSize();
     this->min = band->GetMinimum();
     this->max = band-> GetMaximum();
+    //https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset15GetGeoTransformEPd
+    this->geotransform = (double *) CPLMalloc(sizeof(double)*6);
     this->reference_system = *(file->GetSpatialRef());
 
+    file->GetGeoTransform(this->geotransform);
+    
     this->data = (float *) CPLMalloc(sizeof(float)*width*height);
     CPLErr error = band->RasterIO( GF_Read, 0, 0, width, height,
                     this->data, width, height, GDT_Float32,
@@ -44,3 +49,45 @@ float HeightMap::get_pixel(int x, int y)
     //std::cout << "  offset: " << offset << " " << x << ","<< y;
     return *(this->data + offset);
 }
+
+void HeightMap::blur(float standard_deviation)
+{
+    // standard_deviation does not do anything yet. This is currently a simple box blur
+    int kernel_height = 5;
+    int kernel_width = 5;
+    int kernel_size = kernel_height*kernel_width;
+    float* kernel = (float*) CPLMalloc(sizeof(float)*kernel_size);
+    for (int i=0; i<kernel_size; i++)
+    {
+        *(kernel + i) = 1.0;
+    }
+
+    float* blurred = (float *) CPLMalloc(sizeof(float)*this->width*this->height);
+    #pragma omp parallel
+    {
+    #pragma omp for
+    for (int i = 0; i < this->height * this->width; i++)
+    {
+        float blurred_pixel = 0;
+        for (int j = 0; j < kernel_size; j++)
+        {
+            int x = j%kernel_width - kernel_width/2 + i%this->width;
+            int y = j/kernel_width - kernel_height/2 + i/this->width;
+            if (x<0 || x>=this->width)
+            {
+                x = 0;
+            }
+            
+            if (y<0 || y>=this->height)
+            {
+                y = 0;
+            }
+            blurred_pixel += this->get_pixel(x, y);
+        }
+        *(blurred + i) = blurred_pixel/float(kernel_size);
+    }
+    }
+    free(this->data);
+    this->data = blurred;
+    blurred = nullptr;
+}

src/contour_creator.hh

@@ -9,6 +9,8 @@ class HeightMap
 {
     public:
         float* data;
+        double* geotransform; //!< Six double buffer for storing the affine transformations
+        // https://gdal.org/api/gdaldataset_cpp.html#_CPPv4N11GDALDataset15GetGeoTransformEPd
         int width; //!< width of image in pixels
         int height; //!< height of image in pixels
         float min; //!< Minimum value in image
@@ -17,6 +19,7 @@ class HeightMap
 
         HeightMap(const char* filepath);
         float get_pixel(int x,int y);
+        void blur(float standard_deviation);
 };
 
 /**
@@ -26,10 +29,24 @@ class CellMap
 {
     public:
         uint8_t* cells; //!< pointer to the first cell in the array. uint8_t is a 8 bit unsigned integer 
+        double* geotransform; //!< Six double buffer for storing the affine transformations
         int width; //!< width of image in cells
         int height; //!< height of image in cells
         OGRSpatialReference reference_system;
-
+        CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system, double* geotransform);
-        CellMap(int width, int height, uint8_t* cells, OGRSpatialReference reference_system);
         int get_cell(int x,int y);
 };
+
+class Point
+{
+    public:
+        int x;
+        int y;
+        int mscase;
+        bool allocated = false;
+        Point(int x, int y, int mscase){
+            this -> x = x;
+            this -> y = y;
+            this -> mscase = mscase;
+        };
+};

src/main.cpp

@@ -1,60 +1,191 @@
 #include "contour_creator.hh"
+#include <algorithm>
+#include <array>
 #include <gdal/ogr_api.h>
 #include "ogrsf_frmts.h"
+#include <gdal/ogr_core.h>
+#include <gdal/ogr_feature.h>
+#include <gdal/ogr_geometry.h>
 #include <iostream>
 #include <ostream>
+#include <sys/types.h>
+#include <tuple>
 #include <vector>
 #include <cstdint>
 #include <gdal/gdal.h>
 #include "gdal/gdal_priv.h"
 #include <gdal/gdal_frmts.h>
 #include <omp.h>
+#include "argh.h"
 
-CellMap produce_cellmap(HeightMap* heightmap, float z)
+std::vector<Point> produce_cellmap(HeightMap* heightmap, float z)
 {
-    int length = (heightmap->width-1)*(heightmap->height-1);
+    int length = (heightmap->width-1)*(heightmap->height-1); // Defining length of vector
     uint8_t *cells = (uint8_t *) CPLMalloc(sizeof(uint8_t)*length);
+    std::vector<Point> points; // Initiating a vector of points
     for (int i = 0; i<length; i++) {
         int y = i/(heightmap->width-1);
         int x = i%(heightmap->width-1);
-        uint8_t result = (heightmap->get_pixel(x,y)>z) + 
+        uint8_t result = (heightmap->get_pixel(x,y)>z)*8 + 
-        (heightmap->get_pixel(x+1,y)>z)*2 +
+        (heightmap->get_pixel(x+1,y)>z)*4 +
-        (heightmap->get_pixel(x+1,y+1)>z)*4 +
+        (heightmap->get_pixel(x+1,y+1)>z)*2 +
-        (heightmap->get_pixel(x,y+1)>z)*8;
+        (heightmap->get_pixel(x,y+1)>z);
-        *(cells + i) = result;
+        if (result != 0 and result != 15 ) {
+            points.push_back(Point(x, y, result));
+        };
     }
-    
+    return points;
-    return CellMap(heightmap->width-1, heightmap->height-1, cells, heightmap->reference_system);
-
-
 }
 
-std::vector<CellMap> vector_cellmap(HeightMap* heightmap, int interval)
+bool is_in(int value, std::array<int, 8> array)
+{
+    return std::binary_search(array.begin(), array.end(), value);
+}
+
+std::vector<std::vector<Point>> create_lines(HeightMap* heightmap, int interval)
 {
     int num_contours = (heightmap->max - heightmap->min)/interval;
-    std::vector<CellMap> vector_contours;
+    std::vector<std::vector<Point>> vector_contours;
-    omp_set_num_threads(12);
+
+    std::array<int,8> north{4, 5,6,7,8,9,10,11};
+    std::array<int,8> south{1,2,5,6,9,10,13,14};
+    std::array<int,8> east{2,3,4,5,10,11,12,13};
+    std::array<int,8> west{1,3,5,7,8,10,12,14};
 
     #pragma omp parallel
     {
-        std::vector<CellMap> vec_private;
+        std::vector<std::vector<Point>> vec_private;
         #pragma omp for
         for (int i = 1; i <= num_contours; i++) 
         {
-            vec_private.push_back(produce_cellmap(heightmap, heightmap->min + interval*i));
+            auto points = produce_cellmap(heightmap, heightmap->min + interval*i);
-            std::cout << "Execute thread " << omp_get_num_threads() << " ";
+            std::vector<Point> line;
+            int x;
+            int y;
+            int points_allocated = 0;
+            x = points[0].x;
+            y = points[0].y;
+            int current_case = points[0].mscase;
+            points[0].allocated = true;
+            line.push_back(points[0]);
+            for (int j = 0; j < points.size(); j++){
+                for (int k = 0; k < points.size(); k++){
+                    Point* candidate = &points[k];
+                    if (!candidate->allocated) {
+                        if (candidate->x == x +1 && candidate->y == y /*&& is_in(current_case, east)*/) {
+                            candidate->allocated = true;
+                            x = candidate->x;
+                            y = candidate->y;
+                            current_case = candidate->mscase;
+                            line.push_back(*candidate);
+                            points_allocated++;
+                            break;
+                        }
+                        else if (candidate->x == x -1 && candidate->y == y /*&& is_in(current_case, west)*/) {
+                            x = candidate->x;
+                            y = candidate->y;
+                            current_case = candidate->mscase;
+                            candidate->allocated = true;
+                            line.push_back(*candidate);
+                            points_allocated++;
+                            break;
+                        }
+                        else if (candidate->x == x && candidate->y == y + 1 /*&& is_in(current_case, north)*/) {
+                            x = candidate->x;
+                            y = candidate->y;
+                            current_case = candidate->mscase;
+                            candidate->allocated = true;
+                            line.push_back(*candidate);
+                            points_allocated++;
+                            break;
+                        }
+                        else if (candidate->x == x && candidate->y == y - 1 /*&& is_in(current_case, south)*/) {
+                            x = candidate->x;
+                            y = candidate->y;
+                            current_case = candidate->mscase;
+                            candidate->allocated = true;
+                            line.push_back(*candidate);
+                            points_allocated++;
+                            break;
+                        }
+                    }
+                }
+                if (j > points_allocated)
+                {
+                    vec_private.push_back(line);
+                    line.clear();
+                    //std::cout << points.size() << " ";
+                    for (int k = 0; k < points.size(); k++){
+                        Point* candidate = &points[k];
+                        if (!candidate->allocated)
+                        {
+                            line.push_back(*candidate);
+                            candidate->allocated = true;
+                            x = candidate->x;
+                            y = candidate->y;
+                            current_case = candidate->mscase;
+                            points_allocated++;
+                            break;
+                        }
+                    }
+                }
+            }
         }
         #pragma omp critical
-        vector_contours.insert(vector_contours.end(), vec_private.begin(), vec_private.end());
 
+        vector_contours.insert(vector_contours.end(), vec_private.begin(), vec_private.end());
     }
         return vector_contours;
 }
 
-void write_output_file(std::vector<CellMap> cellmaps, const char *filepath)
+std::tuple<double, double> local_to_projected(double* geotransform, int x, int y)
+{
+    return {
+        geotransform[1] * x + geotransform[2] * y +
+        geotransform[1] * 0.5 + geotransform[2] * 0.5 + geotransform[0],
+        geotransform[4] * x + geotransform[5] * y +
+        geotransform[4] * 0.5 + geotransform[5] * 0.5 + geotransform[3]
+        };
+}
+
+constexpr std::tuple<double, double, double, double> marching_squares_lookup(int mcase)
+{
+    switch (mcase) {
+        case 1:
+            return {0, 0.5, 0.5, 0};
+        case 14:
+            return {0.5, 0, 0, 0.5};
+        case 2: 
+            return {0.5, 0, 1, 0.5};
+        case 13:
+            return {1, 0.5, 0.5, 0};
+        case 3:
+            return {0, 0.5, 1, 0.5};
+        case 12:
+            return {1, 0.5, 0, 0.5};
+        case 4:
+            return {0.5, 1, 1, 0.5};
+        case 11:
+            return {1, 0.5, 0.5, 1};
+        case 5: 
+        case 10:
+            return {0, 0, 0, 0}; //FIXME
+        case 6: 
+            return {0.5, 1, 0.5, 0};
+        case 9:
+            return {0.5, 0, 0.5, 1};
+        case 7: 
+            return {0, 0.5, 0.5, 1};
+        case 8:
+            return {0.5, 0.5, 0, 0.5};
+
+        
+    };
+}
+void write_output_file(std::vector<std::vector<Point>> all_points, const char *filepath, HeightMap* heightmap)
 {
 
-    const char *pszDriverName = "ESRI Shapefile";
+    const char *pszDriverName = "GeoJSON";
     GDALDriver *poDriver;
 
     GDALAllRegister();
@@ -77,7 +208,7 @@ void write_output_file(std::vector<CellMap> cellmaps, const char *filepath)
 
     OGRLayer *poLayer;
 
-    poLayer = poDS->CreateLayer( "contours", &(cellmaps[0]).reference_system, wkbPoint, NULL );
+    poLayer = poDS->CreateLayer( "contours", &heightmap->reference_system, wkbLineString, NULL );
     if( poLayer == NULL )
     {
         printf( "Layer creation failed.\n" );
@@ -93,81 +224,101 @@ void write_output_file(std::vector<CellMap> cellmaps, const char *filepath)
         printf( "Creating Name field failed.\n" );
         exit( 1 );
     }
-    for (int j = 0; j < cellmaps.size(); j++)
+    int width = heightmap->width -1;
+    int height = heightmap->height -1;
+    int size = width * height;
+
+    for (int j = 0; j < all_points.size(); j++)
     {
-        CellMap* cellmap = &cellmaps[j];
+        OGRFeature *feature;
-        for (int i = 0; i<cellmap->height*cellmap->width; i++) {
+        OGRLineString *geometry = new OGRLineString();
-            if (*(cellmap->cells + i) != 0 && *(cellmap->cells + i) != 15)
+        feature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
+        feature->SetField( "Name", j );
+
+        std::vector<Point> points = all_points[j];
+        
+        bool left_to_right = true;
+        
+        if ((points[0].y - points[1].y) < 0)
+            left_to_right = true;
+        else
+            left_to_right = false;
+
+        for (int k = 0; k < points.size(); k++)
         {
-                int x_int = i%cellmap->width;
+            bool left_to_right = true;
-                int y_int = cellmap->height*cellmap->width - i/cellmap->width;
-                double x = double(x_int);
-                double y = double(y_int);
         
-                OGRFeature *poFeature;
+            if ((points[k].x - points[k-1].x) < 0)
+                left_to_right = false;
+            bool top_to_bottom = true;
+            if ((points[k].y - points[k-1].y) < 0)
+                left_to_right = false;
 
-                poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() );
-                poFeature->SetField( "Name", *(cellmap->cells + i) );
 
-                //OGRSpatialReference local;
 
-                //auto poCT = OGRCreateCoordinateTransformation( &local, &cellmap->reference_system );
+            int x_raw = points[k].x;
+            int y_raw = points[k].y;
+            auto [x, y] = local_to_projected(heightmap->geotransform, x_raw, y_raw);
 
-                /*
+            auto [x1, y1, x2, y2] = marching_squares_lookup(points[k].mscase);
-                if( poCT == NULL || !poCT->Transform( 1, &x, &y ) )
-                    printf( "Transformation failed.\n" );
-                */
-                OGRPoint pt;
-                pt.setX( x );
-                pt.setY( y );
             
-                poFeature->SetGeometry( &pt );
+            if (left_to_right or !top_to_bottom)
+            {
+                geometry->setPoint(geometry->getNumPoints(),x + (x1/6) ,y + (y1/6));
+                geometry->setPoint(geometry->getNumPoints(),x + (x2/6) ,y + (y2/6));
+            }
+            else {
+                geometry->setPoint(geometry->getNumPoints(),x + (x2/6) ,y + (y2/6));
+                geometry->setPoint(geometry->getNumPoints(),x + (x1/6) ,y + (y1/6));
+            }
+        }
 
-                if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE )
+        if ( feature->SetGeometry(geometry) != OGRERR_NONE)
+        {
+            printf( "Failed to set geometry.\n" );
+            exit( 1 );
+        }
+        OGRGeometryFactory::destroyGeometry(geometry);
+        if( poLayer->CreateFeature( feature ) != OGRERR_NONE )
         {
             printf( "Failed to create feature in shapefile.\n" );
             exit( 1 );
         }
-
+        OGRFeature::DestroyFeature( feature );
-                OGRFeature::DestroyFeature( poFeature );
     }
-        }
-    }
-
     GDALClose( poDS );
-    /*
-    OGRDataSourceH datasource = OGR_Dr_CreateDataSource(OGRGetDriverByName("GeoJSON"), "contour.geojson", new char*);
-    OGRSpatialReference reference_system = cellmap->reference_system;
-    OGRLayerH layer = OGR_DS_CreateLayer(datasource, "contour", &reference_system, wkbLineString25D, new char*);
-    
-    auto feature = OGR_F_Create(OGR_L_GetLayerDefn(static_cast<OGRLayerH>(layer)));
-
-    OGRGeometryH geometry = OGR_G_CreateGeometry(wkbLineString25D);
-    */
 }
 
 int main(int argc, const char* argv[])
 {
-    const char* filepath = argv[1];
+    argh::parser cmdl(argv);
-    HeightMap map(filepath); 
 
-    std::cout << "x: " << map.width << "    y: " << map.height << "\n";
+    if (cmdl[{ "-h", "--help" }])
-    std::cout << "max: " << map.max << "    min: " << map.min << "\n";
-    
-    auto cellmap = produce_cellmap(&map, 40);
-    /*
-    for (int y = 0; y < cellmap.height; y++)
     {
-        for (int x = 0; x < cellmap.width; x++)
+        std::cout << "Usage:\n"
-        {
+            << "contour_creator [OPTIONS] <input_file>\n"
-            if (cellmap.get_cell(x, y) && cellmap.get_cell(x, y)!=15) {
+            << "-o; --output <FILENAME.geojson> - File to write output to (Default: contours.geojson)\n"
-                std::cout << x << ","<< y << "=" << cellmap.get_cell(x, y) << " ";
+            << "-i; --interval <int> - Set the interval between contours (Default: 5)\n"
+            << "-b; --blur - Blur the image\n"
+            << "--stats - Print statistical information about the heightmap\n";
+            exit(0);
     }
-        }
+    int interval;
-       //std::cout << "\n";
+    cmdl({"-i", "--interval"}, 5) >> interval;
-    }
+    if (interval <= 0)
-    */
+        std::cerr << "Interval must be valid positive integer!" << "\n";
     
-    auto cellmaps = vector_cellmap(&map, 5);
+    std::string output_file;
-    write_output_file(cellmaps, "out.shp");
+    cmdl({"-o", "--output"}, "contours.geojson") >> output_file;
+
+    std::string input_file;
+    cmdl(1) >> input_file;
+    
+
+    HeightMap map(input_file.c_str());
+    if (cmdl[{"-b", "--blur"}])
+        map.blur(0.8);
+    auto lines = create_lines(&map, interval);
+    write_output_file(lines, output_file.c_str(), &map);
+    std::cout << "Contours written to " << output_file << " 🗺️\n";
 }