From 390232c4f0bcb6ec5a209e9fc92ecd6894ad818b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Rostislav=20L=C3=A1n?= <lanrosta@seznam.cz>
Date: Thu, 23 Feb 2023 12:58:24 +0100
Subject: [PATCH] Added curved fingerprint model generation, refactored most of
 the code.

---
 src/filters.py |  85 ++++++---
 src/main.py    | 485 +++++++++++++++++++++++++++++++------------------
 2 files changed, 368 insertions(+), 202 deletions(-)

diff --git a/src/filters.py b/src/filters.py
index d241f8c..e226822 100644
--- a/src/filters.py
+++ b/src/filters.py
@@ -8,6 +8,7 @@ import numpy as np
 import cv2 as cv
 from skimage import filters as skiflt
 from skimage import restoration as skirest
+from scipy import signal as sig
 
 
 # Parent class for all the filters
@@ -33,8 +34,8 @@ class convolve(filter):
         kernel = np.array(params["kernel"]) if params["kernel"] else np.ones(
             (ksize, ksize), np.float32) / np.sqrt(ksize)
 
-        print("with params: ksize: " +
-              str(ksize) + " kernel: \n" + str(kernel))
+        #print("with params: ksize: " +
+        #      str(ksize) + " kernel: \n" + str(kernel))
         self.img = cv.filter2D(self.img, -1, kernel)
 
 
@@ -56,8 +57,8 @@ class blur(filter):
             anchor = (-1, -1)
         ksize = int(params["ksize"]) if params["ksize"] else 3
 
-        print("with params: ksize: " +
-              str(ksize) + " anchor: " + str(anchor))
+        #print("with params: ksize: " +
+        #      str(ksize) + " anchor: " + str(anchor))
         self.img = cv.blur(self.img, ksize=(ksize, ksize), anchor=anchor)
 
 
@@ -72,13 +73,15 @@ class gaussian(filter):
         sigmaX = float(params["sigmaX"]) if params["sigmaX"] else 0
         sigmaY = float(params["sigmaY"]) if params["sigmaY"] else 0
 
-        print("with params: ksize: " + str(ksize) +
-              " sigmaX: " + str(sigmaX) + " sigmaY: " + str(sigmaY))
+        
+        #print("with params: ksize: " + str(ksize) +
+        #      " sigmaX: " + str(sigmaX) + " sigmaY: " + str(sigmaY))
         self.img = cv.GaussianBlur(self.img, (ksize, ksize), sigmaX, sigmaY)
 
 
 class median(filter):
