C++ CByteImage类代码示例

void computeDSI(CByteImage im1,       // source (reference) image
		CByteImage im2,       // destination (match) image
		MRF::CostVal *&dsi,   // computed cost volume
		int nD,               // number of disparities
		int birchfield,       // use Birchfield/Tomasi costs
		int squaredDiffs,     // use squared differences
		int truncDiffs)       // truncated differences
{
    CShape sh = im1.Shape();
    int width = sh.width, height = sh.height, nB = sh.nBands;
    dsi = new MRF::CostVal[width * height * nD];

    int nColors = __min(3, nB);

    // worst value for sumdiff below 
    int worst_match = nColors * (squaredDiffs ? 255 * 255 : 255);
    // truncation threshold - NOTE: if squared, don't multiply by nColors (Eucl. dist.)
    int maxsumdiff = squaredDiffs ? truncDiffs * truncDiffs : nColors * abs(truncDiffs);
    // value for out-of-bounds matches
    int badcost = __min(worst_match, maxsumdiff);

    int dsiIndex = 0;
    for (int y = 0; y < height; y++) {
	for (int x = 0; x < width; x++) {
	    uchar *pix1 = &im1.Pixel(x, y, 0);
	    for (int d = 0; d < nD; d++) {
		int x2 = x-d;
		int dsiValue;
                if (x2 >= 0 && d < nD) { // in bounds
		    uchar *pix2 = &im2.Pixel(x2, y, 0);
                    int sumdiff = 0;
                    for (int b = 0; b < nColors; b++) {
			int diff = 0;
			if (birchfield) {
			    // Birchfield/Tomasi cost
			    int im1c = pix1[b];
			    int im1l = x == 0?   im1c : (im1c + pix1[b - nB]) / 2;
			    int im1r = x == width-1? im1c : (im1c + pix1[b + nB]) / 2;
			    int im2c = pix2[b];
			    int im2l = x2 == 0?   im2c : (im2c + pix2[b - nB]) / 2;
			    int im2r = x2 == width-1? im2c : (im2c + pix2[b + nB]) / 2;
			    int min1 = __min(im1c, __min(im1l, im1r));
			    int max1 = __max(im1c, __max(im1l, im1r));
			    int min2 = __min(im2c, __min(im2l, im2r));
			    int max2 = __max(im2c, __max(im2l, im2r));
			    int di1 = __max(0, __max(im1c - max2, min2 - im1c));
			    int di2 = __max(0, __max(im2c - max1, min1 - im2c));
			    diff = __min(di1, di2);
			} else {
			    // simple absolute difference
			    int di = pix1[b] - pix2[b];
			    diff = abs(di);
			}
			// square diffs if requested (Birchfield too...)
			sumdiff += (squaredDiffs ? diff * diff : diff);
                    }
		    // truncate diffs
		    dsiValue = __min(sumdiff, maxsumdiff);
                } else { // out of bounds: use maximum truncated cost
		    dsiValue = badcost;
		}
		//int x0=-140, y0=-150;
		//if (x==x0 && y==y0)
		//    printf("dsi(%d,%d,%2d)=%3d\n", x, y, d, dsiValue); 

		// The cost of pixel p and label l is stored at dsi[p*nLabels+l]
		dsi[dsiIndex++] = dsiValue;
	    }
	}
    }
    //exit(1);
}

