/src/astc-encoder/Source/astcenc_compute_variance.cpp

Source
// SPDX-License-Identifier: Apache-2.0
// ----------------------------------------------------------------------------
// Copyright 2011-2026 Arm Limited
//
// Licensed under the Apache License, Version 2.0 (the "License"); you may not
// use this file except in compliance with the License. You may obtain a copy
// of the License at:
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
// License for the specific language governing permissions and limitations
// under the License.
// ----------------------------------------------------------------------------

#if !defined(ASTCENC_DECOMPRESS_ONLY)

/**
 * @brief Functions to calculate variance per component in a NxN footprint.
 *
 * We need N to be parametric, so the routine below uses summed area tables in order to execute in
 * O(1) time independent of how big N is.
 *
 * The addition uses a Brent-Kung-based parallel prefix adder. This uses the prefix tree to first
 * perform a binary reduction, and then distributes the results. This method means that there is no
 * serial dependency between a given element and the next one, and also significantly improves
 * numerical stability allowing us to use floats rather than doubles.
 */

#include "astcenc_internal.h"

#include <cassert>

/**
 * @brief Generate a prefix-sum array using the Brent-Kung algorithm.
 *
 * This will take an input array of the form:
 *     v0, v1, v2, ...
 * ... and modify in-place to turn it into a prefix-sum array of the form:
 *     v0, v0+v1, v0+v1+v2, ...
 *
 * @param d      The array to prefix-sum.
 * @param items  The number of items in the array.
 * @param stride The item spacing in the array; i.e. dense arrays should use 1.
 */
static void brent_kung_prefix_sum(
  vfloat4* d,
  size_t items,
  size_t stride
) {
  if (items < 2)
    return;

  size_t lc_stride = 2;
  size_t log2_stride = 1;

  // The reduction-tree loop
  do {
    size_t step = lc_stride >> 1;
    size_t start = lc_stride - 1;
    size_t iters = items >> log2_stride;

    vfloat4 *da = d + (start * stride);
    ptrdiff_t ofs = -static_cast<ptrdiff_t>(step * stride);
    size_t ofs_stride = stride << log2_stride;

    while (iters)
    {
      *da = *da + da[ofs];
      da += ofs_stride;
      iters--;
    }

    log2_stride += 1;
    lc_stride <<= 1;
  } while (lc_stride <= items);

  // The expansion-tree loop
  do {
    log2_stride -= 1;
    lc_stride >>= 1;

    size_t step = lc_stride >> 1;
    size_t start = step + lc_stride - 1;
    size_t iters = (items - step) >> log2_stride;

    vfloat4 *da = d + (start * stride);
    ptrdiff_t ofs = -static_cast<ptrdiff_t>(step * stride);
    size_t ofs_stride = stride << log2_stride;

    while (iters)
    {
      *da = *da + da[ofs];
      da += ofs_stride;
      iters--;
    }
  } while (lc_stride > 2);
}