-    ''' Median blur filter from OpenCV.
+    ''' Median blur filter from scikit-image.
+    Using this over opencv version as that one is limited to 5x5 kernel.
     '''
     def __init__(self, img):
         super().__init__(img)
@@ -86,8 +89,8 @@ class median(filter):
     def apply(self, params):
         ksize = int(params["ksize"]) if params["ksize"] else 3
 
-        print("with params: ksize: " + str(ksize))
-        self.img = cv.medianBlur(np.float32(self.img), ksize)
+        #print("with params: ksize: " + str(ksize))
+        self.img = skiflt.median(self.img, footprint=np.ones((ksize, ksize)))
 
 
 class bilateral(filter):
@@ -102,8 +105,9 @@ class bilateral(filter):
         sigmaColor = int(params["sigmaColor"]) if params["sigmaColor"] else 75
         sigmaSpace = int(params["sigmaSpace"]) if params["sigmaSpace"] else 75
 
-        print("with params: d: " + str(d) + " sigmaColor: " +
-              str(sigmaColor) + " sigmaSpace: " + str(sigmaSpace))
+        #print("with params: d: " + str(d) + " sigmaColor: " +
+        #      str(sigmaColor) + " sigmaSpace: " + str(sigmaSpace))
+        self.img = np.uint8(self.img)
         self.img = cv.bilateralFilter(self.img, d, sigmaColor, sigmaSpace)
 
 
@@ -119,8 +123,8 @@ class denoise(filter):
         sWS = int(params["searchWindowSize"]
                   ) if params["searchWindowSize"] else 21
 
-        print("with params: h: " + str(h) +
-              " tWS: " + str(tWS) + " sWS: " + str(sWS))
+        #print("with params: h: " + str(h) +
+        #      " tWS: " + str(tWS) + " sWS: " + str(sWS))
         self.img = np.uint8(self.img)
         self.img = cv.fastNlMeansDenoising(
             self.img, h, tWS, sWS)
@@ -146,9 +150,9 @@ class denoise_bilateral(filter):
                           ) if params["channelAxis"] else None
         iterations = int(params["iterations"]) if params["iterations"] else 1
 
-        print("with params: sigma_color: " + str(sigmaColor) +
-              " sigma_spatial: " + str(sigmaSpace) + " channel_axis: " +
-              str(channelAxis) + " iterations: " + str(iterations))
+        #print("with params: sigma_color: " + str(sigmaColor) +
+        #      " sigma_spatial: " + str(sigmaSpace) + " channel_axis: " +
+        #      str(channelAxis) + " iterations: " + str(iterations))
 
         for i in range(iterations):
             self.img = skirest.denoise_bilateral(
@@ -172,8 +176,8 @@ class denoise_tv_chambolle(filter):
                           ) if params["channelAxis"] else None
         iterations = int(params["iterations"]) if params["iterations"] else 1
 
-        print("with params: weight: " + str(weight) +
-              " channel_axis: " + str(channelAxis) + " iterations: " + str(iterations))
+        #print("with params: weight: " + str(weight) +
+        #      " channel_axis: " + str(channelAxis) + " iterations: " + str(iterations))
         for i in range(iterations):
             self.img = skirest.denoise_tv_chambolle(
                 self.img, weight=weight, channel_axis=channelAxis)
@@ -190,7 +194,7 @@ class sharpen(filter):
         kernel = np.matrix(params["kernel"]) if params["kernel"] else np.array(
             [[0, -1, 0], [-1, 5, -1], [0, -1, 0]])
 
-        print("with params: kernel: \n" + str(kernel))
+        #print("with params: kernel: \n" + str(kernel))
         self.img = cv.filter2D(self.img, ddepth=-1, kernel=kernel)
 
 
@@ -211,8 +215,8 @@ class unsharp_mask(filter):
         blurred = cv.medianBlur(self.img, ksize)
         lap = cv.Laplacian(blurred, cv.CV_32F)
 
-        print("with params: strength: " +
-              str(strength) + " ksize: " + str(ksize))
+        #print("with params: strength: " +
+        #      str(strength) + " ksize: " + str(ksize))
         self.img = blurred - strength*lap
 
 
@@ -235,8 +239,8 @@ class unsharp_mask_scikit(filter):
                           ) if params["channelAxis"] else None
 
         #self.img = cv.cvtColor(self.img, cv.COLOR_GRAY2RGB)
-        print("with params: radius: " +
-              str(radius) + " amount: " + str(amount))
+        #print("with params: radius: " +
+        #      str(radius) + " amount: " + str(amount))
         self.img = skiflt.unsharp_mask(
             self.img, radius=radius, amount=amount, channel_axis=channelAxis)
         #self.img = cv.cvtColor(self.img, cv.COLOR_RGB2GRAY)
@@ -265,8 +269,39 @@ class morph(filter):
         else:
             anchor = (-1, -1)
 
-        print("with params: kernel: \n" + str(kernel) + " anchor: " +
-              str(anchor) + " iterations: " + str(iterations) + " op: " + str(op))
+        #print("with params: kernel: \n" + str(kernel) + " anchor: " +
+        #      str(anchor) + " iterations: " + str(iterations) + " op: " + str(op))
         self.img = cv.morphologyEx(
             np.uint8(self.img), op=op, kernel=kernel,
             anchor=anchor, iterations=iterations)
+
+
+class gabor(filter):
+    ''' Gabor filter from OpenCV.
+
+    Performs Gabor filtering on the image.
+    '''
+
+    def __init__(self, img):
+        super().__init__(img)
+
+    # TODO: not working properly
+    def apply(self, params):
+        ksize = int(params["ksize"]) if params["ksize"] else 31
+        sigma = float(params["sigma"]) if params["sigma"] else 10.0
+        theta = params["theta"] if params["theta"] else [0,np.pi/16,np.pi-np.pi/16]
+        lambd = float(params["lambd"]) if params["lambd"] else 10.0
+        gamma = float(params["gamma"]) if params["gamma"] else 0.02
+        psi = float(params["psi"]) if params["psi"] else 0.0
+
+        filters = []
+
+        for i in range(len(theta)):
+            g_kernel = cv.getGaborKernel(ksize=(ksize, ksize), sigma=sigma, theta=theta[i], lambd=lambd, gamma=gamma, psi=psi)
+            g_kernel = g_kernel / 1.5 * g_kernel.sum()
+            filters.append(g_kernel)
+
+        tmp = np.zeros_like(self.img)
+        for i in range(len(filters)):
+            tmp = cv.filter2D(self.img, -1, kernel=filters[i])
+            self.img += np.maximum(self.img, tmp)
diff --git a/src/main.py b/src/main.py
index 623794b..2dd9429 100644
--- a/src/main.py
+++ b/src/main.py
@@ -14,8 +14,8 @@ import numpy as np
 import matplotlib.pyplot as plt
 #from PIL import Image
 import cv2 as cv
-
 from stl import mesh
+import math
 
 # Import custom image filter library
 import filters as flt
@@ -25,96 +25,141 @@ class app:
     def __init__(self):
         # Parse arguments from command line
         self.parse_arguments()
-        self.params = {}
+
+        # List and dict for filters and corresponding parameters
         self.filters = []
-        
-        # Parse configuration from json file
+        self.params = {}
+
+        # Parse configuration from json config file
         if self.args.config:
-            self.config_file = self.args.config[0]
-            self.preset_name = self.args.config[1]
+            self.config_file, self.preset_name = self.args.config
             self.config = json.load(open(self.config_file))
             self.parse_conf()
 
-        # If no config file given, expect filters in command line
-        else:
-            if not self.args.filters:
-                print("No filters given, saving original image")
-
+        elif self.args.filters:
             print("No config file given, using command line arguments")
             i = 0
 
+            # Otherwise expect filters from command line
             for filter in self.args.filters:
                 if filter.find('=') == -1:
-                    # if no '=' in filter, it is a new filter
+                    # if no '=' char in filter, it is a new filter name
                     self.filters.append(filter)
                     i += 1
-                    self.params[i] = {} # create empty dict for params
+                    self.params[i] = {}  # create empty dict for params
                 else:
                     # else it's a parameter for current filter
                     key, value = filter.split('=')
                     self.params[i][key] = value
                 self.parse_params(self.params[i])
 
-        if self.args.stl_file and len(self.args.stl_file) == 3:
-            self.stl_file = self.args.stl_file[0]
-            self.height_line = float(self.args.stl_file[1])
-            self.height_base = float(self.args.stl_file[2])
-            self.mode = "2d"
-
-        elif self.args.stl_file and len(self.args.stl_file) == 2:
-            self.stl_file = self.args.stl_file[0]
-            self.height_line = float(self.args.stl_file[1])
-            self.mode = "3d"
-    
         else:
-            print("No STL file given, saving image only")
-            exit(1)
+            print("No filters given, saving original image")
 
         self.input_file = self.args.input_file
         self.output_file = self.args.output_file
         self.dpi = self.args.dpi
         self.mirror = True if self.args.mirror else False
+
         if exists(self.input_file):
-            self.run()
+            self.run_filtering()
+
         else:
-            print("Input file does not exist", file=sys.stderr)
-            exit(1)
 
-    def run(self):
-        # read as numpy.array
-        self.img = cv.imread(
-            self.input_file, cv.IMREAD_GRAYSCALE).astype(np.uint8)
+            self.error_exit("Input file " + self.input_file +
+                            " does not exist")
 
-        self.width = self.img.shape[1]
-        self.height = self.img.shape[0]
-        self.print_size(self.img.shape)
-        fig = plt.figure(figsize=(self.width/self.dpi, self.height/self.dpi),
-                         frameon=False, dpi=self.dpi)
+        if self.args.stl_file:
 
-        ax = plt.Axes(fig, [0., 0., 1., 1.])
-        ax.set_axis_off()
-        fig.add_axes(ax)
+            # Get stl filename
+            self.stl_file = self.args.stl_file[0]
 
-        if self.mirror is True:
-            self.mirror_image()
+            # Get mode and model parameters
+            if self.args.planar:
+                self.mode = "2d"
 
-        # Apply all filters and save image
-        self.apply_filter()
-        self.save_image(fig, ax)
-        plt.close()
-        if self.args.stl_file:
-            self.make_model()
+                if len(self.args.stl_file) < 3:
+                    self.height_base = 10
+                    self.height_line = 2
+                    print(
+                        "Warning: Too few arguments, using default values (10mm base, 2mm lines)")
+                else:
+                    self.height_line = float(self.args.stl_file[1])
+                    self.height_base = float(self.args.stl_file[2])
+                    print("Base height:", self.height_base,
+                          "mm, lines depth/height:", self.height_line, "mm")
+
+            else:
+                self.mode = "3d"
+
+                if len(self.args.stl_file) < 4:
+                    self.height_line = 2
+                    self.curv_rate_x = 0.5
+                    self.curv_rate_y = 0.5
+                    print(
+                        "Warning: Too few arguments, using default values (2mm lines, curvature 0.5 on x, 0.5 on y)")
+                else:
+                    self.height_line = float(self.args.stl_file[1])
+                    self.curv_rate_x = float(
+                        self.args.stl_file[2])  # finger depth
+                    self.curv_rate_y = float(
+                        self.args.stl_file[3])  # finger depth
+                    # self.curv_rate_x = float(self.args.stl_file[2]) # excentricity x
+                    # self.curv_rate_y = float(self.args.stl_file[3]) # excentricity y
+                    print("Line height:", self.height_line,"mm, x axis curvature:", self.curv_rate_x,
+                         ", y axis curvature:", self.curv_rate_y)
+
+            print(self.mode, "mode selected")
+            self.run_stl()
+
+    def parse_arguments(self):
+        '''Parse arguments from command line using argparse library.
+        '''
+
+        parser = ap.ArgumentParser(prog='main.py',
+                                   description='Program for processing a 2D image into 3D fingerprint.',
+                                   usage='%(prog)s [-h] [-m | --mirror | --no-mirror] [-p] input_file output_file dpi ([-c config_file preset | --config config_file preset] | [filters ...]) [-s stl_file | --stl stl_file height_line height_base | --stl_file stl_file height_line curv_rate_x curv_rate_y]')
+
+        # positional arguments
+        parser.add_argument("input_file", type=str, help="input file path")
+        parser.add_argument("output_file", type=str, help="output file path")
+        parser.add_argument("dpi", type=int, help="dpi of used scanner")
+
+        # boolean switch argument
+        parser.add_argument('-m', '--mirror', type=bool, action=ap.BooleanOptionalAction,
+                            help="switch to mirror input image")
+
+        # another boolean switch argument, this time with value, name of the new file and dimensions
+        parser.add_argument('-s', '--stl_file', type=str, nargs='*',
+                            help="create stl model from processed image")
+
+        # another boolean switch argument, this enables 2d mode
+        parser.add_argument('-p', '--planar', type=bool, action=ap.BooleanOptionalAction,
+                            help="make stl shape planar instead of curved one")
+
+        # configuration file containing presets, preset name
+        # pair argument - give both or none
+        parser.add_argument('-c', '--config', nargs=2,
+                            help='pair: name of the config file with presets, name of the preset')
+
+        # array of unknown length, all filter names saved inside
+        parser.add_argument('filters', type=str, nargs='*',
+                            help="list of filter names and their parameters in form [filter_name1 param1=value param2=value filter_name2 param1=value...]")
+
+        self.args = parser.parse_args()
 
     def parse_params(self, params):
-        ''' Parse parameters of filters.
-        Set to None if not given.
-        They are later set in the filter method.
+        '''Parse parameters of filters. Set to None if not given.
+        They are later set to default values in the filter method apply.
         '''
 
