pyview/view/idl_translation_core/ViewCalcData.py

import numpy as np
import scipy.ndimage as sci
import logging
from view.python_core.areas import read_area_file

# solution from https://stackoverflow.com/a/49480246
from .bleach_correction import fitlogdecay, get_bleach_weights

if __package__ is None or __package__ == '':  # if Master_Hannah2018.py is importing this file
    pass
else:  # if we are importing this file from a package, e.g., using view.idl_translation_core.ViewOverview
    pass
#import scipy.optimize as sco


def CalcSigMaster(flag, p1):
    #PRO CalcSigMaster, raw, dark, p, method, relative_changes

    method1 = flag["LE_CalcMethod"] #keep this renaming for the lines below
#    if ((method1 gt  1000) AND (method1 lt  2000)) THEN method1= 1000
#    if ((method1 ge  3000) AND (method1 lt  4000)) THEN method1= 3000
#    if ((method1 ge  4000) AND (method1 lt  5000)) THEN method1= 4000
#    if ((method1 ge  5000) AND (method1 lt  6000)) THEN method1= 5000
#    if ((method1 ge 10000) AND (method1 lt 19000)) THEN method1=10000
#
#; Master Procedure to calculate relative changes from rawdata for all odors
#; several different Methods can be used

# addMeasurement = fix(flag[view_MultiExp])
################################################
# these numbers do not match the IDL numbers any more
    
    if method1 == 0:
        p1 = CalcSigAll0(p1, flag) #absolute values / 1000, effectively raw data
        # in IDL, this was method 3
    elif method1 == 1:
        p1 = CalcSigAll1(p1, flag) #tbw
#    elif method1 == 2:
#        p1 = CalcSigAll2(p1, flag) #tbw
    elif method1 == 3: #simple deltaF, no thrills
        p1 = CalcSigAll3(p1, flag) #
    elif method1 == 4: #simple ratio, no thrills
        p1 = CalcSigAll4(p1, flag) #pure ratio, no thrills
        # in IDL, this was 15
#9: BEGIN
#		;calculates deltaF, F0 are initial frames to be set in program
#		;includes median bleach correction
#		;and baseline shift to inital values to be set in program
#       CalcSigAll9, raw, dark, p, relative_changes
#  END
#10: BEGIN
#		;calculates RATIOS for fura, and deltaF into air control memory
#		;includes median bleach correction
#		;and baseline shift to initial values, to be set in program
#       CalcSigAll10, raw, dark, p, relative_changes
#  END
#11: BEGIN
#		;calculates RATIOS for fura, and deltaF into air control memory
#		;no baseline shift
#       CalcSigAll11, raw, dark, p, relative_changes
#  END
#12: BEGIN
#		;calculates RATIOS for fura into odour 2, and deltaF 340 into air control memory, deltaF 380 into odour 1
#		;baseline shift
#       CalcSigAll12, raw, dark, p, relative_changes
#  END
#13: BEGIN
#		;calculates RATIOS for fura into odour 2, and deltaF 340 into air control memory, deltaF 380 into odour 1
#		;baseline shift
#       CalcSigAll13, raw, dark, p, relative_changes
#  END
#14: BEGIN
#		;calculates RATIOS for fura into odour 2, and deltaF 340 into air control memory, deltaF 380 into odour 1
#		;baseline shift
#		;constant background correction
#		;unsharp mask correction
#       CalcSigAll14, raw, dark, p, relative_changes
#  END
#15: BEGIN
#		;calculates RATIOS for fura into buffer 1, no other treatment/filter/shift/etc whatsoever
#       CalcSigAll15, raw, relative_changes
#  END
#19: BEGIN
#		;same as 9 but only central region 0.2-0.8 in x and y
#		;calculates deltaF, F0 are initial frames to be set in program
#		;includes median bleach correction
#		;and baseline shift to inital values to be set in program
#       CalcSigAll9, raw, dark, p, relative_changes
#  END
#29: BEGIN
#		;same as 9 but left and right half frames are median corrected separatedly
#       CalcSigAll29, raw, dark, p, relative_changes
#  END
#30: BEGIN
#		;calculates FURA ratio and CaGr-deltaF signals from dual sets (Silke 2001)
#       CalcSigAll30, raw, dark, p, relative_changes
#  END
#39: BEGIN
#		;same as 9 but only AL perimeter is considered
#       CalcSigAll39, raw, dark, p, relative_changes
#  END
#1000: BEGIN
#		; method for substracting a smoothed background (smooth function in IDL)- test version
#		CalcSigAll1000, raw, dark, p, relative_changes, method
#  END
    elif (method1 >= 3000 & method1 < 6000): 
        p1,flag = CalcSigAll3000(p1, flag) #tbw
