- Timestamp:
- Jun 3, 2015, 4:10:26 AM (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c
r38340 r38355 18 18 // we do not try to support compression of an image with an associated table (not valid) 19 19 20 # define ESCAPE (A){ fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; }20 # define ESCAPE { fprintf (stderr, "error in %s @ line %d\n", __func__, __LINE__); goto escape; } 21 21 22 22 int gfits_compress_image (Header *header, Matrix *matrix, FTable *ftable, int *Ztile, char *zcmptype) { … … 46 46 47 47 // creates an empty (Naxes = 2, Naxis[i] = 0) table header with XTENSION = BINTABLE, EXTNAME = COMPRESSED_IMAGE 48 if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE (A);48 if (!gfits_create_table_header (theader, "BINTABLE", "COMPRESSED_IMAGE")) ESCAPE; 49 49 50 50 // by using "P" format, we are limited to 32bit pointers (2GB heap) 51 51 // XXX how is the "B" format used? 52 if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE (A);53 if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE (A);52 if (!gfits_define_bintable_column (theader, "1PB(0)", "COMPRESSED_DATA", "compressed image data", "none", 1.0, 0.0)) ESCAPE; 53 if (!gfits_delete (theader, "TUNIT1", 1)) ESCAPE; 54 54 55 55 // allocates the default data array (ftable->buffer) 56 if (!gfits_create_table (theader, ftable)) ESCAPE (A);56 if (!gfits_create_table (theader, ftable)) ESCAPE; 57 57 58 58 // copy original header to output header (XXX this needs to be finished) 59 if (!gfits_copy_keywords_compress (header, theader)) ESCAPE (A);59 if (!gfits_copy_keywords_compress (header, theader)) ESCAPE; 60 60 61 61 // this allocates the pointer data array for the full set of tiles (NOT the heap) 62 62 char *tmpbuffer = NULL; 63 63 ALLOCATE_ZERO (tmpbuffer, char, Ntile); 64 if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE (A);64 if (!gfits_set_bintable_column (theader, ftable, "COMPRESSED_DATA", tmpbuffer, Ntile)) ESCAPE; 65 65 free (tmpbuffer); 66 66 67 67 // add ZIMAGE from output header 68 if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE (A);68 if (!gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE)) ESCAPE; 69 69 70 70 // define compression-specific keywords, update header as needed. 71 if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE (A);72 73 if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE (A);71 if (!gfits_modify (theader, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE; 72 73 if (!gfits_modify_alt (theader, "ZSIMPLE", "%t", 1, TRUE)) ESCAPE; 74 74 // gfits_modify (header, "ZTENSION", "%s", 1, exttype); 75 75 76 76 // supply the Z* header values from the source image header 77 if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE (A);78 79 if (!gfits_modify (theader, "ZNAXIS", "%d", 1, header[0].Naxes)) ESCAPE (A);77 if (!gfits_modify (theader, "ZBITPIX", "%d", 1, header[0].bitpix)) ESCAPE; 78 79 if (!gfits_modify (theader, "ZNAXIS", "%d", 1, header[0].Naxes)) ESCAPE; 80 80 81 81 char keyword[11]; 82 82 for (i = 0; i < header[0].Naxes; i++) { 83 83 snprintf (keyword, 10, "ZNAXIS%d", i + 1); 84 if (!gfits_modify (theader, keyword, OFF_T_FMT, 1, header[0].Naxis[i])) ESCAPE (A);84 if (!gfits_modify (theader, keyword, OFF_T_FMT, 1, header[0].Naxis[i])) ESCAPE; 85 85 } 86 86 87 87 for (i = 0; i < header->Naxes; i++) { 88 88 snprintf (keyword, 10, "ZTILE%d", i + 1); 89 if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE (A);90 } 91 92 if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE (A);93 if (!gfits_modify (theader, "BZERO", "%lf", 1, header[0].bzero)) ESCAPE (A);89 if (!gfits_modify (theader, keyword, "%d", 1, ztile[i])) ESCAPE; 90 } 91 92 if (!gfits_modify (theader, "BSCALE", "%lf", 1, header[0].bscale)) ESCAPE; 93 if (!gfits_modify (theader, "BZERO", "%lf", 1, header[0].bzero)) ESCAPE; 94 94 95 95 // other keywords to define: … … 128 128 129 129 VarLengthColumn zdef; 130 if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE(A); 131 132 // allocate the working buffers to the max tile size 130 if (!gfits_varlength_column_define (ftable, &zdef, 1)) ESCAPE; 131 132 // pixel sizes and tile sizes: 133 // * We have Ntile tiles, each of Nx * Ny * Nz ... tile pixels 134 // -> A pixel in the tile (in the image) has size tile_pixsize = f(BITPIX) 135 // * We copy the pixels in the tile to a continuous buffer to be compressed. 136 // During this copy, we can scale the data (eg, floats -> ints), in which 137 // case the raw buffer may have a different pixel size (and type) than the 138 // image tile 139 // -> A pixel in the raw buffer is f(zcmptype,BITPIX): 140 // -- GZIP_1 : raw_pixsize = tile_pixsize (no scaling is performed) 141 // -- GZIP_2 : raw_pixsize = tile_pixsize (no scaling is performed) 142 // -- RICE_1 : raw pixels must be integers. raw_pixsize is user specified? (ZNAMEn,ZVALn : BYTEPIX = 1,2,4,8) 143 // -- PLIO_1 : raw_pixsize = 2 bytes 144 // -- HCOMPRESS_1 : raw pixels must be integers. 145 146 // size (in pixels) of the largest tile 133 147 int max_tile_size = gfits_imtile_maxsize (matrix, ztile); 134 148 135 // size of a pixel in the raw image tile: 136 int raw_pixsize = abs(header[0].bitpix) / 8; 137 138 // size of a pixel in the raw image tile: 139 int tile_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); 140 149 // size of a pixel in the image tile -- this is probably not needed 150 int tile_pixsize = abs(header[0].bitpix) / 8; 151 152 // size of a pixel in the raw pixel buffer (before compression) 153 int raw_pixsize = gfits_uncompressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); 154 int raw_bitpix = gfits_uncompressed_data_bitpix (zcmptype, header[0].bitpix, NULL, NULL, 0); 155 int cmp_pixsize = gfits_compressed_data_pixsize (zcmptype, header[0].bitpix, NULL, NULL, 0); 156 157 fprintf (stderr, "raw_pixsize: %d, cmp_pixsize: %d, tile_pixsize: %d, raw_bitpix: %d\n", 158 raw_pixsize, cmp_pixsize, tile_pixsize, raw_bitpix); 159 160 // allocate the buffer for compression work 141 161 int Nzdata_alloc = raw_pixsize*max_tile_size + 100; 142 162 ALLOCATE (raw, char, raw_pixsize*max_tile_size); … … 144 164 145 165 // init the otile[] counters 146 if (!gfits_imtile_start (matrix, otile)) ESCAPE (A);166 if (!gfits_imtile_start (matrix, otile)) ESCAPE; 147 167 148 168 // compress the data : copy into a tile, compress the tile, then add to the output table … … 150 170 for (i = 0; i < Ntile; i++) { 151 171 152 // size of the current tile 172 // size of the current tile in pixels 153 173 int Nraw = gfits_tile_size (matrix, otile, ztile); 154 174 155 175 // copy the raw pixels from their native matrix locations to the temporary output buffer 156 if (!gfits_collect_data (matrix, raw, Nraw, raw_pixsize, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE(A); 157 176 // for float -> int scaling by zscale, zzero may be applied 177 if (!gfits_collect_data (matrix, raw, Nraw, raw_bitpix, otile, oblank, ztile, zblank, zscale, zzero)) ESCAPE; 178 179 if (0 && (i == 0)) { 180 int k; 181 fprintf (stderr, "cmp mat: "); 182 for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", matrix->buffer[k]); } 183 fprintf (stderr, "\n"); 184 fprintf (stderr, "cmp raw: "); 185 for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", raw[k]); } 186 fprintf (stderr, "\n"); 187 } 188 189 // gzip compresses bytes which are in BIG-ENDIAN order 158 190 if (!strcasecmp(zcmptype, "GZIP_1")) { 159 if (!gfits_byteswap_zdata (raw, Nraw, raw_pixsize)) ESCAPE(A); 160 } 161 162 // XXX the tile must not be > 2GB 163 164 // optname, optvalue = NULL, Noptions = 0 191 if (!gfits_byteswap_zdata (raw, Nraw, raw_pixsize)) ESCAPE; 192 } 193 if (0 && (i == 0)) { 194 int k; 195 fprintf (stderr, "cmp swp: "); 196 for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", raw[k]); } 197 fprintf (stderr, "\n"); 198 } 199 200 // optname, optvalue = NULL, Noptions = 0 : defaults for RICE_1, HCOMPRESS_1 165 201 166 202 int Nzdata = Nzdata_alloc; // available space, replaced with actual output size on compression 167 if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, tile_pixsize)) ESCAPE(A); 168 169 if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE(A); 170 171 // XXX need to activate 172 // if (FALSE) gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize) ESCAPE(A); 203 if (!gfits_compress_data (zdata, &Nzdata, zcmptype, NULL, NULL, 0, raw, Nraw, raw_pixsize, ztile[1], ztile[0])) ESCAPE; 204 if (0 && (i == 0)) { 205 int k; 206 fprintf (stderr, "cmp dat: "); 207 for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", zdata[k]); } 208 fprintf (stderr, "\n"); 209 } 210 211 if (!gfits_varlength_column_add_data (ftable, zdata, Nzdata, i, &zdef)) ESCAPE; 212 213 // all other compression modes require swapping after compression 214 // (compresssion & decompression operate on native ENDIAN) 215 if (strcasecmp(zcmptype, "GZIP_1")) { 216 if (!gfits_byteswap_zdata (zdata, Nzdata, cmp_pixsize)) ESCAPE; 217 } 173 218 174 219 // update the otile[] counters, carrying to the next dimension if needed … … 176 221 } 177 222 178 if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE(A); 223 if (!gfits_varlength_column_finish (ftable, &zdef)) ESCAPE; 224 if (0) { 225 int k; 226 fprintf (stderr, "cmp tbl: "); 227 for (k = 0; k < 32; k++) { fprintf (stderr, "0x%02hhx ", ftable->buffer[k]); } 228 fprintf (stderr, "\n"); 229 } 179 230 180 231 free (raw); … … 209 260 // raw_pixsize is the bytes / pixel for the input matrix 210 261 // place the raw image bytes for the current tile into the tile buffer 211 int gfits_collect_data (Matrix *matrix, char *raw, int Nraw, int raw_pixsize, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) { 262 // * otile is the counter for the current output tile 263 // * ztile gives the size of the i-th dimension of the current tile 264 // * raw_bitpix defines the size and type of the raw buffer 265 int gfits_collect_data (Matrix *matrix, char *raw, int Nraw, int raw_bitpix, int *otile, int oblank, int *ztile, int zblank, float zscale, float zzero) { 212 266 213 267 int i, j, start, offset, coord, Nline; … … 246 300 // we need to set up switches for all of the possible combinations: 247 301 // this macro is used at the inner switch to run the actual loop 248 // note we copy Optr to Iptr 249 # define SCALE_AND_DIST(TYPE, SIZE) { \ 250 TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ 251 for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ 252 if (*Optr == oblank) { \ 253 *Iptr = zblank; \ 302 # define SCALE_AND_DIST_INT(TYPE, SIZE) { \ 303 TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ 304 for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ 305 if (*TILEptr == oblank) { \ 306 *RAWptr = zblank; \ 254 307 } else { \ 255 * Iptr = (*Optr - zzero) / zscale; \308 *RAWptr = (*TILEptr - zzero) / zscale; \ 256 309 } } } 257 310 258 311 // we need to set up switches for all of the possible combinations: 259 312 // this macro is used at the inner switch to run the actual loop 260 // note we copy Optr to Iptr 261 # define SCALE_AND_DIST_FLOAT(TYPE, SIZE) { \ 262 TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ 263 for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ 264 if (!isfinite(*Optr)) { \ 265 *Iptr = zblank; \ 313 # define SCALE_AND_DIST_FLOAT(TYPE, SIZE) { \ 314 TYPE *TILEptr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ 315 for (j = 0; j < Ztile[0]; j++, TILEptr++, RAWptr++) { \ 316 if (!isfinite(*TILEptr)) { \ 317 *RAWptr = zblank; \ 266 318 } else { \ 267 * Iptr = (*Optr - zzero) / zscale; \319 *RAWptr = (*TILEptr - zzero) / zscale; \ 268 320 } } } 269 321 270 // this macro sets up the outer switch and calls above macro with all output bitpix options 271 # define SETUP_INSIZE(TYPE, SIZE) { \ 272 TYPE *Iptr = (TYPE *) &raw[i*SIZE*Ztile[0]]; \ 322 // this macro sets up the outer switch and calls above the macro with all RAW buffer bitpix options 323 // 324 # define SETUP_RAWSIZE(TYPE, SIZE) { \ 325 TYPE *RAWptr = (TYPE *) &raw[i*SIZE*Ztile[0]]; \ 273 326 switch (matrix->bitpix) { \ 274 case 8: \ 275 SCALE_AND_DIST (char, 1); \ 276 break; \ 277 case 16: \ 278 SCALE_AND_DIST (short, 2); \ 279 break; \ 280 case 32: \ 281 SCALE_AND_DIST (int, 4); \ 282 break; \ 283 case -32: \ 284 SCALE_AND_DIST_FLOAT (float, 4); \ 285 break; \ 286 case -64: \ 287 SCALE_AND_DIST_FLOAT (double, 8); \ 288 break; \ 289 default: \ 290 abort(); \ 327 case 8: SCALE_AND_DIST_INT (char, 1); break; \ 328 case 16: SCALE_AND_DIST_INT (short, 2); break; \ 329 case 32: SCALE_AND_DIST_INT (int, 4); break; \ 330 case -32: SCALE_AND_DIST_FLOAT (float, 4); break; \ 331 case -64: SCALE_AND_DIST_FLOAT (double, 8); break; \ 332 default: abort(); \ 291 333 } } 292 334 293 335 // loop over lines in the tile 294 336 for (i = 0; i < Nline; i++) { 295 switch (raw_pixsize) { 296 case 1: 297 SETUP_INSIZE (char, 1); 298 break; 299 case 2: 300 SETUP_INSIZE (short, 2); 301 break; 302 case 4: 303 SETUP_INSIZE (int, 4); 304 break; 305 case 8: 306 SETUP_INSIZE (double, 8); 307 break; 308 default: 309 abort(); 337 switch (raw_bitpix) { 338 case 8: SETUP_RAWSIZE (char, 1); break; 339 case 16: SETUP_RAWSIZE (short, 2); break; 340 case 32: SETUP_RAWSIZE (int, 4); break; 341 case -32: SETUP_RAWSIZE (float, 4); break; 342 case -64: SETUP_RAWSIZE (double, 8); break; 343 default: abort(); 310 344 } 311 345
Note:
See TracChangeset
for help on using the changeset viewer.