void evaldisp(CFloatImage disp, CFloatImage gtdisp, CByteImage mask, float badthresh, int maxdisp, int rounddisp)
{
    CShape sh = gtdisp.Shape();
    CShape sh2 = disp.Shape();
    CShape msh = mask.Shape();
    int width = sh.width, height = sh.height;
    int width2 = sh2.width, height2 = sh2.height;
    int scale = width / width2;

    if ((!(scale == 1 || scale == 2 || scale == 4))
	|| (scale * width2 != width)
	|| (scale * height2 != height)) {
	printf("   disp size = %4d x %4d\n", width2, height2);
	printf("GT disp size = %4d x %4d\n", width,  height);
	throw CError("GT disp size must be exactly 1, 2, or 4 * disp size");
    }

    int usemask = (msh.width > 0 && msh.height > 0);
    if (usemask && (msh != sh))
	throw CError("mask image must have same size as GT\n");

    int n = 0;
    int bad = 0;
    int invalid = 0;
    float serr = 0;
    for (int y = 0; y < height; y++) {
	for (int x = 0; x < width; x++) {
	    float gt = gtdisp.Pixel(x, y, 0);
	    if (gt == INFINITY) // unknown
		continue;
	    float d = scale * disp.Pixel(x / scale, y / scale, 0);
		float maxd = scale * maxdisp; // max disp range
	    int valid = (d != INFINITY && d <= maxd * 1000);
	    if (valid) {
		d = __max(0, __min(maxd, d)); // clip disps to max disp range
	    }
	    if (valid && rounddisp)
		d = round(d);
	    float err = fabs(d - gt);
	    if (usemask && mask.Pixel(x, y, 0) != 255) { // don't evaluate pixel
	    } else {
		n++;
		if (valid) {
		    serr += err;
		    if (err > badthresh) {
			bad++;
		    }
		} else {// invalid (i.e. hole in sparse disp map)
		    invalid++;
		}
	    }
	}
    }
    float badpercent =  100.0*bad/n;
    float invalidpercent =  100.0*invalid/n;
    float totalbadpercent =  100.0*(bad+invalid)/n;
    float avgErr = serr / (n - invalid); // CHANGED 10/14/2014 -- was: serr / n
    //printf("mask  bad%.1f  invalid  totbad   avgErr\n", badthresh);
    printf("%4.1f  %6.2f  %6.2f   %6.2f  %6.2f\n",   100.0*n/(width * height),
	   badpercent, invalidpercent, totalbadpercent, avgErr);
}

Feature 
HOGFeatureExtractor::operator()(const CByteImage& img_) const
{
	/******** BEGIN TODO ********/
	// Compute the Histogram of Oriented Gradients feature
	// Steps are:
	// 1) Compute gradients in x and y directions. We provide the 
	//    derivative kernel proposed in the paper in _kernelDx and
	//    _kernelDy.
	// 2) Compute gradient magnitude and orientation
	// 3) Add contribution each pixel to HOG cells whose
	//    support overlaps with pixel. Each cell has a support of size
	//    _cellSize and each histogram has _nAngularBins.
	// 4) Normalize HOG for each cell. One simple strategy that is
	//    is also used in the SIFT descriptor is to first threshold
	//    the bin values so that no bin value is larger than some
	//    threshold (we leave it up to you do find this value) and
	//    then re-normalize the histogram so that it has norm 1. A more 
	//    elaborate normalization scheme is proposed in Dalal & Triggs
	//    paper but we leave that as extra credit.
	// 
	// Useful functions:
	// convertRGB2GrayImage, TypeConvert, WarpGlobal, Convolve

	
	int xCells = ceil(1.*img_.Shape().width  / _cellSize);
	int yCells = ceil(1.*img_.Shape().height / _cellSize);

	CFloatImage HOGHist(xCells, yCells, _nAngularBins);
	HOGHist.ClearPixels();

	CByteImage gray(img_.Shape());
	CFloatImage grayF(img_.Shape().width, img_.Shape().height, 1);

	convertRGB2GrayImage(img_, gray);
	TypeConvert(gray, grayF);

	CFloatImage diffX( img_.Shape()), diffY( img_.Shape());

	Convolve(grayF, diffX, _kernelDx);
	Convolve(grayF, diffY, _kernelDy);
	
	CFloatImage grad(grayF.Shape()), grad2(grayF.Shape());
	CFloatImage angl(grayF.Shape()), angl2(grayF.Shape());
	
	for (int y = 0; y <grayF.Shape().height; y++){
		for (int x = 0; x<grayF.Shape().width; x++) {
			grad2.Pixel(x,y,0) = (diffX.Pixel(x,y,0) * diffX.Pixel(x,y,0) +
							     diffY.Pixel(x,y,0) * diffY.Pixel(x,y,0));
			angl2.Pixel(x,y,0) = atan(diffY.Pixel(x,y,0) / abs(diffY.Pixel(x,y,0)));
		}
	}

	// Bilinear Filter
	ConvolveSeparable(grad2, grad, ConvolveKernel_121,ConvolveKernel_121,1);
	ConvolveSeparable(angl2, angl, ConvolveKernel_121,ConvolveKernel_121,1);
	//WriteFile(diffX, "angle.tga");
	//WriteFile(diffY, "angleG.tga");

	for (int y = 0; y <grayF.Shape().height; y++){
		for (int x = 0; x<grayF.Shape().width; x++) {
			// Fit in the bins
			int a = angl.Pixel(x,y,0) / 3.14 * (_nAngularBins) + _nAngularBins/2;		
			// Histogram
			HOGHist.Pixel(floor(1.*x / _cellSize),
						  floor(1.*y / _cellSize),
						  a) += grad.Pixel(x,y,0);
		}
	}
	
	// Normalization 
	float threshold = 0.7;
	for (int y = 0; y < yCells; y++){
		for (int x = 0; x < xCells; x++){
			float total = 0;
			for (int a = 0; a < _nAngularBins; a++) {
				if (HOGHist.Pixel(x,y,a) > threshold)
					HOGHist.Pixel(x,y,a) = threshold;
				// Sum for normalization
				total += HOGHist.Pixel(x,y,a);
			}
			for (int a = 0;a< _nAngularBins; a++) {
				HOGHist.Pixel(x,y,a) /= total;
				
			}
		}
	}
	
	return HOGHist;
	/******** END TODO ********/

}

