I have written a procedure to find the background radiation in an image. I use sigmacutting to remove pixels with high values - often objects; imfit = image ; this is the array used by ;sfit, to find a polynomial equation ;representing the ;background im2 = image ; Image to be returned after removing the background radiation. im3 = image ; Median-filter on im2 to remove bad pixels. im2=median(im2,4) ; Same on imfit imfit=median(imfit,5) Here comes the problem: res0 = sfit(imfit, 0) ; First polynomial fit, of degree 0. sigma0 = stddev(imfit - res0) ; Standard deviation of the residual. im3 = im2 - res0 ; Image with background corresponding to a polynomial fit of degree 0, removed. ; Removes high and low values in imfit imfit[where(imfit gt (res0+3*sigma0) or imfit lt (res0-3*sigma0))]=res0[where(imfit gt (res0+3*sigma0) or imfit lt (res0-3*sigma0))] This has worked fine before, but I have dumped into an image where it does not work. I get the error message: % Attempt to subscript RES0 with is out of range. The thing is that it usually works, but not here. I can't find any differences between the other images either. Is the where-sentence wrong? Typing; help, imfit help, res0 help, sigma0 Gives: IMFIT FLOAT = Array[128, 96] RES0 DOUBLE = Array[128, 96] SIGMA0 DOUBLE = 0.30962909 Furthermore I do like this: degree = 1 res1 = sfit(imfit, degree) ; Second polynomial fit, of degree 1. sigma1 = stddev(imfit - res1) ; Standard deviation of the residual. ; Minimizing the standard deviation of the residual. ; Checks if sigma1 is smaller than sigma0. Finds a new image with ; background, corresponding to a polynomial fit of degree +1, removed ; if sigma1 < sigma0. i=1 while sigma1 lt sigma0 do begin res0=res1 ; Previous polynomial fit. sigma0 = sigma1 ; Previous standard deviation of the residual. im3 = im2 - res0 ; New image with previous polynomial fit removed. i=i+1 tv,res0&wait,5 help, imfit help, res0 help, sigma0 ; Removes high and low values in imfit imfit[where(imfit gt (res0+3*sigma0) or imfit lt (res0-3*sigma0))]=res0[where(imfit gt (res0+3*sigma0) or imfit lt (res0-3*sigma0))] print, res0[where(imfit gt (res0+3*sigma0) or imfit lt (res0-3*sigma0))] ; - actually this always gives the same error output... degree = degree + 1 ; New degree of the polynomial fit. res1 = sfit(imfit, degree) ; New polynomial fit. sigma1 = stddev(imfit - res1) ; New standard deviation of the residual. endwhile Morten