/* See header for documentation. */
void compute_pixel_region_variance(
  astcenc_contexti& ctx,
  const pixel_region_args& arg
) {
  // Unpack the memory structure into local variables
  const astcenc_image* img = arg.img;
  astcenc_swizzle swz = arg.swz;
  bool have_z = arg.have_z;

  size_t size_x = arg.size_x;
  size_t size_y = arg.size_y;
  size_t size_z = arg.size_z;

  size_t offset_x = arg.offset_x;
  size_t offset_y = arg.offset_y;
  size_t offset_z = arg.offset_z;

  size_t alpha_kernel_radius = arg.alpha_kernel_radius;

  float*   input_alpha_averages = ctx.input_alpha_averages;
  vfloat4* work_memory = arg.work_memory;

  // Compute memory sizes and dimensions that we need
  size_t kernel_radius = alpha_kernel_radius;
  size_t kerneldim = 2 * kernel_radius + 1;
  size_t kernel_radius_xy = kernel_radius;
  size_t kernel_radius_z = have_z ? kernel_radius : 0;

  size_t padsize_x = size_x + kerneldim;
  size_t padsize_y = size_y + kerneldim;
  size_t padsize_z = size_z + (have_z ? kerneldim : 0);
  size_t sizeprod = padsize_x * padsize_y * padsize_z;

  int zd_start = have_z ? 1 : 0;

  vfloat4 *varbuf1 = work_memory;
  vfloat4 *varbuf2 = work_memory + sizeprod;

  // Scaling factors to apply to Y and Z for accesses into the work buffers
  size_t yst = padsize_x;
  size_t zst = padsize_x * padsize_y;

  // Scaling factors to apply to Y and Z for accesses into result buffers
  size_t ydt = img->dim_x;
  size_t zdt = img->dim_x * img->dim_y;

  // Macros to act as accessor functions for the work-memory
  #define VARBUF1(z, y, x) varbuf1[z * zst + y * yst + x]
  #define VARBUF2(z, y, x) varbuf2[z * zst + y * yst + x]

  // Load N and N^2 values into the work buffers
  if (img->data_type == ASTCENC_TYPE_U8)
  {
    // Swizzle data structure 4 = ZERO, 5 = ONE
    uint8_t data[6];
    data[ASTCENC_SWZ_0] = 0;
    data[ASTCENC_SWZ_1] = 255;

    for (size_t z = zd_start; z < padsize_z; z++)
    {
      size_t z_src = (z - zd_start) + offset_z;
      if (z_src <= kernel_radius_z)
      {
        z_src = 0;
      }
      else
      {
        z_src -= kernel_radius_z;
      }

      z_src = astc::min(z_src, static_cast<size_t>(img->dim_z - 1));

      uint8_t* data8 = static_cast<uint8_t*>(img->data[z_src]);

      for (size_t y = 1; y < padsize_y; y++)
      {
        size_t y_src = (y - 1) + offset_y;
        if (y_src <= kernel_radius_xy)
        {
          y_src = 0;
        }
        else
        {
          y_src -= kernel_radius_xy;
        }

        y_src = astc::min(y_src, static_cast<size_t>(img->dim_y - 1));

        for (size_t x = 1; x < padsize_x; x++)
        {
          size_t x_src = (x - 1) + offset_x;
          if (x_src <= kernel_radius_xy)
          {
            x_src = 0;
          }
          else
          {
            x_src -= kernel_radius_xy;
          }

          x_src = astc::min(x_src, static_cast<size_t>(img->dim_x - 1));

          data[0] = data8[(4 * img->dim_x * y_src) + (4 * x_src    )];
          data[1] = data8[(4 * img->dim_x * y_src) + (4 * x_src + 1)];
          data[2] = data8[(4 * img->dim_x * y_src) + (4 * x_src + 2)];
          data[3] = data8[(4 * img->dim_x * y_src) + (4 * x_src + 3)];

          uint8_t r = data[swz.r];
          uint8_t g = data[swz.g];
          uint8_t b = data[swz.b];
          uint8_t a = data[swz.a];

          vfloat4 d = vfloat4 (r * (1.0f / 255.0f),
                               g * (1.0f / 255.0f),
                               b * (1.0f / 255.0f),
                               a * (1.0f / 255.0f));

          VARBUF1(z, y, x) = d;
          VARBUF2(z, y, x) = d * d;
        }
      }
    }
  }
  else if (img->data_type == ASTCENC_TYPE_F16)
  {
    // Swizzle data structure 4 = ZERO, 5 = ONE (in FP16)
    uint16_t data[6];
    data[ASTCENC_SWZ_0] = 0;
    data[ASTCENC_SWZ_1] = 0x3C00;

    for (size_t z = zd_start; z < padsize_z; z++)
    {
      size_t z_src = (z - zd_start) + offset_z;
      if (z_src <= kernel_radius_z)
      {
        z_src = 0;
      }
      else
      {
        z_src -= kernel_radius_z;
      }

      z_src = astc::min(z_src, static_cast<size_t>(img->dim_z - 1));

      uint16_t* data16 = static_cast<uint16_t*>(img->data[z_src]);

      for (size_t y = 1; y < padsize_y; y++)
      {
        size_t y_src = (y - 1) + offset_y;
        if (y_src <= kernel_radius_xy)
        {
          y_src = 0;
        }
        else
        {
          y_src -= kernel_radius_xy;
        }

        y_src = astc::min(y_src, static_cast<size_t>(img->dim_y - 1));

        for (size_t x = 1; x < padsize_x; x++)
        {
          size_t x_src = (x - 1) + offset_x;
          if (x_src <= kernel_radius_xy)
          {
            x_src = 0;
          }
          else
          {
            x_src -= kernel_radius_xy;
          }

          x_src = astc::min(x_src, static_cast<size_t>(img->dim_x - 1));

          data[0] = data16[(4 * img->dim_x * y_src) + (4 * x_src    )];
          data[1] = data16[(4 * img->dim_x * y_src) + (4 * x_src + 1)];
          data[2] = data16[(4 * img->dim_x * y_src) + (4 * x_src + 2)];
          data[3] = data16[(4 * img->dim_x * y_src) + (4 * x_src + 3)];

          vint4 di(data[swz.r], data[swz.g], data[swz.b], data[swz.a]);
          vfloat4 d = float16_to_float(di);

          VARBUF1(z, y, x) = d;
          VARBUF2(z, y, x) = d * d;
        }
      }
    }
  }
  else // if (img->data_type == ASTCENC_TYPE_F32)
  {
    assert(img->data_type == ASTCENC_TYPE_F32);

    // Swizzle data structure 4 = ZERO, 5 = ONE (in FP16)
    float data[6];
    data[ASTCENC_SWZ_0] = 0.0f;
    data[ASTCENC_SWZ_1] = 1.0f;

    for (size_t z = zd_start; z < padsize_z; z++)
    {
      size_t z_src = (z - zd_start) + offset_z;
      if (z_src <= kernel_radius_z)
      {
        z_src = 0;
      }
      else
      {
        z_src -= kernel_radius_z;
      }

      z_src = astc::min(z_src, static_cast<size_t>(img->dim_z - 1));

      float* data32 = static_cast<float*>(img->data[z_src]);

      for (size_t y = 1; y < padsize_y; y++)
      {
        size_t y_src = (y - 1) + offset_y;
        if (y_src <= kernel_radius_xy)
        {
          y_src = 0;
        }
        else
        {
          y_src -= kernel_radius_xy;
        }

        y_src = astc::min(y_src, static_cast<size_t>(img->dim_y - 1));

        for (size_t x = 1; x < padsize_x; x++)
        {
          size_t x_src = (x - 1) + offset_x;
          if (x_src <= kernel_radius_xy)
          {
            x_src = 0;
          }
          else
          {
            x_src -= kernel_radius_xy;
          }

          x_src = astc::min(x_src, static_cast<size_t>(img->dim_x - 1));

          data[0] = data32[(4 * img->dim_x * y_src) + (4 * x_src    )];
          data[1] = data32[(4 * img->dim_x * y_src) + (4 * x_src + 1)];
          data[2] = data32[(4 * img->dim_x * y_src) + (4 * x_src + 2)];
          data[3] = data32[(4 * img->dim_x * y_src) + (4 * x_src + 3)];

          float r = data[swz.r];
          float g = data[swz.g];
          float b = data[swz.b];
          float a = data[swz.a];

          vfloat4 d(r, g, b, a);

          VARBUF1(z, y, x) = d;
          VARBUF2(z, y, x) = d * d;
        }
      }
    }
  }

  // Pad with an extra layer of 0s; this forms the edge of the SAT tables
  vfloat4 vbz = vfloat4::zero();
  for (size_t z = 0; z < padsize_z; z++)
  {
    for (size_t y = 0; y < padsize_y; y++)
    {
      VARBUF1(z, y, 0) = vbz;
      VARBUF2(z, y, 0) = vbz;
    }

    for (size_t x = 0; x < padsize_x; x++)
    {
      VARBUF1(z, 0, x) = vbz;
      VARBUF2(z, 0, x) = vbz;
    }
  }

  if (have_z)
  {
    for (size_t y = 0; y < padsize_y; y++)
    {
      for (size_t x = 0; x < padsize_x; x++)
      {
        VARBUF1(0, y, x) = vbz;
        VARBUF2(0, y, x) = vbz;
      }
    }
  }

  // Generate summed-area tables for N and N^2; this is done in-place, using
  // a Brent-Kung parallel-prefix based algorithm to minimize precision loss
  for (size_t z = zd_start; z < padsize_z; z++)
  {
    for (size_t y = 1; y < padsize_y; y++)
    {
      brent_kung_prefix_sum(&(VARBUF1(z, y, 1)), padsize_x - 1, 1);
      brent_kung_prefix_sum(&(VARBUF2(z, y, 1)), padsize_x - 1, 1);
    }
  }

  for (size_t z = zd_start; z < padsize_z; z++)
  {
    for (size_t x = 1; x < padsize_x; x++)
    {
      brent_kung_prefix_sum(&(VARBUF1(z, 1, x)), padsize_y - 1, yst);
      brent_kung_prefix_sum(&(VARBUF2(z, 1, x)), padsize_y - 1, yst);
    }
  }

  if (have_z)
  {
    for (size_t y = 1; y < padsize_y; y++)
    {
      for (size_t x = 1; x < padsize_x; x++)
      {
        brent_kung_prefix_sum(&(VARBUF1(1, y, x)), padsize_z - 1, zst);
        brent_kung_prefix_sum(&(VARBUF2(1, y, x)), padsize_z - 1, zst);
      }
    }
  }

  // Compute a few constants used in the variance-calculation.
  float alpha_kdim = static_cast<float>(2 * alpha_kernel_radius + 1);
  float alpha_rsamples;

  if (have_z)
  {
    alpha_rsamples = 1.0f / (alpha_kdim * alpha_kdim * alpha_kdim);
  }
  else
  {
    alpha_rsamples = 1.0f / (alpha_kdim * alpha_kdim);
  }

  // Use the summed-area tables to compute variance for each neighborhood
  if (have_z)
  {
    for (size_t z = 0; z < size_z; z++)
    {
      size_t z_src = z + kernel_radius_z;
      size_t z_dst = z + offset_z;
      size_t z_low  = z_src - alpha_kernel_radius;
      size_t z_high = z_src + alpha_kernel_radius + 1;

      for (size_t y = 0; y < size_y; y++)
      {
        size_t y_src = y + kernel_radius_xy;
        size_t y_dst = y + offset_y;
        size_t y_low  = y_src - alpha_kernel_radius;
        size_t y_high = y_src + alpha_kernel_radius + 1;

        for (size_t x = 0; x < size_x; x++)
        {
          size_t x_src = x + kernel_radius_xy;
          size_t x_dst = x + offset_x;
          size_t x_low  = x_src - alpha_kernel_radius;
          size_t x_high = x_src + alpha_kernel_radius + 1;

          // Summed-area table lookups for alpha average
          float vasum = (  VARBUF1(z_high,  y_low,  x_low).lane<3>()
                         - VARBUF1(z_high,  y_low, x_high).lane<3>()
                         - VARBUF1(z_high, y_high,  x_low).lane<3>()
                         + VARBUF1(z_high, y_high, x_high).lane<3>()) -
                        (  VARBUF1(z_low,   y_low,  x_low).lane<3>()
                         - VARBUF1(z_low,   y_low, x_high).lane<3>()
                         - VARBUF1(z_low,  y_high,  x_low).lane<3>()
                         + VARBUF1(z_low,  y_high, x_high).lane<3>());

          size_t out_index = z_dst * zdt + y_dst * ydt + x_dst;
          input_alpha_averages[out_index] = (vasum * alpha_rsamples);
        }
      }
    }
  }
  else
  {
    for (size_t y = 0; y < size_y; y++)
    {
      size_t y_src = y + kernel_radius_xy;
      size_t y_dst = y + offset_y;
      size_t y_low  = y_src - alpha_kernel_radius;
      size_t y_high = y_src + alpha_kernel_radius + 1;

      for (size_t x = 0; x < size_x; x++)
      {
        size_t x_src = x + kernel_radius_xy;
        size_t x_dst = x + offset_x;
        size_t x_low  = x_src - alpha_kernel_radius;
        size_t x_high = x_src + alpha_kernel_radius + 1;

        // Summed-area table lookups for alpha average
        float vasum = VARBUF1(0, y_low,  x_low).lane<3>()
                    - VARBUF1(0, y_low,  x_high).lane<3>()
                    - VARBUF1(0, y_high, x_low).lane<3>()
                    + VARBUF1(0, y_high, x_high).lane<3>();

        size_t out_index = y_dst * ydt + x_dst;
        input_alpha_averages[out_index] = (vasum * alpha_rsamples);
      }
    }
  }
}