#		; calculates delta F / F (3000) or FURA (4000) or raw data (5000)
#		; family of methods: first position after the 3 gives the following options
#		; setting:    0 1 2 3 4 5 6 7 8 9
#		;alPerimeter  + + + + - - - - - - select central area or alperimeter for bleach corr
#		;excludeStim  - + - + - + - + - - exclude stimulus time + 3 sec from bleach corr
#		;addNoise     - - + + - - + + - - add high frequency noise to bleach corr
#		; decimals and units give a radius in micrometers for correction of spread light
#		CalcSigAll3000,  method

#8000: BEGIN
#		CalcSigAll8000, method ;
#		;8000 is for multi-focal short ratio measurements, collapses to one depth taking the maximum
#  END
#9000: BEGIN
#		CalcSigAll9000, method ;
#		;9000 is for multi-focal short ratio measurements, collapses to one depth taking the mean
#  END
#10000: BEGIN
#		CalcSigAll10000, method ;10000 has the common blocks
#		;10000 is for multi-focal short ratio measurements, with 3D correction and many sub-options
#  END
    else: print('ViewCalcData.CalcDataMaster: flag.CalcMethod not implemented')
    return 



def CalcSigAll3000(p1, flag):
    #settings values
    FlagRatioSetting   = int(flag["LE_CalcMethod"] / 1000) # first digit. int is like floor
    RatioSetting 	= (FlagRatioSetting == 4)  # is this a FURA measurement?
    startbackground 		   = flag["LE_StartBackground"] # background ab frame ..-
    prestim_endbackground 	= flag["LE_PrestimEndBackground"] # background till n frames before stimulus
    subtractBeginning 		= (startbackground >= 0)


    #run CalcSigAll3000 for the single wavelength in p1.raw1
    p1, flag = CalcSigAll3000_raw(p1, flag) #run the full program on raw1
    # p1.sig is already the perfect calculation if not a ratio
    if RatioSetting: #fura measurement - so far we worked on raw1, i.e. 340
        # run the same thing for p1.raw2(380), therefore do some copying
        p2 = p1.copy() #make a flat copy into p2
        # raw1 is 340, raw2 is 380
        p2.raw1 = p2.raw2 #copy second wavelength into the first slot
        p2, flag = CalcSigAll3000_raw(p2, flag) #run the full program on raw2 (called raw1 here)
    # now calculate the ratio, put into p1.
    # 340/380 is p1/p2
    #but sig1 is deltaF/F, i.e. not F/F0, therefore add 1
        p1.sig1 = (p1.sig1 +1)/(p2.sig1 + 1) # times multiplyFactor for percentages, not defined in python
        #p1.raw1 = p2.raw1 #CalcSigAll3000 modifies raw1 also
        p1.metadata.bleachpar2 = p2.metadata.bleachpar #keep track of the bleach parameters
#	;since this was the ratio of the deltaF, we have to multiply with the inverse of the backframes to get the true ratio
#	;we did the delta F first because otherwise the log-bleach-correction would not work
        backBackFrame = p1.metadata.backframe/p2.metadata.backframe
        for i in range(p1.metadata.frames): #remove the effect of calculating F/F0
            p1.sig1[:,:,i] = p1.sig1[:,:,i]  * backBackFrame
        print('ViewCalcData/CalcSigAll3000: calculated ratio based on method 4xyy ',flag["LE_CalcMethod"])

#;subtract beginning
    if subtractBeginning:

        (startbackground, endbackground) = p1.metadata.background_frames

        backbackframe = np.mean(p1.sig1[:,:,startbackground:endbackground], axis=2) #takes the average for selected time frames
        for i in range(p1.metadata.frames):
            p1.sig1[:,:,i] -= backbackframe 
    return p1, flag # no need to return p1, since it is modified anyway



def CalcSigAll3000_raw(p1, flag):

#;method 4xyy, where x is the setting (see table below), and y is the radius for scatter correction
#3000 for deltaF, 4000 for Fura
#5000 for raw data

#
#;fixed settings: need to go into flags, all of them!
#
#lightThresholdFlag =0 ;not well tested yet, implemented for 3D-2photon.Nov 1007
#;switch OFF until exported to an external flag
#
    noiseseconds   	= 1 #; high frequency cutoff in seconds for lamp noise reduction, when set
#rawDataCopy    		= 0 ; puts raw data into odour set 0, if OFF puts deltaF/F of SLAVE into odour set 0
#;set rawDataCopy to 2 to put bleach corrected raw data into raw1
#;with new '5000' family of loading, do not use any more
#; calculate nr of last frame that will be used for background average
#multiplyfactor 		= 100.0   ; set to 100 to get percentage values
    lightmeanvalue	= 1000.0 #; raw data are first scaled to this mean light value
    cutedge 			= 0.2 #; border area not to be considered for bleach correction when not using area
    bleachmean 		= 0 #; use mean over time and frames as trace to do the bleach correction with; else median

    rawdatadelete	   = flag["VIEW_DeleteRawData"] != 0 # ; deletes raw data at the end, saves disk space
    logexcludeseconds = flag["LELog_ExcludeSeconds"] #]) ; how many seconds after stimulus onset to exclude, when set
    bleachstartframe 	= flag["LE_BleachStartFrame"] #] ); set to 2 to exclude 2 frames at the beginning from log-fit
