;Apply threshold to image img; use the list of non-masked pixels
;(according to mask): non_masked_img, and update it by removing the 
;pixels whose value are above the threshold
function apply_img_threshold, img, non_masked_img, threshold
  above = where(img GT threshold, count)
  if(count ne 0) then begin
                                ;get pixels subscripts that are below
                                ;threshold but are not already masked
     diff = setdifference(non_masked_img, above)
     if(size(diff, /dimensions) GT 1) then begin
        res = diff
     endif else begin
        res = non_masked_img
     endelse
  endif
  return, res
end