+        # TODO: possibly too bloated, sending all possible params to each filter
         possible_params = {"h", "searchWindowSize", "templateWindowSize",
                            "ksize", "kernel", "sigmaX", "sigmaY",
                            "sigmaColor", "sigmaSpace", "d", "anchor", "iterations",
-                           "op", "strength", "amount", "radius", "weight", "channelAxis"}
+                           "op", "strength", "amount", "radius", "weight", "channelAxis",
+                           "theta", "sigma", "lambda", "gamma", "psi"}
+
         for key in possible_params:
             if params.get(key) is None:
                 params[key] = None
@@ -122,168 +167,157 @@ class app:
                 params[key] = params[key]
 
     def parse_conf(self):
-        ''' Parse configuration file if one was given and store filters with their parameters
+        '''Parse configuration file if one was given. 
+        Store filters and their parameters.
         '''
 
+        # Find preset in config file
         if self.preset_name in self.config:
             filter_array = self.config[self.preset_name]
+            # Iterate over filters in preset, store them and their parameters
             for i, filter in enumerate(range(len(filter_array)), start=1):
                 self.filters.append(filter_array[filter]["name"])
                 self.params[i] = {}
                 for attribute, value in filter_array[filter].items():
+                    # Filter name isn't needed in here
                     if attribute != "name":
                         self.params[i][attribute] = value
                 self.parse_params(self.params[i])
             print("Loaded preset: " + self.preset_name +
                   " from file: " + self.config_file)
         else:
-            print("Preset not found", file=sys.stderr)
+            self.error_exit("Preset not found")
 
-    def parse_arguments(self):
-        ''' Parse arguments from command line
+    def error_exit(self, message):
+        '''Print error message and exit.
         '''
 
