def lacosmiciteration(self, verbose = None):
        """
        Performs one iteration of the L.A.Cosmic algorithm.
        It operates on self.cleanarray, and afterwards updates self.mask by adding the newly detected
        cosmics to the existing self.mask. Cleaning is not made automatically ! You have to call
        clean() after each iteration.
        This way you can run it several times in a row to to L.A.Cosmic "iterations".
        See function lacosmic, that mimics the full iterative L.A.Cosmic algorithm.
        
        Returns a dict containing
            - niter : the number of cosmic pixels detected in this iteration
            - nnew : among these, how many were not yet in the mask
            - itermask : the mask of pixels detected in this iteration
            - newmask : the pixels detected that were not yet in the mask
    
        If findsatstars() was called, we exclude these regions from the search.
    
        """
        
        if verbose == None:
            verbose = self.verbose

        if verbose:
            print "Convolving image with Laplacian kernel ..."
        
        # We subsample, convolve, clip negative values, and rebin to original size
        subsam = subsample(self.cleanarray)
        conved = signal.convolve2d(subsam, laplkernel, mode="same", boundary="symm")
        cliped = conved.clip(min=0.0)
        #cliped = np.abs(conved) # unfortunately this does not work to find holes as well ...
        lplus = rebin2x2(cliped)
        
        if verbose:
            print "Creating noise model ..."
            
        # We build a custom noise map, so to compare the laplacian to
        m5 = ndimage.filters.median_filter(self.cleanarray, size=5, mode='mirror')
        # We keep this m5, as I will use it later for the interpolation.
        m5clipped = m5.clip(min=0.00001) # As we will take the sqrt
        noise = (1.0/self.gain) * np.sqrt(self.gain*m5clipped + self.readnoise*self.readnoise)

        if verbose:
           print "Calculating Laplacian signal to noise ratio ..."

        # Laplacian signal to noise ratio :
        s = lplus / (2.0 * noise) # the 2.0 is from the 2x2 subsampling
        # This s is called sigmap in the original lacosmic.cl
        
        # We remove the large structures (s prime) :
        sp = s - ndimage.filters.median_filter(s, size=5, mode='mirror')
         
        if verbose:
            print "Selecting candidate cosmic rays ..."
            
         # Candidate cosmic rays (this will include stars + HII regions)
        candidates = sp > self.sigclip    
        nbcandidates = np.sum(candidates)
        
        if verbose:
            print "  %5i candidate pixels" % nbcandidates
         
         # At this stage we use the saturated stars to mask the candidates, if available :
        if self.satstars is not None:
            if verbose:
                 print "Masking saturated stars ..."
                 candidates = np.logical_and(np.logical_not(self.satstars), candidates)
                 nbcandidates = np.sum(candidates)
        
            if verbose:
                print "  %5i candidate pixels not part of saturated stars" % nbcandidates
         
        if verbose:
            print "Building fine structure image ..."
            
         # We build the fine structure image :
        m3 = ndimage.filters.median_filter(self.cleanarray, size=3, mode='mirror')
        m37 = ndimage.filters.median_filter(m3, size=7, mode='mirror')
        f = m3 - m37
        # In the article that's it, but in lacosmic.cl f is divided by the noise...
        # Ok I understand why, it depends on if you use sp/f or L+/f as criterion.
        # There are some differences between the article and the iraf implementation.
        # So I will stick to the iraf implementation.
        f = f / noise
        f = f.clip(min=0.01) # as we will divide by f. like in the iraf version.
        
        if verbose:
            print "Removing suspected compact bright objects ..."
            
        # Now we have our better selection of cosmics :
        cosmics = np.logical_and(candidates, sp/f > self.objlim)
        # Note the sp/f and not lplus/f ... due to the f = f/noise above.
        
        nbcosmics = np.sum(cosmics)
        
        if verbose:
            print "  %5i remaining candidate pixels" % nbcosmics
        
        # What follows is a special treatment for neighbors, with more relaxed constains.
        
        if verbose:
            print "Finding neighboring pixels affected by cosmic rays ..."
            
        # We grow these cosmics a first time to determine the immediate neighborhod  :
        growcosmics = np.cast['bool'](signal.convolve2d(np.cast['float32'](cosmics), growkernel, mode="same", boundary="symm"))
        
        # From this grown set, we keep those that have sp > sigmalim
        # so obviously not requiring sp/f > objlim, otherwise it would be pointless
        growcosmics = np.logical_and(sp > self.sigclip, growcosmics)
        
        # Now we repeat this procedure, but lower the detection limit to sigmalimlow :
            
        finalsel = np.cast['bool'](signal.convolve2d(np.cast['float32'](growcosmics), growkernel, mode="same", boundary="symm"))
        finalsel = np.logical_and(sp > self.sigcliplow, finalsel)
        
        # Again, we have to kick out pixels on saturated stars :
        if self.satstars is not None:
            if verbose:
                print "Masking saturated stars ..."
            finalsel = np.logical_and(np.logical_not(self.satstars), finalsel)
             
        nbfinal = np.sum(finalsel)
        
        if verbose:
            print "  %5i pixels detected as cosmics" % nbfinal
        
        # Now the replacement of the cosmics...
        # we outsource this to the function clean(), as for some purposes the cleaning might not even be needed.
        # Easy way without masking would be :
        #self.cleanarray[finalsel] = m5[finalsel]
        
        # We find how many cosmics are not yet known :
        newmask = np.logical_and(np.logical_not(self.mask), finalsel)
        nbnew = np.sum(newmask)
        
        # We update the mask with the cosmics we have found :
        self.mask = np.logical_or(self.mask, finalsel)
    
        # We return
        # (used by function lacosmic)
        
        return {"niter":nbfinal, "nnew":nbnew, "itermask":finalsel, "newmask":newmask}