Open mohawk2 opened 2 months ago
Also, reading back a written-out-compressed byte
array is wrong in at least one circumstance (RICE_1 is documented as lossless, https://www.irjet.net/archives/V3/i2/IRJET-V3I2267.pdf), and also can't handle 3D (inc RGB) data:
--- a/IO/FITS/t/fits.t
+++ b/IO/FITS/t/fits.t
@@ -345,6 +345,19 @@ ok all(approx $m51, $m51_2), 'read back written-out bintable FITS file' or diag
$m51->wfits($fname, {compress=>1});
$m51_2 = rfits($fname);
ok all(approx $m51, $m51_2), 'read back written-out compressed FITS file' or diag "got:", $m51_2->info;
+$m51_2->hdrcpy(1);
+$m51_2 = $m51_2->dummy(2,3)->sever;
+$m51_2->hdr->{NAXIS} = 3;
+$m51_2->hdr->{NAXIS3} = 3;
+eval {$m51_2->wfits($fname, {compress=>1})};
+is $@, '', 'wfits can compress RGB image';
+if (eval 'use PDL::IO::Pic; 1') {
+my $shape = rpic(cfile(($fs->updir)x2, qw(Libtmp Transform Cartography earth_height-2048x1024.jpg)));
+$shape->wfits($fname, {compress=>1});
+my $shape2 = rfits($fname);
+ok all($shape == $shape2), 'rfits compress of byte OK'
+ or diag 'first non-matching indices', +($shape != $shape2)->whichND->slice(',0');
+}
}
}
The current CFITSIO rcomp
etc is in https://github.com/HEASARC/cfitsio/blob/develop/ricecomp.c, vs our copy https://github.com/PDLPorters/pdl/blob/master/Libtmp/Compression/ricecomp.c
Some issues on the HEASARC repo remind me that astropy has FITS capability (https://docs.astropy.org/en/stable/io/fits/index.html) so we can use that to compare behaviour.
The above-linked commit fixes the decompression error. This was in fact a compression error, which was apparently fixed in CFITSIO in the meantime so I have just updated us to using their latest.
The above-linked commit fixes the 3D compression problem, by using the same strategy it looks like fpack
etc do: clump higher data dims.
See https://github.com/astropy/astropy/issues/6778 for how to do complex numbers (spoiler: use the 0th dim for real/imaginary).