docs: def backgroundCMD(self, catalog, mode='cloud-in-cells', weights=None):
          """
          Generate an empirical background model in color-magnitude space.
          
          INPUTS:
              catalog: Catalog object
          OUTPUTS:
              background
          """
  
          # Select objects in annulus
          cut_annulus = self.roi.inAnnulus(catalog.lon,catalog.lat)
          color = catalog.color[cut_annulus]
          mag   = catalog.mag[cut_annulus]
  
          # Units are (deg^2)
          solid_angle = ugali.utils.binning.take2D(self.solid_angle_cmd, color, mag,
                                                   self.roi.bins_color, self.roi.bins_mag)
  
          # Weight each object before binning
          # Divide by solid angle and bin size in magnitudes to get number density 
          # [objs / deg^2 / mag^2]
          if weights is None:
              number_density = (solid_angle*self.roi.delta_color*self.roi.delta_mag)**(-1)
          else:
              number_density = weights*(solid_angle*self.roi.delta_color*self.roi.delta_mag)**(-1)
  
          mode = str(mode).lower()
          if mode == 'cloud-in-cells':
              # Apply cloud-in-cells algorithm
              cmd_background = ugali.utils.binning.cloudInCells(color,mag,
                                                                [self.roi.bins_color,self.roi.bins_mag],
                                                                weights=number_density)[0]
          elif mode == 'bootstrap':
              # Not implemented
              raise ValueError("Bootstrap mode not implemented")
              mag_1_array = catalog.mag_1
              mag_2_array = catalog.mag_2
  
              catalog.mag_1 + (catalog.mag_1_err * np.random.normal(0, 1., len(catalog.mag_1)))
              catalog.mag_2 + (catalog.mag_2_err * np.random.normal(0, 1., len(catalog.mag_2)))
  
          elif mode == 'histogram':
              # Apply raw histogram
              cmd_background = np.histogram2d(mag,color,bins=[self.roi.bins_mag,self.roi.bins_color],
                                              weights=number_density)[0]
  
          elif mode == 'kde':
              # Gridded kernel density estimator
              logger.warning("### KDE not implemented properly")
              cmd_background = ugali.utils.binning.kernelDensity(color,mag,
                                                                 [self.roi.bins_color,self.roi.bins_mag],
                                                                 weights=number_density)[0]
          elif mode == 'uniform':
              logger.warning("### WARNING: Uniform CMD")
              hist = np.histogram2d(mag,color,bins=[self.roi.bins_mag,self.roi.bins_color], weights=number_density)[0]
              cmd_background = np.mean(hist)*np.ones(hist.shape)
              observable = (self.solid_angle_cmd > self.minimum_solid_angle)
              cmd_background *= observable
              return cmd_background
          else:
              raise ValueError("Unrecognized mode: %s"%mode)
          ## Account for the objects that spill out of the observable space
          ## But what about the objects that spill out to red colors??
          #for index_color in range(0, len(self.roi.centers_color)):
          #    for index_mag in range(0, len(self.roi.centers_mag)):
          #        if self.solid_angle_cmd[index_mag][index_color] < self.minimum_solid_angle:
          #            cmd_background[index_mag - 1][index_color] += cmd_background[index_mag][index_color]
          #            cmd_background[index_mag][index_color] = 0.
          #            break
  
          cmd_area = self.solid_angle_cmd*self.roi.delta_color*self.roi.delta_mag # [deg^2 * mag^2]
  
          # ADW: This accounts for leakage to faint magnitudes
          # But what about the objects that spill out to red colors??
          # Maximum obsevable magnitude index for each color (uses the fact that
          # np.argmin returns first minimum (zero) instance found.
          # NOTE: More complicated maps may have holes causing problems
  
          observable = (self.solid_angle_cmd > self.minimum_solid_angle)
          index_mag = observable.argmin(axis=0) - 1
          index_color = np.arange(len(self.roi.centers_color))
          # Add the cumulative leakage back into the last bin of the CMD
          leakage = (cmd_background * ~observable).sum(axis=0)
          cmd_background[[index_mag,index_color]] += leakage
          # Zero out all non-observable bins
          cmd_background *= observable
  
          # Avoid dividing by zero by setting empty bins to the value of the 
          # minimum filled bin of the CMD. This choice is arbitrary and 
          # could be replaced by a static minimum, some fraction of the 
          # CMD maximum, some median clipped minimum, etc. However, should 
          # be robust against outliers with very small values.
          min_cmd_background = max(cmd_background[cmd_background > 0.].min(),
                                   1e-4*cmd_background.max())
          cmd_background[observable] = cmd_background[observable].clip(min_cmd_background)
  
          ### # ADW: This is a fudge factor introduced to renormalize the CMD
          ### # to the number of input stars in the annulus. While leakage
          ### # will still smooth the distribution, it shouldn't result in 
          ### fudge_factor = len(mag) / float((cmd_background*cmd_area).sum())
          ### cmd_background *= fudge_factor
  
          return cmd_background