Extracting a mask after doing operations with two fields

[pablomac]

Hello!

I am running simulations to estimate electric fields induced with non-invasive brain stimulation. When using the GUI, I can obtain the field resulting from the interference of two electric fields using the TI Max Modulation function. However, this function is no available on the Python API, so I had to implement my own. The implication of this is that, when I try to extract the fields for a certain part of the brain using the FieldMaskingFilter, I cannot because the fields I calculate are not connected to the necessary structure.

So for example, if I wanted to extract a mask for the simulation of one electrode pair, I would do:
sim_ext = sim.Results()
em_sensor_extractor = sim_ext["Overall Field"]
em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All"
em_sensor_extractor.Normalization.Normalize = True
em_sensor_extractor.Normalization.NewReferenceValue = 0.001
s4l_document.AllAlgorithms.Add(em_sensor_extractor)

inputs = em_sensor_extractor.Outputs["EM E(x,y,z,f0)"]
--Define the mask--
field_masking_filter = analysis.core.FieldMaskingFilter(inputs=inputs)
field_masking_filter.UpdateAttributes()

And then just set materials to False or True depending on the region of interest. What I need is to extract the mask after performing some calculations, so I would like something like:
TImax = TI(EF1, EF2)
--Define the mask--
field_masking_filter = analysis.core.FieldMaskingFilter(inputs=TImax)
field_masking_filter.UpdateAttributes()

So the problem is that I don't know how to create a data structure using my calculated TImax vector.

Could someone please help me out?

Thanks a lot!

[halder]

Hi,
Currently, TI MAM is not exposed in this API. I would recommend getting the Efield sensors, write a function to calculate your TImax. Then create another function where you find entities corresponding to tissue names you want to mask and apply mask filter to calculated input. Please let me know if you have further questions regarding this.

Thank you,
Arjama

Example

def modulation_envelope(e_field_1, e_field_2, dir_vector=None):
E_plus = np.add(e_field_1, e_field_2) # Create the sum of the E fields

if dir_vector is None:
    envelope = np.zeros(e_field_1.shape[0])

    ## Calculate the angles between the two fields for each vector
    dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
    cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
    angles = np.arctan2(cross_angle, dot_angle)
    
    # Flip the direction of the electric field if the angle between the two is greater or equal to 90 degrees
    e_field_2 = np.where(np.broadcast_to(angles >= np.pi/2., (3, e_field_2.shape[0])).T, -e_field_2, e_field_2)
    
    # Recalculate the angles
    dot_angle = np.einsum('ij,ij->i', e_field_1, e_field_2)
    cross_angle = np.linalg.norm(np.cross(e_field_1, e_field_2), axis=1)
    angles = np.arctan2(cross_angle, dot_angle)

    E_minus = np.subtract(e_field_1, e_field_2) # Create the difference of the E fields
    
    # Condition to have two times the E2 field amplitude
    max_condition = np.linalg.norm(e_field_2, axis=1) < np.linalg.norm(e_field_1, axis=1)*np.cos(angles)
    e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), max_condition, False)

    # Condition to have two times the E1 field amplitude
    max_condition_2 = np.linalg.norm(e_field_1, axis=1) < np.linalg.norm(e_field_2, axis=1)*np.cos(angles)
    e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), max_condition_2, False)
    
    # Double magnitudes
    envelope = np.where(e1_gr_e2, 2.0*np.linalg.norm(e_field_2, axis=1), envelope) # 2E2 (First case)
    envelope = np.where(e2_gr_e1, 2.0*np.linalg.norm(e_field_1, axis=1), envelope) # 2E1 (Second case)
    
    # Cross product
    # Redifine e1_gr_e2 to match the area ouside where the max_condition is met
    e1_gr_e2 = np.where(np.linalg.norm(e_field_1, axis=1) > np.linalg.norm(e_field_2, axis=1), np.logical_not(max_condition), False)
    envelope = np.where(e1_gr_e2, 2.0*(np.linalg.norm(np.cross(e_field_2, E_minus), axis=1)/np.linalg.norm(E_minus, axis=1)), envelope) # (First case)
    
    e2_gr_e1 = np.where(np.linalg.norm(e_field_2, axis=1) > np.linalg.norm(e_field_1, axis=1), np.logical_not(max_condition_2), False)
    envelope = np.where(e2_gr_e1, 2.0*(np.linalg.norm(np.cross(-e_field_1, -E_minus), axis=1)/np.linalg.norm(-E_minus, axis=1)), envelope) # (Second case)
else:
    envelope = np.abs(np.abs(np.dot(E_plus, dir_vector)) - np.abs(np.dot(E_minus, dir_vector)))
return envelope

def getEfield(file):
electrode = np.loadtxt(file)
xyz = np.array([electrode[:,0],electrode[:,1],electrode[:,2]]).T
Efield = np.array([electrode[:,3],electrode[:,5],electrode[:,7]]).T
return xyz,Efield

def _findEntities(lm): # finds entities corresponding to tissue names, so I can set them in masking filter
gg=[]
for ii in model.AllEntities():
if(ii.Name in lm):
gg.append(ii)
return gg

tissues = []# created list of tissues I want as first masking filter

Adding a new FieldMaskingFilter

inputs = [TIMAM]
field_masking_filter = analysis.core.FieldMaskingFilter(inputs=inputs)

how to uncheck the "invalidate masked values box"

field_masking_filter.UseNaN = True

clearing selected entities

field_masking_filter.SetAllMaterials(False)
g2=_findEntities(tissues) # finding the entities based on material names in my list above
field_masking_filter.SetEntities(g2)# enableing the entities in first masking filterfield_masking_filter.UpdateAttributes()
document.AllAlgorithms.Add(field_masking_filter)

[pablomac]

Dear Arjama,

Thanks a lot for your reply and for your very detailed example. The problem I was referring to is in the line where you add TIMAM to the mask.
If I run:
inputs = [TIMAM]
field_masking_filter = analysis.core.FieldMaskingFilter(inputs=inputs)

where TIMAM = modulation_envelope(e_field_1, e_field_2, dir_vector=None)

I get the error:
field_masking_filter = analysis.core.FieldMaskingFilter(inputs=[TImax])
File "C:\Program Files\Sim4Life_8.0.0.15165\Python\lib\site-packages\s4l_v1_api\analysiswrappers.py", line 231, in init
self.Inputs[id].Connect(input)
File "C:\Program Files\Sim4Life_8.0.0.15165\Python\lib\site-packages\s4l_v1_api\analysiswrappers.py", line 647, in Connect
raise Exception("Can't connect ports")
Exception: Can't connect ports

So I think my problem is the creation of TIMAM. I Imagine it is not the TI envelope, but rather the TI Envelope inside some data structure that can be connected to the masking filter?

Thanks again!