ebcotrateallocator.java

jpeg2000算法实现

	int precinctIdxA[][][];
	Coord xys[],xyInc;
        SubbandAn sb;
        float threshold;
        BitOutputBuffer hBuff = null;
        byte[] bBuff = null;
	int numLvls;

	// Max number of decomposition levels in each tile
	numLvls = mrl[t][0];
	for(int c=0; c < numComps; c++)
	    if(mrl[t][c]>numLvls)
		numLvls = mrl[t][c];
	numLvls++;
	
	precinctIdxA = new int[numComps][numLvls][numLayers];
	
	// Coord[] xys = packetEnc.getSotEotArrayMax(t, c);
	xys = packetEnc.getSotEotArrayMax(t, 0);
	x0 = xys[0].x;
	y0 = xys[0].y;
	x1 = xys[1].x;
	y1 = xys[1].y;
	
	// Coord xyInc = packetEnc.getIncArrayMax(t, c);
	xyInc = packetEnc.getIncArrayMax(t, 0);
	x_inc = xyInc.x;
	y_inc = xyInc.y;
	
	for(int r=rs ; r<re ; r++){// Loop on resolution levels
	    
	    // Loop on packets
	    for(int yr=y0; yr<y1; yr+=y_inc ){
		for(int xr=x0; xr<x1; xr+=x_inc ){
		    
		    for(int c=cs ; c<ce ; c++){ // Loop on each component
			
			// If no more decomposition levels for this
			// component
			if(r>mrl[t][c]) continue;
			
			// set boolean sopUsed here (SOP markers)
			sopUsed = ((String)encSpec.sops.
				   getTileDef(t)).equalsIgnoreCase("on");
			// set boolean ephUsed here (EPH markers)
			ephUsed = ((String)encSpec.
				   ephs.getTileDef(t)).equalsIgnoreCase("on");
			
			sb = src.getSubbandTree(t,c);
			numLvls = sb.resLvl+1;
			for(int i=numLvls-1; i>r; i--){
			    sb = sb.subb_LL;
			}
			
			xyInc = packetEnc.getIncArray(t,c,r);
			x_inc_rl = xyInc.x;
			y_inc_rl = xyInc.y;
			xys = packetEnc.getSotEotArray(t,c,r);
			x0_rl = xys[0].x;
			y0_rl = xys[0].y;
			
			if( ( (xr==x0) || (xr%x_inc_rl==0) ) &&
			    ( (yr==y0) || (yr%y_inc_rl==0) ) ){
			    
			    for(int l=lys[c][r] ;l<lye ;l++){ //loop on layers
				
                                if(precinctIdxA[c][r][l]>=
                                   maxNumPrec[t][c][r].x*maxNumPrec[t][c][r].y)
                                    continue;

				precinctIdx = precinctIdxA[c][r][l];
				
				threshold = layers[l].rdThreshold;

				findTruncIndices(l,c,r,t,sb,threshold,
						 precinctIdx);
				
				hBuff = packetEnc.
				    encodePacket(l+1,c,r,t,cblks[t][c][r],
						 truncIdxs[t][l][c][r],
						 hBuff, bBuff, precinctIdx);
				
				if( packetEnc.isPacketWritable() ){
				    bsWriter.
					writePacketHead(hBuff.getBuffer(),
							hBuff.getLength(), 
							false,sopUsed,ephUsed);
				    bsWriter.
					writePacketBody(packetEnc.
							getLastBodyBuf(),
							packetEnc.
							getLastBodyLen(),
                                                        false,packetEnc.
                                                        isROIinPkt(),packetEnc.
                                                        getROILen());
				}
				precinctIdxA[c][r][l]++;
			    } // End loop on layers
			} // test packet
		    } // end loop on components
		}
	    } // end loop on precincts
	} // end loop on resolution levels
    }

    /**
     * This function implements the rate-distortion optimization algorithm.
     * It saves the state of any previously generated bit-stream layers and
     * then simulate the formation of a new layer in the bit stream as often
     * as necessary to find the smallest rate-distortion threshold such that
     * the total number of bytes required to represent the layer does not
     * exceed `maxBytes' minus `prevBytes'.  It then restores the state of any
     * previously generated bit-stream layers and returns the threshold.
     *
     * @param layerIdx The index of the current layer
     *
     * @param fmaxt The maximum admissible slope value. Normally the threshold
     * slope of the previous layer.
     *
     * @param maxBytes The maximum number of bytes that can be written. It
     * includes the length of the current layer bistream length and all the
     * previous layers bit streams.
     *
     * @param prevBytes The number of bytes of all the previous layers.
     *
     * @return The value of the slope threshold.
     * */
    private float optimizeBitstreamLayer (int layerIdx, float fmaxt,
                                          int maxBytes, int prevBytes) 
        throws IOException {

        int numTiles;         // The total number of tiles
        int numComp;          // The total number of components
        int numLvls;          // The total number of resolution levels
        int c, n, t;          // Component, rsolution and tile indexes
        int actualBytes;      // Actual number of bytes for a layer
        float fmint;          // Minimum of the current threshold interval
        float ft;             // Current threshold
        SubbandAn sb;         // Current subband
        BitOutputBuffer hBuff;// The packet head buffer
        byte[] bBuff;         // The packet body buffer
        int sidx;             // The index in the summary table
        boolean sopUsed;      // Should SOP markers be used ?
        boolean ephUsed;      // Should EPH markers be used ?
        int precinctIdx;      // Precinct index for current packet
        int x0, y0, x1, y1;   // Coordinates of tile-component
        int x_inc, y_inc;     // Steps for precinct loop 
        int x_inc_rl, y_inc_rl;   // Steps for precinct loop
        int x0_rl, y0_rl;         // Tile-component coordinates in res level
        Coord xys[], xyInc;       // Start coordinates and increment

        // Save the packet encoder state
        packetEnc.save();

        numTiles = src.getNumTiles();
        numComp = src.getNumComps();
        hBuff = null;
        bBuff = null;
        
        // Estimate the minimum slope to start with from the summary
        // information in 'RDSlopesRates'. This is a real minimum since it
        // does not include the packet head overhead, which is always
        // non-zero.

        // Look for the summary entry that gives 'maxBytes' or more data
        for (sidx=RD_SUMMARY_SIZE-1; sidx > 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 (t=0; t < numTiles; t++){
		for (c=0; c < numComp; c++) {
		    // set boolean sopUsed here (SOP markers)
		    sopUsed = ((String)encSpec.sops.getTileDef(t)).
			equalsIgnoreCase("on");
		    // set boolean ephUsed here (EPH markers)
		    ephUsed = ((String)encSpec.ephs.getTileDef(t)).
			equalsIgnoreCase("on");
		    
		    // Get LL subband
                    sb = (SubbandAn) src.getSubbandTree(t,c);
                    numLvls = sb.resLvl + 1;
		    sb = (SubbandAn) sb.getSubbandByIdx(0,0);
		    //loop on resolution levels
		    for(int r=0; r<numLvls; r++) { 
                        precinctIdx = 0;

			// Start-end of tile-component
			xys = packetEnc.getSotEotArrayMax(t,c);
			x0 = xys[0].x;
			y0 = xys[0].y;
			x1 = xys[1].x;
			y1 = xys[1].y;
			
			xyInc = packetEnc.getIncArrayMax(t,c);
			x_inc = xyInc.x;
			y_inc = xyInc.y;
			
			xyInc = packetEnc.getIncArray(t,c,r);
			x_inc_rl = xyInc.x;
			y_inc_rl = xyInc.y;
			
			xys = packetEnc.getSotEotArray(t,c,r);
			x0_rl = xys[0].x;
			y0_rl = xys[0].y;

			for(int yr=y0 ; yr<y1 ; yr+=y_inc )
			    for(int xr=x0 ; xr<x1 ; xr+=x_inc )
				if ( ( (xr==x0) || (xr%x_inc_rl==0) ) &&
				     ( (yr==y0) || (yr%y_inc_rl==0) ) ){

                                    if(precinctIdx>=maxNumPrec[t][c][r].x*
                                       maxNumPrec[t][c][r].y)
                                        continue;

                                    findTruncIndices(layerIdx,c,r,t,sb,ft,
                                                     precinctIdx);
                                    hBuff = packetEnc.
                                        encodePacket(layerIdx+1,c,r,t,
                                                     cblks[t][c][r], 
                                                     truncIdxs[t][layerIdx]
                                                     [c][r],
                                                     hBuff,bBuff,
                                                     precinctIdx);
                                    bBuff = packetEnc.getLastBodyBuf();
                                    actualBytes += bsWriter.
                                        writePacketHead(hBuff.getBuffer(),
                                                        hBuff.getLength(), 
                                                        true, sopUsed,ephUsed);
                                    actualBytes += bsWriter.
                                        writePacketBody(bBuff,
                                                        packetEnc.
                                                        getLastBodyLen(), 
                                                        true,packetEnc.
                                                        isROIinPkt(),packetEnc.
                                                        getROILen());
                                    precinctIdx++;
                                } // end loop on precincts
			sb = sb.parent;
		    }
		}
	    } // End loop on tiles

            // Move the interval bounds according to simulation result
            if (actualBytes > maxBytes) {
                // '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
            packetEnc.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
        // proporti