# frames before bleachstartframe are ignored for bleach correction
    loginitialfactor  = flag["LELog_InitialFactor"]#] ) ; weight factor for curve stretch before stimulus onset
# frames after bleachstartframe and before stimulus-logexccludeseconds are valued with this weight

# for example: curve 5 Hz, 100 Frames long, odor at frame 21
# LogExcludeSeconds=5, BleachStartFrame=3, LogInitialFactor=2
# LogFitting will be made using frames 3:20 and 46:99, but 3:20 will count double
# prestim_endbackground 	(see below) is also used for bleach correction

#variable settings
    flagratiosetting   = int(flag["LE_CalcMethod"] / 1000) # first digit. int is like floor
    myMethod       = int(flag["LE_CalcMethod"] / 100) % 10 #second digit. contains settings for bleach correction
    alperimeter    = myMethod in [0,1,2,3]      #;for 0,1,2,3
    excludestimulus= myMethod in [1,3,5,7] # for 1,3,5,7 (,9)
    addnoise       = myMethod in [2,3,6,7] # for 2,3,6,7
    nobleach       = myMethod in [9]
    airbleach      = myMethod in [8]
#AirBleachAddNoise = 0 ; local setting  - option not implemented yet
#;IF AirBleach THEN alPerimeter = 1 ;use this line for bleach in AL perimeter only

#; setting:    0 1 2 3 4 5 6 7 8 9
#;alPerimeter  + + + + - - - - C - ; with 8 bleaching is in coordinates only (variable: CoorPerimeter = AirBleach)
#;excludeStim  - + - + - + - + - -
#;addNoise     - - + + - - + + - -
#;no bleach    - - - - - - - - - +
#;air bleach   - - - - - - - - + - ;correct with bleach parameters taken from air trial

#;ratiosetting of 3: deltaF/F
#;ratiosetting of 4: make ratio of odor buffer 0 / buffer1
#    ratiosetting 	= (FlagRatioSetting == 4) #not used because the fura is controlled from outside
    # in here, I only work on the main data raw1.
    returnrawdata 	= (flagratiosetting == 5)
    #returnrawdata returns the bleach corrected data to raw1

#external settings
    shrinkfactor 	= flag["LE_ShrinkFaktor"]
#Memory_3Dim				= fix(flag[VIEW_No4Darray]) ; use arrays with only one odor-buffer to save space
    startbackground 		= max(flag["LE_StartBackground"],0) # background ab frame ..- at least 0
#    subtractbeginning 	= (startbackground >= 0)
    prestim_endbackground 	= flag["LE_PrestimEndBackground"] # background till n frames before stimulus
#firstOdor 				= fix(flag[LE_FirstBuffer])
#   radius 					= flag.RM_Radius # necessary for coordinate analysis
# the next lines not translated to Python, because they relate to the way the data was put into the sig variable
# if RatioSetting 		THEN 	firstOdor = 0 ; 0 contains master
#if (firstOdor eq -1) 	then begin ;settings as befor the introduction of LE_usefirstbuffer
#	IF (total(raw1(*,*,*,0)) eq 0) then firstOdor = 1 ELSE firstodor = 0
#endIF

    (startbackground, endbackground) = p1.metadata.background_frames


#SETTING should I do unsharpmask?
    smoothvalue  = int(flag["LE_CalcMethod"] % 100) # last two digits, radius of filter in um
    unsharpmask  = smoothvalue > 0

    if airbleach:
        print('ViewCalcData.CalcSigAll3000: airBleach not implemented yet in Python')
#	IF RatioSetting Then Begin
#		airBleach = 0
#		print, 'CalcSigAll3000: AirBleach not implemented with Ratio dyes'
#	endIF
#	IF (firstOdor gt 0) Then Begin
#		;airBleach = 0
#		firstodor = 0
#		print, 'CalcSigAll3000: using AirBleach with flag[LE_FirstBuffer] greater 0 (???)'
#	endIF
#	IF (total(raw1(*,*,*,0)) eq 0) then begin
#		airBleach = 0
#		firstodor = 1
#		print, 'CalcSigAll3000: AirBleach not useful without data in buffer 0 (air is missing)'
#	endIF
#	IF AirBleachAddNoise THEN begin
#		print, '************CalcSigAll3000 - local flag set to addNoise with AirBleachAddNoise'
#		addNoise = 1
#	endIF

#loadct, SO_MV_colortable



#
#;define variables
#sig1 = fltarr(p1.format_x, p1.format_y, p1.frames, p1.odors+1)
#backframe        = fltarr(p1.format_x, p1.format_y)
# Python: will use p1.sig1 and p1.backframe

###### start calculation: create maskframe, either via alperimeter, or via cutedge
    if alperimeter:
#	;get selected AL area
#        if flag.RM_differentViews > 0: #not implemented for now - 
#   		;special section Tilman
#    		;areaFileName = flag[stg_OdorMaskPath]+strmid(Flag[stg_ReportTag],0,14)+p1.viewLabel+' '+strmid(Flag[stg_ReportTag],14)+'.area'
#    		;print, 'CalcSigAll3000: special Tilman section for area file!!'
#    	;standard section
#            areafilename = flag.STG_OdormaskPath+flag.STG_ReportTag+str(p1.viewLabel)+'.Area'
#        else:
#         maskframe = IDL.restore_maskframe(flag) #restores maskframe from file
        try:
            area_filename = flag.get_existing_area_filepath()

        except FileNotFoundError as fnfe:
            maskframe = np.ones((p1.metadata.format_x, p1.metadata.format_y))
        else:
            maskframe = read_area_file(area_filename)

#	;now AL perimeter is in variable maskframe
        #apply shrinkfactor
        if (shrinkfactor not in [0,1]): #both 0 and 1 mean 'no shrink'
            # copy from ViewLoadData.loadTILL
            shrink = 1/shrinkfactor # flag of 2 means half the size in each dimension
            import scipy.ndimage as scind
        # the command used in IDL for shrink was REBIN
        #raw1(*,*,i,1) = rebin(image(0:p1.format_x*shrinkFactor-1, 0:p1.format_y*shrinkFactor-1) ,p1.format_x, p1.format_y) ;
        # the following shrinks in x and y, not in t, with polinomial order 1. 
        # "nearest" is for treatment at the borders
            maskframe = scind.zoom(maskframe,(shrink,shrink), mode='nearest')
          
        #check for consistent size,
#	;safety check for unequal size
#	;remove last row(s)/column(s) in case of original size incompatible with shrink factor
#	;e.g. shrinkfactor 2 with format_x 173, limit mask to 0..172
#	maskFrame = maskFrame(0:p1.format_x*shrinkfactor-1,0:p1.format_y*shrinkfactor-1 )
#	;resize maskFrame considering shrinkfactor
#	maskSize = size(maskFrame)
#	maskFrame = rebin(maskFrame,maskSize[1]/shrinkfactor,maskSize[2]/shrinkfactor)
#	maskSize = size(maskFrame)
        if (maskframe.shape != (p1.metadata.format_x,p1.metadata.format_y)): # not correct size
            print('ViewCalcData.Calcsigall3000: AREA file has wrong size - maybe make with binning unequal 1?')
            return #jump out of the program
    else: # there is no alPerimeter, i.e. no .Area file
#	;select central region of AL for correction
        maskframe = np.zeros((p1.metadata.format_x, p1.metadata.format_y), dtype=int)
        maskframe[round(p1.metadata.format_x*cutedge):round(p1.metadata.format_x*(1-cutedge)),round(p1.metadata.format_y*cutedge):round(p1.metadata.format_y*(1-cutedge))]=1
#	maskFrame[p1.format_x*CutEdge:p1.format_x*(1-CutEdge),p1.format_y*CutEdge:p1.format_y*(1-CutEdge)]=1
#	;now AL perimeter is in variable maskframe
        print('CalcSigAll3000 - reference area for F and for Bleach is centre with edges cut at ',cutedge)

#########  CoorPerimeter #################
# the following CoorPerimeter section appears not in use - not translated to python yet
# I only found the flag in AirBleach.
# The idea is: get the bleach function only from within the coordinates of .coor
# easy to add if needed
#;CoorPerimeter
#;If CoorPerimeter is set, select the area delimited by coordinates for bleach correction maskframe
#;section above becomes obsolete, because maskframe is overwritten
#IF CoorPerimeter THEN begin
#	print, 'Using coordinates for bleach control (calcsigall3000)'
#	;load coordinates
#	;get file name of coordinate file (.coor)
#	;Get Information about odours from external file
#	   	IF fix(flag[RM_differentViews]) THEN begin
#    			GlomeruliListFile =  flag[stg_OdorMaskPath]+flag[stg_ReportTag]+p1.viewLabel+'.coor'
#    		endIF else begin
#    			GlomeruliListFile =  flag[stg_OdorMaskPath]+flag[stg_ReportTag]+'.coor'
#    	endELSE
#     		;reportFile = flag[stg_OdorReportPath]+flag[stg_ReportTag]+'.expGlo'
#	   	shiftMaskeX = p1.shiftX
#	   	shiftMaskeY = p1.shiftY
#	;empty maskframe
#	maskFrame = bytarr(p1.format_x, p1.format_y)
#   	;read coordinates
#    readListe, GlomeruliNum, liste, GlomeruliListFile, column4=separateLayers
#   	glomerulinum = fix(glomerulinum)
#   	;got the liste of coordinates & identities
#    ;shiftMaskeX, shiftMaskeY
#   	liste(0,*) = liste(0,*) + shiftMaskeX
#   	liste(1,*) = liste(1,*) + shiftMaskeY
#    ;consider shrinkfactor
#    IF (shrinkFactor gt 1) then begin
#			liste(0,*) = liste(0,*)/shrinkFactor
#			liste(1,*) = liste(1,*)/shrinkFactor
#    endIF
#    ;now get infos about each coordinate
#	for glo=0, GlomeruliNum-1 do begin
#					x = liste(0,glo)
#			    	y = liste(1,glo)
#			    	xborders = [x-radius,x+radius]
#			    	IF ((xborders(0) lt 0) OR (xborders(1) gt (p1.format_x-1))) THEN begin
#			    		xborders = (xborders > 0 ) < (p1.format_x-1)
#			    		print, 'x coordinates clipped: ',x, xborders,' in glo ',glo
#			    	ENDif
#			    	yborders = [y-radius,y+radius]
#			    	IF ((yborders(0) lt 0) OR (yborders(1) gt (p1.format_y-1))) THEN begin
#			    		yborders = (yborders > 0 ) < (p1.format_y-1)
#			    		print, 'y coordinates clipped: ',y, yborders,' in glo ',glo
#			    	ENDif
#			    ;now set this area to 1 in maskFrame
#				maskFrame(xborders(0):xborders(1),yborders(0):yborders(1))=1
#		endfor
#endIF ;CoorPerimeter

#positions = where(maskframe)
#;positions contains the pixels to be considered for bleach correction
# in python, do not use np.multiply(raw1,maskframe) 

#dividend   = endbackground-startbackground+1     ;find how many frames for F
#there are no multiple odors in Python, but there may be raw1 and raw2
#         
#;first run through data for light scattering correction, and setting mean to 10000
#for od = firstOdor, p1.odors  do begin
#
#	;calculate F (float)
#    backframe  = total(raw1(*,*,startbackground:endbackground,od),3,/DOUBLE) ;makes a sum of the third dimension (frames)
#	backframe  = float(backframe) / dividend
############### background frame, corresponds to F0

    p1.metadata.backframe = np.mean(p1.raw1[:,:,startbackground:endbackground], axis=2)
    print('ViewCalcData/CalcSigAll3000_raw: background start/stop is frame: ',startbackground, endbackground)
    #takes the average for selected time frames
    
#;calculate intensity correction, only within maskframe
#    backvalue = np.mean(np.multiply(p1.backframe,maskframe))
    backvalue = np.mean(p1.metadata.backframe[maskframe != 0])
#	backValue  = mean(backframe(positions), /double) ; average light values in the backframe area
    if (backvalue == 0): backvalue = 1 #safety net to avoid dividing by 0
    backfactor = lightmeanvalue/backvalue
# apply backFactor outside loop
    p1.sig1 = p1.raw1 * backfactor
########### now sig1 was created, as light intensity corrected copy of raw1
    #       so from now on work on sig1
    
#            buffer[:,:,i] = sci.gaussian_filter(buffer[:,:,i], filterSpaceSize)
#scipy.ndimage.filters.gaussian_filter(input, sigma, order=0, output=None, mode='reflect', cval=0.0, truncate=4.0)

#	; correct data for scattered light
    if unsharpmask: #light scatter correction
        smoothfactor	= flag["VIEW_ScatterLightFactor"]
        # adjust for pixel size - since smoothvalue is given in um
        smX 		= smoothvalue / p1.metadata.pixelsizex
        smY 		= smoothvalue / p1.metadata.pixelsizey
        smoothvalue = (smX+smY)/2.0
        if (smX != smY): print('unequal pixel size not implemented yet - averaging x and y value')
        print('Scattered light correction with ',smoothvalue, ' pixels for ',p1.metadata.pixelsizex*smoothvalue,' microns. Factor:', smoothfactor)
        print('CalcSigAll3000: a negative factor means scattered light correction, a positive factor a spatial high-pass filter')

        for i in range(0,p1.sig1.shape[2]): # go through frames
            # copy a frame, smooth it, copy back
            oneframe = p1.sig1[:,:,i]
#	for i= 0, p1.frames-1 do begin ;correct data first, write corrected data into final dest
#		;dark-frame correction # dark frame deprecated
#		oneFrame = float(raw1(*,*,i,od) - dark1)
#		;intensity correction # done outside the frame loop
#		oneFrame = oneFrame * backFactor
#		;light scatter correction
#		IF unsharpMask THEN begin
            if smoothfactor > 0:
#			IF (SmoothFactor gt 0) THEN begin
#				oneFrame = oneFrame + SmoothFactor*(oneFrame - smooth(oneFrame, smoothValue, /EDGE_TRUNCATE))
                #I use sci.gaussian_filter for the IDL smooth command
                #mode='nearest' extends the last value as a constant
                oneframe = oneframe + smoothfactor*(oneframe - sci.gaussian_filter(oneframe, smoothvalue, mode='nearest'))
#			endIF
            if smoothfactor < 0:
#			IF (SmoothFactor lt 0) THEN begin
#			;correction Jan 05: the upper line amplifies noise - in fact, rather than a scattered light correction
#			;it is a spatial high-pass filter!
#			;the next line smooths the unsharp image
#			;for backward compatibility, this has to be indicated by a NEGATIVE smoothFactor
#				oneFrame = oneFrame - SmoothFactor*smooth((oneFrame - smooth(oneFrame, smoothValue, /EDGE_TRUNCATE)),smoothvalue, /edge_truncate)
                oneframe = oneframe - smoothfactor*sci.gaussian_filter((oneframe - sci.gaussian_filter(oneframe, smoothvalue, mode='nearest')), smoothvalue, mode='nearest')
#			endIF
#		endIF
#		;write corrected data temp. into sig1
#		sig1(*,*,i, od) = oneFrame
            p1.sig1[:,:,i] = oneframe
#	endFOR
#endfor   ; odornr

#;sig1 now contains the scatter-light corrected raw data, with mean at lightmeanvalue
            # copy into raw1, or delete raw1
    if rawdatadelete:
        p1.raw1 = 0 #save memory space
    else:
        p1.raw1 = p1.sig1.copy() # raw now contains scatter corrected "raw" data
#IF rawDataDelete then begin
#	Raw1 = 0
#	print, 'CalcSigAll3000: raw data deleted'
#endIF else begin
#	;copy back into raw
#	raw1 = sig1  #why back into raw?
#endELSE
    

#;make a copy of this corrected data into the experimental foto
#;only if the foto does not yet exist or is empty
    if not hasattr(p1, "foto1"):
        p1.foto1 = np.mean(p1.sig1[:, :, startbackground :endbackground], axis=2)
    elif np.all(p1.foto1 == 0):
        p1.foto1 = np.mean(p1.sig1[:, :, startbackground :endbackground], axis=2)


## lightThresholdFlag apparently was not implemented yet well
#IF (lightThresholdFlag eq 1) then begin
#	;apply light threshold to data
#	;that is: where there is no dye, there cannot be a signal
#	;threshold calculated from foto
#	lightThreshold = 800
#	sig1 = sig1 > lightThreshold
#endIF ;flag



#;calculate deltaF for each buffer
    p1.sig1 = calc_deltaF(p1.sig1, refRange=[startbackground, endbackground])


#;run through deltaF data, bleach and noise correction for each 'odour', i.e. wavelength separately
    if nobleach:
        print('No logarighmic bleach correction in CalcSigAll3000')
    else:
#$
#	;bleach correction!
#	;necessary data structures
#	;excluded frames at the beginning are defined in fitlogdecay
#	correctArray    = fltarr(p1.frames)
#	AllcorrectArray = fltarr(p1.odors+1,p1.frames)
#	AllOdourArray   = fltarr(p1.odors+1,p1.frames)

#calculate correctarray as mean or median
    #convert flag
# - limit this to the mask
#        if bleachmean == 1:
#            correctarray = np.mean(p1.sig1, axis=(0,1))
#        else:
#            correctarray = np.median(p1.sig1, axis=(0,1))
        #not good in terms of syntax
        correctarray = np.zeros(p1.metadata.frames)
        positions = maskframe != 0
        if bleachmean == 1:
            for i in range(p1.metadata.frames):
                correctarray[i] = np.mean(p1.sig1[:,:,i][positions])
        else:
            for i in range(p1.metadata.frames):
                correctarray[i] = np.median(p1.sig1[:,:,i][positions])



##############NOW WORK ON THE BLEACH CORRECTION; that is: calculate the weight array.

        weights = get_bleach_weights(flags=flag, p1_metadata=p1_metadata, exclude_stimulus=excludestimulus)
############
        print(weights)
        print(prestim_endbackground, logexcludeseconds, p1.metadata.frequency, stimulus_on, loginitialfactor, bleachstartframe)
        if airbleach:
            print('CalcSigAll3000: ********* AirBleach not implemented in Python yet')
#				;calculate log function only for AIR
#				IF (od eq 0) THEN fitlogdecay, correctArray, Airdetract, weights
#				detract = Airdetract
#				print, 'LogBleach done on reference measurement!'
        else:
            measurement_label = flag.get_default_measurement_label()
            (detract, opt_parms) = fitlogdecay(correctarray, weights, not flag["VIEW_batchmode"], measurement_label)
            p1.metadata.bleachpar = opt_parms
            print('CalcSigAll3000: Logbleach correction done')
#IDL: pro FitLogDecay, lineIn, fittedOut, weights, verbose ;, noConverge

        if addnoise:
#			;now add high frequency noise (lamp noise) to the fitted exponential
#            dummy = detract[0:BleachStartFrame] #;do not smooth before bleachStartFrame
# $$ implement smooth in the next line, defined above
#            detract = detract + correctarray - smooth(correctarray, p1.frequency*noiseSeconds)
#$$ or gaussian filter buffer[:,:,i] = sci.gaussian_filter(buffer[:,:,i], filterSpaceSize)
            detract = detract + correctarray - sci.gaussian_filter(correctarray, p1.metadata.frequency*noiseseconds)
            print('CalcSigAll3000: addnoise done')
#            detract[0:BleachStartFrame] = dummy[:]

#		;subtract these values
        for i in range(p1.metadata.frames):
            p1.sig1[:,:,i] -= detract[i]

#		;MONITOR the action taken
#		allcorrectArray(od,*)=detract(*) ;this is just for plotting the monitor
#		allOdourArray(od,*)  =correctArray(*) ;this is just for plotting the monitor

##;calculate bleach corrected raw data and put it into raw
            # discontinued, because implemented in 5000 family with returnrawdata
#        if (rawdatacopy == 2):
#            for i in range(p1.frames-1):
##		raw1(*,*,i,*) = (sig1(*,*,i,*) * backFrame(*,*,*))/multiplyfactor ;raw data
#                p1.raw1[:,:,i] = (p1.sig1[:,:,i] * p1.backframe[:,:])/multiplyfactor #
#                # this inverts the deltaF/F ;raw data
##	endFor
        print('CalcSigAll3000: bleach corrected deltaF data in sig1')
#;calculate bleach corrected raw data and put it into sig
        if returnrawdata:# setting is 5000
            for i in range(p1.metadata.frames): # do begin
                # do the inverse of deltaF/F
                p1.raw1[:,:,i] = ((p1.sig1[:,:,i]+1) * p1.metadata.backframe) # /multiplyfactor# ;raw data
            #this only works because backframe uses the same frames as deltaF/F
            print('CalcSigAll3000: bleach corrected raw data calculated')
#        else: #endIF ELSE begin
#  ;shift back to mean value of MultiplyFactor
#  ;'temporary' saves memory space
#            p1.sig1 += multiplyfactor
#            print('CalcSigAll3000: check this line - why would MultiplyFactor be added??')

# in python, plot the curves not implemented (yet?)
#  ;plot the curves used in log-correction
#
#	window, 16, xsize=512, ysize=512
#	minplot = floor(min([allcorrectarray(firstodor:p1.odors,*),allodourarray(firstodor:p1.odors,*)]))
#	maxplot = ceil(max([allcorrectarray(firstodor:p1.odors,*),allodourarray(firstodor:p1.odors,*)]))
#	;minplot = 98.0
#	;maxplot = 102.0
#	plot, correctarray(*), yrange=[minplot,maxplot], color=245, thick=3, ystyle=1, /nodata
#	;plot weight values as curve and bottom line
#	oplot, minplot+weights(*), thick=2, color=255 & correctarray(*)=minplot & oplot, correctarray(*), thick=1, color=255
#	;plot log fits
#	for od = firstOdor, p1.odors  do begin
#		oplot, allcorrectArray(od,*), thick=3, color=240-30*od
#		oplot, allOdourArray(od,*), thick=1,  color=240-30*od
#	end
#	;plot stimulus bar
#	xaxis     = FINDgen(p1.frames)
#	stimlow   = fltarr(p1.frames)
#	stimlow(*)= minplot
#	oplot, xaxis(p1.stimulus_on:p1.stimulus_end), stimlow, thick=5, color = 245
#	;second stimulus; if no second stimulus is given, this should be from 0 to 0
#	IF ((p1.stim2on gt 0) AND (p1.stim2off gt 0) AND (p1.stim2on lt p1.frames-1) AND (p1.stim2on  lt p1.frames-1)) THEN begin
#		oplot, xaxis(p1.stim2on:p1.stim2off), stimlow, thick=5, color = 245
#	endIF

#endELSE;nobleach
    return p1, flag #  no return, p1 is modified already. Try return instead now (8/19)


