Wave Function Collapse implementation in pure GDScript

import_plugin_wfc_overlapping_model_2d_params.gd

@tool
extends EditorImportPlugin

func _get_importer_name() -> String:
	return "wfc_overlapping_model_2d_params.dreamer_wfc.apples"

func _get_visible_name() -> String:
	return "WFC Overlapping Model Params"

func _get_recognized_extensions() -> PackedStringArray:
	return ["png", "bmp", "hdr", "jpg", "jpeg", "svg", "tga", "exr", "webp", "image"]

func _get_priority() -> float:
	return 1.0

func _get_save_extension() -> String:
	return "res"

func _get_resource_type() -> String:
	return "Resource"

func _get_preset_count() -> int:
	return 1

func _get_preset_name(preset_index: int) -> String:
	return "Default"

func _get_import_options(path: String, preset_index: int) -> Array[Dictionary]:
	return [
		{"name": "periodic", "default_value": true},
		{"name": "pattern_size", "default_value": 3, "property_hint": PROPERTY_HINT_RANGE, "hint_string": "2,4,1"},
		{"name": "pattern_symmetry", "default_value": WfcOverlappingModel2DParams.SYMMETRY_EVERY_WHICHAWAY, "property_hint": PROPERTY_HINT_ENUM, "hint_string": "Fixed:1,Mirror:3,Every Whichaway:255"},
	]

func _get_option_visibility(path: String, option_name: StringName, options: Dictionary) -> bool:
	return true

func _get_import_order() -> int:
	return 0

func _import(source_file: String, save_path: String, options: Dictionary, platform_variants: Array[String], gen_files: Array[String]) -> int:
	var image := Image.load_from_file(ProjectSettings.globalize_path(source_file))
	var params := WfcOverlappingModel2DParams.create_from_image(image, options.periodic, options.pattern_size, options.pattern_symmetry)
	var filename := save_path + "." + _get_save_extension()
	return ResourceSaver.save(params, filename)

LICENSE.txt

MIT License

Copyright (c) 2024 Apples

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.

wfc_model_2d.gd

class_name WfcModel2D
extends RefCounted
## An implementation of Wave Function Collapse for 2D grids.
## Typically used with a params object, such as [WfcOverlappingModel2DParams].
## 
## Example usage:
## [codeblock]
## # Creating a params object
## var image: Image = preload("res://wfc_input.png")
## var params := WfcOverlappingModel2DParams.create_from_image(image, true, 3, WfcOverlappingModel2DParams.SYMMETRY_EVERY_WHICHAWAY)
## 
## # Initializing the algorithm
## var wfc := WfcModel2D.new(params)
## 
## # Running the algorithm
## wfc.prepare(100, 100)
## var rng := RandomNumberGenerator.new()
## var run_result := wfc.run(rng, 8000)
## while run_result == WfcModel2D.RUN_RESULT_UNFINISHED:
##     await get_tree().process_frame
##     run_result = wfc.run(rng, 8000)
## match run_result:
##     WfcModel2D.RUN_RESULT_CONTRADICTION:
##         print("Algorithm failed. Please retry.")
##     WfcModel2D.RUN_RESULT_SUCCESS:
##         print("Success!")
## 
## # Get successful result
## assert(run_result == WfcModel2D.RUN_RESULT_SUCCESS)
## var result: Image = params.get_result(wfc)
## [/codeblock]

## List of DIRECTIONS for navigation (W, N, E, S)
const DIRECTIONS: Array[Vector2i] = [Vector2i(-1, 0), Vector2i(0, 1), Vector2i(1, 0), Vector2i(0, -1)]

## Maps a direction index to the index of the OPPOSITE direction
const OPPOSITE: Array[int] = [2, 3, 0, 1]

enum {
	## The model is completely uninitizalized. No operations are permitted other than [method prepare].
	RUN_STATE_UNINITIALIZED,
	## The model is prepared. Valid operations are [method constrain_cell] and [method run].
	RUN_STATE_PREPARED,
	## The model is currently running. Valid operations are [method run].
	RUN_STATE_RUNNING,
	## The model has completed its run. Call [method prepare] to reset.
	RUN_STATE_COMPLETE,
}

