def setzscale(self, z1="auto", z2="auto", nsig=3, samplesizelimit = 10000, border=300):
        """
        We set z1 and z2, according to different algorithms or arguments.
        
        For both z1 and z2, give either :
        
            - "auto" (default automatic, different between z1 and z2)
            - "ex" (extrema)
            - "flat" ("sigma-cuts" around median value, well-suited for flatfields)
            - numeric value like 1230.34
        
        nsig is the number of sigmas to be rejected (used by auto z1 + both flats)
        
        samplesizelimit is the maximum number of pixels to compute statistics on.
        If your image is larger then samplesizelimit, I will use only samplesizelimit pixels of it.
        
        If your image is 3 times border in width and height, I will skip border pixels around the image before
        doing calculations. This is made to get rid of the overscan and prescan etc.
        So you can basically leave this at 300, it will only affect images wider then 900 pixels.
        (300 happens to be a safe value for many telescopes.)
        You can put border = 0 to deactivate this feature.
        
        If you give nothing, the cutoff will not be changed.
        You should set the z scale directly after cropping the image.
        
        """
        
        if self.pilimage != None:
            raise RuntimeError, "Cannot set z scale anymore, PIL image already exists !"
        
        if self.numpyarray.shape[0] > 3 * border and self.numpyarray.shape[1] > 3 * border:
            if border > 0:
                if self.verbose :
                    print "For the stats I will leave a border of %i pixels" % border
                calcarray = self.numpyarray[border:-border, border:-border].copy()
            else:
                calcarray = self.numpyarray.copy()
        else:
            calcarray = self.numpyarray.copy()
            if self.verbose:
                print "Image is too small for a border of %i" % (border)
        
        
        # Starting with the simple possibilities :
        if z1 == "ex" :
            self.z1 = np.min(calcarray)
            if self.verbose:
                print "Setting ex z1 to %f" % self.z1
            
        if z2 == "ex":
            self.z2 = np.max(calcarray)
            if self.verbose:
                print "Setting ex z2 to %f" % self.z2
            
        if type(z1) == type(0) or type(z1) == type(0.0):
            self.z1 = z1
            if self.verbose:
                print "Setting z1 to %f" % self.z1
            
        if type(z2) == type(0) or type(z2) == type(0.0):
            self.z2 = z2
            if self.verbose:
                print "Setting z2 to %f" % self.z2
        
        # Now it gets a little more sophisticated.
        if z1 == "auto" or z2 == "auto" or z1 == "flat" or z2 == "flat":
            
            # To speed up, we do not want to do statistics on the full image if it is large.
            # So we prepare a small random sample of pixels.
            
            calcarray.shape = calcarray.size # We flatten the 2D array
            if calcarray.size > samplesizelimit :
                #selectionindices = np.random.random_integers(low = 0, high = calcarray.size - 1, size=samplesizelimit)
                selectionindices = np.linspace(0, calcarray.size-1, samplesizelimit).astype(np.int)
                statsel = calcarray[selectionindices]
            else :
                statsel = calcarray
            
                        
            #nbrofbins = 10 + int(np.log10(calcarray.size)*10.0)
            #print "Building histogram with %i bins" % nbrofbins
            #nbrofbins = 100
            #hist = np.histogram(statsel, bins=nbrofbins, range=(self.z1, self.z2), normed=False, weights=None, new=True)
            
            medianlevel = np.median(statsel)
            firststd = np.std(statsel)
            
            if z1 == "auto" :
            
                # 2 sigma clipping (quick and dirty star removal) :
                nearskypixvals = statsel[np.logical_and(statsel > medianlevel - 2*firststd, statsel < medianlevel + 2*firststd)]
                
                skylevel = np.median(nearskypixvals)
                secondstd = np.std(nearskypixvals)
                
                if self.verbose :
                    print "Sky level at %f +/- %f" % (skylevel, secondstd)
                
                self.z1 = skylevel - nsig*secondstd
                if self.verbose :
                    print "Setting auto z1 to %f, nsig = %i" % (self.z1, nsig)
                                
                
                
            if z2 == "auto" :
                # Here we want to reject a percentage of high values...
                sortedstatsel = np.sort(statsel) 
                n = round(0.9995 * statsel.size)
                self.z2 = sortedstatsel[n]
                if self.verbose :
                    print "Setting auto z2 to %f" % self.z2 
                
            if z1 == "flat" :
                
                # 5 sigma clipping to get rid of cosmics :
                nearflatpixvals = statsel[np.logical_and(statsel > medianlevel - 5*firststd, statsel < medianlevel + 5*firststd)]
                
                flatlevel = np.median(nearflatpixvals)
                flatstd = np.std(nearflatpixvals)
    
        
                self.z1 = flatlevel - nsig*flatstd
    
                if self.verbose :
                    print "Setting flat z1 : %f, nsig = %i" % (self.z1, nsig)

            if z2 == "flat" : # symmetric to z1
                
                # 5 sigma clipping to get rid of cosmics :
                nearflatpixvals = statsel[np.logical_and(statsel > medianlevel - 5*firststd, statsel < medianlevel + 5*firststd)]
                
                flatlevel = np.median(nearflatpixvals)
                flatstd = np.std(nearflatpixvals)
    
        
                self.z2 = flatlevel + nsig*flatstd
    
                if self.verbose :
                    print "Setting flat z2 : %f, nsig = %i" % (self.z2, nsig)