def calc_deltaF(imgData, refRange):
    '''    
    Input: 3D matrix x,y,t, refRange
    calculates a crude deltaF.
    F0 is the entire movie (default), or the interval refRange (e.g. [14,21])
    :param refRange: list of length 2, frame numbers (indices) of the starting and ending frame of the range to use for
    baseline fluorescence calculation. Note that the ending frame is included in baseline calculation.
    '''

    referenceF = np.mean(imgData[:, :, refRange[0]: refRange[1] + 1], axis=2)

    dead_mask = referenceF == 0
    if np.any(dead_mask):
        # keep pixels that averages 0 over the reference range unmodified. these are most likely pixels dead at 0.
        referenceF[dead_mask] = 1

        logging.getLogger("VIEW").warning("F0 for some pixels were 0, not normalizing these pixels!")

    # convert to txy because framewise division is easier
    imgData_txy = np.moveaxis(imgData, source=-1, destination=0)

    deltaFdata_txy = (imgData_txy / referenceF) - 1

    # convert back
    deltaFdata = np.moveaxis(deltaFdata_txy, source=0, destination=-1)

    logging.getLogger("VIEW").info(
        f"Baseline fluorescence was calculated by averaging frames numbered {refRange[0]}-{refRange[1]}")

    return deltaFdata


def CalcSigAll0(p1, flag):
#pro CalcSigAll3, raw, dark, p, absolute_values
#; Procedure not to calculate any signals
#; Author: Giovanni, Jasdan
#;parameters: (from loadexp.pro)
#; raw        : 4 dimensional array (x,y,frames,odors)
#;              odor=0:   background/air
#;              odor=1... odors from experiment
#; dark       : darkframe of ccd camera
#; p          : parameterset
#;results:
#;relative_changes : 4 dimensional float array with signals for all odor and background
#absolute_values  = fltarr(p.format_x, p.format_y, p.frames, p.odors+1)
#backframe        = fltarr(p.format_x, p.format_y)
#for odor = 0, p.odors  do begin
#	backframe(*) = float(dark)
#	for i= 0, p.frames-1 do begin
#   	  absolute_values(*,*,i, odor) = raw(*,*,i,odor) - backframe(*,*)
#	endfor
#endfor   ; odornr
#absolute_values = absolute_values  / 1000
    
    p1.sig1 = p1.raw1.astype(float) / 1000
    print('ViewCalcData/CalcSigAll0: Signals calculated (method 0). sig1 = raw1/1000')
#end ; of program
    return p1


def CalcSigAll1(p1, flag):
    # pro CalcSigAll3, raw, dark, p, absolute_values
    # ; Procedure not to calculate any signals
    # ; Author: Giovanni, Jasdan
    # ;parameters: (from loadexp.pro)
    # ; raw        : 4 dimensional array (x,y,frames,odors)
    # ;              odor=0:   background/air
    # ;              odor=1... odors from experiment
    # ; dark       : darkframe of ccd camera
    # ; p          : parameterset
    # ;results:
    # ;relative_changes : 4 dimensional float array with signals for all odor and background
    # absolute_values  = fltarr(p.format_x, p.format_y, p.frames, p.odors+1)
    # backframe        = fltarr(p.format_x, p.format_y)
    # for odor = 0, p.odors  do begin
    #	backframe(*) = float(dark)
    #	for i= 0, p.frames-1 do begin
    #   	  absolute_values(*,*,i, odor) = raw(*,*,i,odor) - backframe(*,*)
    #	endfor
    # endfor   ; odornr
    # absolute_values = absolute_values  / 1000

    p1.sig1 = p1.raw2.astype(float) / 1000
    print('ViewCalcData/CalcSigAll1: Signals calculated (method 0). sig1 = raw2/1000')
    # end ; of program
    return p1


def CalcSigAll3(p1, flag):
    '''
    Calculate deltaF, with no thrills at all
    in IDL, this was ???
    '''

    p1.sig1 = calc_deltaF(p1.raw1, refRange=p1.metadata.background_frames)

    logging.getLogger("VIEW").info(
        f'CalcSigAll3: delta F with background: {p1.metadata.background_frames}')
    return p1


def CalcSigAll4(p1, flag):
    '''
    Calculate the ratio, with no thrills at all
    in IDL, this was CalcSigAll15
    '''
#pro CalcSigAll15, raw, relative_changes
#deleteRaw = 1
#multiplyFactor = 100 ; factor to multiply data with when making ratio
#relative_changes = raw
#IF deleteRAW then raw = 0
#relative_changes = float(relative_changes)
#;calculate RATIOS
#relative_changes(*,*,*, 1) = (multiplyFactor*relative_changes(*,*,*, 0)) / relative_changes(*,*,*, 1)
#relative_changes(*,*,*, 0) = 0 ;delete buffer 0
#print, 'plain ratios in buffer 1'
#end ; of program
    p1.sig1 = p1.raw1 / p1.raw2
    print('ViewCalcData/CalcSigAll4: Signals calculated (method 4). Fura, no thrills')
    return p1