improved bad pixel correction, added reading of bad_pixel file; added output...

improved bad pixel correction, added reading of bad_pixel file; added output control (verbose level); apply_calib_on_raw now checks for inf in gain and offset;small bugfixes

IR_config_constants.py

@@ -10,6 +10,7 @@ archivepath="http://archive-webapi.ipp-hgw.mpg.de/Test/raw/W7X/"
 calibpath="\\\\sv-e4-fs-1\\E4-Mitarbeiter\\E4 Diagnostics\\QIR\\Calibrations\\"
 IRCamColdframes_fittingpath=calibpath+"IRCamColdframes_fitting\\"
 IRCamRefImagespath = calibpath+'IRCamReferenceImages\\'
+IRCAMBadPixels_path = calibpath + 'IRCAMBadPixels\\'
 
 #%% Dictionaries
 portcamdict={

IR_image_tools.py

@@ -14,12 +14,54 @@ Tools for:
 """
 import numpy as np
 import matplotlib.pyplot as plt
-from IR_config_constants import portcamdict,IRCamRefImagespath
+from IR_config_constants import portcamdict,IRCamRefImagespath,IRCAMBadPixels_path
 import h5py
 from os.path import join, basename
 import glob
+import datetime
 
 
+def get_OP_by_time(time_ns=None, shot_no=None, program_str=None):
+    '''Derives operation phase (OP) of W7-X based on either:
+       a nanosacond time stamp, a MDSplus style shot no. or a program ID.
+       IN:
+          time_ns      - integer of nanosecond time stamp, 
+                         e.g. 1511972727249834301 (OPTIONAL)
+          shot_no      - integer of MDSplus style shot number, 
+                         e.g. 171207022 (OPTIONAL)
+          program_str  - string of CoDaQ ArchiveDB style prgram number or date, 
+                         e.g. '20171207.022' or '20171207' (OPTIONAL)
+       RETURN:
+          conn         - MDSplus connection object, to be used in e.g. 1511972727249834301
+                         read_MDSplus_image_simple(), read_MDSplus_metadata()
+   '''
+    # derive operation phase (OP) from time as nanosecond time stamp or string
+    if time_ns is not None:
+        dateOP = datetime.datetime.utcfromtimestamp(time_ns/1e9)
+    elif shot_no is not None:
+        dateOP = datetime.datetime.strptime(str(shot_no)[:6], '%y%m%d')
+    elif program_str is not None:
+        dateOP = datetime.datetime.strptime(program_str[:8], '%Y%m%d')
+    else:
+        raise Exception('get_OP_by_time: ERROR! neither time, shot no. or program ID provided')
+        
+    if dateOP.year == 2017:
+        if dateOP.month>8 and dateOP.month<12:
+            return "OP1.2a"
+        elif dateOP.month==8 and dateOP.day>28:
+            return "OP1.2a"
+        elif dateOP.month==12 and dateOP.day<8:
+            return "OP1.2a"
+        else:
+            return None        
+    elif dateOP.year == 2018:
+        return "OP1.2b"
+    elif dateOP.year <= 2016 and dateOP.year >= 2015:
+        if (dateOP.year == 2016 and dateOP.month <= 3) or (dateOP.year == 2015 and dateOP.month == 12):
+            return "OP1.1"
+        else:
+            return None
+
 
 def bestimmtheitsmaß_general(data,fit):
     R=0
@@ -105,6 +147,32 @@ def check_backgroundframe(backgroundframe,threshold=50):
         valid=False    
     return valid,np.mean(dataset)
     
+def read_bad_pixels_from_file(port, shot_no=None, program=None):
+    '''Reads bad pixels stored in *.bpx file on E4 server.
+       Requires one of the optional arguments shot_no or program.
+        IN
+            port            - integer of port no of camera
+            shot_no         - integer of MDSplus style shot number, e.g. 171207022 (OPTIONAL)
+            program         - string of CoDaQ ArchiveDB style prgram number or date, 
+                              e.g. '20171207.022' or '20171207' (OPTIONAL)
+        OUT
+            bad_pixle_list  - list of tuples (row,column) of pixel coordinates 
+                              as integer
+    '''
+    if shot_no is not None:    
+        OP = get_OP_by_time(shot_no=shot_no)
+    elif program is not None:
+        OP = get_OP_by_time(program_str=program)
+    else:
+        raise Exception('read_bad_pixels_from_file: ERROR! Need either shot no. or program string.')
+        
+    port_name = 'AEF{0}'.format(port)
+    bad_pixel_file = 'badpixel_{0}.bpx'.format(portcamdict[OP][port_name][6:])
+    data = np.genfromtxt(IRCAMBadPixels_path+bad_pixel_file, dtype=int)
+    bad_pixle_list = list(zip(data[:,1], data[:,0]))
+    return bad_pixle_list
+    
+
 def find_outlier_pixels(frame,tolerance=3,worry_about_edges=True,plot_it=False):
     # This function finds the bad pixels in a 2D dataset. 
     # Tolerance is the number of standard deviations used for cutoff.
@@ -116,7 +184,6 @@ def find_outlier_pixels(frame,tolerance=3,worry_about_edges=True,plot_it=False):
     mean = np.mean(difference)
     if plot_it:
         
-        
         fig = plt.figure()
         fig.suptitle('find_outlier_pixels: histogram')
         plt.hist(difference.ravel(),50,log=True,histtype='stepfilled')
@@ -144,108 +211,190 @@ def find_outlier_pixels(frame,tolerance=3,worry_about_edges=True,plot_it=False):
    
     return bad_pixels
 
+
 def correct_images(images,badpixels):
+    print('correct_images: New routine restore_bad_pixels() is used and can be called directly. Check out "help(restore_bad_pixels)"')
     if type(badpixels)!=int:
-        for i in range(len(images)):
-            images[i]=(restore_pixels(images[i],np.invert(badpixels==1))).astype(np.float32)
+        if type(images) == list:
+            # return corrected images also as list of 2D arrays
+            images = restore_bad_pixels(images, np.invert(badpixels==1)).astype(np.float32)
+            images = list(images)
+        else:
+            # keep shape
+            images = restore_bad_pixels(images, np.invert(badpixels==1)).astype(np.float32)
+        
+#        for i in range(len(images)):
+#            images[i]=(restore_pixels(images[i],np.invert(badpixels==1))).astype(np.float32)
         print("done")
     return images
 
-def restore_pixels(frame, bad_pixel):# code from Peter from JET
+
+def restore_bad_pixels(frames, bad_pixel, by_list=True, check_neighbours=True, plot_it=False, verbose=0):
+    """Restore bad pixel by interpolation of adjacent pixels. Optionally make 
+       sure that adjacent pixels are not bad (time consuming). Default is to use 
+       a list of bad pixels and a for loop. For many bad pixels consider using 
+       the optinal alternative using a bad pixel mask.
+        IN:
+            frames              - either list of frames as 2D numpy array, 
+                                  or 3D numpy array (frame number, n_rows, n_cols),
+                                  or 2D numpy array (n_rows, n_cols)
+            bad_pixel           - either list of tuples of bad pixel coordinates,
+                                  or mask of pixel status (good=True, bad=False)
+            by_list             - boolean of whether to use a list and a for loop (True),
+                                  or to use a mask of bad pixel and array operations (False)
+                                  (OPTIONAL: if not provided, True (list) is default)
+            check_neighbours    - boolean of whether to check if neighbours of a bad pixel
+                                  are not bad either before computing a mean
+                                  (works only in list mode!)
+                                  (OPTIONAL: if not provided, check is on)
+            plot_it             - boolean to decide whether to plot intermediate
+                                  results or not
+                                  (OPTIONAL: if not provided, switched off)
+            verbose             - integer of feedback level (amount of prints)
+                                  (OPTIONAL: if not provided, only ERROR output)
+        RETURN:
+            frames              - 3D numpy array (frame number, n_rows, n_cols) of 
+                                  corrected frames
+    """
+
+    # make sure frames is correctly shaped
+    if type(frames) == list:
+        frames = np.array(frames)
+        frame_shape = 'list'
+    else:
+        if len(np.shape(frames)) == 2:
+            frames = np.array([frames])
+            frame_shape = '2D'
+        elif len(np.shape(frames)) == 3:
+            frame_shape = '3D'
+            pass
+        else:
+            raise Exception('restore_bad_pixels: ERROR! Unexpected shape of frames.')
+    frame_dtype = frames.dtype
+#    frames = frames.astype(float)
+    n_frames, n_rows, n_cols = np.shape(frames)
+    if plot_it:
+        start_frame = np.copy(frames[0])
+    
     # make sure bad pixel are provided as mask and list  
     if type(bad_pixel) is list:
         blist = bad_pixel
-        bmask = np.ones(frame.shape,dtype=bool)
-        for pix in bad_pixel:
+        bmask = np.ones([n_rows, n_cols],dtype=bool)
+        for pix in blist:
             bmask[pix] = False
         bmask = np.invert(bmask)
     else:
-        bmask = np.invert(bad_pixel)
-        x,y = np.where(bmask)
-        blist = list(zip(x,y))
-#    print('restore_pixels() working on: ',blist)
-    # prepare internal result frames
-    resframe = frame.astype(float)
-    # -----------------------------------
-    # restore by list
-    #
-    height = np.shape(frame)[0]
-    width  = np.shape(frame)[1]
-    for pos in blist:
-        # assume four neighbouring, non-bad pixels exist
-        n_neighbours = 4
-        # -- top neighbour --
-        if pos[0] == 0:
-            # out of bounds
-            top = 0
-            n_neighbours -= 1
-        else:
-            i_t = pos[0] - 1
-            while (bmask[i_t,pos[1]] & (i_t >= 0)):
-                i_t -= 1
-            if i_t >= 0:
-                top = frame[i_t,pos[1]]
-            else:
-                # out of bounds
-                top = 0
-                n_neighbours -= 1
-            #print('top: ',i_t, top)    
-        # -- bottom neighbour --
-        if pos[0] == height-1:
-            # out of bounds
-            bottom = 0
-            n_neighbours -= 1
+        if np.shape(bad_pixel)[0] == n_rows and np.shape(bad_pixel)[1] == n_cols:
+            bmask = np.invert(bad_pixel)            
+            x,y = np.where(bmask)
+            blist = list(zip(x,y))
         else:
-            i_b = pos[0] + 1
-            while (bmask[i_b,pos[1]] & (i_b <= height-2)):
-                i_b += 1
-            if i_b >= 0:
-                bottom = frame[i_b,pos[1]]
-            else:
-                # out of bounds
-                bottom = 0
-                n_neighbours -= 1
-            #print('bottom: ',i_b, bottom)
-        # -- left neighbour --
-        if pos[1] == 0:
-            # out of bounds
-            left = 0
-            n_neighbours -= 1
-        else:
-            i_l = pos[1] - 1
-            while (bmask[pos[0],i_l] & (i_l >= 0)):
-                i_l -= 1
-            if i_l >= 0:
-                left = frame[pos[0],i_l]
-            else:
-                # out of bounds
-                left = 0
-                n_neighbours -= 1
-            #print('left: ',i_l, left)    
-        # -- right neighbour --
-        if pos[1] == width-1:
-            # out of bounds
-            right = 0
-            n_neighbours -= 1
-        else:
-            i_r = pos[1] + 1
-            while (bmask[pos[0],i_r] & (i_r <= width-2)):
-                i_r += 1
-            if i_r >= 0:
-                right = frame[pos[0],i_r]
+            raise Exception('restore_bad_pixels: ERROR! bad_pixel in bad shape {0}'.format(np.shape(bad_pixel)))
+            
+    if verbose > 0:
+        print('restore_bad_pixels: {0} bad pixels to be restored: {1} ... '.format(len(blist), blist[:3]))    
+    
+    # expand frame by rows and columns of zeros to simplify treatment of edges
+    frames = np.dstack([np.zeros([n_frames,n_rows], dtype=frame_dtype), frames, np.zeros([n_frames,n_rows], dtype=frame_dtype)])
+    frames = np.hstack([np.zeros([n_frames,1,n_cols+2], dtype=frame_dtype), frames, np.zeros([n_frames,1,n_cols+2], dtype=frame_dtype)])
+    bmask = np.vstack([np.zeros([1,n_cols], dtype=bool), bmask, np.zeros([1,n_cols], dtype=bool)])
+    bmask = np.hstack([np.zeros([n_rows+2,1], dtype=bool), bmask, np.zeros([n_rows+2,1], dtype=bool)])
+    
+    # define number of neighbours (up to 4) ina an array of expanded frame size
+    n_neighbours = np.ones([n_frames, n_rows+2, n_cols+2])*4
+    n_neighbours[:,1,:] = 3
+    n_neighbours[:,-2,:] = 3
+    n_neighbours[:,:,1] = 3
+    n_neighbours[:,:,-2] = 3
+    n_neighbours[:,1,1] = 2
+    n_neighbours[:,1,-2] = 2
+    n_neighbours[:,-2,1] = 2
+    n_neighbours[:,-2,-2] = 2
+    
+    if by_list:
+        # ===== correct bad pixels using the list of bad pixels =====
+        #
+        
+        for pos in blist:
+            # Note:
+            # pos points to real frame coordinates, while bmask, n_neighbours have been expanded!
+            
+            if check_neighbours:
+                # takes only neighbours that are not bad
+                pos_l = np.where(bmask[pos[0]+1,:pos[1]+1]==False)[0]
+                if len(pos_l) != 0:
+                    pos_l = pos_l[-1]
+                else: 
+                    pos_l = pos[1]+1
+                pos_r = np.where(bmask[pos[0]+1,pos[1]+1:]==False)[0]
+                if len(pos_r) != 0:
+                    pos_r = pos_r[0] + pos[1]+1
+                else: 
+                    pos_r = pos[1]+2
+                pos_t = np.where(bmask[:pos[0]+1,pos[1]+1]==False)[0]
+                if len(pos_t) != 0:
+                    pos_t = pos_t[-1]
+                else: 
+                    pos_t = pos[0]+1
+                pos_b = np.where(bmask[pos[0]+1:,pos[1]+1]==False)[0]
+                if len(pos_b) != 0:
+                    pos_b = pos_b[0] + pos[0]+1
+                else: 
+                    pos_b = pos[0]+2
             else:
-                # out of bounds
-                right = 0
-                n_neighbours -= 1
-            #print('right: ',i_r, right)                
-        # averaging
-        if n_neighbours > 0:
-            #print('original value: ',frame[pos[0],pos[1]])
-            resframe[pos[0],pos[1]] = (top + bottom + left + right)/n_neighbours
-            #print('average of ',n_neighbours,' neighbours: ',frame1[pos[0],pos[1]])
-        else:
-            print('ERROR: no adjacent pixel found!')
-    return resframe
+                # insensitive to neighbours being bad as well!
+                pos_l = pos[1]
+                pos_r = pos[1]+2
+                pos_t = pos[0]
+                pos_b = pos[0]+2
+            average = (frames[:,pos[0]+1,pos_l].astype(float) + 
+                       frames[:,pos[0]+1,pos_r].astype(float) + 
+                       frames[:,pos_t,pos[1]+1].astype(float) + 
+                       frames[:,pos_b,pos[1]+1].astype(float)) / n_neighbours[:,pos[0]+1,pos[1]+1]
+            frames[:,pos[0]+1,pos[1]+1] = average.astype(frame_dtype)
+        frames = frames[:,1:-1,1:-1]
+        
+    else:
+        # ======= correct bad pixels using the bad pixel mask =======
+        #
+        # (insensitive to neighbours being bad as well!)
+   
+        # prepare mask arrays for neighbours by shifting it to left, right, top and bottom
+        bmask_l = np.hstack([bmask[:,1:], np.zeros([n_rows+2,1], dtype=bool)])
+        bmask_r = np.hstack([np.zeros([n_rows+2,1], dtype=bool), bmask[:,:-1]])
+        bmask_t = np.vstack([bmask[1:,:], np.zeros([1,n_cols+2], dtype=bool)])
+        bmask_b = np.vstack([np.zeros([1,n_cols+2], dtype=bool), bmask[:-1,:]])
+        
+    
+        # -----------------------------------
+        # restore by mask
+        #
+        frames[:,bmask] = ( (frames[:,bmask_l].astype(float) + 
+                             frames[:,bmask_r].astype(float) + 
+                             frames[:,bmask_t].astype(float) + 
+                             frames[:,bmask_b].astype(float)) / n_neighbours[:,bmask] ).astype(frame_dtype)
+        frames = frames[:,1:-1,1:-1]
+    
+    # plot comparison
+    if plot_it:
+        plt.figure()
+        plt.title('bad pixel correction of first frame')
+        m = np.mean(start_frame)
+        s = np.std(start_frame)
+        plt.imshow(start_frame, vmin=m-s, vmax=m+s)
+        plt.colorbar()
+        x,y = zip(*blist)
+        plt.scatter(y,x, marker='o', s=5, c='r', linewidths=1)
+        plt.tight_layout()
+        plt.show()
+
+    if frame_shape == 'list':
+        frames = list(frames)
+    elif frame_shape == '2D' and len(np.shape(frames))==3:
+        frames = frames[0]
+        
+    return frames
+
 
 def generate_new_hot_image(cold,reference_cold,reference_hot):
     if cold==None or reference_cold==None or reference_hot==None:

downloadversionIRdata.py

@@ -142,7 +142,7 @@ def read_restdb_old(request_url):
         t=np.array(signal_list['dimensions'])
         return True, t, signal0
 
-def download_LUT(port,time,exposure=0,emissivity=0,camera_filter=0,version=0):
+def download_LUT(port,time,exposure=0,emissivity=0,camera_filter=0,version=0, verbose=0):
     """
     download_LUT(camera,port,time,exposure=0,emissivity=0,camera_filter=0,version=1):
         time in ns
@@ -165,7 +165,8 @@ def download_LUT(port,time,exposure=0,emissivity=0,camera_filter=0,version=0):
         raise Exception
     if version==0:
         version=get_latest_version(stream+"DATASTREAM")
-        print("LUT V"+str(version)+" is used")
+        if verbose>0:
+            print("LUT V"+str(version)+" is used")
     #time=int(fu.TimeToNs([2017,9,26],[8,0,0,0]))
     LUTpar=read_restdb_old(larchivepath+"PARLOG/V"+str(version)+"/_signal.json?from="+str(time-10)+"&upto="+str(time+20))
     if LUTpar[0]:
@@ -177,12 +178,12 @@ def download_LUT(port,time,exposure=0,emissivity=0,camera_filter=0,version=0):
             del LUTpar, LUTs
             return True,LUT
         else:
-            print("unable to download the LUTs")
+            print("Warning: unable to download the LUTs")
             del LUTpar, LUTs
             return False,0
     else:
         del LUTpar
-        print("unable to find LUTs, check your request")
+        print("Warning: unable to find LUTs, check your request")
         return False,0
 
 def download_NUC_by_program(port,program,exposure,version=0):
@@ -373,7 +374,7 @@ def download_raw_images_by_times(port,starttime,stoptime,version=0,intervalSize=
             return False, 0,-1
 
 if fastDL:
-    def download_raw_images_by_program_via_png(port,program,time_s=0,version=0,threads=1):
+    def download_raw_images_by_program_via_png(port,program,time_s=0,version=0,threads=1, verbose=0):
 #        prog=AKF_1.get_program_from_PID(program)
         try:
             t_program = AKF_2.get_program_from_to(program)
@@ -439,7 +440,8 @@ if fastDL:
                 jobs = []
                 out_q=multiprocessing.Queue()
                 for i in range(threads):
-                    print("Start Thread ",i+1)                
+                    if verbose>0:
+                        print("Start Thread ",i+1)                
                     p = multiprocessing.Process(target=download_raw_images_png_by_times_thread, args=(port,tim[intervalls[i]:intervalls[i+1]],out_q,i,version,))
                     jobs.append(p)
                     p.start()
@@ -448,7 +450,8 @@ if fastDL:
                     resultdict.append(out_q.get())
                 for p in jobs:
                     p.join()
-                print("all threads are done")
+                if verbose>0:
+                    print("all threads are done")
                 order=[]
                 for ele in resultdict:
                     order.append(ele[0])
@@ -918,7 +921,7 @@ def get_temp_from_raw_by_program_V1(portnr,program,time_s=0,emi=0.8,divertorpart
 #    prog=AKF_1.get_program_from_PID(program)
     try:
         t_program = AKF_2.get_program_from_to(program)
-        prog =AKF_2.get_program_list(t_program[0], t_program[1])
+        prog = AKF_2.get_program_list(t_program[0], t_program[1])
 #    if prog[0]:
 #        t0=prog[1]['trigger']['0'][0]
         t1=prog[0]['trigger']['1'][0]
@@ -1253,18 +1256,26 @@ def get_nuced_raw_by_program_fullthreads(portnr,program,time_s=0,emi=0.8,T_versi
         else:
             return success,np.array(times),images,valid        
 
-def apply_calib_on_raw(images,background,LUT,refT=28.5,gain=0,offset=0,gain_error=0,offset_error=0,fullbackground=False,give_ERROR=False):
+def apply_calib_on_raw(images,background,LUT,refT=28.5,gain=0,offset=0,gain_error=0,offset_error=0,fullbackground=False,give_ERROR=False, verbose=0):
     try:
         #images=np.array(raw,dtype=np.uint16)
 #        del raw
 #        images=images.swapaxes(1,2)
         if type(gain)!=int and type(offset)!=int:
-            print(datetime.datetime.now(),"NUCing")
+            if verbose>0:
+                print(datetime.datetime.now(),"NUCing")
+            # eliminate bad offset and gain points
+            offset[offset==np.inf] = 0
+            offset[offset==-np.inf] = 0
+            gain[gain==np.inf] = 1
+            gain[gain==-np.inf] = 1
+            # try to apply NUC
             if give_ERROR:
                 images,error_images=apply_NUC(images,gain,offset,gain_error,offset_error,give_ERROR)
             else:
                 images=apply_NUC(images,gain,offset)
-        print(datetime.datetime.now(),"background treatment")
+        if verbose>0:
+            print(datetime.datetime.now(),"background treatment")
 #        if fullbackground:
         
         for i in range(len(images)):
@@ -1274,7 +1285,8 @@ def apply_calib_on_raw(images,background,LUT,refT=28.5,gain=0,offset=0,gain_erro
 #                images[i]=images[i]-background
 #        images=np.array(images.clip(min=0),dtype=np.uint16)
         images=[np.round(im.clip(min=0)).astype(np.uint16) for im in images]
-        print(datetime.datetime.now(),"applying LUT")
+        if verbose>0:
+            print(datetime.datetime.now(),"applying LUT")
         LUT=np.array([LUT[1],LUT[2]])
         if give_ERROR:
             terror=[]#np.zeros(np.shape(images))
@@ -1319,7 +1331,7 @@ def apply_NUC(images,gain,offset,gain_error=0,offset_error=0,give_error=False):
         print(E)
         return 0
         
-def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivity=0.8):
+def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivity=0.8, verbose=0):
 #    prog=AKF_1.get_program_from_PID(program)
     try:
         t_program = AKF_2.get_program_from_to(program)
@@ -1344,7 +1356,7 @@ def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivit
             t_exp=int(expo_DL[2][0])
             del expo_DL
             time=int(TimeToNs([2017,9,26],[8,0,0,0]))   
-            LUT_DL=download_LUT(port,time,t_exp,emissivity,cfilter,version)
+            LUT_DL=download_LUT(port,time,t_exp,emissivity,cfilter,version, verbose=verbose-1)
             if LUT_DL[0]:
                 LUT=LUT_DL[1]
                 del LUT_DL
@@ -1352,6 +1364,8 @@ def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivit
                 if back_DL[0]:
                     background=back_DL[2]
                     if Temp_V==1:
+                        if verbose>0:
+                            print('use temperature calibration version 1')
                         backtime=back_DL[1]
                         backtime=backtime.tolist()
                         divertorpart="all"   
@@ -1365,14 +1379,17 @@ def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivit
                             else:
                                 raise Exception("Unable to find thermocouple data")
                     elif Temp_V==2:
+                        if verbose>0:
+                            print('use temperature calibration version 2')
                         frame=background.copy()
                         background=get_average_background_recangle(port,background)
                         refT=28.5
                     elif Temp_V==3:
+                        if verbose>0:
+                            print('use temperature calibration version 3')
                         frame=background.copy()
                         background=get_average_background_recangle(port,background)
                         refT=28.5
-                        print("test of V3")
 #                        back_off=estimate_offset(port,program)                      
 #                        background=(background-back_off)/(back_emissivity)+back_off
                     else:
@@ -1432,7 +1449,12 @@ def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivit
                                 offset_error=np.array(NUC_DL[1][5])
                                 badpixels=np.array(NUC_DL[1][3],dtype=np.ubyte)
                                 if np.max(badpixels)==0:
-                                    badpixels=find_badpixels(port,gain,offset,niterations=10,tolerance=10)
+                                    if verbose>0:
+                                        print(datetime.datetime.now(),"Scanning for bad pixel")
+                                    # initial list from config files
+                                    init_bp_list = IR_tools.read_bad_pixels_from_file(port, program=program)
+                                    # find more bad pixel
+                                    badpixels = find_badpixels(port, gain, offset, init_bp_list=init_bp_list, niterations=10, tolerance=10, verbose=verbose-1)
                     else:
                         gain=0
                         offset=0
@@ -1447,23 +1469,33 @@ def get_calib_data(port,program,emissivity=0.8,Temp_V=2,version=0,back_emissivit
             raise Exception("no exposure time found")
         return background,LUT,refT,gain,offset,badpixels,t_exp,cfilter,gain_error,offset_error
     except:
-        print("Program was not found")
-        return 0,0,0,0,0
+        print("Warning: Program was not found")
+        return 0,0,0,0,0,0,0,0,0,0
 
-def find_badpixels(port,gain,offset,niterations=3,tolerance=10,plot_it=False):
-    badpixels=np.zeros(np.shape(gain))
+def find_badpixels(port, gain, offset, init_bp_list=None, niterations=3,
+                   tolerance=10, plot_it=False, verbose=0):
+    badpixels = np.zeros(np.shape(gain))
 #    certain_bads=np.zeros(np.shape(gain))
 #    gainmax=12#100
 #    badpixels+=(gain>gainmax)
 #    badpixels+=(gain<0.1)
     FOV=get_FOV_mask(port)
+    # take initial bad pixels into account
+    if init_bp_list is not None:
+        if verbose>0:
+            print("use {0} initial bad pixels from file".format(len(init_bp_list)))
+        for pix in init_bp_list:
+            badpixels[pix] = 1
+        gain = IR_tools.restore_bad_pixels(gain*FOV, init_bp_list)
+    
     last_number=0
     finished=False
     n=0
     while (n<=niterations and not finished):#len(badlist)>=last_number):      
         badlist=IR_tools.find_outlier_pixels(gain*FOV,plot_it=plot_it,tolerance=tolerance)
-        gain=IR_tools.restore_pixels(gain*FOV,badlist)
-        print("number of found bad pixels: ",len(badlist))
+        gain=IR_tools.restore_bad_pixels(gain*FOV,badlist)
+        if verbose>0:
+            print("number of found bad pixels: ",len(badlist))
         n+=1
         if len(badlist)>=last_number and n>2:
             finished=True