def voxelize(self, points):
        if not self.voxelized:
            # compute the boundary of the 3D points
            Xmin = np.nanmin(points[0,:]) - self.margin
            Xmax = np.nanmax(points[0,:]) + self.margin
            Ymin = np.nanmin(points[1,:]) - self.margin
            Ymax = np.nanmax(points[1,:]) + self.margin
            Zmin = np.nanmin(points[2,:]) - self.margin
            Zmax = np.nanmax(points[2,:]) + self.margin
            self.min_x = Xmin
            self.min_y = Ymin
            self.min_z = Zmin
            self.max_x = Xmax
            self.max_y = Ymax
            self.max_z = Zmax

            # step size
            self.step_x = (Xmax-Xmin) / self.grid_size
            self.step_y = (Ymax-Ymin) / self.grid_size
            self.step_z = (Zmax-Zmin) / self.grid_size
            self.voxelized = True

        # compute grid indexes
        indexes = np.zeros_like(points, dtype=np.float32)
        indexes[0,:] = np.floor((points[0,:] - self.min_x) / self.step_x)
        indexes[1,:] = np.floor((points[1,:] - self.min_y) / self.step_y)
        indexes[2,:] = np.floor((points[2,:] - self.min_z) / self.step_z)

        # crash the grid indexes
        # grid_indexes = indexes[0,:] * self.grid_size * self.grid_size + indexes[1,:] * self.grid_size + indexes[2,:]
        # I = np.isnan(grid_indexes)
        # grid_indexes[I] = -1
        # grid_indexes = grid_indexes.reshape(self.height, self.width).astype(np.int32)

        return indexes

def normalizeFloatImage3(floatImage):

    mn_0 = np.nanmin(np.nanmin(floatImage[:, :,1]))
    mx_0 = np.nanmax(np.nanmax(floatImage[:, :, 1]))

    rows = floatImage.shape[0]
    cols = floatImage.shape[1]
    mn_0 = 1000
    mx_0 =-1000
    for r in range(rows):
        for c in range(cols):
            if floatImage[r,c, 2] <= 0:
                mn_0 = mn_0 if mn_0 <=  floatImage[r,c, 1 ] else floatImage[r,c, 1 ]
                mx_0 = mx_0 if mx_0  >=  floatImage[r,c, 1 ] else floatImage[r,c, 1 ]

    fctr_0 = 255.0/(mx_0 - mn_0)

    for r in range(rows):
        for c in range(cols):
            if floatImage[r,c, 2]  <= 0:
                floatImage[r,c, 1] = (floatImage[r,c, 1] - mn_0) * fctr_0
            else:
                floatImage[r,c, 1]=0

    print "mn_0: ", mn_0, "  mx_0: ", mx_0

def rescale_to_depth_image(original_image, opencv_image):
    nan_max = np.nanmax(original_image.pixels[..., 2])
    nan_min = np.nanmin(original_image.pixels[..., 2])
    depth_pixels = opencv_image.astype(np.float) / 255.0
    depth_pixels = rescale(depth_pixels, 0.0, 1.0, nan_min, nan_max)
    depth_pixels[np.isclose(np.nanmin(depth_pixels), depth_pixels)] = np.nan
    return depth_pixels

    def _get_Tp_limits(self):
        """Get the limits for the graphs in temperature and pressure, based on 
        SI units: [Tmin, Tmax, pmin, pmax]"""
        T_lo,T_hi,P_lo,P_hi = self.limits
        Ts_lo,Ts_hi = self._get_sat_bounds(CoolProp.iT)
        Ps_lo,Ps_hi = self._get_sat_bounds(CoolProp.iP)

        if T_lo is None:            T_lo  = 0.0
        elif T_lo < self.ID_FACTOR: T_lo *= Ts_lo
        if T_hi is None:            T_hi  = 1e6
        elif T_hi < self.ID_FACTOR: T_hi *= Ts_hi
        if P_lo is None:            P_lo  = 0.0
        elif P_lo < self.ID_FACTOR: P_lo *= Ps_lo
        if P_hi is None:            P_hi  = 1e10
        elif P_hi < self.ID_FACTOR: P_hi *= Ps_hi

        try: T_lo = np.nanmax([T_lo, self._state.trivial_keyed_output(CoolProp.iT_min)])
        except: pass
        try: T_hi = np.nanmin([T_hi, self._state.trivial_keyed_output(CoolProp.iT_max)])
        except: pass
        try: P_lo = np.nanmax([P_lo, self._state.trivial_keyed_output(CoolProp.iP_min)])
        except: pass
        try: P_hi = np.nanmin([P_hi, self._state.trivial_keyed_output(CoolProp.iP_max)])
        except: pass

        return [T_lo,T_hi,P_lo,P_hi]