-        parser = ap.ArgumentParser(prog='main.py',
-                                   description='Program for processing a 2D image into 3D fingerprint.',
-                                   usage='%(prog)s [-h] [-m | --mirror | --no-mirror] input_file output_file dpi \
-                                   ([-c config_file preset | --config config_file preset] | [filters ...]) \
-                                   [-s stl_file | --stl_file stl_file depth_total depth_line]')
+        print("ERROR:", message, file=sys.stderr)
+        exit(1)
 
-        # positional arguments
-        parser.add_argument("input_file", type=str, help="location with input file")
-        parser.add_argument("output_file", type=str, help="output file location")
-        parser.add_argument("dpi", type=int, help="scanner dpi")
-
-        # boolean switch argument
-        parser.add_argument('-m', "--mirror", help="mirror input image", type=bool, action=ap.BooleanOptionalAction)
+#------------------------- FILTERING -------------------------#
 
-        # another boolean switch argument, this time with value, name of the new file and dimensions
-        parser.add_argument('-s', '--stl_file', type=str, nargs='*',
-                            help="make planar model from processed image", required=False)
+    def run_filtering(self):
+        '''Load from input file, store as numpy.array, 
+        process image using filters and save to output file.
+        '''
 
-        # file with configuration containing presets, new preset name
-        # pair argument - give both or none
-        parser.add_argument('-c', '--config', nargs=2,
-                            help='pair: name of the config file with presets, name of the preset')
+        self.img = cv.imread(
+            self.input_file, cv.IMREAD_GRAYSCALE).astype(np.uint8)
 
-        # array of unknown length, all filter names saved inside
-        parser.add_argument('filters', type=str, nargs='*',
-                            help="list of filter names and their parameters in form [filter_name1 param1=value1 param2=value2 filter_name2 param1=value1...]")
+        self.height, self.width = self.img.shape
+        print("Height: " + str(self.height) + " px and width: "
+              + str(self.width) + " px", file=sys.stderr)
+        size = (self.width/self.dpi, self.height/self.dpi)
+        fig = plt.figure(figsize=size, frameon=False, dpi=self.dpi)
 
-        self.args = parser.parse_args()
+        ax = plt.Axes(fig, [0., 0., 1., 1.])
+        ax.set_axis_off()
+        fig.add_axes(ax)
 
-    def filter_factory(self, filter_name):
-        ''' Selects filter method of filters library.
-        '''
+        if self.mirror is True:
+            self.mirror_image()
 
-        print("Applying " + filter_name + " filter ", end='')
-        return getattr(flt, filter_name)
+        # Apply all filters and save image
+        self.apply_filters()
+        self.save_image(fig, ax)
+        plt.close()
 
     def mirror_image(self):
-        ''' Mirror image when mirroring is needed,
-        should be used only if we want a positive model
-         '''
+        '''Mirror image using opencv, should be used if we want a positive model.
+        '''
 
         print("Mirroring image", file=sys.stderr)
         self.img = cv.flip(self.img, 1)  # 1 for vertical mirror
 
-    def apply_filter(self):
-        ''' Apply filters to image.
-
-        Apply the filters one by one, if none were given, just save original image output.
+    def apply_filters(self):
+        '''Apply filters to image one by one.
+        In case none were given, pass and save original image to the output file.
         '''
 
         if len(self.filters) == 0:
-            # No filter given, just save the image
             pass
         else:
             # Apply all filters
             for i, filter_name in enumerate(self.filters):
-                filter = self.filter_factory(filter_name)
+                # Get filter class from filter.py, use the apply method
+                filter = getattr(flt, filter_name)
                 filter.apply(self, self.params[i+1])
 
     def save_image(self, fig, ax):
-        ''' Save processed image.
-        Colormap set to grayscale to avoid color mismatch.
+        '''Save processed image to the output file.
         '''
 
-        print("Saving image to ", self.output_file, file=sys.stderr)
+        print("Saving image to", self.output_file, file=sys.stderr)
+        # Colormap must be set to grayscale to avoid color mismatch.
         ax.imshow(self.img, cmap="gray")
-        fig.savefig(fname=self.output_file, dpi='figure')
+        fig.savefig(fname=self.output_file, dpi=self.dpi)
 
-    def print_size(self, size):
-        print("Image of height: " + str(size[0]) + 
-            " px and width: " + str(size[1]) + " px", file=sys.stderr)
+#------------------------- STL GENERATION -------------------------#
 
