jj2000.j2k.entropy.encoder.EBCOTRateAllocator Maven / Gradle / Ivy
Show all versions of jai-imageio-jpeg2000 Show documentation
/*
* $RCSfile: EBCOTRateAllocator.java,v $
* $Revision: 1.1 $
* $Date: 2005/02/11 05:02:08 $
* $State: Exp $
*
* Class: EBCOTRateAllocator
*
* Description: Generic interface for post-compression
* rate allocator.
*
*
*
* COPYRIGHT:
*
* This software module was originally developed by Raphaël Grosbois and
* Diego Santa Cruz (Swiss Federal Institute of Technology-EPFL); Joel
* Askelöf (Ericsson Radio Systems AB); and Bertrand Berthelot, David
* Bouchard, Félix Henry, Gerard Mozelle and Patrice Onno (Canon Research
* Centre France S.A) in the course of development of the JPEG2000
* standard as specified by ISO/IEC 15444 (JPEG 2000 Standard). This
* software module is an implementation of a part of the JPEG 2000
* Standard. Swiss Federal Institute of Technology-EPFL, Ericsson Radio
* Systems AB and Canon Research Centre France S.A (collectively JJ2000
* Partners) agree not to assert against ISO/IEC and users of the JPEG
* 2000 Standard (Users) any of their rights under the copyright, not
* including other intellectual property rights, for this software module
* with respect to the usage by ISO/IEC and Users of this software module
* or modifications thereof for use in hardware or software products
* claiming conformance to the JPEG 2000 Standard. Those intending to use
* this software module in hardware or software products are advised that
* their use may infringe existing patents. The original developers of
* this software module, JJ2000 Partners and ISO/IEC assume no liability
* for use of this software module or modifications thereof. No license
* or right to this software module is granted for non JPEG 2000 Standard
* conforming products. JJ2000 Partners have full right to use this
* software module for his/her own purpose, assign or donate this
* software module to any third party and to inhibit third parties from
* using this software module for non JPEG 2000 Standard conforming
* products. This copyright notice must be included in all copies or
* derivative works of this software module.
*
* Copyright (c) 1999/2000 JJ2000 Partners.
* */
package jj2000.j2k.entropy.encoder;
import java.awt.Point;
import java.io.IOException;
import jj2000.j2k.codestream.Markers;
import jj2000.j2k.codestream.PrecInfo;
import jj2000.j2k.codestream.ProgressionType;
import jj2000.j2k.codestream.writer.BitOutputBuffer;
import jj2000.j2k.codestream.writer.CodestreamWriter;
import jj2000.j2k.codestream.writer.PktEncoder;
import jj2000.j2k.entropy.Progression;
import jj2000.j2k.util.FacilityManager;
import jj2000.j2k.util.MathUtil;
import jj2000.j2k.util.MsgLogger;
import jj2000.j2k.util.ProgressWatch;
import jj2000.j2k.wavelet.analysis.SubbandAn;
import com.github.jaiimageio.jpeg2000.impl.J2KImageWriteParamJava;
/**
* This implements the EBCOT post compression rate allocation algorithm. This
* algorithm finds the most suitable truncation points for the set of
* code-blocks, for each layer target bitrate. It works by first collecting
* the rate distortion info from all code-blocks, in all tiles and all
* components, and then running the rate-allocation on the whole image at
* once, for each layer.
*
* This implementation also provides some timing features. They can be
* enabled by setting the 'DO_TIMING' constant of this class to true and
* recompiling. The timing uses the 'System.currentTimeMillis()' Java API
* call, which returns wall clock time, not the actual CPU time used. The
* timing results will be printed on the message output. Since the times
* reported are wall clock times and not CPU usage times they can not be added
* to find the total used time (i.e. some time might be counted in several
* places). When timing is disabled ('DO_TIMING' is false) there is no penalty
* if the compiler performs some basic optimizations. Even if not the penalty
* should be negligeable.
*
* @see PostCompRateAllocator
*
* @see CodedCBlkDataSrcEnc
*
* @see jj2000.j2k.codestream.writer.CodestreamWriter
* */
public class EBCOTRateAllocator extends PostCompRateAllocator {
/** Whether to collect timing information or not: false. Used as a compile
* time directive. */
private final static boolean DO_TIMING = false;
/** The wall time for the initialization. */
private long initTime;
/** The wall time for the building of layers. */
private long buildTime;
/** The wall time for the writing of layers. */
private long writeTime;
/**
* 5D Array containing all the coded code-blocks:
*
*
* - 1st index: tile index
* - 2nd index: component index
* - 3rd index: resolution level index
* - 4th index: subband index
* - 5th index: code-block index
*
**/
private CBlkRateDistStats cblks[][][][][];
/**
* 6D Array containing the indices of the truncation points. It actually
* contains the index of the element in CBlkRateDistStats.truncIdxs that
* gives the real truncation point index.
*
*
* - 1st index: tile index
* - 2nd index: layer index
* - 3rd index: component index
* - 4th index: resolution level index
* - 5th index: subband index
* - 6th index: code-block index
*
**/
private int truncIdxs[][][][][][];
/**
* Maximum number of precincts :
*
*
* - 1st dim: tile index.
* - 2nd dim: component index.
* - 3nd dim: resolution level index.
*
*/
private Point numPrec[][][];
/** Array containing the layers information. */
private EBCOTLayer layers[];
/** The log of 2, natural base */
private static final double LOG2 = Math.log(2);
/** The normalization offset for the R-D summary table */
private static final int RD_SUMMARY_OFF = 24;
/** The size of the summary table */
private static final int RD_SUMMARY_SIZE = 64;
/** The relative precision for float data. This is the relative tolerance
* up to which the layer slope thresholds are calculated. */
private static final float FLOAT_REL_PRECISION = 1e-4f;
/** The precision for float data type, in an absolute sense. Two float
* numbers are considered "equal" if they are within this precision. */
private static final float FLOAT_ABS_PRECISION = 1e-10f;
/**
* Minimum average size of a packet. If layer has less bytes than the this
* constant multiplied by number of packets in the layer, then the layer
* is skipped.
* */
private static final int MIN_AVG_PACKET_SZ = 32;
/**
* The R-D summary information collected from the coding of all
* code-blocks. For each entry it contains the accumulated length of all
* truncation points that have a slope not less than
* '2*(k-RD_SUMMARY_OFF)', where 'k' is the entry index.
*
* Therefore, the length at entry 'k' is the total number of bytes of
* code-block data that would be obtained if the truncation slope was
* chosen as '2*(k-RD_SUMMARY_OFF)', without counting the overhead
* associated with the packet heads.
*
*
This summary is used to estimate the relation of the R-D slope to
* coded length, and to obtain absolute minimums on the slope given a
* length.
**/
private int RDSlopesRates[];
/** Packet encoder. */
private PktEncoder pktEnc;
/** The layer specifications */
private LayersInfo lyrSpec;
/** The maximum slope accross all code-blocks and truncation points. */
private float maxSlope;
/** The minimum slope accross all code-blocks and truncation points. */
private float minSlope;
/**
* Initializes the EBCOT rate allocator of entropy coded data. The layout
* of layers, and their bitrate constraints, is specified by the 'lyrs'
* parameter.
*
* @param src The source of entropy coded data.
*
* @param lyrs The layers layout specification.
*
* @param writer The bit stream writer.
*
* @see ProgressionType
* */
public EBCOTRateAllocator(CodedCBlkDataSrcEnc src, LayersInfo lyrs,
CodestreamWriter writer,
J2KImageWriteParamJava wp) {
super(src,lyrs.getTotNumLayers(),writer,wp);
int minsbi, maxsbi;
int i;
SubbandAn sb, sb2;
Point ncblks = null;
// If we do timing create necessary structures
if (DO_TIMING) {
// If we are timing make sure that 'finalize' gets called.
System.runFinalizersOnExit(true);
// The System.runFinalizersOnExit() method is deprecated in Java
// 1.2 since it can cause a deadlock in some cases. However, here
// we use it only for profiling purposes and is disabled in
// production code.
initTime = 0L;
buildTime = 0L;
writeTime = 0L;
}
// Save the layer specs
lyrSpec = lyrs;
//Initialize the size of the RD slope rates array
RDSlopesRates = new int[RD_SUMMARY_SIZE];
//Get number of tiles, components
int nt = src.getNumTiles();
int nc = getNumComps();
//Allocate the coded code-blocks and truncation points indexes arrays
cblks = new CBlkRateDistStats[nt][nc][][][];
truncIdxs = new int[nt][numLayers][nc][][][];
int cblkPerSubband; // Number of code-blocks per subband
int mrl; // Number of resolution levels
int l; // layer index
int s; //subband index
// Used to compute the maximum number of precincts for each resolution
// level
int tx0, ty0, tx1, ty1; // Current tile position in the reference grid
int tcx0, tcy0, tcx1, tcy1; // Current tile position in the domain of
// the image component
int trx0, try0, trx1, try1; // Current tile position in the reduced
// resolution image domain
int xrsiz, yrsiz; // Component sub-sampling factors
Point tileI = null;
Point nTiles = null;
int xsiz,ysiz,x0siz,y0siz;
int xt0siz,yt0siz;
int xtsiz,ytsiz;
int cb0x = src.getCbULX();
int cb0y = src.getCbULY();
src.setTile(0,0);
for (int t=0; ttrx0) {
numPrec[t][c][r].x = (int)Math.ceil((trx1-cb0x)/twoppx)
- (int)Math.floor((trx0-cb0x)/twoppx);
} else {
numPrec[t][c][r].x = 0;
}
if (try1>try0) {
numPrec[t][c][r].y = (int)Math.ceil((try1-cb0y)/twoppy)
- (int)Math.floor((try0-cb0y)/(double)twoppy);
} else {
numPrec[t][c][r].y = 0;
}
minsbi = (r==0) ? 0 : 1;
maxsbi = (r==0) ? 1 : 4;
cblks[t][c][r] = new CBlkRateDistStats[maxsbi][];
for(l=0; l=0; n--) {
layers[n] = new EBCOTLayer();
}
minlsz = 0; // To keep compiler happy
for( int t=0 ; t totenclength)
nextbytes = totenclength;
}
else {
nextbytes = 1;
}
loopnlyrs = lyrSpec.getExtraLayers(i)+1;
ls = Math.exp(Math.log((double)nextbytes/basebytes)/loopnlyrs);
layers[n].optimize = true;
for (l = 0; l < loopnlyrs; l++) {
newbytes = (int)basebytes - lastbytes - ho;
if (newbytes < minlsz) { // Skip layer (too small)
basebytes *= ls;
numLayers--;
continue;
}
lastbytes = (int)basebytes - ho;
layers[n].maxBytes = lastbytes;
basebytes *= ls;
n++;
}
i++; // Goto next optimization point
}
// Ensure minimum size of last layer (this one determines overall
// bitrate)
n = numLayers-2;
nextbytes = (int) (lyrSpec.getTotBitrate()*np) - ho;
newbytes = nextbytes - ((n>=0) ? layers[n].maxBytes : 0);
while (newbytes < minlsz) {
if (numLayers == 1) {
if (newbytes <= 0) {
throw new
IllegalArgumentException("Overall target bitrate too "+
"low, given the current "+
"bit stream header overhead");
}
break;
}
// Delete last layer
numLayers--;
n--;
newbytes = nextbytes - ((n>=0) ? layers[n].maxBytes : 0);
}
// Set last layer to the overall target bitrate
n++;
layers[n].maxBytes = nextbytes;
layers[n].optimize = true;
// Re-initialize progression order changes if needed Default values
Progression[] prog1,prog2;
prog1 = (Progression[])wp.getProgressionType().getDefault();
int nValidProg = prog1.length;
for(int prg=0; prgnumLayers){
prog1[prg].lye = numLayers;
}
}
if(nValidProg==0)
throw new Error("Unable to initialize rate allocator: No "+
"default progression type has been defined.");
// Tile specific values
for(int t=0; tnumLayers){
prog1[prg].lye = numLayers;
}
}
if(nValidProg==0)
throw new Error("Unable to initialize rate allocator: No "+
"default progression type has been defined.");
}
} // End loop on tiles
if (DO_TIMING) initTime += System.currentTimeMillis()-stime;
}
/**
* This method gets all the coded code-blocks from the EBCOT entropy coder
* for every component and every tile. Each coded code-block is stored in
* a 5D array according to the component, the resolution level, the tile,
* the subband it belongs and its position in the subband.
*
* For each code-block, the valid slopes are computed and converted
* into the mantissa-exponent representation.
* */
private void getAllCodeBlocks() {
int numComps, numTiles, numBytes;
int c, r, t, s, sidx, k;
int slope;
SubbandAn subb;
CBlkRateDistStats ccb = null;
Point ncblks = null;
int last_sidx;
float fslope;
long stime = 0L;
maxSlope = 0f;
minSlope = Float.MAX_VALUE;
//Get the number of components and tiles
numComps = src.getNumComps();
numTiles = src.getNumTiles();
SubbandAn root,sb;
int cblkToEncode = 0;
int nEncCblk = 0;
ProgressWatch pw = FacilityManager.getProgressWatch();
//Get all coded code-blocks Goto first tile
src.setTile(0,0);
for (t=0; t=0; k--) {
fslope = ccb.truncSlopes[k];
if (fslope > maxSlope) maxSlope = fslope;
if (fslope < minSlope) minSlope = fslope;
sidx = getLimitedSIndexFromSlope(fslope);
for (; sidx > last_sidx; sidx--) {
RDSlopesRates[sidx] +=
ccb.truncRates[ccb.truncIdxs[k]];
}
last_sidx = getLimitedSIndexFromSlope(fslope);
}
//Fills code-blocks array
cblks[t][c][r][s][(ccb.m*ncblks.x)+ccb.n] = ccb;
ccb = null;
if(DO_TIMING) initTime += System.currentTimeMillis()-stime;
}
}
if(pw!=null) {
pw.terminateProgressWatch();
}
//Goto next tile
if(t=numLayers-1 ) {
throw new IllegalArgumentException("The first and the"+
" last layer "+
"thresholds"+
" must be optimized");
}
rdThreshold = estimateLayerThreshold(maxBytes,layers[l-1]);
}
for(int t=0; tmrlc[c]) continue;
lys[c][r] = lye;
}
} // End loop on progression
} // End loop on tiles
if (DO_TIMING) writeTime += System.currentTimeMillis()-stime;
}
/**
* Write a piece of bit stream according to the
* RES_LY_COMP_POS_PROG progression mode and between given bounds
*
* @param t Tile index.
*
* @param rs First resolution level index.
*
* @param re Last resolution level index.
*
* @param cs First component index.
*
* @param ce Last component index.
*
* @param lys First layer index for each component and resolution.
*
* @param lye Index of the last layer.
* */
public void writeResLyCompPos(int t,int rs,int re,int cs,int ce,
int lys[][],int lye) throws IOException {
boolean sopUsed; // Should SOP markers be used ?
boolean ephUsed; // Should EPH markers be used ?
int nc = src.getNumComps();
int[] mrl = new int[nc];
SubbandAn sb;
float threshold;
BitOutputBuffer hBuff = null;
byte[] bBuff = null;
int nPrec = 0;
// Max number of resolution levels in the tile
int maxResLvl = 0;
for(int c=0; cmaxResLvl) maxResLvl = mrl[c];
}
int minlys; // minimum layer start index of each component
for(int r=rs; rmaxResLvl) continue;
minlys = 100000;
for(int c=cs; c=lys[c].length) continue;
if(lmrl[c]) continue;
nPrec = numPrec[t][c][r].x*numPrec[t][c][r].y;
for(int p=0; pr; i--) {
sb = sb.subb_LL;
}
threshold = layers[l].rdThreshold;
findTruncIndices(l,c,r,t,sb,threshold,p);
hBuff = pktEnc.encodePacket(l+1,c,r,t,cblks[t][c][r],
truncIdxs[t][l][c][r],
hBuff,bBuff,p);
if(pktEnc.isPacketWritable()) {
bsWriter.writePacketHead(hBuff.getBuffer(),
hBuff.getLength(),
false,sopUsed,ephUsed);
bsWriter.writePacketBody(pktEnc.getLastBodyBuf(),
pktEnc.getLastBodyLen(),
false,pktEnc.isROIinPkt(),
pktEnc.getROILen());
}
} // End loop on precincts
} // End loop on components
} // End loop on layers
} // End loop on resolution levels
}
/**
* Write a piece of bit stream according to the
* LY_RES_COMP_POS_PROG progression mode and between given bounds
*
* @param t Tile index.
*
* @param rs First resolution level index.
*
* @param re Last resolution level index.
*
* @param cs First component index.
*
* @param ce Last component index.
*
* @param lys First layer index for each component and resolution.
*
* @param lye Index of the last layer.
* */
public void writeLyResCompPos(int t,int rs,int re,int cs,int ce,
int[][] lys,int lye) throws IOException {
boolean sopUsed; // Should SOP markers be used ?
boolean ephUsed; // Should EPH markers be used ?
int nc = src.getNumComps();
int mrl;
SubbandAn sb;
float threshold;
BitOutputBuffer hBuff = null;
byte[] bBuff = null;
int nPrec = 0;
int minlys = 100000; // minimum layer start index of each component
for(int c=cs; cmrl) continue;
if(r>=lys[c].length) continue;
if(lr; i--) {
sb = sb.subb_LL;
}
threshold = layers[l].rdThreshold;
findTruncIndices(l,c,r,t,sb,threshold,p);
hBuff = pktEnc.encodePacket(l+1,c,r,t,cblks[t][c][r],
truncIdxs[t][l][c][r],
hBuff,bBuff,p);
if(pktEnc.isPacketWritable()) {
bsWriter.writePacketHead(hBuff.getBuffer(),
hBuff.getLength(),
false,sopUsed,ephUsed);
bsWriter.writePacketBody(pktEnc.getLastBodyBuf(),
pktEnc.getLastBodyLen(),
false,pktEnc.isROIinPkt(),
pktEnc.getROILen());
}
} // end loop on precincts
} // end loop on components
} // end loop on resolution levels
} // end loop on layers
}
/**
* Write a piece of bit stream according to the
* COMP_POS_RES_LY_PROG progression mode and between given bounds
*
* @param t Tile index.
*
* @param rs First resolution level index.
*
* @param re Last resolution level index.
*
* @param cs First component index.
*
* @param ce Last component index.
*
* @param lys First layer index for each component and resolution.
*
* @param lye Index of the last layer.
* */
public void writePosCompResLy(int t,int rs,int re,int cs,int ce,
int[][] lys,int lye) throws IOException {
boolean sopUsed; // Should SOP markers be used ?
boolean ephUsed; // Should EPH markers be used ?
int nc = src.getNumComps();
int mrl;
SubbandAn sb;
float threshold;
BitOutputBuffer hBuff = null;
byte[] bBuff = null;
// Computes current tile offset in the reference grid
Point nTiles = src.getNumTiles(null);
Point tileI = src.getTile(null);
int x0siz = src.getImgULX();
int y0siz = src.getImgULY();
int xsiz = x0siz + src.getImgWidth();
int ysiz = y0siz + src.getImgHeight();
int xt0siz = src.getTilePartULX();
int yt0siz = src.getTilePartULY();
int xtsiz = src.getNomTileWidth();
int ytsiz = src.getNomTileHeight();
int tx0 = (tileI.x==0) ? x0siz : xt0siz+tileI.x*xtsiz;
int ty0 = (tileI.y==0) ? y0siz : yt0siz+tileI.y*ytsiz;
int tx1 = (tileI.x!=nTiles.x-1) ? xt0siz+(tileI.x+1)*xtsiz : xsiz;
int ty1 = (tileI.y!=nTiles.y-1) ? yt0siz+(tileI.y+1)*ytsiz : ysiz;
// Get precinct information (number,distance between two consecutive
// precincts in the reference grid) in each component and resolution
// level
PrecInfo prec; // temporary variable
int p; // Current precinct index
int gcd_x = 0; // Horiz. distance between 2 precincts in the ref. grid
int gcd_y = 0; // Vert. distance between 2 precincts in the ref. grid
int nPrec = 0; // Total number of found precincts
int[][] nextPrec = new int [ce][]; // Next precinct index in each
// component and resolution level
int minlys = 100000; // minimum layer start index of each component
int minx = tx1; // Horiz. offset of the second precinct in the
// reference grid
int miny = ty1; // Vert. offset of the second precinct in the
// reference grid.
int maxx = tx0; // Max. horiz. offset of precincts in the ref. grid
int maxy = ty0; // Max. vert. offset of precincts in the ref. grid
for(int c=cs; cmrl) continue;
if (r=0; p--) {
prec = pktEnc.getPrecInfo(t,c,r,p);
if(prec.rgulx!=tx0) {
if(prec.rgulxmaxx) maxx = prec.rgulx;
}
if(prec.rguly!=ty0){
if(prec.rgulymaxy) maxy = prec.rguly;
}
if(nPrec==0) {
gcd_x = prec.rgw;
gcd_y = prec.rgh;
} else {
gcd_x = MathUtil.gcd(gcd_x,prec.rgw);
gcd_y = MathUtil.gcd(gcd_y,prec.rgh);
}
nPrec++;
} // precincts
} // resolution levels
} // components
if(nPrec==0) {
throw new Error("Image cannot have no precinct");
}
int pyend = (maxy-miny)/gcd_y+1;
int pxend = (maxx-minx)/gcd_x+1;
int y = ty0;
int x = tx0;
for(int py=0; py<=pyend; py++) { // Vertical precincts
for(int px=0; px<=pxend; px++) { // Horiz. precincts
for(int c=cs; cmrl) continue;
if(nextPrec[c][r] >=
numPrec[t][c][r].x*numPrec[t][c][r].y) {
continue;
}
prec = pktEnc.getPrecInfo(t,c,r,nextPrec[c][r]);
if((prec.rgulx!=x) || (prec.rguly!=y)) {
continue;
}
for(int l=minlys; l=lys[c].length) continue;
if(lr; i--) {
sb = sb.subb_LL;
}
threshold = layers[l].rdThreshold;
findTruncIndices(l,c,r,t,sb,threshold,
nextPrec[c][r]);
hBuff = pktEnc.encodePacket(l+1,c,r,t,
cblks[t][c][r],
truncIdxs[t][l][c][r],
hBuff,bBuff,
nextPrec[c][r]);
if(pktEnc.isPacketWritable()) {
bsWriter.writePacketHead(hBuff.getBuffer(),
hBuff.getLength(),
false,sopUsed,
ephUsed);
bsWriter.writePacketBody(pktEnc.
getLastBodyBuf(),
pktEnc.
getLastBodyLen(),
false,
pktEnc.isROIinPkt(),
pktEnc.getROILen());
}
} // layers
nextPrec[c][r]++;
} // Resolution levels
} // Components
if(px!=pxend) {
x = minx+px*gcd_x;
} else {
x = tx0;
}
} // Horizontal precincts
if(py!=pyend) {
y = miny+py*gcd_y;
} else {
y = ty0;
}
} // Vertical precincts
// Check that all precincts have been written
for(int c=cs; cmrl) continue;
if(nextPrec[c][r]mrl) continue;
nextPrec[c] = new int[mrl+1];
if (r=0; p--) {
prec = pktEnc.getPrecInfo(t,c,r,p);
if(prec.rgulx!=tx0) {
if(prec.rgulxmaxx) maxx = prec.rgulx;
}
if(prec.rguly!=ty0){
if(prec.rgulymaxy) maxy = prec.rguly;
}
if(nPrec==0) {
gcd_x = prec.rgw;
gcd_y = prec.rgh;
} else {
gcd_x = MathUtil.gcd(gcd_x,prec.rgw);
gcd_y = MathUtil.gcd(gcd_y,prec.rgh);
}
nPrec++;
} // precincts
} // resolution levels
} // components
if(nPrec==0) {
throw new Error("Image cannot have no precinct");
}
int pyend = (maxy-miny)/gcd_y+1;
int pxend = (maxx-minx)/gcd_x+1;
int y;
int x;
for(int c=cs; cmrl) continue;
if(nextPrec[c][r] >=
numPrec[t][c][r].x*numPrec[t][c][r].y) {
continue;
}
prec = pktEnc.getPrecInfo(t,c,r,nextPrec[c][r]);
if((prec.rgulx!=x) || (prec.rguly!=y)) {
continue;
}
for(int l=minlys; l=lys[c].length) continue;
if(lr; i--) {
sb = sb.subb_LL;
}
threshold = layers[l].rdThreshold;
findTruncIndices(l,c,r,t,sb,threshold,
nextPrec[c][r]);
hBuff = pktEnc.encodePacket(l+1,c,r,t,
cblks[t][c][r],
truncIdxs[t][l][c][r],
hBuff,bBuff,
nextPrec[c][r]);
if(pktEnc.isPacketWritable()) {
bsWriter.writePacketHead(hBuff.getBuffer(),
hBuff.getLength(),
false,sopUsed,
ephUsed);
bsWriter.writePacketBody(pktEnc.
getLastBodyBuf(),
pktEnc.
getLastBodyLen(),
false,
pktEnc.isROIinPkt(),
pktEnc.getROILen());
}
} // Layers
nextPrec[c][r]++;
} // Resolution levels
if(px!=pxend) {
x = minx+px*gcd_x;
} else {
x = tx0;
}
} // Horizontal precincts
if(py!=pyend) {
y = miny+py*gcd_y;
} else {
y = ty0;
}
} // Vertical precincts
} // components
// Check that all precincts have been written
for(int c=cs; cmrl) continue;
if(nextPrec[c][r]mrl) continue;
if (r=0; p--) {
prec = pktEnc.getPrecInfo(t,c,r,p);
if(prec.rgulx!=tx0) {
if(prec.rgulxmaxx) maxx = prec.rgulx;
}
if(prec.rguly!=ty0){
if(prec.rgulymaxy) maxy = prec.rguly;
}
if(nPrec==0) {
gcd_x = prec.rgw;
gcd_y = prec.rgh;
} else {
gcd_x = MathUtil.gcd(gcd_x,prec.rgw);
gcd_y = MathUtil.gcd(gcd_y,prec.rgh);
}
nPrec++;
} // precincts
} // resolution levels
} // components
if(nPrec==0) {
throw new Error("Image cannot have no precinct");
}
int pyend = (maxy-miny)/gcd_y+1;
int pxend = (maxx-minx)/gcd_x+1;
int x,y;
for(int r=rs; rmrl) continue;
if(nextPrec[c][r] >=
numPrec[t][c][r].x*numPrec[t][c][r].y) {
continue;
}
prec = pktEnc.getPrecInfo(t,c,r,nextPrec[c][r]);
if((prec.rgulx!=x) || (prec.rguly!=y)) {
continue;
}
for(int l=minlys; l=lys[c].length) continue;
if(lr; i--) {
sb = sb.subb_LL;
}
threshold = layers[l].rdThreshold;
findTruncIndices(l,c,r,t,sb,threshold,
nextPrec[c][r]);
hBuff = pktEnc.encodePacket(l+1,c,r,t,
cblks[t][c][r],
truncIdxs[t][l][c][r],
hBuff,bBuff,
nextPrec[c][r]);
if(pktEnc.isPacketWritable()) {
bsWriter.writePacketHead(hBuff.getBuffer(),
hBuff.getLength(),
false,sopUsed,
ephUsed);
bsWriter.writePacketBody(pktEnc.
getLastBodyBuf(),
pktEnc.
getLastBodyLen(),
false,
pktEnc.isROIinPkt(),
pktEnc.getROILen());
}
} // layers
nextPrec[c][r]++;
} // Components
if(px!=pxend) {
x = minx+px*gcd_x;
} else {
x = tx0;
}
} // Horizontal precincts
if(py!=pyend) {
y = miny+py*gcd_y;
} else {
y = ty0;
}
} // Vertical precincts
} // Resolution levels
// Check that all precincts have been written
for(int c=cs; cmrl) continue;
if(nextPrec[c][r] 0; sidx--) {
if (RDSlopesRates[sidx] >= maxBytes) {
break;
}
}
// Get the corresponding minimum slope
fmint = getSlopeFromSIndex(sidx);
// Ensure that it is smaller the maximum slope
if (fmint>=fmaxt) {
sidx--;
fmint = getSlopeFromSIndex(sidx);
}
// If we are using the last entry of the summary, then that
// corresponds to all the data, Thus, set the minimum slope to 0.
if (sidx<=0) fmint = 0;
// We look for the best threshold 'ft', which is the lowest threshold
// that generates no more than 'maxBytes' code bytes.
// The search is done iteratively using a binary split algorithm. We
// start with 'fmaxt' as the maximum possible threshold, and 'fmint'
// as the minimum threshold. The threshold 'ft' is calculated as the
// middle point of 'fmaxt'-'fmint' interval. The 'fmaxt' or 'fmint'
// bounds are moved according to the number of bytes obtained from a
// simulation, where 'ft' is used as the threshold.
// We stop whenever the interval is sufficiently small, and thus
// enough precision is achieved.
// Initialize threshold as the middle point of the interval.
ft = (fmaxt+fmint)/2f;
// If 'ft' reaches 'fmint' it means that 'fmaxt' and 'fmint' are so
// close that the average is 'fmint', due to rounding. Force it to
// 'fmaxt' instead, since 'fmint' is normally an exclusive lower
// bound.
if (ft <= fmint) ft = fmaxt;
do {
// Get the number of bytes used by this layer, if 'ft' is the
// threshold, by simulation.
actualBytes = prevBytes;
src.setTile(0,0);
for (int t=0; tmaxBytes) {
// 'ft' is too low and generates too many bytes, make it the
// new minimum.
fmint = ft;
} else {
// 'ft' is too high and does not generate as many bytes as we
// are allowed too, make it the new maximum.
fmaxt = ft;
}
// Update 'ft' for the new iteration as the middle point of the
// new interval.
ft = (fmaxt+fmint)/2f;
// If 'ft' reaches 'fmint' it means that 'fmaxt' and 'fmint' are
// so close that the average is 'fmint', due to rounding. Force it
// to 'fmaxt' instead, since 'fmint' is normally an exclusive
// lower bound.
if (ft <= fmint) ft = fmaxt;
// Restore previous packet encoder state
pktEnc.restore();
// We continue to iterate, until the threshold reaches the upper
// limit of the interval, within a FLOAT_REL_PRECISION relative
// tolerance, or a FLOAT_ABS_PRECISION absolute tolerance. This is
// the sign that the interval is sufficiently small.
} while (ft < fmaxt*(1f-FLOAT_REL_PRECISION) &&
ft < (fmaxt-FLOAT_ABS_PRECISION));
// If we have a threshold which is close to 0, set it to 0 so that
// everything is taken into the layer. This is to avoid not sending
// some least significant bit-planes in the lossless case. We use the
// FLOAT_ABS_PRECISION value as a measure of "close" to 0.
if (ft <= FLOAT_ABS_PRECISION) {
ft = 0f;
} else {
// Otherwise make the threshold 'fmaxt', just to be sure that we
// will not send more bytes than allowed.
ft = fmaxt;
}
return ft;
}
/**
* This function attempts to estimate a rate-distortion slope threshold
* which will achieve a target number of code bytes close the
* `targetBytes' value.
*
* @param targetBytes The target number of bytes for the current layer
*
* @param lastLayer The previous layer information.
*
* @return The value of the slope threshold for the estimated layer
* */
private float estimateLayerThreshold(int targetBytes,
EBCOTLayer lastLayer) {
float log_sl1; // The log of the first slope used for interpolation
float log_sl2; // The log of the second slope used for interpolation
float log_len1; // The log of the first length used for interpolation
float log_len2; // The log of the second length used for interpolation
float log_isl; // The log of the interpolated slope
float log_ilen; // Log of the interpolated length
float log_ab; // Log of actual bytes in last layer
int sidx; // Index into the summary R-D info array
float log_off; // The log of the offset proportion
int tlen; // The corrected target layer length
float lthresh; // The threshold of the last layer
float eth; // The estimated threshold
// In order to estimate the threshold we base ourselves in the summary
// R-D info in RDSlopesRates. In order to use it we must compensate
// for the overhead of the packet heads. The proportion of overhead is
// estimated using the last layer simulation results.
// NOTE: the model used in this method is that the slope varies
// linearly with the log of the rate (i.e. length).
// NOTE: the model used in this method is that the distortion is
// proprotional to a power of the rate. Thus, the slope is also
// proportional to another power of the rate. This translates as the
// log of the slope varies linearly with the log of the rate, which is
// what we use.
// 1) Find the offset of the length predicted from the summary R-D
// information, to the actual length by using the last layer.
// We ensure that the threshold we use for estimation actually
// includes some data.
lthresh = lastLayer.rdThreshold;
if (lthresh > maxSlope) lthresh = maxSlope;
// If the slope of the last layer is too small then we just include
// all the rest (not possible to do better).
if (lthresh < FLOAT_ABS_PRECISION) return 0f;
sidx = getLimitedSIndexFromSlope(lthresh);
// If the index is outside of the summary info array use the last two,
// or first two, indexes, as appropriate
if (sidx >= RD_SUMMARY_SIZE-1) sidx = RD_SUMMARY_SIZE-2;
// Get the logs of the lengths and the slopes
if (RDSlopesRates[sidx+1] == 0) {
// Pathological case, we can not use log of 0. Add
// RDSlopesRates[sidx]+1 bytes to the rates (just a crude simple
// solution to this rare case)
log_len1 = (float)Math.log((RDSlopesRates[sidx]<<1)+1);
log_len2 = (float)Math.log(RDSlopesRates[sidx]+1);
log_ab = (float)Math.log(lastLayer.actualBytes+
RDSlopesRates[sidx]+1);
}
else {
log_len1 = (float)Math.log(RDSlopesRates[sidx]);
log_len2 = (float)Math.log(RDSlopesRates[sidx+1]);
log_ab = (float)Math.log(lastLayer.actualBytes);
}
log_sl1 = (float)Math.log(getSlopeFromSIndex(sidx));
log_sl2 = (float)Math.log(getSlopeFromSIndex(sidx+1));
log_isl = (float)Math.log(lthresh);
log_ilen = log_len1 +
(log_isl-log_sl1)*(log_len1-log_len2)/(log_sl1-log_sl2);
log_off = log_ab - log_ilen;
// Do not use negative offsets (i.e. offset proportion larger than 1)
// since that is probably a sign that our model is off. To be
// conservative use an offset of 0 (i.e. offset proportiojn 1).
if (log_off < 0) log_off = 0f;
// 2) Correct the target layer length by the offset.
tlen = (int)(targetBytes/(float)Math.exp(log_off));
// 3) Find, from the summary R-D info, the thresholds that generate
// lengths just above and below our corrected target layer length.
// Look for the index in the summary info array that gives the largest
// length smaller than the target length
for (sidx = RD_SUMMARY_SIZE-1; sidx>=0 ; sidx--) {
if (RDSlopesRates[sidx] >= tlen) break;
}
sidx++;
// Correct if out of the array
if (sidx >= RD_SUMMARY_SIZE) sidx = RD_SUMMARY_SIZE-1;
if (sidx <= 0) sidx = 1;
// Get the log of the lengths and the slopes that are just above and
// below the target length.
if (RDSlopesRates[sidx] == 0) {
// Pathological case, we can not use log of 0. Add
// RDSlopesRates[sidx-1]+1 bytes to the rates (just a crude simple
// solution to this rare case)
log_len1 = (float)Math.log(RDSlopesRates[sidx-1]+1);
log_len2 = (float)Math.log((RDSlopesRates[sidx-1]<<1)+1);
log_ilen = (float)Math.log(tlen+RDSlopesRates[sidx-1]+1);
}
else {
// Normal case, we can safely take the logs.
log_len1 = (float)Math.log(RDSlopesRates[sidx]);
log_len2 = (float)Math.log(RDSlopesRates[sidx-1]);
log_ilen = (float)Math.log(tlen);
}
log_sl1 = (float)Math.log(getSlopeFromSIndex(sidx));
log_sl2 = (float)Math.log(getSlopeFromSIndex(sidx-1));
// 4) Interpolate the two thresholds to find the target threshold.
log_isl = log_sl1 +
(log_ilen-log_len1)*(log_sl1-log_sl2)/(log_len1-log_len2);
eth = (float)Math.exp(log_isl);
// Correct out of bounds results
if (eth > lthresh) eth = lthresh;
if (eth < FLOAT_ABS_PRECISION) eth = 0f;
// Return the estimated threshold
return eth;
}
/**
* This function finds the new truncation points indices for a packet. It
* does so by including the data from the code-blocks in the component,
* resolution level and tile, associated with a R-D slope which is larger
* than or equal to 'fthresh'.
*
* @param layerIdx The index of the current layer
*
* @param compIdx The index of the current component
*
* @param lvlIdx The index of the current resolution level
*
* @param tileIdx The index of the current tile
*
* @param subb The LL subband in the resolution level lvlIdx, which is
* parent of all the subbands in the packet. Except for resolution level 0
* this subband is always a node.
*
* @param fthresh The value of the rate-distortion threshold
* */
private void findTruncIndices(int layerIdx, int compIdx, int lvlIdx,
int tileIdx, SubbandAn subb,
float fthresh, int precinctIdx) {
int minsbi, maxsbi, b, bIdx, n;
Point ncblks = null;
SubbandAn sb;
CBlkRateDistStats cur_cblk;
PrecInfo prec = pktEnc.getPrecInfo(tileIdx,compIdx,lvlIdx,precinctIdx);
Point cbCoord;
sb = subb;
while(sb.subb_HH != null) {
sb = sb.subb_HH;
}
minsbi = (lvlIdx==0) ? 0 : 1;
maxsbi = (lvlIdx==0) ? 1 : 4;
int yend,xend;
sb = (SubbandAn)subb.getSubbandByIdx(lvlIdx, minsbi);
for(int s=minsbi; sIf the value to return is lower than 0, 0 is returned. If it is
* larger than the maximum table index, then the maximum is returned.
*
* @param slope The slope value
*
* @return The index for the summary table of the slope.
* */
private static int getLimitedSIndexFromSlope(float slope) {
int idx;
idx = (int)Math.floor(Math.log(slope)/LOG2) + RD_SUMMARY_OFF;
if (idx < 0) {
return 0;
}
else if (idx >= RD_SUMMARY_SIZE) {
return RD_SUMMARY_SIZE-1;
}
else {
return idx;
}
}
/**
* Returns the minimum slope value associated with a summary table
* index. This minimum slope is just 2^(index-RD_SUMMARY_OFF).
*
* @param index The summary index value.
*
* @return The minimum slope value associated with a summary table index.
* */
private static float getSlopeFromSIndex(int index) {
return (float)Math.pow(2,(index-RD_SUMMARY_OFF));
}
}