clara/bluenoise-raytracer
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

included targetrdf

Browse source
This commit is contained in:
CDaut 2024-05-21 12:16:51 +02:00
parent 39ac0cf6c7
commit 3ea1becd6f
Signed by: clara
GPG key ID: 223391B52FAD4463
9 changed files with 1276 additions and 2 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.idea/.name generated
View file

@ -1 +1 @@
vk_raytracing_tutorial
targetrdf

4
.idea/misc.xml generated
View file

@ -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>

22
targetrdf/CMakeLists.txt Normal file
View file

@ -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})

108
targetrdf/README.md Normal file
View file

@ -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

399
targetrdf/common.cc Normal file
View file

@ -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);
}

138
targetrdf/common.hh Normal file
View file

@ -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

138
targetrdf/param.cc Normal file
View file

@ -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");
}
}

53
targetrdf/param.hh Normal file
View file

@ -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

414
targetrdf/targetrdf.cc Normal file
View file

@ -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);
}