To be honest, I finished this implementation many months ago (in March, I think), but I've been busy with other things that I haven't had time to write about it yet.

Originally, this idea came to me thanks to a friend, and I already had code for image stippling, but it has a couple of issues:

1. It's not my own implementation (I coded it, but it was based on codetrain work)
2. It is in JavaScript

but putting those issues aside, I wondered: If I consider each of these stippling points as a node, what would be the best way to draw edges between them so that it resembles to the original image?

My image stippling code already assigns a weight to each node depending of how important it's in the whole image (thanks to voronoid stippling), so I could assume that those nodes should also have important edges (?).

This way I could assign a different width to each edge, I don't call it weight because it's the weight of the node and that weight is merely used for visualization purposes, so in this context edge width and node weight means the same.

so, since edges doesn't have a weight... I want to connect all the nodes in a graph keeping the whole covered distance as short as possible and that it's also a cycle (this condition will be explained later).

If I rephrase it a bit:

Given a list of nodes and the distances between each pair of nodes, what is the shortest possible route that visits each node exactly once and returns to the origin node?

that is the Travelling Salesman Problem, which is NP-hard :(... but before of giving up, let me explain how this could be approched if we relax the conditions a bit.

oh and btw, this will be in Rust, but I'm not so good at it yet, so... enjoy.

TSP Art

Although TSP is a NP-hard problem, I don't really care so much about this condition: shortest possible route. I prefer a shorter route, but that is not mandatory in this context, what really matters most is:

1. Visit each node exactly once
2. Returns to the origin node
3. No cross edges

last condition is important and not originally part of the TSP, but we require a planar graph so anyway the condition of shortest possible route will probably imply another kind of algorithms with this new condition.

I'll use two different strategies for the tour generation, one is the cheapest insertion algorithm that I'll use only when the amount of nodes is less than 2048, otherwise I'll use a greedy strategy, both of them uses the 2-opt algorithm to remove cross edges.

My greedy algorithm is pretty similar to the cheapest insertion algorithm, but instead of trying with every node and calculate the cost of that insertion, I'll skip some nodes and taking in consideration only the nodes that minimizes the euclidean distance as an heuristic.

impl TourStrategy for GreedyStrategy {
  fn build_tour(&self, points: &Vec<(f32, f32)>, hull_points: &Vec<Point<f32>>) -> Vec<usize> {
    let mut hull_indices: Vec<usize> = hull_points
      .into_iter()
      .map(|p| {
        points
          .iter()
          .position(|&pt| pt.0 == p.x() && pt.1 == p.y())
          .unwrap()
      })
      .collect();

    if hull_indices.len() > 1 && hull_indices.first() == hull_indices.last() {
      hull_indices.pop();
    }
    let mut tour: Vec<usize> = hull_indices.clone();
    let mut in_tour: Vec<bool> = vec![false; points.len()];
    for &idx in &tour {
      in_tour[idx] = true;
    }

    let mut kdtree = KdTree::<f32, 2>::with_capacity(points.len());
    for (idx, &(x, y)) in points.iter().enumerate() {
      if !in_tour[idx] {
        kdtree.add(&[x, y], idx as u64);
      }
    }

    while tour.len() < points.len() {
      let mut best_increase: f32 = f32::INFINITY;
      let mut best_p = None;
      let mut best_k = None;

      let step: usize = (tour.len() >> 8).max(1).min(4);
      let max_neighbors: usize = if step > 4 { 8 } else { 3 };
      for k in (0..tour.len()).step_by(step) {
        let next_k: usize = (k + 1) % tour.len();
        let i: usize = tour[k];
        let j: usize = tour[next_k];
        let p1: (f32, f32) = points[i];
        let p2: (f32, f32) = points[j];

        let mid_x: f32 = (p1.0 + p2.0) / 2.0;
        let mid_y: f32 = (p1.1 + p2.1) / 2.0;

        let neighbors =
          kdtree.nearest_n::<SquaredEuclidean>(&[mid_x, mid_y], max_neighbors);
        for neighbor in neighbors {
          let p: usize = neighbor.item as usize;
          if in_tour[p] {
            continue;
          }

          let increase: f32 = distance(points[i], points[p])
            + distance(points[p], points[j])
            - distance(points[i], points[j]);

          if increase < best_increase {
            best_increase = increase;
            best_p = Some(p);
            best_k = Some(k);
          }
        }
      }

      if let (Some(p), Some(k)) = (best_p, best_k) {
        let next_k: usize = (k + 1) % tour.len();

        if next_k == 0 {
          tour.push(p);
        } else {
          tour.insert(next_k, p);
        }

        in_tour[p] = true;

        let neighbors =
          kdtree.nearest_n::<SquaredEuclidean>(&[points[p].0, points[p].1], 1);
        for neighbor in neighbors.iter() {
          let idx = neighbor.item as usize;
          if idx == p {
            kdtree.remove(&[points[p].0, points[p].1], p as u64);
            break;
          }
        }
      } else {
        let remaining: Vec<usize> =
          (0..points.len()).filter(|&idx| !in_tour[idx]).collect();

        if !remaining.is_empty() {
          let p: usize = remaining[0];
          tour.push(p);
          in_tour[p] = true;

          kdtree.remove(&[points[p].0, points[p].1], p as u64);
        } else {
          break;
        }
      }
      if tour.len() % 32 == 0 {
        two_opt(&points, &mut tour);
      }
    }

    two_opt(&points, &mut tour);
    tour
  }
}

Fourier transform

This is actually not part of TSP art, and I didn't saw to anyone else doing something similar but, now that the code is able to find a closed path wihtout cross edges, it's also possible to describe this discrete path with Fourier epicycles. It means, you can now describe a path that draws your image or at least, kinda, using just Fourier.

pub fn compute_fourier_series(
  points: &[(f32, f32)],
  num_epicycles: usize,
) -> (Complex<f32>, Vec<Epicycle>) {
  let n_points: usize = points.len();
  let mut buffer: Vec<Complex<f32>> = points
    .iter()
    .map(|&(x, y)| Complex { re: x, im: y })
    .collect();
  let mut planner = FftPlanner::new();
  let fft = planner.plan_fft_forward(n_points);
  fft.process(&mut buffer);

  let c_0 = buffer[0] / n_points as f32;

  let mut coeffs: Vec<(i32, f32, f32)> = (1..n_points)
    .map(|k| {
      let freq = if k <= n_points / 2 {
        k as i32
      } else {
        k as i32 - n_points as i32
      };
      let coeff = buffer[k] / n_points as f32;
      let radius = coeff.norm();
      let phase = coeff.arg();
      (freq, radius, phase)
    })
    .collect();

  coeffs.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap());
  let num_epicycles = num_epicycles.min(coeffs.len());
  let selected_coeffs = &coeffs[0..num_epicycles];

  let epicycles: Vec<Epicycle> = selected_coeffs
    .iter()
    .map(|&(freq, radius, phase)| Epicycle {
      radius,
      freq,
      phase,
    })
    .collect();

  (c_0, epicycles)
}

pub fn compute_position(c_0: Complex<f32>, epicycles: &[Epicycle], t: f32) -> (f32, f32) {
  let mut x = c_0.re;
  let mut y = c_0.im;
  for epicycle in epicycles {
    let angle = 2.0 * PI * epicycle.freq as f32 * t + epicycle.phase;
    x += epicycle.radius * angle.cos();
    y += epicycle.radius * angle.sin();
  }
  (x, y)
}

Additional

Other assets I've generated with this code