def clean(self, mask = None, verbose = None):
        """
        Given the mask, we replace the actual problematic pixels with the masked 5x5 median value.
        This mimics what is done in L.A.Cosmic, but it's a bit harder to do in python, as there is no
        readymade masked median. So for now we do a loop...
        Saturated stars, if calculated, are also masked : they are not "cleaned", but their pixels are not
        used for the interpolation.
        
        We will directly change self.cleanimage. Instead of using the self.mask, you can supply your
        own mask as argument. This might be useful to apply this cleaning function iteratively.
        But for the true L.A.Cosmic, we don't use this, i.e. we use the full mask at each iteration.

        """
        if verbose == None:
            verbose = self.verbose
        if mask == None:
            mask = self.mask
            
        if verbose:
            print "Cleaning cosmic affected pixels ..."
        
        # So... mask is a 2D array containing False and True, where True means "here is a cosmic"
        # We want to loop through these cosmics one by one.
        cosmicindices = np.argwhere(mask)
        # This is a list of the indices of cosmic affected pixels.
        #print cosmicindices
        
        # We put cosmic ray pixels to np.Inf to flag them :
        self.cleanarray[mask] = np.Inf
        
        # Now we want to have a 2 pixel frame of Inf padding around our image.
        w = self.cleanarray.shape[0]
        h = self.cleanarray.shape[1]
        padarray = np.zeros((w+4,h+4))+np.Inf
        padarray[2:w+2,2:h+2] = self.cleanarray.copy() # that copy is important, we need 2 independent arrays
        
        # The medians will be evaluated in this padarray, skipping the np.Inf.
        # Now in this copy called padarray, we also put the saturated stars to np.Inf, if available :
        if self.satstars is not None:
            padarray[2:w+2,2:h+2][self.satstars] = np.Inf
            # Viva python, I tested this one, it works...
        
        # A loop through every cosmic pixel :
        for cosmicpos in cosmicindices:
            x = cosmicpos[0]
            y = cosmicpos[1]
            cutout = padarray[x:x+5, y:y+5].ravel() # remember the shift due to the padding !
            #print cutout
            # Now we have our 25 pixels, some of them are np.Inf, and we want to take the median
            goodcutout = cutout[cutout != np.Inf]
            #print np.alen(goodcutout)
            
            if np.alen(goodcutout) >= 25 :
                # This never happened, but you never know ...
                raise RuntimeError, "Mega error in clean !"
            elif np.alen(goodcutout) > 0 :
                replacementvalue = np.median(goodcutout)
            else :    
                # i.e. no good pixels : Shit, a huge cosmic, we will have to improvise ...
                print "OH NO, I HAVE A HUUUUUUUGE COSMIC !!!!!"
                replacementvalue = self.guessbackgroundlevel()
            
            # We update the cleanarray,
            # but measure the medians in the padarray, so to not mix things up...
            self.cleanarray[x, y] = replacementvalue
            
        # That's it.
        if verbose:
            print "Cleaning done"

        # FYI, that's how the LACosmic cleaning looks in iraf :
        """
        imarith(outmask,"+",finalsel,outmask)
        imreplace(outmask,1,lower=1,upper=INDEF) # ok so outmask = 1 are the cosmics
        imcalc(outmask,inputmask,"(1.-10000.*im1)",verb-)
        imarith(oldoutput,"*",inputmask,inputmask)
        median(inputmask,med5,5,5,zloreject=-9999,zhi=INDEF,verb-)
        imarith(outmask,"*",med5,med5)
        if (i>1) imdel(output)
        imcalc(oldoutput//","//outmask//","//med5,output,"(1.-im2)*im1+im3",verb-)
        
        # =
        
        merging to full mask
        inputmask = 1.0 - 10000.0 * finalsel # So this is 1.0, but cosmics are very negative
        inputmask = oldoutput * inputmask # orig image, with very negative cosmics
        med5 = median of inputmask, but rejecting these negative cosmics
        # i dunno how to do this in python -> had to do the loop
        med5 = finalsel * med5 # we keep only the cosmics of this median
        # actual replacement :
        output = (1.0 - outmask)*oldoutput + med5 # ok 
        """
Given the mask, we replace the actual problematic pixels with the masked 5x5 median value.
        This mimics what is done in L.A.Cosmic, but it's a bit harder to do in python, as there is no
        readymade masked median. So for now we do a loop...
        Saturated stars, if calculated, are also masked : they are not "cleaned", but their pixels are not
        used for the interpolation.
        
        We will directly change self.cleanimage. Instead of using the self.mask, you can supply your
        own mask as argument. This might be useful to apply this cleaning function iteratively.
        But for the true L.A.Cosmic, we don't use this, i.e. we use the full mask at each iteration.