def normalize_RGB_image(RGB, M):
  # Assumes that {RGB} is an image with positive and negative samples
  # and arbitrary mean, with shape {(ny,nx,nc)} where {nc} (number of
  # channels).
  # 
  # Then the procedure scales and shifts the samples in each channel, so
  # that those samples where the mask image {M} is nonzero are
  # normalized. The image {M} must be an array with shape {(ny,nx)}
  # whose values are either 0 or 1.
  #
  # If {ball} is true, the normalization makes the selected samples have
  # zero mean and unit variance. If {ball} is false, the normalization
  # ensures that those samples span the {{0 _ 1]} interval, apart from a
  # small percentage; and clips those which lie outside that interval.
  #
  # The result is a {numpy} array {RGB} of floats with shape
  # {(ny,nx,nc)} with float samples in {[0 _ 1]}.

  ny, nx, nc = numpy.shape(RGB)
  assert numpy.shape(M) == (ny, nx)
      
  # Scale and shift each channel:
  for ic in range(nc):
    cvals = [ \
      RGB[iy,ix,ic] \
      for ix in range(nx) \
      for iy in range(ny) \
      if M[iy,ix] != 0 \
    ]
    
    if ball:
      savg = numpy.mean(cvals)
      sdev = numpy.std(cvals) # New {numpy} allows {mean = savg} arg.
      # Normalize to zero mean,unit variance, no clipping:
      for iy in range(ny):
        for ix in range(nx):
          RGB[iy,ix,ic] = (RGB[iy,ix,ic] - savg)/sdev 
    else:
      # Normalize to {95%} of samples spanning {[0_1]}, with clipping:
      smin = numpy.quantile(cvals, 0.025)
      smax = numpy.quantile(cvals, 0.975)
      sys.stderr.write(f"quantile range(channel {ic}):")
      sys.stderr.write(f" [{smin:7.4f} _ {smax:7.4f}]")
      smid = (smin + smax)/2
      srad = (smax - smin)/2

      nbig = 0
      fsmp_min = +inf
      fsmp_max = -inf
      for iy in range(ny):
        for ix in range(nx):
          fsmp = (RGB[iy,ix,ic] - smid)/srad 
          fsmp = 0.5*(1 + fsmp)
          if fsmp < 0.0 or fsmp > 1.0:
            nbig += 1
            fsmp_min = min(fsmp_min, fsmp)
            fsmp_max = max(fsmp_max, fsmp)
          fsmp = min(1.0, max(0.0, fsmp))
          RGB[iy,ix,ic] = fsmp 
      sys.stderr.write(f"channel {ic}: {nbig:4d} out of range")
      sys.stderr.write(f" actual range [{fsmp_min:7.4f} _ {fsmp_max:7.4f}]\n")
  return RGB
  # ----------------------------------------------------------------------