help scripting...

[SylvainT]

Hi all,

I'm trying to script some little things, and I need some help to go a bit further...

For now, my script aims at opening a T1, a mesh, do a CutMesh fusion. I'd like to control, through my script, the orientation of the cutting plane in that fusion (by default it seems to be axial): first, I'd like to change it from axial to coronal for example; second I'd like to make it oblique, and to control quantitatively the angle. How can I do this?

-------------
# load image and mesh
croppedIm = a.loadObject( ... );
brainMesh = a.loadObject( ... )

# perform CutMesh fusion
cutMeshFusion = a.fusionObjects([brainMesh, croppedIm], "FusionCutMeshMethod")

# view the fusion in a 3D window
w3 = a.createWindow("3D")
w3.addObjects(cutMeshFusion)
-------------

Then, I go on in this script by adding the same image into that same window (w3.addObjects(croppedIm))... Similar wish: by default, the image show up in axial orientation. First, I'd like to change it to another orientation (coronal or sagital). Second, I'd like to make it oblique with a given angle...

Thanks in advance for your help,

Sylvain

ps: I looked in the reference of the anatomist api, but without success... if you can point me to where I can find these kind of info on my own, that'd be great

[riviere]

Hi Sylvain,
You can use the SliceParams command to set the cut mesh plane. You may specify the plane either using a quaternion (rotation from the axial plane), or more simply by giving the plane equation coefficients:

Code: Select all

a.execute( 'SliceParams', objects=[cutMeshFusion], plane=[-0.5, -0.5, 0, 100] )

To change the orientation of a view, you can use the Camera command:

Code: Select all

a.execute( 'Camera', windows=[w3], view_quaternion=[0.5277269, 0.0827091999999, 0.009927049999, 0.845319] )

or change the slice orientation using the same command:

Code: Select all

a.execute( 'Camera', windows=[w3], slice_quaternion=[0.5277269, 0.0827091999999, 0.009927049999, 0.845319] )

You can get the current view orientations using the getInfos() method:

Code: Select all

>>> i = w3.getInfos()
>>> print i
{'boundingbox_max': [152.773, 143.08699999999999, 171.98400000000001],
 'boundingbox_min': [-152.773, -143.08699999999999, -171.98400000000001],
 'geometry': [868, 481, 455, 463],
 'group': 0,
 'objects': [2, 4],
 'observer_position': [129.74299999999999,
                       129.74299999999999,
                       84.825000000000003],
 'position': [119.818, 115.765, 87.524000000000001, 0],
 'referential': 1,
 'slice_quaternion': [0.52772699999999995,
                      0.082709199999999997,
                      0.0099270499999999998,
                      0.84531900000000004],
 'type': 'AWindow',
 'view_quaternion': [0.53561000000000003,
                     0.0154371,
                     -0.018251,
                     0.84412699999999996],
 'view_size': [384, 386],
 'windowType': '3D',
 'zoom': 1}

It's useful to copy/paste quaternion values (because I admit interpretation of those quaternions are not easy at first sight).

Denis

[SylvainT]

Hi Denis, thanks for your answer!

Hi Sylvain,
You can use the SliceParams command to set the cut mesh plane. You may specify the plane either using a quaternion (rotation from the axial plane), or more simply by giving the plane equation coefficients:

Code: Select all

a.execute( 'SliceParams', objects=[cutMeshFusion], plane=[-0.5, -0.5, 0, 100] )

The SliceParams command works, thanks!

To change the orientation of a view, you can use the Camera command:

Code: Select all

a.execute( 'Camera', windows=[w3], view_quaternion=[0.5277269, 0.0827091999999, 0.009927049999, 0.845319] )

or change the slice orientation using the same command:

Code: Select all

a.execute( 'Camera', windows=[w3], slice_quaternion=[0.5277269, 0.0827091999999, 0.009927049999, 0.845319] )

The Camera command does not work (neither with view_quaternion or slice_quaternion, or none of the other keywords I get from getInfos() )... I get the following output message (which is still cryptic to me coz I'm not yet very familiar with python... it doesn't seem to be an error message per say... I dunno what that means...), and nothing changes in the display:

Code: Select all

In [109]: a.execute( 'Camera', windows=[w3], slice_quaternion=[0.5277269, 0.0827091999999, 0.009927049999, 0.845319] )
Out[109]: <anatomist.cpp.RegularCommand object at 0x94e53ec>

Cheers,

Sylvain

[Dominique Geffroy]

Hi,

Indeed, there was a bug with the Camera command in the 3.2.0 version. It is corrected in the 3.2.1 version.

Dominique

[SylvainT]

Hi Dominique, thanks, that solved it

I'm a bit lost in the definition of the slice_ or view_ quaternions of the camera command, but I'm working on it...

