Adapted from Moreland's InterpolateColor. This uses the full
    transformations to compute a color mapping at a set of sampled points.

def moreland_interpolate_sampled (srgb1, srgb2, refwhite=CIELAB_D65,
                                  nsamples=DEFAULT_SAMPLE_POINTS):
    """Adapted from Moreland's InterpolateColor. This uses the full
    transformations to compute a color mapping at a set of sampled points."""

    msh1, msh2 = srgb_to_msh (np.asarray ([srgb1, srgb2],
                                          dtype=np.float), refwhite)

    raddiff = msh1[H] - msh2[H]
    while raddiff > np.pi:
        raddiff -= 2 * np.pi
    while raddiff < -np.pi:
        raddiff += 2 * np.pi
    raddiff = np.abs (raddiff)

    x = np.linspace (0, 1, nsamples).reshape ((nsamples, 1))
    x = np.repeat (x, 3, 1)

    if msh1[S] <= 0.05 or msh2[S] <= 0.05 or raddiff < np.pi/3:
        # Colors are too close together to comfortably put white in between.
        # Our interpolation won't have a control point, and won't actually be
        # divergent.

        if msh1[S] < 0.05 and msh2[S] > 0.05:
            msh1[H] = moreland_adjusthue (msh1, msh1[M])
        elif msh2[S] < 0.05 and msh1[S] > 0.05:
            msh2[H] = moreland_adjusthue (msh2, msh2[M])

        samples = np.empty ((4, nsamples))

        msh = (1 - x) * msh1 + x * msh2
        samples[0] = x[:,0]
        samples[1:4] = msh_to_srgb (msh, refwhite).T
    else:
        # Colors are not too close together -- we can add a white control
        # point in the middle, and do two interpolations joined piecewise. We
        # then use 2*nsamples-1 (not actually nsamples -- shhh) samples for
        # the spline fit

        msh3 = msh2
        msh2a = np.asarray ([np.max ([msh1[M], msh3[M], 88]), 0, 0])
        msh2b = msh2a.copy ()

        if msh1[S] < 0.05 and msh2a[S] > 0.05:
            msh1[H] = moreland_adjusthue (msh2a, msh1[M])
        elif msh2a[S] < 0.05 and msh1[S] > 0.05:
            msh2a[H] = moreland_adjusthue (msh1, msh2a[M])

        if msh2b[S] < 0.05 and msh3[S] > 0.05:
            msh2b[H] = moreland_adjusthue (msh3, msh2b[M])
        elif msh3[S] < 0.05 and msh2b[S] > 0.05:
            msh3[H] = moreland_adjusthue (msh2b, msh3[M])

        samples = np.empty ((4, 2*nsamples-1))

        msh = (1 - x) * msh1 + x * msh2a
        samples[0,:nsamples] = 0.5 * x[:,0]
        samples[1:4,:nsamples] = msh_to_srgb (msh, refwhite).T

        msh = (1 - x) * msh2b + x * msh3
        samples[0,nsamples-1:] = 0.5 * x[:,0] + 0.5
        samples[1:4,nsamples-1:] = msh_to_srgb (msh, refwhite).T

    return samples