-    def make_model(self):
-        '''After processing image, make a lithophane from it.
+    def run_stl(self):
+        '''Make heightmap, create mesh and save as stl file.
         '''
 
-        print("Making heighthmap", file=sys.stderr)
         self.prepare_heightmap()
 
         if self.mode == "2d":
-            print("Converting to stl format", file=sys.stderr)
             self.make_stl_planar()
-            plt.show()
-            print(f"Saving lithophane to ", self.stl_file, file=sys.stderr)
-            self.save_stl_2d()
+
         elif self.mode == "3d":
-            self.map_image_to_3d()
-            plt.show()
-            self.save_stl_3d()
+            self.make_stl_curved()
+
         else:
-            print("Mode not supported", file=sys.stderr)
-            exit(1)      
+            self.error_exit("Mode not supported")
+
+        plt.show()
+        print(f"Saving model to ", self.stl_file, file=sys.stderr)
+        self.save_stl()
 
     def prepare_heightmap(self):
-        ''' Create numpy meshgrid.
-        Modify image values to get usable depth values.
+        '''Modify image values to get usable height/depth values.
+        Check validity of dimension parameters.
         '''
 
+        # TODO: redo, too complicated, add extra params, redo checks
+
         if self.img.dtype == np.float32 or self.img.dtype == np.float64:
             print("Converting to uint8", file=sys.stderr)
             self.img = self.img * 255
             self.img = self.img.astype(np.uint8)
+            print("Creating mesh", file=sys.stderr)
 
         if self.mode == "2d":
             if self.height_base <= 0:
-                print("Depth of plate height must be positive", file=sys.stderr)
-                exit(1)
+                self.error_exit("Depth of plate height must be positive")
 
             if self.height_line + self.height_base <= 0:
-                print("Line depth must be less than plate thickness", file=sys.stderr)
-                exit(1)
-        
-            print("Base height:", self.height_base,
-                "mm, lines depth/height:", self.height_line, "mm")
+                self.error_exit("Line depth must be less than plate thickness")
 
             # Transform image values to get a heightmap
-            if self.height_line < 0:
-                self.img = (self.height_base + (1 - self.img/255)
-                            * self.height_line)
-            else:
-                self.img = (self.height_base + (1 - self.img/255)
-                            * self.height_line)
+            self.img = (self.height_base + (1 - self.img/255)
+                        * self.height_line)
 
         if self.mode == "3d":
-            #TODO add some checks and print info
-            pass
+            # TODO check curvature values and print info
+            # TODO: copy pasta code, remove
+            # Transform image values to get a heightmap
+            self.img = (1 - self.img/255) * self.height_line
 
-    def add_faces(self, faces, c):
+    def append_faces(self, faces, c):
         # Function to add faces to the list
         faces.append([c, c + 1, c + 2])
         faces.append([c + 1, c + 3, c + 2])
         return c + 4
 
     def make_stl_planar(self):
-        ''' Create mesh from meshgrid.
+        '''Create mesh from meshgrid.
         Create vertices from meshgrid, add depth values from image.
         Create faces from vertices. Add vectors and faces to the model.
 
         From wikipedia.org/wiki/STL_(file_format):
         ascii stl format consists of repeating structures:
-        
+
         facet normal ni nj nk # normal vector
             outer loop
                 vertex v1x v1y v1z # vertex 1
@@ -299,7 +333,7 @@ class app:
 
         self.meshgrid_2d = np.meshgrid(x, y)
 
-        # Add the image matrix to the 2D meshgrid and create 1D array of 3D points
+        # Add the image matrix to the 2D meshgrid and create 1D array of 3D pointsd
         vertex_arr = np.vstack(list(map(np.ravel, self.meshgrid_2d))).T
         z = (self.img / 10).reshape(-1, 1)
         vertex_arr = np.concatenate((vertex_arr, z), axis=1)
@@ -320,42 +354,42 @@ class app:
                 vertices.append([vertex_arr[i+1][j]])
                 vertices.append([vertex_arr[i+1][j+1]])
 
-                count = self.add_faces(faces, count)
+                count = self.append_faces(faces, count)
 
         # Add faces for the backside of the lithophane
         null_arr = np.copy(vertex_arr)
         for i in range(self.height):
             for j in range(self.width):
                 null_arr[i][j][2] = 0
-    
+
         # Back side faces
         for i in range(self.height - 1):
-            for j in range(self.width - 1): 
+            for j in range(self.width - 1):
 
                 vertices.append([null_arr[i][j]])
                 vertices.append([null_arr[i+1][j]])
                 vertices.append([null_arr[i][j+1]])
                 vertices.append([null_arr[i+1][j+1]])
 
-                count = self.add_faces(faces, count)
+                count = self.append_faces(faces, count)
 
         # Horizontal side faces
-        for j in range(self.height - 1):
-            vertices.append([vertex_arr[j][0]])
-            vertices.append([vertex_arr[j+1][0]])
-            vertices.append([null_arr[j][0]])
-            vertices.append([null_arr[j+1][0]])
+        for i in range(self.height - 1):
+            vertices.append([vertex_arr[i][0]])
+            vertices.append([vertex_arr[i+1][0]])
+            vertices.append([null_arr[i][0]])
+            vertices.append([null_arr[i+1][0]])
 