I found another problem, with the sliceParams command of my cut mesh fusion this time...
I can define a sagital ( plane=[1, 0, 0, -48] ) or coronal ( plane=[0, 1, 0, -48] ) cutting plane, but when I try an axial one (with tons of attemps: plane=[0, 0, -1, 48], plane=[0, 0, 1, -48], plane=[0, 0, -1, 20], plane=[0, 0, 1, -20], plane=[0, 0, 1, 20] etc.), I get error messages (several hundreds of them):

Code: Select all

planarMeshSortPoints: Open polygon
.. closing open polygon
planarMeshSortPoints: Open polygon
.. closing open polygon
planarMeshSortPoints: Open polygon
.. closing open polygon
etc.

and a unique error message at the end:

Code: Select all

duplicate points cleared

and everything's frozen after that! (I have to kill the ipython process to be able to escape...)

Any clue? If needed, I can send you the files...

Cheers,

Sylvain

[SylvainT]

...hmmmmmmmmm, actually, after having killed the process the first times it happened, I didn't kill it the last time... and a few minutes later, it unfreezed and said:

Code: Select all

12126 points processed
5286 vertices / 1762 polygons in planar mesh

so it must have been computing something... the result is totally blank on the display of my 3D window (as if my cutting plane doesn't intersect the mesh, which is totally possible coz I didn't adjust the values); still, it has to be some kind of bug, coz, as I said in my previous message, this doesn't happen for coronal and sagittal planes (even when the defined plane doesn't intersect the mesh)

[riviere]

Hi Sylvain,
The algorithm that re-meshes the planar polygon of the cut mesh has some failures when the mesh to cut has some problems, especially when it contains crossing polygons. In that case the algorithm normally just stops and you end up with an unfinished polygon (the volume texture is not complete), sometimes it is completely empty. However it generally doesn't take longer to compute, so I don't know why it takes so long. Did you check the CPU usage (using the top comand or a graphical perfmeter) ?
Denis

[SylvainT]

Hi Denis,
I just tried it again to check the CPU usage, and it did stay at 100% during those several minutes, so it seems that it's actually computing something...
Sylvain

[riviere]

Really ? Then it looks like a real bug... I can have a look at your data to see if I can reproduce the problem, yes.
Denis

[riviere]

Hi Sylvain,
I had a look at your data, and of course you are right: the cut mesh algorithm hangs. I know there are problems on meshes which contain intersecting polygons, but I thought they were detected and the algorithm aborts when it detects such problems. But apparently all cases are not detected...
When I developped the plane remeshing algorithm, I had ideas on how to handle such cases, but I never finished it (since it worked on "good" meshes), and now I don't remember the details of it - it would take me several days to get into it again and it would require calm and concentration, which I almost never have in our lab, so I can't promise I will look at it soon... sorry...
Denis

[SylvainT]

Hello Denis, thanks for looking at the data... It's not really the answer I was hoping for, but I understand the situation

Here are two extra ideas/questions that might help me at this point:

- is there a way to "fix" a mesh with intersecting polygons? some resampling could do, for example, even at a lower resolution (smaller number of nodes) since I do not need a super high precision on this project...

- do you have a list of the tools known to potentially create such "bad" meshes, so that I can avoid using them to create mine? or know to insure good properties for the resulting mesh, so that I only use those?

Cheers,

Sylvain

[riviere]

Hi Sylvain,

The intersecting polygons problem is due (I think) to the decimation/smoothing step of the mesher, which does not check this kind of overlaps when "pushing" vertices.
I haven't investigated a lot on how to solve the problem: I have tried once a silly "brute force" exaustive test which had prohibitive execution times, so I gave up at that moment and didn't think about it since then... it should be possible to optimize it using kdtrees for instance, but it is a bit of work in any way. What should be done is to detect every triangles intersection, and split such triangles so that the intersection is materialized by an edge of the mesh. This would (normally) prevent the "cut mesh" algorithm in Anatomist to fail, but would not make the meshes nicer.
- "Resampling" to a lower number of nodes/polygons might actually worsen the problem because vertices would go easily further from their original location.
- Lowering of disabling the decimation/smoothing may help, inversely. I am not sure if completely disabling it would guarantee a good behaviour or not - I have not tried, to be true.

Denis

[SylvainT]

- Lowering of disabling the decimation/smoothing may help, inversely. I am not sure if completely disabling it would guarantee a good behaviour or not - I have not tried, to be true.

OK, I'll try to look into this... Thanks Denis!

[riviere]

Hi Sylvain,

Just to inform you that I have finally completely changed the mesh cut remeshing algorithm, it will be replaced using an OpenGL (GLU) tesselation routine, which takes the intersecting polygons problem into account. It also seems to work faster...
This will be in the next version 4.1 of Anatomist (not in 4.0.2 that we released yesterday).

Denis

[SylvainT]

excellent, looking forward to it so that I can try to use it in my code!!

if you want me to beta-test, let me know!

cheers,

sylvain