Compositor: Add Anti-Aliasing node

This is an implementation of Enhanced Subpixel Morphological Antialiasing (SMAA)

The algorithm was proposed by:
  Jorge Jimenez, Jose I. Echevarria, Tiago Sousa, Diego Gutierrez

This node provides only SMAA 1x mode, so the operation will be done with no spatial
multisampling nor temporal supersampling. See Patch for comparisons.

The existing AA operation seems to be used only for binary images by some other nodes.
Using SMAA for binary images needs no important parameter such as "threshold", so we
perhaps can switch the operation to SMAA, though that changes existing behavior.

Notes:
1. The program code assumes the screen coordinates are DirectX style that the
   vertical direction is upside-down, so "top" and "bottom" actually represent bottom
   and top, respectively.

Thanks for Habib Gahbiche (zazizizou) to polish and finalize this patch.

Reviewed By: jbakker

Differential Revision: https://developer.blender.org/D2411

4 files changed, 1243 insertions, 0 deletions

diff --git a/extern/CMakeLists.txt b/extern/CMakeLists.txt
index 235c2fa931a..824b2fb0e9c 100644
--- a/extern/CMakeLists.txt
+++ b/extern/CMakeLists.txt
@@ -109,3 +109,7 @@ endif()
 if(WITH_MOD_FLUID)
   add_subdirectory(mantaflow)
 endif()
