Context Definitions
The outer surface depicted in your Wikipedia image represents the volume bounded by the Colour Matching Functions (CMFS). The CMFS are defined by the CIE as follows:
tristimulus values of monochromatic stimuli of equal radiant power
Another name for the equal radiant power illuminant in the CIE system is Equal Energy or E.
The CIE XYZ tristimulus values are themselves defined by the CIE:
amounts of the 3 reference colour stimuli, in a given trichromatic system, required to match the colour of the stimulus considered
This is directly related to Colour Matching Experiments from Maxwell and Wright and Guild (1931) and probably out of the scope of this answer.
Outer Surface Generation
The outer surface is likely built using square wave spectral power distributions (SPDs), i.e. slices of the the Equal Energy illuminant. Bruce Lindbloom has an explanation on this page.
The methodology is to convert the square wave SPDs with various bandwidth to CIE XYZ tristimulus values, this will form the outer surface.
Assuming 5 bins, a first set of SPDs would be like that:
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
The second one:
1 1 0 0 0
0 1 1 0 0
0 0 1 1 0
0 0 0 1 1
1 0 0 0 1
The third:
1 1 1 0 0
0 1 1 1 0
0 0 1 1 1
1 0 0 1 1
1 1 0 0 1
etc...
Here a Python snippet using Colour to generate the points for the outer surface:
import colour
import numpy as np
class NearestNeighbourInterpolator(colour.KernelInterpolator):
def __init__(self, *args, **kwargs):
kwargs['kernel'] = colour.kernel_nearest_neighbour
super(NearestNeighbourInterpolator, self).__init__(*args, **kwargs)
def generate_square_waves(samples):
square_waves = []
square_waves_basis = np.tril(np.ones((samples, samples)))[0:-1, :]
for i in range(samples):
square_waves.append(np.roll(square_waves_basis, i))
return np.vstack((np.zeros(samples), np.vstack(square_waves),
np.ones(samples)))
def XYZ_outer_surface(samples):
XYZ = []
wavelengths = np.linspace(colour.DEFAULT_SPECTRAL_SHAPE.start,
colour.DEFAULT_SPECTRAL_SHAPE.end, samples)
for wave in generate_square_waves(samples):
spd = colour.SpectralPowerDistribution(
wave, wavelengths).align(
colour.DEFAULT_SPECTRAL_SHAPE,
interpolator=NearestNeighbourInterpolator)
XYZ.append(colour.spectral_to_XYZ(spd))
return np.array(XYZ).reshape(len(XYZ), -1, 3)
# 43 is picked as number of samples to have integer wavelengths.
colour.write_image(XYZ_outer_surface(43), 'CIE_XYZ_outer_surface.exr')
And the output: