problem in LIME raytracer

[migueldvb]

In our research group we have identified a problem in the LIME raytracer
causing the central pixel fluxes in a ~1/r^2 distribution to be overestimated.
For the same model with different pixel sizes but the same physical scale
the fluxes inside the central 0.2'' square region differ considerably (calculated
by the sum statistic): the image with finer pixel sampling has 3x less flux
than the other one. This could explain the discrepancy we're seeing in the
visibility curves when comparing LIME results with ALMA observations.
In order to solve this issue we need to modify the ray tracer to allow finer
sampling in the central pixels and tried the two alternatives below without success.

There is an antialias option in the documentation, which performs stochastic
supersampling. We have tried with values of 10 and 100, but this did
not fully resolve the problem. It seems that this algorithm is currently not
implemented in the raytracing code.

We have also modified the value of maxNumRaysPerPixel in the raytrace function
from 20 to 0 to indicate that there is no maximum number of rays per pixel.
However this did not solve the problem either.

[mcordiner]

The problem arises when there is a strongly varying number density inside a given pixel (in this case, the central pixel(s). Using log(r) grid sampling results in a much higher density of grid points for the central, high-density regions. Given that the LIME raytracer sampling is chosen based on the cell density (one for each cell, up to the limit given by maxNumRaysPerPixel), the resulting average flux in the central pixel becomes biased by the high-density region towards r=0.

I wrote a fix for raytrace.c that implements a regular cartesian sampling of the central pixel(s) (similar to what RATRAN can do), which seems to solve the problem. I'd be happy to share this fix if anyone wants to try it.

[allegroLeiden]

Sorry not to reply to this until now. LIME support in the immediate future is problematic. I'm nearly at the end of my contract at Leiden and have had to devote my time to other things for some time. I am about to spend a few weeks on LIME, to try to fix some of the outstanding issues and to finalize the whole python/CASA interface additions. @mcordiner that sounds like a useful alternative but priority for me is to address #261 and then merge my python PR.

[allegroLeiden]

I've wondered for some time if the 'roughness' visible in the central pixels means that the sampling is too sparse there. Traditional LIME used a cockamamie algorithm bodged up to avoid the slowness of the rejection method of choosing random point locations to follow a prescribed density function. This could easily result in undersampling in the centre and is also poorly adapted to models which are not dominated by a central density maximum. par->samplingAlgorithm=1 uses my alternate algorithm, which is fast but which also follows the prescribed point density fairly exactly; however there is no guarantee either that this gives a good sampling density at all places. The only way to really decide on this is if one could look at the maximum optical depth (MOD) per link/edge between grid points. Ideally this ought to be of order unity. MOD<<1 means the grid points are denser than they need to be; MOD>>1 means they are too sparse. I have had plans to write a diagnostic for this but this has run against the same wall as with the rest of my wish list for LIME: lack of available man-hours.

[mcordiner]

In my opinion, the new raytracing algorithm could be improved by taking into account the physical (sky-projected) area associated with each ray, and then weighting the resulting pixel flux accordingly. To implement this (assuming a one-ray-per-grid-cell approach), you'd just need to know the sky-projected area of each cell.

Also, I guess it wouldn't be that hard to implement an option for RATRAN-style raytracing, that should give good results for the case of models dominated by a single maximum at the image centre, and allow easier benchmarking against that code.

[allegroLeiden]

I just had a deeper look at this and I confess I am a bit puzzled. What happens in raytrace is as follows:

• All internal grid points are projected to the image plane. (Some others are added on the locus of the circular projected model rim.) This forms the set of rays.
• Each ray is followed through the model and the flux at each channel calculated for it.
• Each ray's flux is added to the image pixel the ray passes through. (For every channel.) The tally of the rays passing through each pixel is kept.
• Pixels with more than 1 ray have their cumulative values divided by the number of rays passing through that pixel. I.e., this is just an average value. All other pixels have their values obtained by interpolation.

The average performed is I suppose a biased one, but if the model properties vary markedly within a pixel then any averaging scheme (e.g. sampling on a cartesian grid within the pixel, or ditto random- or quasi-random sampling) will be biased in that it will not reflect such variations well.

Mathematically one should not expect that the flux found by summing all pixels should be independent of the pixel sizes. A trivial inequality illustrating why not is

(a+b+c+d+e)/5 != [(a+b+c)/3 + (d+e)/2]/2

The code which implements the algorithm is not very safe, in that the tally of rays per pixel and the summing of flux are performed in different places. I'll make this more robust and see if that makes a difference.

The thing I notice, at least with the example model, is not so much over-estimation as the relatively non-uniform nature of the solution across the image. This leads to a 'speckled' appearance of the central pixels, i.e. those for which averaging is performed. I can't see any solution to this apart from increasing the number of grid points. I'll try this as well.

[allegroLeiden]

BTW @migueldvb you should have set maxNumRaysPerPixel to a value <0 to flag the desire for no upper bound.

[allegroLeiden]

I upload 4 images, which were made as follows:

• a1: standard model.
• a2: standard model, image at twice the resolution.
• b1: standard model, but with 8 times as many grid points (thus average edge lengths half that of standard model).
• b2: same as b1 but with twice the resolution.

All are displayed with the same brightness scale. All show the +500 m/s channel.

[allegroLeiden]

My take-home:

• The in-pixel sampling gives quite a reasonable representation of the actual solution, in that there are no marked or unexpected results from doubling the resolution.
• Halving the average distance between grid points does, as expected, remove a great deal of the spatial inhomogeneity of the solution.

So I can't really see an issue. If any of you have a model which displays one, I'd be happy to take a look at it.

[mcordiner]

I agree that the current raytracing implementation works well in many cases, as long as there aren't any sharp features (containing lots of grid points) with sizes much less than a single image pixel. Your above test images appear quite smoothly varying. However, in our problem case, there is a strong density increase towards the image center, that is well sampled by the logarithmic radial grid, but poorly sampled in terms of image pixels. This is actually for a model of a solar system body, so the size of the central density enhancement is ~ 1km, however, I can imagine protostellar envelopes with ~1/r*2 density profiles where a similar situation could occur, but of course, on a much larger size scale (and at at similarly greater distance). It's impractical for us to have images 10^5x10^5 pixels in size, so the central regions need a robust supersampling algorithm, or else their fluxes get severely overestimated.

The attached image shows a zoom of the central portion of the model, and the statistics are taken inside the magenta box. The box is the same physical size (in arcseconds) in both images, but the image on the right has 10 times smaller pixels sizes. The region on the left has a total flux (sum) 1.4 times the total flux of the image on the right, which is quite a severe discrepancy. I can email the model file if you're interested in having a look.

[allegroLeiden]

Ok I think I have a better appreciation of the problem now. You want to have more accurate flux estimates, which I guess is only reasonable. I can think of 3 candidate algorithms for improving matters, 2 suggested by @mcordiner and 1 which occurred to me:

1. Cartesian grid.
2. Weighting by cell area.
3. 2^D-tree sampling.

The 1st is probably the simplest to implement, but with this one would risk undersampling the flux, since it is easy to imagine a situation in which a sharp maximum falls between points. However the point about comparing to RATRAN is well made.

Weighting by cell area: these would have to be Voronoi cells. I guess we can get these without too much extra graft since qhull has to be called anyway for the sparse-pixel interpolation, although we get the Delaunay duals for this. Personally I am very leery of having anything more to do with qhull whatsoever, but since we already have exposure to this package, I suppose in for a penny, in for a pound. However, another objection to this scheme is that flux does not necessarily follow grid point density.

The 2^D-tree thing is already coded up in Lime (in module tree_random.c) as an alternative to the original manky hard-wired grid-and-smooth algorithm for selecting grid points. I tried to make it general, not hard-wired for D=3. It can be used not only to quickly generate random points following an arbitrary distribution but in principle also to estimate the integral of that distribution. This method would probably provide the most accurate fluxes per pixel of the three, although decent accuracy would probably require more samples (thus running time) than most people are willing to afford. In other words it may be too much gun for the bear.

Now the crunch is that I am very reluctant to do any more substantial work on Lime. I've just merged a giant PR with fixes to the CASA interface, and since that was what they mostly hired me to do, and since I have had very little time to pursue my own research interests during the last 3.5 years, I am looking forward now to at last being able to give them some attention during the small amount of time I have left in this contract. So if one of you wants to submit a PR for a changed algorithm in the way you desire, I'll have a look at it, and merge it if the code meets reasonable standards. That's about the fairest I can say.

[christianbrinch]

A historical note: LIME used to have a RATRAN style raytracer in the early versions. In fact, the raytracer was a simple translation of Michiels SKY code. It worked mostly fine, but we encountered the exact same problem as the one described here, namely, that the center pixel(s) were off in some cases. SKY improved on this by introducing oversampling, which helped but did not solve the issue. Kees Dullemond had experienced the same issue with the raytracer in RADMC-3D and he tried to solve it with a clever subsampling scheme, which unfortunately could become extremely slow.
Kees and I devised a solution which uses gridpoint based raytracing rather than the old raster based raytracing. A variation of this method (which is described in this paper: http://adsabs.harvard.edu/abs/2014MNRAS.440.3285B, see fig 6 for an illustration of the problem: there is a hole in the center of the raster image of the disk because the four central rays miss the peak.) was implemented in LIME some 5-6 years ago by me and iirc later on also in RADMC-3D by Kees. It works quite well, so I don't know why Michiel decided to go back to the old raster method. It is slower, it does not conserve the flux properly and it has aliasing issues.

It is true that grid point based raytracing rely critically on proper sampling of the model. However, it will give you an exact mapping of the flux associated with the sampling of the model - pixel based raytracing will not. The idea is that if you sample the model properly, the flux will only change linearly between gridpoints and therefore, by triangulating the rays, it is possible to interpolate the flux between the rays (not generally possible for pixels) and achieve arbitrary image resolution, even with relatively few rays.

[mcordiner]

To me, a decent approach would be to continue with the current algorithm, but do a weighted average over the rays per pixel based on the projected area of each grid cell (instead of the current, unweighted average). However, I'm not at all familiar with what Qhull can do. If it can return the sky-projected cell areas that would be ideal, or alternatively, we can calculate them using the actual cell areas and their normal vectors, if those are available?

Christian - which version of LIME implements the Brinch and Dullemond method? If it works well, surely it would be a good idea to make it available again (perhaps as an option in the latest version).

[christianbrinch]

Certainly 1.4 and 1.5. Not sure about 1.6. Probably not.

[allegroLeiden]

Direct generation of visibilities a la Brinch and Dullemond is a good idea, but AFAIR this was never fully implemented in LIME.

[allegroLeiden]

Just looked at the code. In the version of 1.4 which I got from Jes Jorgensen, a sky-plane grid image is produced. There's no direct generation of visibilities. The algorithm was essentially as follows:

• Rays are passed through each grid point and projected onto the image plane.
• Qhull is called a second time, in 2D, for these projected points.
• A sky-plane i.e. orthogonal gridded image is constructed via what looks like an interpolation scheme.

This algorithm was taken offline in version 1.5, although the code was kept, since it seemed a useful approach which however needed some work. A version of the algorithm, with simpler and more transparent interpolation, and with separate treatment of low-density vs. high-density pixels, was reinstated in 1.6. This has been retained with only minor changes since that time.