+
+if (WITH_COMPOSITOR)
+	add_subdirectory(smaa_areatex)
+endif()
diff --git a/extern/smaa_areatex/CMakeLists.txt b/extern/smaa_areatex/CMakeLists.txt
new file mode 100644
index 00000000000..2386b0e7b79
--- /dev/null
+++ b/extern/smaa_areatex/CMakeLists.txt
@@ -0,0 +1,26 @@
+# ***** BEGIN GPL LICENSE BLOCK *****
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+# The Original Code is Copyright (C) 2017, Blender Foundation
+# All rights reserved.
+#
+# The Original Code is: all of this file.
+#
+# Contributor(s): IRIE Shinsuke
+#
+# ***** END GPL LICENSE BLOCK *****
+
+add_executable(smaa_areatex smaa_areatex.cpp)
diff --git a/extern/smaa_areatex/README.blender b/extern/smaa_areatex/README.blender
new file mode 100644
index 00000000000..9c409142ae8
--- /dev/null
+++ b/extern/smaa_areatex/README.blender
@@ -0,0 +1,5 @@
+Project: smaa-cpp
+URL: https://github.com/iRi-E/smaa-cpp
+License: MIT
+Upstream version: 0.4.0
+Local modifications:
diff --git a/extern/smaa_areatex/smaa_areatex.cpp b/extern/smaa_areatex/smaa_areatex.cpp
new file mode 100644
index 00000000000..971706fd64d
--- /dev/null
+++ b/extern/smaa_areatex/smaa_areatex.cpp
@@ -0,0 +1,1208 @@
+/**
+ * Copyright (C) 2016-2017 IRIE Shinsuke
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+/*
+ * smaa_areatex.cpp  version 0.4.0
+ *
+ * This is a part of smaa-cpp that is an implementation of
+ * Enhanced Subpixel Morphological Antialiasing (SMAA) written in C++.
+ *
+ * This program is C++ rewrite of AreaTex.py included in the original
+ * SMAA ditribution:
+ *
+ *   https://github.com/iryoku/smaa/tree/master/Scripts
+ */
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+#include <cmath>
+
+/*------------------------------------------------------------------------------*/
+/* Type Definitions */
+
+class Int2;
+class Dbl2;
+
+class Int2 {
+public:
+	int x, y;
+
+	Int2() { this->x = this->y = 0; }
+	Int2(int x) { this->x = this->y = x; }
+	Int2(int x, int y) { this->x = x; this->y = y; }
+
+	operator Dbl2();
+
+	Int2 operator + (Int2 other) { return Int2(x + other.x, y + other.y); }
+	Int2 operator * (Int2 other) { return Int2(x * other.x, y * other.y); }
+};
+
+class Dbl2 {
+public:
+	double x, y;
+
+	Dbl2() { this->x = this->y = 0.0; }
+	Dbl2(double x) { this->x = this->y = x; }
+	Dbl2(double x, double y) { this->x = x; this->y = y; }
+
+	Dbl2 apply(double (* func)(double)) { return Dbl2(func(x), func(y)); }
+
+	operator Int2();
+
+	Dbl2 operator + (Dbl2 other) { return Dbl2(x + other.x, y + other.y); }
+	Dbl2 operator - (Dbl2 other) { return Dbl2(x - other.x, y - other.y); }
+	Dbl2 operator * (Dbl2 other) { return Dbl2(x * other.x, y * other.y); }
+	Dbl2 operator / (Dbl2 other) { return Dbl2(x / other.x, y / other.y); }
+	Dbl2 operator += (Dbl2 other) { return Dbl2(x += other.x, y += other.y); }
+	bool operator == (Dbl2 other) { return (x == other.x && y == other.y); }
+};
+
+Int2::operator Dbl2() { return Dbl2((double)x, (double)y); }
+Dbl2::operator Int2() { return Int2((int)x, (int)y); }
+
+/*------------------------------------------------------------------------------*/
+/* Data to Calculate Areatex */
+
+/* Texture sizes: */
+/* (it's quite possible that this is not easily configurable) */
+static const int SUBSAMPLES_ORTHO = 7;
+static const int SUBSAMPLES_DIAG  = 5;
+static const int MAX_DIST_ORTHO_COMPAT = 16;
+static const int MAX_DIST_ORTHO = 20;
+static const int MAX_DIST_DIAG  = 20;
+static const int TEX_SIZE_ORTHO = 80; /* 16 * 5 slots = 80 */
+static const int TEX_SIZE_DIAG  = 80; /* 20 * 4 slots = 80 */
+
+/* Number of samples for calculating areas in the diagonal textures: */
+/* (diagonal areas are calculated using brute force sampling) */
+static const int SAMPLES_DIAG = 30;
+
+/* Maximum distance for smoothing u-shapes: */
+static const int SMOOTH_MAX_DISTANCE = 32;
+
+/*------------------------------------------------------------------------------*/
+/* Offset Tables */
+
+/* Offsets for subsample rendering */
+static const double subsample_offsets_ortho[SUBSAMPLES_ORTHO] = {
+	0.0,    /* 0 */
+	-0.25,  /* 1 */
+	0.25,   /* 2 */
+	-0.125, /* 3 */
+	0.125,  /* 4 */
+	-0.375, /* 5 */
+	0.375   /* 6 */
+};
+
+static const Dbl2 subsample_offsets_diag[SUBSAMPLES_DIAG] = {
+	{ 0.00,   0.00},  /* 0 */
+	{ 0.25,  -0.25},  /* 1 */
+	{-0.25,   0.25},  /* 2 */
+	{ 0.125, -0.125}, /* 3 */
+	{-0.125,  0.125}  /* 4 */
+};
+
+/* Mapping offsets for placing each pattern subtexture into its place */
+enum edgesorthoIndices
+{
+	EDGESORTHO_NONE_NONE = 0,
+	EDGESORTHO_NONE_NEGA = 1,
+	EDGESORTHO_NONE_POSI = 2,
+	EDGESORTHO_NONE_BOTH = 3,
+	EDGESORTHO_NEGA_NONE = 4,
+	EDGESORTHO_NEGA_NEGA = 5,
+	EDGESORTHO_NEGA_POSI = 6,
+	EDGESORTHO_NEGA_BOTH = 7,
+	EDGESORTHO_POSI_NONE = 8,
+	EDGESORTHO_POSI_NEGA = 9,
+	EDGESORTHO_POSI_POSI = 10,
+	EDGESORTHO_POSI_BOTH = 11,
+	EDGESORTHO_BOTH_NONE = 12,
+	EDGESORTHO_BOTH_NEGA = 13,
+	EDGESORTHO_BOTH_POSI = 14,
+	EDGESORTHO_BOTH_BOTH = 15,
+};
+
+static const Int2 edgesortho_compat[16] = {
+	{0, 0}, {0, 1}, {0, 3}, {0, 4}, {1, 0}, {1, 1}, {1, 3}, {1, 4},
+	{3, 0}, {3, 1}, {3, 3}, {3, 4}, {4, 0}, {4, 1}, {4, 3}, {4, 4}
+};
+
+static const Int2 edgesortho[16] = {
+	{0, 0}, {0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, {1, 3},
+	{2, 0}, {2, 1}, {2, 2}, {2, 3}, {3, 0}, {3, 1}, {3, 2}, {3, 3}
+};
+
+enum edgesdiagIndices
+{
+	EDGESDIAG_NONE_NONE = 0,
+	EDGESDIAG_NONE_VERT = 1,
+	EDGESDIAG_NONE_HORZ = 2,
+	EDGESDIAG_NONE_BOTH = 3,
+	EDGESDIAG_VERT_NONE = 4,
+	EDGESDIAG_VERT_VERT = 5,
+	EDGESDIAG_VERT_HORZ = 6,
+	EDGESDIAG_VERT_BOTH = 7,
+	EDGESDIAG_HORZ_NONE = 8,
+	EDGESDIAG_HORZ_VERT = 9,
+	EDGESDIAG_HORZ_HORZ = 10,
+	EDGESDIAG_HORZ_BOTH = 11,
+	EDGESDIAG_BOTH_NONE = 12,
+	EDGESDIAG_BOTH_VERT = 13,
+	EDGESDIAG_BOTH_HORZ = 14,
+	EDGESDIAG_BOTH_BOTH = 15,
+};
+
+static const Int2 edgesdiag[16] = {
+	{0, 0}, {0, 1}, {0, 2}, {0, 3}, {1, 0}, {1, 1}, {1, 2}, {1, 3},
+	{2, 0}, {2, 1}, {2, 2}, {2, 3}, {3, 0}, {3, 1}, {3, 2}, {3, 3}
+};
+
+/*------------------------------------------------------------------------------*/
+/* Miscellaneous Utility Functions */
+
+/* Linear interpolation: */
+static Dbl2 lerp(Dbl2 a, Dbl2 b, double p)
+{
+	return a + (b - a) * Dbl2(p);
+}
+
+/* Saturates a value to [0..1] range: */
+static double saturate(double x)
+{
+	return 0.0 < x ? (x < 1.0 ? x : 1.0) : 0.0;
+}
+
+/*------------------------------------------------------------------------------*/
+/* Horizontal/Vertical Areas */
+
+class AreaOrtho {
+	double m_data[SUBSAMPLES_ORTHO][TEX_SIZE_ORTHO][TEX_SIZE_ORTHO][2];
+	bool m_compat;
+	bool m_orig_u;
+public:
+	AreaOrtho(bool compat, bool orig_u) : m_compat(compat), m_orig_u(orig_u) {}
+
+	double *getData() { return (double *)&m_data; }
+	Dbl2 getPixel(int offset_index, Int2 coords) {
+		return Dbl2(m_data[offset_index][coords.y][coords.x][0],
+			    m_data[offset_index][coords.y][coords.x][1]);
+	}
+
+	void areaTex(int offset_index);
+private:
+	void putPixel(int offset_index, Int2 coords, Dbl2 pixel) {
+		m_data[offset_index][coords.y][coords.x][0] = pixel.x;
+		m_data[offset_index][coords.y][coords.x][1] = pixel.y;
+	}
+
+	Dbl2 smoothArea(double d, Dbl2 a1, Dbl2 a2);
+	Dbl2 makeQuad(int x, double d, double o);
+	Dbl2 area(Dbl2 p1, Dbl2 p2, int x);
+	Dbl2 calculate(int pattern, int left, int right, double offset);
+};
+
+/* Smoothing function for small u-patterns: */
+Dbl2 AreaOrtho::smoothArea(double d, Dbl2 a1, Dbl2 a2)
+{
+	Dbl2 b1 = (a1 * Dbl2(2.0)).apply(sqrt) * Dbl2(0.5);
+	Dbl2 b2 = (a2 * Dbl2(2.0)).apply(sqrt) * Dbl2(0.5);
+	double p = saturate(d / (double)SMOOTH_MAX_DISTANCE);
+	return lerp(b1, a1, p) + lerp(b2, a2, p);
+}
+
+/* Smoothing u-patterns by quadratic function: */
+Dbl2 AreaOrtho::makeQuad(int x, double d, double o)
+{
+	double r = (double)x;
+
+	/* fmin() below is a trick to smooth tiny u-patterns: */
+	return Dbl2(r, (1.0 - fmin(4.0, d) * r * (d - r) / (d * d)) * o);
+}
+
+/* Calculates the area under the line p1->p2, for the pixel x..x+1: */
+Dbl2 AreaOrtho::area(Dbl2 p1, Dbl2 p2, int x)
+{
+	Dbl2 d = p2 - p1;
+	double x1 = (double)x;
+	double x2 = x1 + 1.0;
+
+	if ((x1 >= p1.x && x1 < p2.x) || (x2 > p1.x && x2 <= p2.x)) { /* inside? */
+		double y1 = p1.y + (x1 - p1.x) * d.y / d.x;
+		double y2 = p1.y + (x2 - p1.x) * d.y / d.x;
+
+		if ((copysign(1.0, y1) == copysign(1.0, y2) ||
+		     fabs(y1) < 1e-4 || fabs(y2) < 1e-4)) { /* trapezoid? */
+			double a = (y1 + y2) / 2.0;
+			if (a < 0.0)
+				return Dbl2(fabs(a), 0.0);
+			else
+				return Dbl2(0.0, fabs(a));
+		}
+		else { /* Then, we got two triangles: */
+			double x = p1.x - p1.y * d.x / d.y, xi;
+			double a1 = x > p1.x ? y1 * modf(x, &xi) / 2.0 : 0.0;
+			double a2 = x < p2.x ? y2 * (1.0 - modf(x, &xi)) / 2.0 : 0.0;
+			double a = fabs(a1) > fabs(a2) ? a1 : -a2;
+			if (a < 0.0)
+				return Dbl2(fabs(a1), fabs(a2));
+			else
+				return Dbl2(fabs(a2), fabs(a1));
+		}
+	}
+	else
+		return Dbl2(0.0, 0.0);
+}
+
+/* Calculates the area for a given pattern and distances to the left and to the */
+/* right, biased by an offset: */
+Dbl2 AreaOrtho::calculate(int pattern, int left, int right, double offset)
+{
+	Dbl2 a1, a2;
+
+	/*
+	 * o1           |
+	 *      .-------´
+	 * o2   |
+	 *
+	 *      <---d--->
+	 */
+	double d = (double)(left + right + 1);
+
+	double o1 = 0.5 + offset;
+	double o2 = 0.5 + offset - 1.0;
+
+	switch (pattern) {
+		case EDGESORTHO_NONE_NONE:
+		{
+			/*
+			 *
+			 *    ------
+			 *
+			 */
+			return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_POSI_NONE:
+		{
+			/*
+			 *
+			 *   .------
+			 *   |
+			 *
+			 * We only offset L patterns in the crossing edge side, to make it
+			 * converge with the unfiltered pattern 0 (we don't want to filter the
+			 * pattern 0 to avoid artifacts).
+			 */
+			if (left <= right)
+				return area(Dbl2(0.0, o2), Dbl2(d / 2.0, 0.0), left);
+			else
+				return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_NONE_POSI:
+		{
+			/*
+			 *
+			 *    ------.
+			 *          |
+			 */
+			if (left >= right)
+				return area(Dbl2(d / 2.0, 0.0), Dbl2(d, o2), left);
+			else
+				return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_POSI_POSI:
+		{
+			/*
+			 *
+			 *   .------.
+			 *   |      |
+			 */
+			if (m_orig_u) {
+				a1 = area(Dbl2(0.0, o2), Dbl2(d / 2.0, 0.0), left);
+				a2 = area(Dbl2(d / 2.0, 0.0), Dbl2(d, o2), left);
+				return smoothArea(d, a1, a2);
+			}
+			else
+				return area(makeQuad(left, d, o2), makeQuad(left + 1, d, o2), left);
+			break;
+		}
+		case EDGESORTHO_NEGA_NONE:
+		{
+			/*
+			 *   |
+			 *   `------
+			 *
+			 */
+			if (left <= right)
+				return area(Dbl2(0.0, o1), Dbl2(d / 2.0, 0.0), left);
+			else
+				return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_BOTH_NONE:
+		{
+			/*
+			 *   |
+			 *   +------
+			 *   |
+			 */
+			return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_NEGA_POSI:
+		{
+			/*
+			 *   |
+			 *   `------.
+			 *          |
+			 *
+			 * A problem of not offseting L patterns (see above), is that for certain
+			 * max search distances, the pixels in the center of a Z pattern will
+			 * detect the full Z pattern, while the pixels in the sides will detect a
+			 * L pattern. To avoid discontinuities, we blend the full offsetted Z
+			 * revectorization with partially offsetted L patterns.
+			 */
+			if (fabs(offset) > 0.0) {
+				a1 = area(Dbl2(0.0, o1), Dbl2(d, o2), left);
+				a2 = area(Dbl2(0.0, o1), Dbl2(d / 2.0, 0.0), left);
+				a2 += area(Dbl2(d / 2.0, 0.0), Dbl2(d, o2), left);
+				return (a1 + a2) / Dbl2(2.0);
+			}
+			else
+				return area(Dbl2(0.0, o1), Dbl2(d, o2), left);
+			break;
+		}
+		case EDGESORTHO_BOTH_POSI:
+		{
+			/*
+			 *   |
+			 *   +------.
+			 *   |      |
+			 */
+			return area(Dbl2(0.0, o1), Dbl2(d, o2), left);
+			break;
+		}
+		case EDGESORTHO_NONE_NEGA:
+		{
+			/*
+			 *          |
+			 *    ------´
+			 *
+			 */
+			if (left >= right)
+				return area(Dbl2(d / 2.0, 0.0), Dbl2(d, o1), left);
+			else
+				return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_POSI_NEGA:
+		{
+			/*
+			 *          |
+			 *   .------´
+			 *   |
+			 */
+			if (fabs(offset) > 0.0) {
+				a1 = area(Dbl2(0.0, o2), Dbl2(d, o1), left);
+				a2 = area(Dbl2(0.0, o2), Dbl2(d / 2.0, 0.0), left);
+				a2 += area(Dbl2(d / 2.0, 0.0), Dbl2(d, o1), left);
+				return (a1 + a2) / Dbl2(2.0);
+			}
+			else
+				return area(Dbl2(0.0, o2), Dbl2(d, o1), left);
+			break;
+		}
+		case EDGESORTHO_NONE_BOTH:
+		{
+			/*
+			 *          |
+			 *    ------+
+			 *          |
+			 */
+			return Dbl2(0.0, 0.0);
+			break;
+		}
+		case EDGESORTHO_POSI_BOTH:
+		{
+			/*
+			 *          |
+			 *   .------+
+			 *   |      |
+			 */
+			return area(Dbl2(0.0, o2), Dbl2(d, o1), left);
+			break;
+		}
+		case EDGESORTHO_NEGA_NEGA:
+		{
+			/*
+			 *   |      |
+			 *   `------´
+			 *
+			 */
+			if (m_orig_u) {
+				a1 = area(Dbl2(0.0, o1), Dbl2(d / 2.0, 0.0), left);
+				a2 = area(Dbl2(d / 2.0, 0.0), Dbl2(d, o1), left);
+				return smoothArea(d, a1, a2);
+			}
+			else
+				return area(makeQuad(left, d, o1), makeQuad(left + 1, d, o1), left);
+			break;
+		}
+		case EDGESORTHO_BOTH_NEGA:
+		{
+			/*
+			 *   |      |
+			 *   +------´
+			 *   |
+			 */
+			return area(Dbl2(0.0, o2), Dbl2(d, o1), left);
+			break;
+		}
+		case EDGESORTHO_NEGA_BOTH:
+		{
+			/*
+			 *   |      |
+			 *   `------+
+			 *          |
+			 */
+			return area(Dbl2(0.0, o1), Dbl2(d, o2), left);
+			break;
+		}
+		case EDGESORTHO_BOTH_BOTH:
+		{
+			/*
+			 *   |      |
+			 *   +------+
+			 *   |      |
+			 */
+			return Dbl2(0.0, 0.0);
+			break;
+		}
+	}
+
+	return Dbl2(0.0, 0.0);
+}
+
+/*------------------------------------------------------------------------------*/
+/* Diagonal Areas */
+
+class AreaDiag {
+	double m_data[SUBSAMPLES_DIAG][TEX_SIZE_DIAG][TEX_SIZE_DIAG][2];
+	bool m_numeric;
+	bool m_orig_u;
+public:
+	AreaDiag(bool numeric, bool orig_u) : m_numeric(numeric), m_orig_u(orig_u) {}
+
+	double *getData() { return (double *)&m_data; }
+	Dbl2 getPixel(int offset_index, Int2 coords) {
+		return Dbl2(m_data[offset_index][coords.y][coords.x][0],
+			    m_data[offset_index][coords.y][coords.x][1]);
+	}
+
+	void areaTex(int offset_index);
+private:
+	void putPixel(int offset_index, Int2 coords, Dbl2 pixel) {
+		m_data[offset_index][coords.y][coords.x][0] = pixel.x;
+		m_data[offset_index][coords.y][coords.x][1] = pixel.y;
+	}
+
+	double area1(Dbl2 p1, Dbl2 p2, Int2 p);
+	Dbl2 area(Dbl2 p1, Dbl2 p2, int left);
+	Dbl2 areaTriangle(Dbl2 p1L, Dbl2 p2L, Dbl2 p1R, Dbl2 p2R, int left);
+	Dbl2 calculate(int pattern, int left, int right, Dbl2 offset);
+};
+
+/* Calculates the area under the line p1->p2 for the pixel 'p' using brute */
+/* force sampling: */
+/* (quick and dirty solution, but it works) */
+double AreaDiag::area1(Dbl2 p1, Dbl2 p2, Int2 p)
+{
+	if (p1 == p2)
+		return 1.0;
+
+	double xm = (p1.x + p2.x) / 2.0, ym = (p1.y + p2.y) / 2.0;
+	double a = p2.y - p1.y;
+	double b = p1.x - p2.x;
+	int count = 0;
+
+	for (int ix = 0; ix < SAMPLES_DIAG; ix++) {
+		double x = (double)p.x + (double)ix / (double)(SAMPLES_DIAG - 1);
+		for (int iy = 0; iy < SAMPLES_DIAG; iy++) {
+			double y = (double)p.y + (double)iy / (double)(SAMPLES_DIAG - 1);
+			if (a * (x - xm) + b * (y - ym) > 0.0) /* inside? */
+				count++;
+		}
+	}
+	return (double)count / (double)(SAMPLES_DIAG * SAMPLES_DIAG);
+}
+
+/* Calculates the area under the line p1->p2: */
+/* (includes the pixel and its opposite) */
+Dbl2 AreaDiag::area(Dbl2 p1, Dbl2 p2, int left)
+{
+	if (m_numeric) {
+		double a1 = area1(p1, p2, Int2(1, 0) + Int2(left));
+		double a2 = area1(p1, p2, Int2(1, 1) + Int2(left));
+		return Dbl2(1.0 - a1, a2);
+	}
+
+	/* Calculates the area under the line p1->p2 for the pixel 'p' analytically */
+	Dbl2 d = p2 - p1;
+	if (d.x == 0.0)
+		return Dbl2(0.0, 1.0);
+
+	double x1 = (double)(1 + left);
+	double x2 = x1 + 1.0;
+	double ymid = x1;
+	double xtop = p1.x + (ymid + 1.0 - p1.y) * d.x / d.y;
+	double xmid = p1.x + (ymid       - p1.y) * d.x / d.y;
+	double xbot = p1.x + (ymid - 1.0 - p1.y) * d.x / d.y;
+
+	double y1 = p1.y + (x1 - p1.x) * d.y / d.x;
+	double y2 = p1.y + (x2 - p1.x) * d.y / d.x;
+	double fy1 = y1 - floor(y1);
+	double fy2 = y2 - floor(y2);
+	int iy1 = (int)floor(y1 - ymid);
+	int iy2 = (int)floor(y2 - ymid);
+
+	if (iy1 <= -2) {
+		if (iy2 == -1)
+			return Dbl2(1.0 - (x2 - xbot) * fy2 * 0.5, 0.0);
+		else if (iy2 == 0)
+			return Dbl2((xmid + xbot) * 0.5 - x1, (x2 - xmid) * fy2 * 0.5);
+		else if (iy2 >= 1)
+			return Dbl2((xmid + xbot) * 0.5 - x1, x2 -  (xtop + xmid) * 0.5);
+		else /* iy2 < -1 */
+			return Dbl2(1.0, 0.0);
+	}
+	else if (iy1 == -1) {
+		if (iy2 == -1)
+			return Dbl2(1.0 - (fy1 + fy2) * 0.5, 0.0);
+		else if (iy2 == 0)
+			return Dbl2((xmid - x1) * (1.0 - fy1) * 0.5, (x2 - xmid) * fy2 * 0.5);
+		else if (iy2 >= 1)
+			return Dbl2((xmid - x1) * (1.0 - fy1) * 0.5, x2 - (xtop + xmid) * 0.5);
+		else /* iy2 < -1 */
+			return Dbl2(1.0 - (xbot - x1) * fy1 * 0.5, 0.0);
+	}
+	else if (iy1 == 0) {
+		if (iy2 == -1)
+			return Dbl2((x2 - xmid) * (1.0 - fy2) * 0.5, (xmid - x1) * fy1 * 0.5);
+		else if (iy2 == 0)
+			return Dbl2(0.0, (fy1 + fy2) * 0.5);
+		else if (iy2 >= 1)
+			return Dbl2(0.0, 1.0 - (xtop - x1) * (1.0 - fy1) * 0.5);
+		else /* iy2 < -1 */
+			return Dbl2(x2 - (xmid + xbot) * 0.5, (xmid - x1) * fy1 * 0.5);
+	}
+	else { /* iy1 > 0 */
+		if (iy2 == -1)
+			return Dbl2((x2 - xtop) * (1.0 - fy2) * 0.5, (xtop + xmid) * 0.5 - x1);
+		else if (iy2 == 0)
+			return Dbl2(0.0, 1.0 - (x1 - xtop) * (1.0 - fy2) * 0.5);
+		else if (iy2 >= 1)
+			return Dbl2(0.0, 1.0);
+		else /* iy2 < -1 */
+			return Dbl2(x2 - (xmid + xbot) * 0.5, (xtop + xmid) * 0.5 - x1);
+	}
+}
+
+/* Calculate u-patterns using a triangle: */
+Dbl2 AreaDiag::areaTriangle(Dbl2 p1L, Dbl2 p2L, Dbl2 p1R, Dbl2 p2R, int left)
+{
+	double x1 = (double)(1 + left);
+	double x2 = x1 + 1.0;
+
+	Dbl2 dL = p2L - p1L;
+	Dbl2 dR = p2R - p1R;
+	double xm = ((p1L.x * dL.y / dL.x - p1L.y) - (p1R.x * dR.y / dR.x - p1R.y)) / (dL.y / dL.x - dR.y / dR.x);
+
+	double y1 = (x1 < xm) ? p1L.y + (x1 - p1L.x) * dL.y / dL.x : p1R.y + (x1 - p1R.x) * dR.y / dR.x;
+	double y2 = (x2 < xm) ? p1L.y + (x2 - p1L.x) * dL.y / dL.x : p1R.y + (x2 - p1R.x) * dR.y / dR.x;
+
+	return area(Dbl2(x1, y1), Dbl2(x2, y2), left);
+}
+
+/* Calculates the area for a given pattern and distances to the left and to the */
+/* right, biased by an offset: */
+Dbl2 AreaDiag::calculate(int pattern, int left, int right, Dbl2 offset)
+{
+	Dbl2 a1, a2;
+
+	double d = (double)(left + right + 1);
+
+	/*
+	 * There is some Black Magic around diagonal area calculations. Unlike
+	 * orthogonal patterns, the 'null' pattern (one without crossing edges) must be
+	 * filtered, and the ends of both the 'null' and L patterns are not known: L
+	 * and U patterns have different endings, and we don't know what is the
+	 * adjacent pattern. So, what we do is calculate a blend of both possibilites.
+	 */
+	switch (pattern) {
+		case EDGESDIAG_NONE_NONE:
+		{
+			/*
+			 *
+			 *         .-´
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   ´
+			 *
+			 */
+			a1 = area(Dbl2(1.0, 1.0), Dbl2(1.0, 1.0) + Dbl2(d), left); /* 1st possibility */
+			a2 = area(Dbl2(1.0, 0.0), Dbl2(1.0, 0.0) + Dbl2(d), left); /* 2nd possibility */
+			return (a1 + a2) / Dbl2(2.0); /* Blend them */
+			break;
+		}
+		case EDGESDIAG_VERT_NONE:
+		{
+			/*
+			 *
+			 *         .-´
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   |
+			 *   |
+			 */
+			a1 = area(Dbl2(1.0, 0.0) + offset, Dbl2(0.0, 0.0) + Dbl2(d), left);
+			a2 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d), left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_NONE_HORZ:
+		{
+			/*
+			 *
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   ´
+			 *
+			 */
+			a1 = area(Dbl2(0.0, 0.0), Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0), Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_VERT_HORZ:
+		{
+			/*
+			 *
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   |
+			 *   |
+			 */
+			if (m_orig_u)
+				return area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			else
+				return areaTriangle(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d),
+						    Dbl2(0.0, 0.0), Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			break;
+		}
+		case EDGESDIAG_HORZ_NONE:
+		{
+			/*
+			 *
+			 *         .-´
+			 *       .-´
+			 *     .-´
+			 * ----´
+			 *
+			 *
+			 */
+			a1 = area(Dbl2(1.0, 1.0) + offset, Dbl2(0.0, 0.0) + Dbl2(d), left);
+			a2 = area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d), left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_BOTH_NONE:
+		{
+			/*
+			 *
+			 *         .-´
+			 *       .-´
+			 *     .-´
+			 * --.-´
+			 *   |
+			 *   |
+			 */
+			a1 = area(Dbl2(1.0, 1.0) + offset, Dbl2(0.0, 0.0) + Dbl2(d), left);
+			a2 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d), left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_HORZ_HORZ:
+		{
+			/*
+			 *
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 * ----´
+			 *
+			 *
+			 */
+			return area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			break;
+		}
+		case EDGESDIAG_BOTH_HORZ:
+		{
+			/*
+			 *
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 * --.-´
+			 *   |
+			 *   |
+			 */
+			a1 = area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_NONE_VERT:
+		{
+			/*
+			 *         |
+			 *         |
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   ´
+			 *
+			 */
+			a1 = area(Dbl2(0.0, 0.0), Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0), Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_VERT_VERT:
+		{
+			/*
+			 *         |
+			 *         |
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   |
+			 *   |
+			 */
+			return area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			break;
+		}
+		case EDGESDIAG_NONE_BOTH:
+		{
+			/*
+			 *         |
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   ´
+			 *
+			 */
+			a1 = area(Dbl2(0.0, 0.0), Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0), Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_VERT_BOTH:
+		{
+			/*
+			 *         |
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 *   .-´
+			 *   |
+			 *   |
+			 */
+			a1 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_HORZ_VERT:
+		{
+			/*
+			 *         |
+			 *         |
+			 *       .-´
+			 *     .-´
+			 * ----´
+			 *
+			 *
+			 */
+			if (m_orig_u)
+				return area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			else
+				return areaTriangle(Dbl2(1.0, 1.0) + offset, Dbl2(2.0, 1.0) + Dbl2(d),
+						    Dbl2(1.0, 0.0), Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			break;
+		}
+		case EDGESDIAG_BOTH_VERT:
+		{
+			/*
+			 *         |
+			 *         |
+			 *       .-´
+			 *     .-´
+			 * --.-´
+			 *   |
+			 *   |
+			 */
+			a1 = area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_HORZ_BOTH:
+		{
+			/*
+			 *         |
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 * ----´
+			 *
+			 *
+			 */
+			a1 = area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+		case EDGESDIAG_BOTH_BOTH:
+		{
+			/*
+			 *         |
+			 *         .----
+			 *       .-´
+			 *     .-´
+			 * --.-´
+			 *   |
+			 *   |
+			 */
+			a1 = area(Dbl2(1.0, 1.0) + offset, Dbl2(1.0, 1.0) + Dbl2(d) + offset, left);
+			a2 = area(Dbl2(1.0, 0.0) + offset, Dbl2(1.0, 0.0) + Dbl2(d) + offset, left);
+			return (a1 + a2) / Dbl2(2.0);
+			break;
+		}
+	}
+
+	return Dbl2(0.0, 0.0);
+}
+
+/*------------------------------------------------------------------------------*/
+/* Main Loops */
+
+void AreaOrtho::areaTex(int offset_index)
+{
+	double offset = subsample_offsets_ortho[offset_index];
+	int max_dist = m_compat ? MAX_DIST_ORTHO_COMPAT : MAX_DIST_ORTHO;
+
+	for (int pattern = 0; pattern < 16; pattern++) {
+		Int2 e = Int2(max_dist) * (m_compat ? edgesortho_compat : edgesortho)[pattern];
+		for (int left = 0; left < max_dist; left++) {
+			for (int right = 0; right < max_dist; right++) {
+				Dbl2 p = calculate(pattern, left * left, right * right, offset);
+				Int2 coords = e + Int2(left, right);
+
+				putPixel(offset_index, coords, p);
+			}
+		}
+	}
+	return;
+}
+
+void AreaDiag::areaTex(int offset_index)
+{
+	Dbl2 offset = subsample_offsets_diag[offset_index];
+
+	for (int pattern = 0; pattern < 16; pattern++) {
+		Int2 e = Int2(MAX_DIST_DIAG) * edgesdiag[pattern];
+		for (int left = 0; left < MAX_DIST_DIAG; left++) {
+			for (int right = 0; right < MAX_DIST_DIAG; right++) {
+				Dbl2 p = calculate(pattern, left, right, offset);
+				Int2 coords = e + Int2(left, right);
+
+				putPixel(offset_index, coords, p);
+			}
+		}
+	}
+	return;
+}
+
+/*------------------------------------------------------------------------------*/
+/* Write File to Specified Location on Disk */
+
+/* C/C++ source code (arrays of floats) */
+static void write_double_array(FILE *fp, const double *ptr, int length, const char *array_name, bool quantize)
+{
+	fprintf(fp, "static const float %s[%d] = {", array_name, length);
+
+	for (int n = 0; n < length; n++) {
+		if (n > 0)
+			fprintf(fp, ",");
+		fprintf(fp, (n % 8 != 0) ? " " : "\n\t");
+
+		if (quantize)
+			fprintf(fp, "%3d / 255.0", (int)(*(ptr++) * 255.0));
+		else
+			fprintf(fp, "%1.8lf", *(ptr++));
+	}
+
+	fprintf(fp, "\n};\n");
+}
+
+static void write_csource(AreaOrtho *ortho, AreaDiag *diag, FILE *fp, bool subsampling, bool quantize)
+{
+	fprintf(fp, "/* This file was generated by smaa_areatex.cpp */\n");
+
+	fprintf(fp, "\n/* Horizontal/Vertical Areas */\n");
+	write_double_array(fp, ortho->getData(),
+			   TEX_SIZE_ORTHO * TEX_SIZE_ORTHO * 2 * (subsampling ? SUBSAMPLES_ORTHO : 1),
+			   "areatex", quantize);
+
+	fprintf(fp, "\n/* Diagonal Areas */\n");
+	write_double_array(fp, diag->getData(),
+			   TEX_SIZE_DIAG * TEX_SIZE_DIAG * 2 * (subsampling ? SUBSAMPLES_DIAG : 1),
+			   "areatex_diag", quantize);
+}
+
+/* .tga File (RGBA 32bit uncompressed) */
+static void write_tga(AreaOrtho *ortho, AreaDiag *diag, FILE *fp, bool subsampling)
+{
+	int subsamples = subsampling ? SUBSAMPLES_ORTHO : 1;
+	unsigned char header[18] = {0, 0,
+				    2,   /* uncompressed RGB */
+				    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+				    32,  /* 32bit */
+				    8};  /* 8bit alpha, left to right, bottom to top */
+
+	/* Set width and height */
+	header[12] = (TEX_SIZE_ORTHO + TEX_SIZE_DIAG)      & 0xff;
+	header[13] = ((TEX_SIZE_ORTHO + TEX_SIZE_DIAG) >> 8) & 0xff;
+	header[14] = (subsamples * TEX_SIZE_ORTHO)      & 0xff;
+	header[15] = ((subsamples * TEX_SIZE_ORTHO) >> 8) & 0xff;
+
+	/* Write .tga header */
+	fwrite(header, sizeof(unsigned char), sizeof(header) / sizeof(unsigned char), fp);
+
+	/* Write pixel data  */
+	for (int i = subsamples - 1; i >= 0; i--) {
+		for (int y = TEX_SIZE_ORTHO - 1; y >= 0; y--) {
+			for (int x = 0; x < TEX_SIZE_ORTHO; x++) {
+				Dbl2 p = ortho->getPixel(i, Int2(x, y));
+				fputc(0, fp);                            /* B */
+				fputc((unsigned char)(p.y * 255.0), fp); /* G */
+				fputc((unsigned char)(p.x * 255.0), fp); /* R */
+				fputc(0, fp);                            /* A */
+			}
+
+			for (int x = 0; x < TEX_SIZE_DIAG; x++) {
+				if (i < SUBSAMPLES_DIAG) {
+					Dbl2 p = diag->getPixel(i, Int2(x, y));
+					fputc(0, fp);                            /* B */
+					fputc((unsigned char)(p.y * 255.0), fp); /* G */
+					fputc((unsigned char)(p.x * 255.0), fp); /* R */
+					fputc(0, fp);                            /* A */
+				}
+				else {
+					fputc(0, fp);
+					fputc(0, fp);
+					fputc(0, fp);
+					fputc(0, fp);
+				}
+			}
+		}
+	}
+}
+
+/* .raw File (R8G8 raw data) */
+static void write_raw(AreaOrtho *ortho, AreaDiag *diag, FILE *fp, bool subsampling)
+{
+	int subsamples = subsampling ? SUBSAMPLES_ORTHO : 1;
+
+	/* Write pixel data  */
+	for (int i = 0; i < subsamples; i++) {
+		for (int y = 0; y < TEX_SIZE_ORTHO; y++) {
+			for (int x = 0; x < TEX_SIZE_ORTHO; x++) {
+				Dbl2 p = ortho->getPixel(i, Int2(x, y));
+				fputc((unsigned char)(p.x * 255.0), fp); /* R */
+				fputc((unsigned char)(p.y * 255.0), fp); /* G */
+			}
+
+			for (int x = 0; x < TEX_SIZE_DIAG; x++) {
+				if (i < SUBSAMPLES_DIAG) {
+					Dbl2 p = diag->getPixel(i, Int2(x, y));
+					fputc((unsigned char)(p.x * 255.0), fp); /* R */
+					fputc((unsigned char)(p.y * 255.0), fp); /* G */
+				}
+				else {
+					fputc(0, fp);
+					fputc(0, fp);
+				}
+			}
+		}
+	}
+}
+
+static int generate_file(AreaOrtho *ortho, AreaDiag *diag, const char *path, bool subsampling, bool quantize, bool tga, bool raw)
+{
+	FILE *fp = fopen(path, tga ? "wb" : "w");
+
+	if (!fp) {
+		fprintf(stderr, "Unable to open file: %s\n", path);
+		return 1;
+	}
+
+	fprintf(stderr, "Generating %s\n", path);
+
+	if (tga)
+		write_tga(ortho, diag, fp, subsampling);
+	else if (raw)
+		write_raw(ortho, diag, fp, subsampling);
+	else
+		write_csource(ortho, diag, fp, subsampling, quantize);
+
+	fclose(fp);
+
+	return 0;
+}
+
+int main(int argc, char **argv)
+{
+	bool subsampling = false;
+	bool quantize = false;
+	bool tga = false;
+	bool raw = false;
+	bool compat = false;
+	bool numeric = false;
+	bool orig_u = false;
+	bool help = false;
+	char *outfile = NULL;
+	int status = 0;
+
+	for (int i = 1; i < argc; i++) {
+		char *ptr = argv[i];
+		if (*ptr++ == '-' && *ptr != '\0') {
+			char c;
+			while ((c = *ptr++) != '\0') {
+				if (c == 's')
+					subsampling = true;
+				else if (c == 'q')
+					quantize = true;
+				else if (c == 't')
+					tga = true;
+				else if (c == 'r')
+					raw = true;
+				else if (c == 'c')
+					compat = true;
+				else if (c == 'n')
+					numeric = true;
+				else if (c == 'u')
+					orig_u = true;
+				else if (c == 'h')
+					help = true;
+				else {
+					fprintf(stderr, "Unknown option: -%c\n", c);
+					status = 1;
+					break;
+				}
+			}
+		}
+		else if (outfile) {
+			fprintf(stderr, "Too much file names: %s, %s\n", outfile, argv[i]);
+			status = 1;
+		}
+		else
+			outfile = argv[i];
+
+		if (status != 0)
+			break;
+	}
+
+	if (status == 0 && !help && !outfile) {
+		fprintf(stderr, "File name was not specified.\n");
+		status = 1;
+	}
+
+	if (status != 0 || help) {
+		fprintf(stderr, "Usage: %s [OPTION]... OUTFILE\n", argv[0]);
+		fprintf(stderr, "Options:\n");
+		fprintf(stderr, "    -s    Calculate data for subpixel rendering\n");
+		fprintf(stderr, "    -q    Quantize data to 256 levels\n");
+		fprintf(stderr, "    -t    Write TGA image instead of C/C++ source\n");
+		fprintf(stderr, "    -r    Write R8G8 raw image instead of C/C++ source\n");
+		fprintf(stderr, "    -c    Generate compatible orthogonal data that subtexture size is 16\n");
+		fprintf(stderr, "    -n    Numerically calculate diagonal data using brute force sampling\n");
+		fprintf(stderr, "    -u    Process orthogonal / diagonal U patterns in older ways\n");
+		fprintf(stderr, "    -h    Print this help and exit\n");
+		fprintf(stderr, "File name OUTFILE usually should have an extension such as .c, .h, or .tga,\n");
+		fprintf(stderr, "except for a special name '-' that means standard output.\n\n");
+		fprintf(stderr, "Example:\n");
+		fprintf(stderr, "  Generate TGA file exactly same as AreaTexDX10.tga bundled with the\n");
+		fprintf(stderr, "  original implementation:\n\n");
+		fprintf(stderr, "  $ smaa_areatex -stcnu AreaTexDX10.tga\n\n");
+		return status;
+	}
+
+	AreaOrtho *ortho = new AreaOrtho(compat, orig_u);
+	AreaDiag *diag = new AreaDiag(numeric, orig_u);
+
+	/* Calculate areatex data */
+	for (int i = 0; i < (subsampling ? SUBSAMPLES_ORTHO : 1); i++)
+		ortho->areaTex(i);
+
+	for (int i = 0; i < (subsampling ? SUBSAMPLES_DIAG : 1); i++)
+		diag->areaTex(i);
+
+	/* Generate .tga, .raw, or C/C++ source file, or write the data to stdout */
+	if (strcmp(outfile, "-") != 0)
+		status = generate_file(ortho, diag, outfile, subsampling, quantize, tga, raw);
+	else if (tga)
+		write_tga(ortho, diag, stdout, subsampling);
+	else if (raw)
+		write_raw(ortho, diag, stdout, subsampling);
+	else
+		write_csource(ortho, diag, stdout, subsampling, quantize);
+
+	delete ortho;
+	delete diag;
+
+	return status;
+}
+
+/* smaa_areatex.cpp ends here */