enum {
	## The run completed and was a success. The solution is stored in [member output].
	RUN_RESULT_SUCCESS,
	## The run completed and was a failute. No solution was found.
	RUN_RESULT_CONTRADICTION,
	## The run is unfinished.
	RUN_RESULT_UNFINISHED,
}

## For each possible assignment value and direction, provides a list of allowed adjacent values.
var legal_adjacencies: Array[PackedInt32Array] # [t * 4 + d][]

## The occurrence rate of each assignment value. Used to calculate entropy.
var weights: PackedFloat64Array # [t]

## The current state of the model. Must not be modified by the user.
var run_state: int = RUN_STATE_UNINITIALIZED

## The width of the output solution.
var output_size_x: int

## The height of the output solution.
var output_size_y: int

## If true, the output solution will be seamlessly tileable.
var is_output_periodic: bool

## The output produced by [method run], row major order.
var output: PackedInt32Array # [y][x] => [y * output_size_x + x]

enum {
	_OBSERVE_RESULT_PROPAGATE,
	_OBSERVE_RESULT_ALL_DONE,
	_OBSERVE_RESULT_CONTRADICTION,
}

enum {
	_PROPAGATE_RESULT_UNFINISHED,
	_PROPAGATE_RESULT_DONE,
	_PROPAGATE_RESULT_CONTRADICTION,
}

var _num_values: int # amount of t values

var _wave: PackedByteArray # [y][x][t] => [(y * output_size_x + x) * _num_values + t]
var _wave_length: int # output_size_x * output_size_y

var _num_compatible_from_dir: PackedInt32Array # [(x, y)][t][OPPOSITE[d]] => [((y * output_size_x + x) * _num_values + t) * 4 + d]

var _propagate_stack: PackedInt32Array # interleaved tuples (i, t)
var _propagate_stack_top: int

var _weight_log_weights: PackedFloat64Array # [t]

var _num_possible_assignments: PackedInt32Array # [y][x] => [y * output_size_x + x]

# These variables are used to cache and quickly compute entropy values
var _starting_entropy: float # maximum possible entropy
var _sum_of_weights: float
var _sum_of_weight_log_weights: float
var _sums_of_weights: PackedFloat64Array # [y][x] => [y * output_size_x + x]
var _sums_of_weight_log_weights: PackedFloat64Array # [y][x] => [y * output_size_x + x]

# Shannon entropy for each cell, measured in nats
var _entropies: PackedFloat64Array # [y][x] => [y * output_size_x + x]

## [param params]: initial parameters to use.
func _init(params: WfcModel2DParams = null):
	if params:
		legal_adjacencies = params.legal_adjacencies.duplicate(true)
		weights = params.weights.duplicate()

## Resets the current state. Must be called before each new run.
func prepare(p_output_size_x: int, p_output_size_y: int, p_is_output_periodic: bool = false) -> void:
	if run_state == RUN_STATE_UNINITIALIZED:
		_initialize_wfc()
	output_size_x = p_output_size_x
	output_size_y = p_output_size_y
	is_output_periodic = p_is_output_periodic
	_clear()
	run_state = RUN_STATE_PREPARED

func assign_values(x: int, y: int, values: Array[int]) -> void:
	var i := x + y * output_size_x
	for t in _num_values:
		if not t in values and _wave[i * _num_values + t]:
			_ban(i, t)

func ban_values(x: int, y: int, values: Array[int]) -> void:
	var i := x + y * output_size_x
	for t in values:
		if _wave[i * _num_values + t]:
			_ban(i, t)

