void compute_ideal_weights_for_decimation(
	const endpoints_and_weights& eai_in,
	endpoints_and_weights& eai_out,
	const decimation_info& di,
	float* dec_weight_ideal_value,
	float* dec_weight_ideal_sig
) {
	unsigned int texel_count = di.texel_count;
	unsigned int weight_count = di.weight_count;

	promise(texel_count > 0);
	promise(weight_count > 0);

	// This function includes a copy of the epw from eai_in to eai_out. We do it here because we
	// want to load the data anyway, so we can avoid loading it from memory twice.
	eai_out.ep = eai_in.ep;
	eai_out.is_constant_weight_error_scale = eai_in.is_constant_weight_error_scale;

	// Ensure that the end of the output arrays that are used for SIMD paths later are filled so we
	// can safely run SIMD elsewhere without a loop tail. Note that this is always safe as weight
	// arrays always contain space for 64 elements
	unsigned int weight_count_simd = round_up_to_simd_multiple_vla(weight_count);
	for (unsigned int i = weight_count; i < weight_count_simd; i++)
	{
		dec_weight_ideal_value[i] = 0.0f;
	}

	// If we have a 1:1 mapping just shortcut the computation - clone the weights into both the
	// weight set and the output epw copy.

	// Transfer enough to also copy zero initialized SIMD over-fetch region
	unsigned int texel_count_simd = round_up_to_simd_multiple_vla(texel_count);
	if (texel_count == weight_count)
	{
		for (unsigned int i = 0; i < texel_count_simd; i++)
		{
			// Assert it's an identity map for valid texels, and last valid value for any overspill
			assert(((i < texel_count) && (i == di.weight_texel[0][i])) ||
			       ((i >= texel_count) && (texel_count - 1 == di.weight_texel[0][i])));
			dec_weight_ideal_value[i] = eai_in.weights[i];
			dec_weight_ideal_sig[i] = eai_in.weight_error_scale[i];

			eai_out.weights[i] = eai_in.weights[i];
			eai_out.weight_error_scale[i] = eai_in.weight_error_scale[i];
		}

		return;
	}
	// If we don't have a 1:1 mapping just clone the weights into the output epw copy and then do
	// the full algorithm to decimate weights.
	else
	{
		for (unsigned int i = 0; i < texel_count_simd; i++)
		{
			eai_out.weights[i] = eai_in.weights[i];
			eai_out.weight_error_scale[i] = eai_in.weight_error_scale[i];
		}
	}

	// Otherwise compute an estimate and perform single refinement iteration
	alignas(ASTCENC_VECALIGN) float infilled_weights[BLOCK_MAX_TEXELS];

	// Compute an initial average for each decimated weight
	bool constant_wes = eai_in.is_constant_weight_error_scale;
	vfloat weight_error_scale(eai_in.weight_error_scale[0]);

	// This overshoots - this is OK as we initialize the array tails in the
	// decimation table structures to safe values ...
	for (unsigned int i = 0; i < weight_count; i += ASTCENC_SIMD_WIDTH)
	{
		// Start with a small value to avoid div-by-zero later
		vfloat weight_weight(1e-10f);
		vfloat initial_weight = vfloat::zero();

		// Accumulate error weighting of all the texels using this weight
		vint weight_texel_count(di.weight_texel_count + i);
		unsigned int max_texel_count = hmax(weight_texel_count).lane<0>();
		promise(max_texel_count > 0);

		for (unsigned int j = 0; j < max_texel_count; j++)
		{
			// Not all lanes may actually use j texels, so mask out if idle
			vmask active = weight_texel_count > vint(j);

			vint texel(di.weight_texel[j] + i);
			texel = select(vint::zero(), texel, active);

			vfloat weight = loada(di.weights_flt[j] + i);
			weight = select(vfloat::zero(), weight, active);

			if (!constant_wes)
			{
				weight_error_scale = gatherf(eai_in.weight_error_scale, texel);
			}

			vfloat contrib_weight = weight * weight_error_scale;

			weight_weight += contrib_weight;
			initial_weight += gatherf(eai_in.weights, texel) * contrib_weight;
		}

		storea(weight_weight, dec_weight_ideal_sig + i);
		storea(initial_weight / weight_weight, dec_weight_ideal_value + i);
	}

	// Populate the interpolated weight grid based on the initital average
	// Process SIMD-width texel coordinates at at time while we can. Safe to
	// over-process full SIMD vectors - the tail is zeroed.
	for (unsigned int i = 0; i < texel_count; i += ASTCENC_SIMD_WIDTH)
	{
		vfloat weight = bilinear_infill_vla(di, dec_weight_ideal_value, i);
		storea(weight, infilled_weights + i);
	}

	// Perform a single iteration of refinement
	// Empirically determined step size; larger values don't help but smaller drops image quality
	constexpr float stepsize = 0.25f;
	constexpr float chd_scale = -WEIGHTS_TEXEL_SUM;

	for (unsigned int i = 0; i < weight_count; i += ASTCENC_SIMD_WIDTH)
	{
		vfloat weight_val = loada(dec_weight_ideal_value + i);

		// Accumulate error weighting of all the texels using this weight
		// Start with a small value to avoid div-by-zero later
		vfloat error_change0(1e-10f);
		vfloat error_change1(0.0f);

		// Accumulate error weighting of all the texels using this weight
		vint weight_texel_count(di.weight_texel_count + i);
		unsigned int max_texel_count = hmax(weight_texel_count).lane<0>();
		promise(max_texel_count > 0);

		for (unsigned int j = 0; j < max_texel_count; j++)
		{
			// Not all lanes may actually use j texels, so mask out if idle
			vmask active = weight_texel_count > vint(j);

			vint texel(di.weight_texel[j] + i);
			texel = select(vint::zero(), texel, active);

			vfloat contrib_weight = loada(di.weights_flt[j] + i);
			contrib_weight = select(vfloat::zero(), contrib_weight, active);

			if (!constant_wes)
			{
 				weight_error_scale = gatherf(eai_in.weight_error_scale, texel);
			}

			vfloat scale = weight_error_scale * contrib_weight;
			vfloat old_weight = gatherf(infilled_weights, texel);
			vfloat ideal_weight = gatherf(eai_in.weights, texel);

			error_change0 += contrib_weight * scale;
			error_change1 += (old_weight - ideal_weight) * scale;
		}


		vfloat step = (error_change1 * chd_scale) / error_change0;
		step = clamp(-stepsize, stepsize, step);

		// Update the weight; note this can store negative values.
		storea(weight_val + step, dec_weight_ideal_value + i);
	}
}