def calc_norm_summary_tables(accuracy_tbl, time_tbl):
    """
    Calculate normalized performance/ranking summary, as numpy
    matrices as usual for convenience, and matrices of additional
    statistics (min, max, percentiles, etc.)

    Here normalized means relative to the best which gets a 1, all
    others get the ratio resulting from dividing by the performance of
    the best.
    """
    # Min across all minimizers, i.e. for each fit problem what is the lowest chi-squared and the lowest time
    min_sum_err_sq = np.nanmin(accuracy_tbl, 1)
    min_runtime = np.nanmin(time_tbl, 1)

    # create normalised tables
    norm_acc_rankings = accuracy_tbl / min_sum_err_sq[:, None]
    norm_runtimes = time_tbl / min_runtime[:, None]

    summary_cells_acc = np.array([np.nanmin(norm_acc_rankings, 0),
                                  np.nanmax(norm_acc_rankings, 0),
                                  stats.nanmean(norm_acc_rankings, 0),
                                  stats.nanmedian(norm_acc_rankings, 0)
                                  ])

    summary_cells_runtime = np.array([np.nanmin(norm_runtimes, 0),
                                      np.nanmax(norm_runtimes, 0),
                                      stats.nanmean(norm_runtimes, 0),
                                      stats.nanmedian(norm_runtimes, 0)
                                      ])

    return norm_acc_rankings, norm_runtimes, summary_cells_acc, summary_cells_runtime

def plotLL(fname='out4.npy'):
    plt.figure()
    h= np.linspace(0,1,21)
    g= np.linspace(0,1,21)
    m=np.linspace(0,2,21)
    d=np.linspace(0,2,21)
    out=np.load(fname)
    print np.nanmax(out),np.nanmin(out)
    rang=np.nanmax(out)-np.nanmin(out)
    maxloc= np.squeeze(np.array((np.nanmax(out)==out).nonzero()))
    H,G=np.meshgrid(h,g)
    print maxloc
    for mm in range(m.size/2):
        for dd in range(d.size/2):
            plt.subplot(10,10,(9-mm)*10+dd+1)
            plt.pcolormesh(h,g,out[:,:,mm*2,dd*2].T,
                           vmax=np.nanmax(out),vmin=np.nanmax(out)-rang/4.)
            plt.gca().set_xticks([])
            plt.gca().set_yticks([])
            if mm==maxloc[2]/2 and dd==maxloc[3]/2:
                plt.plot(h[maxloc[0]],g[maxloc[1]],'ow',ms=8)
            if dd==0:
                print mm,dd
                plt.ylabel('%.1f'%m[mm*2])
            if mm==0: plt.xlabel('%.1f'%d[dd*2])
    plt.title(fname[:6])