func save_state() -> Dictionary:
	return {
		run_state = run_state,
		output_size_x = output_size_x,
		output_size_y = output_size_y,
		is_output_periodic = is_output_periodic,
		_num_values = _num_values,
		_weight_log_weights = _weight_log_weights,
		_sum_of_weights = _sum_of_weights,
		_sum_of_weight_log_weights = _sum_of_weight_log_weights,
		_starting_entropy = _starting_entropy,
		_wave_length = _wave_length,
		_wave = _wave.duplicate(),
		_num_compatible_from_dir = _num_compatible_from_dir.duplicate(),
		output = output.duplicate(),
		_num_possible_assignments = _num_possible_assignments.duplicate(),
		_sums_of_weights = _sums_of_weights.duplicate(),
		_sums_of_weight_log_weights = _sums_of_weight_log_weights.duplicate(),
		_entropies = _entropies.duplicate(),
		_propagate_stack = _propagate_stack.duplicate(),
		_propagate_stack_top = _propagate_stack_top,
	}

func restore_state(state: Dictionary) -> void:
	run_state = state.run_state
	output_size_x = state.output_size_x
	output_size_y = state.output_size_y
	is_output_periodic = state.is_output_periodic
	_num_values = state._num_values
	_weight_log_weights = state._weight_log_weights
	_sum_of_weights = state._sum_of_weights
	_sum_of_weight_log_weights = state._sum_of_weight_log_weights
	_starting_entropy = state._starting_entropy
	_wave_length = state._wave_length
	_wave = state._wave.duplicate()
	_num_compatible_from_dir = state._num_compatible_from_dir.duplicate()
	output = state.output.duplicate()
	_num_possible_assignments = state._num_possible_assignments.duplicate()
	_sums_of_weights = state._sums_of_weights.duplicate()
	_sums_of_weight_log_weights = state._sums_of_weight_log_weights.duplicate()
	_entropies = state._entropies.duplicate()
	_propagate_stack = state._propagate_stack.duplicate()
	_propagate_stack_top = state._propagate_stack_top

## Runs the WFC algorithm and attempts to find a solution.
##
## [param rng] is the RNG used by all random variables.
## [param time_limit] is the run time limit, measured in microseconds. A value of -1 means unlimited.
##
## Returns one of [constant RUN_RESULT_SUCCESS], [constant RUN_RESULT_CONTRADICTION],
## or [constant RUN_RESULT_UNFINISHED].
##
## If successful, the solution will be stored in [member output].
##
## If a [param time_limit] is given, will only run iterations of the algorithm while the time
## elapsed is less than [param time_limit]. If the time limit is hit before the algorithm is
## complete, will return [constant RUN_RESULT_UNFINISHED]. If [method run] is called again while
## unfinished, the samerng must be passed in again to ensure determinism.
func run(rng := RandomNumberGenerator.new(), time_limit: int = -1) -> int:
	match run_state:
		RUN_STATE_UNINITIALIZED:
			push_error("Tried to call run() for the first time without being prepared. Call prepare() first.")
			return RUN_RESULT_CONTRADICTION
		RUN_STATE_COMPLETE:
			push_error("Tried to call run() after already being complete. Call prepare() to reset for another run.")
			return RUN_RESULT_CONTRADICTION
	
	run_state = RUN_STATE_RUNNING
	
	var start_time := Time.get_ticks_usec()
	
	# Start off with a propagate to clear out any pre-applied constraints.
	# Usually this is needed when the user calls constrain_cell().
	# Also needed when propagation was previously left unfinished due to the time_limit.
	match _propagate(time_limit, start_time):
		_PROPAGATE_RESULT_UNFINISHED:
			# If _propagate was left unfinished, consider the run unfinished.
			return RUN_RESULT_UNFINISHED
		_PROPAGATE_RESULT_CONTRADICTION:
			# If a contradiction was found, end the run.
			run_state = RUN_STATE_COMPLETE
			return RUN_RESULT_CONTRADICTION
	
	# Iterate until the time_limit is reached.
	while time_limit == -1 or Time.get_ticks_usec() - start_time < time_limit:
		match _observe(rng):
			_OBSERVE_RESULT_PROPAGATE:
				# If not done yet, _observe will tell us we need to _propagate.
				match _propagate(time_limit, start_time):
					_PROPAGATE_RESULT_UNFINISHED:
						return RUN_RESULT_UNFINISHED
					_PROPAGATE_RESULT_CONTRADICTION:
						run_state = RUN_STATE_COMPLETE
						return RUN_RESULT_CONTRADICTION
			_OBSERVE_RESULT_CONTRADICTION:
				# If a contradiction was found, end the run.
				run_state = RUN_STATE_COMPLETE
				return RUN_RESULT_CONTRADICTION
			_OBSERVE_RESULT_ALL_DONE:
				# Finished successfully. No _propagate necessary.
				# Fill out the output result.
				for i: int in range(_wave_length):
					for t: int in range(_num_values):
						if _wave[i * _num_values + t]:
							output[i] = t
							break
				run_state = RUN_STATE_COMPLETE
				return RUN_RESULT_SUCCESS
	
	# Unfinished due to reaching the time_limit.
	return RUN_RESULT_UNFINISHED

