Extracting a mask after doing operations with two fields

Hello!
I am running simulations to estimate electric fields induced with noninvasive 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!

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,
ArjamaExample
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 fieldsif 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,Efielddef _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 ggtissues = []# 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) 
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\sitepackages\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\sitepackages\s4l_v1_api\analysiswrappers.py", line 647, in Connect
raise Exception("Can't connect ports")
Exception: Can't connect portsSo 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!