stringart/main.cpp

627 lines
22 KiB
C++
Raw Permalink Blame History

This file contains ambiguous Unicode characters!

This file contains ambiguous Unicode characters that may be confused with others in your current locale. If your use case is intentional and legitimate, you can safely ignore this warning. Use the Escape button to highlight these characters.

/**
* \file main.cpp
* \brief implementation of the string art algorithm
* \author GrumpyDeveloper (Sascha Nitsch)
* \copyright 2023 Sascha Nitsch
* Licensed under GPL3 or later license
* https://contentnation.net/en/grumpydevelop/stringart
* SPDX-FileCopyrightText: 2023 Sascha Nitsch (@grumpydevelop@contentnation.net) https://contentnation.net/en/grumpydevelop
* SPDX-License-Identifier: GPL-3.0-or-later
*
*/
// compile with
/// g++ --std=c++20 -march=native`Magick++-config --cxxflags --cppflags` -Wall -Werror -o main main.cpp -O3 `Magick++-config --ldflags --libs`
#include <math.h>
#include <Magick++.h>
#include <string.h>
#include <sys/time.h>
#include <inttypes.h>
#include <semaphore>
#include <barrier>
#include <vector>
#include <unordered_map>
#include <mutex>
#include <thread>
#include <condition_variable>
// which basic nail placement algorithms should be used
// #define grid
#define multicircle
// #define circle
class Main {
private:
/// indicator, which threads are busy
uint64_t m_busyFlag;
/// number of threads
int16_t m_numThreads;
/// worker threads
std::vector<std::thread> m_worker;
/// result from checkLine thread
int64_t* m_lineResult;
/// next target to process in thread pool
int16_t m_nextTarget;
/// result pool lock
std::mutex m_resultLock;
/// mutex for work notification
std::mutex m_workMutex;
/// work notification condition
std::condition_variable m_workNotification;
/// mutex for finished motification
std::mutex m_finishedMutex;
/// finished notification condition
std::binary_semaphore m_finishedNotification{0};
/// sync point for worker threads and main
std::barrier<void(*)()> m_syncPoint;
/// last position (number of nail)
int16_t m_lastPosition = 0;
/// tasks left
int32_t m_tasksLeft;
/// running flag for threads
bool m_running;
/// list of nails
const std::vector<uint32_t>& m_nails;
/// definition of a point for our drawn line vector
struct Point {
/// x coordinate
uint16_t x;
/// y coordinate
uint16_t y;
/// color value
uint8_t color;
};
/// typedef for a list of points tha make a line
typedef std::vector<Point> td_pointsInLine;
/// the line from src to dst, key = (src << 16) + dst
typedef std::unordered_map<uint32_t, td_pointsInLine> td_linesFromSource;
/// our line storage
td_linesFromSource m_linesFromSource;
/// internal image width
uint16_t m_imgWidth;
/// current state of image
int16_t* m_currentState;
/// desired target state
uint8_t* m_targetState;
/// penalty duplication factor
float m_duplicateFactor;
/// map of used paths
uint16_t *m_usedpaths;
/// number of nails
int16_t m_numberOfNails;
/// the actual weight function to calculate how off we are to the target
/// \param value current value
/// \param target desired target
/// \retval distance to target
inline int64_t weightFunction(int16_t value, int16_t target) const {
return (value - target) * (value - target);
}
/// swaps two numbers
/// \param a first number
/// \param b second number
inline void swap(int16_t* a , int16_t* b) {
int16_t temp = *a;
*a = *b;
*b = temp;
}
/// return floating part of number
/// \param x number to process
/// \retval the data behind the dot
inline float fPartOfNumber(float x) {
if (x > 0) {
return x - floor(x);
}
return x - (floor(x) + 1);
}
/// add given point to vector if col is > 0
/// \param pil pointer to vector
/// \param x x coordinate
/// \param y y coordingte
/// \param col color
inline void pb(td_pointsInLine* pil, uint16_t x, uint16_t y, uint8_t col) {
if (col > 0) {
pil->push_back({x, y, col});
}
}
/// draw line using Xiaolin Wus line algorithm
/// \param x0 source x coordinate
/// \param y0 source y coordinate
/// \param x1 destination x coordinate
/// \param y1 destination y coordinate
/// \param color color to draw
td_pointsInLine drawAALine(int16_t x0 , int16_t y0 , int16_t x1 , int16_t y1, uint8_t color) {
td_pointsInLine pil;
bool steep = abs(y1 - y0) > abs(x1 - x0);
// swap the co-ordinates if slope > 1 or we
// draw backwards
if (steep) {
swap(&x0, &y0);
swap(&x1, &y1);
}
if (x0 > x1) {
swap(&x0, &x1);
swap(&y0, &y1);
}
// compute the slope
float dx = x1 - x0;
float dy = y1 - y0;
float gradient = dy / dx;
if (dx == 0.0) {
gradient = 1;
}
int16_t xpxl1 = x0;
int16_t xpxl2 = x1;
float intersectY = y0;
// main loop
if (steep) {
int16_t x;
for (x = xpxl1 ; x <= xpxl2 ; ++x) {
// pixel coverage is determined by fractional
// part of y co-ordinate
pb(&pil, static_cast<uint16_t>(intersectY), static_cast<uint16_t>(x), static_cast<uint8_t>(color * (1 - fPartOfNumber(intersectY))));
if (intersectY >= 1) {
pb(&pil, static_cast<uint16_t>(intersectY - 1), static_cast<uint16_t>(x), static_cast<uint8_t>(color * fPartOfNumber(intersectY)));
}
intersectY += gradient;
}
} else {
int16_t x;
for (x = xpxl1 ; x <= xpxl2 ; ++x) {
// pixel coverage is determined by fractional
// part of y co-ordinate
pb(&pil, static_cast<uint16_t>(x), static_cast<uint16_t>(intersectY), static_cast<uint8_t>(color * (1 - fPartOfNumber(intersectY))));
if (intersectY >= 1) {
pb(&pil, static_cast<uint16_t>(x), static_cast<uint16_t>(intersectY - 1), static_cast<uint8_t>(color * fPartOfNumber(intersectY)));
}
intersectY += gradient;
}
}
return pil;
}
float checkLine(int16_t target) const {
if (m_lastPosition == target) return INT64_MAX;
/// the diff on current lastPosition -> target
int64_t testDiff = 0;
uint16_t src = std::min(m_lastPosition, target);
uint16_t dst = std::max(m_lastPosition, target);
td_linesFromSource::const_iterator lttIter = m_linesFromSource.find((src << 16) + dst);
// calculate difference to target
// for each point
td_pointsInLine::const_iterator pilIter = lttIter->second.begin();
while (pilIter != lttIter->second.end()) {
uint16_t x = (*pilIter).x;
uint16_t y = (*pilIter).y;
uint8_t sub = (*pilIter).color;
uint32_t index = y * m_imgWidth + x;
int16_t cur = m_currentState[index];
int16_t goal = m_targetState[index];
// subtract previous error
testDiff -= weightFunction(cur, goal);
cur -= sub;
// add new error
testDiff += weightFunction(cur, goal);
++pilIter;
}
float duplicatePenalty = (m_duplicateFactor != 1) ?
pow(m_duplicateFactor, m_usedpaths[std::min(m_lastPosition, target)* m_numberOfNails + std::max(m_lastPosition, target)])
: 1;
return testDiff * duplicatePenalty;
}
void checkLines(int16_t tid) {
bool first = true;
while (m_running) {
m_syncPoint.arrive_and_wait();
// wait for notification
{
std::unique_lock<std::mutex> lock(m_workMutex);
m_busyFlag &= ~(1 << tid);
if (!m_busyFlag) {
if (!first) {
m_finishedNotification.release();
} else {
first = false;
}
}
m_workNotification.wait(lock);
m_busyFlag |= (1 << tid);
}
// did we wake up because of shutdown?
if (!m_running) break;
int16_t tasksLeft = 0;
int16_t target = 0;
do { // as long as there is work
{
std::unique_lock<std::mutex> lock(m_resultLock);
target = m_nextTarget++;
tasksLeft = m_tasksLeft--;
}
if (target < m_numberOfNails) {
m_lineResult[target] = checkLine(target);
}
} while (tasksLeft > 0);
}
}
public:
/// \brief constructor
/// \param nail vector with nail positions
/// \param duplicateFactor duplication penalty factor
Main(const std::vector<uint32_t>& nails, float duplicateFactor, int16_t numThreads) : m_syncPoint(numThreads + 1, [](){}), m_nails(nails) {
m_duplicateFactor = duplicateFactor;
m_numberOfNails = nails.size();
m_lineResult = reinterpret_cast<int64_t*>(malloc(sizeof(int64_t) * m_numberOfNails));
for (uint16_t i = 0; i < m_numberOfNails; ++i) {
m_lineResult[i] = INT64_MAX;
}
// a lookup of used paths to count repeats
m_usedpaths = reinterpret_cast<uint16_t*>(malloc(m_numberOfNails * m_numberOfNails * sizeof(uint16_t)));
bzero(m_usedpaths, m_numberOfNails * m_numberOfNails * sizeof(uint16_t));
m_nextTarget = 0;
m_imgWidth = 0;
m_running = true;
m_targetState = NULL;
m_currentState = NULL;
m_tasksLeft = 0;
m_numThreads = numThreads;
m_busyFlag = 0;
m_worker.reserve(numThreads);
for (uint16_t i = 0; i < numThreads; ++i) {
m_worker.push_back(std::thread(&Main::checkLines, this, i));
}
// wait until all threads are there
m_syncPoint.arrive_and_wait();
}
/// destructor
~Main() {
m_running = false;
{
std::unique_lock<std::mutex> lock(m_workMutex);
m_workNotification.notify_all();
}
for (uint16_t i = 0; i < m_numThreads; ++i) {
m_worker[i].join();
}
free(m_lineResult);
}
/// \brief main function
/// \param resolutionX X resolution of internal image
/// \param resolutionY Y resolution of internal image
/// \param maxIter maximal number of iterations to run
/// \param lineColor line color to use
int run(const char* imageName, Magick::Image* img, uint16_t resolutionX, uint16_t resolutionY, int16_t requestedNumberOfNails, uint16_t maxIter, uint8_t lineColor) {
printf("res: %ix%i nails: %i maxIter: %i duplicatePenalty %.1f color: %i\n", resolutionX, resolutionY, m_numberOfNails, maxIter, m_duplicateFactor, lineColor);
// for time measurement
struct timeval tv1, tv2;
gettimeofday(&tv1, NULL);
for (uint16_t src = 0; src < m_numberOfNails; ++src) {
for (uint16_t dst = src + 1; dst < m_numberOfNails; ++dst) {
td_pointsInLine pointsInLine = drawAALine(m_nails[src] >> 16, m_nails[src] & 0xFFFF, m_nails[dst] >> 16, m_nails[dst] & 0xFFFF, lineColor);
m_linesFromSource.insert(std::make_pair((src << 16) + dst, pointsInLine));
}
}
uint32_t channels = img->channels();
m_imgWidth = img->columns();
uint16_t imgHeight = img->rows();
printf("target image %i x %i x %i\n", m_imgWidth, imgHeight, channels);
MagickCore::Quantum *pixels = img->getPixels(0, 0, m_imgWidth, imgHeight);
m_targetState = reinterpret_cast<uint8_t*>(malloc(m_imgWidth * imgHeight));
if (channels == 1) { // monochrome image
for (uint32_t i = 0; i < m_imgWidth * imgHeight; ++i) {
m_targetState[i] = pixels[i] >> 8;
}
} else if (channels == 2) { // color + alpha?
for (uint32_t i = 0; i < m_imgWidth * imgHeight; ++i) {
m_targetState[i] = pixels[i << 1] >> 8;
}
} else { // RGB or RGBA
for (uint32_t i = 0; i < m_imgWidth * imgHeight; ++i) {
m_targetState[i] = ((pixels[i*channels] + pixels[i*channels + 1] + pixels[i*channels + 2]) / 3) >> 8;
}
}
#ifdef DEBUGIMG
FILE* debugfh = fopen("debug.pnm", "wb");
fprintf(debugfh, "P5\n# debug\n%i %i\n255\n", realWidth, realHeight);
fwrite(targetState, 1, realWidth * realHeight, debugfh);
fclose(debugfh);
#endif
// thread path
std::vector<uint16_t> path;
// add start position
path.push_back(0);
/// last thread end position
m_lastPosition = 0;
/// storage for the current state (all previous drawn threads)
m_currentState = reinterpret_cast<int16_t*>(malloc(m_imgWidth * imgHeight * 2));
// clear states
uint32_t widthXheight = m_imgWidth * imgHeight;
for (uint32_t i = 0; i < widthXheight; ++i) {
m_currentState[i] = 255;
}
// current iteration
uint32_t iter = 0;
// list of used nails with their counter
uint8_t usedPins[m_numberOfNails] = {0};
// number of continous jump tries if we got stuck
uint16_t jumps = 0;
/// total diff from currentState to targetState
int64_t totalDiff = 0;
// calculate inital difference
for (uint32_t i = 0; i < widthXheight; ++i) {
totalDiff += weightFunction(255, m_targetState[i]);
}
printf("start diff %li\n", totalDiff);
while ((iter < maxIter) && jumps*2 < m_numberOfNails) {
++iter;
#ifdef SANITYCHECK
int64_t sanity = 0;
for (uint32_t Z = 0; Z < widthXheight; ++Z) {
int16_t cur = currentState[Z];
int16_t goal = targetState[Z];
sanity += weightFunction(cur, goal);
}
if (sanity != totalDiff) {
printf("%i: total: %li, sanity: %li, diff: %li\n", iter, totalDiff, sanity, sanity - totalDiff);
}
#endif
/// current best difference
int64_t bestDiff = INT64_MAX;
int16_t bestTarget = -1;
{
std::unique_lock<std::mutex> lock(m_resultLock);
m_nextTarget = 0;
m_tasksLeft = m_numberOfNails;
}
// notify threads
{
std::unique_lock<std::mutex> lock(m_workMutex);
m_workNotification.notify_all();
}
m_syncPoint.arrive_and_wait();
// wait for threads, results are in m_lineResults
m_finishedNotification.acquire();
// search best result
for (int16_t target = 0; target < m_numberOfNails; ++target) {
if (target == m_lastPosition) continue;
if (m_lineResult[target] < bestDiff) {
bestTarget = target;
bestDiff = m_lineResult[target];
}
}
if (bestDiff < 0) {
// apply current best
td_linesFromSource::const_iterator lttIter = m_linesFromSource.find((std::min(m_lastPosition, bestTarget) << 16) + std::max(m_lastPosition, bestTarget));
td_pointsInLine::const_iterator pilIter = lttIter->second.begin();
while (pilIter != lttIter->second.end()) {
uint16_t x = (*pilIter).x;
uint16_t y = (*pilIter).y;
int16_t sub = (*pilIter).color;
uint32_t index = y * m_imgWidth + x;
m_currentState[index] -= sub;
++pilIter;
}
} else {
// we got worse, jump to random place to continue
// printf("j %3i %3i -> %3i(%4i, %4i) bestDiff %8li (%li) iter %i path %i\n", jumps, lastPosition, bestTarget, pos[bestTarget] >> 16, pos[bestTarget] & 0xFFFF, realBestDiff, totalDiff - realBestDiff, iter, m_usedpaths[std::min(m_lastPosition, bestTarget)* m_numberOfNails + std::max(m_lastPosition, bestTarget)]]);
if (jumps) { // undo last jump, was not working anyway
--usedPins[m_lastPosition];
path.pop_back();
}
// select next target randomly (kind of, intentially producing the same numbers)
bestTarget = random() % m_numberOfNails;
path.push_back(bestTarget);
m_lastPosition = bestTarget;
++usedPins[bestTarget];
++jumps;
--iter;
continue;
}
/// reset jump counter (faster to always set)
jumps = 0;
if (bestTarget < 0) {
printf("no best\n");
break;
}
// update path map
++m_usedpaths[std::min(m_lastPosition, bestTarget)* m_numberOfNails + std::max(m_lastPosition, bestTarget)];
// add new stop
path.push_back(bestTarget);
// update used pins
++usedPins[bestTarget];
// update diff
totalDiff += bestDiff;
// progress report
if (iter % 100 == 0) {
printf("best %4i -> %4i(%4i, %4i) diff %12li iter %5i path %i\n", m_lastPosition, bestTarget, m_nails[bestTarget] >> 16, m_nails[bestTarget] & 0xFFFF,
totalDiff, iter, m_usedpaths[std::min(m_lastPosition, bestTarget) * m_numberOfNails + std::max(m_lastPosition, bestTarget)]);
}
// set new start position
m_lastPosition = bestTarget;
// update current state from best map
}
printf("size %li\n", path.size());
// we are done, create output svg
FILE* fh = fopen("map.svg", "wb");
fprintf(fh, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 %i %i\">\n<rect width=\"%i\" height=\"%i\" fill=\"#ffffff\" />\n<g style=\"fill:none;stroke:#000000;stroke-opacity:%.2f;stroke-width:1\">\n", m_imgWidth, imgHeight, m_imgWidth, imgHeight, lineColor / 255.0);
uint32_t counter = 0;
for (uint16_t i : path) {
if ((counter & 255) == 0) {
fprintf(fh, "<path d=\"M%i %i", m_nails[i] >> 16, m_nails[i] & 0xffff);
} else {
fprintf(fh, "L%i %i", m_nails[i] >> 16, m_nails[i] & 0xffff);
}
if ((counter & 255) == 255) {
fprintf(fh, "\" />\n");
}
++counter;
}
if ((counter & 255) != 0) {
fprintf(fh, "\" />\n");
}
gettimeofday(&tv2, NULL);
float timeNeeded = (tv2.tv_sec - tv1.tv_sec) + (tv2.tv_usec - tv1.tv_usec) / 1000000.0;
fprintf(fh, "</g>\n");
fprintf(fh, "<text x=\"30\" y=\"10\" style=\"font-weight:bold;font-size:60px;font-family:'DejaVu Serif'\"><tspan x=\"30\" y=\"70\">%s %i %i</tspan><tspan x=\"70\" y=\"150\"> %i %i %.2f %i</tspan><tspan x=\"30\" y=\"%i\">%i nails, %li paths, %.1f sec</tspan></text>", imageName, resolutionX, resolutionY, requestedNumberOfNails, maxIter, m_duplicateFactor, lineColor, imgHeight - 30, m_numberOfNails, path.size(), timeNeeded);
fprintf(fh, "</svg>");
fclose(fh);
// cleanup
free(m_targetState);
free(m_currentState);
free(m_usedpaths);
path.clear();
m_linesFromSource.clear();
return 0;
}
};
/// \brief main entry point
/// \param argc number of command line arguments
/// \param argv command line arguments
int main(int argc, char* argv[]) {
if (argc != 8) {
printf("usage: %s <image name> <resolution x> <resolution y> <number of nails> <max number of iterations> <penalty for duplicate path usage> <lineColor>\n", argv[0]);
return 1;
}
// copy command line data to easier to use variables
const char* imageName = argv[1];
uint16_t resolutionX = atoi(argv[2]);
uint16_t resolutionY = atoi(argv[3]);
uint16_t requestedNumberOfNails = atoi(argv[4]);
uint16_t maxIter = atoi(argv[5]);
float duplicateFactor = atof(argv[6]);
uint8_t lineColor = atoi(argv[7]);
/// nail positions (x << 16) + y
std::vector<uint32_t> nails;
// initialize image magick, load image and resize to target coordinates
Magick::InitializeMagick(NULL);
Magick::Image img(imageName);
if (img.depth() != 8) {
printf("only 8 bit images supported\n");
return 1;
}
img.sample(Magick::Geometry(resolutionX, resolutionY));
// fix potential size differences between requested and delivered size
uint16_t imgWidth = img.columns();
uint16_t imgHeight = img.rows();
// position nails
#ifdef circle
for (uint16_t i = 0; i < requestedNumberOfNails; ++i) {
float x = sin(2.0 * M_PI * i / requestedNumberOfNails) * (imgWidth-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / requestedNumberOfNails) * (imgHeight-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
#endif
#ifdef multicircle
uint16_t count = requestedNumberOfNails/1.5;
for (uint16_t i = 0; i < count; ++i) {
float x = sin(2.0 * M_PI * i / count) * (imgWidth-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / count) * (imgHeight-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
uint16_t width = imgWidth / 1.2 - 1;
uint16_t height = imgHeight / 1.2 - 1;
count = requestedNumberOfNails / 1.5;
for (uint16_t i = 0; i < count; ++i) {
float x = sin(2.0 * M_PI * i / count) * (width-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / count) * (height-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
width = imgWidth / 1.5 - 1;
height = imgHeight / 1.5 - 1;
count = requestedNumberOfNails / 2;
for (uint16_t i = 0; i < count; ++i) {
float x = sin(2.0 * M_PI * i / count) * (width-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / count) * (height-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
width = imgWidth / 2 - 1;
height = imgHeight / 2 - 1;
count = requestedNumberOfNails / 3;
for (uint16_t i = 0; i < count; ++i) {
float x = sin(2.0 * M_PI * i / count) * (width-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / count) * (height-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
width = imgWidth / 3 - 1;
height = imgHeight / 3 - 1;
count = requestedNumberOfNails / 4;
for (uint16_t i = 0; i < count; ++i) {
float x = sin(2.0 * M_PI * i / count) * (width-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / count) * (height-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
width = imgWidth / 5 - 1;
height = imgHeight / 5 - 1;
count = requestedNumberOfNails / 6;
for (uint16_t i = 0; i < count; ++i) {
float x = sin(2.0 * M_PI * i / count) * (width-1) / 2.0 + imgWidth / 2.0;
float y = cos(2.0 * M_PI * i / count) * (height-1) / 2.0 + imgHeight / 2.0;
nails.push_back((static_cast<uint32_t>(floor(x)) << 16) + static_cast<uint16_t>(floor(y)));
}
nails.push_back((static_cast<uint32_t>(floor(imgWidth / 2.0)) << 16) + static_cast<uint16_t>(floor(imgHeight / 2.0)));
#endif
#ifdef grid
uint8_t sq_pins = sqrt(requestedNumberOfNails);
float distX = static_cast<float>(imgWidth - 1) / (sq_pins - 1);
float distY = static_cast<float>(imgHeight - 1) / (sq_pins - 1);
for (uint16_t y = 0; y < sq_pins; ++y) {
for (uint16_t x = 0; x < sq_pins; ++x) {
nails.push_back((static_cast<uint32_t>(floor(distX * x)) << 16) + static_cast<uint16_t>(floor(distY * y)));
}
}
#endif
Main m(nails, duplicateFactor, std::min(std::thread::hardware_concurrency(), 64U));
m.run(imageName, &img, resolutionX, resolutionY, requestedNumberOfNails, maxIter, lineColor);
Magick::TerminateMagick();
}