# Initializes the constant state of the model.
func _initialize_wfc() -> void:
	assert(run_state == RUN_STATE_UNINITIALIZED)
	assert(weights.size() == legal_adjacencies.size() / DIRECTIONS.size())
	_num_values = weights.size()
	
	_weight_log_weights.resize(_num_values)
	_sum_of_weights = 0
	_sum_of_weight_log_weights = 0
	
	for t: int in range(0, _num_values):
		_weight_log_weights[t] = weights[t] * log(weights[t])
		_sum_of_weights += weights[t]
		_sum_of_weight_log_weights += _weight_log_weights[t]
	
	_starting_entropy = log(_sum_of_weights) - _sum_of_weight_log_weights / _sum_of_weights

# Clears and resets the per-run state of the model.
func _clear() -> void:
	_wave_length = output_size_x * output_size_y
	_wave.resize(_wave_length * _num_values)
	_num_compatible_from_dir.resize(_wave_length * _num_values * 4)
	output.resize(_wave_length)
	
	_num_possible_assignments.resize(_wave_length)
	_sums_of_weights.resize(_wave_length)
	_sums_of_weight_log_weights.resize(_wave_length)
	_entropies.resize(_wave_length)
	
	_propagate_stack.resize(_wave_length * _num_values * 2)
	_propagate_stack_top = 0
	
	for i: int in range(_wave_length):
		for t: int in range(_num_values):
			_wave[i * _num_values + t] = 1
			for d: int in range(4):
				var p_i: int = t * 4 + OPPOSITE[d]
				_num_compatible_from_dir[(i * _num_values + t) * 4 + d] = legal_adjacencies[p_i].size()
		
		_num_possible_assignments[i] = weights.size()
		_sums_of_weights[i] = _sum_of_weights
		_sums_of_weight_log_weights[i] = _sum_of_weight_log_weights
		_entropies[i] = _starting_entropy
		output[i] = -1

