[GRASS-dev] Re: [GRASS-user] low radiance values after i.atcor and i.topo.corr

[Glynn Clements]

Markus Neteler wrote:

> >>  r.mapcalc b10=lsat7_2000_10/10000.
> >
> > [snip]
> >
> >> Error: unreliable clustering
> >> try a smaller initial number of clusters
> >> Warning: Removed a singular subsignature; number 1; -1 remain
> ...
> > Could it be that the definition of ZERO in subcluster.c:
> >
> >        #define ZERO 1e-10
> >
> > is too high when the inputs are scaled down?
> 
> Apparently yes.
> 
> > For an NxN matrix, scaling all of the the elements down by a factor of
> > K will scale the determinant by K^N.
> >
> > AFAICT[1], the matrix elements are roughly proportional to the squares
> > of the input values, so if you're scaling the values down by 1e-5, the
> > matrix elements are scaled down by 1e-10, and with 7 bands the
> > determinant would be scaled down by 1e-70.
> >
> > [1] Which isn't saying much; it's not like I actually understand the
> > algorithm being used here.
> 
> Documentation and a potentially more advanced version I found here:
> http://cobweb.ecn.purdue.edu/~bouman/software/segmentation/
> 
> > Can you try dividing by 10, 100, etc rather than 10000, to see if
> > there's a threshold at which the accuracy starts to be an issue?
> >
> > And/or make ZERO much smaller.
> 
> For /10000. it works with ZERO = 0.0 only. The result is reasonable.
> 
> Keeping ZERO = 0.0 and dividing by /10. only produces the same result
> (but the Rissanen output differs), same for /100. Differences with r.mapcalc
> are NULL. Same for /100000.

In retrospact, that's to be expected. Even dividing by 10 would (IIUC)
result in the determinant being scaled by 1e-14.

> Looks like solved with ZERO patch to be submitted?

Yes, but it may not be trivial. Changing ZERO to 0.0 would essentially
disable the check for a singular matrix (I don't know whether that's
acceptable in practice).

Ideally, ZERO should be a variable whose magnitude varies with the
data. One option would be to compute e.g. the RMS of the entire data
as it is read in. Another option would be to compute the RMS of the
matrix elements within compute_constants().

That still leaves the issue of determining a constant factor. Do you
have any integer maps which produce singular matrices for lower
numbers of subclasses, so that a reasonable value can be determined?

Actually, I've just noticed a couple of issues with G_ludcmp() (in
lib/gmath/lu.c), which is used to perform the inversion.

First, it adjusts its calculations according the range of the
elements:

    for (i=0;i<n;i++) 
    {
        big=0.0;
        for (j=0;j<n;j++)
            if ((temp=fabs(a[i][j])) > big)
		big=temp;
        if (big == 0.0)
	{
	    *d = 0.0;
	    return 0; /* Singular matrix  */
	}
        vv[i]=1.0/big;
    }

Second has its own constant:

	#define TINY 1.0e-20;

This may make the result useless if the matrix elements are themselves
relatively small.

So, it may be necessary to pre-scale the matrix elements by some
factor, then scale the inverse (and determinant) to compensate.

-- 
Glynn Clements <glynn at gclements.plus.com>