static float prepare_error_weight_block(
	const astcenc_context& ctx,
	const astcenc_image& image,
	const block_size_descriptor& bsd,
	const image_block& blk,
	error_weight_block& ewb
) {
	unsigned int idx = 0;
	bool any_mean_stdev_weight =
	    ctx.config.v_rgb_mean != 0.0f || ctx.config.v_rgb_stdev != 0.0f || \
	    ctx.config.v_a_mean != 0.0f || ctx.config.v_a_stdev != 0.0f;

	vfloat4 color_weights(ctx.config.cw_r_weight,
	                      ctx.config.cw_g_weight,
	                      ctx.config.cw_b_weight,
	                      ctx.config.cw_a_weight);

	// This works because HDR is imposed globally at compression time
	unsigned int rgb_lns = blk.rgb_lns[0];
	unsigned int a_lns = blk.alpha_lns[0];
	vint4 use_lns(rgb_lns, rgb_lns, rgb_lns, a_lns);
	vmask4 lns_mask = use_lns != vint4::zero();

	promise(bsd.xdim > 0);
	promise(bsd.ydim > 0);
	promise(bsd.zdim > 0);

	for (unsigned int z = 0; z < bsd.zdim; z++)
	{
		for (unsigned int y = 0; y < bsd.ydim; y++)
		{
			for (unsigned int x = 0; x < bsd.xdim; x++)
			{
				unsigned int xpos = x + blk.xpos;
				unsigned int ypos = y + blk.ypos;
				unsigned int zpos = z + blk.zpos;

				if (xpos >= image.dim_x || ypos >= image.dim_y || zpos >= image.dim_z)
				{
					ewb.error_weights[idx] = vfloat4(1e-11f);
				}
				else
				{
					vfloat4 derv(65535.0f);

					// Compute derivative if we have any use of LNS
					if (any(lns_mask))
					{
						vfloat4 data = blk.texel(idx);
						vint4 datai = lns_to_sf16(float_to_int(data));

						vfloat4 dataf = float16_to_float(datai);
						dataf = max(dataf, 6e-5f);

						vfloat4 data_lns1 = dataf * 1.05f;
						data_lns1 = float_to_lns(data_lns1);

						vfloat4 data_lns2 = dataf;
						data_lns2 = float_to_lns(data_lns2);

						vfloat4 divisor_lns = dataf * 0.05f;

						// Clamp derivatives between 1/32 and 2^25
						float lo = 1.0f / 32.0f;
						float hi = 33554432.0f;
						vfloat4 derv_lns = clamp(lo, hi, (data_lns1 - data_lns2) / divisor_lns);
						derv = select(derv, derv_lns, lns_mask);
					}

					// Compute error weight
					vfloat4 error_weight(ctx.config.v_rgb_base,
					                     ctx.config.v_rgb_base,
					                     ctx.config.v_rgb_base,
					                     ctx.config.v_a_base);

					unsigned int ydt = image.dim_x;
					unsigned int zdt = image.dim_x * image.dim_y;

					if (any_mean_stdev_weight)
					{
						vfloat4 avg = ctx.input_averages[zpos * zdt + ypos * ydt + xpos];
						avg = max(avg, 6e-5f);
						avg = avg * avg;

						vfloat4 variance = ctx.input_variances[zpos * zdt + ypos * ydt + xpos];
						variance = variance * variance;

						float favg = hadd_rgb_s(avg) * (1.0f / 3.0f);
						float fvar = hadd_rgb_s(variance) * (1.0f / 3.0f);

						float mixing = ctx.config.v_rgba_mean_stdev_mix;
						avg.set_lane<0>(favg * mixing + avg.lane<0>() * (1.0f - mixing));
						avg.set_lane<1>(favg * mixing + avg.lane<1>() * (1.0f - mixing));
						avg.set_lane<2>(favg * mixing + avg.lane<2>() * (1.0f - mixing));

						variance.set_lane<0>(fvar * mixing + variance.lane<0>() * (1.0f - mixing));
						variance.set_lane<1>(fvar * mixing + variance.lane<1>() * (1.0f - mixing));
						variance.set_lane<2>(fvar * mixing + variance.lane<2>() * (1.0f - mixing));

						vfloat4 stdev = sqrt(max(variance, 0.0f));

						vfloat4 scalea(ctx.config.v_rgb_mean, ctx.config.v_rgb_mean, ctx.config.v_rgb_mean, ctx.config.v_a_mean);
						avg = avg * scalea;

						vfloat4 scales(ctx.config.v_rgb_stdev, ctx.config.v_rgb_stdev, ctx.config.v_rgb_stdev, ctx.config.v_a_stdev);
						stdev = stdev * scales;

						error_weight = error_weight + avg + stdev;
						error_weight = 1.0f / error_weight;
					}

					if (ctx.config.flags & ASTCENC_FLG_MAP_NORMAL)
					{
						// Convert from 0 to 1 to -1 to +1 range.
						float xN = ((blk.data_r[idx] * (1.0f / 65535.0f)) - 0.5f) * 2.0f;
						float yN = ((blk.data_a[idx] * (1.0f / 65535.0f)) - 0.5f) * 2.0f;

						float denom = 1.0f - xN * xN - yN * yN;
						denom = astc::max(denom, 0.1f);
						denom = 1.0f / denom;
						error_weight.set_lane<0>(error_weight.lane<0>() * (1.0f + xN * xN * denom));
						error_weight.set_lane<3>(error_weight.lane<3>() * (1.0f + yN * yN * denom));
					}

					if (ctx.config.flags & ASTCENC_FLG_USE_ALPHA_WEIGHT)
					{
						float alpha_scale;
						if (ctx.config.a_scale_radius != 0)
						{
							alpha_scale = ctx.input_alpha_averages[zpos * zdt + ypos * ydt + xpos];
						}
						else
						{
							alpha_scale = blk.data_a[idx] * (1.0f / 65535.0f);
						}

						alpha_scale = astc::max(alpha_scale, 0.0001f);

						alpha_scale *= alpha_scale;
						error_weight.set_lane<0>(error_weight.lane<0>() * alpha_scale);
						error_weight.set_lane<1>(error_weight.lane<1>() * alpha_scale);
						error_weight.set_lane<2>(error_weight.lane<2>() * alpha_scale);
					}

					error_weight = error_weight * color_weights;
					error_weight = error_weight * ctx.deblock_weights[idx];

					// When we loaded the block to begin with, we applied a transfer function and
					// computed the derivative of the transfer function. However, the error-weight
					// computation so far is based on the original color values, not the
					// transfer-function values. As such, we must multiply the error weights by the
					// derivative of the inverse of the transfer function, which is equivalent to
					// dividing by the derivative of the transfer function.

					error_weight = error_weight / (derv * derv * 1e-10f);
					ewb.error_weights[idx] = error_weight;
				}
				idx++;
			}
		}
	}

	// Small bias to avoid divide by zeros and NaN propagation later
	vfloat4 texel_weight_sum(1e-17f);
	vfloat4 error_weight_sum(1e-17f);

	int texels_per_block = bsd.texel_count;
	for (int i = 0; i < texels_per_block; i++)
	{
		texel_weight_sum += ewb.error_weights[i] * blk.texel(i);
		error_weight_sum += ewb.error_weights[i];

		float wr = ewb.error_weights[i].lane<0>();
		float wg = ewb.error_weights[i].lane<1>();
		float wb = ewb.error_weights[i].lane<2>();
		float wa = ewb.error_weights[i].lane<3>();

		ewb.texel_weight_r[i] = wr;
		ewb.texel_weight_g[i] = wg;
		ewb.texel_weight_b[i] = wb;
		ewb.texel_weight_a[i] = wa;

		ewb.texel_weight_rg[i] = (wr + wg) * 0.5f;
		ewb.texel_weight_rb[i] = (wr + wb) * 0.5f;
		ewb.texel_weight_gb[i] = (wg + wb) * 0.5f;

		ewb.texel_weight_gba[i] = (wg + wb + wa) * 0.333333f;
		ewb.texel_weight_rba[i] = (wr + wb + wa) * 0.333333f;
		ewb.texel_weight_rga[i] = (wr + wg + wa) * 0.333333f;
		ewb.texel_weight_rgb[i] = (wr + wg + wb) * 0.333333f;

		ewb.texel_weight[i] = (wr + wg + wb + wa) * 0.25f;
	}

	ewb.block_error_weighted_rgba_sum = texel_weight_sum;
	ewb.block_error_weight_sum = error_weight_sum;

	return hadd_s(error_weight_sum);
}