# Attempts to assign a value to the cell with the lowest entropy.
func _observe(rng: RandomNumberGenerator) -> int:
	var lowest_seen_entropy: float = 1.0e+4 # Initialized to a Big Number (tm).
	var unobserved_cell: int = -1
	
	# Find the next unobserved index based on entropy.
	for i: int in range(_wave_length):
		@warning_ignore("integer_division")
		var v := Vector2i(i % output_size_x, i / output_size_x)
		
		# If any cell has no remaining possible assignments, we've created a contradiction.
		var num_possible_assignments := _num_possible_assignments[i]
		if num_possible_assignments == 0:
			return _OBSERVE_RESULT_CONTRADICTION
		
		# Entropy check.
		var entropy: float = _entropies[i]
		if num_possible_assignments > 1 and entropy <= lowest_seen_entropy:
			# Add a small random value to the entropy just to shake things up.
			entropy += 1.0e-6 * rng.randf()
			if entropy < lowest_seen_entropy:
				lowest_seen_entropy = entropy
				unobserved_cell = i
	
	# If there are no unobserved indices left, we're completely done.
	if unobserved_cell == -1:
		return _OBSERVE_RESULT_ALL_DONE
	
	# Choose a random _num_values from the cell's possible assignments.
	
	var possible_assignments := PackedInt32Array() # [i] => t
	var possible_assignment_weights := PackedFloat64Array() # [i] => weights[t]
	var sum_of_possible_assignment_weights: float = 0.0
	
	for t: int in range(_num_values):
		if _wave[unobserved_cell * _num_values + t]:
			var w := weights[t]
			possible_assignments.append(t)
			possible_assignment_weights.append(w)
			sum_of_possible_assignment_weights += w
	
	var weight_threshold := rng.randf_range(0.0, sum_of_possible_assignment_weights)
	
	var chosen_assignment: int = 0
	
	var partial_sum: float = 0;
	for i: int in range(possible_assignment_weights.size()):
		partial_sum += possible_assignment_weights[i]
		if partial_sum >= weight_threshold:
			chosen_assignment = possible_assignments[i]
			break
	
	# Ban all assignments for this cell except for the chosen one.
	for t: int in range(_num_values):
		if (_wave[unobserved_cell * _num_values + t] != int(t == chosen_assignment)):
			_ban(unobserved_cell, t)
	
	return _OBSERVE_RESULT_PROPAGATE

# Propagates wave state changes from _ban.
func _propagate(time_limit: int, start_time: int) -> int:
	while _propagate_stack_top > 0 and (time_limit == -1 or Time.get_ticks_usec() - start_time < time_limit):
		var i1: int = _propagate_stack[_propagate_stack_top - 2]
		var t1: int = _propagate_stack[_propagate_stack_top - 1]
		_propagate_stack_top -= 2
		
		@warning_ignore("integer_division")
		var v1 := Vector2i(i1 % output_size_x, i1 / output_size_x)
		
		for d: int in range(4):
			var v2 := v1 + DIRECTIONS[d]
			
			if not is_output_periodic:
				if not Rect2i(0, 0, output_size_x, output_size_y).has_point(v2):
					continue
			else:
				# If the output is periodic, we need to propagate to every cell and also wrap around
				# the edges of the output.
				v2.x = posmod(v2.x, output_size_x)
				v2.y = posmod(v2.y, output_size_y)
			
			var i2: int = v2.x + v2.y * output_size_x
			var p_i1: int = t1 * 4 + d
			
			for t2: int in legal_adjacencies[p_i1]:
				var comp_i := (i2 * _num_values + t2) * 4 + d
				_num_compatible_from_dir[comp_i] -= 1
				if _num_compatible_from_dir[comp_i] == 0:
					_ban(i2, t2)
					if _num_possible_assignments[i2] == 0:
						return _PROPAGATE_RESULT_CONTRADICTION
	
	if _propagate_stack_top > 0:
		return _PROPAGATE_RESULT_UNFINISHED
	
	return _PROPAGATE_RESULT_DONE

# Removes the specified assignment from a cell's possibilities. Queues up a change for propagation.
func _ban(i: int, t: int) -> void:
	var wave_index := i * _num_values + t
	assert(_wave[wave_index])
	_wave[wave_index] = 0
	
	for d: int in range(4):
		_num_compatible_from_dir[wave_index * 4 + d] = 0
	
	# Push onto the propagation _propagate_stack.
	_propagate_stack[_propagate_stack_top] = i
	_propagate_stack[_propagate_stack_top + 1] = t
	_propagate_stack_top += 2
	
	_num_possible_assignments[i] -= 1
	_sums_of_weights[i] -= weights[t]
	_sums_of_weight_log_weights[i] -= _weight_log_weights[t]
	
	var sum: float = _sums_of_weights[i]
	_entropies[i] = log(sum) - _sums_of_weight_log_weights[i] / sum

wfc_model_2d_params.gd

