~ubuntu-branches/ubuntu/maverick/pdl/maverick

 

       {  /* convenience block for bounds handling and data weighting */

                  if(sum > 4) /* 2 pixels -- four half-widths */

                    i = ndims; /* exit early if this pixel is too far outside the footprint of the ideal point */

       if(bound_ok)

         }/* end of bound_ok check */

 

       }  /* end of data collection convenience block */

     } while(i<ndims);  /* end of total data accumulation loop */

        my $tosize;

        ($tosize = $osize->where($osize == 0)) .= 1.0;

          $out->hdr->{"CRPIX$d"} = 1 + ($out->dim($d-1)-1)/2 ;

  my $method = _opt($opt,['m','method','Method'],undef);

  my $idx = $me->invert($ndc->double);

  my $integrate = ($method =~ m/^[jghJGH]/) if defined($method);

A simplistic table inversion routine is included.  This means that

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

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

destroyed.

  $a = float rfits('m51.fits');

  $lookup = xvals(128,128) -> cat( 10**(yvals(128,128)/50) * 256/10**2.55 );

  $me->{odim} = $source->dim($source->ndims-1);

     my($data,$p,$table,$scale,$offset) = @_;

 

    $data = pdl($data) 

      unless ((ref $data) && (UNIVERSAL::isa($data,'PDL')));

      $main::foo = $data;

    my($a)= ($table

             ->interpND(float($data) * $scale + $offset,

      ($a->reorder( $dnd..($dnd + $table->ndims - $data->dim(0)-1)

  $me->{func} = sub {my($data,$p) = @_;  &$lookup_func($data,$p,$p->{table},$p->{scale},$p->{offset})};

      my $data = shift;

        if($me->{idim} != $me->{odim}) {

         barf("t_lookup: can't calculate an inverse of a projection operation! (idim != odim)");

        }

       print "t_lookup: Warning, table inversion is only weakly supported.  \n   I'll try to invert it using a pretty boneheaded algorithm...\n  (If it takes too long, consider writing a faster algorithm!)\n   Calculating inverse table (will be cached)...\n" if($PDL::verbose || $PDL::debug || $PDL::Transform::debug);

        my $itable = zeroes($p->{table});

        my $minvals = $p->{table}->clump($me->{idim})->minimum;

        my $maxvals = $p->{table}->clump($me->{idim})->maximum;

        

        # Scale so that the range runs from 0 through the top pixel in the table

        my $scale = (  pdl(  $itable->dims  )->(0:-2)-1  ) / 

                    (($maxvals - $minvals)+ (($maxvals-$minvals) == 0));

        my $offset = - ($minvals/$scale);

 

        $p->{iscale} = $scale;

        $p->{ioffset} = $offset;

                      

        my $enl_scale = $p->{'enl_scale'} || 10;

        my $smallcoords = ndcoords((pdl($enl_scale * 2 + 1)->at(0)) x $me->{idim})/ $enl_scale - 1;

 

        # $oloop runs over (point, index) for all points in the output table, in 

        # $p->{table} output space

        $oloop = ndcoords($itable->mv(-1,0)->((0)))->

                        mv(0,-1)->

                        clump($itable->ndims-1);  # oloop: (pixel, index)

        $oloop->mv(-1,0) -= $offset;

        $oloop->mv(-1,0) /= $scale;

 

        # The Right Thing to do here is to take the outer product of $itable and $otable, then collapse 

        # to find minimum distance.  But memory demands for that would be HUGE.  Instead we do an 

        # elementwise calculation.

 

        print "t_lookup: inverting ".$oloop->dim(0)." points...\n" if($PDL::verbose || $PDL::debug || $PDL::Transform::debug);

        my $pt = $p->{table}->mv(-1,0); # pt runs (index, x,y,...)

 

        my $itable_flattened = zeroes($oloop);

 

        for $i(0..$oloop->dim(0)-1) {

 

            my $olp = $oloop->(($i));                # olp runs (index)

            my $diff = ($pt - $olp);                 # diff runs (index, x, y, ...)

            my $r2 = ($diff * $diff)->sumover;       # r2 runs (x,y,...)

            my $c = whichND($r2==$r2->min)->(:,(0)); # c runs (index)

 

            # Now zero in on the neighborhood around the point of closest approach.

            my $neighborhood = $p->{table}->interpND($c + $smallcoords,{method=>'linear',bound=>'t'})->

                     mv(-1,0); # neighborhood runs (index, dx, dy,...)

            $diff = $neighborhood - $olp;        # diff runs (index, dx, dy, ...)

            $r2 = ($diff * $diff)->sumover;      # r2 runs (dx,dy,...)

            my $dc = $smallcoords->mv(0,-1)->indexND(0+whichND($r2==$r2->min)->(:,(0)));

            

 

            my $coord = $c + $dc;

            # At last, we've found the best-fit to an enl_scale'th of an input-table pixel.

            # Back it out to input-science coordinates, and stuff it in the inverse table.

            $itable_flattened->(($i)) .= $coord;

            

            print " $i..." if( ($i%1000==0) && ( $PDL::verbose || $PDL::debug || $PDL::Transform::debug));

        }

        $itable->clump($itable->ndims-1) .= $itable_flattened;

        $itable->mv(-1,0) -= $p->{offset};

        $itable->mv(-1,0) /= $p->{scale};

 

        $p->{itable} = $itable;

      &$lookup_func($data,$p, $p->{itable},$p->{iscale},$p->{ioffset}) ;

  pop @_ if (($#_ % 2 ==0) && !defined($_[-1]));

  #suppresses a warning if @_ has an odd number of elements and the

  #last is undef

 

    # If CPi_j aren't present, and N=2, then try using CROTA or

        # Rotation matrix rotates counterclockwise to get from sci to pixel coords

        # (detector has been rotated ccw, according to FITS standard)

      $cpm .= pdl( [cos($cr), sin($cr)],[-sin($cr),cos($cr)] );

#If there are no CDELTs, then we assume they are all 1.0,

#as in PDL::Graphics::PGPLOT::Window::_FITS_tr.

    $cd+=1 if (all($cd==0));