[GRASS-dev] Re: [SoC] gdalwarp OpenCL Performance (Week 9)

[Hamish]

[cc'ing the world, hope that's ok]

Seth wrote:
> The main problem isn't how
> computationally intensive com_par() is, it's that if I want
> to rewrite the calculate() main loop into OpenCL, then I
> need to rewrite everything that is called by it. When I run
> the calculate() OpenCL kernel on the GPU, it can't go back
> and call the pj_do_proj() on the host CPU.
> 
> So we either need to rewrite com_par() to pull pj_do_proj()
> outside of the main calculate() loop, or we rewrite all of
> pj_do_proj() into OpenCL. The former is difficult, but I
> believe the latter is effectively impossible.

I concur.

> Thus we need to figure a way for getting pj_do_proj() outside
> of the main calculate() loop and passing all it's results into
> com_par() in an array. Then we can copy the array to the GPU,
> and run the kernel without interruption. Does this make sense?

I've got another idea: all it's doing is setting up a triangle
in Cartesian space to make the trigonometry work easier. I don't
think it is really necessary to do a full transform back to the
map projection to accomplish that, if we neglect the ellipsoid
we can just set east and north axes to nearly the same scale by *cos(lat).

While that might not be as nice as a real map projection, it
actually avoids the 'convergence angle'* problem, which may be
mucking up the right-triangle assumption anyway.

[*] convergence angle is the difference between grid north and
true north at a point, for spearfish this is about 0.9 degrees.
  (in grass: `g.region -n`)


try the attached patch against trunk which removes the 
pj_do_proj() call and replaces it with a cheap estimation. The
resulting stepsinangle and stepcosangle vars seem within about
0.2-0.4 of the pj_do_proj() version, but after factoring in no
convergence angle error I couldn't tell you how much better or
worse the result is, or vs. the old sunAzimuthAngle method, or
how much this even matters. In my tests I got perfectly identical
output maps before and after the change. Perhaps I was not using
the right combo of options to trigger that, but I was using
shading and no horizon maps, so I'm pretty sure I use that code
fully. ?

bonus: the patch makes r.sun in Mode 2 run 16% faster.

(fwiw I don't really follow why the grid resolution (stepxy) is
in that calc.. shrug)


> Honestly, if com_par() was particularly computationally
> intensive, that would only be a reason to put as much of it
> as possible on the GPU. The CPU is slow. :)

fwiw my tests were in the spearfish sample dataset:

# DEM (mode 1)
g.region rast=elevation.dem
r.sun -s elevation.dem glob_rad=global_radDay180.14 day=180 time=14

r.colors -e global_radDay180.14 col=grey
d.shadedmap rel=elevation.dem drape=global_radDay180.14 br=30


# Gaussian mound (mode 2)
# r.surf.volcano is from the grass addons repo.
g.region -d
r.surf.volcano out=gauss method=gaussian kurtosis=1
r.sun -s elevin="gauss" glob_rad="rad.global.30minT" day=180 step=0.5

# check for deviations between old and new code:
r.mapcalc "diff = rad.global.ORIG - rad.global.NEW"
r.univar diff


> Besides r.sun, my next project will be getting my warping
> code integrated into GRASS's warping & interpolation
> modules. Depending on how my PhD goes, I may be putting
> effort into other GRASS modules this fall, though. ;) What's
> related to land cover classification if I'm writing a new
> classifier module already?

umm, pass. beyond my area of expertise. [cc list]
 
> Once I get some appropriate 4+ band imagery, I'll be making
> some tweaks for the CPU code to use SSE also. Debugging
> OpenCL is hell, but it's pretty powerful & versatile
> once you get it running. (Looks like I have some 6-band
> imagery, I know what I'm doing today...)

>From the GRASS download page you can find the Spearfish Imagery
extra mapset, and also the new North Carolina sample dataset.
Both include multiband LANDSAT data. Also the NC dataset has
geotiff versions available for download if you want to push it
straight into gdalwarp. otherwise just use `r.out.gdal` to
create the geotiffs from the grass data.

Also you can go to
 http://grass.osgeo.org/wiki/Global_datasets#LANDSAT
where there is a link to where you can access the USGS LANDSAT
archive. (requests may take a few days to process)


regards,
Hamish



      
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rsun_no_pj.diff
Type: text/x-patch
Size: 2457 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-dev/attachments/20100729/a81f20da/rsun_no_pj.bin