"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Libtmp/Transform/transform.pd" between
PDL-2.074.tar.gz and PDL-2.075.tar.gz

About: PDL (Perl Data Language) aims to turn perl into an efficient numerical language for scientific computing (similar to IDL and MatLab).

transform.pd  (PDL-2.074):transform.pd  (PDL-2.075)
skipping to change at line 38 skipping to change at line 38
PDL::Transform is a convenient way to represent coordinate PDL::Transform is a convenient way to represent coordinate
transformations and resample images. It embodies functions mapping transformations and resample images. It embodies functions mapping
R^N -> R^M, both with and without inverses. Provision exists for R^N -> R^M, both with and without inverses. Provision exists for
parametrizing functions, and for composing them. You can use this parametrizing functions, and for composing them. You can use this
part of the Transform object to keep track of arbitrary functions part of the Transform object to keep track of arbitrary functions
mapping R^N -> R^M with or without inverses. mapping R^N -> R^M with or without inverses.
The simplest way to use a Transform object is to transform vector The simplest way to use a Transform object is to transform vector
data between coordinate systems. The L</apply> method data between coordinate systems. The L</apply> method
accepts a PDL whose 0th dimension is coordinate index (all other accepts a PDL whose 0th dimension is coordinate index (all other
dimensions are threaded over) and transforms the vectors into the new dimensions are broadcasted over) and transforms the vectors into the new
coordinate system. coordinate system.
Transform also includes image resampling, via the L</map> method. Transform also includes image resampling, via the L</map> method.
You define a coordinate transform using a Transform object, then use You define a coordinate transform using a Transform object, then use
it to remap an image PDL. The output is a remapped, resampled image. it to remap an image PDL. The output is a remapped, resampled image.
You can define and compose several transformations, then apply them You can define and compose several transformations, then apply them
all at once to an image. The image is interpolated only once, when all at once to an image. The image is interpolated only once, when
all the composed transformations are applied. all the composed transformations are applied.
skipping to change at line 145 skipping to change at line 145
called with the input coordinate, and the "params" hash. This called with the input coordinate, and the "params" hash. This
springboarding is done via explicit ref rather than by subclassing, springboarding is done via explicit ref rather than by subclassing,
for convenience both in coding new transforms (just add the for convenience both in coding new transforms (just add the
appropriate sub to the module) and in adding custom transforms at appropriate sub to the module) and in adding custom transforms at
run-time. Note that, if possible, new C<func>s should support run-time. Note that, if possible, new C<func>s should support
L<inplace|PDL::Core/inplace> operation to save memory when the data are flagged L<inplace|PDL::Core/inplace> operation to save memory when the data are flagged
inplace. But C<func> should always return its result even when inplace. But C<func> should always return its result even when
flagged to compute in-place. flagged to compute in-place.
C<func> should treat the 0th dimension of its input as a dimensional C<func> should treat the 0th dimension of its input as a dimensional
index (running 0..N-1 for R^N operation) and thread over all other input index (running 0..N-1 for R^N operation) and broadcast over all other input
dimensions. dimensions.
=item inv =item inv
Ref to an inverse method that reverses the transformation. It must Ref to an inverse method that reverses the transformation. It must
accept the same "params" hash that the forward method accepts. This accept the same "params" hash that the forward method accepts. This
key can be left undefined in cases where there is no inverse. key can be left undefined in cases where there is no inverse.
=item idim, odim =item idim, odim
skipping to change at line 621 skipping to change at line 621
OtherPars=>'pdl *in; pdl *out; pdl *map; SV *boundary; SV *method; OtherPars=>'pdl *in; pdl *out; pdl *map; SV *boundary; SV *method;
long big; double blur; double sv_min; char flux; SV *bv', long big; double blur; double sv_min; char flux; SV *bv',
Code => <<'+==EOD_map_c_code==', Code => <<'+==EOD_map_c_code==',
/* /*
* Pixel interpolation & averaging code * Pixel interpolation & averaging code
* *
* Calls a common coordinate-transformation block (see following hdr) * Calls a common coordinate-transformation block (see following hdr)
* that isn't dependent on the type of the input variable. * that isn't dependent on the type of the input variable.
* *
* The inputs are SVs to avoid hassling with threadloops; threading * The inputs are SVs to avoid hassling with broadcastloops; broadcasting
* is handled internally. To simplify the threading business, any * is handled internally. To simplify the broadcasting business, any
* thread dimensions should all be collapsed to a single one by the * broadcast dimensions should all be collapsed to a single one by the
* perl front-end. * perl front-end.
* *
*/ */
pdl *in = $COMP(in); pdl *in = $COMP(in);
pdl *out = $COMP(out); pdl *out = $COMP(out);
pdl *map = $COMP(map); pdl *map = $COMP(map);
PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */ PDL_Indx ndims = map->ndims -1; /* Number of dimensions we're working in */
PDL_Double tmp[3*ndims*ndims+ndims]; /* Workspace for prefrobnication */ PDL_Double tmp[3*ndims*ndims+ndims]; /* Workspace for prefrobnication */
PDL_Indx ovec[ndims]; /* output pixel loop vector */ PDL_Indx ovec[ndims]; /* output pixel loop vector */
PDL_Indx ivec[ndims]; /* input pixel loop vector */ PDL_Indx ivec[ndims]; /* input pixel loop vector */
PDL_Indx ibvec[ndims]; /* input pixel base offset vector */ PDL_Indx ibvec[ndims]; /* input pixel base offset vector */
PDL_Double dvec[ndims]; /* Residual vector for linearization */ PDL_Double dvec[ndims]; /* Residual vector for linearization */
PDL_Double tvec[ndims]; /* Temporary floating-point vector */ PDL_Double tvec[ndims]; /* Temporary floating-point vector */
PDL_Double acc[in->dims[ndims]]; /* Threaded accumulator */ PDL_Double acc[in->dims[ndims]]; /* Broadcasted accumulator */
PDL_Double wgt[in->dims[ndims]]; /* Threaded weight accumulator */ PDL_Double wgt[in->dims[ndims]]; /* Broadcasted weight accumulator */
PDL_Double wgt2[in->dims[ndims]]; /* Threaded weight accumulator for badval fin PDL_Double wgt2[in->dims[ndims]]; /* Broadcasted weight accumulator for badval
ding */ finding */
char bounds[ndims]; /* Boundary condition packed string */ char bounds[ndims]; /* Boundary condition packed string */
PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */ PDL_Indx index_stash[ndims]; /* Stash to store the opening index of dim sample scans */
char method; /* Method identifier (gets one of 'h','g') */ char method; /* Method identifier (gets one of 'h','g') */
PDL_Long big; /* Max size of input footprint for each pix */ PDL_Long big; /* Max size of input footprint for each pix */
PDL_Double blur; /* Scaling of filter */ PDL_Double blur; /* Scaling of filter */
PDL_Double sv_min; /* minimum singular value */ PDL_Double sv_min; /* minimum singular value */
char flux; /* Flag to indicate flux conservation */ char flux; /* Flag to indicate flux conservation */
PDL_Double *map_ptr; PDL_Double *map_ptr;
PDL_Long i, j; PDL_Long i, j;
$GENERIC() badval = SvNV($COMP(bv)); $GENERIC() badval = SvNV($COMP(bv));
skipping to change at line 1076 skipping to change at line 1076
default: default:
{ {
char buf[80]; char buf[80];
sprintf(buf,"This can't happen: method='%c'",method); sprintf(buf,"This can't happen: method='%c'",method);
$CROAK("%s", buf); $CROAK("%s", buf);
} }
} }
{ /* convenience block -- accumulate the current point into the weight ed sum. */ { /* convenience block -- accumulate the current point into the weight ed sum. */
/* This is more than simple assignment because we have our own expli cit poor */ /* This is more than simple assignment because we have our own expli cit poor */
/* man's threadloop here, so we accumulate each threaded element sep arately. */ /* man's broadcastloop here, so we accumulate each broadcasted eleme nt separately. */
$GENERIC() *dat = (($GENERIC() *)(in->data)) + i_off; $GENERIC() *dat = (($GENERIC() *)(in->data)) + i_off;
PDL_Indx max = out->dims[ndims]; PDL_Indx max = out->dims[ndims];
for( i=0; i < max; i++ ) { for( i=0; i < max; i++ ) {
if( (badval==0) || (*dat != badval) ) { if( (badval==0) || (*dat != badval) ) {
acc[i] += *dat * alpha; acc[i] += *dat * alpha;
dat += in->dimincs[ndims]; dat += in->dimincs[ndims];
wgt[i] += alpha; wgt[i] += alpha;
} }
wgt2[i] += alpha; // Accumulate what weight we would have with no bad values wgt2[i] += alpha; // Accumulate what weight we would have with no bad values
} }
skipping to change at line 1496 skipping to change at line 1496
This value lets you set the lower limit of the transformation's This value lets you set the lower limit of the transformation's
singular values in the hanning and gaussian methods, limiting the singular values in the hanning and gaussian methods, limiting the
minimum radius of influence associated with each output pixel. Large minimum radius of influence associated with each output pixel. Large
numbers yield smoother interpolation in magnified parts of the image numbers yield smoother interpolation in magnified parts of the image
but don't affect reduced parts of the image. but don't affect reduced parts of the image.
=item big, Big (default = 0.5) =item big, Big (default = 0.5)
This is the largest allowable input spot size which may be mapped to a This is the largest allowable input spot size which may be mapped to a
single output pixel by the hanning and gaussian methods, in units of single output pixel by the hanning and gaussian methods, in units of
the largest non-thread input dimension. (i.e. the default won't let the largest non-broadcast input dimension. (i.e. the default won't let
you reduce the original image to less than 5 pixels across). This places you reduce the original image to less than 5 pixels across). This places
a limit on how long the processing can take for pathological transformations. a limit on how long the processing can take for pathological transformations.
Smaller numbers keep the code from hanging for a long time; larger numbers Smaller numbers keep the code from hanging for a long time; larger numbers
provide for photometric accuracy in more pathological cases. Numbers larer provide for photometric accuracy in more pathological cases. Numbers larer
than 1.0 are silly, because they allow the entire input array to be compressed than 1.0 are silly, because they allow the entire input array to be compressed
into a region smaller than a single pixel. into a region smaller than a single pixel.
Wherever an output pixel would require averaging over an area that is too Wherever an output pixel would require averaging over an area that is too
big in input space, it instead gets NaN or the bad value. big in input space, it instead gets NaN or the bad value.
skipping to change at line 1911 skipping to change at line 1911
############################## ##############################
## Anti-aliasing code: ## Anti-aliasing code:
## Condition the input and call the pixelwise C interpolator. ## Condition the input and call the pixelwise C interpolator.
## ##
barf("PDL::Transform::map: Too many dims in transformation\n") barf("PDL::Transform::map: Too many dims in transformation\n")
if($in->ndims < $idx->ndims-1); if($in->ndims < $idx->ndims-1);
#################### ####################
## Condition the threading -- pixelwise interpolator only threads ## Condition the broadcasting -- pixelwise interpolator only broadcasts
## in 1 dimension, so squish all thread dimensions into 1, if necessary ## in 1 dimension, so squish all broadcast dimensions into 1, if necessary
my @iddims = $idx->dims; my @iddims = $idx->dims;
my $in2 = $in->ndims == $#iddims my $in2 = $in->ndims == $#iddims
? $in->dummy(-1,1) ? $in->dummy(-1,1)
: $in->reorder($nd..$in->ndims-1, 0..$nd-1) : $in->reorder($nd..$in->ndims-1, 0..$nd-1)
->clump($in->ndims - $nd) ->clump($in->ndims - $nd)
->mv(0,-1); ->mv(0,-1);
#################### ####################
# Allocate the output array # Allocate the output array
my $o2 = PDL->new_from_specification($in2->type, my $o2 = PDL->new_from_specification($in2->type,
skipping to change at line 2338 skipping to change at line 2338
primitive; if you want more, try composing the linear transform with primitive; if you want more, try composing the linear transform with
this one. this one.
The prescribed values in the lookup table are treated as The prescribed values in the lookup table are treated as
pixel-centered: that is, if your input array has N elements per row pixel-centered: that is, if your input array has N elements per row
then valid data exist between the locations (-0.5) and (N-0.5) in then valid data exist between the locations (-0.5) and (N-0.5) in
lookup pixel space, because the pixels (which are numbered from 0 to lookup pixel space, because the pixels (which are numbered from 0 to
N-1) are centered on their locations. N-1) are centered on their locations.
Lookup is done using L<interpND|PDL::Primitive/interpnd>, so the boundary condit ions Lookup is done using L<interpND|PDL::Primitive/interpnd>, so the boundary condit ions
and threading behaviour follow from that. and broadcasting behaviour follow from that.
The indexed-over dimensions come first in the table, followed by a The indexed-over dimensions come first in the table, followed by a
single dimension containing the column vector to be output for each single dimension containing the column vector to be output for each
set of other dimensions -- ie to output 2-vectors from 2 input set of other dimensions -- ie to output 2-vectors from 2 input
parameters, each of which can range from 0 to 49, you want an index parameters, each of which can range from 0 to 49, you want an index
that has dimension list (50,50,2). For the identity lookup table that has dimension list (50,50,2). For the identity lookup table
you could use C<cat(xvals(50,50),yvals(50,50))>. you could use C<cat(xvals(50,50),yvals(50,50))>.
If you want to output a single value per input vector, you still need If you want to output a single value per input vector, you still need
that last index threading dimension -- if necessary, use C<dummy(-1,1)>. that last index broadcasting dimension -- if necessary, use C<dummy(-1,1)>.
The lookup index scaling is: out = lookup[ (scale * data) + offset ]. The lookup index scaling is: out = lookup[ (scale * data) + offset ].
A simplistic table inversion routine is included. This means that A simplistic table inversion routine is included. This means that
you can (for example) use the C<map> method with C<t_lookup> transformations. you can (for example) use the C<map> method with C<t_lookup> transformations.
But the table inversion is exceedingly slow, and not practical for tables But the table inversion is exceedingly slow, and not practical for tables
larger than about 100x100. The inversion table is calculated in its larger than about 100x100. The inversion table is calculated in its
entirety the first time it is needed, and then cached until the object is entirety the first time it is needed, and then cached until the object is
destroyed. destroyed.
skipping to change at line 2464 skipping to change at line 2464
if($data->dim(0) > $me->{idim}) { if($data->dim(0) > $me->{idim}) {
barf("Too many dims (".$data->dim(0).") for your table (".$me->{idim}.")\n "); barf("Too many dims (".$data->dim(0).") for your table (".$me->{idim}.")\n ");
}; };
my($x)= ($table my($x)= ($table
->interpND(float($data) * $scale + $offset, ->interpND(float($data) * $scale + $offset,
$p->{interpND_opt} $p->{interpND_opt}
) )
); );
# Put the index dimension (and threaded indices) back at the front of # Put the index dimension (and broadcasted indices) back at the front of
# the dimension list. # the dimension list.
my($dnd) = $data->ndims - 1; my($dnd) = $data->ndims - 1;
return ($x -> ndims > $data->ndims - 1) ? return ($x -> ndims > $data->ndims - 1) ?
($x->reorder( $dnd..($dnd + $table->ndims - $data->dim(0)-1) ($x->reorder( $dnd..($dnd + $table->ndims - $data->dim(0)-1)
, 0..$data->ndims-2 , 0..$data->ndims-2
) )
) : $x; ) : $x;
}; };
$me->{func} = sub {my($data,$p) = @_; &$lookup_func($data,$p,$p->{table},$p-> {scale},$p->{offset})}; $me->{func} = sub {my($data,$p) = @_; &$lookup_func($data,$p,$p->{table},$p-> {scale},$p->{offset})};
 End of changes. 10 change blocks. 
16 lines changed or deleted 16 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)