/******************* TO DO *********************
* AccumulateBlend:
*	INPUT:
*		img: a new image to be added to acc
*		acc: portion of the accumulated image where img is to be added
*       M: the transformation mapping the input image 'img' into the output panorama 'acc'
*		blendWidth: width of the blending function (horizontal hat function;
*	    try other blending functions for extra credit)
*	OUTPUT:
*		add a weighted copy of img to the subimage specified in acc
*		the first 3 band of acc records the weighted sum of pixel colors
*		the fourth band of acc records the sum of weight
*/
static void AccumulateBlend(CByteImage& img, CFloatImage& acc, CTransform3x3 M, float blendWidth)
{
    // BEGIN TODO
    // Fill in this routine
	// get shape of acc and img
	CShape sh = img.Shape();
    int width = sh.width;
    int height = sh.height;
	CShape shacc = acc.Shape();
    int widthacc = shacc.width;
    int heightacc = shacc.height;
	
	// get the bounding box of img in acc
	int min_x, min_y, max_x, max_y;
	ImageBoundingBox(img, M, min_x, min_y, max_x, max_y);

	CVector3 p;
	double newx, newy;

	// Exposure Compensation
	double lumaScale = 1.0;
	double lumaAcc = 0.0;
	double lumaImg = 0.0;
	int cnt = 0;

	for (int ii = min_x; ii < max_x; ii++)
		for (int jj = min_y; jj < max_y; jj++)
		{
			// flag: current pixel black or not
			bool flag = false;
			p[0] = ii; p[1] = jj; p[2] = 1;
			p = M.Inverse() * p;
			newx = p[0] / p[2];
			newy = p[1] / p[2];
			// If in the overlapping region
			if (newx >=0 && newx < width && newy >=0 && newy < height)
			{
				if (acc.Pixel(ii,jj,0) == 0 &&
					acc.Pixel(ii,jj,1) == 0 &&
					acc.Pixel(ii,jj,2) == 0)
					flag = true;
				if (img.PixelLerp(newx,newy,0) == 0 &&
					img.PixelLerp(newx,newy,1) == 0 &&
					img.PixelLerp(newx,newy,2) == 0)
					flag = true;
				if (!flag)
				{
					// Compute Y using RGB (RGB -> YUV)
					lumaAcc = 0.299 * acc.Pixel(ii,jj,0) +
							   0.587 * acc.Pixel(ii,jj,1) +
							   0.114 * acc.Pixel(ii,jj,2);
					lumaImg = 0.299 * img.PixelLerp(newx,newy,0) +
							   0.587 * img.PixelLerp(newx,newy,1) +
							   0.114 * img.PixelLerp(newx,newy,2);
					
					if (lumaImg != 0)
					{
						double scale = lumaAcc / lumaImg;
						if (scale > 0.5 && scale < 2)
						{
							lumaScale += lumaAcc / lumaImg;
							cnt++;
						}
					}
				}
			}
		}

	if (cnt != 0)
		lumaScale = lumaScale / (double)cnt;
	else lumaScale = 1.0;

	// add every pixel in img to acc, feather the region withing blendwidth to the bounding box,
	// pure black pixels (caused by warping) are not added
	double weight;
	
	for (int ii = min_x; ii < max_x; ii++)
		for (int jj = min_y; jj < max_y; jj++)
		{
			p[0] = ii; p[1] = jj; p[2] = 1;
			p = M.Inverse() * p;
			newx = p[0] / p[2];
			newy = p[1] / p[2];
			if ((newx >= 0) && (newx < width-1) && (newy >= 0) && (newy < height-1))
			{
				weight = 1.0;
				if ( (ii >= min_x) && (ii < (min_x+blendWidth)) )
//.........这里部分代码省略.........

void ReadFileTGA(CByteImage& img, const char* filename)
{
    // Open the file and read the header
    FILE *stream = fopen(filename, "rb");
    if (stream == 0)
        throw CError("ReadFileTGA: could not open %s", filename);
    CTargaHead h;
    if (fread(&h, sizeof(CTargaHead), 1, stream) != 1)
	    throw CError("ReadFileTGA(%s): file is too short", filename);

    // Throw away the image descriptor
    if (h.idLength > 0)
    {
        char* tmp = new char[h.idLength];
        int nread = fread(tmp, sizeof(uchar), h.idLength, stream);
        delete tmp;   // throw away this data
        if (nread != h.idLength)
	        throw CError("ReadFileTGA(%s): file is too short", filename);
    }
    bool isRun = (h.imageType & 8) != 0;
    bool reverseRows = (h.descriptor & TargaScreenOrigin) != 0;
    int fileBytes = (h.pixelSize + 7) / 8;

    // Read the colormap
    uchar colormap[TargaCMapSize][TargaCMapBands];
    int cMapSize = 0;
    bool grayRamp = false;
    if (h.colorMapType == 1)
    {
        cMapSize = (h.cMapLength[1] << 8) + h.cMapLength[0];
        if (h.cMapBits != 24)
            throw CError("ReadFileTGA(%s): only 24-bit colormap currently supported", filename);
	    int l = fileBytes * cMapSize;
        if (l > TargaCMapSize * TargaCMapBands)
	        throw CError("ReadFileTGA(%s): colormap is too large", filename);
	    if (fread(colormap, sizeof(uchar), l, stream) != l)
	        throw CError("ReadFileTGA(%s): could not read the colormap", filename);

        // Check if it's just a standard gray ramp
	int i;
        for (i = 0; i < cMapSize; i++) {
            for (int j = 0; j < TargaCMapBands; j++)
                if (colormap[i][j] != i)
                    break;
        }
        grayRamp = (i == cMapSize);    // didn't break out too soon
    }
    bool isGray = 
        h.imageType == TargaRawBW || h.imageType == TargaRunBW ||
        grayRamp &&
        (h.imageType == TargaRawColormap || h.imageType == TargaRunColormap);
    bool isRaw = h.imageType == TargaRawBW || h.imageType == TargaRawRGB ||
        h.imageType == TargaRawRGB && isGray;

    // Determine the image shape
    CShape sh(h.width, h.height, (isGray) ? 1 : 4);
    
    // Allocate the image if necessary
    img.ReAllocate(sh, false);

    // Construct a run-length code reader
    CTargaRLC rlc(! isRaw);

    // Read in the rows
    for (int y = 0; y < sh.height; y++)
    {
        int yr = reverseRows ? sh.height-1-y : y;
        uchar* ptr = (uchar *) img.PixelAddress(0, yr, 0);
        if (fileBytes == sh.nBands && isRaw)
        {
            // Special case for raw image, same as destination
            int n = sh.width*sh.nBands;
    	    if (fread(ptr, sizeof(uchar), n, stream) != n)
    	        throw CError("ReadFileTGA(%s): file is too short", filename);
        }
        else
        {
            // Read one pixel at a time
            for (int x = 0; x < sh.width; x++, ptr += sh.nBands)
            {
                uchar* buf = rlc.getBytes(fileBytes, stream);
                if (fileBytes == 1 && sh.nBands == 1)
                {
                    ptr[0] = buf[0];
                }
                else if (fileBytes == 1 && sh.nBands == 4)
                {
                    for (int i = 0; i < 3; i++)
                        ptr[i] = (isGray) ? buf[0] : colormap[buf[0]][i];
                    ptr[3] = 255;   // full alpha;
                }
                else if ((fileBytes == 3 || fileBytes == 4) && sh.nBands == 4)
                {
                    int i;
                    for (i = 0; i < fileBytes; i++)
                        ptr[i] = buf[i];
                    if (i == 3) // missing alpha channel
                        ptr[3] = 255;   // full alpha;
                }
                else
//.........这里部分代码省略.........

// TODO: the following function should go somewhere else, perhaps in ImageIO.cpp
// Make sure the image has the smallest number of bands before writing.
// That is, if it's 4 bands with full alpha, reduce to 3 bands.  
// If it's 3 bands with constant colors, make it 1-band.
CByteImage removeRedundantBands(CByteImage img)
{
    CShape sh = img.Shape();
	int w = sh.width, h = sh.height, nB = sh.nBands;
	int x, y;
	if (nB < 3)
		return img;

	// check if full alpha if alpha channel present
	bool fullAlpha = true;
	if (nB == 4) {
		for (y = 0; y < h && fullAlpha; y++) {
			uchar *pix = &img.Pixel(0, y, 0);
			for (x = 0; x < w; x++) {
				if (pix[3] != 255) {
					fullAlpha = false;
					break;
				}
				pix += nB;
			}
		}
	}
	if (!fullAlpha)
		return img;

	// check for equal colors
	bool equalColors = true;
	for (y = 0; y < h && equalColors; y++) {
		uchar *pix = &img.Pixel(0, y, 0);
		for (x = 0; x < w; x++) {
			if (pix[0] != pix[1] ||
				pix[0] != pix[2] ||
				pix[1] != pix[2]) {
					equalColors = false;
					break;
				}
				pix += nB;
		}
	}
	// at this point, if nB == 4 we can reduce to at least 3 bands,
	// and if equalColors we can reduce to 1 band.
	if (! equalColors && nB < 4)
		return img;

	int newNB = equalColors ? 1 : 3;

	if (DEBUG_ImageIOpng)
		fprintf(stderr, "reducing from %d to %d bands\n", nB, newNB);

	CShape sh2(w, h, newNB);
	CByteImage img2(sh2);
	
	for (y = 0; y < h; y++) {
		uchar *pix = &img.Pixel(0, y, 0);
		uchar *pix2 = &img2.Pixel(0, y, 0);
		for (x = 0; x < w; x++) {
			for (int b = 0; b < newNB; b++) {
				pix2[b] = pix[b];
			}
			pix += nB;
			pix2 += newNB;
		}
	}

	return img2;
}

void WriteFilePNG(CByteImage img, const char* filename)
{
	img = removeRedundantBands(img);

    CShape sh = img.Shape();
    int width = sh.width, height = sh.height, nBands = sh.nBands;

	// Make sure the image has the smallest number of bands before writing.
	// That is, if it's 4 bands with full alpha, reduce to 3 bands.  
	// If it's 3 bands with constant colors, make it 1-band.

    FILE *stream = fopen(filename, "wb");
    if (stream == 0)
        throw CError("WriteFilePNG: could not open %s", filename);

    png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL,
      (png_error_ptr)pngfile_error, (png_error_ptr)NULL);

	if (!png_ptr) {
        fclose(stream);
		throw CError("WriteFilePNG: error creating png structure");
	}

    info_ptr = png_create_info_struct(png_ptr);
    if (!info_ptr) {
        fclose(stream);
		throw CError("WriteFilePNG: error creating png structure");
    }

	png_init_io(png_ptr, stream);

	int bits = 8;
	int colortype =
		nBands == 1 ? PNG_COLOR_TYPE_GRAY :
		nBands == 3 ? PNG_COLOR_TYPE_RGB :
	                  PNG_COLOR_TYPE_RGB_ALPHA;
	png_set_IHDR(png_ptr, info_ptr, width, height, 
		bits, colortype,
		PNG_INTERLACE_NONE, 
		PNG_COMPRESSION_TYPE_DEFAULT, 
		PNG_FILTER_TYPE_DEFAULT);

	// write the file header information
	png_write_info(png_ptr, info_ptr);

	// swap the BGR pixels in the DiData structure to RGB
	png_set_bgr(png_ptr);

	//  allocate a vector of row pointers
	std::vector<uchar *> rowPtrs;
	rowPtrs.resize(height);
	for (int y = 0; y<height; y++)
		rowPtrs[y] = &img.Pixel(0, y, 0);

	// write the whole image
	png_write_image(png_ptr, &rowPtrs[0]);

	// write the additional chunks to the PNG file (not really needed)
	png_write_end(png_ptr, info_ptr);

	png_destroy_write_struct(&png_ptr, (png_infopp) NULL);

    fclose (stream);
}

void ReadFilePNG(CByteImage& img, const char* filename)
{
    // open the PNG input file
    FILE *stream = fopen(filename, "rb");
    if (stream == 0)
        throw CError("ReadFilePNG: could not open %s", filename);

    // first check the eight byte PNG signature
    png_byte pbSig[8];
    fread(pbSig, 1, 8, stream);
	if (!png_check_sig(pbSig, 8)) {
        fclose(stream);
        throw CError("ReadFilePNG: invalid PNG signature");
	}

    // create the two png(-info) structures
    png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL,
      (png_error_ptr)pngfile_error, (png_error_ptr)NULL);

	if (!png_ptr) {
        fclose(stream);
		throw CError("ReadFilePNG: error creating png structure");
	}

	info_ptr = png_create_info_struct(png_ptr);
	if (!info_ptr)
	{
		png_destroy_read_struct(&png_ptr, NULL, NULL);
        fclose(stream);
		throw CError("ReadFilePNG: error creating png structure");
	}

	png_init_io(png_ptr, stream);
	png_set_sig_bytes(png_ptr, 8);

	// read all PNG info up to image data
	png_read_info(png_ptr, info_ptr);

	// get width, height, bit-depth and color-type
	int width, height, bits, colorType, nBands;

	png_uint_32 pwidth, pheight;

	png_get_IHDR(png_ptr, info_ptr, 
		     &pwidth, &pheight,
		     //(png_uint_32 *)&width, (png_uint_32 *)&height,
		     &bits, &colorType, NULL, NULL, NULL);

	width = pwidth;
	height = pheight;

	nBands = (int)png_get_channels(png_ptr, info_ptr);

	if (DEBUG_ImageIOpng)
	fprintf(stderr, " w=%d, h=%d, %2d bits, %s, nB=%d",
		width, height, bits, 
		colorType == PNG_COLOR_TYPE_GRAY ? "gray" :
		colorType == PNG_COLOR_TYPE_PALETTE ? "plt " :
		colorType == PNG_COLOR_TYPE_RGB ? "rgb " :
		colorType == PNG_COLOR_TYPE_RGB_ALPHA ? "rgba" :
		colorType == PNG_COLOR_TYPE_GRAY_ALPHA ? "gr-a" : "??? ",
		nBands);


	// get rid of lower-order byte in 16-bit images
	// TODO: could allow this and read in IntImage in this case...
	if (bits == 16)
		png_set_strip_16(png_ptr);

	// change palette color into RGB
	if (colorType == PNG_COLOR_TYPE_PALETTE)
		png_set_expand(png_ptr);

	// want at least 8 bits
	if (bits < 8)
		png_set_expand(png_ptr);

	// if there is a transparent palette entry, create alpha channel
	if (png_get_valid(png_ptr, info_ptr, PNG_INFO_tRNS))
		png_set_expand(png_ptr);

	// make gray images with alpha channel into RGBA -- TODO: or just ignore alpha?
	if (colorType == PNG_COLOR_TYPE_GRAY_ALPHA)
		// colorType == PNG_COLOR_TYPE_GRAY       // but leave gray images alone
		png_set_gray_to_rgb(png_ptr);

	// set the background color to draw transparent and alpha images over.
	// only needed for gray images with alpha 
	if (colorType == PNG_COLOR_TYPE_GRAY_ALPHA ||
		colorType == PNG_COLOR_TYPE_GRAY) {
		png_color_16 *pBackground;
		if (png_get_bKGD(png_ptr, info_ptr, &pBackground))
			png_set_background(png_ptr, pBackground, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0);
		}

	// if required set gamma conversion
	// this seems to cause problems, so let's just leave gamma alone.
	//double gamma;
	//if (png_get_gAMA(png_ptr, info_ptr, &gamma)) {
	// //fprintf(stderr, "\n reading gamma %lf\n", gamma);
//.........这里部分代码省略.........

int BlendPairs(int argc, const char *argv[])
{
    // Blend a sequence of images given the pairwise transformations
    if (argc < 5)
    {
        printf("usage: %s pairlist.txt outimg.tga blendWidth\n", argv[1]);
        return -1;
    }
    const char *pairlist= argv[2];
    const char *outfile = argv[3];
    float blendWidth    = (float) atof(argv[4]);

    // Open the list of image pairs
    FILE *stream = fopen(pairlist, "r");
    if (stream == 0)
        throw CError("%s: could not open the file %s", argv[1], pairlist);

    // Construct the list of images and translations
    CImagePositionV ipList;
    char line[1024], infile1[1024], infile2[1024];
    // float rel_t[2];
    CTransform3x3 M;

    int n;
    for (n = 0; fgets(line, 1024, stream); n++)
    {
        // Compute the position from the PREVIOUS displacement
        CImagePosition ip;

		if (n == 0) {
			ip.position = CTransform3x3::Translation(0.0, 0.0); /* Ident matrix */
			// ip.position = CTransform3x3::Rotation(45.0); /* Ident matrix */
		} else {
			ip.position = M * ipList[n-1].position;
		}

        // for (int k = 0; k < 2; k++)
        //    ip.position[k] = (n > 0) ? ipList[n-1].position[k] - rel_t[k] : 0.0f;

        // Read the line and push the image onto the list
        // if (sscanf(line, "%s %s %f %f", infile1, infile2,
        //           &rel_t[0], &rel_t[1]) != 4)
        //    throw CError("%s: error reading %s", argv[1], pairlist);

#if 0
        if (sscanf(line, "%s %s %lf %lf", infile1, infile2,
			&(M[0][2]), &(M[1][2])) != 4)
			throw CError("%s: error reading %s\n", argv[1], pairlist);
#else
        if (sscanf(line, "%s %s %lf %lf %lf %lf %lf %lf %lf %lf %lf", infile1, infile2,
			&(M[0][0]), &(M[0][1]), &(M[0][2]), 
			&(M[1][0]), &(M[1][1]), &(M[1][2]), 
			&(M[2][0]), &(M[2][1]), &(M[2][2])) != 11)
			throw CError("%s: error reading %s\n", argv[1], pairlist);
#endif

		// The coordinate system of the ImageLib images has y reflected -- so let's fix that!

		// M[1][2] = -M[1][2];

		ip.imgName = std::string(infile1);
		CByteImage imgTmp;
		ReadFile(imgTmp, infile1);
		ip.img = CByteImage(imgTmp.Shape());
		FlipImage(imgTmp, ip.img);

		ipList.push_back(ip);
	}

	// Read in the last image
	CImagePosition ip;

	// for (int k = 0; k < 2; k++)
	//     ip.position[k] = ipList[n-1].position[k] - rel_t[k];
	ip.position = M * ipList[n-1].position;

	ip.imgName = std::string(infile2);
	CByteImage imgTmp;
	ReadFile(imgTmp, infile2);
	ip.img = CByteImage(imgTmp.Shape());
	FlipImage(imgTmp, ip.img);

	ipList.push_back(ip);
	fclose(stream);

	CByteImage resultTmp = BlendImages(ipList, blendWidth);

	// Flip again
	CByteImage result(resultTmp.Shape());
	FlipImage(resultTmp, result);

	WriteFile(result, outfile);
	return 0;
}