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 |