class_name WfcModel2DParams
extends Resource
## Provides storage for initial parameters to a [WfcModel2D].
## It is expected that this class will be extended to implement various forms
## of Wave Function Collapse, e.g. [WfcOverlappingModel2DParams]

## For each possible assignment value and direction, provides a list of allowed adjacent values.
## Indexed as t * 4 + d, where t is the assignment value, and d is the direction index.
@export var legal_adjacencies: Array[PackedInt32Array] # [t * 4 + d][]

## The occurrence rate of each assignment value. Used to calculate entropy.
@export var weights: PackedFloat64Array # [t]

wfc_overlapping_model_2d_params.gd

class_name WfcOverlappingModel2DParams
extends WfcModel2DParams
## Implements the overlapped pattern model of WFC.
## It is expected to be used as part of the import process for an input image,
## but can also be created at runtime for dynamic inputs.
##
## The WFC model expects abstract values labelled 0..N, where N is the total number of values.
## Here, we analyze an input inage for patterns, and assign each pattern an index value.
## These index values are what WFC will assign to cells, and we can then convert each cell
## to a color by examining its pattern after WFC has assigned a value to all cells.
##
## Normally, for non-tiling outputs, we could reduce the output size of the WFC model
## by pattern_size - 1, and then infer the color values of the resulting border by examining
## nearby patterns. This would be a significant optimization.
## To implement this, the WFC model's run() function would need to be modified to only operate on a
## specific region of the output, and this would need to be specified by the params object.
## Propagation still needs to occur on the entire output before run() is called, but while run()
## is running, propagation only needs to work in the specified output region.

enum {
	## The original orientation of the pattern from the input.
	SYMMETRY_INPUT = 1 << 0,
	## The pattern rotated 90 degrees counterclockwise.
	SYMMETRY_ROT_90 = 1 << 1,
	## The pattern rotated 180 degrees counterclockwise.
	SYMMETRY_ROT_180 = 1 << 2,
	## The pattern rotated 270 degrees counterclockwise.
	SYMMETRY_ROT_270 = 1 << 3,
	## The pattern mirrored on the Y axis.
	SYMMETRY_HFLIP = 1 << 4,
	## The pattern mirrored on the Y axis, and then rotated 90 degrees counterclockwise.
	SYMMETRY_HFLIP_ROT_90 = 1 << 5,
	## The pattern mirrored on the Y axis, and then rotated 180 degrees counterclockwise.
	SYMMETRY_HFLIP_ROT_180 = 1 << 6,
	## The pattern mirrored on the Y axis, and then rotated 270 degrees counterclockwise.
	SYMMETRY_HFLIP_ROT_270 = 1 << 7,
	## [constant SYMMETRY_INPUT] | [constant SYMMETRY_HFLIP].
	SYMMETRY_MIRROR = SYMMETRY_INPUT | SYMMETRY_HFLIP,
	## All symmetries.
	SYMMETRY_EVERY_WHICHAWAY = 0b11111111,
}

## The input image from which the patterns were gathered.
@export var input: Image

## If true, the input image will be tiled.
@export var is_input_periodic: bool

## Width and height of the patterns.
@export var pattern_size: int

## A set of bits, each corresponding to a symmetry which is enabled if the bit is set.
## Each enabled symmetry will be used to rotate or mirror each pattern from the input.
@export_flags("Input", "Rotate CCW 90", "Rotate CCW 180", "Rotate CCW 270", "Flip Horizontal", "Flip and Rotate 90", "Flip and Rotate 180", "Flip and Rotate 270")
var pattern_symmetry: int = SYMMETRY_EVERY_WHICHAWAY

## A list of all the patterns.
## A pattern is a [member pattern_size] x [member pattern_size] array,
## where each element is an index of the [member colors] array.
## Each pattern's index is the value used in the WFC model cell values.
@export var patterns: Array[PackedByteArray]

## A list of the colors from the input image.
## Each color's index is the value used in the pattern arrays.
@export var colors: PackedColorArray

