1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471 |
- # -*- coding: utf-8 -*-
- """
- Created on Thu May 31 17:34:29 2018
- @author: Giovanni Galizia
- """
- import pandas as pd
- import scipy as sp
- from scipy import ndimage
- import numpy as np
- import matplotlib.pyplot as plt
- import os
- from view.python_core.p1_class.filters import apply_filter
- from view.python_core.paths import get_existing_raw_data_filename
- from view.python_core.io import read_tif_2Dor3D
- import tifffile
- import pathlib as pl
- import logging
- def load_pst(filename):
- """
- read tillvision based .pst files as uint16.
- """
- # filename can have an extension (e.g. .pst), or not
- # reading stack size from inf
- #inf_path = os.path.splitext(filename)[0] + '.inf'
- #this does not work for /data/030725bR.pst\\dbb10F, remove extension by hand,
- #assuming it is exactly 3 elements
- if filename[-4] == '.':
- filename = filename[:-4] #reomove extension
- meta = {}
- with open(filename+'.inf','r') as fh:
- # fh.next()
- for line in fh.readlines():
- try:
- k,v = line.strip().split('=')
- meta[k] = v
- except:
- pass
- # reading stack from pst
- shape = sp.int32((meta['Width'],meta['Height'],meta['Frames']))
- expected_units = np.prod(shape)
- assert pl.Path(f"{filename}.pst").stat().st_size >= 2 * expected_units, \
- f"Expected atleast {2 * expected_units} bytes in {filename}.pst. Found {pl.Path(filename).stat().st_size}"
- raw = sp.fromfile(filename+'.pst', dtype='int16', count=expected_units)
- data = sp.reshape(raw,shape,order='F')
- # was swapping x, y axes; commented out to retain original order
- # data = data.swapaxes(0,1)
- # data is upside down as compared to what we see in TillVision
- data = np.flip(data, axis=1)
- data = data.astype('uint16')
- return data
- def read_lsm(path):
- """ takes a path to a lsm file, reads the file with the tifffile lib and
- returns a np array
- """
- data = tifffile.imread(path)
- Data_cut = data[0,0,:,:,:] # empirical ...
- Data_cut_rot = sp.swapaxes(Data_cut,0,2)
- Data_cut_rot_flip = sp.flip(Data_cut_rot, axis=1)
- return Data_cut_rot_flip
- def loadTILL(flag, p1):
- filename = get_existing_raw_data_filename(flags=flag, p1_metadata=p1, extension=".pst")
- #if not exist try FID format version oct18
- # if still not exist try FID format old
- print('ViewLoadData.loadTILL. Loading data :',filename)
- # write data directly into p1 structure: p1.Raw1
- p1["raw1"] = load_pst(filename.rstrip(".pst"))
- # saving the complete filename of raw data for reference
- p1["full_raw_data_path"] = filename
- # use shrinkfactor to shrink
- if (flag["LE_ShrinkFaktor"] not in [0,1]): #both 0 and 1 mean 'no shrink'
- shrink = 1/flag["LE_ShrinkFaktor"] # 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
- p1.raw1 = scind.zoom(p1.raw1,(shrink,shrink,1), order=1, mode='nearest')
- # put shape into p structure
- datasize = p1.raw1.shape
- p1.format_x = datasize[0]
- p1.format_y = datasize[1]
- p1.frames = datasize[2]
- # IDL: MedianCorrection, oneOdour, fix(flag(csm_median)); mean filter added in pyVIEW
- # apply median filter first, then mean filter, depending on flags
- p1.raw1 = apply_filter(matrix_in=p1.raw1, view_flags=flag, filter_type="median")
- p1.raw1 = apply_filter(matrix_in=p1.raw1, view_flags=flag, filter_type="mean")
- # calculate foto
- p1.foto1 = p1.raw1[:,:,1:6].mean(axis = 2) #average of frames 1:6 p1.raw1[:,:,1:6]
- # IDL: foto1(i,j) = mean(raw1(i,j,1:5,0))
- #
- # calculate darkfoto: deprecated
- #
- if (flag["LE_AskForAir"] != 0):
- print('ViewLoadData.loadTILL. Loading Air not implemented yet')
- # use gettillloginfo to get the info for p1.control
- # load that measurement p_air["raw1"] = load_pst(filename)
- # check that shape is equal
- #
- return p1 #p1 contains all data for this measurement
- def loadFURA(flag, p1):
- # loads FURA measurements
- # here I deviate from IDL. I just call loadTILL two times (in IDL, that code was effectively duplicated)
- # but I play around with the p1 structure.
- p2 = p1.copy() # use p2 for the second wavelength.
- p2.dbb1 = p2.dbb2 # BUT loadTILL does not know, so swap the dbb info
- # load wavelength 1
- p1 = loadTILL(flag, p1)
- # load wavelength 2
- p2 = loadTILL(flag, p2)
- # copy data of wavelength 2 into 1, since they belong together
- p1["raw2"] = p2.raw1
- p1["foto2"] = p2.foto1
- if (flag["LE_AskForAir"] != 0):
- print('ViewLoadData.loadFURA. Loading Air not implemented yet')
- # in fact, it was never implemented in IDL either
- return p1 #p1 contains all data for this measurement
- def loadZeiss(flag, p1):
- filename = get_existing_raw_data_filename(flags=flag, p1_metadata=p1, extension=".lsm")
- print('ViewLoadData.loadZeiss: Loading data :',filename)
- #;load zeiss data set
- print('loadZeiss: loading Zeiss dataset ', filename)
- dataIn = read_lsm(filename)
- #;define data sets for IDL-View
- p1.raw1 = dataIn #;for old empty air format
- # saving the complete filename of raw data for reference
- p1["full_raw_data_path"] = filename
- datasize = p1.raw1.shape
- p1.format_x = datasize[0]
- p1.format_y = datasize[1]
- p1.frames = datasize[2]
- #;median correction
- p1.raw1 = MedianCorrection(p1.raw1, flag)
- # ;define foto as the mean of the first six frames
- p1.foto1 = p1.raw1[:, :, 1:6].mean(axis=2)
- print('loadZeiss: data loaded')
- return p1 #loadZeiss
- def load_VIEW_tif(flag, p1):
- filename = get_existing_raw_data_filename(flags=flag, p1_metadata=p1, extension=".tif")
- logging.getLogger("VIEW").info(f"Loading VIEW-tif from: {filename}")
- try:
- p1.raw1 = read_tif_2Dor3D(filename)
- except Exception as e:
- raise IOError(f"Error reading {filename}")
- # saving the complete filename of raw data for reference
- p1["full_raw_data_path"] = filename
- datasize = p1.raw1.shape
- assert len(datasize) == 3, \
- f"{filename} cannot be loaded as a VIEW-tif as it has {len(datasize)} dimensions (3 expected)"
- p1.format_x = datasize[0]
- p1.format_y = datasize[1]
- p1.frames = datasize[2]
- # ;median correction
- p1.raw1 = MedianCorrection(p1.raw1, flag)
- # ;define foto as the mean of the first six frames
- p1.foto1 = p1.raw1[:, :, 1:6].mean(axis=2)
- return p1
- def create_raw_data667(p1_metadata, peaksignal):
- '''
- 'this is the backup copy of the old create_raw_data666, that contained only square glomeruli
- '''
- print('ViewLoadData/LoadRaw667: generating test data with many fixed settings!')
- stimonset = p1_metadata.stimulus_on
- stimoffset = p1_metadata.stimulus_end
- ### for playing around start after #:
- # peaksignal, stimonset, stimoffset = 10, 25, 35 # - is in command line now
- # local flags
- lightlevel = 1000 # baseline value of each pixel.
- # at the setup, you would set exposure time and light level accordingly
- randomseed = 0 # set to None for irreproducible random numbers
- # add chipnoise
- darknoise = True
- darknoise_mean = 100 # mean of background values in the chip
- darknoise_amplitude = 20 # max amplitude of background chip noise
- # I assume chip noise to be random over time
- #
- xSize, ySize, tSize = 172, 130, 80 # 120, 99, 80
- rampwidth = 10 # width of the ramp, left 3x, right 2x
- stepsize = 13 # must be a divisor of ySze
- # glomerular locations are (make sure to exclude rampwidth)
- glo1 = (50, 70, 50, 70) # x,x - y,y
- glo2 = (60, 80, 60, 80) # x,x - y,y
- glo3 = (65, 75, 45, 55)
- glo4 = (70, 90, 30, 50)
- glo5 = (50, 65, 30, 45)
- glo6 = (70, 90, 50, 70) # off response
- glo7 = (50, 70, 90, 110) # biphasic response
- glo8 = (75, 95, 90, 110) # negative response
- shotnoise = True # add photon noise
- # add bad pixels on the CCD chip, i.e. pixels with value 0
- badpixels = True
- badpixels_clip = 0.98 # 0.9 for 10% bad pixels
- #
- # add global bleaching
- bleach = True
- bleach_percent = 10 # percentage bleaching
- bleach_tau = tSize / 5 # assume decay to 37% of percentage over this length
- # add lamp noise
- lampnoise = 0.1 # value in percent. 0 for no noise. Value is PEAK noise
- # add scattered light
- scatterlight = 0 # 2 # unit is gaussian pixels; 0 for no scattered light
- # change code, do exclude lateral step image
- # TODO
- # - add light scatter (smooth): modify
- # - add movement
- #
- def applytimecourse(matrix, glo_location, oneline):
- '''
- applies oneline as a factor (length: tSize)
- to square window (glomerulus) defined with upperleft, lowerright
- glo_location is (x_left,x_right,y_down, y_up)
- (y_down < y_up; but safety net included)
- '''
- x1 = np.min((glo_location[0], glo_location[1]))
- x2 = np.max((glo_location[0], glo_location[1]))
- y1 = np.min((glo_location[2], glo_location[2]))
- y2 = np.max((glo_location[3], glo_location[3]))
- matrix[x1:x2, y1:y2, :] = matrix[x1:x2, y1:y2, :] * oneline.reshape(1, 1, tSize)
- return matrix
- # set random numbers
- np.random.seed(randomseed)
- # create still image
- # with average light level as visible to the experimenter
- # therefore subtract darknoise_mean from lightlevel
- stillimage = np.full((xSize, ySize), (lightlevel - darknoise_mean))
- # create a ramp from 0 to 2*lightlevel
- oneline = np.linspace(0, 2 * lightlevel, ySize).reshape(1, ySize)
- stillimage[0:rampwidth, :] = oneline # first band, going up
- stillimage[rampwidth:2 * rampwidth, :] = np.flip(oneline) # first band, going up
- # create a ramp in the center, in order to avoid border effects
- oneline = np.linspace(0, 2 * lightlevel, ySize // 2).reshape(1, ySize // 2)
- stillimage[2 * rampwidth:np.int(2.5 * rampwidth),
- ySize // 4:ySize // 4 + ySize // 2] = oneline # first band, going up
- stillimage[np.int(2.5 * rampwidth):3 * rampwidth, ySize // 4:ySize // 4 + ySize // 2] = np.flip(
- oneline) # first band, going up
- # create a bit of a photo:.
- # background is darker
- stillimage[3 * rampwidth:xSize - rampwidth, 0:ySize] = (lightlevel - darknoise_mean) * 0.7
- # large area (antennal lobe) is average
- stillimage[40:95, 15:110] = (lightlevel - darknoise_mean)
- # some glomeruli are brighter, to unequal extent
- thisglo = glo1
- stillimage[thisglo[0]:thisglo[1], thisglo[2]:thisglo[3]] = (lightlevel - darknoise_mean) * 1.3
- thisglo = glo2
- stillimage[thisglo[0]:thisglo[1], thisglo[2]:thisglo[3]] = (lightlevel - darknoise_mean) * 1.5
- thisglo = glo3
- stillimage[thisglo[0]:thisglo[1], thisglo[2]:thisglo[3]] = (lightlevel - darknoise_mean) * 1.1
- thisglo = glo4
- stillimage[thisglo[0]:thisglo[1], thisglo[2]:thisglo[3]] = (lightlevel - darknoise_mean) * 1.8
- thisglo = glo5
- stillimage[thisglo[0]:thisglo[1], thisglo[2]:thisglo[3]] = (lightlevel - darknoise_mean) * 1.3
- thisglo = glo6
- stillimage[thisglo[0]:thisglo[1], thisglo[2]:thisglo[3]] = (lightlevel - darknoise_mean) * 1.3
- if badpixels:
- # bad pixels are fixed in space - therefore create in an image
- noise = np.random.rand(xSize, ySize) > badpixels_clip # values in [0-1]
- stillimage[noise] = 0
- # assumes that from now on everything is multiplicative, not additive
- # check result
- # plt.imshow(stillimage[:,:].T, origin='lower')
- # broadcast still image to time domain
- # movie memory structure
- raw666 = np.zeros((xSize, ySize, tSize))
- stillimage = stillimage.reshape(xSize, ySize, 1)
- raw666 = raw666 + stillimage # use broadcasting
- # now go into the time domain:
- # add bleaching
- if bleach:
- fullframe = (0, xSize, 0, ySize) # apply bleaching to full frame
- oneline = np.arange(tSize)
- oneline = (1 - bleach_percent / 200) + (bleach_percent / 100) * np.exp(-oneline / bleach_tau)
- raw666 = applytimecourse(raw666, fullframe, oneline)
- # plt.plot(oneline)
- if lampnoise: # if lampnoise == 0, it is like false.
- fullframe = (0, xSize, 0, ySize) # apply bleaching to full frame
- oneline = 1 + (lampnoise / 100) * np.random.rand(tSize)
- raw666 = applytimecourse(raw666, fullframe, oneline)
- # plt.plot(oneline)
- # add glomeruli
- ###sawtooth
- oneline = np.zeros(tSize)
- oneline[stimonset:stimoffset] = np.linspace(0, peaksignal / 100, stimoffset - stimonset)
- oneline[stimoffset:] = np.linspace(peaksignal / 100, 0, tSize - stimoffset)
- oneline = oneline + 1 # baseline factor is 1
- raw666 = applytimecourse(raw666, glo1, oneline)
- raw666 = applytimecourse(raw666, glo2, oneline)
- ###realistic:
- tau1 = stimoffset - stimonset
- tau2 = 2 * tau1
- c = 1
- oneline = np.ones(tSize)
- response = np.arange(tSize - stimonset)
- response = (-1) * np.exp(-response / tau1) + c * np.exp(-response / tau2)
- response = ((peaksignal / 100) / np.max(response)) * response # scale to peaksignal
- oneline[stimonset:] = oneline[stimonset:] + response
- raw666 = applytimecourse(raw666, glo3, oneline)
- raw666 = applytimecourse(raw666, glo4, oneline)
- raw666 = applytimecourse(raw666, glo5, oneline)
- # late response
- oneline = np.ones(tSize)
- oneline[stimoffset:] = oneline[stimoffset:] + response[:(tSize - stimoffset)]
- raw666 = applytimecourse(raw666, glo6, oneline)
- # biphasic response
- oneline = np.ones(tSize)
- oneline[stimonset:] = oneline[stimonset:] + response
- oneline[stimoffset:] = oneline[stimoffset:] - 2 * response[:(tSize - stimoffset)]
- raw666 = applytimecourse(raw666, glo7, oneline)
- # negative response, builds on factors defined in realistic
- oneline = np.ones(tSize)
- oneline[stimonset:] = oneline[stimonset:] - response
- raw666 = applytimecourse(raw666, glo8, oneline)
- # plt.plot(oneline)
- # add scattered light, as gaussian filter
- if scatterlight:
- raw666[3 * rampwidth:xSize - 2 * rampwidth, :, :] = sp.ndimage.gaussian_filter(
- raw666[3 * rampwidth:xSize - 2 * rampwidth, :, :], sigma=[scatterlight, scatterlight, 0])
- if shotnoise: # add noise, random and proportional to square root of signal
- # photons are poisson distributed
- # noisier = lambda p : np.random.poisson(p)
- # raw666 = noisier(raw666)
- raw666 = np.random.poisson(raw666) # one line solution! Python is so cool...
- if darknoise:
- # add chip noise.
- # This comes at the end, because it is not influenced by light
- noise = np.random.rand(xSize, ySize, tSize) # values in [0-1]
- noise = noise * darknoise_amplitude + darknoise_mean - darknoise_amplitude / 2
- raw666 = raw666 + noise
- # now, AFTER the noise, get a ramp of totally clean signals for calibration
- # create a step function with stepsize steps, also to 2*lightlevel
- # this area has no signals in time - only differrent background, no noise
- stepfunction = np.repeat(np.arange(stepsize), ySize / stepsize)
- stepfunction = stepfunction * (2 * lightlevel / np.max(stepfunction))
- # shift by darknoise_mean to avoid having values of 0, which cause trouble when dividing in CalcSig
- stepfunction = stepfunction + darknoise_mean
- stepfunction = stepfunction.reshape(1, ySize, 1)
- raw666[-rampwidth:, :, :] = stepfunction # set to uniform background level
- # create a step function with stepsize steps
- stepfunction = np.repeat(np.arange(stepsize), ySize / stepsize) - stepsize // 2 # negative and positive signals
- # calibrate to peaksignal max both ways
- stepfunction = ((peaksignal / 100) / np.max(stepfunction)) * stepfunction
- # values are given in percent
- stepfunction = (stepfunction + 1) * lightlevel
- raw666[-2 * rampwidth:-rampwidth, :, :] = lightlevel # set to uniform background level
- stepfunction = stepfunction.reshape(1, ySize, 1)
- raw666[-2 * rampwidth:-rampwidth, :, stimonset:stimoffset] = stepfunction
- # plt.plot(raw666[-2*rampwidth,:,10])
- # plt.plot(raw666[-2*rampwidth,:,25])
- # plt.show()
- ##
- # camera is 12 bit, therefore clip to 4095, and integer
- raw666 = np.clip(raw666.astype(int), 0, 4095)
- # plt.imshow(raw666[:,:,30].T, origin='lower')
- # plt.show()
- # plt.plot(raw666[100,10,:])
- # plt.plot(raw666[100,60,:])
- # plt.plot(raw666[100,80,:])
- # plt.plot(raw666[80,60,:])
- # plt.plot(raw666[60,60,:])
- # plt.show()
- return raw666 # end of create_raw_data667
- def create_raw_data666(p1_metadata, peaksignal):
- '''
- Giovanni, October 2020, based on previous 666 for square glomeruli, which is now called 667
- Add: glomeruli with round shapes, because CaImAn detects round shapes
- Add: movement (for now just lateral movements)
- Add: frames out of focus (simulated with smoothing)
- Add: unequal bleaching (created by adding non-bleaching base fluorescence)
- -> variable noBleach_lightlevel
- move bad pixels to the end (after artificial movement)
-
- Parameters
- ----------
- p1_metadata : TYPE
- contains all parameters, used are stimulus_on and stimulus_end
- peaksignal : TYPE real
- like odor concentration, i.e. how strong the response in the fictive glomeruli.
- Returns
- -------
- matrix : TYPE
- corresponds to p1.raw, to be copied there by caller
- '''
- print('ViewLoadData/LoadRaw666: generating test data with many fixed settings!')
-
- addRightBand = False
- #right band has fixed value squares, no noise, no change over time.
- addLeftBand = False
- # lateral band left is a continuus up and down, including noise
- stimulus_pulse_start_end_frames = p1_metadata["pulsed_stimuli_handler"].get_pulse_start_end_frames()
- stimonset = stimulus_pulse_start_end_frames[0][0]
- stimoffset = stimulus_pulse_start_end_frames[0][1]
- ### for playing around start after #:
- #####
- #peaksignal, stimonset, stimoffset = 10, 25, 35 # - is in command line now
- # local flags
- lightlevel = 1000 # baseline value of each pixel.
- # at the setup, you would set exposure time and light level accordingly
- noBleach_lightlevel = 500 # add this base lavel after fictive bleaching,
- # the effect is that bleaching cannot be completely corrected by our standard procedure,
- # because it will not have the same parameters in the deltaF/F data across all glomeruli.
- randomseed = 0 # set to None for irreproducible random numbers
- # add chipnoise
- darknoise = True
- darknoise_mean = 100 # mean of background values in the chip
- darknoise_amplitude = 20 # max amplitude of background chip noise
-
- # subtract the constant parts from the light level,
- # because in the experiment we set the light to the sum of them all
- brain_lightlevel = lightlevel - darknoise_mean - noBleach_lightlevel#
- # I assume chip noise to be random over time
- #
- xSize, ySize, tSize = 172, 130, 80 # 120, 99, 80
- rampwidth = 10 # width of the ramp, 3 left, 2 right, i.e.
- # usable range therefore is, in x, 30-152, center is 92
- stepsize = 13 # must be a divisor of ySze
- # glomerular locations for square glomeruli are (make sure to exclude rampwidth)
- # coordinates are x from left to right
- # y from bottom to top
- # glo1 = (50, 70, 90, 110) # x,x - y,y
- # glo2 = (60, 80, 60, 80) # x,x - y,y
- # glo3 = (65, 75, 45, 55)
- # glo4 = (70, 90, 30, 50)
- # glo5 = (50, 65, 30, 45)
- # glo6 = (70, 90, 50, 70) # off response
- # glo7 = (50, 70, 90, 110) # biphasic response
- # glo8 = (75, 95, 90, 110) # negative response
- # glomerular locations for circle glomeruli, syntax centerX, centerY, radius
- # x from left to right, y from bottom to top
- circle1left = (50, 65, 20)
- circle2 = (60, 85, 15)
- circle3 = (70, 97, 10)
- circle4top = (92, 105, 18)
- circle5 = (110, 85, 12)
- circle6right= (120, 65, 15)
- circle7 = (110, 45, 10)
- circle8bottom=(92, 25, 19)
- circle9 = (70, 32, 10)
- circle0 = (60, 50, 16)
- circleCenter= (90, 65, 15)
- shotnoise = True # add photon noise
- # add bad pixels on the CCD chip, i.e. pixels with value 0
- badpixels = True
- badpixels_clip = 0.98 # 0.9 for 10% bad pixels
- #
- # add global bleaching
- bleach = True
- bleach_percent = 10 # percentage bleaching
- bleach_tau = tSize / 5 # assume decay to 37% of percentage over this length
- # add lamp noise
- lampnoise = 0.1 # value in percent. 0 for no noise. Value is PEAK noise
- # add scattered light
- scatterlight = 2 # 2 # unit is gaussian pixels; 0 for no scattered light
- # smoothBackground
- smoothBackground = 5
- # smoooth some frames, to mimic out of focus movement
- smoothOutOfFocus = 5
- smoothOutOfFocus_frames = [4,5,34,35,36,37,50,51,52,53,53,55,70,71,72]
- # add movement
- AddMovement = True
- AddMovementSize = 3 # in pixels, x and y direction, gaussian distribution
- # change code, do exclude lateral step image
- # TODO
- # - add light scatter (smooth): modify
- # - add movement
- #
- def applytimecourse_square_glo(matrix, square_glo, oneline):
- '''
- applies oneline as a factor (length: tSize)
- to square window (glomerulus) defined with upperleft, lowerright
- glo_location is (x_left,x_right,y_down, y_up)
- (y_down < y_up; but safety net included)
- '''
- x1 = np.min((square_glo[0], square_glo[1]))
- x2 = np.max((square_glo[0], square_glo[1]))
- y1 = np.min((square_glo[2], square_glo[2]))
- y2 = np.max((square_glo[3], square_glo[3]))
- matrix[x1:x2, y1:y2, :] = matrix[x1:x2, y1:y2, :] * oneline.reshape(1, 1, tSize)
- return matrix
- # functions to create circular regions, and apply time courses to them
- # def create_circular_mask(h, w, center=None, radius=None):
- # # from
- # # https://stackoverflow.com/questions/44865023/how-can-i-create-a-circular-mask-for-a-numpy-array
- # if center is None: # use the middle of the image
- # center = (int(w/2), int(h/2))
- # if radius is None: # use the smallest distance between the center and image walls
- # radius = min(center[0], center[1], w-center[0], h-center[1])
-
- # Y, X = np.ogrid[:h, :w]
- # dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)
-
- # mask = dist_from_center <= radius
- # return mask
- def create_circular_mask(x, y, glomerulus):
- # glomerulus is (x, y, radius)
- # x from left to right, y from bottom to top
- Y, X = np.ogrid[:x, :y]
- # why x/y are inverted in the glomerulus coordinates (rows/columns instead of x/y) evades me
- dist_from_center = np.sqrt((X - glomerulus[1])**2 + (Y-glomerulus[0])**2)
- mask = dist_from_center <= glomerulus[2]
- return mask # this is a boolean 2D array.
- def applytimecourse_mask(matrix, mask, oneline):
- '''
- applies oneline as a factor (length: tSize)
- to a boolean glomerulus (2D) in mask
- '''
- #uses broadcasting to apply 2D mask
- matrix[mask,] = matrix[mask,] * oneline.reshape(1, 1, tSize)
- return matrix
- def applytimecourse_glomerulus(matrix, glomerulus, oneline):
- mask = create_circular_mask(matrix.shape[0], matrix.shape[1], glomerulus)
- matrix = applytimecourse_mask(matrix, mask, oneline)
- return matrix
- # set random numbers
- np.random.seed(randomseed)
- # create still image
- # with average light level as visible to the experimenter
- # therefore subtract darknoise_mean from lightlevel (has been done avoe)
- stillimage = np.full((xSize, ySize), brain_lightlevel, dtype='float64')
-
- # create a bit of a photo:.
- # background is darker
- stillimage[:, :] = brain_lightlevel * 0.7
- # large circle in the middle (antennal lobe) is average
- circlesize = min(xSize//2, ySize//2)
- position = (xSize//2, ySize//2, circlesize)
- stillimage[create_circular_mask(xSize, ySize, position)] = brain_lightlevel
- # two concentric circles inside, brighter
- position = (xSize//2, ySize//2, circlesize//1.5)
- stillimage[create_circular_mask(xSize, ySize, position)] = brain_lightlevel * 1.5
- position = (xSize//2, ySize//2, circlesize//2)
- stillimage[create_circular_mask(xSize, ySize, position)] = brain_lightlevel * 2
- # some glomeruli have different brightness, to unequal extent
- stillimage[create_circular_mask(xSize, ySize, circle1left)] *= 1.7
- stillimage[create_circular_mask(xSize, ySize, circle2)] *= 0.9
- stillimage[create_circular_mask(xSize, ySize, circle3)] *= 2.0
- stillimage[create_circular_mask(xSize, ySize, circle4top)] *= 0.7
- stillimage[create_circular_mask(xSize, ySize, circle5)] *= 1.4
- stillimage[create_circular_mask(xSize, ySize, circle6right)] *= 0.9
- stillimage[create_circular_mask(xSize, ySize, circle7)] *= 1.5
- stillimage[create_circular_mask(xSize, ySize, circle8bottom)] *= 0.7
- stillimage[create_circular_mask(xSize, ySize, circle9)] *= 1.9
- stillimage[create_circular_mask(xSize, ySize, circle0)] *= 0.8
- stillimage[create_circular_mask(xSize, ySize, circleCenter)] *= 2.0
- # smoooth image, sigma is smoothBackground
- stillimage = sp.ndimage.gaussian_filter(stillimage, smoothBackground, mode='nearest')
- if addLeftBand:
- # left side, overwriting everything here: create four ramps
- # create a ramp from 0 to 2*lightlevel
- oneline = np.linspace(0, 2 * lightlevel, ySize).reshape(1, ySize)
- stillimage[0:rampwidth, :] = oneline # first band, going up
- stillimage[rampwidth:2 * rampwidth, :] = np.flip(oneline) # second band, going down
- # create a ramp in the center, in order to avoid border effects
- oneline = np.linspace(0, 2 * lightlevel, ySize // 2).reshape(1, ySize // 2)
- stillimage[2 * rampwidth:np.int(2.5 * rampwidth),
- ySize // 4 : ySize // 4 + ySize // 2] = oneline # third half band, going up
- stillimage[np.int(2.5 * rampwidth):3 * rampwidth,
- ySize // 4:ySize // 4 + ySize // 2] = np.flip(oneline) # fourth half band, going down
-
- # check result while debugging
- # plt.imshow(stillimage[:,:].T, origin='lower')
- # broadcast still image to time domain
- # movie memory structure
- raw666 = np.zeros((xSize, ySize, tSize))
- stillimage = stillimage.reshape(xSize, ySize, 1)
- raw666 += stillimage # use broadcasting
- # now go into the time domain:
- # add bleaching
- if bleach:
- fullframe = (0, xSize, 0, ySize) # apply bleaching to full frame
- oneline = np.arange(tSize)
- oneline = (1 - bleach_percent / 200) + (bleach_percent / 100) * np.exp(-oneline / bleach_tau)
- # calculate deltaF/f first - avoid values of 0, therefore add 10
- stillimage += 10
- raw666 += 10
- raw666 /= stillimage
- raw666 = applytimecourse_square_glo(raw666, fullframe, oneline)
- # calculate back to "original" data
- raw666 *= stillimage
- raw666 += np.min(raw666) # avoid negative values
- # plt.plot(oneline)
- if lampnoise: # if lampnoise == 0, it is like false.
- fullframe = (0, xSize, 0, ySize) # apply lampnoise to full frame
- oneline = 1 + (lampnoise / 100) * np.random.rand(tSize)
- raw666 = applytimecourse_square_glo(raw666, fullframe, oneline)
- # plt.plot(oneline)
- # add glomeruli, using different time courses
- ###sawtooth
- oneline = np.zeros(tSize)
- oneline[stimonset:stimoffset] = np.linspace(0, peaksignal / 100, stimoffset - stimonset)
- oneline[stimoffset:] = np.linspace(peaksignal / 100, 0, tSize - stimoffset)
- oneline = oneline + 1 # baseline factor is 1
- #raw666 = applytimecourse_square_glo(raw666, glo1, oneline)
- # raw666 = applytimecourse_square_glo(raw666, glo2, oneline)
- raw666 = applytimecourse_glomerulus(raw666, circle1left, oneline)
- ###realistic:
- tau1 = stimoffset - stimonset
- tau2 = 2 * tau1
- c = 1
- oneline = np.ones(tSize)
- response = np.arange(tSize - stimonset)
- response = (-1) * np.exp(-response / tau1) + c * np.exp(-response / tau2)
- response = ((peaksignal / 100) / np.max(response)) * response # scale to peaksignal
- oneline[stimonset:] += response
- # raw666 = applytimecourse_square_glo(raw666, glo3, oneline)
- # raw666 = applytimecourse_square_glo(raw666, glo4, oneline)
- # raw666 = applytimecourse_square_glo(raw666, glo5, oneline)
- # apply to circle glomerulus
- raw666 = applytimecourse_glomerulus(raw666, circle2, oneline)
- oneline = (oneline -1)*1.5 +1
- raw666 = applytimecourse_glomerulus(raw666, circle3, oneline)
- # late response
- oneline = np.ones(tSize)
- oneline[stimoffset:] = oneline[stimoffset:] + response[:(tSize - stimoffset)]
- # raw666 = applytimecourse_square_glo(raw666, glo6, oneline)
- # apply to circle glomerulus
- raw666 = applytimecourse_glomerulus(raw666, circle4top, oneline)
- oneline = (oneline -1)*1.5 +1
- raw666 = applytimecourse_glomerulus(raw666, circle5, oneline)
- # slightly delayed
- oneline = np.ones(tSize)
- oneline[stimonset+5:] += response[:-5]
- raw666 = applytimecourse_glomerulus(raw666, circleCenter, oneline)
- # biphasic response
- oneline = np.ones(tSize)
- oneline[stimonset:] = oneline[stimonset:] + response
- oneline[stimoffset:] = oneline[stimoffset:] - 2 * response[:(tSize - stimoffset)]
- # raw666 = applytimecourse_square_glo(raw666, glo7, oneline)
- # apply to circle glomerulus
- raw666 = applytimecourse_glomerulus(raw666, circle6right, oneline)
- oneline = (oneline -1)*1.5 +1
- raw666 = applytimecourse_glomerulus(raw666, circle7, oneline)
- oneline = (oneline -1)*0.6 +1
- raw666 = applytimecourse_glomerulus(raw666, circle0, oneline)
- # negative response, builds on factors defined in realistic
- oneline = np.ones(tSize)
- oneline[stimonset:] = oneline[stimonset:] - response
- # raw666 = applytimecourse_square_glo(raw666, glo8, oneline)
- # apply to circle glomerulus
- raw666 = applytimecourse_glomerulus(raw666, circle8bottom, oneline)
- oneline = (oneline -1)*1.5 +1
- raw666 = applytimecourse_glomerulus(raw666, circle9, oneline)
- # plt.plot(oneline)
- # add scattered light, as gaussian filter
- if scatterlight:
- raw666[3*rampwidth:xSize - 2*rampwidth, :, :] = sp.ndimage.gaussian_filter(
- raw666[3 * rampwidth:xSize - 2 * rampwidth, :, :], sigma=[scatterlight, scatterlight, 0])
- if shotnoise: # add noise, random and proportional to square root of signal
- # photons are poisson distributed
- # noisier = lambda p : np.random.poisson(p)
- # raw666 = noisier(raw666)
- raw666 = np.random.poisson(raw666) # one line solution! Python is so cool...
- ##movement. I want movement that is subpixel, therefore each frame is rescaled 10fold.
- # using for-loop to save memory space
- if AddMovement:
- # create random movement array in x and y. AddMovementSize
- # movement is slow, therefore random numbers first on array a tenth of the size
- # then blown up by 10, and limited to tSize elements
- x_movement = np.random.rand(tSize//10 + 1)*AddMovementSize
- x_movement = ndimage.zoom(x_movement, 10)[:tSize]
- y_movement = np.random.rand(tSize//10 + 1)*AddMovementSize
- y_movement = ndimage.zoom(y_movement, 10)[:tSize]
- # shift image by image
- for i in range(tSize):
- raw666[3*rampwidth:xSize - 2*rampwidth,:,i] = ndimage.shift(
- raw666[3*rampwidth:xSize - 2*rampwidth,:,i], (x_movement[i], y_movement[i]),
- mode='nearest')
- if darknoise:
- # add chip noise.
- # This comes at the end, because it is not influenced by light
- noise = np.random.rand(xSize, ySize, tSize) # values in [0-1]
- noise = noise * darknoise_amplitude + darknoise_mean - darknoise_amplitude / 2
- raw666 = raw666 + noise
- # add background fluorescence level that is not affected by any of the above manipulations
- # simulating fluorescence that is not from the reporter
- # and that is not affected by bleaching
- raw666 += noBleach_lightlevel
-
- # now add movement. Smooth to mimic vertical movement
- # smoooth some frames, to mimit out of focus movement
- # smoothOutOfFocus = 5
- # smoothOutOfFocus_frames = [4,5,34,35,36,37,50,51,52,53,53,55,70,71,72]
- if smoothOutOfFocus: #only if not 0
- # syntax as for scatterlight
- raw666[3 * rampwidth:xSize - 2 * rampwidth, :, smoothOutOfFocus_frames] = sp.ndimage.gaussian_filter(
- raw666[3 * rampwidth:xSize - 2 * rampwidth, :, smoothOutOfFocus_frames], sigma=[scatterlight, scatterlight, 0])
-
- if badpixels: ##move to after movement creation
- # bad pixels are fixed in space - therefore equal for third dimension
- noise = np.random.rand(xSize, ySize) > badpixels_clip # values in [0-1]
- raw666[noise,:] = 0
- # now, AFTER the noise, get a ramp of totally clean signals for calibration
- if addRightBand:
- # create a step function with stepsize steps, also to 2*lightlevel
- # this area has no signals in time - only differrent background, no noise
- stepfunction = np.repeat(np.arange(stepsize), ySize / stepsize)
- stepfunction = stepfunction * (2 * lightlevel / np.max(stepfunction))
- # shift by darknoise_mean to avoid having values of 0, which cause trouble when dividing in CalcSig
- stepfunction = stepfunction + darknoise_mean
- stepfunction = stepfunction.reshape(1, ySize, 1)
- raw666[-rampwidth:, :, :] = stepfunction # set to uniform background level
-
- # create a step function with stepsize steps
- stepfunction = np.repeat(np.arange(stepsize), ySize / stepsize) - stepsize // 2 # negative and positive signals
- # calibrate to peaksignal max both ways
- stepfunction = ((peaksignal / 100) / np.max(stepfunction)) * stepfunction
- # values are given in percent
- stepfunction = (stepfunction + 1) * lightlevel
- raw666[-2 * rampwidth:-rampwidth, :, :] = lightlevel # set to uniform background level
- stepfunction = stepfunction.reshape(1, ySize, 1)
- raw666[-2 * rampwidth:-rampwidth, :, stimonset:stimoffset] = stepfunction
- # plt.plot(raw666[-2*rampwidth,:,10])
- # plt.plot(raw666[-2*rampwidth,:,25])
- # plt.show()
- ##
- # camera is 12 bit, therefore clip to 4095, and integer
- raw666 = np.clip(raw666.astype(int), 0, 4095)
- # plt.imshow(raw666[:,:,30].T, origin='lower')
- # plt.show()
- # plt.plot(raw666[100,10,:])
- # plt.plot(raw666[100,60,:])
- # plt.plot(raw666[100,80,:])
- # plt.plot(raw666[80,60,:])
- # plt.plot(raw666[60,60,:])
- # plt.show()
- return raw666
- def MovementCorrection(MatrixIN, flag):
- #;input: a 3-dim-matrix MatrixIN
- #;output: another matrix, with each frame shifted in x and y so that the movements are corrected
- #;this output has been removed - program should crash when called in old mode
- #;output: shiftArray contains how much each frame was shifted
- #;options
- def old_RollIt(frames, shiftArray, neighbourSearch, MatrixIN, matrixOUT, correlation, maxShiftX, maxShiftY):
- # uses shiftArray, matrix etc.
- # does not work yet (26.7.2018)
- bestXshift = shiftArray[0,frames]
- bestYshift = shiftArray[1,frames]
- # ;take values from previous shift
- for sx in range( (-1)*neighbourSearch, neighbourSearch+1):
- for sy in range( (-1)*neighbourSearch, neighbourSearch+1):
- shift1x = shiftArray[0,frames] + sx
- shift1y = shiftArray[1,frames] + sy
- if (shift1x > maxShiftX) : shift1x = maxShiftX
- if (shift1y > maxShiftY) : shift1y = maxShiftY
- if (shift1x < -maxShiftX): shift1x = -maxShiftX
- if (shift1y < -maxShiftY): shift1y = -maxShiftY
- result = np.corrcoef(MatrixIN[ (maxShiftX+shift1x):(maxShiftX+smallMsizeX+shift1x),
- (maxShiftY+shift1y):(maxShiftY+smallMsizeY+shift1y),
- frames].flatten(),
- MatrixIN[ (maxShiftX):(maxShiftX+smallMsizeX),
- (maxShiftY):(maxShiftY+smallMsizeY),
- compareFrame].flatten())
- # result = correlate(matrixIN(maxShiftX+shift1x:(maxShiftX+smallMsizeX+shift1x-1),maxShiftY+shift1y:(maxShiftY+smallMsizeY+shift1Y-1),frames), $
- # matrixIN(maxShiftX:(maxShiftX+smallMsizeX-1),MaxShiftY:(MaxShiftY+smallMsizeY-1),compareFrame) )
- #
- if (result[1][0] > correlation): # THEN begin ; found a better shift
- bestXshift = shift1x
- bestYshift = shift1y
- correlation= result[1][0]
- if (bestXshift > maxShiftX): print('*****Movementcorrection: shiftX exceeds positive limit!')
- if (bestXshift < -maxShiftX): print('*****Movementcorrection: shiftX exceeds negative limit!')
- if (bestYshift > maxShiftY): print('*****Movementcorrection: shiftY exceeds positive limit!')
- if (bestYshift < -maxShiftY): print('*****Movementcorrection: shiftY exceeds negative limit!')
- shiftArray[0,frames] = bestXshift
- shiftArray[1,frames] = bestYshift
- shiftArray[2,frames] = int(correlation*10)
- # the IDL program shifted by the negative amount??
- matrixOUT[:,:,frames] = np.roll(matrixOUT[:,:,frames], bestXshift, axis=0)
- matrixOUT[:,:,frames] = np.roll(matrixOUT[:,:,frames], bestYshift, axis=1)
- return (shiftArray, matrixOUT, correlation) #Rollit within MovementCorrection
- def RollIt(frames, compareFrame, neighbourSearch, MatrixIN, MatrixOUT, shiftArray):
- #simplified version - shifts always from central 0 position
- #this solution is slow, because all shifts are evaluated every time
- # in IDL (old_RollIt), neighbour search meant that only few additional shifts were tested each time
- # shift matrix by all pixel values up to neighbourSearch
- correlation = 0
- for sx in range( (-1)*neighbourSearch, neighbourSearch+1):
- for sy in range( (-1)*neighbourSearch, neighbourSearch+1):
- # for each shift, calculate correlation
- shiftMatrix = np.roll( MatrixIN[:,:,frames], sx, axis=0)
- shiftMatrix = np.roll(shiftMatrix, sy, axis=1)
- #calculate correlation, avoiding the border region
- result = np.corrcoef(shiftMatrix[ (neighbourSearch):(-neighbourSearch),
- (neighbourSearch):(-neighbourSearch)].flatten(),
- MatrixIN[ (neighbourSearch):(-neighbourSearch),
- (neighbourSearch):(-neighbourSearch),
- compareFrame].flatten())
- if (result[1][0] > correlation): # THEN begin ; found a better shift
- bestXshift = sx
- bestYshift = sy
- correlation= result[1][0]
- shiftArray[0,frames] = bestXshift
- shiftArray[1,frames] = bestYshift
- shiftArray[2,frames] = int(correlation*10)
- # the IDL program shifted by the negative amount??
- MatrixOUT[:,:,frames] = np.roll( MatrixIN[:,:,frames], bestXshift, axis=0)
- MatrixOUT[:,:,frames] = np.roll(MatrixOUT[:,:,frames], bestYshift, axis=1)
- return (shiftArray, MatrixOUT) #Rollit within MovementCorrection
- # local settings for MovementCorrection
- maxShift = 0.10 #;proportion of maximal shift allowed
- compareFrame = 5 #;which frame to take as standard view (i.e. not moved)
- plotShift = True #; plot the shift result
- logFile = True #;
- # logFileName = flag.STG_OdorReportPath +flag.STG_ReportTag +'_Shift.log'
- # ; if set to 0 then all possible movements are checked
- (sizeX, sizeY, sizeZ) = MatrixIN.shape
- neighbourSearch = min([int(sizeX * maxShift),int(sizeY * maxShift)])
- #;define ShiftArray which contains all shift information
- shiftArray = np.zeros((3, sizeZ), dtype=int) #;dim0:shiftX,dim1:shiftY, dim3:corr*100
- #;define outputmatrix
- MatrixOUT = MatrixIN.copy()
- # ;up from compareframe
- for frames in range(compareFrame+1, sizeZ):
- (shiftArray, MatrixOUT) = RollIt(frames, compareFrame, neighbourSearch, MatrixIN, MatrixOUT, shiftArray)
- # ;down from compareframe
- for frames in range(compareFrame-1, -1, -1):
- (shiftArray, MatrixOUT) = RollIt(frames, compareFrame, neighbourSearch, MatrixIN, MatrixOUT, shiftArray)
- #
- #;suboptimal approach: go through all possibel shifts
- # implemented differently from IDL now, by increasing neighbourSearch to very high values
- ######## done movement correction
- if plotShift:
- # evenly sampled time at 200ms intervals
- t = np.arange(0., sizeZ, 1) #maybe replace 1 with framerate, to have time
- # use p1.frequency or the like
- # red dashes, blue squares and green triangles
- plt.plot(t, shiftArray[0,:], 'r--', t, shiftArray[1,:], 'b--', t, shiftArray[2,:], 'g^')
- plt.show()
- # window, 0
- # ;loadCT, 39
- # range = max([maxshiftX,maxshiftY])
- # plot, shiftArray(0,*), yrange=[(-1)*range,range],color=140
- # oplot, shiftarray(1,*), color=254
- # oplot, shiftarray(2,*), color=255
- if logFile:
- print('ViewLoadData.MovementCorrection: logFile not implemented - still to do')
- #no need for log file, shifts are written in other program
- # openW, 10, logFileName, /APPEND
- # printF, 10, string(9b)+string(9b)+'Followes report shift for measurement ',p1.ex_name,'*',p1.experiment
- # for i=0,sizeZ-1 do begin #
- # printF, 10, strtrim(string(shiftArray(0,i)),2),+string(9b)+strtrim(string(shiftArray(1,i)),2) #
- # endFOR
- # printF, 10, string(9b)+string(9b)+'****************end of************** ',p1.ex_name
- # close, 10
- return (MatrixOUT, shiftArray) #MovementCorrection
- def ReadWriteMovementValues(flagWrite, MovementList, p1, flag):
- #;reads or writes movement values for within one measurement
- #;uses standard file name conventions
- #;write or read data? readFlag now boolean
- readFlag = (flagWrite == 'read')
- #any other string means: write.
- #filename for the list
- MoveFile = os.path.join(flag["STG_OdorInfoPath"],flag["STG_ReportTag"])+'.moveList'
- #what is the current measurement??
- #???inputlabel = fix(inputOdorFile)
- #debug MoveFile
- #MoveFile = 'C:/Users/Giovanni Galizia/Documents/Production/2017_Alja_PostOdor/Filme_data/data/ORN_Glomeruli/al_100312_e.moveList'
- if readFlag:
- # ;read the file
- if os.path.isfile(MoveFile):
- print('ReadWriteMovementValues: open info file for reading movement:', MoveFile)
- # format of movelist:
- # first column is the name of the measurement. For each measurement there are six rows.
- # first three rows relate to odor 0, second three rows to odor 1. Take odor 1 here.
- # the second column is the row label: 0 for x movement, 1 for y movement, 2 for "quality", i.e. 0,1,2 ... repeat this...
- # the third column is the odor flag, generally 0,0,0,1,1,1, ...repeat this...
- # HERE: take row 4 for x movement, and row 5 for y movement, IGNORE the rest for now.
- # movelist is tab-separated
- moveList_df = pd.read_csv(MoveFile, sep='\t', header=None)
- # which is the measurement I am interested in?
- measuNum = p1.messungszahl
- # I only want the rows that have measuNum in the first column
- moveList_df = moveList_df.loc[moveList_df.iloc[:,0] == measuNum]
- # remove the first two columns, and keep only fourth and fifth row
- #moveList_df = moveList_df.iloc[3:5,3:] #take rows 3-5, i.e. second three rows
- moveList_df = moveList_df.iloc[-3:,3:p1.metadata.frames+3] #take last three rows,
- # and only as many columns as frames starting at 3
- #;definition of movement list in the master is:
- #;movementList = intarr(3,p1.frames,p1.odors+1); x/y/quality, frames, data(odor/wavelength)
- # because in IDL there were, generally, two odors, the first fictive, the second not
- # in most moveLists there will be 6 lines, with identical values in the first 3 and secont 3
- # or with 0 in the first 3 lines.
- # therefore, here in Phython I take the LAST 3 lines
- # this way, in the future, I can write 3 lines only.
- # will go wrong if several odors are used (e.g. separate movementlist in Fura? CHECK)
- MovementList = moveList_df.values # report matrix only
- else: #write file
- print('ReadWriteMovementValues: open info file for writing ', MoveFile)
- # add the first three columns to movementList matrix
- # first column: p1.messungszahl
- # second column 0,1,2 (for x/y/quality)
- # third column 0 (for odor - not used in Python, yet)
- out_info = np.array([[p1.messungszahl,0,0],[p1.messungszahl,1,0],[p1.messungszahl,2,0]])
- out_move = pd.DataFrame(np.append(out_info, MovementList, axis=1)) #also convert to dataframe
- if os.path.isfile(MoveFile): #file exists
- print('ReadWriteMovementValues: appending to old file')
- # load previous data
- old_move = pd.read_csv(MoveFile, sep='\t', header=None)
- out_move = pd.concat([old_move,out_move]) #add this one to the existing list
- # now write the new movementList to file.
- # use dataframe format for my lazyness - np.savetxt would be leaner
- # do not add header line (also IDL did not do that any more)
- out_move.to_csv(MoveFile, sep='\t', header=False, index=False)
- #movementList is just the shiftArray, movementList_w has the first three columns added
- return MovementList #
- #end; ReadWriteMovementValues
- def MovementCorrectionMaster(p1, flag):
- #
- #;old movement correction (giovanni, beware, does not work for all datasets)
- #;or new movement correction (mathias, yet to be tested as of April 1st, 2003), with flag setting above 10
- #
- #; recalculate movement (flag set to 1), or recalculate and save the shifts (flag set to 2), or use saved shifts (flag set to 3)
- MovementFlag = flag.CSM_Movement
- ShrinkFactor = flag.LE_ShrinkFaktor #
- # movement values are stored WITHOUT shrink factor
- MathiasCorrection = (MovementFlag > 10)
- if MathiasCorrection: MovementFlag = MovementFlag - 10
- #;either recalculate movement (flag set to 1), or recalculate and save the shifts (flag set to 2), or use saved shifts (flag set to 3)
- if (MovementFlag == 1): print('ViewLoadData.MovementCorrectionMaster: on the spot')
- if (MovementFlag == 2): print('ViewLoadData.MovementCorrectionMaster: on the spot, saving values')
- if (MovementFlag == 3): print('ViewLoadData.MovementCorrectionMaster: using saved values')
- #IF (movementFlag eq 4) THEN print, 'MovementCorrectionMaster: using saved values, being equal for both sets (master/slave)'
- #IF (movementFlag eq 5) THEN print, 'Extracting Shifts from Movement List - ignoring movements'
- #
- #;check movement on monitor?
- # MonitorMovement = 1 # ;(movementFlag eq 2)
- #
- #;define movement value array
- # odors = 1 # or copy p1.odors, if that exists
- # MovementList = np.zeros((3, p1.frames, odors), dtype=int) # I take 2 odors for backwards compatibility
- # MovementList = np.zeros(3, p1.frames, 1, dtype=int) # I take 2 odors for backwards compatibility
- # the two odors are: 0 for control (empty) or fura1, 1 for measurement or fura2
- #intarr(3,p1.frames,p1.odors+1); x/y/quality, frames, data(odor/wavelength)
- MovementList = np.zeros((3, p1.metadata.frames), dtype=int) # I would take 2 odors for backwards compatibility
- ####################
- # Problem here: odors is abolished in this Python setting, but for FURA I don't have a solution yet
- # the thing is: I cannot do a movement correction separately for the two wavelengths, without alignem them,
- # because they might move astray.
- # For new data, movement correction should use new tools anyway.
- #
- #
- # ;go through raw data, do movement correction, for each 'odorset' separatedly (e.g. also for the two wavelengths in FURA)
- if MovementFlag in [1,2]:
- # IF ((MovementFlag ge 1) AND (MovementFlag le 2)) THEN begin
- # for i=0, p1.odors-1 do begin
- # IF (total(raw1[*,*,*,i]) gt 0) THEN begin
- if MathiasCorrection:
- print('MovementCorrectionMaster - MathiasCorrection still to do')
- # MovementCorrectionMathias()
- # ;MovementCorrectionMathias, raw1[*,*,*,i], movementList[*,*,i] ;in ImageALMathias
- # movementList[0:1,*,i] = MovementCorrectionMathias(raw1[*,*,*,i]) ;eine Funktion??
- # ;MovementCorrectionMathias gibt nur x/y, nicht qualität 'raus
- else:
- #MovementCorrection, raw1[*,*,*,i], movementList[*,*,i] # MatrixIN, shiftArray, flag
- print('MovementCorrectionMaster: now calling MovementCorrection - be patient!')
- (MatrixOUT, MovementList) = MovementCorrection(p1.raw1, flag)
- print('MovementCorrectionMaster: now done with MovementCorrection.')
- if MovementFlag in [2]:
- MovementList = ReadWriteMovementValues('write', MovementList, p1, flag)
- if MovementFlag in [3]:
- # ;load movementList for flag set to 3
- # IF ((MovementFlag ge 3) AND (MovementFlag le 4)) THEN begin
- # ReadWriteMovementValues, 'read', movementList
- # print, 'Reading movement information - values are corrected for shrinkraktor! (MovementCorrectionMaster)'
- # IF (ShrinkFactor ge 2) then begin
- # movementList[0,*,*] = fix(movementList[0,*,*] / shrinkfactor)
- # movementList[1,*,*] = fix(movementList[1,*,*] / shrinkfactor)
- # endIF
- # endIF
- MovementList = ReadWriteMovementValues('read', MovementList, p1, flag)
- if ShrinkFactor > 1:
- MovementList[0:2,:] = (MovementList[0:2,:] / ShrinkFactor).astype(int)
- # now apply to data
- MatrixOUT = p1.raw1.copy
- MatrixOUT[:,:,frames] = np.roll(MatrixOUT[:,:,frames], bestXshift, axis=0)
- MatrixOUT[:,:,frames] = np.roll(MatrixOUT[:,:,frames], bestYshift, axis=1)
- if MovementFlag in [4,5]:
- print('MovementCorrectionMaster: flag 4,5 not implemented yet in phython')
- # ;write info from one odor into the other one
- # IF (MovementFlag eq 4) THEN begin
- # ;which one is the odor that has information?
- # IF total(movementList[*,*,0]) eq 0 then begin
- # ;0 is the empty odor
- # movementList[*,*,0] = movementList[*,*,1]
- # endIF else begin
- # ;1 is the empty odor
- # movementList[*,*,1] = movementList[*,*,0]
- # endELSE
- # endIF
- # IF (movementFlag eq 5) THEN begin
- # ReadWriteMovementValues, 'read', movementList
- # ;no movement correction, but extract shifts
- # referenceImage = fix(flag[LE_StartBackground])
- # referenceOdorBuffer = 1
- # p1.shiftX = movementList[0,referenceImage,referenceOdorBuffer]
- # p1.shiftY = movementList[1,referenceImage,referenceOdorBuffer]
- # monitorMovement = 0 ; no need to see the movement on screen
- # print, 'Movement correction OFF, shifts taken from moveList frame ',referenceimage,'; and odorBuffer ',referenceOdorBuffer
- # ENDIF else begin
- # ;apply movementList to raw data
- # for i=0, p1.odors do begin
- # for j=0, p1.frames-1 do begin
- # raw1[*,*,j,i] = shift(raw1[*,*,j,i],movementList[0,j,i],movementList[1,j,i])
- # endFor
- # endFOR
- # endELSE ;movementFlag 5
- #
- # ;check movement with screen output for flag set to 2
- # ;show corrected data
- # IF monitorMovement then begin
- # ;part1: plot shifts and quality
- # window, 2*p1.odors, xsize=250, ysize=250
- # loadCT, 39
- # ;schleife über odors und x/y/qualität
- # rangeTop = max(MovementList(0:1,*,*))
- # rangeBot = min(MovementList(0:1,*,*))
- # plot, movementList[2,*,0], yrange=[rangeBot-1,rangeTop+1],color=255
- # oplot, movementList[1,*,0], color=254
- # oplot, movementList[0,*,0], color=140
- #
- # ;part 2, plot x/y intersection over time
- # for i=0, p1.odors-1 do begin
- # IF (total(raw1[*,*,*,i]) gt 0) THEN begin
- # showMovement, raw1[*,*,*,i], i+p1.odors ;opens one window for each 'odor'
- # endIF
- # endFOR
- # endIF
- #
- #
- #
- #now overwrite p1.raw1
- p1.raw1 = MatrixOUT
- #end; MovementCorrectionMaster
- return p1 # end movementcorrectionmaster
- # def loadDataMaster(flag, p1):
- # # ;correct flag setting
- # #IF (fix(flag[LE_ShrinkFaktor]) eq 0) THEN flag[LE_ShrinkFaktor]=1
- # if flag["LE_ShrinkFaktor"] == 0:
- # print("ViewLoadData.loaddatamaster: flag.le_shrinkfactor == 0, add correction ")
- # #;that has become necessary for pixel display and coordinate calculation
- # #;in the loading routines 0 and 1 were already equivalent
- # #;Giovanni, FEb.05
- # setup = flag["LE_loadExp"]
- # #;flag to clip dark pixels to a minimum value, Gio Sept 2015
- # #;default value 0, i.e. no clipping.
- # # ClipPixels = flag.le_ClipPixels
- # #
- # # ;load experiment
- #
- # # try to read dbb1 value as the absolute path of a VIEW-tiff. If there are any expected exceptions, abort.
- # try:
- # p1 = load_VIEW_tif(flag, p1)
- # setup = 30
- # except FileNotFoundError as fnfe:
- # pass
- # except AssertionError as ae:
- # if str(ae).startswith("Error reading"):
- # pass
- # except IOError as ioe:
- # if str(ioe).endswith("(3 expected)"):
- # pass
- #
- # if (setup == 0):
- # print("ViewLoadData.loaddatamaster: setup is 0, old Berlin system - load not implemented in python yet")
- # # inFile = flag[stg_DataPath]+inputMark
- # # loadexp, InFile, raw1, foto1, dark1, p1
- # # parts = str_sep(inputMark, dirsign) ; separate parts of path+name
- # # info = size(parts) ; get info on variable PARTS, info(1) holds size of array
- # # experiment = parts(info(1)-1) ; name is at last position of array PARTS
- # # p1.experiment = experiment
- # # endIF
- #
- # elif (setup == 3):
- # p1 = loadTILL(flag, p1)
- # # IF setup eq 3 THEN begin
- # # IF fix(flag[VIEW_No4Darray]) THEN begin
- # # loadTILL_3Dim, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # # endIF else begin
- # # loadTILL, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # # endELSE
- # # endIF
- # elif (setup == 4):
- # p1 = loadFURA(flag, p1)
- # # IF setup eq 4 THEN loadFURA, inputMark, measu2, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # # ;setup 4 is compatible with VIEW_no4Darray
- # #
- # # IF setup eq 5 THEN loadFURA_and_CaGr, inputMark, measu2, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # #
- # # IF setup eq 6 THEN loadZseries, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # #
- # # IF setup eq 7 THEN begin
- # # loadTILL, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # # shuffleMultipleLayers, 3 ; 3 for three layers, recorded in one film, to be splitted as if 3 odors
- # # endIF
- # #
- # # ;load till photonics setup, with split image of two wavelength,
- # # ;treat as FURA data
- # # IF setup eq 8 THEN loadMicroImager, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # #
- # #
- # # IF setup eq 13 THEN loadTILLbyte, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # #
- # # IF setup eq 14 THEN loadIDLraw, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # #
- # #
- # # IF setup eq 17 THEN begin
- # # print, 'loading using setup type 17'
- # # loadTILL3Layer, inputMark, measu2, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # # ;loads measurements done on threelayers in 3 separate dbb files, and glues them side by side
- # # ; Inga, Aug. 2014, Konstanz
- # # endIF
- # #
- # # ;setup 9 loads single wavelength data, using the DBB-file specified in the measu2 column of the .lst file (TILL data)
- # # ;usefult to load only the slave, or only the 380nm measurements in FURA data
- # # IF setup eq 9 THEN loadTILL, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 2
- # #
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # IF setup eq 15 THEN loadOlympusR, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # #
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # IF setup eq 20 THEN loadZeiss, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # elif (setup == 20):
- # p1 = loadZeiss(flag, p1)
- # #
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # ;21 reads a .tif file (the .lst may contain any extension, that is ignored and replaced by .tif
- # # IF setup eq 21 THEN Load3DTiff, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # # ;22 reads a _affine.tif
- # # IF setup eq 22 THEN Load3DTiff, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # # ;23 reads a _full.tif
- # # IF setup eq 23 THEN Load3DTiff, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # #
- # #
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # IF setup eq 26 THEN load_3DArray_Raw, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # #
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # IF setup eq 30 THEN load_redshirtimaging, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst', 1
- # #
- # #
- # # ;load format from Olympus Imaging system, Andy March 2010
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # IF setup eq 40 THEN loadOlympusR, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # # ;inputmark can be the file name, or the line tag number in the .lst file .
- # # IF setup eq 41 THEN loadOlympusR_ratio, inputMark, flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.lst'
- # elif setup == 30: # internal tif format
- # pass # should have been read above
- # elif (setup == 665): #generate test data, use odor concentration as max value
- # p1 = LoadRaw666(flag, p1, peaksignal = -10)
- # #generate test data, peaksignal -10
- # elif (setup == 666):
- # p1 = LoadRaw666(flag, p1, peaksignal = p1.odor_nr)
- # # this allows to use the column odor concentration in a list file to generate different signal sizes
- # elif (setup == 667): #generate test data
- # p1 = LoadRaw666(flag, p1, peaksignal = 10)
- # #generate test data, peaksignal 10
- # else:
- #
- # print("ViewLoadData.loaddatamaster: setup load not implemented in python yet")
- # raise NotImplementedError
- # #
- # #;SlaveCorrection ##NOT IMPLEMENTED IN PYTHON - from old mirror system
- # #;if the image of the slave is shifted/rotated/distorted, correct it
- # #IF ((setup eq 4) OR (setup eq 8)) THEN begin
- # # ;the correction parameters are in file ....slaveMod
- # # ;in the directory with the .inf files
- # # ;a file 'all.slaveMod' is active for ALL files in that /listen directory
- # # ;a file 'animal.slaveMod' overwrites 'all.slaveMod' should both be present
- # # modifyFile1 = flag[stg_OdorInfoPath]+flag[stg_ReportTag]+'.slaveMod'
- # # modifyFile2 = flag[stg_OdorInfoPath]+'all'+'.slaveMod'
- # # existFile1 = existFile(modifyFile1)
- # # existFile2 = existFile(modifyFile2)
- # # IF (existFile1 OR existFile2) THEN begin ;modify if there is the necessary information
- # # ;read that information
- # # IF existFile1 THEN modifyFile = modifyFile1 ELSE modifyFile = modifyFile2
- # # openR, unit, modifyFile, /GET_LUN
- # # line = ''
- # # readF, unit, line
- # # readF, unit, line ;the second line contains the modify information
- # # free_lun, unit ;close the file
- # # print, 'LoadDataMaster: Modifying SLAVE according to information in file: ',modifyFile
- # # print, 'Modify parameters are: ', line
- # # ;line is: shiftX shiftY mod rotAng magni rotCentX rotCentY
- # # shiftArray = str_sep(line, string(9b))
- # # ;the slave to be modified is the dataset raw(*,*,*,1)
- # # ;first the shift, then the rotation
- # # for i=0, p1.frames-1 do begin
- # # b = float(reform(raw1(*,*,i,1))) ; take that frame, remove spurious dimensions
- # # b = shift(b,shiftArray(0),shiftArray(1))
- # # ;modify this slice - after shift
- # # b = rot(b, shiftArray(3),shiftArray(4)/100.0,shiftArray(5),shiftArray(6), /interp, /pivot)
- # # ;Result = ROT( A, Angle, [Mag, X0, Y0] [, /INTERP] [, CUBIC=value{-1 to 0}] [, MISSING=value] [, /PIVOT] )
- # # ;rotate around point x0/y0 by Angle clockwise, Magnify by Mag
- # # raw1(*,*,i,1) = fix(b)
- # # endFor
- # # endIF
- # #endIF
- # #
- # #;movementCorrection
- # #;the settings in movementCorrection affect the shift correction - therefore, movement correction preceeds shift correction
- # if (flag["CSM_Movement"] > 0):
- # p1 = MovementCorrectionMaster(p1, flag)
- # #endIF ; movementCorrection
- # #
- # #
- # #;shift correction: if flag is 0: no shift here, but coordinates are shifted
- # #; if flag is 1: shift here , coordinates are not shifted
- # #; if flag is 2: shift entry in p1.shift is ignored altogether (e.g. when movement correction contains shift information)
- # #$ DataShift not translated yet
- # #IF (fix(flag[CSM_DataShift]) eq 1) THEN begin
- # # print, 'LoadDataMaster: shift Data by ',p1.shiftX, p1.shiftY
- # # ;correct shifts
- # # p1.shiftX = (-1)*p1.shiftX / fix(flag[LE_ShrinkFaktor])
- # # p1.shiftY = (-1)*p1.shiftY / fix(flag[LE_ShrinkFaktor])
- # # ;shift the data instead of shifting the masks
- # # raw1 = shift(raw1,p1.shiftX,p1.shiftY,0,0)
- # # ;put the shifts to 0, to avoid that the masks will be shifted
- # # p1.shiftX = 0
- # # p1.shiftY = 0
- # #endIF
- # #IF (fix(flag[CSM_DataShift]) eq 2) THEN begin
- # # print, 'LoadDataMaster: shift information ignored'
- # # ;put the shifts to 0, to avoid that the masks will be shifted
- # # p1.shiftX = 0
- # # p1.shiftY = 0
- # #endIF
- # #
- # #
- # #
- # #IF (fix(flag[VIEW_No4Darray]) eq 1) THEN begin
- # # IF ((setup eq 4) OR (setup eq 3) OR (setup eq 8)) THEN begin
- # # ;for these options specific care has been taken
- # # end else begin
- # # ;remove the odor buffer 0
- # # ;watchout with FURA data
- # # raw1 = raw1(*,*,*,1:p1.odors)
- # # p1.odors = p1.odors-1
- # # odor = (odor - 1)>0
- # # IF (fix(flag[LE_UseFirstBuffer]) ge 1) THEN flag[LE_UseFirstBuffer] = flag[LE_UseFirstBuffer] - 1
- # # print, '*****************LoadDataMaster: deleted odor buffer 1, only 0. '
- # # endELSE
- # #endIF
- # #
- # ## skip NOT IMPLEMENTED in python (yet?)
- # #;skip frames if requested
- # #SF_upFront = fix(flag(CSM_SkipFrmUpFront))
- # #SF_atBack = fix(flag(CSM_SkipFrmAtBack))
- # #SkipFrames = (SF_upFront ne 0) OR (SF_atBack ne 0)
- # #IF SkipFrames THEN begin
- # # ;if the numbers are negative, they are relative to Stimulus ONSET, otherwise to FIRST FRAME,
- # # IF SF_upFront lt 0 THEN begin
- # # SK_firstFrame = p1.stimulus_on + SF_upFront
- # # endIF else begin
- # # SK_firstFrame = SF_upFront
- # # endELSE
- # # ;if the numbers are negative, they are relative to Stimulus ONSET, otherwise to FIRST FRAME,
- # # IF SF_atBack lt 0 THEN begin
- # # SK_lastFrame = p1.stimulus_on - SF_atBack
- # # endIF else begin
- # # SK_lastFrame = SF_atBack
- # # endELSE
- # # IF (SF_atBack eq 0) THEN begin
- # # SK_lastFrame = p1.frames -1
- # # endIF
- # # ;error messages
- # # IF (SK_firstFrame lt 0) then begin
- # # print, 'skipping frames with wrong parameters in LoadDataMaster.pro: SK_firstFrame = ',SK_firstframe
- # # stop
- # # endIF
- # # IF (SK_lastFrame ge p1.frames) then begin
- # # print, 'skipping frames with wrong parameters in LoadDataMaster.pro: SK_lastFrame = ',SK_lastFrame
- # # stop
- # # endIF
- # # SK_Frames = SK_lastFrame - SK_firstFrame + 1
- # # IF SK_Frames lt 1 THEN begin
- # # print, 'skipping frames with wrong parameters in LoadDataMaster.pro: SK_Frames = ',SK_Frames
- # # stop
- # # endIF
- # # ;now cut the excess frames
- # # raw1 = raw1(*,*,SK_firstFrame:SK_lastFrame,*)
- # # print, 'Frames skipped. Kept frames from ',SK_firstFrame,' to ',SK_lastFrame
- # # ;correct parameters
- # # p1.frames = SK_Frames
- # # p1.stimulus_on = P1.stimulus_on - SK_firstFrame
- # # p1.stimulus_end = p1.stimulus_end - SK_firstFrame
- # # p1.stimulus_del = p1.stimulus_del - SK_firstFrame*p1.trial_ticks
- # # p1.stim2ON = p1.stim2ON - SK_firstFrame
- # # p1.stim2OFF = p1.stim2OFF - SK_firstFrame
- # # p1.stim2_del = p1.stim2_del - SK_firstFrame*p1.trial_ticks
- # #endIF ;SkipFrames
- # #
- # #IF ClipPixels gt 0 THEN begin
- # # print, 'LoadDataMaster: Clip data to ', ClipPixels
- # # raw1 = raw1 > ClipPixels
- # #endIF
- # #
- # #IF ClipPixels lt 0 THEN begin
- # # ;clips pixels to the positive value, and creates a mask file that can be used for overviews
- # # ;example of how to run: set to -300, all pixels <300 will be set to 300, and a mask for the remaining pixels is saved
- # # ;for faster performance, the second time the flag can be set to 300, because the mask file already exists
- # # ClipPixels = (-1)*ClipPixels
- # # print, 'LoadDataMaster: Clip data to ', ClipPixels
- # # raw1 = raw1 > ClipPixels
- # # ;create mask
- # # IF (fix(flag[VIEW_No4Darray]) eq 1) THEN dataPlace = 0 ELSE dataPlace = 1
- # # copyData = reform(raw1(*,*,*,dataPlace)) ;take the measurement
- # # copyData = copyData gt ClipPixels ;gives 1 for all pixels gt clip
- # # ;because we are not working on single frames, average for the time course
- # # copyData = total(copyData, 3);sums across frames, now floating point
- # # copyData = copyData / p1.frames;
- # # copyData = copyData gt 0.5 ; take all pixels where for at least half of the frames
- # # ; the value was larger than ClipPixels
- # # maskFrame = byte(copyData) ;convert to byte
- # # ;no save this mask to file
- # # ;file name taken from singleOverview.pro
- # # areaFileName = flag[stg_OdorMaskPath]+flag[stg_ReportTag]+'.Area'
- # # print, 'LoadDataMaster.pro: saving area file to ',areaFileName
- # # save, filename=areaFileName, maskFrame
- # # copyData = 0B ; empty memory space
- # # maskFrame = 0B
- # #endIF ; ClipPixels lt 0
- # #
- # #
- # data_sampling_period = pd.Timedelta(f"{p1.trial_ticks}ms")
- # p1.pulsed_stimuli_handler.initialize_stimulus_offset(flag["mv_correctStimulusOnset"], data_sampling_period)
- #
- # print('ViewLoadData.LoadDataMaster: Data loaded, setup type ',setup)
- # return # no return any more. p1
|