stringart/main.cpp

/**
 * \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++ -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 <vector>
#include <unordered_map>
#include "inttypes.h"

// which basic nail placement algorithms should be used
// #define grid
#define multicircle
// #define circle

class Main {
 private:
  /// last position (number of nail)
  int16_t m_lastPosition = 0;

  /// definition of a point for our dwarn 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;

  uint16_t *m_usedpaths;

  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) {
    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 Wu’s 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) {
    /// 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);
    /*if (lttIter == m_linesFromSource.end()) {
      printf("itt fail %i %i\n", src, dst);
      abort();
    }*/
    // 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;
  }

 public:
  /// \brief main function
  /// \param resolutionX X resolution of internal image
  /// \param resolutionY Y resolution of internal image
  /// \param nail vector with nail positions
  /// \param maxIter maximal number of iterations to run
  /// \param duplicateFactor duplication penality factor
  /// \param lineColor line color to use
  int run(const char* imageName, Magick::Image* img, uint16_t resolutionX, uint16_t resolutionY, int16_t requestedNumberOfNails, std::vector<uint32_t> nails, uint16_t maxIter, float duplicateFactor, uint8_t lineColor) {
    m_duplicateFactor = duplicateFactor;
    int16_t numberOfNails = nails.size();
    printf("res: %ix%i nails: %i maxIter: %i duplicatePenalty %.1f color: %i\n", resolutionX, resolutionY, numberOfNails, maxIter, duplicateFactor, lineColor);

    // for time measurement
    struct timeval tv1, tv2;
    gettimeofday(&tv1, NULL);


    for (uint16_t src = 0; src < numberOfNails; ++src) {
      for (uint16_t dst = src + 1; dst < numberOfNails; ++dst) {
        td_pointsInLine pointsInLine = drawAALine(nails[src] >> 16, nails[src] & 0xFFFF, nails[dst] >> 16, 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);
    m_numberOfNails = numberOfNails;
    // a lookup of used paths to count repeats
    m_usedpaths = reinterpret_cast<uint16_t*>(malloc(numberOfNails * numberOfNails * sizeof(uint16_t)));
    bzero(m_usedpaths, numberOfNails * numberOfNails * sizeof(uint16_t));

    /// 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));
    /// temp storage to save (current) best version, will be continously updated
    int16_t* bestState = 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;
      bestState[i] = 255;
    }
    // current iteration
    uint32_t iter = 0;
    // list of used nails with their counter
    uint8_t usedPins[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 %li\n", totalDiff);
    while ((iter < maxIter) && jumps*2 < 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;
      /// compensated diff includes penality when reusing paths
      int64_t compensatedBestDiff = INT64_MAX;
      int16_t bestTarget = -1;
      // printf("source %i\n", m_lastPosition); fflush(stdout);
      for (int16_t target = 0; target < numberOfNails; ++target) {
        if (target == m_lastPosition) continue;
        int64_t diff = checkLine(target);
        if (diff < bestDiff) {
          bestTarget = target;
          bestDiff = diff;
        }
      }
      // printf("bestTarget %i diff %li\n", bestTarget, bestDiff);
      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));
        /*if (lttIter == m_linesFromSource.end()) {
          printf("iter fail %i %i\n", m_lastPosition, bestTarget);
          abort();
        }*/
        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;
          int16_t cur = m_currentState[index];
          cur -= sub;
          bestState[index] = cur;
          ++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() % 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];
      // progress report
      if (iter % 100 == 0) {
        printf("best %4i -> %4i(%4i, %4i) diff %9li (%12li) iter %5i path %i\n", m_lastPosition, bestTarget, nails[bestTarget] >> 16, nails[bestTarget] & 0xFFFF, compensatedBestDiff, totalDiff + bestDiff, iter, m_usedpaths[std::min(m_lastPosition, bestTarget) * m_numberOfNails + std::max(m_lastPosition, bestTarget)]);
      }
      // set new start position
      m_lastPosition = bestTarget;
      // update diff
      totalDiff += bestDiff;
      // update current state from best map
      memcpy(m_currentState, bestState, widthXheight * 2);
    }
    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", nails[i] >> 16, nails[i] & 0xffff);
      } else {
        fprintf(fh, "L%i %i", nails[i] >> 16, 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, duplicateFactor, lineColor, imgHeight - 30, numberOfNails, path.size(), timeNeeded);
    fprintf(fh, "</svg>");
    fclose(fh);
    // cleanup
    free(m_targetState);
    free(bestState);
    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 numberOfNails = 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 < numberOfNails; ++i) {
      float x = sin(2.0 * M_PI * i / numberOfNails) * (imgWidth-1) / 2.0 + imgWidth / 2.0;
      float y = cos(2.0 * M_PI * i / numberOfNails) * (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 = numberOfNails/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 = numberOfNails/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 = numberOfNails/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 = numberOfNails/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 = numberOfNails/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 = numberOfNails / 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(numberOfNails);
  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;
  m.run(imageName, &img, resolutionX, resolutionY, numberOfNails, nails, maxIter, duplicateFactor, lineColor);

  Magick::TerminateMagick();

}