Re: [unclassified] contourc

[Alexander Barth]

On 10/7/07, Kai Habel <address@hidden> wrote:

Alexander Barth schrieb:
> Hi all,
> I seems that there are several implementations for computing contours:
> the __contourc__.cc originally from plplot and the script
> implementation in contourf.m from octplot. For me, the implementation
> in contourf.m works more reliably than the __contourc__.cc
> implementation. __contourc__.cc produces sometimes discontinuous
> contour lines. As contourf.m is written, it would be very easy to
> extract the contour generation code and make a new contourc.m to
> replace __contourc__.cc.
>
> I would also be interested to implement arbitrary quadrilateral meshes.
> What do you think?
>
> Thanks,
> Alex
>
Hello Alexander, Hello all,

as far as i know there are two implementations for computing contours.
The 'old' __contourc__.cc (before April 07) from octplot (originating
from plplot etc.) and a improved 'new' version (April' 07) also from
octplot to support filled contours. This was required (as you have
noticed), because the 'old' version generates in some cases broken contours.
The 'old' version has also been included in the octave tree.

If you talk about the octplot/contourf.m implemetation, this does not
include any contour finding as it relies on contourc.m/__contourc__.cc
for this purpose. The combination of contourf.m and __contourc__.cc from
octplot should work pretty well.

Unfortunately the 'new' version of __contourc__.cc (also to be found in
the octplot svn) has some problems as well. So, I have a third version
of __contourc__.cc - a major rewrite, which is almost finished. I plan
to post it for discussion on the maintainers list very soon. (If all
works out well, maybe this week).

Regarding your last question, is this related to this thread [1] or a
more general approach?

Kai

[1]
http://www.nabble.com/contour-plots-only-work-for-meshgrid-%28uniform%29-spaces-and-not-other-spaces.-tf4541422.html

The approach that I had in mind for general quadrilateral data is indeed the same as described in the thread you sent me. The only difference is that I would compute the contours in index space (i,j) and then map the contours in (x,y) space. Basically something like:

    [ii, jj] = meshgrid(1 : columns(Z),1 : rows(Z));
    # call to contourc for "meshgridded" data
    c_ij  = contourc(ii, jj, Z, lvl);

    c = zeros(size(c_ij));
    loop over segments    
       k = start_index : end_index
       c(1,k) = interp2(ii,jj,X,c_ij(1,k),c_ij(2,k));
       c(2,k) = interp2(ii,jj,Y,c_ij(1,k),c_ij(2,k));

This is very similar to the approach of Peter. The only case where the two approaches differ is when the coordinate system is singular (if I understand the patch of Peter correctly). For example:

[r,theta] = ndgrid(0:10, 0:2*pi);
x = r .* cos(theta);
y = r .* sin(theta);

in which case all y(:,1) will be zero.
But I haven't test the patch of Peter and I might be wrong.

Thanks
Alex