gcx-discussion

[Matei C. <cm...@us...>]

Update of /cvsroot/gcx/gcx/src/ccd
In directory sfp-cvs-1.v30.ch3.sourceforge.com:/tmp/cvs-serv9885/src/ccd

Modified Files:
	Makefile.am badpix.c ccd.h ccd_frame.c dslr.c sources.c warp.c 
Added Files:
	jpeg.c tiff.c 
Log Message:
Merge my modifications to CVS

Index: dslr.c
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/dslr.c,v
retrieving revision 1.14
retrieving revision 1.15
diff -C2 -d -r1.14 -r1.15
*** dslr.c	14 Sep 2009 14:38:43 -0000	1.14
--- dslr.c	7 Mar 2013 20:49:33 -0000	1.15
***************
*** 980,986 ****
  	frame->h = endian_to_host_16(byteorder, mrw->prd->ccd_size_y);
  
- 	frame->pix_size = 4;
- 	frame->pix_format = PIX_FLOAT;
- 
  	if (alloc_frame_data(frame)) {
  		err_printf("read_mrw_file: cannot allocate data for frame\n");
--- 980,983 ----
***************
*** 1067,1073 ****
  	strncpy(frame->name, raw->filename, 255);
  
  	if (endian_to_host_16(byteorder, mrw->prd->ccd_size_x) >
  	    endian_to_host_16(byteorder, mrw->prd->img_size_x)) {
- 
  		frame->x_skip = (endian_to_host_16(byteorder, mrw->prd->ccd_size_x) -
  				 endian_to_host_16(byteorder, mrw->prd->img_size_x)) / 2;
--- 1064,1070 ----
  	strncpy(frame->name, raw->filename, 255);
  
+ #if 0
  	if (endian_to_host_16(byteorder, mrw->prd->ccd_size_x) >
  	    endian_to_host_16(byteorder, mrw->prd->img_size_x)) {
  		frame->x_skip = (endian_to_host_16(byteorder, mrw->prd->ccd_size_x) -
  				 endian_to_host_16(byteorder, mrw->prd->img_size_x)) / 2;
***************
*** 1086,1089 ****
--- 1083,1087 ----
  		fits_add_keyword(frame, "CCDSKIP2", strbuf);
  	}
+ #endif
  
  	if (raw->date_obs)
***************
*** 1437,1443 ****
  #endif
  
- 	frame->pix_size = 4;
- 	frame->pix_format = PIX_FLOAT;
- 
  	/* save color coefficients */
  	frame->rmeta.wbr  = cr2->wb_gain_r;
--- 1435,1438 ----

Index: ccd.h
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/ccd.h,v
retrieving revision 1.38
retrieving revision 1.39
diff -C2 -d -r1.38 -r1.39
*** ccd.h	28 Sep 2011 06:05:47 -0000	1.38
--- ccd.h	7 Mar 2013 20:49:33 -0000	1.39
***************
*** 62,74 ****
  };
  
- // values describing the geometry of a frame from a ccd chip
- struct ccd_geometry {
- 	unsigned w; /* total size of the frame */
- 	unsigned h;
- 	int x_skip; /* number of pixels skipped at the beginning of each line
- 		     * (ref: first active pixel) */
- 	int y_skip; /* number of lines skipped at the beginning of each frame */
- };
- 
  // structure describing a ccd exposure
  struct exp_data {
--- 62,65 ----
***************
*** 174,179 ****
  	int ref_count;
  	char *filename;         // filename of this map, if any
- 	unsigned int x_skip;	// from frame geometry when the map was created
- 	unsigned int y_skip;
  	unsigned int bin_x;
  	unsigned int bin_y;
--- 165,168 ----
***************
*** 186,213 ****
  #define FRAME_NAME_SIZE 256
  
! /*  the big one: struct ccd_frame is the internal representation of a frame or  */
! /*  fits file in cx. Every function that operates on images uses a struct ccd_frame */
! 
! /*  the frames' data areas are allocated to hold DEFAULT_PIX_SIZE bytes/pixel  */
! /*  (usually 4). Functions that process frames in ccd_frame.c and x11ops */
! /*  expect the frames to be float. frame_to_float can be used for the trasnsformation */
! 
! /*  pix_size shows the amount of space allocated for the pixels; if smaller
!  *  pixels are used (e.g. byte/short), pix_size is not affected */
! 
! #define DEFAULT_PIX_SIZE sizeof(float)
! 
! // pixel formats
! #define PIX_FLOAT 1
! #define PIX_BYTE 2
! #define PIX_SHORT 3 	// 16-bit using out native endianess
! #define PIX_16LE 4	// 16-bit little endian
! #define PIX_16BE 5	// 16-bit big endian
! 
! #ifndef LITTLE_ENDIAN
! #ifndef BIG_ENDIAN
! #define LITTLE_ENDIAN
! #endif
! #endif
  
  struct raw_metadata {
--- 175,179 ----
  #define FRAME_NAME_SIZE 256
  
! #define PIXEL_SIZE sizeof(float)
  
  struct raw_metadata {
***************
*** 235,250 ****
  
  
  struct ccd_frame {
! 	int ref_count;
! 	unsigned loaded; /* how much of the frame's data is loaded (for incremental
! 			  * loads and readout - in bytes */
! /*   how many need this frame. When ref_count is 1, a release_frame */
! /*   deallocates the frame. The frame is created with ref_count=1 */
  	int w;		// width of frame in pixels
  	int h;		// height of frame in pixels
! 	int x_skip; 	// number of pixels skipped at the beginning of each line
! 			// (ref: first active pixel)
! 	int y_skip; 	// number of lines skipped at the beginning of each frame
  	unsigned magic; // an unique number identifying the frame (science, dark, flatfield etc)
  	struct exp_data exp;
  	struct im_stats stats;
--- 201,217 ----
  
  
+ /*  the big one: struct ccd_frame is the internal representation of a frame or  */
+ /*  fits file in gcx. Every function that operates on images uses a struct ccd_frame */
+ 
+ /*  the frames' data areas are allocated to hold floats  */
  struct ccd_frame {
! 	int ref_count; //   how many need this frame. When ref_count is 1, a release_frame
! 	               //   deallocates the frame. The frame is created with ref_count=1
! 
  	int w;		// width of frame in pixels
  	int h;		// height of frame in pixels
! 
  	unsigned magic; // an unique number identifying the frame (science, dark, flatfield etc)
+ 
  	struct exp_data exp;
  	struct im_stats stats;
***************
*** 254,266 ****
  			 * normally, it's the amount of data (in bytes)
  			 * in the array */
! 	int pix_size;	// bytes/pixel
! 	int pix_format;	// format pixels are stored in
! 	void *dat;	// the data array
  	FITS_row *var;	/* malloced array of all unrecognized header lines */
  	int nvar;	/* number of var[] */
  	char name[256];	// name of frame
! 	void *rdat;     // rgb image planes
! 	void *gdat;
! 	void *bdat;
  };
  
--- 221,231 ----
  			 * normally, it's the amount of data (in bytes)
  			 * in the array */
! 	float *dat;	// the data array
  	FITS_row *var;	/* malloced array of all unrecognized header lines */
  	int nvar;	/* number of var[] */
  	char name[256];	// name of frame
! 	float *rdat;     // rgb image planes
! 	float *gdat;
! 	float *bdat;
  };
  
***************
*** 274,289 ****
  #define FILE_FRAME	0x0000008
  
- // bits describing the operations on the frame
- 
- #define FRAME_MATH	0x0000010	// frame resulting from computation between frames
- #define FRAME_BADPIX	0x0000020	// frame had had the bad pixels replaced
- #define FRAME_BIASC	0x0000040	// bias was substracted from the frame
- #define FRAME_FLATC	0x0000080	// frame was flatfield corrected
- #define FRAME_STACK	0x0000100	// frame results from a stak op
- 
- #define FRAME_UNCERTAIN 0x1000000	// frame containins uncertain data
- 
  /* frame image structure */
- 
  #define FRAME_HAS_CFA     0x0100000 /* frame contains CFA data */
  #define FRAME_VALID_RGB   0x0200000 /* rgb data was generated from bayer */
--- 239,243 ----
***************
*** 298,302 ****
  
  // bad pixel fixing methods
- 
  #define BADPIX_MEDIAN 1
  #define BADCOLUMN_AVG 0
--- 252,255 ----
***************
*** 577,580 ****
--- 530,535 ----
  int sub_frames (struct ccd_frame *fr, struct ccd_frame *fr1);
  
+ int overscan_correction(struct ccd_frame *fr, double pedestal, int x, int y, int w, int h);
+ 
  
  struct ccd_frame *read_image_file(char *filename, char *ungz, int force_unsigned,
***************
*** 596,599 ****
--- 551,563 ----
  extern int set_color_field(struct ccd_frame *fr);
  
+ // from tiff.c
+ extern int tiff_filename(char *filename);
+ extern struct ccd_frame *read_tiff_file(char *filename);
+ 
+ // from jpeg.c
+ extern int jpeg_filename(char *filename);
+ extern struct ccd_frame *read_jpeg_file(char *filename);
+ 
+ 
  // from badpix.c
  extern int save_bad_pix(struct bad_pix_map *map);
***************
*** 676,682 ****
--- 640,655 ----
  // from warp.c
  
+ enum {
+ 	PAR_RESAMPLE_NEAREST,
+ 	PAR_RESAMPLE_BILINEAR,
+ 	PAR_RESAMPLE_BSPLINE,
+ 	PAR_RESAMPLE_CATMULL,
+ 	PAR_RESAMPLE_MITCHELL
+ };
+ 
  extern int make_shift_ctrans(struct ctrans *ct, double dx, double dy);
  extern int make_roto_translate(struct ctrans *ct, double dx, double dy, double xs, double ys, double rot);
  extern int shift_frame(struct ccd_frame *fr, double dx, double dy);
+ extern int shift_scale_rotate_frame(struct ccd_frame *fr, double dx, double dy, double ds, double dt, int resampling);
  extern int linear_x_shear(struct ccd_frame *in, struct ccd_frame *out, double a, double c);
  extern int linear_y_shear(struct ccd_frame *in, struct ccd_frame *out, double b, double c,
***************
*** 693,700 ****
  
  #endif
- 
- 
- 
- 
- 
- 
--- 666,667 ----

Index: sources.c
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/sources.c,v
retrieving revision 1.21
retrieving revision 1.22
diff -C2 -d -r1.21 -r1.22
*** sources.c	28 Sep 2011 06:05:47 -0000	1.21
--- sources.c	7 Mar 2013 20:49:33 -0000	1.22
***************
*** 50,54 ****
  
  #define RING_MAX 3000
! #define MAX_STAR_CANDIDATES 16384	/* max number of candidates before we find the first star
  					 */
  
--- 50,54 ----
  
  #define RING_MAX 3000
! #define MAX_STAR_CANDIDATES 65536	/* max number of candidates before we find the first star
  					 */
  
***************
*** 102,109 ****
  		dy2 = (iy - yc) * (iy - yc);
  		for (ix = xs; ix < xe; ix++) {
- 			v = get_pixel_luminence(fr, ix, iy);
  			if (((r = (ix - xc) * (ix - xc) + dy2)) 
  			    < rsq1 || r > rsq2)
  				continue;
  			if (v < min || v > max) {
  				nskipped ++;
--- 102,110 ----
  		dy2 = (iy - yc) * (iy - yc);
  		for (ix = xs; ix < xe; ix++) {
  			if (((r = (ix - xc) * (ix - xc) + dy2)) 
  			    < rsq1 || r > rsq2)
  				continue;
+ 
+ 			v = get_pixel_luminence(fr, ix, iy);
  			if (v < min || v > max) {
  				nskipped ++;
***************
*** 611,615 ****
  
  //	d3_printf("extract_stars: frame size is %dx%d\n", fr->w, fr->h);
- //	d3_printf("extract_stars: frame pixel format is %d [%d]\n", fr->pix_format, fr->pix_size);
  
  	minpk = fr->stats.cavg + 2 * sigmas * fr->stats.csigma;
--- 612,615 ----
***************
*** 647,653 ****
  			}
  // check that we have a few connected pixels above the cut
! 			ring_stats(fr, 1.0*x, 1.0*y, 0, 2, 
  				   QUAD1|QUAD2|QUAD3|QUAD4, rsn, skycut, HUGE);
! 			if (rsn->used < NCONN) {
  //				d3_printf("only %d connected pixels found\n", rsn->used);
  				continue;
--- 647,653 ----
  			}
  // check that we have a few connected pixels above the cut
! 			ret = ring_stats(fr, 1.0*x, 1.0*y, 0, 2,
  				   QUAD1|QUAD2|QUAD3|QUAD4, rsn, skycut, HUGE);
! 			if (ret < 0 || rsn->used < NCONN) {
  //				d3_printf("only %d connected pixels found\n", rsn->used);
  				continue;

Index: Makefile.am
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/Makefile.am,v
retrieving revision 1.7
retrieving revision 1.8
diff -C2 -d -r1.7 -r1.8
*** Makefile.am	26 Jul 2011 08:41:54 -0000	1.7
--- Makefile.am	7 Mar 2013 20:49:33 -0000	1.8
***************
*** 10,14 ****
  	edb.c aphot.c worldpos.c sources.c \
  	warp.c recipy.c errlog.c use_dcraw.c \
! 	ccd.h dslr.h
  
  CLEANFILES = *~
--- 10,14 ----
  	edb.c aphot.c worldpos.c sources.c \
  	warp.c recipy.c errlog.c use_dcraw.c \
! 	ccd.h dslr.h tiff.c jpeg.c
  
  CLEANFILES = *~

--- NEW FILE: jpeg.c ---
/*******************************************************************************
  Copyright(c) 2013 Matei Conovici. All rights reserved.

  This program is free software; you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by the Free
  Software Foundation; either version 2 of the License, or (at your option)
  any later version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
  more details.

  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 59
  Temple Place - Suite 330, Boston, MA  02111-1307, USA.

  The full GNU General Public License is included in this distribution in the
  file called LICENSE.

  Contact Information: mco...@gm...
*******************************************************************************/

/*  jpeg.c -- handle JPEG files */

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <setjmp.h>

#include "config.h"
#include "ccd.h"

#ifdef HAVE_LIBJPEG

#include <jpeglib.h>

int jpeg_filename(char *filename)
{
	char *p;

	p = strchr(filename, '.');
	if (!p)
		return -1;

	p++;
	if (!strcasecmp(p, "jpeg") || !strcasecmp(p, "jpg"))
		return 0;

	return -1;
}

static void err_exit(j_common_ptr cinfo)
{
	longjmp(*(jmp_buf *)cinfo->client_data, 1);
}

static int scanline_to_frame(struct ccd_frame *fr, JSAMPLE *buf, int row, int ncolors)
{
	float *dp;
	int color, i;

	for (color = 0; color < ncolors; color++) {
		dp = fr->dat;
		if (ncolors != 1) {
			switch (color) {
			case 0:
				dp = fr->rdat;
				break;
			case 1:
				dp = fr->gdat;
				break;
			case 2:
				dp = fr->bdat;
				break;
			default:
				return 0;
			}
		}

		dp += row * fr->w;

		for (i = 0; i < fr->w; i++) {
			*dp++ = (float) *((JSAMPLE *) buf + i * ncolors + color);
		}
	}

	return 0;
}


struct ccd_frame *read_jpeg_file(char *filename)
{
	FILE *fin;
	struct ccd_frame *fr = NULL;
	struct jpeg_decompress_struct cinfo;
	struct jpeg_error_mgr jerr;
	int i;
	JSAMPLE *buf = NULL;
	jmp_buf errbuf;

	if ((fin = fopen(filename, "rb")) == NULL)
		return NULL;

	cinfo.err = jpeg_std_error(&jerr);

	/* override exit fn */
	cinfo.err->error_exit = err_exit;
	cinfo.client_data = &errbuf;

	/* have you seen my axe ? I have some error recovery to do */
	if (setjmp(errbuf)) {
		release_frame(fr);
		fr = NULL;
		goto out;
	}

	jpeg_create_decompress(&cinfo);
	jpeg_stdio_src(&cinfo, fin);

	if (jpeg_read_header(&cinfo, TRUE) != JPEG_HEADER_OK)
		goto out;

	jpeg_start_decompress(&cinfo);

	if (cinfo.output_components != 1 && cinfo.output_components != 3)
		goto out;

	fr = new_frame(cinfo.output_width, cinfo.output_height);
	if (cinfo.output_components == 3) {
		free(fr->dat);
		alloc_frame_rgb_data(fr);
		fr->magic |= FRAME_VALID_RGB;
	}

	if ((buf = malloc(cinfo.output_width * cinfo.output_components * sizeof(JSAMPLE))) == NULL) {
		release_frame(fr);
		fr = NULL;
		goto out;
	}

	for (i = 0; i < cinfo.output_height; i++) {
		jpeg_read_scanlines(&cinfo, &buf, 1);
		scanline_to_frame(fr, buf, i, cinfo.output_components);
	}

out:
	free(buf);
	fclose(fin);
	jpeg_destroy_decompress(&cinfo);

	return fr;

}

#endif

Index: badpix.c
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/badpix.c,v
retrieving revision 1.22
retrieving revision 1.23
diff -C2 -d -r1.22 -r1.23
*** badpix.c	27 Sep 2009 15:19:48 -0000	1.22
--- badpix.c	7 Mar 2013 20:49:33 -0000	1.23
***************
*** 122,127 ****
  
  error:
- 	map->x_skip = fr->x_skip;
- 	map->y_skip = fr->y_skip;
  	map->bin_x = fr->exp.bin_x;
  	map->bin_y = fr->exp.bin_y;
--- 122,125 ----
***************
*** 154,159 ****
  
  	fprintf(fp, "pixels %d\n", map->pixels);
- 	fprintf(fp, "x_skip %d\n", map->x_skip);
- 	fprintf(fp, "y_skip %d\n", map->y_skip);
  	fprintf(fp, "bin_x %d\n", map->bin_x);
  	fprintf(fp, "bin_y %d\n", map->bin_y);
--- 152,155 ----
***************
*** 191,200 ****
  	if (ret != 1)
  		goto bad_format;
- 	ret = fscanf(fp, " x_skip %d", &map->x_skip);
- 	if (ret != 1)
- 		goto bad_format;
- 	ret = fscanf(fp, " y_skip %d", &map->y_skip);
- 	if (ret != 1)
- 		goto bad_format;
  	ret = fscanf(fp, " bin_x %d", &map->bin_x);
  	if (ret != 1)
--- 187,190 ----
***************
*** 521,531 ****
  		/* regular BW image */
  		for (i = 0; i < map->pixels; i++) {
- 			if (fr->exp.bin_y != 0 && fr->exp.bin_x != 0) {
- 				frx = map->pix[i].x - fr->x_skip / fr->exp.bin_x;
- 				fry = map->pix[i].y - fr->y_skip / fr->exp.bin_y;
- 			} else {
  				frx = map->pix[i].x;
  				fry = map->pix[i].y;
- 			}
  
  			if (frx > 1 && frx < fr->w - 2 && fry > 1 && fry < fr->h - 2) {
--- 511,516 ----
***************
*** 551,553 ****
  	return 0;
  }
- 
--- 536,537 ----

Index: ccd_frame.c
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/ccd_frame.c,v
retrieving revision 1.49
retrieving revision 1.50
diff -C2 -d -r1.49 -r1.50
*** ccd_frame.c	3 Dec 2011 21:48:07 -0000	1.49
--- ccd_frame.c	7 Mar 2013 20:49:33 -0000	1.50
***************
*** 23,30 ****
  
  /*  ccd_frame.c: frame operation functions */
- /*  $Revision$ */
- /*  $Date$ */
- 
- /* Many functions here assume we have float frames - this must be fixed */
  
  #define _GNU_SOURCE
--- 23,26 ----
***************
[...1166 lines suppressed...]
- 	case PIX_16LE:
- 		op = (float *)fr->dat + all - 1;
- 		ip = fr->dat + all * 2 - 2;
- 		for (i=0; i< all; i++) {
- 			*op-- = (256.0 * (* (unsigned char *) (ip+1))
- 				 + 1.0 * (* (unsigned char *) (ip)));
- 			ip -= 2;
- 		}
- 		fr->pix_format = PIX_FLOAT;
- 		return 0;
- #endif /* LITTLE_ENDIAN */
- 		default:
- 		err_printf("cannot convert unknown format %d to float\n");
- 		return -1;
- 	}
- 
- }
  
  /* get a double field
--- 1399,1402 ----

--- NEW FILE: tiff.c ---
/*******************************************************************************
  Copyright(c) 2013 Matei Conovici. All rights reserved.

  This program is free software; you can redistribute it and/or modify it
  under the terms of the GNU General Public License as published by the Free
  Software Foundation; either version 2 of the License, or (at your option)
  any later version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
  more details.

  You should have received a copy of the GNU General Public License along with
  this program; if not, write to the Free Software Foundation, Inc., 59
  Temple Place - Suite 330, Boston, MA  02111-1307, USA.

  The full GNU General Public License is included in this distribution in the
  file called LICENSE.

  Contact Information: mco...@gm...
*******************************************************************************/

/*  tiff.c -- handle TIFF files */


#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>

#include "config.h"
#include "ccd.h"

#ifdef HAVE_LIBTIFF

#include <tiffio.h>


int tiff_filename(char *filename)
{
	char *p;

	p = strchr(filename, '.');
	if (!p)
		return -1;

	p++;
	if (!strcasecmp(p, "tiff") || !strcasecmp(p, "tif"))
		return 0;

	return -1;
}

static int scanline_separate_to_frame(struct ccd_frame *fr,
				      tdata_t buf, uint32_t row,
				      uint16_t color, uint16_t ncolors,
				      uint16_t fmt, uint16_t bpp)
{
	float *dp = fr->dat;
	int i;

	if (ncolors != 1) {
		switch (color) {
		case 0:
			dp = fr->rdat;
			break;
		case 1:
			dp = fr->gdat;
			break;
		case 2:
			dp = fr->bdat;
			break;
		default:
			// ignore alpha
			return 0;
		}
	}

	dp += row * fr->w;

	for (i = 0; i < fr->w; i++) {
		switch (fmt) {
		case SAMPLEFORMAT_INT:
			switch (bpp) {
			case 8:
				*dp++ = (float) *((int8_t *) buf + i);
				break;
			case 16:
				*dp++ = (float) *((int16_t *) buf + i);
				break;
			case 32:
				*dp++ = (float) *((int32_t *) buf + i);
				break;
			default:
				return -1;
			}
			break;

		case SAMPLEFORMAT_IEEEFP:
			if (bpp != 32)
				return -1;

			*dp++ = *((float *) buf + i);
			break;

		case SAMPLEFORMAT_UINT:
		default:
			switch (bpp) {
			case 8:
				*dp++ = (float) *((uint8_t *) buf + i);
				break;
			case 16:
				*dp++ = (float) *((uint16_t *) buf + i);
				break;
			case 32:
				*dp++ = (float) *((uint32_t *) buf + i);
				break;
			default:
				return -1;
			}
		}
	}

	return 0;
}

static int scanline_contig_to_frame(struct ccd_frame *fr,
				     tdata_t buf, uint32_t row, uint16_t ncolors,
				     uint16_t fmt, uint16_t bpp)
{
	float *dp = fr->dat;
	int i, color;

	for (color = 0; color < ncolors; color++) {
		if (ncolors != 1) {
			switch (color) {
			case 0:
				dp = fr->rdat;
				break;
			case 1:
				dp = fr->gdat;
				break;
			case 2:
				dp = fr->bdat;
				break;
			default:
				// ignore alpha
				return 0;
			}
		}

		dp += row * fr->w;

		for (i = 0; i < fr->w; i++) {
			switch (fmt) {
			case SAMPLEFORMAT_INT:
				switch (bpp) {
				case 8:
					*dp++ = (float) *((int8_t *) buf + i * ncolors + color);
					break;
				case 16:
					*dp++ = (float) *((int16_t *) buf + i * ncolors + color);
					break;
				case 32:
					*dp++ = (float) *((int32_t *) buf + i * ncolors + color);
					break;
				default:
					return -1;
				}
				break;

			case SAMPLEFORMAT_IEEEFP:
				if (bpp != 32)
					return -1;

				*dp++ = *((float *) buf + i * ncolors + color);
				break;

			case SAMPLEFORMAT_UINT:
			default:
				switch (bpp) {
				case 8:
					*dp++ = (float) *((uint8_t *) buf + i * ncolors + color);
					break;
				case 16:
					*dp++ = (float) *((uint16_t *) buf + i * ncolors + color);
					break;
				case 32:
					*dp++ = (float) *((uint32_t *) buf + i * ncolors + color);
					break;
				default:
					return -1;
				}
			}
		}
	}

	return 0;
}



struct ccd_frame *read_tiff_file(char *filename)
{
	struct ccd_frame *fr = NULL;
	TIFF *tif;
	uint32_t w, h, row;
	uint16_t nsamples, config, bpp, fmt, color;
	tdata_t buf;
	int ret = 0;

	if ((tif = TIFFOpen(filename, "r")) == NULL)
		return NULL;

	ret += TIFFGetField(tif, TIFFTAG_IMAGEWIDTH, &w);
	ret += TIFFGetField(tif, TIFFTAG_IMAGELENGTH, &h);
	ret += TIFFGetField(tif, TIFFTAG_SAMPLESPERPIXEL, &nsamples);
	ret += TIFFGetField(tif, TIFFTAG_PLANARCONFIG, &config);

	/* we understand only the simplest of TIFF files */
	if (ret != 4)
		goto out;

	fmt = SAMPLEFORMAT_UINT;
	TIFFGetField(tif, TIFFTAG_SAMPLEFORMAT, &fmt);

	bpp = 16;
	TIFFGetField(tif, TIFFTAG_BITSPERSAMPLE, &bpp);

	if (nsamples != 1 && nsamples != 3 && nsamples != 4)
		goto out;

	/* FIXME: should attempt to extract some metadata, like date/time, EXIF data, etc */

	fr = new_frame(w, h);
	if (nsamples != 1) {
		free(fr->dat);
		alloc_frame_rgb_data(fr);
		fr->magic |= FRAME_VALID_RGB;
	}

	buf = _TIFFmalloc(TIFFScanlineSize(tif));
	if (config == PLANARCONFIG_CONTIG) {
		/* packed RGB(A) */
		for (row = 0; row < h; row++) {
			if (TIFFReadScanline(tif, buf, row, 0) != 1)
				goto out_err;

			if (scanline_contig_to_frame(fr, buf, row, nsamples, fmt, bpp))
				goto out_err;
		}
	} else if (config == PLANARCONFIG_SEPARATE) {
		/* separate planes */
		for (color = 0; color < nsamples; color++) {
			for (row = 0; row < h; row++) {
				if (TIFFReadScanline(tif, buf, row, color) != 1)
					goto out_err;

				if (scanline_separate_to_frame(fr, buf, row, color, nsamples, fmt, bpp))
					goto out_err;
			}
		}
	} else
		goto out_err;

	_TIFFfree(buf);

out:
	TIFFClose(tif);

	return fr;

out_err:
	_TIFFfree(buf);
	TIFFClose(tif);
	release_frame(fr);

	return NULL;
}

#endif

Index: warp.c
===================================================================
RCS file: /cvsroot/gcx/gcx/src/ccd/warp.c,v
retrieving revision 1.23
retrieving revision 1.24
diff -C2 -d -r1.23 -r1.24
*** warp.c	19 Aug 2009 18:40:19 -0000	1.23
--- warp.c	7 Mar 2013 20:49:33 -0000	1.24
***************
*** 274,278 ****
  		dpi = get_color_plane(fr, plane_iter);
  		dpo = get_color_plane(new, plane_iter);
! 		memcpy(dpo, dpi, all * fr->pix_size);
  	}
  	fr->stats.statsok = 0;
--- 274,278 ----
  		dpi = get_color_plane(fr, plane_iter);
  		dpo = get_color_plane(new, plane_iter);
! 		memcpy(dpo, dpi, all * PIXEL_SIZE);
  	}
  	fr->stats.statsok = 0;
***************
*** 330,334 ****
  }
  
- /*
  // compute the coordinate transform
  static void do_ctrans(struct ctrans *ct, double u, double v, double *x, double *y)
--- 330,333 ----
***************
*** 360,364 ****
  	}
  }
! */
  // do linear shear transform in the x direction
  // x = c ( u - a * v )
--- 359,363 ----
  	}
  }
! 
  // do linear shear transform in the x direction
  // x = c ( u - a * v )
***************
*** 665,668 ****
--- 664,935 ----
  }
  
+ 
+ /*
+ static void coord_transform(double dx, double dy, double ds, double dt, double x, double y, double *xx, double *yy)
+ {
+ 	*xx = dx + (x * ds - dx) * cos(dt) + (y * ds - dy) * sin(dt);
+ 	*yy = dy - (x * ds - dx) * sin(dt) + (y * ds - dy) * cos(dt);
+ }
+ */
+ 
+ #if 0
+ 	ct->u0 = dx;
+ 	ct->v0 = dy;
+ 	sa = sin(-rot);
+ 	ca = cos(-rot);
+ 	ct->a[1][0] = 1.0 / xs * ca;
+ 	ct->a[0][1] = - 1.0 / xs * sa;
+ 	ct->b[0][1] = 1.0 / ys * ca;
+ 	ct->b[1][0] = 1.0 / ys * sa;
+ 	return 0;
+ #endif
+ 
+ struct cmatrix {
+ 	double d[3][3];
+ 	double i[3][3];
+ };
+ 
+ void make_cmatrix(struct cmatrix *mat, double dx, double dy, double ds, double dt)
+ {
+ 	int i, j;
+ 	double st, ct;
+ 	double det;
+ 
+ 	st = sin(dt);
+ 	ct = cos(dt);
+ 
+ 	for (i = 0; i < 3; i++)
+ 		for (j = 0; j < 3; j++)
+ 			mat->d[i][j] = (i == j) ? 1.0 : 0.0;
+ 
+ 	mat->d[0][2] = dx;
+ 	mat->d[1][2] = dy;
+ 
+ 	mat->d[0][0] = ct * ds;
+ 	mat->d[0][1] = - st * ds;
+ 	mat->d[1][0] = st * ds;
+ 	mat->d[1][1] = ct * ds;
+ 
+ 	det = mat->d[0][0] * mat->d[1][1] - mat->d[0][1] * mat->d[1][0];
+ 	if (fabs(det) < 10e-30) {
+ 		err_printf("singular matrix\n");
+ 		return;
+ 	}
+ 
+ 	mat->i[0][0] = 1.0 / det * mat->d[1][1];
+ 	mat->i[0][1] = - 1.0 / det * mat->d[0][1];
+ 	mat->i[0][2] = 1.0 / det * (mat->d[0][1] * mat->d[1][2] - mat->d[0][2] * mat->d[1][1]);
+ 
+ 	mat->i[1][0] = - 1.0 / det * mat->d[1][0];
+ 	mat->i[1][1] = 1.0 / det * mat->d[0][0];
+ 	mat->i[1][2] = 1.0 / det * (mat->d[0][2] * mat->d[1][0] - mat->d[0][0] * mat->d[1][2]);
+ 
+ 	mat->i[2][0] = 0.0;
+ 	mat->i[2][1] = 0.0;
+ 	mat->i[2][2] = 1.0 / det * (mat->d[0][0] * mat->d[1][1] - mat->d[0][1] * mat->d[1][0]);
+ 
+ }
+ 
+ void do_cmatrix_d(struct cmatrix *mat, double u, double v, double *x, double *y)
+ {
+ 	*x = mat->d[0][0] * u + mat->d[0][1] * v + mat->d[0][2];
+ 	*y = mat->d[1][0] * u + mat->d[1][1] * v + mat->d[1][2];
+ }
+ 
+ static inline void do_cmatrix_i(struct cmatrix *mat, double u, double v, double *x, double *y)
+ {
+ 	*x = mat->i[0][0] * u + mat->i[0][1] * v + mat->i[0][2];
+ 	*y = mat->i[1][0] * u + mat->i[1][1] * v + mat->i[1][2];
+ }
+ 
+ 
+ static inline double cubic_weight(double x, double a, double b)
+ {
+ 	double v = fabs(x);
+ 	double q = 0.0;
+ 
+ 	if (v < 1) {
+ 		q = (-6*a - 9*b + 12) * v * v *v +
+ 		    (6*a + 12*b - 18) * v * v -
+ 		    2*b + 6;
+ 	} else if (v < 2) {
+ 		q = (-6*a - b) * v * v * v +
+ 		    (30*a + 6*b) * v * v +
+ 		    (-48*a - 12*b) * v +
+ 		    24*a + 8*b;
+ 	}
+ 
+ 	return 1.0 / 6.0 * q;
+ }
+ 
+ static inline double interpolate_pixel_cubic(float *data, int w, int h, double x, double y, double a, double b)
+ {
+ 	int x0, y0;
+ 	double p, q;
+ 	int i, j, u, v;
+ 
+ 	if (x < 3 || x >= w - 3 || y < 3 || y >= h - 3)
+ 		return 0.0;
+ 
+ 	x0 = (int) floor(x);
+ 	y0 = (int) floor(y);
+ 
+ 	// n = extent of kernel
+ 	// q = 0
+ 	// for j = 0 ... 2n-1 do // iterate over 2n line
+ 	//     let v <- vj = |y0| + j - n + 1
+ 	//     let p <- 0
+ 	//     for i = 0 ... 2n-1 do // iterate over 2n cols
+ 	//         let u <- ui = |x0| + i - n +1
+ 	//         let p <- p + I(u,v) * w(x0-u)
+ 	//     done
+ 	//     let q = q + p * w(y0-v)
+ 	// done
+ 
+ 	q = 0.0;
+ 	for (j = 0; j < 4; j++) {
+ 		p = 0.0;
+ 		v = y0 + j - 1;
+ 		for (i = 0; i < 4; i++) {
+ 			u = x0 + i - 1;
+ 			p += data[v * w + u] * cubic_weight(x - (double) u, a, b);
+ 		}
+ 
+ 		q += p * cubic_weight(y - (double) v, a, b);
+ 	}
+ 
+ 
+ 	return q;
+ }
+ 
+ int shift_scale_rotate_frame(struct ccd_frame *fr, double dx, double dy, double ds, double dt, int resampling)
+ {
+ 	struct cmatrix cm;
+ 	int i, j;
+ 	int xx, yy;
+ 	double x, y;
+ 	double a = 0.0, b = 0.0, c, d;
+ 	int cubic = 1;
+ 	int index = 0;
+ 	float *ddat = NULL, *sdat = NULL;
+ 	int plane_iter;
+ 
+ 	make_cmatrix(&cm, dx, dy, ds, dt);
+ 
+ 	struct ccd_frame *nfr = new_frame_fr(fr, fr->w, fr->h);
+ 
+ 	switch (resampling) {
+ 
+ 	case PAR_RESAMPLE_NEAREST:
+ 		for (i = 0; i < nfr->h; i++) {
+ 			for (j = 0; j < nfr->w; j++) {
+ 				do_cmatrix_i(&cm, j, i, &x, &y);
+ 
+ 				yy = (int) floor(y + 0.5);
+ 				xx = (int) floor(x + 0.5);
+ 
+ 				if (yy >= 0 && yy < fr->h && xx >= 0 && xx < fr->w) {
+ 					plane_iter = 0;
+ 					while ((plane_iter = color_plane_iter(fr, plane_iter))) {
+ 						sdat = get_color_plane(fr, plane_iter);
+ 						ddat = get_color_plane(nfr, plane_iter);
+ 
+ 						ddat[index] = sdat[yy*fr->w + xx];
+ 					}
+ 					index++;
+ 
+ 				} else {
+ 					index++;
+ 				}
+ 			}
+ 		}
+ 		cubic = 0;
+ 		break;
+ 
+ 	case PAR_RESAMPLE_BILINEAR:
+ 		for (i = 0; i < nfr->h; i++) {
+ 			for (j = 0; j < nfr->w; j++) {
+ 				do_cmatrix_i(&cm, j, i, &x, &y);
+ 
+ 				if (x < 1 || x > fr->w - 1 || y < 1 || y > fr->h - 1)
+ 					index++;
+ 				else {
+ 					a = x - floor(x);
+ 					b = 1 - a;
+ 					c = y - floor(y);
+ 					d = 1 - c;
+ 
+ 					int xx = (int) floor(x);
+ 					int yy = (int) floor(y);
+ 
+ 					plane_iter = 0;
+ 					while ((plane_iter = color_plane_iter(fr, plane_iter))) {
+ 						sdat = get_color_plane(fr, plane_iter);
+ 						ddat = get_color_plane(nfr, plane_iter);
+ 
+ 
+ 						ddat[index] = b * d * sdat[yy * fr->w + xx] +
+ 							      a * d * sdat[yy * fr->w + xx + 1] +
+ 							      b * c * sdat[(yy + 1) * fr->w + xx] +
+ 							      a * c * sdat[(yy + 1) * fr->w + xx + 1];
+ 					}
+ 					index++;
+ 				}
+ 			}
+ 		}
+ 		cubic = 0;
+ 		break;
+ 
+ 	case PAR_RESAMPLE_BSPLINE:	     /* Cubic B-spline */
+ 		a = 1.0; b = 0.0;
+ 		break;
+ 
+ 	case PAR_RESAMPLE_CATMULL:	     /* Catmull-Rom */
+ 		a = 0.0; b = 0.5;
+ 		break;
+ 
+ 	case PAR_RESAMPLE_MITCHELL:	     /* Mitchell-Netravali */
+ 		a = 1/3.0; b = 1/3.0;
+ 		break;
+ 
+ 	default:
+ 		return -1;
+ 	}
+ 
+ 	if (cubic) {
+ 		for (i = 0; i < nfr->h; i++) {
+ 			for (j = 0; j < nfr->w; j++) {
+ 				do_cmatrix_i(&cm, j, i, &x, &y);
+ 
+ 				/* a = 1.0, b = 0.0     -> Cubic B-spline
+ 				   a = 0.0, b = 0.5     -> Catmull-Rom
+ 				   a = 1/3.0, b = 1/3.0 -> Mitchell-Netravali
+ 				*/
+ 				plane_iter = 0;
+ 				while ((plane_iter = color_plane_iter(fr, plane_iter))) {
+ 					sdat = get_color_plane(fr, plane_iter);
+ 					ddat = get_color_plane(nfr, plane_iter);
+ 
+ 					ddat[index] = interpolate_pixel_cubic(sdat, fr->w, fr->h, x, y, a, b);
+ 				}
+ 				index++;
+ 			}
+ 		}
+ 	}
+ 
+ 	fr->stats.statsok = 0;
+ 
+ 	free(fr->dat);  fr->dat  = nfr->dat;  nfr->dat = NULL;
+ 	free(fr->rdat); fr->rdat = nfr->rdat; nfr->rdat = NULL;
+ 	free(fr->gdat); fr->gdat = nfr->gdat; nfr->gdat = NULL;
+ 	free(fr->bdat); fr->bdat = nfr->bdat; nfr->bdat = NULL;
+ 
+ 	release_frame(nfr);
+ 
+ 	return 0;
+ 
+ }
+ 
+ 
  // create a gaussian filter kernel of given sigma
  // requires a prealloced table of floats of suitable size (size*size)
***************
*** 710,712 ****
  	return 0;
  }
- 
--- 977,978 ----