def draw_hmap_old(hmap, yvals, fname=None):
    """
    Plot a matrix as a heat map and write an image file.
    :param hmap: Heat map matrix.
    :param yvals: Heat map Y labels (e.g. amino acid names).
    :param fname: Destination image file.
    """
    if np.nanmax(hmap) > abs(np.nanmin(hmap)):
        vmax = np.nanmax(hmap)
        vmin = -np.nanmax(hmap)
    else:
        vmax = abs(np.nanmin(hmap))
        vmin = np.nanmin(hmap)
    fig = plt.figure()
    plt.figure(figsize=(20,10))
    plt.imshow(hmap, cmap='RdBu', interpolation = 'nearest',aspect='auto',vmin = vmin ,vmax = vmax )
    plt.xlim(0, hmap.shape[1])
    plt.ylim(0, hmap.shape[0])
    ax = plt.gca()
    fig.set_facecolor('white')
    ax.set_xlim((-0.5, hmap.shape[1] -0.5))
    ax.set_ylim((-0.5, hmap.shape[0] -0.5))
    ax.set_yticks([x for x in xrange(0, hmap.shape[0])])
    ax.set_yticklabels(yvals)
    ax.set_xticks(range(0,76,5))
    ax.set_xticklabels(range(2,76,5)+['STOP'])
    ax.set_ylabel('Residue')
    ax.set_xlabel('Ub Sequence Position')
    cb = plt.colorbar()
    cb.set_clim(vmin=vmin, vmax=vmax)
    cb.set_label('Relative Fitness')
    if fname is not None:
        plt.savefig(fname, bbox_inches='tight')
    return fig

def TestPlot(fig=None):
    A = numpy.array([1,2,3,4,2,5,8,3,2,3,5,6])
    B = numpy.array([8,7,3,6,4,numpy.nan,9,3,7,numpy.nan,2,4])
    C = numpy.array([6,3,4,7,2,1,1,7,8,4,3,2])
    D = numpy.array([5,2,4,5,3,8,2,5,3,5,6,8])
    
    # A work around to get the histograms overplotted with each other to overlap correctly;
    histrangelist = [(numpy.nanmin(A),numpy.nanmax(A)),(numpy.nanmin(B),numpy.nanmax(B)),
                (numpy.nanmin(C),numpy.nanmax(C)),(numpy.nanmin(D),numpy.nanmax(D))]
    
    data = numpy.array([A,B,C,D])
    labels = ['A','3','C','D']

    fig = GridPlot(data,labels=labels, no_tick_labels=True, color='black', 
                    hist=True, histbins=3, histloc='tl', histrangelist=histrangelist, fig=None) 
    
    # Data of note to plot in different color
    A2 = numpy.array([1,2,3,4])
    B2 = numpy.array([8,7,3,6])
    C2 = numpy.array([6,3,4,7])
    D2 = numpy.array([5,2,4,5])
    data2 = numpy.array([A2,B2,C2,D2])
    
    fig = GridPlot(data2,labels=labels, no_tick_labels=True, color='red', 
                hist=True, histbins=3, histloc='tr', histrangelist=histrangelist, fig=fig) 
    
    return fig

def bin_fit(x, y, buckets=3):
     
    assert buckets in [3,25]

    xstd=np.nanstd(x)
    
    if buckets==3:
        binlimits=[np.nanmin(x), -xstd/2.0,xstd/2.0 , np.nanmax(x)]
    elif buckets==25:
    
        steps=xstd/4.0
        binlimits=np.arange(-xstd*3.0, xstd*3.0, steps)
    
        binlimits=[np.nanmin(x)]+list(binlimits)+[np.nanmax(x)]
    
    fit_y=[]
    err_y=[]
    x_values_to_plot=[]
    for binidx in range(len(binlimits))[1:]:
        lower_bin_x=binlimits[binidx-1]
        upper_bin_x=binlimits[binidx]

        x_values_to_plot.append(np.mean([lower_bin_x, upper_bin_x]))

        y_in_bin=[y[idx] for idx in range(len(y)) if x[idx]>=lower_bin_x and x[idx]<upper_bin_x]

        fit_y.append(np.nanmedian(y_in_bin))
        err_y.append(np.nanstd(y_in_bin))

    ## no zeros
    

    return (binlimits, x_values_to_plot, fit_y, err_y)

