It's histogram you are looking for, not index. Actually, indexND example pretty clear shows that it does something else. In your case, some selected elements from all-zeroes piddle became incremented from 0 to 1, that's all. Rather (demo output looks correct to me):

use strict;
use warnings;
use feature 'say';
use PDL;
use PDL::IO::Image;
use open IO => 'raw';

my $data = `convert -size 5x5 radial-gradient:red-green bmp:-`;
my $im = PDL::IO::Image-> new_from_file( \$data )
                       -> pixels_to_pdl
                       -> long;

my $flat = $im-> copy;
$flat-> slice( [],[],2 )-> inplace-> shiftleft( 16, 0 );
$flat-> slice( [],[],1 )-> inplace-> shiftleft( 8, 0 );

my $cube = $flat-> mv( 2, 0 )
                -> sumover
                -> flat
                -> hist( 0, 0x1000000, 1 )
                -> reshape( 256, 256, 256 );

### Demo

say $im;
say my $uniq_reds   = $im-> slice( [],[],0 )-> uniq;
say my $uniq_greens = $im-> slice( [],[],1 )-> uniq;
say $cube-> dice( $uniq_reds, $uniq_greens, 0 );


__END__

[
 [
  [  0   0   0   0   0]
  [  0  75 128  75   0]
  [  0 128 255 128   0]
  [  0  75 128  75   0]
  [  0   0   0   0   0]
 ]
 [
  [128 128 128 128 128]
  [128  91  64  91 128]
  [128  64   0  64 128]
  [128  91  64  91 128]
  [128 128 128 128 128]
 ]
 [
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
  [0 0 0 0 0]
 ]
]

[0 75 128 255]
[0 64 91 128]

[
 [
  [ 0  0  0  1]
  [ 0  0  4  0]
  [ 0  4  0  0]
  [16  0  0  0]
 ]
]
[download]

vr++ correctly identifies this is a histogram type of problem. This note is just about an alternative approach I stumbled on very recently.

The new PDL::IO::STL, in the imminent PDL 2.079, uses PDL::VectorValued if installed, to deduplicate its set of colour-vectors, then point the triangles' vertex-indices at the correct entry in the vertex-vector. It uses vsearchvec of PDL::VectorValued::Utils for this. Summarised from https://github.com/PDLPorters/pdl/blob/32a99b563358f7ad3ee26fba22fb24b5c02380e8/IO/STL/STL.pm#L150-L158:

use PDL::VectorValued::Utils;
# $pdl is ncoords=3,nvertices=3,ntriangles
my $uniqv = $pdl->uniqvec; # ncoords=3,nunique
my $face_indices = vsearchvec($pdl, $uniqv); # nvertices=3,ntriangles
[download]

From this, if you histogram the vertex-indices to get a count of each, you can then turn this into a pair of ndarrays ($coords, $count), and do

$cube->indexND($coords) .= $count;
[download]

This will work because then you are only writing into each cell once. As the "CAVEAT" in range of PDL::Slices notes:

It's quite possible to select multiple ranges that intersect. In that case, modifying the ranges doesn't have a guaranteed result in the original PDL -- the result is an arbitrary choice among the valid values. For some things that's OK; but for others it's not. In particular, this doesn't work:

pdl> $photon_list = new PDL::RandVar->sample(500)->reshape(2,250)*10
pdl> $histogram = zeroes(10,10)
pdl> $histogram->range($photon_list,1)++;  #not what you wanted
[download]

The reason is that if two photons land in the same bin, then that bin doesn't get incremented twice. (That may get fixed in a later version...)

Another alternative approach would be to try PDL::NDBin.