Index: /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c =================================================================== --- /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c (revision 38323) +++ /branches/eam_branches/ohana.20150429/src/libfits/matrix/F_compress_M.c (revision 38323) @@ -0,0 +1,372 @@ +# include +# include +# include + +# define ESCAPE { \ + if (ztile != NULL) free (ztile); \ + if (optname != NULL) { \ + for (j = 0; j < Noptions; j++) { \ + free (optname[j]); \ + free (optvalue[j]); \ + } \ + free (optname); \ + free (optvalue); \ + } \ + if (out != NULL) free (out); \ + if (otile != NULL) free (otile); \ + if (ntile != NULL) free (ntile); \ + return (FALSE); } + +# define MOD_KEYWORD(ZNAME,NAME,TYPE,OUT,VALUE) { \ + gfits_modify (theader, ZNAME, TYPE, 1, OUT); \ + gfits_delete (theader, NAME, 1); } + +// the user needs to specify the various compression options: +// ztile[i] +// zcmptype + +int gfits_compress_image (Header *header, Matrix *matrix, FTable *ftable) { + + off_t Nzdata, Nzrows, zcol; + int i, j, status, zimage, Nout, max_tile_size; + char cmptype[80]; + char zaxis[10], naxis[10], key[10], word[81], exttype[81], checksum[81], datasum[81]; + int Noptions, NOPTIONS, zblank, oblank; + VarLengthColumn zdef; + float zscale, zzero; + + int zdata_pixsize, odata_pixsize; + + int *ztile = NULL; + int *otile = NULL; + int *ntile = NULL; + char **optname = NULL; + char **optvalue = NULL; + char *out = NULL; + char *zdata = NULL; + + Noptions = 0; + + // is ZIMAGE present? + // NOTE target of %t must be int length + status = gfits_scan_alt (ftable->header, "ZIMAGE", "%t", 1, &zimage); + if (!status || !zimage) ESCAPE; + + // copy original header to output header + gfits_copy_header (header, theader); + + // add ZIMAGE from output header + gfits_modify_alt (theader, "ZIMAGE", "%t", 1, TRUE); + + // translate image-specific keywords to Z versions (updates the Z version, removes the old version) + MOD_KEYWORD ("ZBITPIX", "BITPIX", "%d", header->bitpix); + MOD_KEYWORD ("ZNAXIS", "NAXIS", "%d", header->Naxes); + + for (i = 0; i < header->Naxes; i++) { + snprintf (zaxis, 10, "ZNAXIS%d", i + 1); + snprintf (naxis, 10, "NAXIS%d", i + 1); + MOD_KEYWORD (zaxis, naxis, OFF_T_FMT, header->Naxis[i]); + } + + for (i = 0; i < header->Naxes; i++) { + snprintf (key, 10, "ZTILE%d", i + 1); + if (!gfits_modify (theader, key, "%d", 1, ztile[i])) ESCAPE; + } + + // define compression-specific keywords, update header as needed. + if (!gfits_modify (header, "ZCMPTYPE", "%s", 1, zcmptype)) ESCAPE; + + // search for algorithm-specific keywords. these are used to control compression options + // note the difference between the keywords (1 indexed) and the variables (0 indexed) + for (i = 0; i < Noptions; i++) { + snprintf (key, 10, "ZNAME%d", i + 1); + if (!gfits_modify (header, key, "%s", 1, optname[i])) break; + + snprintf (key, 10, "ZVAL%d", i + 1); + if (!gfits_modify (header, key, "%s", 1, optvalue[i])) break; + } + +# if (0) + // XXX get ZMASKCMP + int zsimple, have_zsimple; + int have_ztension; + + // NOTE target of %t must be int length + // XXX how do I specify this?? + // have_zsimple = gfits_scan_alt (header, "ZSIMPLE", "%t", 1, &zsimple); + // have_ztension = gfits_scan (header, "ZTENSION", "%s", 1, exttype); + + // this is a very bogus case: we cannot have both keywords + if (have_zsimple && have_ztension) ESCAPE; + + // if neither are present, we have an image that is not really following the standard: + // assume it is a PHU + if (!have_zsimple && !have_ztension) { + header->simple = TRUE; + gfits_extended_to_primary (header, header->simple, "Image data"); + + MOD_KEYWORD_ALT ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); + MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); + } + + // have_zsimple : image comes from a PHU + if (have_zsimple) { + header->simple = zsimple; + gfits_delete (header, "ZSIMPLE", 1); + gfits_extended_to_primary (header, header->simple, "Image data"); + + MOD_KEYWORD_ALT ("ZEXTEND", "EXTEND", "%t", &header->extend, header->extend); + MOD_KEYWORD_ALT ("ZBLOCKED", "BLOCKED", "%t", &header->extend, header->extend); + } + + // have_ztension : image comes from an extension + if (have_ztension) { + gfits_delete (header, "ZTENSION", 1); + gfits_modify_extended (header, exttype, "Image extension"); + + // we may have an uncompressed PCOUNT / GCOUNT value, otherwise set to 0,1 + if (gfits_scan (header, "ZPCOUNT", "%d", 1, &header->pcount)) { + gfits_delete (header, "ZPCOUNT", 1); + gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + } else { + header->pcount = 0; + gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + } + if (gfits_scan (header, "ZGCOUNT", "%d", 1, &header->gcount)) { + gfits_delete (header, "ZGCOUNT", 1); + gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); + } else { + header->pcount = 1; + gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); + } + } else { + header->pcount = 0; + header->gcount = 1; + gfits_modify (header, "PCOUNT", "%d", 1, header->pcount); + gfits_modify (header, "GCOUNT", "%d", 1, header->gcount); + } + + MOD_KEYWORD ("ZHECKSUM", "CHECKSUM", "%s", checksum, checksum); + MOD_KEYWORD ("ZDATASUM", "DATASUM", "%s", datasum, datasum); +# endif + + zscale = 1; + gfits_scan (header, "ZSCALE", "%f", 1, &zscale); + + zblank = 32767; + gfits_scan (header, "ZBLANK", "%d", 1, &zblank); + + oblank = 32767; + gfits_scan (header, "BLANK", "%d", 1, &oblank); + + zzero = 0; + gfits_scan (header, "ZZERO", "%f", 1, &zzero); + + // XXX do not try to support compression of an image with an associated table (not valid) + + // create a table with a column called COMPRESSED_DATA + // output format may be 1PB, 1PI, 1PJ + + // find the COMPRESSED_DATA column (format should be 1PB, 1PI, 1PJ) + for (i = 1; TRUE; i++) { + snprintf (key, 10, "TTYPE%d", i); + if (!gfits_scan (ftable->header, key, "%s", 1, word)) ESCAPE; + if (!strcmp (word, "COMPRESSED_DATA")) break; + } + zcol = i; + + if (!gfits_varlength_column_define (ftable, &zdef, zcol)) ESCAPE; + gfits_delete (header, "TFIELDS", 1); + snprintf (key, 10, "TTYPE"OFF_T_FMT, zcol); + gfits_delete (header, key, 1); + snprintf (key, 10, "TFORM"OFF_T_FMT, zcol); + gfits_delete (header, key, 1); + + // counters for tile numbers + ALLOCATE (otile, int, matrix->Naxes); + ALLOCATE (ntile, int, matrix->Naxes); + max_tile_size = 1; + for (i = 0; i < matrix->Naxes; i++) { + otile[i] = 0; + ntile[i] = (matrix->Naxis[i] % ztile[i]) ? (matrix->Naxis[i] / ztile[i] + 1) : (matrix->Naxis[i] / ztile[i]); + max_tile_size *= ztile[i]; + + // ztile[i] is the default (or max) tile size in the i-th dimension + // ntile[i] is the number of tiles in the i-th dimension + // otile[i] is the current output tile counter in the i-th dimension + } + + // this takes place in three major steps: + // 1) rawdata -> tiledata (copy a tile chunk to the tiledata buffer) + // 2) compress tiledata -> zdata + // 3) add zdata to table data heap + + // heapdata.pixsize : zdef.format + // zdata.pixsize : depends on compression + // odata.pixsize : depends on compression + // idata.pixsize : header.bitpix + + // size of an element in the vartable heap section + zdata_pixsize = gfits_vartable_heap_pixsize (zdef.format); + + // size of a pixel in the final image (not needed) + // idata_pixsize = abs(header[0].bitpix) / 8; + + // size of a pixel in the output from the decompression routine + raw_pixsize = gfits_uncompressed_data_pixsize (cmptype, header[0].bitpix, optname, optvalue, Noptions) + ALLOCATE (raw_buffer, char, raw_pixsize*max_tile_size); + ALLOCATE (zdata, char, raw_pixsize*max_tile_size); // + + // compress the data : copy into a tile, compress the tile, then add to the output table + // each tile -> 1 row of the output table + int done = FALSE; + while (!done) { + + // size of the current tile + Nraw = gfits_tile_size (matrix, otile, ztile); + + // copy the raw pixels from their native matrix locations to the temporary output buffer + if (!gfits_collect_data (matrix, raw, Nraw, raw_pixsize, otile, oblank, ztile, zblank, zscale, zzero)) return (FALSE); + + // XXX the tile must not be > 2GB + + Nzdata = raw_pixsize*max_tile_size; // available space, replaced with actual output size on compression + if (!gfits_compress_data (zdata, &Nzdata, cmptype, optname, optvalue, Noptions, raw, Nraw, raw_pixsize)) return (FALSE); + + if (!gfits_varlength_column_addrow (ftable, zdata, Nzdata, &zdef)) return FALSE; + + // XXX need to activate + if (FALSE) gfits_byteswap_zdata (zdata, Nzdata, zdata_pixsize) return (FALSE); + + // update the tile counters, carrying to the next dimension if needed + for (j = 0; j < matrix->Naxes; j++) { + otile[j] ++; + if (otile[j] == ntile[j]) { + if (j == matrix->Naxes - 1) { + done = TRUE; + break; + } + otile[j] = 0; + } else { + break; + } + } + } + return (TRUE); +} + +// raw_pixsize is the bytes / pixel for the input matrix +// place the raw image bytes for the current tile into the tile buffer +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) { + + int i, j, start, offset, coord, Nline; + int *counter = NULL; + int *Ztile = NULL; + + ALLOCATE (counter, int, matrix->Naxes); // counter for current row in tile to copy + ALLOCATE (Ztile, int, matrix->Naxes); // true sizes of this tile (in pixels) + + for (i = 0; i < matrix->Naxes; i++) { + counter[i] = 0; + Ztile[i] = MIN ((matrix->Naxis[i] - otile[i]*ztile[i]), ztile[i]); + } + + // number of lines in the tile (in pixels) + Nline = 1; + for (i = 1; i < matrix->Naxes; i++) { + Nline *= Ztile[i]; + } + + // confirm the tile size: + assert (Nraw == Ztile[0]*Nline); + + // set the starting point of the tile in the matrix buffer + // start = otile[0]*ztile[0] + otile[1]*ztile[1]*Naxis[0] + otile[2]*ztile[2]*Naxis[0]*Naxis[1] + ...; + // start = otile[0]*ztile[0] + Naxis[0]*(otile[1]*ztile[1] + Naxis[1]*(otile[2]*ztile[2] + ...)); + start = otile[matrix->Naxes-1]*ztile[matrix->Naxes-1]; + for (i = matrix->Naxes - 2; i >= 0; i--) { + coord = otile[i]*ztile[i]; + start = start*matrix->Naxis[i] + coord; + } + + // pixel offset in output array relative to tile start + offset = 0; + + // we need to set up switches for all of the possible combinations: + // this macro is used at the inner switch to run the actual loop + // note we copy Optr to Iptr +# define SCALE_AND_DIST(TYPE, SIZE) { \ + TYPE *Optr = (TYPE *) &matrix->buffer[SIZE*(offset + start)]; \ + for (j = 0; j < Ztile[0]; j++, Iptr++, Optr++) { \ + if (*Optr == oblank) { \ + *Iptr = zblank; \ + } else { \ + *Iptr = (*Optr - zzero) / zscale; \ + } } } + + // this macro sets up the outer switch and calls above macro with all output bitpix options +# define SETUP_INSIZE(TYPE, SIZE) { \ + TYPE *Iptr = (TYPE *) &raw[i*SIZE*Ztile[0]]; \ + switch (matrix->bitpix) { \ + case 8: \ + SCALE_AND_DIST (char, 1); \ + break; \ + case 16: \ + SCALE_AND_DIST (short, 2); \ + break; \ + case 32: \ + SCALE_AND_DIST (int, 4); \ + break; \ + case -32: \ + SCALE_AND_DIST (float, 4); \ + break; \ + case -64: \ + SCALE_AND_DIST (double, 8); \ + break; \ + default: \ + abort(); \ + } } + + // loop over lines in the tile + for (i = 0; i < Nline; i++) { + switch (raw_pixsize) { + case 1: + SETUP_INSIZE (char, 1); + break; + case 2: + SETUP_INSIZE (short, 2); + break; + case 4: + SETUP_INSIZE (int, 4); + break; + default: + abort(); + } + + // update the counters, carrying to the next dimension if needed + for (j = 1; j < matrix->Naxes; j++) { + counter[j] ++; + if (counter[j] == Ztile[j]) { + counter[j] = 0; + } else { + break; + } + } + if (j == matrix->Naxes) assert (i == Nline - 1); // we should be done here... + + // Naxes = 3 + // offset = counter[1]*matrix->Naxis[0] + counter[2]*matrix->Naxis[0]*matrix->Naxis[1] + + // offset = matrix->Naxis[0]*(counter[1] + matrix->Naxis[1]*(counter[2] + matrix->Naxis[2]*...)) + + // determine the offset of the next line relative to the start position + offset = counter[matrix->Naxes - 1]; + for (j = matrix->Naxes - 2; j >= 0; j--) { + offset = offset*matrix->Naxis[j] + counter[j]; + } + } + + free (counter); + free (Ztile); + + return (TRUE); +}