## Creates a parameter set from an input image using the overlapping technique.
static func create_from_image(input: Image, is_input_periodic: bool, pattern_size: int,  pattern_symmetry: int) -> WfcOverlappingModel2DParams:
	var result := WfcOverlappingModel2DParams.new()
	result.input = input
	result.is_input_periodic = is_input_periodic
	result.pattern_size = pattern_size
	result.pattern_symmetry = pattern_symmetry
	
	var input_width := input.get_width()
	var input_height := input.get_height()
	
	# First we need to convert the input image into an indexed byte array.
	# Each color will be added to result.colors, and its index in that array
	# will be used as the value in the sample array.
	
	var sample := PackedByteArray() # [x + y * input_width]
	sample.resize(input_width * input_height)
	
	for i: int in range(sample.size()):
		var color := input.get_pixel(i % input_width, i / input_width)
		var k := result.colors.find(color)
		if k == -1:
			assert(result.colors.size() < 256)
			result.colors.append(color)
			k = result.colors.size() - 1
		sample[i] = k
	
	# Maps patterns to their indices.
	# Mostly used for de-duplication.
	var pattern_indices: Dictionary = {} # { [pattern: PackedByteArray]: int }
	
	var xmax: int = input_width if is_input_periodic else input_width - pattern_size + 1
	var ymax: int = input_height if is_input_periodic else input_height - pattern_size + 1
	
	# Loop over each pixel and make patterns.
	# A pattern is formed by using that pixel position as the upper-left corner
	# of a pattern_size x pattern_size square.
	
	for y in range(0, ymax):
		for x in range(0, xmax):
			var current_patterns: Array[PackedByteArray] = []
			
			var make_pattern := func make_pattern(dx: int, dy: int) -> int:
				return sample[(x + dx) % input_width + (y + dy) % input_height * input_width]
			
			# Each of these entries represents a symmetry,
			# in the same order as the SYMMETRY_ enum values.
			current_patterns.append(_pattern(make_pattern, pattern_size)) # SYMMETRY_INPUT
			current_patterns.append(_rotate(current_patterns[0], pattern_size)) # SYMMETRY_ROT_90
			current_patterns.append(_rotate(current_patterns[1], pattern_size)) # SYMMETRY_ROT_180
			current_patterns.append(_rotate(current_patterns[2], pattern_size)) # SYMMETRY_ROT_270
			current_patterns.append(_reflect(current_patterns[0], pattern_size)) # SYMMETRY_HFLIP
			current_patterns.append(_reflect(current_patterns[3], pattern_size)) # SYMMETRY_HFLIP_ROT_90
			current_patterns.append(_reflect(current_patterns[2], pattern_size)) # SYMMETRY_HFLIP_ROT_180
			current_patterns.append(_reflect(current_patterns[1], pattern_size)) # SYMMETRY_HFLIP_ROT_270
			
			for k in range(8):
				# Skip symmetries not flagged in pattern_symmetry
				if not (pattern_symmetry & (1 << k)):
					continue
				
				var p := current_patterns[k]
				if pattern_indices.has(p):
					# If a pattern has been seen already, we need to increment its weight.
					# This is the only situation where weights can be anything other than 1.
					var index: int = pattern_indices[p]
					result.weights[index] = result.weights[index] + 1
				else:
					result.patterns.append(p)
					result.weights.append(1.0)
					pattern_indices[p] = result.patterns.size() - 1
	
	# After all patterns have been discovered, we must construct their legal_adjacencies.
	# Patterns are legally adjacent if their overlapping regions contain the same values.
	# Overlapping regions are produced by shifting one of the patterns by a single pixel
	# in any cardinal direction.
	
	result.legal_adjacencies.resize(result.patterns.size() * 4)
	
	# t is the index of the pattern which we are calculating adjacencies for.
	for t: int in range(result.patterns.size()):
		for d: int in range(4):
			var p := result.legal_adjacencies[t * 4 + d]
			# t2 is the index of the pattern which might be adjacent to t.
			for t2: int in range(result.patterns.size()):
				if _agrees(result.patterns[t], result.patterns[t2], WfcModel2D.DIRECTIONS[d], pattern_size):
					p.append(t2)
	
	return result