-            count = self.add_faces(faces, count)
+            count = self.append_faces(faces, count)
 
             max = self.width - 1
 
-            vertices.append([vertex_arr[j+1][max]])
-            vertices.append([vertex_arr[j][max]])
-            vertices.append([null_arr[j+1][max]])
-            vertices.append([null_arr[j][max]])
+            vertices.append([vertex_arr[i+1][max]])
+            vertices.append([vertex_arr[i][max]])
+            vertices.append([null_arr[i+1][max]])
+            vertices.append([null_arr[i][max]])
 
-            count = self.add_faces(faces, count)
+            count = self.append_faces(faces, count)
 
         # Vertical side faces
         for j in range(self.width - 1):
@@ -364,7 +398,7 @@ class app:
             vertices.append([null_arr[0][j+1]])
             vertices.append([null_arr[0][j]])
 
-            count = self.add_faces(faces, count)
+            count = self.append_faces(faces, count)
 
             max = self.height - 1
 
@@ -373,41 +407,60 @@ class app:
             vertices.append([null_arr[max][j]])
             vertices.append([null_arr[max][j+1]])
 
-            count = self.add_faces(faces, count)
+            count = self.append_faces(faces, count)
 
         # Convert to numpy arrays
         faces = np.array(faces)
         vertices = np.array(vertices)
 
         # Create the mesh - vertices.shape (no_faces, 3, 3)