/* See header for documentation. */
size_t init_compute_averages(
  const astcenc_image& img,
  size_t alpha_kernel_radius,
  const astcenc_swizzle& swz,
  avg_args& ag
) {
  size_t size_x = img.dim_x;
  size_t size_y = img.dim_y;
  size_t size_z = img.dim_z;

  // Compute maximum block size and from that the working memory buffer size
  size_t kernel_radius = alpha_kernel_radius;
  size_t kerneldim = 2 * kernel_radius + 1;

  bool have_z = size_z > 1;
  size_t max_blk_size_xy = have_z ? 16u : 32u;

  size_t max_blk_size_z_lim = have_z ? 16u : 1u;
  size_t max_blk_size_z = astc::min(size_z, max_blk_size_z_lim);

  size_t max_padsize_xy = max_blk_size_xy + kerneldim;
  size_t max_padsize_z = max_blk_size_z + (have_z ? kerneldim : 0);

  // Perform block-wise averages calculations across the image
  // Initialize fields which are not populated until later
  ag.arg.size_x = 0;
  ag.arg.size_y = 0;
  ag.arg.size_z = 0;
  ag.arg.offset_x = 0;
  ag.arg.offset_y = 0;
  ag.arg.offset_z = 0;
  ag.arg.work_memory = nullptr;

  ag.arg.img = &img;
  ag.arg.swz = swz;
  ag.arg.have_z = have_z;
  ag.arg.alpha_kernel_radius = alpha_kernel_radius;

  ag.img_size_x = size_x;
  ag.img_size_y = size_y;
  ag.img_size_z = size_z;
  ag.blk_size_xy = max_blk_size_xy;
  ag.blk_size_z = max_blk_size_z;
  ag.work_memory_size = 2 * max_padsize_xy * max_padsize_xy * max_padsize_z;

  // The parallel task count
  size_t tasks_z = astc::get_block_count_safe(size_z, max_blk_size_z);
  size_t tasks_y = astc::get_block_count_safe(size_y, max_blk_size_xy);
  return tasks_z * tasks_y;
}