## Gets the indices of all patterns which have the given [param color] as their primary color.
## The primary color of a pattern is the color of its upper-left pixel.
## [param color]: the color.
func get_patterns_for_color(color: Color) -> Array[int]:
	var color_value := colors.find(color)
	if color_value == -1:
		return []
	
	var result: Array[int] = []
	for t: int in range(patterns.size()):
		if patterns[t][0] == color_value:
			result.append(t)
	return result

## Gets the indices of all patterns which contain the given [param color] anywhere in them.
## [param color]: the color.
func get_patterns_containing_color(color: Color) -> Array[int]:
	var color_value := colors.find(color)
	assert(color_value != -1)
	if color_value == -1:
		return []
	
	var result: Array[int] = []
	for t: int in range(patterns.size()):
		if patterns[t].has(color_value):
			result.append(t)
	return result

## Converts the output of a model to an image.
## The [param model] must have completed a run and be in [constant WfcModel2D.RUN_STATE_COMPLETE].
## If the run was a success, this returns an image the same size as the model output,
## with the primary colors from the patterns which were assigned to each cell.
## If the run was unsuccessful, this returns null.
##
## [param model]: the model. Must have been [method WfcModel2D.run].
func get_result(model: WfcModel2D) -> Image:
	assert(model.run_state == WfcModel2D.RUN_STATE_COMPLETE)
	if model.run_state != WfcModel2D.RUN_STATE_COMPLETE:
		return null
	if model.output[0] < 0:
		return null
	var result := Image.create(model.output_size_x, model.output_size_y, false, Image.FORMAT_RGBA8)
	for y in range(model.output_size_y):
		var dy: int = 0 if y < model.output_size_y - pattern_size + 1 else pattern_size - 1
		for x in range(model.output_size_x):
			var dx: int = 0 if x < model.output_size_x - pattern_size + 1 else pattern_size - 1
			result.set_pixel(x, y, colors[patterns[model.output[x - dx + (y - dy) * model.output_size_x]][dx + dy * pattern_size]])
	return result

# Constructs a pattern using a mapping function.
# [param f]: called for each position in the pattern, func f(x: int, y: int) -> int
static func _pattern(f: Callable, pattern_size: int) -> PackedByteArray:
	var result := PackedByteArray()
	result.resize(pattern_size * pattern_size)
	for y in range(0, pattern_size):
		for x in range(0, pattern_size):
			result[x + y * pattern_size] = f.call(x, y)
	return result

# Constructs a copy of the givne pattern, rotated 90 degrees counterclockwise.
static func _rotate(p: PackedByteArray, pattern_size: int) -> PackedByteArray:
	return _pattern(func (x: int, y: int): return p[pattern_size - 1 - y + x * pattern_size], pattern_size)

# Constructs a copy of the givne pattern, mirrored across the Y axis.
static func _reflect(p: PackedByteArray, pattern_size: int) -> PackedByteArray:
	return _pattern(func (x: int, y: int): return p[pattern_size - 1 - x + y * pattern_size], pattern_size)

# Determines if the overlap of the two patterns has matching pixels.
static func _agrees(p1: PackedByteArray, p2: PackedByteArray, dv: Vector2i, pattern_size: int) -> bool:
	var xmin: int = 0 if dv.x < 0 else dv.x
	var xmax: int = (dv.x + pattern_size) if dv.x < 0 else pattern_size
	var ymin: int = 0 if dv.y < 0 else dv.y
	var ymax: int = (dv.y + pattern_size) if dv.y < 0 else pattern_size
	for y in range(ymin, ymax):
		for x in range(xmin, xmax):
			if p1[x + pattern_size * y] != p2[x - dv.x + pattern_size * (y - dv.y)]:
				return false
	return true