diff --git a/.idea/.name b/.idea/.name
index a4d4a59..67a735f 100644
--- a/.idea/.name
+++ b/.idea/.name
@@ -1 +1 @@
-vk_raytracing_tutorial
\ No newline at end of file
+targetrdf
\ No newline at end of file
diff --git a/.idea/misc.xml b/.idea/misc.xml
index 0b76fe5..07b89be 100644
--- a/.idea/misc.xml
+++ b/.idea/misc.xml
@@ -3,5 +3,7 @@
   <component name="CMakePythonSetting">
     <option name="pythonIntegrationState" value="YES" />
   </component>
-  <component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$" />
+  <component name="CMakeWorkspace" PROJECT_DIR="$PROJECT_DIR$/targetrdf">
+    <contentRoot DIR="$PROJECT_DIR$" />
+  </component>
 </project>
\ No newline at end of file
diff --git a/targetrdf/CMakeLists.txt b/targetrdf/CMakeLists.txt
new file mode 100644
index 0000000..3925604
--- /dev/null
+++ b/targetrdf/CMakeLists.txt
@@ -0,0 +1,22 @@
+cmake_minimum_required(VERSION 3.20)
+project(targetrdf CXX)
+
+if (CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
+    set(CMAKE_CXX_STANDARD 20)
+endif ()
+
+set(CMAKE_CXX_FLAGS "-O3 -Wall -g")
+
+set(SOURCE_FILES
+        common.cc
+        param.cc
+        targetrdf.cc
+        targetrdf.cc
+)
+
+set(HEADER_FILES
+        common.hh
+        param.hh
+)
+
+add_executable(${PROJECT_NAME} ${SOURCE_FILES} ${HEADER_FILES})
\ No newline at end of file
diff --git a/targetrdf/README.md b/targetrdf/README.md
new file mode 100644
index 0000000..a687423
--- /dev/null
+++ b/targetrdf/README.md
@@ -0,0 +1,108 @@
+This directory contains the source code for the paper
+
+  D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
+  Controlled Aliasing", ACM Trans. Graph., 2013, in press
+
+If you have questions about the paper or the included source code, please 
+email
+
+	Daniel Heck <dheck@gmx.de>
+
+The main program is called "targetrdf", and it attempts to construct a point
+set matching a given "target RDF" (or a given target spectrum). As a simple
+example, here is how to generate a step blue noise point set, as presented in
+the paper:
+
+	./targetrdf --steprp --crit 0.606 --out step.txt
+
+
+  --steprp      
+
+        Chooses a radial power spectrum with a step characteristic
+
+  --crit 0.606  
+
+        Sets the position of the step; 0.606 is the numerical value for v_max,
+        the maximal realizable step position
+
+  --out step.txt
+
+        Writes the final point set to 'step.txt'
+
+
+There are many additional parameters that can be used to influence the
+optimization process. A full list is shown by "targetrdf --help".
+
+
+  --npts n
+  --seed s
+
+        Specifies the number of points to use. The points are initialized to
+        random positions, and the seed for the RNG can be set to different
+        values to generate a family of similar point sets.
+
+  --in file
+
+        Load initial point set from 'file'
+
+  --reference file
+
+        Loads a point set from 'file' and tries to construct a new point set
+        with the same RDF.
+
+  --rdf file
+  --spectrum file
+  --rp-as-rdf file
+
+        Loads an RDF or power spectrum from 'file' and tries to construct a
+        matching point set. 'rp-as-rdf' can be used to generate "dual" point
+        sets by reinterpreting a power spectrum as an RDF.
+
+  --peakrp
+  --peak power
+  --peaksmooth sigma
+  --crit
+
+        These paramters are used to generate the single-peak blue noise point
+        sets discussed in the paper.
+
+  --nbins n
+  --smoothing sigma
+
+        Specifies the number of bins and the amount of smoothing to use when
+        calculating RDFs.
+
+FILE FORMATS
+------------
+
+The default file formats used by 'targetrdf' are very simple.
+
+Point sets:
+
+  - the first line contains the number of points
+  - the remaining lines contain the x/y coordinates of the points 
+    in the unit square
+
+        4096
+        0.162776 0.74937
+        0.0798951 0.880201
+        0.58569 0.782718
+        0.697048 0.380065
+        0.883548 0.747492
+        0.432829 0.344747
+        ...
+
+RDF/Power spectrum:
+
+  - the file contains ordinate/abscissa pairs:
+
+        6.10352e-05 0.00088529
+        0.000183105 0.00107669
+        0.000305176 0.00145942
+        ...
+        0.499695 1.0007
+        0.499817 1.00067
+        0.499939 1.00064
+
+  - it is assumed that ordinate spacing is uniform
+
diff --git a/targetrdf/common.cc b/targetrdf/common.cc
new file mode 100644
index 0000000..baffddf
--- /dev/null
+++ b/targetrdf/common.cc
@@ -0,0 +1,399 @@
+//
+// Source code for the paper
+// 
+// D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
+// Controlled Aliasing", ACM Trans. Graph., 2013, in press
+//
+// Copyright (C) 2012,2013 Daniel Heck and Thomas Schloemer
+//
+#include "common.hh"
+#include <fstream>
+#include <iostream>
+#include <string.h>
+#include <cstdlib>
+
+// Precomputed arrays for Bessel function computation
+#define PI_4    0.7853981633974483096f
+#define DR1     5.7831859629467845211f
+
+static const float JP[5] = {
+    -6.068350350393235E-008f,
+    6.388945720783375E-006f,
+    -3.969646342510940E-004f,
+    1.332913422519003E-002f,
+    -1.729150680240724E-001f
+};
+static const float MO[8] = {
+    -6.838999669318810E-002f,
+    1.864949361379502E-001f,
+    -2.145007480346739E-001f,
+    1.197549369473540E-001f,
+    -3.560281861530129E-003f,
+    -4.969382655296620E-002f,
+    -3.355424622293709E-006f,
+    7.978845717621440E-001f
+};
+static const float PH[8] = {
+    3.242077816988247E+001f,
+    -3.630592630518434E+001f,
+    1.756221482109099E+001f,
+    -4.974978466280903E+000f,
+    1.001973420681837E+000f,
+    -1.939906941791308E-001f,
+    6.490598792654666E-002f,
+    -1.249992184872738E-001f
+};
+
+// Single precision approximation to j0, about twice as fast as the standard
+// implementation. From cephes lib, available at http://www.netlib.org/cephes/
+inline float Polynomial(float x, const float *p, int n) {
+    float y = p[0];
+    for (int i = 1; i < n; ++i)
+        y = y * x + p[i];
+    return y;
+}
+
+inline float j0f(float f) {
+    float x = fabsf(f);
+    if (x <= 2.f) {
+        float z = x*x;
+        if (x < 1.0e-3f)
+            return (1.f - 0.25f*z);
+        return (z - DR1) * Polynomial(z, JP, 5);
+	}
+    float q = 1.f / x;
+    float w = sqrtf(q);
+    float p = w * Polynomial(q, MO, 8);
+    w = q*q;
+    float xn = q * Polynomial(w, PH, 8) - PI_4;
+    return p * cosf(xn + x);
+}
+
+bool HasSuffix(const std::string &s, const std::string &suffix) {
+    int i = s.size() - suffix.size();
+    if (i >= 0)
+        return s.substr(i) == suffix;
+    return false;
+}
+
+Endianness SystemEndianness() {
+    union { uint32_t i; char c[4]; } e;
+    e.i = 1;
+    return e.c[0] == 0 ? BigEndian : LittleEndian;
+}
+
+void SwapEndian4(uint8_t *buf, int nn) {
+    for (int i=0; i<nn; i++) {
+        std::swap(buf[4*i], buf[4*i+3]);
+        std::swap(buf[4*i+1], buf[4*i+2]);
+    }
+}
+
+bool WriteFloats(FILE *fp, const float *p, int n, Endianness endian) {
+    Endianness sysEndian = SystemEndianness();
+    const int nbuf = 256;
+    uint8_t buf[nbuf*sizeof(float)];
+    while (n > 0) {
+        int nn = std::min(n, nbuf);
+        memcpy(buf, p, nn * sizeof(float));
+        if (endian != sysEndian)
+            SwapEndian4(buf, nn);
+        if (fwrite(buf, sizeof(float), nn, fp) != (size_t)nn)
+            return false;
+        p += nn; n -= nn;
+    }
+    return true;
+}
+
+bool ReadFloats(FILE *fp, float *p, int n, Endianness endian) {
+    Endianness sysEndian = SystemEndianness();
+    const int nbuf = 256;
+    uint8_t buf[nbuf*sizeof(float)];
+    while (n > 0) {
+        int nn = std::min(n, nbuf);
+        if (fread(buf, sizeof(float), nn, fp) != (size_t)nn)
+            return false;
+        if (endian != sysEndian)
+            SwapEndian4(buf, nn);
+        memcpy(p, buf, nn*sizeof(float));
+        p += nn; n -= nn;
+    }
+    return true;
+}
+
+void ReadPFM(int &width, int &height, float *&data, const char *fname) {
+    FILE* fp = fopen(fname, "rb");
+    char name[100];
+    float f;
+    fscanf(fp, "%s\n%d %d %f\n", name, &width, &height, &f);
+    printf("f is %f\n", f);
+    data = new float[width * height];
+    fread((char*) data, sizeof(float), width * height, fp);
+    fclose(fp);
+}
+
+void WritePFM(int width, int height, const float *data, const char *fname) {
+    FILE* fp = fopen(fname, "wb");
+    fprintf(fp, "Pf\n");
+    fprintf(fp, "%d %d\n-1.00\n", width, height);
+    fwrite((const char*) data, sizeof(float), width * height, fp);
+    fclose(fp);
+}
+
+
+void LoadPoints(const std::string &fname, std::vector<Point> &points) {
+    int npts = 0;
+    if (HasSuffix(fname, ".rps")) {
+        FILE *fp= fopen(fname.c_str(), "rb");
+        if (!fp) {
+            std::cerr << "Cannot load '" << fname << "'.\n";
+            exit(1);
+        }
+        fseek(fp, 0, SEEK_END);
+        npts = ftell(fp) / 8;
+        points.resize(npts);
+        fseek(fp, 0, SEEK_SET);
+        ReadFloats(fp, &points[0].x, npts*2, LittleEndian);
+        fclose(fp);
+    } else {
+        std::ifstream fp(fname.c_str());
+        if (!fp) {
+            std::cerr << "Cannot load '" << fname << "'.\n";
+            exit(1);
+        }
+        fp >> npts;
+        points.reserve(npts);
+        Point pt;
+        while ((int)points.size() < npts && (fp >> pt.x >> pt.y))
+            points.push_back(pt);
+    }
+}
+
+void WritePoints(const std::string &fname, std::vector<Point> &points) {
+    std::ofstream fp(fname.c_str());
+    if (!fp) {
+        std::cerr << "Cannot create '" << fname << "'.\n";
+        exit(1);
+    }
+    fp << points.size() << "\n";
+    for (unsigned i = 0; i < points.size(); i++)
+        fp << points[i].x << " " << points[i].y << "\n"; 
+}
+
+
+void SavePGM(float *p, int rows, int cols, const char *fname) {
+    FILE *fp = fopen(fname, "wb");
+    if (!fp) {
+        fprintf(stderr, "Cannot open '%s' for writing.\n", fname);
+        exit(1);
+    }
+
+    int maxval = 255;
+    fprintf(fp, "P2\n%d %d\n%d\n", cols, rows, maxval);
+
+    for (int y=0; y<rows; y++) {
+        for (int x=0; x<cols; x++) {
+            int pp = (int)(p[y*cols + x] * maxval);
+            if (pp < 0) pp = 0;
+            if (pp > maxval) pp = maxval;
+            fprintf(fp, "% 3d ", pp);
+        }
+        fprintf(fp, "\n");
+    }
+
+    fclose(fp);
+}
+
+
+// -------------------- Curve --------------------
+
+Curve::Curve(const Curve &c) {
+    x0 = c.x0; 
+    x1 = c.x1;
+    dx = c.dx;
+    y = c.y;
+}
+
+Curve::Curve(int size, float x0, float x1) {
+    y.resize(size);
+    
+    this->x0 = x0;
+    this->x1 = x1;
+    dx = 1.f/size * (x1 - x0);
+}
+
+float Curve::At(float x) const {
+    float xx = (x-x0)/dx;
+    if (xx < 0) xx = 0;
+    int i = (int)floor(xx);
+    int ii = i+1;
+    float a = xx-i;
+    if (i >= (int)y.size()) i = y.size()-1;
+    if (ii >= (int)y.size()) ii = y.size()-1;
+    return (1-a)*y[i] + a * y[ii];
+}
+
+void Curve::Write(const std::string &fname) {
+    std::ofstream os(fname.c_str());
+    for (int j=0; j<size(); j++)
+        os << ToX(j) << " " << y[j] << "\n";
+}
+
+Curve Curve::Read(const std::string &fname) {
+    std::ifstream is(fname.c_str());
+    std::vector<float> x, y;
+    float xx, yy;
+    while (is >> xx >> yy) {
+        x.push_back(xx);
+        y.push_back(yy);
+    }
+    if (x.empty()) {
+        std::cerr << "Could not load RDF from '" << fname << "'.\n";
+        exit(1);
+    }
+
+    int n = x.size();
+    float dx = n > 1 ? x[1] - x[0] : 1;
+    Curve c(n, x[0]-dx/2, x[n-1] + dx/2);
+    for (int i=0; i<n; i++)
+        c.y[i] = y[i];
+    return c;
+}
+
+
+
+void Curve::FilterGauss(const Curve &source, float sigma) {
+    if (sigma == 0){
+        y = source.y;
+        return;
+    }
+  
+    float dd = (float)source.size() / size();
+  
+    for (int i=0; i<size(); i++) {
+        float a = 0.0f, sumw = 0.0f;
+
+        int jmin = std::max(0, (int)floor(dd * (i - 5 * sigma)));
+        int jmax = std::min((int)source.size()-1, (int)ceil(dd * (i + 5*sigma)));
+        for (int j=jmin; j<=jmax; j++) {
+            float d = j/dd-i;
+            float w = exp(-d*d/(2*sigma*sigma));
+            a += source[j] * w;
+            sumw += w;
+        }
+        y[i] = a / sumw;
+    }
+}
+
+float Integrate(const Curve &c, float x0, float x1, float *E) {
+    bool negate = false;
+    if (x0 > x1) {
+        std::swap(x0, x1);
+        negate = true;
+    }
+
+    float T = 0.0f, M = 0.0f;
+    int i0 = c.ToIndex(x0);
+    int i1 = c.ToIndex(x1+c.dx);
+    if (i0 < 0) i0 = 0;
+    if (i1 < 0) i1 = 0;
+    if (i0 >= c.size()) i0 = c.size()-1;
+    if (i1 >= c.size()) i1 = c.size()-1;
+
+    T = 0.5 * (c[i0] + c[i1]);
+    for (int i=i0+2; i<i1; i+=2) {
+        T += c[i];
+        M += c[i-1];
+    }
+    if (E != NULL)
+        *E = c.dx * 2.0f * fabs(T-M);
+    float a = c.dx*(T+M);
+    return negate ? -a : a;
+}
+
+
+// -------------------- RDF Calculation --------------------
+
+void CalcRDF(Curve &c, size_t npts, const float *pts, float (*distfunc)(const float *, const float *)) {
+    std::vector<unsigned long> bins(c.size());
+    for (unsigned j=0; j<npts; j++) {
+        for (unsigned k=j+1; k<npts; k++) {
+            float d = distfunc(&pts[2*j], &pts[2*k]);
+            int idx = c.ToIndex(d);
+            if (0 <= idx && idx < c.size())
+                bins[idx]++;
+//            if (d < c.x1)
+//                bins[c.ToIndex(d)]++;
+        }
+    }
+    const float scale = npts * (npts-1)/2 * M_PI * c.dx * c.dx;
+    for (int j=0; j<c.size(); j++)
+        c[j] = bins[j] / (scale * (2*j+1));
+}
+
+Curve CalcRDF(int numbins, size_t npts, const float *pts, 
+        float sigma, float (*distfunc)(const float *, const float *))
+{
+    Curve c(numbins, 0, 0.5f);
+    CalcRDF(c, npts, pts, distfunc);
+    if (sigma == 0)
+        return c;
+    else
+        return FilterGauss(numbins, c, sigma);
+}
+
+float HannWindow(float x, float xlim) {
+    return x > xlim ? 0.0f : Sqr(cosf(M_PI * x / (2*xlim)));
+}
+
+float BlackmanWindow(float x, float xlim) {
+    return x > xlim ? 0.0f : 0.43 + 0.5*cosf(M_PI*x/xlim) + 0.08*cosf(2*M_PI*x/xlim);
+}
+
+float HammingWindow(float x, float xlim) {
+    return x > xlim ? 0.0f : 0.54f + 0.46f * cosf(M_PI * x/xlim);
+}
+
+float WelchWindow(float x, float xlim) {
+    return x > xlim ? 0.0f : 1 - x*x/(xlim*xlim);
+}
+
+Curve RDF2Power(int npts, const Curve &rdf, int nbins, float x0, float x1) {
+    Curve power(nbins, x0, x1);
+    Curve tmp(rdf);
+    for (int i=0; i<nbins; i++) {
+        const float u0 = power.ToX(i);
+        const float u = 2*M_PI*u0;
+        const float wndsize = rdf.x1 * std::min(0.5f, std::max(0.2f, 4*u0/sqrtf(npts)));
+        for (int j=0; j<tmp.size(); j++) {
+            float x = rdf.ToX(j);
+            float wnd = BlackmanWindow(x, wndsize);
+//            float wnd = HannWindow(x, wndsize);
+            tmp[j] = (rdf[j]-1) * j0f(u*x) * x * wnd;
+        }
+        power[i] = fabs(1 + 2 * M_PI * Integrate(tmp) * npts);
+    }
+    return power;
+}
+
+Curve Power2RDF(int npts, const Curve &power, int nbins, float x0, float x1, float smoothing) {
+    Curve rdf(nbins, x0, x1);
+    Curve tmp(power);
+    for (int i=0; i<rdf.size(); i++) {
+        float u = rdf.ToX(i);
+        for (int j=0; j<tmp.size(); j++) {
+            float x = tmp.ToX(j);
+            float wnd = BlackmanWindow(x, power.x1);
+            tmp[j] = (power[j]-1) * jn(0, 2*M_PI * u * x) * x * wnd;
+        }
+        rdf[i] = 1 + 2 * M_PI * Integrate(tmp) / npts;
+    }
+    for (int i=0; i<rdf.size(); i++) {
+        if (rdf[i] < 0) {
+            printf("Warning: RDF has negative components\n");
+            break;
+        }
+    }
+    return FilterGauss(nbins, rdf, smoothing);
+}
+
diff --git a/targetrdf/common.hh b/targetrdf/common.hh
new file mode 100644
index 0000000..fd928d6
--- /dev/null
+++ b/targetrdf/common.hh
@@ -0,0 +1,138 @@
+//
+// Source code for the paper
+// 
+// D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
+// Controlled Aliasing", ACM Trans. Graph., 2013, in press
+//
+// Copyright (C) 2012,2013 Daniel Heck and Thomas Schloemer
+//
+#ifndef COMMON_HH_INCLUDED
+#define COMMON_HH_INCLUDED
+
+#include <string>
+#include <vector>
+#include <cmath>
+#include <stdint.h>
+
+// -------------------- Math --------------------
+
+template <class T, class U>
+inline T Clamp(T x, U min, U max) {
+    if (x < min)
+	return min;
+    else if (x > max)
+	return max;
+    else
+	return x;
+}
+
+
+// -------------------- Utility --------------------
+
+bool HasSuffix(const std::string &s, const std::string &suffix);
+
+enum Endianness {
+    BigEndian, LittleEndian
+};
+
+Endianness SystemEndianness();
+void SwapEndian4(uint8_t *buf, int n);
+bool WriteFloats(FILE *fp, const float *p, int n, Endianness endian);
+bool ReadFloats(FILE *fp, float *p, int n, Endianness endian);
+
+void ReadPFM(int &width, int &height, float *&data, const char *fname);
+void WritePFM(int width, int height, const float *data, const char *fname);
+
+// Output floating point matrix as PGM image
+void SavePGM(float *p, int rows, int cols, const char *fname);
+
+
+// -------------------- Point --------------------
+
+struct Point {
+    float x, y;
+    Point(float x = 0, float y = 0) : x(x), y(y) {}
+};
+
+void LoadPoints(const std::string &fname, std::vector<Point> &points);
+void WritePoints(const std::string &fname, std::vector<Point> &points);
+
+// -------------------- Curve --------------------
+
+class Curve {
+    std::vector<float> y;
+public:
+    float x0, x1, dx;
+    Curve() { x0 = x1 = dx = 0.0f; }
+    Curve(const Curve &c);
+    Curve(int size, float x0, float x1);
+
+    float &operator[](size_t i) { return y[i]; }
+    const float &operator[](size_t i) const { return y[i]; }
+
+    int size() const { return (int)y.size(); }
+
+    int ToIndex(float x) const {
+        return static_cast<int>((x-x0)/dx);
+    }
+
+    float ToXCenter(int index) const {
+        return x0 + (index+0.5f)*dx;
+    }
+    float ToX(int index) const {
+        return x0 + index*dx;
+    }
+
+    // Evaluate curve at point 'x'. Use linear interpolation if 'x' falls
+    // between two bins.
+    float At(float x) const;
+
+    void Write(const std::string &fname);
+    static Curve Read(const std::string &fname);
+
+    void FilterGauss(const Curve &source, float sigma);
+};
+
+
+inline Curve FilterGauss(int nbins, const Curve &source, float sigma) {
+    Curve c(nbins, source.x0, source.x1);
+    c.FilterGauss(source, sigma);
+    return c;
+}
+
+float Integrate(const Curve &c, float x0, float x1, float *E = NULL);
+
+inline float Integrate(const Curve &c) {
+    return Integrate(c, c.x0, c.x1);
+}
+
+
+// -------------------- Toroidal distance --------------------
+
+inline static float Sqr(float x) {
+    return x*x;
+}
+  
+// Compute toroidal distance between points a and b. Both a and b must
+// lie inside the unit square.
+inline static float DistSqrToroidal(float const *a, float const *b) {
+    float dx = a[0]-b[0], dy = a[1]-b[1];
+    return Sqr(dx + (dx < -.5f) - (dx > .5f)) + 
+        Sqr(dy + (dy < -.5f) - (dy > .5f));
+}
+
+inline static float DistToroidal(float const *a, float const *b) {
+    return sqrtf(DistSqrToroidal(a, b));
+}
+
+// -------------------- RDF calculations --------------------
+
+void CalcRDF(Curve &c, size_t npts, const float *pts, 
+        float (*distfunc)(const float *, const float *) = DistToroidal);
+Curve CalcRDF(int numbins, size_t npts, const float *pts, 
+        float sigma = 0, float (*distfunc)(const float *, const float *) = DistToroidal);
+
+Curve RDF2Power(int npts, const Curve &rdf, int nbins, float x0, float x1);
+Curve Power2RDF(int npts, const Curve &power, int nbins, float x0, float x1, float smoothing);
+
+#endif
diff --git a/targetrdf/param.cc b/targetrdf/param.cc
new file mode 100644
index 0000000..1c1201f
--- /dev/null
+++ b/targetrdf/param.cc
@@ -0,0 +1,138 @@
+//
+// Source code for the paper
+// 
+// D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
+// Controlled Aliasing", ACM Trans. Graph., 2013, in press
+//
+// Copyright (C) 2012,2013 Daniel Heck and Thomas Schloemer
+//
+#include "param.hh"
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <cctype>
+#include <algorithm>
+#include <string>
+
+
+static bool IsBool(std::string str) {
+    std::transform(str.begin(), str.end(), str.begin(), ::tolower);
+    return str == "true" || str == "false";
+}
+
+static bool ToBool(std::string str) {
+    std::transform(str.begin(), str.end(), str.begin(), ::tolower);
+    return str == "true";
+}
+
+
+ParamList::ParamList() {
+}
+
+Param *ParamList::Find(const std::string &key) {
+    for (int i=0; i<(int)list.size(); ++i) {
+        if (list[i].name == key) {
+            list[i].used = true;
+            return &list[i];
+        }
+    }
+    return NULL;
+}
+
+Param *ParamList::Define(const std::string &key, const std::string &dflt) {
+    Param *p = Insert(key);
+    p->value = dflt;
+    p->set = false;
+    return p;
+}
+
+Param *ParamList::Set(const std::string &key, const std::string &val) {
+    Param *p = Insert(key);
+    p->value = val;
+    p->set = true;
+    return p;
+}
+
+Param *ParamList::Insert(const std::string &key) {
+    if (Param *p = Find(key))
+        return p;
+    else {
+        Param param;
+        param.name = key;
+        param.set = false;
+        param.used = false;
+        list.push_back(param);
+        return &list[list.size()-1];
+    }
+}
+
+void ParamList::Parse(int argc, char **argv, std::vector<std::string> &args) {
+    for (int index=1; index < argc; index++) {
+        std::string arg = argv[index];
+        size_t eqpos = arg.find("=");
+        size_t beg = arg.find_first_not_of("-");
+        std::string name = arg.substr(beg, eqpos-beg);
+        if (eqpos != std::string::npos) {
+            std::string val = arg.substr(eqpos+1, std::string::npos);
+            Set(name, val);
+        } else if (!arg.empty() && arg[0] == '-') {
+            Param *p = Find(name);
+            if (p && IsBool(p->value)) {
+                Set(name, "true");
+            } else {
+                std::string next = ++index < argc ? argv[index] : "";
+                Set(name, next);
+            }
+        } else
+            args.push_back(arg);
+    }
+}
+
+float ParamList::GetFloat(const std::string &key) {
+    if (const Param *p = Find(key))
+        return atof(p->value.c_str());
+    abort();
+}
+
+float ParamList::GetFloat(const std::string &key, float dflt) {
+    if (const Param *p = Find(key))
+        return p->set ? atof(p->value.c_str()) : dflt;
+    return dflt;
+}
+
+int ParamList::GetInt(const std::string &key, int dflt) {
+    if (const Param *p = Find(key))
+        return atoi(p->value.c_str());
+    return dflt;
+}
+
+std::string ParamList::GetString(const std::string &key, std::string dflt) {
+    if (const Param *p = Find(key))
+        return p->value;
+    return dflt;
+}
+
+bool ParamList::GetBool(const std::string &key, bool dflt) {
+    if (const Param *p = Find(key))
+        return ToBool(p->value);
+    return dflt;
+}
+
+const Param *ParamList::UnusedOption() const {
+    for (size_t i=0; i<list.size(); i++) {
+        if (!list[i].used && list[i].set)
+            return &list[i];
+    }
+    return NULL;
+}
+
+void ParamList::Print() const {
+    for (size_t i=0; i<list.size(); i++) {
+        const Param *p = &list[i];
+        fprintf(stderr, "    %s", p->name.c_str());
+        if (!p->value.empty()) 
+            fprintf(stderr, " (%s)\n", p->value.c_str());
+        else 
+            fprintf(stderr, "\n");
+    }
+}
diff --git a/targetrdf/param.hh b/targetrdf/param.hh
new file mode 100644
index 0000000..96d04df
--- /dev/null
+++ b/targetrdf/param.hh
@@ -0,0 +1,53 @@
+//
+// Source code for the paper
+// 
+// D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
+// Controlled Aliasing", ACM Trans. Graph., 2013, in press
+//
+// Copyright (C) 2012,2013 Daniel Heck and Thomas Schloemer
+//
+#ifndef PARAM_HH_INCLUDED
+#define PARAM_HH_INCLUDED
+
+#include <exception>
+#include <string>
+#include <vector>
+
+// Command line parsing. Parameters can be specified on the command 
+// line in one of the following forms:
+//
+//  --param value
+//  --param=value
+//  param=value
+
+struct Param {
+    std::string name;
+    std::string value;
+    bool set;                   // was specified on command line
+    bool used;                  // was queried by program
+};
+
+class ParamList {
+    std::vector<Param> list;
+public:
+    ParamList();
+
+    Param *Define(const std::string &key, const std::string &dflt);
+    void Parse(int argc, char **argv, std::vector<std::string> &remaining);
+
+    float GetFloat(const std::string &key);
+    float GetFloat(const std::string &key, float dflt);
+    int GetInt(const std::string &key, int dflt = 0);
+    std::string GetString(const std::string &key, std::string dflt = "");
+    bool GetBool(const std::string &key, bool dflt = false);
+
+    const Param *UnusedOption() const;
+
+    void Print() const;
+private:
+    Param *Set(const std::string &key, const std::string &val);
+    Param *Find(const std::string &key);
+    Param *Insert(const std::string &key);
+};
+
+#endif
diff --git a/targetrdf/targetrdf.cc b/targetrdf/targetrdf.cc
new file mode 100644
index 0000000..1cc72c2
--- /dev/null
+++ b/targetrdf/targetrdf.cc
@@ -0,0 +1,414 @@
+//
+// Source code for the paper
+// 
+// D. Heck and T. Schloemer and O. Deussen, "Blue Noise Sampling with
+// Controlled Aliasing", ACM Trans. Graph., 2013, in press
+//
+// Copyright (C) 2012,2013 Daniel Heck and Thomas Schloemer
+//
+#include "param.hh"
+#include "common.hh"
+#include <iostream>
+#include <fstream>
+#include <cstdlib>
+#include <algorithm>
+#include <cfloat>
+#include <signal.h>
+#include <sys/time.h>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+int nbins = -1;                 // default: number of points
+bool randomize_force = false;
+int maxattempts = 5;
+float Tmax = 1, Tmin = 1e-3, Tfactor = 0.9;
+
+// A flag set on SIGKILL to stop optimization
+bool abort_flag = false;
+
+float CalcGradients(const std::vector<Point> &pts, 
+        const Curve &rdf, const Curve &target,
+        std::vector<Point> &gradient) {
+    Curve force(nbins, 0, 0.5f);
+
+    for (int j=0; j<force.size(); j++)
+        force[j] = force.dx * (rdf[j] - target[j]);
+    for (int j=1; j<force.size(); j++)
+        force[j] += force[j-1];
+    for (int j=1; j<force.size(); j++) {
+        float x = force.ToX(j);
+        force[j] /= (pts.size() * x * x);
+        // if (force[j] < 0) force[j] = 0;
+    }
+    force.Write("force.dat");
+
+    float maxforce = 0.0f;
+    gradient.resize(pts.size());
+    
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+{
+#ifdef _OPENMP
+#pragma omp for schedule(static)
+#endif
+    for (unsigned int i=0; i<pts.size(); i++) {
+        Point grad;
+        for (unsigned int j = 0; j < pts.size(); j++) {
+            if (i == j)
+                continue;
+
+            float dx = pts[j].x - pts[i].x;
+            float dy = pts[j].y - pts[i].y;
+            dx += (dx < -0.5f) - (dx > 0.5f);
+            dy += (dy < -0.5f) - (dy > 0.5f);
+            float dist2 = dx * dx + dy * dy;
+
+            // RDF is only reliable up to 0.5, if we go higher, the periodic
+            // boundary conditions cause anisotropies in the final spectrum
+            if (dist2 > 0.49f * 0.49f || dist2 == 0) continue;
+
+            float dist = sqrtf(dist2);
+            float f = force.At(dist);// / dist;
+            grad.x -= dx * f;
+            grad.y -= dy * f;
+        }
+        float ff = grad.x*grad.x + grad.y*grad.y;
+        if (ff > maxforce)
+            maxforce = ff;
+        gradient[i] = grad;
+    }
+}
+    return sqrtf(maxforce);
+}
+
+std::vector<Point> MovePoints(const std::vector<Point> &pts,
+        const std::vector<Point> &gradient, float stepsize)
+{
+    std::vector<Point> newpts(pts.size());
+    for (unsigned i=0; i<pts.size(); i++) {
+        Point force = gradient[i];
+        if (randomize_force) {
+            force.x *= drand48();
+            force.y *= drand48();
+        }
+        float x = pts[i].x + stepsize * force.x;
+        float y = pts[i].y + stepsize * force.y;
+
+        while (x >= 1) x -= 1;
+        while (y >= 1) y -= 1;
+        while (x < 0) x += 1;
+        while (y < 0) y += 1;
+
+        newpts[i].x = x;
+        newpts[i].y = y;
+    }
+    return newpts;
+}
+
+
+float CalcEnergy(const Curve &a, const Curve &b) {
+    Curve c(a);
+    for (int i=0; i<c.size(); i++) 
+        c[i] = powf(c[i] - b[i], 2.0f);
+    return sqrtf(Integrate(c));
+//    fprintf(stderr, "e = %g\n", x);
+//    return x;
+}
+
+float MainOptimization(const std::vector<Point> &pts,
+        Curve &target,
+        float smoothing,
+        std::vector<Point> &output)
+{
+
+    std::vector<Point> current = pts, best = pts;
+    std::vector<Point> gradient;
+    Curve rdf = CalcRDF(nbins, current.size(), &current[0].x, smoothing);
+    float bestenergy = CalcEnergy(rdf, target);
+    float T = Tmax;             // temperature
+    int attempts = 0;
+    while (!abort_flag && T >= Tmin) {
+        // Calculate gradients and move points
+        float maxgrad = CalcGradients(current, rdf, target, gradient);
+        float stepsize = T / (sqrt(current.size()) * maxgrad);
+        current = MovePoints(current, gradient, stepsize);
+
+        // Perform iteration control
+        rdf = CalcRDF(nbins, current.size(), &current[0].x, smoothing);
+        float energy = CalcEnergy(rdf, target);
+        attempts++;
+        if (energy < bestenergy) {
+            attempts = 0;
+            best = current;
+            bestenergy = energy;
+        } else if (energy > bestenergy * 1.2) {
+            attempts = maxattempts;
+            current = best;
+        }
+        printf("%1.5g %1.5g\t%5.5g\t%g\n", 
+                T, stepsize, energy, bestenergy);
+        if (attempts >= maxattempts) {
+            attempts = 0;
+            T *= Tfactor;
+        }
+    }
+    output = best;
+    return bestenergy;
+}
+
+
+void FunctionStep(float critx, Curve &c, int npts) {
+    float x0 = critx * sqrt(2/(sqrt(3)*npts));
+    for (int i=0; i<c.size(); i++) {
+        c[i] = (c.ToX(i) >= x0) ? 1 : 0;
+    }
+}
+
+void FunctionJinc(float critx, Curve &c, int npts) {
+    float x0 = critx / sqrt(2/(sqrt(3)*npts));
+    for (int i=0; i<c.size(); i++) {
+        float xx0 = c.ToX(i);
+    
+        float a = 0.0;
+        const int N=20;
+        for (int j=0; j<N; j++) {
+            float x = xx0 + j*c.dx/N;
+            if (x < 1e-5) a += 1-M_PI*x0*x0/npts;
+            else {
+                float xx = 2*M_PI*x0*x;
+                a += 1 - 2*M_PI*x0*x0/npts * jn(1, xx) / xx;
+            }
+        }
+        c[i] = a / N;
+    }
+}
+
+void FunctionPeak(float critx, float peaky, Curve &c, int npts) {
+    for (int i=0; i<c.size(); i++) {
+        c[i] = (c.ToX(i) > critx) ? 1 : 0;
+        c[i] = ((int) c.ToX(i) == (int) critx) ? peaky : c[i];
+    }
+}
+
+
+void Usage() {
+    std::cerr << "Usage: targetrdf [options]\n"
+              << "  --help                          show this message\n"
+              << "  --out file                      output file\n"
+              << "  --dda file                      output RDF in DDA form, Zhou et al.-style"
+              << "\nTarget RDF\n"
+              << "  --reference pointset            determine RDF from point set\n"
+              << "  --steprp                        generate step power\n"
+              << "  --steprdf                       generate step RDF\n"
+              << "  --peakrp                        generate peak power\n"
+              << "  --crit freq\n"
+              << "  --peak power                    height of peak\n"
+              << "  --peaksmooth sigma (default 0)  Gaussian applied to peak power\n"
+              << "  --rdf file                      load desired RDF from file\n"
+              << "  --spectrum file                 load desired spectrum from file\n"
+              << "  --rp-as-rdf file                load spectrum from file, use as RDF\n"
+              << "\nOptimization Parameters\n"
+              << "  --nbins n (default #points)     histogram size for RDF\n"
+              << "  --smoothing sigma (default 8)   Gaussian applied to all RDFs\n"
+              << "  --randomize                     randomly scale force\n"
+              << "\nInput RDF\n"
+              << "  --in file                       input point set\n"
+              << "  --npts n (default 4096)         number of random points\n"
+              << "  --seed num                      random seed\n"
+        ;
+}
+
+static void sighandler(int sig) {
+    fprintf(stderr, "Aborting...\n");
+    abort_flag = true;
+}
+
+enum AnalyticCurve { 
+    RDF_Step,
+    RDF_Jinc,
+    RP_Peak
+};
+
+int main(int argc, char **argv) {
+    ParamList params;
+    params.Define("in", "");
+    params.Define("out", "targetrdf.txt");
+    params.Define("dda", "targetdda.pfm");
+    params.Define("reference", "");
+    params.Define("rdf", "");
+    params.Define("rp-as-rdf", "");
+    params.Define("spectrum", "");
+    params.Define("seed", "1");
+    params.Define("nbins", "-1");
+    params.Define("smoothing", "8");
+    params.Define("randomize", "false");
+    params.Define("npts", "4096");
+    params.Define("steprp", "false");
+    params.Define("steprdf", "false");
+    params.Define("peakrp", "false");
+    params.Define("crit", "");
+    params.Define("peak", "");
+    params.Define("peaksmooth", "0");
+    params.Define("help", "false");
+
+    std::vector<std::string> args;
+    params.Parse(argc, argv, args);
+
+    std::string infile = params.GetString("in");
+    std::string outfile = params.GetString("out");
+    std::string ddafile = params.GetString("dda");
+    std::string reference = params.GetString("reference");
+    std::string rdffile = params.GetString("rdf");
+    std::string spectrum = params.GetString("spectrum");
+    std::string rp_as_rdf = params.GetString("rp-as-rdf");
+    nbins = params.GetInt("nbins");
+    randomize_force = params.GetBool("randomize");
+    int npts = params.GetInt("npts");
+    long int seed = params.GetInt("seed");
+    float smoothing = params.GetFloat("smoothing");
+
+    AnalyticCurve analyticCurve = RDF_Jinc;
+    float critfreq = 0.0f, peakpower = 1.0f, peaksmooth = 0.0f;
+
+    if (params.GetBool("steprp")) {
+        analyticCurve = RDF_Jinc;
+        critfreq = params.GetFloat("crit", 0.606f);
+    }
+    if (params.GetBool("steprdf")) {
+        analyticCurve = RDF_Step;
+        critfreq = params.GetFloat("crit", 1.0f);
+    }
+    if (params.GetBool("peakrp")) {
+        analyticCurve = RP_Peak;
+        critfreq = params.GetFloat("crit", 36.5f);
+        peakpower = params.GetFloat("peak", 1.0f);
+        peaksmooth = params.GetFloat("peaksmooth", 0.0f);
+    }
+    bool show_usage = params.GetBool("help");
+
+    if (!args.empty()) {
+        std::cerr << "Unexpected argument '" << args[0] << "'\n";
+        show_usage = true;
+    }
+    if (const Param *p = params.UnusedOption()) {
+        std::cerr << "Unknown option '" << p->name << "'\n";
+        show_usage = true;
+    }
+
+    if (!reference.empty() && !rdffile.empty()) {
+        std::cerr << "--rdf and --reference are mutually exclusive.\n";
+        show_usage = true;
+    }
+  
+    if (show_usage) {
+        Usage();
+        exit(1);
+    }
+
+    std::vector<Point> pts, result;
+    srand48(seed);
+    if (infile.empty()) {
+        fprintf(stderr, "No input file specified, using %d random points\n", npts);
+        pts.resize(npts);
+        for (int i=0; i<npts; i++) {
+            pts[i].x = drand48();
+            pts[i].y = drand48();
+        }
+    } else
+        LoadPoints(infile, pts);
+    result.resize(pts.size());
+
+    if (nbins == -1)
+        nbins = pts.size();
+
+    printf("Using %d bins, smoothing=%g\n", nbins, smoothing);
+
+    // Prepare target RDF
+    Curve target(nbins, 0, 0.5f);
+    if (!reference.empty()) {
+        std::vector<Point> referencepts;
+        LoadPoints(reference, referencepts);
+        target = CalcRDF(nbins, referencepts.size(), &referencepts[0].x, 
+                smoothing);
+    } else if (!rdffile.empty()) {
+        Curve tmp = Curve::Read(rdffile.c_str());
+        target = FilterGauss(nbins, tmp, smoothing);
+    } else if (!spectrum.empty()) {
+        Curve tmp  = Curve::Read(spectrum.c_str());
+        target = Power2RDF(npts, tmp, nbins, 0, 0.5f, smoothing);
+    } else if (!rp_as_rdf.empty()) {
+        Curve tmp = Curve::Read(rp_as_rdf.c_str());
+        for (int i=0; i<nbins; i++) {
+            int j = tmp.ToIndex(target.ToX(i) * npts);
+            if (j >= 0 && j < tmp.size())
+                target[i] = tmp[j];
+            else
+                target[i] = 1;
+        }
+    } else if (analyticCurve == RDF_Jinc) {
+        FunctionJinc(critfreq, target, pts.size());
+        Curve tmp(nbins, 0, npts / 2);
+        tmp = RDF2Power(npts, target, nbins, 0, npts / 2);
+        tmp.Write("target_rp.dat");
+    } else if (analyticCurve == RDF_Step) {
+        FunctionStep(critfreq, target, pts.size());
+    } else if (analyticCurve == RP_Peak) {
+        Curve tmp(nbins, 0, npts / 2);
+        FunctionPeak(critfreq, peakpower, tmp, pts.size());
+        if (peaksmooth > 0.f)
+            tmp = FilterGauss(nbins, tmp, peaksmooth);
+        tmp.Write("target_rp.dat");
+        target = Power2RDF(npts, tmp, nbins, 0, 0.5f, smoothing);
+    }
+    target.Write("target_rdf.dat");
+    printf("Target RDF written to 'target_rdf.dat'\n");
+    
+    // std::string fname("cccvt_dda.pfm");
+    // int width, height;
+    // float* data = NULL;
+    // ReadPFM(width, height, data, fname.c_str());
+    // printf("Read %s of size %dx%d\n", fname.c_str(), width, height);
+    // for (int x = 0; x < width; ++x)
+    //     for (int y = 0; y < height; ++y)
+    //         printf("%d, %d: %f\n", x, y, data[x + y*width]);
+    // if (data) delete[] data;
+    
+    if (!ddafile.empty()) {
+        // 128, 0.125
+        // 256, 0.175
+        // 512, 0.225
+        const int size = 256;
+        const float range = 0.175f;  // Guessed empirically
+        std::vector<float> dda(size * size);
+        for (int x = 0; x < size; ++x) {
+            for (int y = 0; y < size; ++y) {
+                float dx = (x - size/2) * range / (size/2);
+                float dy = (y - size/2) * range / (size/2);
+                float dist = sqrtf(dx*dx + dy*dy);
+                dda[x + y*size] = target.At(dist) / (size/2); // Scaling by size/2 guessed
+            }
+        }
+        fprintf(stderr, "Writing '%s'\n", ddafile.c_str());
+        WritePFM(size, size, &dda[0], ddafile.c_str());
+    }
+
+    fprintf(stderr, "Starting optimization\n");
+    signal(SIGINT, sighandler);
+
+    result = pts;
+    timeval t0, t;
+    
+    gettimeofday(&t0, NULL);
+    float bestd = MainOptimization(pts, target, smoothing, result);
+    gettimeofday(&t, NULL);
+    
+    long long delta = (long long)(t.tv_sec  - t0.tv_sec ) * 1000ll +
+                      (long long)(t.tv_usec - t0.tv_usec) / 1000ll;
+    fprintf(stderr, "- finished: %5.5g, time: %3.fs -\n", bestd,
+            delta / 1000.0);
+    fprintf(stderr, "Writing '%s'\n", outfile.c_str());
+    WritePoints(outfile, result);
+}