Workletize the logistic map.
This commit is contained in:
parent
57e8750b33
commit
eac6c21c84
@ -13,51 +13,90 @@
|
|||||||
#include <vector>
|
#include <vector>
|
||||||
#include <vtkm/cont/DataSetBuilderUniform.h>
|
#include <vtkm/cont/DataSetBuilderUniform.h>
|
||||||
#include <vtkm/io/ImageWriterPNG.h>
|
#include <vtkm/io/ImageWriterPNG.h>
|
||||||
|
#include <vtkm/worklet/WorkletMapField.h>
|
||||||
|
|
||||||
|
|
||||||
|
// The logistic map is xᵢ₊₁ = rxᵢ(1-xᵢ).
|
||||||
|
// If we start this iteration out at (say) x₀ = 0.5, the map has "transients",
|
||||||
|
// which we must iterate away to produce the final image.
|
||||||
|
struct LogisticBurnIn : public vtkm::worklet::WorkletMapField
|
||||||
|
{
|
||||||
|
using ControlSignature = void(WholeArrayOut);
|
||||||
|
using ExecutionSignature = void(_1, WorkIndex);
|
||||||
|
|
||||||
|
template <typename OutputArrayPortalType>
|
||||||
|
VTKM_EXEC void operator()(OutputArrayPortalType& outputArrayPortal, vtkm::Id workIndex) const
|
||||||
|
{
|
||||||
|
vtkm::Id width = outputArrayPortal.GetNumberOfValues();
|
||||||
|
double rmin = 2.9;
|
||||||
|
double r = rmin + (4.0 - rmin) * workIndex / (width - 1);
|
||||||
|
double x = 0.5;
|
||||||
|
// 2048 should be enough iterations to get rid of the transients:
|
||||||
|
int n = 0;
|
||||||
|
while (n++ < 2048)
|
||||||
|
{
|
||||||
|
x = r * x * (1 - x);
|
||||||
|
}
|
||||||
|
outputArrayPortal.Set(workIndex, x);
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
// After burn-in, the iteration is periodic but in general not convergent,
|
||||||
|
// i.e., for large enough i, there exists an integer p > 0 such that
|
||||||
|
// xᵢ₊ₚ = xᵢ for all i.
|
||||||
|
// So color the pixels corresponding to xᵢ, xᵢ₊₁, .. xᵢ₊ₚ.
|
||||||
|
struct LogisticLimitPoints : public vtkm::worklet::WorkletMapField
|
||||||
|
{
|
||||||
|
using ControlSignature = void(WholeArrayIn, WholeArrayOut);
|
||||||
|
using ExecutionSignature = void(_1, _2, WorkIndex);
|
||||||
|
|
||||||
|
template <typename InputArrayPortalType, typename OutputArrayPortalType>
|
||||||
|
VTKM_EXEC void operator()(const InputArrayPortalType& inputArrayPortal,
|
||||||
|
OutputArrayPortalType& outputArrayPortal,
|
||||||
|
vtkm::Id workIndex) const
|
||||||
|
{
|
||||||
|
vtkm::Id width = inputArrayPortal.GetNumberOfValues();
|
||||||
|
double x = inputArrayPortal.Get(workIndex);
|
||||||
|
double rmin = 2.9;
|
||||||
|
double r = rmin + (4.0 - rmin) * workIndex / (width - 1);
|
||||||
|
vtkm::Vec4f orange(1.0, 0.5, 0.0, 0.0);
|
||||||
|
|
||||||
|
// We can't display need more limit points than pixels of height:
|
||||||
|
vtkm::Id limit_points = 0;
|
||||||
|
vtkm::Id height = 1800;
|
||||||
|
while (limit_points++ < height)
|
||||||
|
{
|
||||||
|
vtkm::Id j = vtkm::Round(x * (height - 1));
|
||||||
|
outputArrayPortal.Set(j * width + workIndex, orange);
|
||||||
|
x = r * x * (1 - x);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
|
||||||
int main()
|
int main()
|
||||||
{
|
{
|
||||||
size_t height = 1800;
|
vtkm::Id height = 1800;
|
||||||
size_t width = height * 1.618;
|
vtkm::Id width = height * 1.618;
|
||||||
vtkm::cont::DataSetBuilderUniform dsb;
|
vtkm::cont::DataSetBuilderUniform dsb;
|
||||||
vtkm::cont::DataSet ds = dsb.Create(vtkm::Id2(width, height));
|
vtkm::cont::DataSet ds = dsb.Create(vtkm::Id2(width, height));
|
||||||
|
|
||||||
std::vector<double> x(width, 0.5);
|
vtkm::cont::ArrayHandle<double> x;
|
||||||
|
x.Allocate(width);
|
||||||
|
vtkm::cont::Invoker invoke;
|
||||||
|
invoke(LogisticBurnIn{}, x);
|
||||||
|
|
||||||
double rmin = 2.9;
|
vtkm::cont::ArrayHandle<vtkm::Vec4f> pixels;
|
||||||
for (size_t i = 0; i < width; ++i)
|
pixels.Allocate(width * height);
|
||||||
|
auto wp = pixels.WritePortal();
|
||||||
|
for (vtkm::Id i = 0; i < pixels.GetNumberOfValues(); ++i)
|
||||||
{
|
{
|
||||||
double r = rmin + (4.0 - rmin) * i / (width - 1);
|
wp.Set(i, vtkm::Vec4f(0, 0, 0, 0));
|
||||||
int n = 0;
|
|
||||||
// 2048 should be enough iterations to be "converged";
|
|
||||||
// though of course the iterations actually don't all converge but cycle or are chaotic.
|
|
||||||
while (n++ < 2048)
|
|
||||||
{
|
|
||||||
x[i] = r * x[i] * (1 - x[i]);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
vtkm::Vec4f v(1.0, 0.5, 0.0, 0.0);
|
|
||||||
std::vector<vtkm::Vec4f> pixelValues(width * height, vtkm::Vec4f(0, 0, 0, 0));
|
|
||||||
size_t iterates = 0;
|
|
||||||
// We don't need more iterates than pixels of height,
|
|
||||||
// by the pigeonhole principle.
|
|
||||||
while (iterates++ < height)
|
|
||||||
{
|
|
||||||
for (size_t i = 0; i < width; ++i)
|
|
||||||
{
|
|
||||||
double r = rmin + (4.0 - rmin) * i / (width - 1);
|
|
||||||
double y = x[i];
|
|
||||||
assert(y >= 0 && y <= 1);
|
|
||||||
size_t j = std::round(y * (height - 1));
|
|
||||||
pixelValues[j * width + i] = v;
|
|
||||||
x[i] = r * x[i] * (1 - x[i]);
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
invoke(LogisticLimitPoints{}, x, pixels);
|
||||||
std::string colorFieldName = "pixels";
|
std::string colorFieldName = "pixels";
|
||||||
ds.AddPointField(colorFieldName, pixelValues);
|
ds.AddPointField(colorFieldName, pixels);
|
||||||
std::string filename = "logistic.png";
|
std::string filename = "logistic.png";
|
||||||
vtkm::io::ImageWriterPNG writer(filename);
|
vtkm::io::ImageWriterPNG writer(filename);
|
||||||
writer.WriteDataSet(ds, colorFieldName);
|
writer.WriteDataSet(ds, colorFieldName);
|
||||||
|
Loading…
Reference in New Issue
Block a user