#endif

Coverage Report

Created: 2026-08-08 06:50

Line	Count	Source
1		// SPDX-License-Identifier: Apache-2.0
2		// ----------------------------------------------------------------------------
3		// Copyright 2011-2026 Arm Limited
4		//
5		// Licensed under the Apache License, Version 2.0 (the "License"); you may not
6		// use this file except in compliance with the License. You may obtain a copy
7		// of the License at:
8		//
9		// http://www.apache.org/licenses/LICENSE-2.0
10		//
11		// Unless required by applicable law or agreed to in writing, software
12		// distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
13		// WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
14		// License for the specific language governing permissions and limitations
15		// under the License.
16		// ----------------------------------------------------------------------------
17
18		#if !defined(ASTCENC_DECOMPRESS_ONLY)
19
20		/**
21		* @brief Functions to calculate variance per component in a NxN footprint.
22		*
23		* We need N to be parametric, so the routine below uses summed area tables in order to execute in
24		* O(1) time independent of how big N is.
25		*
26		* The addition uses a Brent-Kung-based parallel prefix adder. This uses the prefix tree to first
27		* perform a binary reduction, and then distributes the results. This method means that there is no
28		* serial dependency between a given element and the next one, and also significantly improves
29		* numerical stability allowing us to use floats rather than doubles.
30		*/
31
32		#include "astcenc_internal.h"
33
34		#include <cassert>
35
36		/**
37		* @brief Generate a prefix-sum array using the Brent-Kung algorithm.
38		*
39		* This will take an input array of the form:
40		* v0, v1, v2, ...
41		* ... and modify in-place to turn it into a prefix-sum array of the form:
42		* v0, v0+v1, v0+v1+v2, ...
43		*
44		* @param d The array to prefix-sum.
45		* @param items The number of items in the array.
46		* @param stride The item spacing in the array; i.e. dense arrays should use 1.
47		*/
48		static void brent_kung_prefix_sum(
49		vfloat4* d,
50		size_t items,
51		size_t stride
52	0	) {
53	0	if (items < 2)
54	0	return;
55
56	0	size_t lc_stride = 2;
57	0	size_t log2_stride = 1;
58
59		// The reduction-tree loop
60	0	do {
61	0	size_t step = lc_stride >> 1;
62	0	size_t start = lc_stride - 1;
63	0	size_t iters = items >> log2_stride;
64
65	0	vfloat4 da = d + (start stride);
66	0	ptrdiff_t ofs = -static_cast<ptrdiff_t>(step * stride);
67	0	size_t ofs_stride = stride << log2_stride;
68
69	0	while (iters)
70	0	{
71	0	da = da + da[ofs];
72	0	da += ofs_stride;
73	0	iters--;
74	0	}
75
76	0	log2_stride += 1;
77	0	lc_stride <<= 1;
78	0	} while (lc_stride <= items);
79
80		// The expansion-tree loop
81	0	do {
82	0	log2_stride -= 1;
83	0	lc_stride >>= 1;
84
85	0	size_t step = lc_stride >> 1;
86	0	size_t start = step + lc_stride - 1;
87	0	size_t iters = (items - step) >> log2_stride;
88
89	0	vfloat4 da = d + (start stride);
90	0	ptrdiff_t ofs = -static_cast<ptrdiff_t>(step * stride);
91	0	size_t ofs_stride = stride << log2_stride;
92
93	0	while (iters)
94	0	{
95	0	da = da + da[ofs];
96	0	da += ofs_stride;
97	0	iters--;
98	0	}
99	0	} while (lc_stride > 2);
100	0	}
101
102		/* See header for documentation. */
103		void compute_pixel_region_variance(
104		astcenc_contexti& ctx,
105		const pixel_region_args& arg
106	0	) {
107		// Unpack the memory structure into local variables
108	0	const astcenc_image* img = arg.img;
109	0	astcenc_swizzle swz = arg.swz;
110	0	bool have_z = arg.have_z;
111
112	0	size_t size_x = arg.size_x;
113	0	size_t size_y = arg.size_y;
114	0	size_t size_z = arg.size_z;
115
116	0	size_t offset_x = arg.offset_x;
117	0	size_t offset_y = arg.offset_y;
118	0	size_t offset_z = arg.offset_z;
119
120	0	size_t alpha_kernel_radius = arg.alpha_kernel_radius;
121
122	0	float* input_alpha_averages = ctx.input_alpha_averages;
123	0	vfloat4* work_memory = arg.work_memory;
124
125		// Compute memory sizes and dimensions that we need
126	0	size_t kernel_radius = alpha_kernel_radius;
127	0	size_t kerneldim = 2 * kernel_radius + 1;
128	0	size_t kernel_radius_xy = kernel_radius;
129	0	size_t kernel_radius_z = have_z ? kernel_radius : 0;
130
131	0	size_t padsize_x = size_x + kerneldim;
132	0	size_t padsize_y = size_y + kerneldim;
133	0	size_t padsize_z = size_z + (have_z ? kerneldim : 0);
134	0	size_t sizeprod = padsize_x * padsize_y * padsize_z;
135
136	0	int zd_start = have_z ? 1 : 0;
137
138	0	vfloat4 *varbuf1 = work_memory;
139	0	vfloat4 *varbuf2 = work_memory + sizeprod;
140
141		// Scaling factors to apply to Y and Z for accesses into the work buffers
142	0	size_t yst = padsize_x;
143	0	size_t zst = padsize_x * padsize_y;
144
145		// Scaling factors to apply to Y and Z for accesses into result buffers
146	0	size_t ydt = img->dim_x;
147	0	size_t zdt = img->dim_x * img->dim_y;
148
149		// Macros to act as accessor functions for the work-memory
150	0	#define VARBUF1(z, y, x) varbuf1[z * zst + y * yst + x]
151	0	#define VARBUF2(z, y, x) varbuf2[z * zst + y * yst + x]
152
153		// Load N and N^2 values into the work buffers
154	0	if (img->data_type == ASTCENC_TYPE_U8)
155	0	{
156		// Swizzle data structure 4 = ZERO, 5 = ONE
157	0	uint8_t data[6];
158	0	data[ASTCENC_SWZ_0] = 0;
159	0	data[ASTCENC_SWZ_1] = 255;
160
161	0	for (size_t z = zd_start; z < padsize_z; z++)
162	0	{
163	0	size_t z_src = (z - zd_start) + offset_z;
164	0	if (z_src <= kernel_radius_z)
165	0	{
166	0	z_src = 0;
167	0	}
168	0	else
169	0	{
170	0	z_src -= kernel_radius_z;
171	0	}
172
173	0	z_src = astc::min(z_src, static_cast<size_t>(img->dim_z - 1));
174
175	0	uint8_t* data8 = static_cast<uint8_t*>(img->data[z_src]);
176
177	0	for (size_t y = 1; y < padsize_y; y++)
178	0	{
179	0	size_t y_src = (y - 1) + offset_y;
180	0	if (y_src <= kernel_radius_xy)
181	0	{
182	0	y_src = 0;
183	0	}
184	0	else
185	0	{
186	0	y_src -= kernel_radius_xy;
187	0	}
188
189	0	y_src = astc::min(y_src, static_cast<size_t>(img->dim_y - 1));
190
191	0	for (size_t x = 1; x < padsize_x; x++)
192	0	{
193	0	size_t x_src = (x - 1) + offset_x;
194	0	if (x_src <= kernel_radius_xy)
195	0	{
196	0	x_src = 0;
197	0	}
198	0	else
199	0	{
200	0	x_src -= kernel_radius_xy;
201	0	}
202
203	0	x_src = astc::min(x_src, static_cast<size_t>(img->dim_x - 1));
204
205	0	data[0] = data8[(4 * img->dim_x * y_src) + (4 * x_src )];
206	0	data[1] = data8[(4 * img->dim_x * y_src) + (4 * x_src + 1)];
207	0	data[2] = data8[(4 * img->dim_x * y_src) + (4 * x_src + 2)];
208	0	data[3] = data8[(4 * img->dim_x * y_src) + (4 * x_src + 3)];
209
210	0	uint8_t r = data[swz.r];
211	0	uint8_t g = data[swz.g];
212	0	uint8_t b = data[swz.b];
213	0	uint8_t a = data[swz.a];
214
215	0	vfloat4 d = vfloat4 (r * (1.0f / 255.0f),
216	0	g * (1.0f / 255.0f),
217	0	b * (1.0f / 255.0f),
218	0	a * (1.0f / 255.0f));
219
220	0	VARBUF1(z, y, x) = d;
221	0	VARBUF2(z, y, x) = d * d;
222	0	}
223	0	}
224	0	}
225	0	}
226	0	else if (img->data_type == ASTCENC_TYPE_F16)
227	0	{
228		// Swizzle data structure 4 = ZERO, 5 = ONE (in FP16)
229	0	uint16_t data[6];
230	0	data[ASTCENC_SWZ_0] = 0;
231	0	data[ASTCENC_SWZ_1] = 0x3C00;
232
233	0	for (size_t z = zd_start; z < padsize_z; z++)
234	0	{
235	0	size_t z_src = (z - zd_start) + offset_z;
236	0	if (z_src <= kernel_radius_z)
237	0	{
238	0	z_src = 0;
239	0	}
240	0	else
241	0	{
242	0	z_src -= kernel_radius_z;
243	0	}
244
245	0	z_src = astc::min(z_src, static_cast<size_t>(img->dim_z - 1));
246
247	0	uint16_t* data16 = static_cast<uint16_t*>(img->data[z_src]);
248
249	0	for (size_t y = 1; y < padsize_y; y++)
250	0	{
251	0	size_t y_src = (y - 1) + offset_y;
252	0	if (y_src <= kernel_radius_xy)
253	0	{
254	0	y_src = 0;
255	0	}
256	0	else
257	0	{
258	0	y_src -= kernel_radius_xy;
259	0	}
260
261	0	y_src = astc::min(y_src, static_cast<size_t>(img->dim_y - 1));
262
263	0	for (size_t x = 1; x < padsize_x; x++)
264	0	{
265	0	size_t x_src = (x - 1) + offset_x;
266	0	if (x_src <= kernel_radius_xy)
267	0	{
268	0	x_src = 0;
269	0	}
270	0	else
271	0	{
272	0	x_src -= kernel_radius_xy;
273	0	}
274
275	0	x_src = astc::min(x_src, static_cast<size_t>(img->dim_x - 1));
276
277	0	data[0] = data16[(4 * img->dim_x * y_src) + (4 * x_src )];
278	0	data[1] = data16[(4 * img->dim_x * y_src) + (4 * x_src + 1)];
279	0	data[2] = data16[(4 * img->dim_x * y_src) + (4 * x_src + 2)];
280	0	data[3] = data16[(4 * img->dim_x * y_src) + (4 * x_src + 3)];
281
282	0	vint4 di(data[swz.r], data[swz.g], data[swz.b], data[swz.a]);
283	0	vfloat4 d = float16_to_float(di);
284
285	0	VARBUF1(z, y, x) = d;
286	0	VARBUF2(z, y, x) = d * d;
287	0	}
288	0	}
289	0	}
290	0	}
291	0	else // if (img->data_type == ASTCENC_TYPE_F32)
292	0	{
293	0	assert(img->data_type == ASTCENC_TYPE_F32);
294
295		// Swizzle data structure 4 = ZERO, 5 = ONE (in FP16)
296	0	float data[6];
297	0	data[ASTCENC_SWZ_0] = 0.0f;
298	0	data[ASTCENC_SWZ_1] = 1.0f;
299
300	0	for (size_t z = zd_start; z < padsize_z; z++)
301	0	{
302	0	size_t z_src = (z - zd_start) + offset_z;
303	0	if (z_src <= kernel_radius_z)
304	0	{
305	0	z_src = 0;
306	0	}
307	0	else
308	0	{
309	0	z_src -= kernel_radius_z;
310	0	}
311
312	0	z_src = astc::min(z_src, static_cast<size_t>(img->dim_z - 1));
313
314	0	float* data32 = static_cast<float*>(img->data[z_src]);
315
316	0	for (size_t y = 1; y < padsize_y; y++)
317	0	{
318	0	size_t y_src = (y - 1) + offset_y;
319	0	if (y_src <= kernel_radius_xy)
320	0	{
321	0	y_src = 0;
322	0	}
323	0	else
324	0	{
325	0	y_src -= kernel_radius_xy;
326	0	}
327
328	0	y_src = astc::min(y_src, static_cast<size_t>(img->dim_y - 1));
329
330	0	for (size_t x = 1; x < padsize_x; x++)
331	0	{
332	0	size_t x_src = (x - 1) + offset_x;
333	0	if (x_src <= kernel_radius_xy)
334	0	{
335	0	x_src = 0;
336	0	}
337	0	else
338	0	{
339	0	x_src -= kernel_radius_xy;
340	0	}
341
342	0	x_src = astc::min(x_src, static_cast<size_t>(img->dim_x - 1));
343
344	0	data[0] = data32[(4 * img->dim_x * y_src) + (4 * x_src )];
345	0	data[1] = data32[(4 * img->dim_x * y_src) + (4 * x_src + 1)];
346	0	data[2] = data32[(4 * img->dim_x * y_src) + (4 * x_src + 2)];
347	0	data[3] = data32[(4 * img->dim_x * y_src) + (4 * x_src + 3)];
348
349	0	float r = data[swz.r];
350	0	float g = data[swz.g];
351	0	float b = data[swz.b];
352	0	float a = data[swz.a];
353
354	0	vfloat4 d(r, g, b, a);
355
356	0	VARBUF1(z, y, x) = d;
357	0	VARBUF2(z, y, x) = d * d;
358	0	}
359	0	}
360	0	}
361	0	}
362
363		// Pad with an extra layer of 0s; this forms the edge of the SAT tables
364	0	vfloat4 vbz = vfloat4::zero();
365	0	for (size_t z = 0; z < padsize_z; z++)
366	0	{
367	0	for (size_t y = 0; y < padsize_y; y++)
368	0	{
369	0	VARBUF1(z, y, 0) = vbz;
370	0	VARBUF2(z, y, 0) = vbz;
371	0	}
372
373	0	for (size_t x = 0; x < padsize_x; x++)
374	0	{
375	0	VARBUF1(z, 0, x) = vbz;
376	0	VARBUF2(z, 0, x) = vbz;
377	0	}
378	0	}
379
380	0	if (have_z)
381	0	{
382	0	for (size_t y = 0; y < padsize_y; y++)
383	0	{
384	0	for (size_t x = 0; x < padsize_x; x++)
385	0	{
386	0	VARBUF1(0, y, x) = vbz;
387	0	VARBUF2(0, y, x) = vbz;
388	0	}
389	0	}
390	0	}
391
392		// Generate summed-area tables for N and N^2; this is done in-place, using
393		// a Brent-Kung parallel-prefix based algorithm to minimize precision loss
394	0	for (size_t z = zd_start; z < padsize_z; z++)
395	0	{
396	0	for (size_t y = 1; y < padsize_y; y++)
397	0	{
398	0	brent_kung_prefix_sum(&(VARBUF1(z, y, 1)), padsize_x - 1, 1);
399	0	brent_kung_prefix_sum(&(VARBUF2(z, y, 1)), padsize_x - 1, 1);
400	0	}
401	0	}
402
403	0	for (size_t z = zd_start; z < padsize_z; z++)
404	0	{
405	0	for (size_t x = 1; x < padsize_x; x++)
406	0	{
407	0	brent_kung_prefix_sum(&(VARBUF1(z, 1, x)), padsize_y - 1, yst);
408	0	brent_kung_prefix_sum(&(VARBUF2(z, 1, x)), padsize_y - 1, yst);
409	0	}
410	0	}
411
412	0	if (have_z)
413	0	{
414	0	for (size_t y = 1; y < padsize_y; y++)
415	0	{
416	0	for (size_t x = 1; x < padsize_x; x++)
417	0	{
418	0	brent_kung_prefix_sum(&(VARBUF1(1, y, x)), padsize_z - 1, zst);
419	0	brent_kung_prefix_sum(&(VARBUF2(1, y, x)), padsize_z - 1, zst);
420	0	}
421	0	}
422	0	}
423
424		// Compute a few constants used in the variance-calculation.
425	0	float alpha_kdim = static_cast<float>(2 * alpha_kernel_radius + 1);
426	0	float alpha_rsamples;
427
428	0	if (have_z)
429	0	{
430	0	alpha_rsamples = 1.0f / (alpha_kdim * alpha_kdim * alpha_kdim);
431	0	}
432	0	else
433	0	{
434	0	alpha_rsamples = 1.0f / (alpha_kdim * alpha_kdim);
435	0	}
436
437		// Use the summed-area tables to compute variance for each neighborhood
438	0	if (have_z)
439	0	{
440	0	for (size_t z = 0; z < size_z; z++)
441	0	{
442	0	size_t z_src = z + kernel_radius_z;
443	0	size_t z_dst = z + offset_z;
444	0	size_t z_low = z_src - alpha_kernel_radius;
445	0	size_t z_high = z_src + alpha_kernel_radius + 1;
446
447	0	for (size_t y = 0; y < size_y; y++)
448	0	{
449	0	size_t y_src = y + kernel_radius_xy;
450	0	size_t y_dst = y + offset_y;
451	0	size_t y_low = y_src - alpha_kernel_radius;
452	0	size_t y_high = y_src + alpha_kernel_radius + 1;
453
454	0	for (size_t x = 0; x < size_x; x++)
455	0	{
456	0	size_t x_src = x + kernel_radius_xy;
457	0	size_t x_dst = x + offset_x;
458	0	size_t x_low = x_src - alpha_kernel_radius;
459	0	size_t x_high = x_src + alpha_kernel_radius + 1;
460
461		// Summed-area table lookups for alpha average
462	0	float vasum = ( VARBUF1(z_high, y_low, x_low).lane<3>()
463	0	- VARBUF1(z_high, y_low, x_high).lane<3>()
464	0	- VARBUF1(z_high, y_high, x_low).lane<3>()
465	0	+ VARBUF1(z_high, y_high, x_high).lane<3>()) -
466	0	( VARBUF1(z_low, y_low, x_low).lane<3>()
467	0	- VARBUF1(z_low, y_low, x_high).lane<3>()
468	0	- VARBUF1(z_low, y_high, x_low).lane<3>()
469	0	+ VARBUF1(z_low, y_high, x_high).lane<3>());
470
471	0	size_t out_index = z_dst * zdt + y_dst * ydt + x_dst;
472	0	input_alpha_averages[out_index] = (vasum * alpha_rsamples);
473	0	}
474	0	}
475	0	}
476	0	}
477	0	else
478	0	{
479	0	for (size_t y = 0; y < size_y; y++)
480	0	{
481	0	size_t y_src = y + kernel_radius_xy;
482	0	size_t y_dst = y + offset_y;
483	0	size_t y_low = y_src - alpha_kernel_radius;
484	0	size_t y_high = y_src + alpha_kernel_radius + 1;
485
486	0	for (size_t x = 0; x < size_x; x++)
487	0	{
488	0	size_t x_src = x + kernel_radius_xy;
489	0	size_t x_dst = x + offset_x;
490	0	size_t x_low = x_src - alpha_kernel_radius;
491	0	size_t x_high = x_src + alpha_kernel_radius + 1;
492
493		// Summed-area table lookups for alpha average
494	0	float vasum = VARBUF1(0, y_low, x_low).lane<3>()
495	0	- VARBUF1(0, y_low, x_high).lane<3>()
496	0	- VARBUF1(0, y_high, x_low).lane<3>()
497	0	+ VARBUF1(0, y_high, x_high).lane<3>();
498
499	0	size_t out_index = y_dst * ydt + x_dst;
500	0	input_alpha_averages[out_index] = (vasum * alpha_rsamples);
501	0	}
502	0	}
503	0	}
504	0	}
505
506		/* See header for documentation. */
507		size_t init_compute_averages(
508		const astcenc_image& img,
509		size_t alpha_kernel_radius,
510		const astcenc_swizzle& swz,
511		avg_args& ag
512	0	) {
513	0	size_t size_x = img.dim_x;
514	0	size_t size_y = img.dim_y;
515	0	size_t size_z = img.dim_z;
516
517		// Compute maximum block size and from that the working memory buffer size
518	0	size_t kernel_radius = alpha_kernel_radius;
519	0	size_t kerneldim = 2 * kernel_radius + 1;
520
521	0	bool have_z = size_z > 1;
522	0	size_t max_blk_size_xy = have_z ? 16u : 32u;
523
524	0	size_t max_blk_size_z_lim = have_z ? 16u : 1u;
525	0	size_t max_blk_size_z = astc::min(size_z, max_blk_size_z_lim);
526
527	0	size_t max_padsize_xy = max_blk_size_xy + kerneldim;
528	0	size_t max_padsize_z = max_blk_size_z + (have_z ? kerneldim : 0);
529
530		// Perform block-wise averages calculations across the image
531		// Initialize fields which are not populated until later
532	0	ag.arg.size_x = 0;
533	0	ag.arg.size_y = 0;
534	0	ag.arg.size_z = 0;
535	0	ag.arg.offset_x = 0;
536	0	ag.arg.offset_y = 0;
537	0	ag.arg.offset_z = 0;
538	0	ag.arg.work_memory = nullptr;
539
540	0	ag.arg.img = &img;
541	0	ag.arg.swz = swz;
542	0	ag.arg.have_z = have_z;
543	0	ag.arg.alpha_kernel_radius = alpha_kernel_radius;
544
545	0	ag.img_size_x = size_x;
546	0	ag.img_size_y = size_y;
547	0	ag.img_size_z = size_z;
548	0	ag.blk_size_xy = max_blk_size_xy;
549	0	ag.blk_size_z = max_blk_size_z;
550	0	ag.work_memory_size = 2 * max_padsize_xy * max_padsize_xy * max_padsize_z;
551
552		// The parallel task count
553	0	size_t tasks_z = astc::get_block_count_safe(size_z, max_blk_size_z);
554	0	size_t tasks_y = astc::get_block_count_safe(size_y, max_blk_size_xy);
555	0	return tasks_z * tasks_y;
556	0	}
557
558		#endif