def plot_nontarget_betas_n_back(t_vols_n_back_beta_1, b_vols_smooth_n_back, in_brain_mask, brain_structure, nice_cmap, n_back):

  beta_index = 1

  b_vols_smooth_n_back[~in_brain_mask] = np.nan
  t_vols_n_back_beta_1[~in_brain_mask] = np.nan
  min_val = np.nanmin(b_vols_smooth_n_back[...,(40,50,60),beta_index])
  max_val = np.nanmax(b_vols_smooth_n_back[...,(40,50,60),beta_index])

  plt.figure()

  for map_index, depth in (((3,2,1), 40),((3,2,3), 50),((3,2,5), 60)):
    plt.subplot(*map_index)
    plt.title("z=%d,%s" % (depth, n_back + "-back nontarget,beta values"))
    plt.imshow(brain_structure[...,depth], alpha=0.5)
    plt.imshow(b_vols_smooth_n_back[...,depth,beta_index], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
    plt.colorbar()
    plt.tight_layout()

  t_min_val = np.nanmin(t_vols_n_back_beta_1[...,(40,50,60)])
  t_max_val = np.nanmax(t_vols_n_back_beta_1[...,(40,50,60)])

  for map_index, depth in (((3,2,2), 40),((3,2,4), 50),((3,2,6), 60)):
    plt.subplot(*map_index)
    plt.title("z=%d,%s" % (depth, n_back + "-back nontarget,t values"))
    plt.imshow(brain_structure[...,depth], alpha=0.5)
    plt.imshow(t_vols_n_back_beta_1[...,depth], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
    plt.colorbar()
    plt.tight_layout()

  plt.savefig(os.path.join(output_filename, "sub011_nontarget_betas_%s_back.png" % (n_back)), format='png', dpi=500)  

    def callback_function(self,
                          optimiser_output,
                          minimise_function_result,
                          was_accepted):
        if not(was_accepted):
            return

        if self.current_success_number == 0:
            # First result
            self.successful_results[0] = minimise_function_result
            self.current_success_number = 1

        elif (minimise_function_result >=
              np.nanmin(self.successful_results) + self.confidence):
            # Reject result
            0

        elif (minimise_function_result >=
              np.nanmin(self.successful_results) - self.confidence):
            # Agreeing result
            self.successful_results[
                self.current_success_number
            ] = minimise_function_result

            self.current_success_number += 1

        elif (minimise_function_result <
              np.nanmin(self.successful_results) - self.confidence):
            # New result
            self.successful_results[0] = minimise_function_result
            self.current_success_number = 1

        if self.current_success_number >= self.n:
            return True

def plot_noise_regressor_betas(b_vols_smooth, t_vols_beta_6_to_9, brain_structure, in_brain_mask, nice_cmap):

  plt.figure()

  min_val = np.nanmin(b_vols_smooth[...,40,(6,7,9)])
  max_val = np.nanmax(b_vols_smooth[...,40,(6,7,9)])

  plt.subplot(3,2,1)
  plt.title("z=%d,%s" % (40, "linear drift,betas"))
  b_vols_smooth[~in_brain_mask] = np.nan
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(b_vols_smooth[...,40,6], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,3)
  plt.title("z=%d,%s" % (40, "quadratic drift,betas"))
  b_vols_smooth[~in_brain_mask] = np.nan
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(b_vols_smooth[...,40,7], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,5)
  plt.title("z=%d,%s" % (40, "second PC,betas"))
  b_vols_smooth[~in_brain_mask] = np.nan
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(b_vols_smooth[...,40,9], cmap=nice_cmap, alpha=0.5, vmin=min_val, vmax=max_val)
  plt.colorbar()
  plt.tight_layout()

  t_vols_beta_6_to_9[0][~in_brain_mask] = np.nan
  t_vols_beta_6_to_9[1][~in_brain_mask] = np.nan
  t_vols_beta_6_to_9[3][~in_brain_mask] = np.nan

  t_min_val = np.nanmin([t_vols_beta_6_to_9[0][...,40], t_vols_beta_6_to_9[1][...,40], t_vols_beta_6_to_9[3][...,40]])
  t_max_val = np.nanmax([t_vols_beta_6_to_9[0][...,40], t_vols_beta_6_to_9[1][...,40], t_vols_beta_6_to_9[3][...,40]])

  plt.subplot(3,2,2)
  plt.title("z=%d,%s" % (40, "linear drift,t values"))
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(t_vols_beta_6_to_9[0][...,40], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,4)
  plt.title("z=%d,%s" % (40, "quadratic drift,t values"))
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(t_vols_beta_6_to_9[1][...,40], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.subplot(3,2,6)
  plt.title("z=%d,%s" % (40, "second PC,t values"))
  plt.imshow(brain_structure[...,40], alpha=0.5)
  plt.imshow(t_vols_beta_6_to_9[3][...,40], cmap=nice_cmap, alpha=0.5, vmin=t_min_val, vmax=t_max_val)
  plt.colorbar()
  plt.tight_layout()

  plt.savefig(os.path.join(output_filename, "sub001_noise_regressors_betas_map.png"), format='png', dpi=500)  

    def __call__(self, transform_xy, x1, y1, x2, y2):
        """
        get extreme values.

        x1, y1, x2, y2 in image coordinates (0-based)
        nx, ny : number of divisions in each axis
        """
        x_, y_ = np.linspace(x1, x2, self.nx), np.linspace(y1, y2, self.ny)
        x, y = np.meshgrid(x_, y_)
        lon, lat = transform_xy(np.ravel(x), np.ravel(y))

        # iron out jumps, but algorithm should be improved.
        # This is just naive way of doing and my fail for some cases.
        # Consider replacing this with numpy.unwrap
        # We are ignoring invalid warnings. They are triggered when
        # comparing arrays with NaNs using > We are already handling
        # that correctly using np.nanmin and np.nanmax
        with np.errstate(invalid='ignore'):
            if self.lon_cycle is not None:
                lon0 = np.nanmin(lon)
                lon -= 360. * ((lon - lon0) > 180.)
            if self.lat_cycle is not None:
                lat0 = np.nanmin(lat)
                lat -= 360. * ((lat - lat0) > 180.)

        lon_min, lon_max = np.nanmin(lon), np.nanmax(lon)
        lat_min, lat_max = np.nanmin(lat), np.nanmax(lat)

        lon_min, lon_max, lat_min, lat_max = \
                 self._adjust_extremes(lon_min, lon_max, lat_min, lat_max)

        return lon_min, lon_max, lat_min, lat_max

    def generate_artist(self):
           container=self.get_container()
           if self.isempty() is False:
               return

           x, y = self._eval_xy() # this handles "use_var"
 
           lp=self.getp("loaded_property") 

           if True:
              x, y = self.getp(("x", "y"))
              if y is None: return
              if x is None: return 
              
              if (x is not None  and
                  y is not None):   
                  self._data_extent=[np.nanmin(x), np.nanmax(x), 
                                     np.nanmin(y), np.nanmax(y)]

                  if len(y.shape) == 1:
                     kywds = self._var["kywds"]
                     args, self._tri =  tri_args(x, y, self._tri) 
                     kywds['mask'] = self.getp('mask')
                     kywds['linestyle'] = self.getp('linestyle')
                     kywds['linewidth'] = self.getp('linewidth')
                     kywds['color'] = self.getp('color')
                     a =  triplot(container, *args, **kywds)
                     self.set_artist(a[0])
                     self._other_artists = a[1:]

           if lp is not None:
              for i in range(0, len(lp)):
                  self.set_artist_property(self._artists[i], lp[i])
              self.delp("loaded_property")
           self.set_rasterized()

    def acquire_data(self, var_name=None, slice_=()):
        if var_name in self._variables:
            vars = [var_name]
        else:
            vars = self._variables

        if not isinstance(slice_, tuple): slice_ = (slice_,)

        for vn in vars:
            var = self._data_array[vn]

            ndims = len(var.shape)
            # Ensure the slice_ is the appropriate length
            if len(slice_) < ndims:
                slice_ += (slice(None),) * (ndims-len(slice_))

            arri = ArrayIterator(var, self._block_size)[slice_]
            for d in arri:
                if d.dtype.char is "S":
                    # Obviously, we can't get the range of values for a string data type!
                    rng = None
                elif isinstance(d, numpy.ma.masked_array):
                    # TODO: This is a temporary fix because numpy 'nanmin' and 'nanmax'
                    # are currently broken for masked_arrays:
                    # http://mail.scipy.org/pipermail/numpy-discussion/2011-July/057806.html
                    dc = d.compressed()
                    if dc.size == 0:
                        rng = None
                    else:
                        rng = (numpy.nanmin(dc), numpy.nanmax(dc))
                else:
                    rng = (numpy.nanmin(d), numpy.nanmax(d))
                yield vn, arri.curr_slice, rng, d

        return

 def test_threshold_filter_nan(self):
     src = self.make_src(nan=True)
     self.e.add_source(src)
     threshold = Threshold()
     self.e.add_filter(threshold)
     self.assertEqual(np.nanmin(src.scalar_data), np.nanmin(threshold.outputs[0].point_data.scalars.to_array()))
     self.assertEqual(np.nanmax(src.scalar_data), np.nanmax(threshold.outputs[0].point_data.scalars.to_array()))

def get_region_boxes(sp, reg2sp):
  x = np.arange(0, sp.shape[1])
  y = np.arange(0, sp.shape[0])
  xv, yv = np.meshgrid(x, y)
  maxsp = np.max(sp)
  sp1=sp.reshape(-1)-1
  xv = xv.reshape(-1)
  yv = yv.reshape(-1)
  spxmin = accum.my_accumarray(sp1,xv, maxsp, 'min')
  spymin = accum.my_accumarray(sp1,yv, maxsp, 'min')
  spxmax = accum.my_accumarray(sp1,xv, maxsp, 'max')
  spymax = accum.my_accumarray(sp1,yv, maxsp, 'max')
  
  Z = reg2sp.astype(float, copy=True)
  Z[reg2sp==0] = np.inf
  xmin = np.nanmin(np.multiply(spxmin.reshape(-1,1), Z),0)
  ymin = np.nanmin(np.multiply(spymin.reshape(-1,1), Z),0)
  xmax = np.amax(np.multiply(spxmax.reshape(-1,1), reg2sp),0)
  ymax = np.amax(np.multiply(spymax.reshape(-1,1), reg2sp), 0)
  xmin[np.isinf(xmin)]=0
  ymin[np.isinf(ymin)]=0
  

  boxes = np.hstack((xmin.reshape(-1,1), ymin.reshape(-1,1), xmax.reshape(-1,1), ymax.reshape(-1,1)))
  return boxes 

def classify(request):
    C = json.loads(request.POST["C"])
    try:
        features, labels = get_multi_features(request)
    except ValueError as e:
        return HttpResponse(json.dumps({"status": e.message}))
    try:
        kernel = get_kernel(request, features)
    except ValueError as e:
        return HttpResponse(json.dumps({"status": e.message}))
    
    learn = "No"  
    values=[]

    try:
        domain = json.loads(request.POST['axis_domain'])
        x, y, z = svm.classify_svm(sg.GMNPSVM, features, labels, kernel, domain, learn, values, C, False)
    except Exception as e:
        return HttpResponse(json.dumps({"status": repr(e)}))

#    z = z + np.random.rand(*z.shape) * 0.01
	
    z_max = np.nanmax(z)
    z_min = np.nanmin(z)
    z_delta = 0.1*(np.nanmax(z)-np.nanmin(z))
    data = {"status": "ok",
            "domain": [z_min-z_delta, z_max+z_delta],
            "max": z_max+z_delta,
            "min": z_min-z_delta,
            "z": z.tolist()}

    return HttpResponse(json.dumps(data))

    def _axes_domain(self, nx=None, ny=None, background_patch=None):
        """Returns x_range, y_range"""
        DEBUG = False

        transform = self._crs_transform()

        ax_transform = self.axes.transAxes
        desired_trans = ax_transform - transform

        nx = nx or 30
        ny = ny or 30
        x = np.linspace(1e-9, 1 - 1e-9, nx)
        y = np.linspace(1e-9, 1 - 1e-9, ny)
        x, y = np.meshgrid(x, y)

        coords = np.concatenate([x.flatten()[:, None], y.flatten()[:, None]], 1)

        in_data = desired_trans.transform(coords)

        ax_to_bkg_patch = self.axes.transAxes - background_patch.get_transform()

        ok = np.zeros(in_data.shape[:-1], dtype=np.bool)
        # XXX Vectorise contains_point
        for i, val in enumerate(in_data):
            # convert the coordinates of the data to the background
            # patches coordinates
            background_coord = ax_to_bkg_patch.transform(coords[i : i + 1, :])
            bkg_patch_contains = background_patch.get_path().contains_point
            if bkg_patch_contains(background_coord[0, :]):
                color = "r"
                ok[i] = True
            else:
                color = "b"

            if DEBUG:
                import matplotlib.pyplot as plt

                plt.plot(coords[i, 0], coords[i, 1], "o" + color, clip_on=False, transform=ax_transform)
        #                plt.text(coords[i, 0], coords[i, 1], str(val), clip_on=False,
        #                         transform=ax_transform, rotation=23,
        #                         horizontalalignment='right')

        inside = in_data[ok, :]
        x_range = np.nanmin(inside[:, 0]), np.nanmax(inside[:, 0])
        y_range = np.nanmin(inside[:, 1]), np.nanmax(inside[:, 1])

        # XXX Cartopy specific thing. Perhaps make this bit a specialisation
        # in a subclass...
        crs = self.crs
        if isinstance(crs, Projection):
            x_range = np.clip(x_range, *crs.x_limits)
            y_range = np.clip(y_range, *crs.y_limits)

            # if the limit is >90 of the full x limit, then just use the full
            # x limit (this makes circular handling better)
            prct = np.abs(np.diff(x_range) / np.diff(crs.x_limits))
            if prct > 0.9:
                x_range = crs.x_limits

        return x_range, y_range

	def show(self,**kwargs):
		display = kwargs.get('display', True)
		show_layers = kwargs.get('show_layers',self.layers)
		try:
			show_layers=sorted(show_layers)
		except TypeError:
			show_layers=[show_layers]
		extent=kwargs.get('extent', 
						max_axis(*tuple(_image.axis for _image in self.image_sorted[self.layers[0]])))
		vmin=kwargs.get('vmin')
		vmax=kwargs.get('vmax')
		fig = plt.figure(figsize=(8, 8*abs((extent[3]-extent[2])*1./(extent[1]-extent[0]))))
		for layer in show_layers:
			for image in self.image_sorted[layer]:
				if layer==show_layers[0] and image==self.image_sorted[layer][0]:
					if not vmin:
						kwargs['vmin']=np.nanmin(image.image)
						vmin=np.nanmin(image.image)
					if not vmax:
						kwargs['vmax']=np.nanmax(image.image)
						vmax=np.nanmax(image.image)
					image.show(hold=True,**kwargs)
				else:
					image.show(hold=True,vmin=vmin,vmax=vmax,scalebar='off',colorbar='off')
		plt.xlim(extent[:2])
		plt.ylim(extent[-2:])
		if display:
			plt.show()
		else:
			return fig