Re: [unclassified] contourc

[Peter A. Gustafson]

On Sunday 07 October 2007 22:49:05 Alexander Barth wrote:
> 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

I tried running the above code with my mapping and it failed with divide by 
zero.  I haven't tried Alex's code, but I admit I did not consider a singular 
case when I provided my mapping.

Pete