#!/usr/bin/perl
#
use Statistics::OLS;
use Math::FFT;
use Math::Stat;
use Math::Spline qw(spline linsearch binsearch);
use Math::Derivative qw(Derivative2);
use Math::Approx;
use Astro::FITS::CFITSIO qw( :longnames :constants );
use PDL;
use PDL::Fit::Polynomial; 
use PDL::Filter::Linear;
use PGPLOT;  # Load PGPLOT module
use PDL::Fit::Gaussian;
use PDL::Slatec;
use PDL::Image2D;
#use PDL::Matrix;

#require "/work1/ssa/perl/MY/my.pl";
#require "/home/sanchez/3DSoft_Ubuntu/R3D/my.pl";
require "/home/sanchez/sda2/code/R3D/my.pl";

if ($#ARGV<2) {
    print "USE: index_seg_cube.pl indices.out SEGMENTATION.fits indeces.cube.FITS\n";
    exit;
}

$infile=$ARGV[0];
$segfile=$ARGV[1];
$outfile=$ARGV[2];

$nseg=0;
$n_index=0;
open(FH,"<$infile");
while($line=<FH>) {
    chop($line);
    @data=split(" ",$line);
    if ($nseg==0) {
	$n_index=($#data+1)/3;
	for ($i=0;$i<$n_index;$i++) {
	    $name[$i]=$data[3*$i];
	    $name[$i+$n_index]="e_".$data[3*$i];
	}
    }
    if ((3*$n_index)==($#data+1)) {
#    print "(3*$n_index==($#data+1)\n";
	for ($i=0;$i<$n_index;$i++) {
	    $INDEX[$i][$nseg]=$data[3*$i+1];
	    $e_INDEX[$i][$nseg]=$data[3*$i+2];
	}
    } else {
	for ($i=0;$i<$n_index;$i++) {
	    $INDEX[$i][$nseg]=0;
	    $e_INDEX[$i][$nseg]=0;
	}
    }
    $nseg++;
}
close(FH);
#exit;


$s_in=rfits($segfile);
($nx,$ny)=$s_in->dims;
$nz=2*$n_index;
$pdl_cube=zeroes($nx,$ny,$nz);




for ($i=0;$i<$nx;$i++) {
    for ($j=0;$j<$ny;$j++) {
	$is=$s_in->at($i,$j);
	$is_out=$is-1;
	if ($is_out>=0) {
	    for ($k=0;$k<$n_index;$k++) {
		$kk=$k+$n_index;
#		print "$i/$nx , $j/$ny $k/$nz $kk/$nz $INDEX[$k][$is_out] $e_INDEX[$k][$is_out]\n";
		$val=$INDEX[$k][$is_out];
		$e_val=$e_INDEX[$k][$is_out];
		set($pdl_cube,($i,$j,$k),$val);
		set($pdl_cube,($i,$j,$kk),$val);
	    }
	}
    }
}


$h = {NAXIS=>3, NAXIS1=>$nx, NAXIS2=>$ny, NAXIS3=>$nz, COMMENT=>"FIT-header"};
$$h{FILENAME} = $outfile;
for ($k=0;$k<$n_index;$k++) {
    $index="INDEX".$k;
    $$h{$index}=$name[$k];
    $kk=$k+$n_index;
    $index="INDEX".$kk;
    $$h{$index}=$name[$kk];
}
$pdl_cube->sethdr($h);
$pdl_cube->wfits($outfile);


exit;