diff --git a/.gitignore b/.gitignore
index eb3602d..6134f93 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,3 +16,4 @@ coverage_output
 .*.swp
 .project
 build/
+perceptualdiff.dSYM/
diff --git a/compare_args.cpp b/compare_args.cpp
index 42956d8..0b0250b 100644
--- a/compare_args.cpp
+++ b/compare_args.cpp
@@ -102,7 +102,7 @@ static void print_help()
 #endif
 }
 
-bool CompareArgs::parse_args(int argc, char **argv)
+bool CompareArgs::parse_args(const int argc, const char **argv)
 {
     if (argc <= 1)
     {
diff --git a/compare_args.h b/compare_args.h
index 9a94b65..88de6cc 100644
--- a/compare_args.h
+++ b/compare_args.h
@@ -35,7 +35,7 @@ class CompareArgs
 
     CompareArgs();
 
-    bool parse_args(int argc, char **argv);
+    bool parse_args(const int argc, const char **argv);
     void print_args() const;
 
     std::shared_ptr<RGBAImage> image_a_;
diff --git a/lpyramid.cpp b/lpyramid.cpp
index cbbbce1..e9e9e14 100644
--- a/lpyramid.cpp
+++ b/lpyramid.cpp
@@ -22,22 +22,7 @@ Place, Suite 330, Boston, MA 02111-1307 USA
 #include <cassert>
 
 
-static std::vector<float> copy(const float *img,
-                               const unsigned int width,
-                               const unsigned int height)
-{
-    const auto max = width * height;
-    std::vector<float> out(max);
-    for (auto i = 0u; i < max; i++)
-    {
-        out[i] = img[i];
-    }
-
-    return out;
-}
-
-
-LPyramid::LPyramid(const float *image,
+LPyramid::LPyramid(const std::vector<float> &image,
                    const unsigned int width, const unsigned int height)
     : width_(width), weight_(height)
 {
@@ -47,7 +32,7 @@ LPyramid::LPyramid(const float *image,
     {
         if (i == 0 or width * height <= 1)
         {
-            levels_[i] = copy(image, width, height);
+            levels_[i] = image;
         }
         else
         {
@@ -65,27 +50,21 @@ void LPyramid::convolve(std::vector<float> &a,
     assert(b.size() > 1);
 
     const float Kernel[] = {0.05f, 0.25f, 0.4f, 0.25f, 0.05f};
-    #pragma omp parallel for
+    #pragma omp parallel for shared(a, b, Kernel)
     for (auto y = 0u; y < weight_; y++)
     {
         for (auto x = 0u; x < width_; x++)
         {
-            auto index = y * width_ + x;
-            a[index] = 0.0f;
+            const auto index = y * width_ + x;
+            auto a_index = 0.0f;
             for (auto i = -2; i <= 2; i++)
             {
                 for (auto j = -2; j <= 2; j++)
                 {
                     int nx = x + i;
                     int ny = y + j;
-                    if (nx < 0)
-                    {
-                        nx = -nx;
-                    }
-                    if (ny < 0)
-                    {
-                        ny = -ny;
-                    }
+                    nx = std::max(nx, -nx);
+                    ny = std::max(ny, -ny);
                     if (nx >= static_cast<long>(width_))
                     {
                         nx = 2 * width_ - nx - 1;
@@ -94,10 +73,11 @@ void LPyramid::convolve(std::vector<float> &a,
                     {
                         ny = 2 * weight_ - ny - 1;
                     }
-                    a[index] +=
+                    a_index +=
                         Kernel[i + 2] * Kernel[j + 2] * b[ny * width_ + nx];
                 }
             }
+            a[index] = a_index;
         }
     }
 }
diff --git a/lpyramid.h b/lpyramid.h
index 6796269..d4ae2a5 100644
--- a/lpyramid.h
+++ b/lpyramid.h
@@ -29,9 +29,9 @@ class LPyramid
 {
 public:
 
-    LPyramid(const float *image, unsigned int width, unsigned int height);
+    LPyramid(const std::vector<float> &image, const unsigned int width, const unsigned int height);
 
-    float get_value(unsigned int x, unsigned int y, unsigned int level) const;
+    float get_value(unsigned int x, const unsigned int y, const unsigned int level) const;
 
 private:
 
diff --git a/metric.cpp b/metric.cpp
index 7842765..4480ce7 100644
--- a/metric.cpp
+++ b/metric.cpp
@@ -51,7 +51,7 @@ static constexpr float to_degrees(const float radians)
 //
 // Returns the threshold luminance given the adaptation luminance.
 // Units are candelas per meter squared.
-static float tvi(float adaptation_luminance)
+static float tvi(const float adaptation_luminance)
 {
     const auto log_a = log10f(adaptation_luminance);
 
@@ -83,7 +83,7 @@ static float tvi(float adaptation_luminance)
 
 // computes the contrast sensitivity function (Barten SPIE 1989)
 // given the cycles per degree (cpd) and luminance (lum)
-static float csf(float cpd, float lum)
+static float csf(const float cpd, const float lum)
 {
     const auto a = 440.f * powf((1.f + 0.7f / lum), -0.2f);
     const auto b = 0.3f * powf((1.0f + 100.0f / lum), 0.15f);
@@ -96,7 +96,7 @@ static float csf(float cpd, float lum)
 * Visual Masking Function
 * from Daly 1993
 */
-static float mask(float contrast)
+static float mask(const float contrast)
 {
     const auto a = powf(392.498f * contrast, 0.7f);
     const auto b = powf(0.0153f * a, 4.f);
@@ -105,12 +105,12 @@ static float mask(float contrast)
 
 
 // convert Adobe RGB (1998) with reference white D65 to XYZ
-static void adobe_rgb_to_xyz(float r, float g, float b,
+static void adobe_rgb_to_xyz(const float r, const float g, const float b,
                              float &x, float &y, float &z)
 {
     // matrix is from http://www.brucelindbloom.com/
-    x = r * 0.576700f + g * 0.185556f + b * 0.188212f;
-    y = r * 0.297361f + g * 0.627355f + b * 0.0752847f;
+    x = r * 0.576700f  + g * 0.185556f  + b * 0.188212f;
+    y = r * 0.297361f  + g * 0.627355f  + b * 0.0752847f;
     z = r * 0.0270328f + g * 0.0706879f + b * 0.991248f;
 }
 
@@ -131,16 +131,17 @@ struct White
 static const White global_white;
 
 
-static void xyz_to_lab(float x, float y, float z, float &L, float &A, float &B)
+static void xyz_to_lab(const float x, const float y, const float z, float &L, float &A, float &B)
 {
     const float epsilon = 216.0f / 24389.0f;
     const float kappa = 24389.0f / 27.0f;
+    const float r[] = {
+        x / global_white.x,
+        y / global_white.y,
+        z / global_white.z
+    };
     float f[3];
-    float r[3];
-    r[0] = x / global_white.x;
-    r[1] = y / global_white.y;
-    r[2] = z / global_white.z;
-    for (unsigned int i = 0; i < 3; i++)
+    for (auto i = 0u; i < 3; i++)
     {
         if (r[i] > epsilon)
         {
@@ -157,7 +158,7 @@ static void xyz_to_lab(float x, float y, float z, float &L, float &A, float &B)
 }
 
 
-static unsigned int adaptation(float num_one_degree_pixels)
+static unsigned int adaptation(const float num_one_degree_pixels)
 {
     auto num_pixels = 1.f;
     auto adaptation_level = 0u;
@@ -176,7 +177,7 @@ static unsigned int adaptation(float num_one_degree_pixels)
 
 bool yee_compare(CompareArgs &args)
 {
-    if ((args.image_a_->get_width() != args.image_b_->get_width()) or
+    if ((args.image_a_->get_width()  != args.image_b_->get_width()) or
         (args.image_a_->get_height() != args.image_b_->get_height()))
     {
         args.error_string_ = "Image dimensions do not match\n";
@@ -203,12 +204,6 @@ bool yee_compare(CompareArgs &args)
     }
 
     // Assuming colorspaces are in Adobe RGB (1998) convert to XYZ.
-    std::vector<float> a_x(dim);
-    std::vector<float> a_y(dim);
-    std::vector<float> a_z(dim);
-    std::vector<float> b_x(dim);
-    std::vector<float> b_y(dim);
-    std::vector<float> b_z(dim);
     std::vector<float> a_lum(dim);
     std::vector<float> b_lum(dim);
 
@@ -222,25 +217,34 @@ bool yee_compare(CompareArgs &args)
         std::cout << "Converting RGB to XYZ\n";
     }
 
-    #pragma omp parallel for
+    const auto gamma = args.gamma_;
+    const auto luminance = args.luminance_;
+
+    #pragma omp parallel for shared(args, a_lum, b_lum, a_a, a_b, b_a, b_b)
     for (auto y = 0u; y < h; y++)
     {
         for (auto x = 0u; x < w; x++)
         {
             const auto i = x + y * w;
-            auto r = powf(args.image_a_->get_red(i) / 255.0f, args.gamma_);
-            auto g = powf(args.image_a_->get_green(i) / 255.0f, args.gamma_);
-            auto b = powf(args.image_a_->get_blue(i) / 255.0f, args.gamma_);
-            adobe_rgb_to_xyz(r, g, b, a_x[i], a_y[i], a_z[i]);
+            const auto a_color_r = powf(args.image_a_->get_red(i) / 255.0f, gamma);
+            const auto a_color_g = powf(args.image_a_->get_green(i) / 255.0f, gamma);
+            const auto a_color_b = powf(args.image_a_->get_blue(i) / 255.0f, gamma);
+            float a_x;
+            float a_y;
+            float a_z;
+            adobe_rgb_to_xyz(a_color_r, a_color_g, a_color_b, a_x, a_y, a_z);
             float l;
-            xyz_to_lab(a_x[i], a_y[i], a_z[i], l, a_a[i], a_b[i]);
-            r = powf(args.image_b_->get_red(i) / 255.0f, args.gamma_);
-            g = powf(args.image_b_->get_green(i) / 255.0f, args.gamma_);
-            b = powf(args.image_b_->get_blue(i) / 255.0f, args.gamma_);
-            adobe_rgb_to_xyz(r, g, b, b_x[i], b_y[i], b_z[i]);
-            xyz_to_lab(b_x[i], b_y[i], b_z[i], l, b_a[i], b_b[i]);
-            a_lum[i] = a_y[i] * args.luminance_;
-            b_lum[i] = b_y[i] * args.luminance_;
+            xyz_to_lab(a_x, a_y, a_z, l, a_a[i], a_b[i]);
+            const auto b_color_r = powf(args.image_b_->get_red(i) / 255.0f, gamma);
+            const auto b_color_g = powf(args.image_b_->get_green(i) / 255.0f, gamma);
+            const auto b_color_b = powf(args.image_b_->get_blue(i) / 255.0f, gamma);
+            float b_x;
+            float b_y;
+            float b_z;
+            adobe_rgb_to_xyz(b_color_r, b_color_g, b_color_b, b_x, b_y, b_z);
+            xyz_to_lab(b_x, b_y, b_z, l, b_a[i], b_b[i]);
+            a_lum[i] = a_y * luminance;
+            b_lum[i] = b_y * luminance;
         }
     }
 
@@ -249,8 +253,8 @@ bool yee_compare(CompareArgs &args)
         std::cout << "Constructing Laplacian Pyramids\n";
     }
 
-    const LPyramid la(&a_lum[0], w, h);
-    const LPyramid lb(&b_lum[0], w, h);
+    const LPyramid la(a_lum, w, h);
+    const LPyramid lb(b_lum, w, h);
 
     const auto num_one_degree_pixels =
         to_degrees(2 *
@@ -284,60 +288,35 @@ bool yee_compare(CompareArgs &args)
     auto pixels_failed = 0u;
     auto error_sum = 0.;
 
-    #pragma omp parallel for reduction(+ : pixels_failed) reduction(+ : error_sum)
+    #pragma omp parallel for reduction(+ : pixels_failed, error_sum) shared(args, a_a, a_b, b_a, b_b, cpd, F_freq)
     for (auto y = 0u; y < h; y++)
     {
         for (auto x = 0u; x < w; x++)
         {
-            const auto index = x + y * w;
-            float contrast[MAX_PYR_LEVELS - 2];
-            float sum_contrast = 0;
+            const auto index = y * w + x;
+            const auto adapt = std::max((la.get_value(x, y, adaptation_level) +
+                                         lb.get_value(x, y, adaptation_level)) * 0.5f,
+                                        1e-5f);
+            auto sum_contrast = 0.f;
+            auto factor = 0.f;
             for (auto i = 0u; i < MAX_PYR_LEVELS - 2; i++)
             {
-                auto n1 =
+                const auto n1 =
                     fabsf(la.get_value(x, y, i) - la.get_value(x, y, i + 1));
-                auto n2 =
+                const auto n2 =
                     fabsf(lb.get_value(x, y, i) - lb.get_value(x, y, i + 1));
-                auto numerator = (n1 > n2) ? n1 : n2;
-                auto d1 = fabsf(la.get_value(x, y, i + 2));
-                auto d2 = fabsf(lb.get_value(x, y, i + 2));
-                auto denominator = (d1 > d2) ? d1 : d2;
-                if (denominator < 1e-5f)
-                {
-                    denominator = 1e-5f;
-                }
-                contrast[i] = numerator / denominator;
-                sum_contrast += contrast[i];
-            }
-            if (sum_contrast < 1e-5)
-            {
-                sum_contrast = 1e-5f;
-            }
-            float F_mask[MAX_PYR_LEVELS - 2];
-            auto adapt = la.get_value(x, y, adaptation_level) +
-                         lb.get_value(x, y, adaptation_level);
-            adapt *= 0.5f;
-            if (adapt < 1e-5)
-            {
-                adapt = 1e-5f;
-            }
-            for (auto i = 0u; i < MAX_PYR_LEVELS - 2; i++)
-            {
-                F_mask[i] = mask(contrast[i] * csf(cpd[i], adapt));
-            }
-            auto factor = 0.f;
-            for (auto i = 0u; i < MAX_PYR_LEVELS - 2; i++)
-            {
-                factor += contrast[i] * F_freq[i] * F_mask[i] / sum_contrast;
-            }
-            if (factor < 1)
-            {
-                factor = 1;
-            }
-            if (factor > 10)
-            {
-                factor = 10;
+                const auto numerator = std::max(n1, n2);
+                const auto d1 = fabsf(la.get_value(x, y, i + 2));
+                const auto d2 = fabsf(lb.get_value(x, y, i + 2));
+                const auto denominator = std::max(std::max(d1, d2), 1e-5f);
+                const auto contrast = numerator / denominator;
+                const auto F_mask = mask(contrast * csf(cpd[i], adapt));
+                factor += contrast * F_freq[i] * F_mask;
+                sum_contrast += contrast;
             }
+            sum_contrast = std::max(sum_contrast, 1e-5f);
+            factor /= sum_contrast;
+            factor = std::min(std::max(factor, 1.f), 10.f);
             const auto delta =
                 fabsf(la.get_value(x, y, 0) - lb.get_value(x, y, 0));
             error_sum += delta;
@@ -359,11 +338,9 @@ bool yee_compare(CompareArgs &args)
                     // Don't do color test at all.
                     color_scale = 0.0;
                 }
-                auto da = a_a[index] - b_a[index];
-                auto db = a_b[index] - b_b[index];
-                da = da * da;
-                db = db * db;
-                const auto delta_e = (da + db) * color_scale;
+                const auto da = a_a[index] - b_a[index];
+                const auto db = a_b[index] - b_b[index];
+                const auto delta_e = (da * da + db * db) * color_scale;
                 error_sum += delta_e;
                 if (delta_e > factor)
                 {
diff --git a/perceptualdiff.cpp b/perceptualdiff.cpp
index 1fd8851..ca5b432 100644
--- a/perceptualdiff.cpp
+++ b/perceptualdiff.cpp
@@ -30,7 +30,7 @@ Place, Suite 330, Boston, MA 02111-1307 USA
 #include <string>
 
 
-int main(int argc, char **argv)
+int main(const int argc, const char **argv)
 {
     CompareArgs args;
 
diff --git a/rgba_image.h b/rgba_image.h
index 3c7b69d..0b0e343 100644
--- a/rgba_image.h
+++ b/rgba_image.h
@@ -38,33 +38,33 @@ class RGBAImage
 {
 public:
 
-    RGBAImage(unsigned int w, unsigned int h, const std::string &name="")
+    RGBAImage(const unsigned int w, const unsigned int h, const std::string &name="")
         : width_(w), weight_(h), name_(name), data_(w * h)
     {
     }
 
-    unsigned char get_red(unsigned int i) const
+    unsigned char get_red(const unsigned int i) const
     {
         return (data_[i] & 0xff);
     }
 
-    unsigned char get_green(unsigned int i) const
+    unsigned char get_green(const unsigned int i) const
     {
         return ((data_[i] >> 8) & 0xff);
     }
 
-    unsigned char get_blue(unsigned int i) const
+    unsigned char get_blue(const unsigned int i) const
     {
         return ((data_[i] >> 16) & 0xff);
     }
 
-    unsigned char get_alpha(unsigned int i) const
+    unsigned char get_alpha(const unsigned int i) const
     {
         return ((data_[i] >> 24) & 0xff);
     }
 
-    void set(unsigned char r, unsigned char g, unsigned char b,
-             unsigned char a, unsigned int i)
+    void set(const unsigned char r, const unsigned char g, const unsigned char b,
+             const unsigned char a, const unsigned int i)
     {
         data_[i] = r | (g << 8) | (b << 16) | (a << 24);
     }
@@ -79,17 +79,17 @@ class RGBAImage
         return weight_;
     }
 
-    void set(unsigned int x, unsigned int y, unsigned int d)
+    void set(const unsigned int x, const unsigned int y, const unsigned int d)
     {
         data_[x + y * width_] = d;
     }
 
-    unsigned int get(unsigned int x, unsigned int y) const
+    unsigned int get(const unsigned int x, const  unsigned int y) const
     {
         return data_[x + y * width_];
     }
 
-    unsigned int get(unsigned int i) const
+    unsigned int get(const unsigned int i) const
     {
         return data_[i];
     }