From 3ea1becd6f6c970e1adfc4222cd6148971ef8077 Mon Sep 17 00:00:00 2001 From: CDaut Date: Tue, 21 May 2024 12:16:51 +0200 Subject: [PATCH] included targetrdf --- .idea/.name | 2 +- .idea/misc.xml | 4 +- targetrdf/CMakeLists.txt | 22 +++ targetrdf/README.md | 108 ++++++++++ targetrdf/common.cc | 399 +++++++++++++++++++++++++++++++++++++ targetrdf/common.hh | 138 +++++++++++++ targetrdf/param.cc | 138 +++++++++++++ targetrdf/param.hh | 53 +++++ targetrdf/targetrdf.cc | 414 +++++++++++++++++++++++++++++++++++++++ 9 files changed, 1276 insertions(+), 2 deletions(-) create mode 100644 targetrdf/CMakeLists.txt create mode 100644 targetrdf/README.md create mode 100644 targetrdf/common.cc create mode 100644 targetrdf/common.hh create mode 100644 targetrdf/param.cc create mode 100644 targetrdf/param.hh create mode 100644 targetrdf/targetrdf.cc 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 @@ - + + + \ 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 + +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 +#include +#include +#include + +// 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 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 &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 &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 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 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 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 bins(c.size()); + for (unsigned j=0; j 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 +#include +#include +#include + +// -------------------- Math -------------------- + +template +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 &points); +void WritePoints(const std::string &fname, std::vector &points); + +// -------------------- Curve -------------------- + +class Curve { + std::vector 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((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 +#include +#include +#include +#include +#include + + +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 &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; iname.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 +#include +#include + +// 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 list; +public: + ParamList(); + + Param *Define(const std::string &key, const std::string &dflt); + void Parse(int argc, char **argv, std::vector &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 +#include +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#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 &pts, + const Curve &rdf, const Curve &target, + std::vector &gradient) { + Curve force(nbins, 0, 0.5f); + + for (int j=0; j 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 MovePoints(const std::vector &pts, + const std::vector &gradient, float stepsize) +{ + std::vector newpts(pts.size()); + for (unsigned i=0; i= 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 &pts, + Curve &target, + float smoothing, + std::vector &output) +{ + + std::vector current = pts, best = pts; + std::vector gradient; + Curve rdf = CalcRDF(nbins, current.size(), ¤t[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(), ¤t[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= x0) ? 1 : 0; + } +} + +void FunctionJinc(float critx, Curve &c, int npts) { + float x0 = critx / sqrt(2/(sqrt(3)*npts)); + for (int i=0; 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 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 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 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= 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 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); +}