On Feb 26, 2007, at 12:10 PM, Jarle Brinchmann wrote:
I'm afraid that's not exactly what I'm looking for. After
re-analyzing the problem and reading chapter 3 of the PDL book, I
realize what I actually need to do is find or write a "primitive"
code
that will project a line into a plane, and then let threading take
care of doing that to each line in my image and give me a cube.
Again, sorry to be so late getting into the fray here. It sounds
like the problem is to take a sparse collection of (N-1)-dimensional
data and interpolate them on to a full N-dimensional data set. For
N=2, the data are linear and the data set is an image.
Here's example code to do it in that case -- the more general case is
almost as trivial but less clear to read, and is left as an exercise :-)
On input, let the data sets be contained in a list of 1-D PDLS,
called @dat. Also, let the transformations to
2-D be defined by a list of 2x2-PDLs (matrices) called @mat and a
list of 2-PDL offsets called @of, and let $nx and $ny contain the
size of the output array onto which you want to interpolate. Then do
this:
# Assemble 2-D coordinates for each data point in each line
undef @coords;
for $i(0..$#dat) {
$coords = ($mat[$i] x ndcoords($data2)->(*1))->((0)) + $of[i];
push(@coords,$coords);
}
# Now make $coords: a 2xM-PDL list of all data point coordinates
$coords = $coords[0]->glue(1,$coords[1..$#coords]);
# Convert all the data points to a single M-PDL
$data = $dat[0]->glue(0,$dat[1..$#dat]);
# ocoords gets a 2x($nx)x($ny) array with the coordinates of each
pixel in the
# output array.
$ocoords = ndcoords($nx,$ny)->clump(1..2);
# odiff gets a 2xMx(nx)x(ny) array with the differences between the
coordinates
# of each pixel in the list and each pixel in the final image.
$odiff = $coords->(:,:,*1) - $ocoords->(:,*1);
# Take 1/R^2: square the differences, collapse by summation, and
take reciprocal
$fact = 1.0/(($odiff*$odiff)->sumover)
# Now multiply by the data and collapse to (nx)x(ny) by summation,
and divide
# pixelwise by the summed factors themselves to get the weighted
average values.
$out = ($fact*$data)->sumover / $fact->sumover;
You can of course use a different weighting function by replacing the
1/R^2 with your favorite analytic expression.
The above works great if you can fit several times (nx)x(ny)xM bytes
in RAM, but as with all threading tricks it trades memory and CPU
cycles for programmer time. For a production run or for large data
sets you're better off with a PP solution.
Cheers,
Craig
...
I
think I can manage that, but I haven't worked with c much, and I'm
not
sure I understand how to use PP completely, but hopefully I'll learn.
Sort of on that note, I was hoping to use the code for, say,
diagonal,
or some other similar code, as an example, but I can't find it.
Core.pm and Slices.pm have references to diagonal and diagonalI, but
I'm not sure where to follow those. Where can I find out about how
the code is structured/where the source is?
Not going to be much help about your actual coding problem - it
does sound indeed that you want to use PP for it. As for finding
the source, this can sometimes be a bit tricky because the manual
pages sometimes lead you astray..
diagonal is a good case in point. When you do help diagonal in the
perldl shell you will be told that it is defined in Core.pm. This
is in Basic/Core/Core.pm in the PDL source tree. A lot of key
functionality is defined there. However the diagonal (and some
other routines) are defined as simple wrappers around a different
function, namely diagonalI. Here then comes the tricky part.. where
is this.
Searching blindly will get you there (something like find . -name
'*.pd' -exec grep -i -H 'diagonalI' {}\;) but who wants to do that
each time. Anyway, first priority is to find a .pd file because
these are the ones that contain the PP definitions you want. And I
find that unless I want to do something hairy I tend to only need
to look at two files for core functionality: Basic/Primitive/
primitive.pd and Basic/Slices/slices.pd - these define many of the
most basic functions and are useful to browse to get some information.
That said, diagonalI is _not_ the function I'd start with if I were
you... In fact, it does sound like given your problem that what you
really want to do is to look at the routines in Image2D/image2d.pd
and ImageND/imagend.pd and get a feel for what you need from these.
Cheers,
Jarle.
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
_______________________________________________
Perldl mailing list
[email protected]
http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