-        self.stl_mesh_2d = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
+        self.stl_lithophane = mesh.Mesh(
+            np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
         for i, face in enumerate(faces):
             for j in range(3):
-                self.stl_mesh_2d.vectors[i][j] = vertices[face[j], :]
-
-    def save_stl_2d(self):
-        ''' Save final mesh to stl file.
-        '''
-
-        self.stl_mesh_2d.save(self.stl_file)
+                self.stl_lithophane.vectors[i][j] = vertices[face[j], :]
 
-    def map_image_to_3d(self):
-        ''' Map fingerprint to finger model.
+    def make_stl_curved(self):
+        '''Map fingerprint to finger model.
         '''
 
+        # TODO: if this is the same as 2D, move to heightmap to reduce duplicate code
         x = np.linspace(0, self.width * 25.4 / self.dpi, self.width)
         y = np.linspace(0, self.height * 25.4 / self.dpi, self.height)
 
-        z1 = np.logspace(0, 10, int(np.ceil(self.width / 2)), base=0.7)
+        self.meshgrid_3d = np.meshgrid(x, y)
+
+        # Method 1 - logspace and logarithmic curve
+        '''z1 = np.logspace(0, 10, int(np.ceil(self.width / 2)), base=0.7)
         z2 = np.logspace(10, 0, int(np.floor(self.width / 2)), base=0.7)
         ztemp = 5*np.concatenate((z1, z2))
 
         z = np.array([])
         for i in range(self.height):
-            z = np.concatenate((z, ztemp * pow(np.log(i+2), -1)))
+            z = np.concatenate((z, ztemp + 25*(((i+50)/20)**(-1/2))))
         z = z.reshape(-1, 1)
 
-        self.meshgrid_3d = np.meshgrid(x, y)
+        self.img = (self.img / 10).reshape(-1, 1)
+        z += self.img'''
+
+        # Method 2 - 2 ellipses
+        z = np.array([])
+        for x in range(self.width):
+            z = np.append(z, np.sqrt(1 - (2*x/self.width - 1)**2)
+                            * (self.curv_rate_x**2))
+        z = np.tile(z, (self.height, 1))
+        for y in range(self.height):
+            new = np.sqrt((1 - ((self.height - y)/self.height)**2)
+                          * (self.curv_rate_y**2))
+            z[y] = z[y] + new
+                                  
+        # TODO: clip responsivelly
+        bottom = z[0][math.floor(self.width/2)]
+        #top = self.curv_rate_x**2 + self.curv_rate_y
+        #np.clip(z, bottom, top, out=z)
+        z = z.reshape(-1, 1)
+        self.img = (self.img / 10).reshape(-1, 1)
+        z += self.img
 
         vertex_arr = np.vstack(list(map(np.ravel, self.meshgrid_3d))).T
         vertex_arr = np.concatenate((vertex_arr, z), axis=1)
@@ -416,41 +469,119 @@ class app:
         count = 0
         vertices = []
         faces = []
+        min_point = 0
+        for i in range(self.height - 1):
+            if vertex_arr[i][0][2] <= bottom:
+                min_point = i
+
+        # Add faces for the backside of the lithophane
+        vec_side = (vertex_arr[self.height-1][0][2] -
+                    vertex_arr[min_point][0][2]) / (self.height - min_point)
+        null_arr = np.copy(vertex_arr)
+        for i in range(self.height):
+            for j in range(self.width):
+                null_arr[i][j][2] = 0
+                #null_arr[i][j][2] = bottom + vec_side * i
+                # for smaller mesh
 
         # Iterate over all vertices, create faces
         for i in range(self.height - 1):
             for j in range(self.width - 1):
-
+                if (vertex_arr[i][j][2] <= null_arr[i][j][2]
+                    or vertex_arr[i+1][j][2] <= null_arr[i+1][j][2]
+                    or vertex_arr[i][j+1][2] <= null_arr[i][j+1][2]
+                    or vertex_arr[i+1][j+1][2] <= null_arr[i+1][j+1][2]):
+                    continue
                 vertices.append([vertex_arr[i][j]])
                 vertices.append([vertex_arr[i][j+1]])
                 vertices.append([vertex_arr[i+1][j]])
                 vertices.append([vertex_arr[i+1][j+1]])
 
-                count = self.add_faces(faces, count)
+                count = self.append_faces(faces, count)
+
+        # Rotated back side faces
+        for i in range(self.height - 1):
+            for j in range(self.width - 1):
+                if (vertex_arr[i][j][2] <= null_arr[i][j][2]):
+                    continue
+
+                vertices.append([null_arr[i][j]])
+                vertices.append([null_arr[i+1][j]])
+                vertices.append([null_arr[i][j+1]])
+                vertices.append([null_arr[i+1][j+1]])
+
+                count = self.append_faces(faces, count)
+
+        # Horizontal side faces
+        for i in range(self.height - 1): # right
+            #if (vertex_arr[i][0][2] < null_arr[i][0][2]):
+            #    continue
+
+            vertices.append([vertex_arr[i][0]])
+            vertices.append([vertex_arr[i+1][0]])
+            vertices.append([null_arr[i][0]])
+            vertices.append([null_arr[i+1][0]])
+
+            count = self.append_faces(faces, count)
+
+        for i in range(self.height - 1): # left
+            max = self.width - 1
+            #if (vertex_arr[i][max][2] < null_arr[i][max][2]):
+            #    continue
+
+            vertices.append([vertex_arr[i+1][max]])
+            vertices.append([vertex_arr[i][max]])
+            vertices.append([null_arr[i+1][max]])
+            vertices.append([null_arr[i][max]])
+
+            count = self.append_faces(faces, count)
 
-        #self.finger_base = mesh.Mesh(np.zeros(, dtype=mesh.Mesh.dtype))
+        # Vertical side faces
+        for j in range(self.width - 1): # top
+            #if (vertex_arr[0][j][2] < null_arr[0][j][2]):
+            #    continue
+
+            vertices.append([vertex_arr[0][j+1]])
+            vertices.append([vertex_arr[0][j]])
+            vertices.append([null_arr[0][j+1]])
+            vertices.append([null_arr[0][j]])
+
+            count = self.append_faces(faces, count)
+
+        for j in range(self.width - 1): # bottom
+            max = self.height - 1
+            #if (vertex_arr[max][j][2] < null_arr[max][j][2]):
+            #    continue
 
-        # linear projection
-            # extrude lines in 1 direction 
-        # cylinder / circular projection
-            # extrude lines in direction of a suitable cylinder
-        # normal projection
-            # extrude lines in the direction of normals of given finger model
+            vertices.append([vertex_arr[max][j]])
+            vertices.append([vertex_arr[max][j+1]])
+            vertices.append([null_arr[max][j]])
+            vertices.append([null_arr[max][j+1]])
+
+            count = self.append_faces(faces, count)
 
         # Convert to numpy arrays
         faces = np.array(faces)
         vertices = np.array(vertices)
 
         # Create the mesh - vertices.shape (no_faces, 3, 3)
-        self.mesh_finger = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
+        self.mesh_finger = mesh.Mesh(
+            np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
         for i, face in enumerate(faces):
             for j in range(3):
                 self.mesh_finger.vectors[i][j] = vertices[face[j], :]
 
-    def save_stl_3d(self):
-        ''' Save final mesh to stl file.
+        # print(self.mesh_finger.normals)
+
+    def save_stl(self):
+        '''Save final mesh to stl file.
         '''
 
-        self.mesh_finger.save(self.stl_file)
+        if self.mode == "3d":
+            self.mesh_finger.save(self.stl_file)
+        else:
+            self.stl_lithophane.save(self.stl_file)
+
 
+# run the application
 image = app()