Code optimization: ray tracing through the 3D grid.

[axpulkki]

Below MatLab profile report for the two main functions. It is very clear that the ray tracing through the volume and associated computations are by far the most significant target for our optimization. We need to implement more intelligent ray tracing algorithm. The algorithm needs to provide the grid points that the ray travels through. One option would be to use https://www.mathworks.com/matlabcentral/fileexchange/26852-a-fast-voxel-traversal-algorithm-for-ray-tracing.

[nasa902neel]

L 75 - L85: map_LOS_2_G_data .... I thought this is location of time crunch.. I will look in more detail

[nasa902neel]

The bresenham algo is also probably what you want:

https://www.mathworks.com/matlabcentral/fileexchange/21057-3d-bresenham-s-line-generation?focused=5102923&tab=function

[axpulkki]

Could you find out which one is faster/better for us? The routine https://www.mathworks.com/matlabcentral/fileexchange/26852-a-fast-voxel-traversal-algorithm-for-ray-tracing or Bresenham? I can then implement the recommended ray racing method in our software.

[nasa902neel]

I am trying to figure out how to implement either into your code.... but I am still struggling to fully understand the meaning of each variable.

Also, I think insignificant difference between the two matlab codes. They look similar.

[axpulkki]

OK. If you hit the wall let me know and I will do the implementation over the weekend.

[nasa902neel]

I thought i would point out that my current NPS_eg_bresen code does not actually implemented the bresenham algo... it actually runs an algo proposed by by Mller and Trumbore (1997) - through a subroutine called TraingleRayIntersection. This subroutine inputs start location and direction of ray.

I think this code is a little more robust...

I have uploaded a ghost routine called Bresenham_line3d... some online notes suggested this might be less robust, so I am currently not using it. I kept the routine as a reminder that we have another option available to test for speed. This subroutine inputs start and end location of ray.

[axpulkki]

OK. So lets go with you NPS_eg_bresen right now. But I will try modify the code so that we can easily swap between different ray tracing algorithms if necessary. So I will essentially use "the standardized" ray tracing call you had in your routine (i.e. IN: start_point, end_point, grid_definition. OUT: grid_indices).

[nasa902neel]

Yes, simply use [OUT: grid_indices] = NPS_eg_bresen [IN: start_point, end_point, grid_definition].

Ignore the ghost code. We can discuss in future if we want to use it... It will be a simple adjustment to my NPS code.

[nasa902neel]

I will be working the late afternoon of sat and most sunday. Sat morning n afternoon text me if needed.. (first is tennis then the science march along the mall).

[axpulkki]

km2rs and related conversion routines used by bresen-code appear to be missing from the repository.

[axpulkki]

But I don't think those conversions are actually needed in the bresen?

[axpulkki]

Neel, I believe your bresen routine fails to run in case the ray does not pierce the cube at all. Could you modify the routine to take into account rays that possibly do not pierce the cube. Also, could we modify the I/O of the function as <IN: start_point, direction, grid_definition. OUT: grid_indices> where direction is vector giving the direction of the ray starting from start_point? I am now implementing our tomography routines with the Amanatide & Woo ray tracing. After those mods inserting also your algorithm in the mix should be plug and play.

[nasa902neel]

Sorry, I have now added the utility routines. I am starting to make adjustment to the nps_eg routine so that the outputcorrectly handles a blank array for inputs that do not go through cube.

[axpulkki]

Learned through profiling that the MatLab sub2ind is horribly slow. Better use your version of the conversion Neel. I have that now in Amanatide & Woo ray tracing.

[nasa902neel]

yeah, I read through the sub2ind code.... it is good as a generalised code for any dimensions. Mine works only for 3D cubes.... so the equation is a single line and quick.

But to be fair, I did not see anything in the sub2ind to jump out as horribly slow, it might just be the parsing of large data cubes?!?

[axpulkki]

But it was shit slow... Another function I am trying to optimize is setdiff. Once I am done with that I think the easy optimizations are mostly completed. Code is already shitloads faster.

[nasa902neel]

great!! ... I could not properly test... as the program crashed for me after 40hours of running.

[axpulkki]

I have now uploaded the routines with new ray tracing implemented. I am performing speed test to check exactly how much improvement we achieved. I will work next restructuring the Tomography_Using_ART.

[nasa902neel]

r u in the office already? I am here.. and will walk to your desk if you are here?

[axpulkki]

I am still home... Should be in the office by 10am.

[axpulkki]

The setdiff (even with LIGHT_setdiff) is still taking perhaps too much time. It might be worthwhile to see how that part of Tomography_Using_ART could be perhaps restructured. We do want to make sure we are not spreading electrons twice to the same grid points.

[nasa902neel]

ok.

[nasa902neel]

I am dnlding your new code now and wil begin testing my stuff with it.

[axpulkki]

For the test, old code ran in 1447 s and the new one in 37 s. So about factor of 40 improvement.

[nasa902neel]

WOW!.. so the bresenham algo's ready did help :) ... I noticed that you simply rewrote everything yourself... so I thought I would go back to my code and update things and see if we get same results independently.

[axpulkki]

Yes, I decided to write everything. But the code should be structured now so that we can easily plug in also other ray tracing algorithms (such as Bresenham) to test which is the best fit for us.

[axpulkki]

Another benchmark: I ran the current sheet computation with FOV 2-9 and 50 arsec resolution coronagraph through 85x61x85 grid and 3 ART iteration rounds over 6 views. The whole thing took 50 minutes. So we can definitely start looking at realistic scenarios with reasonable computation times.