[PD-cvs] externals/pdp/system/image Makefile, 1.2, 1.3 pdp_imageproc_common.c, 1.2, 1.3 pdp_imageproc_mmx.c, 1.2, 1.3 pdp_imageproc_portable.c, 1.2, 1.3 pdp_llconv.c, 1.2, 1.3 pdp_llconv_mmx.c, 1.2, 1.3 pdp_llconv_portable.c, 1.2, 1.3 pdp_resample.c, 1.2, 1.3

Hans-Christoph Steiner eighthave at users.sourceforge.net
Fri Dec 16 02:05:41 CET 2005


Update of /cvsroot/pure-data/externals/pdp/system/image
In directory sc8-pr-cvs1.sourceforge.net:/tmp/cvs-serv6756/system/image

Added Files:
	Makefile pdp_imageproc_common.c pdp_imageproc_mmx.c 
	pdp_imageproc_portable.c pdp_llconv.c pdp_llconv_mmx.c 
	pdp_llconv_portable.c pdp_resample.c 
Log Message:
checking in pdp 0.12.4 from http://zwizwa.fartit.com/pd/pdp/pdp-0.12.4.tar.gz

--- NEW FILE: pdp_llconv.c ---
/*
 *   Pure Data Packet system implementation. : low level format conversion code
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */

/* this file contains low level image conversion code 
   nominated as "the ugliest part of pdp"
   some code is mmx, most is not. */

/* seem's there's some confusion between rgb and bgr formats.
   not that it matters much which is supposed to be the "real"
   rgb or bgr, but opengl and v4l seem to disagree on endianness
   grounds.. */

#include "pdp_llconv.h"
#include "pdp_mmx.h"
#include "pdp_post.h"


/* all symbols are C style */
#ifdef __cplusplus
extern "C"
{
#endif


#define CLAMP8(x) (((x)<0) ? 0 : ((x>255)? 255 : (x)))
#define CLAMP16(x) (((x)<-0x7fff) ? -0x7fff : ((x>0x7fff) ? 0x7fff : (x)))
#define FP(x) ((int)(((float)(x)) * 256.0f))

#define CLAMP CLAMP8

/* some prototypes for functions defined elsewhere */
void llconv_yvu_planar_s16u8(short int *src, unsigned char *dst, unsigned int nbpixels);
void llconv_yuv_planar_u8s16(unsigned char* source, short int *dest, int nbpixels);
void llconv_grey_s16u8(short int *src, unsigned char *dst, unsigned int nbpixels);
void llconv_yvu_planar_u8s16(unsigned char* source, short int *dest, int nbpixels);


static inline int rgb2y(int r, int g, int b)
{
    return (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
}
static inline int rgb2v(int r, int g, int b)
{
    return (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b) + FP(128);
}
static inline int rgb2u(int r, int g, int b)
{
    return -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b) + FP(128);
}


/* swap top to bottom */
static inline void _exchange_row(char *row1, char *row2, int size)
{
    int mask = ~(sizeof(int)-1);
    int *irow1 = (int *)row1;
    int *irow2 = (int *)row2;

    /* transfer words */
    while (size & mask){
	int tmp = *irow1;
	*irow1++ = *irow2;
	*irow2++ = tmp;
	size -= sizeof(int);
    }

    row1 = (char *)irow1;
    row2 = (char *)irow2;

    /* transfer rest bytes */
    while (size){
	int tmp = *row1;
	*row1++ = *row2;
	*row2++ = tmp;
	size--;
    }
}

void pdp_llconv_flip_top_bottom(char *data, int width, int height, int pixelsize)
{
    int linesize = width * pixelsize;
    int i;
    char *row1 = data;
    char *row2 = data + linesize * (height-1);

    if (height <= 1) return;
    if (width  <= 0) return;

    while (row1 < row2){
	_exchange_row(row1, row2, linesize);
	row1 += linesize;
	row2 -= linesize;
    }
}

/* "standard" 8 bit conversion routine */
static void llconv_rgb2yvu(unsigned char* src, unsigned char* dst, int nbpixels)
{
    int r,g,b,y,v,u,i;
    for (i=0; i<nbpixels; i++){
	r = src[0];
	g = src[1];
	b = src[2];

	y = rgb2y(r,g,b);
	v = rgb2v(r,g,b);
	u = rgb2u(r,g,b);

	dst[0] = CLAMP(y>>8);
	dst[1] = CLAMP(v>>8);
	dst[2] = CLAMP(u>>8);

	src += 3;
	dst += 3;
    }
}

static void llconv_yvu16planar2rgbpacked(short int *src, unsigned char *dst, int w, int h)
{

/*
B = 1.164(Y - 16)                   + 2.018(U - 128)
G = 1.164(Y - 16) - 0.813(V - 128) - 0.391(U - 128)
R = 1.164(Y - 16) + 1.596(V - 128)}
*/

    int r,g,b,y,u,v,b1,g1,r1,y1,xoff,yoff;
    int size = w*h;
    int voffset = size;
    int uoffset = size + (size>>2);
    int rgboff;
    int rgbw = w*3;
    int lumoff = 0;
    int chromoff = 0;

    for(yoff=0; yoff<w*h; yoff+=2*w){
	for(xoff=0; xoff<w; xoff+=2){

	    /* calculate offsets */
	    rgboff = 3 * (xoff + yoff);
	    lumoff = xoff + yoff;
	    chromoff =  (xoff >> 1) + (yoff >> 2);
	    
	    /* get uv values */
	    v = src[voffset + chromoff];
	    u = src[uoffset + chromoff];

	    /* calculate chroma contrib for 2x2 pixblock */
	    b1 = FP(2.018) * u;
	    g1 = FP(-0.813) * v + FP(-0.391) * u;
	    r1 = FP(1.596) * v;

	    /* TOP LEFT */

	    /* top left luma contrib */
	    y = src[lumoff] << 1;
	    y1 = FP(1.164) * y;
	    y1 -= FP(16*256);
	    
	    b = (b1 + y1)>>16;
	    g = (g1 + y1)>>16;
	    r = (r1 + y1)>>16; 

	    /* store top left rgb pixel */
	    dst[rgboff+0] = CLAMP8(r);
	    dst[rgboff+1] = CLAMP8(g);
	    dst[rgboff+2] = CLAMP8(b);

	    /* TOP RIGHT */

	    /* top right luma contrib */
	    y = src[lumoff + 1] << 1;
	    y1 = FP(1.164) * y;
	    y1 -= FP(16*256);
	    
	    b = (b1 + y1)>>16;
	    g = (g1 + y1)>>16;
	    r = (r1 + y1)>>16; 

	    /* store top right rgb pixel */
	    dst[rgboff+3] = CLAMP8(r);
	    dst[rgboff+4] = CLAMP8(g);
	    dst[rgboff+5] = CLAMP8(b);


	    /* BOTTOM LEFT */

	    /* bottom left luma contrib */
	    y = src[lumoff+w] << 1;
	    y1 = FP(1.164) * y;
	    y1 -= FP(16*256);
	    
	    b = (b1 + y1)>>16;
	    g = (g1 + y1)>>16;
	    r = (r1 + y1)>>16; 

	    /* store bottom left rgb pixel */
	    dst[rgboff+rgbw+0] = CLAMP8(r);
	    dst[rgboff+rgbw+1] = CLAMP8(g);
	    dst[rgboff+rgbw+2] = CLAMP8(b);

	    /* BOTTOM RIGHT */

	    /* bottom right luma contrib */
	    y = src[lumoff + w + 1] << 1;
	    y1 = FP(1.164) * y;
	    y1 -= FP(16*256);
	    
	    b = (b1 + y1)>>16;
	    g = (g1 + y1)>>16;
	    r = (r1 + y1)>>16; 

	    /* store bottom right rgb pixel */
	    dst[rgboff+rgbw+3] = CLAMP8(r);
	    dst[rgboff+rgbw+4] = CLAMP8(g);
	    dst[rgboff+rgbw+5] = CLAMP8(b);

	}

    }

}



/* common 8 bit rgb -> 16 bit yvu */
inline static void llconv_rgb2yvu_planar16sub_indexed(unsigned char* src, short int* dst, int w, int h, int ir, int ig, int ib, int stride)
{
    int r,g,b,y,v,u,i,j,k;
    int size = w*h;

    int voffset = size;
    int uoffset = size + (size>>2);


    int loffset = w * stride;

    k=0;
    for (j=0; j<w*h; j+=(w<<1)){
	k = stride * j;
	for (i=0; i<w; i+=2){


	    // well, this seems to work... strange though
	    r = src[k+ir];
	    g = src[k+ig];
	    b = src[k+ib];
	    
	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v =  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u = -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[i+j] = CLAMP16(y >> 1);

	    r = src[k+stride+ir];
	    g = src[k+stride+ig];
	    b = src[k+stride+ib];
	    
	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v +=  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u += -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[i+j+1] = CLAMP16(y >> 1);



	    r = src[loffset + k+ir];
	    g = src[loffset + k+ig];
	    b = src[loffset + k+ib];
	    
	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v =  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u = -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[w+i+j] = CLAMP16(y >> 1);

	    r = src[loffset + k+stride+ir];
	    g = src[loffset + k+stride+ig];
	    b = src[loffset + k+stride+ib];
	    
	    k += 2 * stride;

	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v +=  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u += -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[w+i+j+1] = CLAMP16(y >> 1);

	    dst[uoffset+ (i>>1) + (j>>2)] = (CLAMP16(u >> 1));
	    dst[voffset+ (i>>1) + (j>>2)] = (CLAMP16(v >> 1));
	}
    }
}

/* 8 bit rgb to 16 bit planar subsampled yvu */
static void llconv_rgb2yvu_planar16sub(unsigned char* src, short int* dst, int w, int h)
{
    llconv_rgb2yvu_planar16sub_indexed(src,dst,w,h,0,1,2,3);
}

/* 8 bit rgba to 16 bit planar subsampled yvu */
static void llconv_rgba2yvu_planar16sub(unsigned char* src, short int* dst, int w, int h)
{
    llconv_rgb2yvu_planar16sub_indexed(src,dst,w,h,0,1,2,4);
}

/* 8 bit bgr to 16 bit planar subsampled yvu */
static void llconv_bgr2yvu_planar16sub(unsigned char* src, short int* dst, int w, int h)
{
    llconv_rgb2yvu_planar16sub_indexed(src,dst,w,h,2,1,0,3);
}

/* 8 bit bgra to 16 bit planar subsampled yvu */
static void llconv_bgra2yvu_planar16sub(unsigned char* src, short int* dst, int w, int h)
{
    llconv_rgb2yvu_planar16sub_indexed(src,dst,w,h,2,1,0,4);
}


/* 8 bit rgb to 8 bit planar subsampled yvu */
static void llconv_rgb2yvu_planar8sub(unsigned char* src, unsigned char *dst, int w, int h)
{
    int r,g,b,y,v,u,i,j,k;
    int size = w*h;

    int voffset = size;
    int uoffset = size + (size>>2);


    int loffset = w * 3;

    k=0;
    for (j=0; j<w*h; j+=(w<<1)){
	k = 3 * j;
	for (i=0; i<w; i+=2){


	    // well, this seems to work... strange though
	    r = src[k];
	    g = src[k+1];
	    b = src[k+2];
	    
	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v =  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u = -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[i+j] = CLAMP8(y >> 8);

	    r = src[k+3];
	    g = src[k+4];
	    b = src[k+5];
	    
	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v +=  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u += -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[i+j+1] = CLAMP8(y >> 8);



	    r = src[loffset + k];
	    g = src[loffset + k+1];
	    b = src[loffset + k+2];
	    
	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v =  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u = -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[w+i+j] = CLAMP8(y >> 8);

	    r = src[loffset + k+3];
	    g = src[loffset + k+4];
	    b = src[loffset + k+5];
	    
	    k += 6;

	    y =  (FP(0.257) * r) + (FP(0.504) * g) + (FP(0.098) * b) + FP(16);
	    v +=  (FP(0.439) * r) - (FP(0.368) * g) - (FP(0.071) * b);
	    u += -(FP(0.148) * r) - (FP(0.291) * g) + (FP(0.439) * b);

	    dst[w+i+j+1] = CLAMP8(y >> 8);

	    dst[uoffset+ (i>>1) + (j>>2)] = (CLAMP8((u >> 9)+128));
	    dst[voffset+ (i>>1) + (j>>2)] = (CLAMP8((v >> 9)+128));
	}
    }
}


/* these seem to be pretty slow */

static void llconv_yvu2rgb(unsigned char* src, unsigned char* dst, int nbpixels)
{
    int r,g,b,y,v,u,i;
    for (i=0; i<nbpixels; i++){
	y = src[0];
	v = src[1];
	u = src[2];


	b = FP(1.164) * (y - 16)                         + FP(2.018) * (u - 128);
	g = FP(1.164) * (y - 16) - FP(0.813) * (v - 128) - FP(0.391) * (u - 128);
	r = FP(1.164) * (y - 16) + FP(1.596) * (v - 128);

	dst[0] = CLAMP(r>>8);
	dst[1] = CLAMP(g>>8);
	dst[2] = CLAMP(b>>8);

	src += 3;
	dst += 3;
    }
}



/* convert yvu to yuyv */
static void llconv_yvu2yuyv(unsigned char *src, unsigned char *dst, unsigned int nbpixels)
{
    unsigned int y1, y2, u, v, i;

    for (i = 0; i < nbpixels/2; i++){

	y1 = src[0];
	y2 = src[3];
	v = (src[1] + src[4]) >> 1;
	u = (src[2] + src[5]) >> 1;
	dst[0] = y1;
	dst[1] = u;
	dst[2] = y2;
	dst[3] = v;

	src += 6;
	dst += 4;

    }

}



/* convert yuvu packed 8 bit unsigned to yv12 planar 16bit signed */
static void llconv_yuyv_packed_u8s16(unsigned char* ucsource, short int *sidest, unsigned int w, unsigned int h)
{
    unsigned int i, j;
    unsigned int *source = (unsigned int *)ucsource;

    unsigned int *dest = (unsigned int *)sidest;
    unsigned int uoffset = (w*h)>>1;
    unsigned int voffset = (w*h + ((w*h) >> 2)) >> 1;

    for(j=0; j < (h*w)>>1; j +=(w)){
	for(i=0; i< (w>>1); i+=2){
	    unsigned int y,u,v;
	    unsigned int v00, v01, v10, v11;
	    v00 = source[i+j];
	    v01 = source[i+j+1];
	    v10 = source[i+j+(w>>1)];
	    v11 = source[i+j+(w>>1)+1];
	    
	    // save luma
	    dest[i+j]          = ((v00 & 0x00ff00ff) << 7);
	    dest[i+j+1]        = ((v01 & 0x00ff00ff) << 7);
	    dest[i+j+(w>>1)]   = ((v10 & 0x00ff00ff) << 7);
	    dest[i+j+(w>>1)+1] = ((v11 & 0x00ff00ff) << 7);

	    // compute chroma

	    // mask out luma & shift right
	    v00 = (v00 & 0xff00ff00)>>1;
	    v01 = (v01 & 0xff00ff00)>>1;
	    v10 = (v10 & 0xff00ff00)>>1;
	    v11 = (v11 & 0xff00ff00)>>1;
	    
	    // average 2 scan lines
	    v00 += v10;
	    v01 += v11;

	    // combine
	    v = (v01 << 16) | (v00 & 0x0000ffff);
	    u = (v01 & 0xffff0000) | (v00 >> 16);

	    // flip sign bits for u,v
	    u ^= 0x80008000;
	    v ^= 0x80008000;

	    // save chroma
	    dest[uoffset + (i>>1) + (j>>2)] = u;
	    dest[voffset + (i>>1) + (j>>2)] = v;
	}
    }


}

#define CONVERT(x,y) ((x) + ((y)<<16))

void pdp_llconv(void *src, int stype, void *dst, int dtype, int w, int h)
{
    int conversion = CONVERT(stype, dtype);
    void *tmpbuf;

    switch(CONVERT(stype, dtype)){

    case CONVERT( RIF_YVU__P411_U8, RIF_YVU__P411_S16 ):
	llconv_yvu_planar_u8s16((unsigned char*)src, (short int *)dst, w*h);
	break;

    case CONVERT( RIF_YUV__P411_U8, RIF_YVU__P411_S16 ):
	llconv_yuv_planar_u8s16((unsigned char*)src, (short int *)dst, w*h);
	break;

    case CONVERT( RIF_YUYV_P____U8, RIF_YVU__P411_S16 ):
	llconv_yuyv_packed_u8s16((unsigned char*)src, (short int *)dst, w, h);
	break;

    case CONVERT( RIF_RGB__P____U8, RIF_YVU__P411_U8 ):
	llconv_rgb2yvu_planar8sub((unsigned char*) src, (unsigned char*) dst, w, h);
	break;

    case CONVERT( RIF_RGB__P____U8, RIF_YVU__P411_S16 ):
	llconv_rgb2yvu_planar16sub((unsigned char*) src, (short int*) dst, w, h);
	break;

    case CONVERT( RIF_RGBA_P____U8, RIF_YVU__P411_S16 ):
	llconv_rgba2yvu_planar16sub((unsigned char*) src, (short int*) dst, w, h);
	break;

    case CONVERT( RIF_BGR__P____U8, RIF_YVU__P411_S16 ):
	llconv_bgr2yvu_planar16sub((unsigned char*) src, (short int*) dst, w, h);
	break;

    case CONVERT( RIF_BGRA_P____U8, RIF_YVU__P411_S16 ):
	llconv_bgra2yvu_planar16sub((unsigned char*) src, (short int*) dst, w, h);
	break;

    case CONVERT( RIF_YVU__P411_S16, RIF_RGB__P____U8 ):
	llconv_yvu16planar2rgbpacked((short int*) src, (unsigned char*) dst, w, h);
	break;

    case CONVERT( RIF_YVU__P411_S16, RIF_YVU__P411_U8 ):
	llconv_yvu_planar_s16u8((short int*)src, (unsigned char*)dst, w*h);
	break;

    case CONVERT( RIF_GREY______S16, RIF_GREY______U8 ):
	llconv_grey_s16u8((short int*)src, (unsigned char*)dst, w*h);
	break;
    default:
	pdp_post("pdp_llconv: WARNING: no conversion routine defined for (%d)->(%d)", stype, dtype);

    }

}


#ifdef __cplusplus
}
#endif

--- NEW FILE: pdp_llconv_portable.c ---

/*
 *   Pure Data Packet system implementation. : portable low level format conversion code
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */

#define CLAMP(x) (((x)<0) ? 0 : ((x>255)? 255 : (x)))
#define FP(x) ((int)(((float)(x)) * 256.0f))

void pixel_unpack_portable_u8s16_y(unsigned char *src ,short int *dst, unsigned int nbpixels)
{
    unsigned int i;
    for (i=0; i<nbpixels; i++) dst[i] = ((short int)(src[i])) << 7;
}

void pixel_unpack_portable_u8s16_uv(unsigned char *src ,short int *dst, unsigned int nbpixels)
{
    unsigned int i;
    for (i=0; i<nbpixels; i++) dst[i] = (((short int)(src[i])) << 8) ^ 0x8000;
}


void pixel_pack_portable_s16u8_y(short int *src, unsigned char *dst, unsigned int nbpixels)
{
    unsigned int i;
    for (i=0; i<nbpixels; i++) dst[i] = (unsigned char)(CLAMP(src[i]>>7));
}

void pixel_pack_portable_s16u8_uv(short int *src, unsigned char *dst, unsigned int nbpixels)
{
    unsigned int i;
    unsigned short *usrc = (unsigned short *)src;
    for (i=0; i<nbpixels; i++) dst[i] = ((usrc[i]^0x8000)>>8);
}


/* convert greyscale 8 bit unsigned to 16bit signed */
void llconv_grey_s16u8(short int *src, unsigned char *dst, unsigned int nbpixels)
{
    pixel_pack_portable_s16u8_y(src, dst, nbpixels);
}

/* convert yvu planar 411 16 bit signed to 8 bit unsigned */
void llconv_yvu_planar_s16u8(short int *src, unsigned char *dst, unsigned int nbpixels)
{
    pixel_pack_portable_s16u8_y(src, dst, nbpixels);
    pixel_pack_portable_s16u8_uv(src + nbpixels, dst + nbpixels, nbpixels>>1);

}


/* convert yvu planar 411 8 bit unsigned to yv12 planar 16bit signed */
void llconv_yvu_planar_u8s16(unsigned char* source, short int *dest, int nbpixels)
{
    pixel_unpack_portable_u8s16_y(source, dest, nbpixels);
    pixel_unpack_portable_u8s16_uv(&source[nbpixels], &dest[nbpixels], nbpixels>>1);
}

/* convert yuv planar 411 8 bit unsigned to yv12 planar 16bit signed */
void llconv_yuv_planar_u8s16(unsigned char* source, short int *dest, int nbpixels)
{
    pixel_unpack_portable_u8s16_y(source, dest, nbpixels);
    pixel_unpack_portable_u8s16_uv(&source[nbpixels], &dest[nbpixels + (nbpixels>>2)], nbpixels>>2);
    pixel_unpack_portable_u8s16_uv(&source[nbpixels + (nbpixels>>2)], &dest[nbpixels], nbpixels>>2);
}



--- NEW FILE: pdp_imageproc_portable.c ---
/*
 *   Pure Data Packet. portable image processing routines.
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */



#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "pdp_imageproc.h"

/* pdp memory alloc/dealloc prototype */
void *pdp_alloc(int size);
void pdp_dealloc(void *);


// utility stuff
inline static s32 float2fixed(float f)
{
    if (f > 1) f = 1;
    if (f < -1) f = -1;
    f *= 0x7fff;
    return (s32)f;
}



#define CLAMP16(x) (((x) > 0x7fff) ? 0x7fff : (((x) < -0x7fff) ? -0x7fff : (x)))

// add two images
void pdp_imageproc_add_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    int a, b;
    unsigned int i;
    for (i=0; i<width*height; i++){
	a = (int)image[i];
	b = (int)image2[i];
	image[i] = (s16)(CLAMP16(a+b));
    }
    
}

// mul two images
void pdp_imageproc_mul_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    int a, b;
    unsigned int i;
    for (i=0; i<width*height; i++){
	a = (int)image[i];
	b = (int)image2[i];
	image[i] = (s16)((a*b)>>15);
    }
    
}

// mix 2 images
void *pdp_imageproc_mix_new(void){return pdp_alloc(2*sizeof(s32));}
int  pdp_imageproc_mix_nb_stackwords(void) {return PDP_IMAGEPROC_NB_STACKWORDS(2*sizeof(s32));}
void pdp_imageproc_mix_delete(void *x) {pdp_dealloc (x);}
void pdp_imageproc_mix_setleftgain(void *x, float gain)
{
    s32 *d = (s32 *)x;
    d[0] = float2fixed(gain);
}
void pdp_imageproc_mix_setrightgain(void *x, float gain)
{
    s32 *d = (s32 *)x;
    d[1] = float2fixed(gain);
}
void pdp_imageproc_mix_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    s32 *d = (s32 *)x;
    u32 i;
    s32 a,b;

    for(i=0; i<width*height; i++){
	a = (s32)image[i];
	b = (s32)image2[i];
	a = (a*d[0] + b*d[1]) >> 15;
	image[i] = (s16)CLAMP16(a);
    }
	
}


// random mix 2 images
void *pdp_imageproc_randmix_new(void){return pdp_alloc(2*sizeof(s32));;}
int  pdp_imageproc_randmix_nb_stackwords(void) {return PDP_IMAGEPROC_NB_STACKWORDS(2*sizeof(s32));}
void pdp_imageproc_randmix_delete(void *x) {pdp_dealloc(x);}
void pdp_imageproc_randmix_setthreshold(void *x, float threshold)
{
    s32 *d = (s32 *)x;
    if (threshold > 1.0f) threshold = 1.0f;
    if (threshold < 0.0f) threshold = 0.0f;
    d[0] = float2fixed(threshold);
}
void pdp_imageproc_randmix_setseed(void *x, float seed)
{
    s32 *d = (s32 *)x;
    d[1] = float2fixed(seed);
}
void pdp_imageproc_randmix_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    s32 *d = (s32 *)x;
    u32 i;
    s16 r;
    srandom((u32)d[1]);


    for(i=0; i<width*height; i++){
	// get a random val between 0 and 0x7fff
	r = (s16)(random() & 0x7fff);
	if (r < d[0]) image[i] = image2[i];
    }
}


// 3x1 or 1x3 in place convolution
// orientation
void *pdp_imageproc_conv_new(void){return(pdp_alloc(6*sizeof(s32)));}
void pdp_imageproc_conv_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_conv_setmin1(void *x, float val)
{
    s32 *d = (s32 *)x;
    d[0] = float2fixed(val);
}
void pdp_imageproc_conv_setzero(void *x, float val)
{
    s32 *d = (s32 *)x;
    d[1] = float2fixed(val);
}
void pdp_imageproc_conv_setplus1(void *x, float val)
{
    s32 *d = (s32 *)x;
    d[2] = float2fixed(val);
}
void pdp_imageproc_conv_setbordercolor(void *x, float val)
{
    s32 *d = (s32 *)x;
    d[3] = float2fixed(val);
}
void pdp_imageproc_conv_setorientation(void *x, u32 val){((u32 *)x)[4] = val;}
void pdp_imageproc_conv_setnbpasses(void *x, u32 val){((u32 *)x)[5] = val;}

static inline void pdp_imageproc_conv_scanline(void *x, s16 *data, u32 count, s32 stride)
{
    s32 *d = (s32 *)x;
    s32 a,b,c,r;
    u32 i;

    a = d[3]; //border
    b = data[0];
    c = data[stride];

    for(i = 0; i < count-2; i++){
	r = a*d[0] + b*d[1] + c*d[2];
	a = data[0];
	b = data[stride];
	c = data[stride<<1];
	data[0] = (s16)CLAMP16(r>>15);
	data += stride;
    }
    r = a*d[0] + b*d[1] + c*d[2];
    a = data[0];
    b = data[stride];
    c = d[3]; //border
    data[0] = (s16)CLAMP16(r>>15);
    r = a*d[0] + b*d[1] + c*d[2];
    data[stride] = (s16)CLAMP16(r>>15);

}

void pdp_imageproc_conv_process(void *x, u32 width, u32 height, s16 *image)
{
    s32 *d = (s32 *)x;
    u32 i, j;
    u32 orientation = d[4];
    u32 nbp = d[5];
    if (orientation == PDP_IMAGEPROC_CONV_HORIZONTAL){
	for(i=0; i<width*height; i+=width)
	    for(j=0; j<nbp; j++)
		pdp_imageproc_conv_scanline(x, image+i, width, 1);

    }

    if (orientation == PDP_IMAGEPROC_CONV_VERTICAL){
	for(i=0; i<width; i++)
	    for(j=0; j<nbp; j++)
		pdp_imageproc_conv_scanline(x, image+i, height, width);

    }



	
}

// apply a gain to an image
void *pdp_imageproc_gain_new(void){return(pdp_alloc(2*sizeof(s32)));}
void pdp_imageproc_gain_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_gain_setgain(void *x, float gain)
{
    /* convert float to s16 + shift */
    s32 *d = (s32 *)x;
    s32 g;
    int i;
    float sign;
    s32 shift = 0;
    
    sign = (gain < 0) ? -1 : 1;
    gain *= sign;

    /* max shift = 16 */
    for(i=0; i<=16; i++){
	if (gain < 0x4000){
	    gain *= 2;
	    shift++;
	}
	else break;
    }

    gain *= sign;
    g = (s32) gain;

    //g = 0x4000;
    //shift = 14;

    d[0]=g;
    d[1]=shift;
}
void pdp_imageproc_gain_process(void *x, u32 width, u32 height, s16 *image)
{
    s32 *d = (s32 *)x;
    s32 a;
    u32 i;
    for (i=0; i<width*height; i++){
	a = (s32)image[i];
	image[i] = (s16)(CLAMP16((a * d[0]) >> d[1]));
    }
}

// colour rotation for 2 colour planes
void *pdp_imageproc_crot2d_new(void){return pdp_alloc(4*sizeof(s32));}
void pdp_imageproc_crot2d_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_crot2d_setmatrix(void *x, float *matrix)
{
    s32 *d = (s32 *)x;
    d[0] = float2fixed(matrix[0]);
    d[1] = float2fixed(matrix[1]);
    d[2] = float2fixed(matrix[2]);
    d[3] = float2fixed(matrix[3]);

}
void pdp_imageproc_crot2d_process(void *x, s16 *image, u32 width, u32 height)
{
    s32 *d = (s32 *)x;
    u32 i,j;
    s32 a1,a2,c1,c2;

    for(i=0, j=width*height; i<width*height; i++, j++){
	c1 = (s32)image[i];
	c2 = (s32)image[j];
	
	a1 = d[0] * c1;
	a2 = d[1] * c1;
	a1+= d[2] * c2;
	a2+= d[3] * c2;

	a1 >>= 15;
	a2 >>= 15;

	image[i] = (s16)CLAMP16(a1);
	image[j] = (s16)CLAMP16(a2);
    }
}

// biquad and biquad time
typedef struct
{
    s32 ma1;
    s32 ma2;
    s32 b0;
    s32 b1;
    s32 b2;

    s32 u0;
    s32 u1;

    s32 u0_save;
    s32 u1_save;

    u32 nbpasses;
    u32 direction;
} t_bq;
void *pdp_imageproc_bq_new(void){return pdp_alloc(sizeof(t_bq));}
void pdp_imageproc_bq_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_bq_setnbpasses(void *x, u32 i){((t_bq *)x)->nbpasses = i;}
void pdp_imageproc_bq_setdirection(void *x, u32 i){((t_bq *)x)->direction = i;}
void pdp_imageproc_bq_setcoef(void *x, float *coef) // a0,-a1,-a2,b0,b1,b2,u0,u1
{
    s32 *d = (s32 *)x;
    float ia0 = 1.0f / coef[0];

    /* all coefs are s1.14 fixed point */
    /* representing values -2 < x < 2  */
    /* so scale down before using the ordinary s0.15 float->fixed routine */

    ia0 *= 0.5f;

    // coef
    d[0] = float2fixed(ia0*coef[1]); // -a1
    d[1] = float2fixed(ia0*coef[2]); // -a2
    d[2] = float2fixed(ia0*coef[3]); // b0
    d[3] = float2fixed(ia0*coef[4]); // b1
    d[4] = float2fixed(ia0*coef[5]); // b2


    // state to reset too
    d[5] = float2fixed(coef[6]);
    d[6] = float2fixed(coef[7]);

}

#define A1 d[0]
#define A2 d[1]
#define B0 d[2]
#define B1 d[3]
#define B2 d[4]
/*
 	# DIRECT FORM II BIQUAD (from pixel_biquad_s16.s)
 	#
 	# y[k]  = b0 * x[k] + u1[k-1]
 	# u1[k] = b1 * x[k] + u2[k-1] - a1 * y[k]
	# u2[k] = b2 * x[k]           - a2 * y[k]
*/

/* remark A1 and A2 are already negated) */


static inline void pdp_imageproc_bq_scanline(void *x, s16 *data, u32 count, s32 stride)
{

    s32 *d = (s32 *)x;
    s32 u1,u2, xx, yy;

    u32 i;

    u1 = d[7];
    u2 = d[8];

    for(i = 0; i < count; i++){

	xx = (s32)data[0];

	yy = ((B0 * xx)>>14) + u1;
	u1 = ((B1 * xx)>>14) + u2 + ((A1 * yy)>>14);
	u2 = ((B2 * xx)>>14)      + ((A2 * yy)>>14);

	data[0] = (s16)CLAMP16(yy);

	data += stride;

    }

    d[7] = u1;
    d[8] = u2;

}

void pdp_imageproc_bqt_process(void *x, u32 width, u32 height, s16 *image, s16 *state1, s16 *state2)
{
    s32 *d = (s32 *)x;
    u32 i;
    s32 u1, u2, xx, yy;

    for (i=0; i<width*height; i++){

	xx = (s32)image[i];
	u1 = (s32)state1[i];
	u2 = (s32)state2[i];

	yy = ((B0 * xx)>>14) + u1;
	u1 = ((B1 * xx)>>14) + u2 + ((A1 * yy)>>14);
	u2 = ((B2 * xx)>>14)      + ((A2 * yy)>>14);

	image[i] = (s16)CLAMP16(yy);
	state1[i] = (s16)CLAMP16(u1);
	state2[i] = (s16)CLAMP16(u2);
    }
	
	
}

void pdp_imageproc_bq_process(void *x, u32 width, u32 height, s16 *data)
{
    s32 *d = (s32 *)x;
    
    u32 nbp = d[9];
    u32 direction = d[10];
    unsigned int i,j, offset;


    /* VERTICAL */
    offset = (height-1)*width;

    if ((direction & PDP_IMAGEPROC_BIQUAD_TOP2BOTTOM)
	&& (direction &  PDP_IMAGEPROC_BIQUAD_BOTTOM2TOP)){

	for(i=0; i<width; i++){
	    for (j=0; j<nbp; j++){
		pdp_imageproc_bq_scanline(x, data+i, height, width); //T->B
		pdp_imageproc_bq_scanline(x, data+offset+i, height, -width); //B->T
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_TOP2BOTTOM){
	for(i=0; i<width; i++){
	    for (j=0; j<nbp; j++){
		pdp_imageproc_bq_scanline(x, data+i, height, width); //T->B
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_BOTTOM2TOP){
	for(i=0; i<width; i++){
	    for (j=0; j<nbp; j++){
		pdp_imageproc_bq_scanline(x, data+offset+i, height, -width); //B->T
	    }
	}
    }

    /* HORIZONTAL */

    offset = width-1;
    if ((direction & PDP_IMAGEPROC_BIQUAD_LEFT2RIGHT)
	&& (direction & PDP_IMAGEPROC_BIQUAD_RIGHT2LEFT)){

	for(i=0; i<(width*height); i += width){
	    for (j=0; j<nbp; j++){
		pdp_imageproc_bq_scanline(x, data+i, width, 1); //L->R
		pdp_imageproc_bq_scanline(x, data+offset+i, width, -1); //R->L
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_LEFT2RIGHT){
	for(i=0; i<(width*height); i += width){
	    for (j=0; j<nbp; j++){
		pdp_imageproc_bq_scanline(x, data+i, width, 1); //L->R
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_RIGHT2LEFT){
	for(i=0; i<(width*height); i += width){
	    for (j=0; j<nbp; j++){
		pdp_imageproc_bq_scanline(x, data+offset+i, width, -1); //R->L

	    }
	}
    }

}

// produce a random image
// note: random number generator can be platform specific
// however, it should be seeded. (same seed produces the same result)
void *pdp_imageproc_random_new(void){return pdp_alloc(sizeof(s32));}
void pdp_imageproc_random_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_random_setseed(void *x, float seed)
{
    float *f = (float *)x;
    u32 *d = (u32 *)x;
    f[0] = seed;
    srandom(d[0]);
}

void pdp_imageproc_random_process(void *x, u32 width, u32 height, short int *image)
{
    s32 *d = (u32 *)x;
    u32 i;
    s32 r;
    srandom(d[0]);
    for (i=0; i<(width*height); i++) {
	r = random();
	image[i] = r;
    }
    d[0] = random();
}



/* resampling code */
// zoom + rotate

/* bilinear resampling core routine */
/* virtual coordinates are the lowest 16 bits in virt_x and virt_y*/
static inline s32 pdp_resample_bilin(s16 *image, s32 width, s32 height, s32 virt_x, s32 virt_y)
{

    s32 fp_x, fp_y, frac_x, frac_y, f, offset, r_1, r_2;

    //virt_x &= 0xffff;
    //virt_y &= 0xffff;

    fp_x = virt_x * (width - 1);
    fp_y = virt_y * (height - 1);

    frac_x = fp_x & (0xffff);
    frac_y = fp_y & (0xffff);

    offset = (fp_x >> 16) + (fp_y >> 16) * width;
    image += offset;

    f = 0x10000 - frac_x;

    r_1 = ((f * (s32)(image[0])  +  frac_x * (s32)(image[1])))>>16;

    image += width;

    r_2 = ((f * (s32)(image[0])  +  frac_x * (s32)(image[1])))>>16;

    f = 0x10000 - frac_y;

    return ((f * r_1 + frac_y * r_2)>>16);
    
}

typedef struct
{
    float centerx;
    float centery;
    float zoomx;
    float zoomy;
    float angle;
} t_affine_map;


void *pdp_imageproc_resample_affinemap_new(void)
{

    t_affine_map *a  = (t_affine_map *)pdp_alloc(sizeof(t_affine_map));
    a->centerx = 0.5;
    a->centery = 0.5;
    a->zoomx = 1.0;
    a->zoomy = 1.0;
    a->angle = 0.0f;
    return (void *)a;
}

void pdp_imageproc_resample_affinemap_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_resample_affinemap_setcenterx(void *x, float f){((t_affine_map *)x)->centerx = f;}
void pdp_imageproc_resample_affinemap_setcentery(void *x, float f){((t_affine_map *)x)->centery = f;}
void pdp_imageproc_resample_affinemap_setzoomx(void *x, float f){((t_affine_map *)x)->zoomx = f;}
void pdp_imageproc_resample_affinemap_setzoomy(void *x, float f){((t_affine_map *)x)->zoomy = f;}
void pdp_imageproc_resample_affinemap_setangle(void *x, float f){((t_affine_map *)x)->angle = f;}
void pdp_imageproc_resample_affinemap_process(void *x, u32 width, u32 height, s16 *srcimage, s16 *dstimage)
{
    t_affine_map *a = (t_affine_map *)x;
    double izx = 1.0f / (a->zoomx);
    double izy = 1.0f / (a->zoomy);
    double scale = (double)0xffffffff;
    double scalew = scale / ((double)(width - 1));
    double scaleh = scale / ((double)(height - 1));
    double cx = ((double)a->centerx) * ((double)(width - 1));
    double cy = ((double)a->centery) * ((double)(height - 1));
    double angle = a->angle * (-M_PI / 180.0);
    double c = cos(angle);
    double s = sin(angle);

    /* affine x, y mappings in screen coordinates */
    double mapx(double x, double y){return cx + izx * ( c * (x-cx) + s * (y-cy));}
    double mapy(double x, double y){return cy + izy * (-s * (x-cx) + c * (y-cy));}

    u32 colstate_x = (u32)(scalew * mapx(0,0));
    u32 colstate_y = (u32)(scaleh * mapy(0,0));
    u32 rowstate_x = colstate_x;
    u32 rowstate_y = colstate_y;

    u32 row_inc_x = (u32)(scalew * (mapx(1,0)-mapx(0,0)));
    u32 row_inc_y = (u32)(scaleh * (mapy(1,0)-mapy(0,0)));
    u32 col_inc_x = (u32)(scalew * (mapx(0,1)-mapx(0,0)));
    u32 col_inc_y = (u32)(scaleh * (mapy(0,1)-mapy(0,0)));

    u32 i,j;

    for (j=0; j<height; j++){
	for (i=0; i<width; i++){
	    *dstimage++ = pdp_resample_bilin(srcimage, width, height, rowstate_x>>16, rowstate_y>>16);
	    rowstate_x += row_inc_x;
	    rowstate_y += row_inc_y;
	}
	colstate_x += col_inc_x;
	colstate_y += col_inc_y;
	rowstate_x = colstate_x;
	rowstate_y = colstate_y;
    }

}





// polynomials




typedef struct
{
    u32 order;
    u32 nbpasses;
    s32 coefs[0];
} t_cheby;

void *pdp_imageproc_cheby_new(int order)
{
    t_cheby *z;
    int i;
    if (order < 2) order = 2;
    z = (t_cheby *)pdp_alloc(sizeof(t_cheby) + (order + 1) * sizeof(s32));
    z->order = order;
    z->coefs[0] = 0;
    z->coefs[1] = 0x7fff;
    for (i=2; i<=order; i++) z->coefs[i] = 0;
    return z;
}
void pdp_imageproc_cheby_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_cheby_setnbpasses(void *x, u32 n){((t_cheby *)x)->nbpasses = n;}
void pdp_imageproc_cheby_setcoef(void *x, u32 n, float f)
{

    t_cheby *z = (t_cheby *)x;
    if (n <= z->order){
	z->coefs[n] = (s32)(f * 32767.0f); // coefs are in s16.15 format
    }

}
void pdp_imageproc_cheby_process(void *x, u32 width, u32 height, s16 *image)
{

    t_cheby *z = (t_cheby *)x;
    u32 iterations = z->nbpasses;
    u32 i,j,k;
    s32 *c = z->coefs;
    for (j=0; j < (height*width); j++){
	s32 acc = (s32)image[j];
	for (i=0; i<iterations; i++){
	    s32 T2 = 0x7fff; /* 1 */
	    s32 T1 = acc;
	    s32 t;
	    s32 in = acc;
	    acc = c[0] + ((in*c[1])>>15);
	    for (k=2; k<=z->order; k++){
		t = ((T1*in)>>14) - T2; /* T_n = 2 x T_n-1 - T_n-2 */
		T2 = T1;
		T1 = t;
		acc += ((c[k] * t)>>15);
	    }
	}
	image[j] = (s16)(CLAMP16(acc));
    }
}

--- NEW FILE: pdp_llconv_mmx.c ---

/*
 *   Pure Data Packet system implementation. : wrapper for mmx low level format conversion code
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */


#include "pdp_mmx.h"



/* convert greyscale 8 bit unsigned to 16bit signed */
void llconv_grey_s16u8(short int *src, unsigned char *dst, unsigned int nbpixels)
{
    pixel_pack_s16u8_y(src, dst, nbpixels>>3);
}

/* convert yvu planar 411 16 bit signed to 8 bit unsigned */
void llconv_yvu_planar_s16u8(short int *src, unsigned char *dst, unsigned int nbpixels)
{
    pixel_pack_s16u8_y(src, dst, nbpixels>>3);
    pixel_pack_s16u8_uv(src + nbpixels, dst + nbpixels, nbpixels>>4);
}


/* convert yvu planar 411 8 bit unsigned to yv12 planar 16bit signed */
void llconv_yvu_planar_u8s16(unsigned char* source, short int *dest, int nbpixels)
{
    pixel_unpack_u8s16_y(source, dest, nbpixels>>3);
    pixel_unpack_u8s16_uv(&source[nbpixels], &dest[nbpixels], nbpixels>>4);
}

/* convert yuv planar 411 8 bit unsigned to yv12 planar 16bit signed */
void llconv_yuv_planar_u8s16(unsigned char* source, short int *dest, int nbpixels)
{
    pixel_unpack_u8s16_y(source, dest, nbpixels>>3);
    pixel_unpack_u8s16_uv(&source[nbpixels], &dest[nbpixels + (nbpixels>>2)], nbpixels>>5);
    pixel_unpack_u8s16_uv(&source[nbpixels + (nbpixels>>2)], &dest[nbpixels], nbpixels>>5);
}


--- NEW FILE: Makefile ---
include ../../Makefile.config

all: $(PDP_TARGET)

OBJECTS = pdp_llconv.o pdp_resample.o pdp_imageproc_common.o

OBJECTS_MMX = pdp_imageproc_mmx.o pdp_llconv_mmx.o
OBJECTS_PORTABLE =  pdp_imageproc_portable.o pdp_llconv_portable.o




linux_mmx: $(OBJECTS_MMX) $(OBJECTS)

linux: $(OBJECTS_PORTABLE) $(OBJECTS)

darwin: $(OBJECTS_PORTABLE) $(OBJECTS)

clean:
	rm -f *~
	rm -f *.o

--- NEW FILE: pdp_resample.c ---
/*
 *   Pure Data Packet system file. - image resampling routines
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */


#include <string.h>
#include "pdp_resample.h"


/*

efficient bilinear resampling ??
performance: how to eliminate divides? -> virtual coordinates 2^k x 2^k (conf. opengl)

i.e. 16 bit virtual coordinates: easy modular addressing

*/


/* code in this file should go out to be replaced by code in pdp_imageproc */

static s32 pdp_resample_bilin(s16 *image, s32 width, s32 height, s32 virt_x, s32 virt_y)
{

    s32 fp_x, fp_y, frac_x, frac_y, f, offset, r_1, r_2;

    virt_x &= 0xffff;
    virt_y &= 0xffff;

    fp_x = virt_x * (width - 1);
    fp_y = virt_y * (height - 1);

    frac_x = fp_x & (0xffff);
    frac_y = fp_y & (0xffff);

    offset = (fp_x >> 16) + (fp_y >> 16) * width;
    image += offset;

    f = 0x10000 - frac_x;

    r_1 = ((f * (s32)(image[0])  +  frac_x * (s32)(image[1])))>>16;

    image += width;

    r_2 = ((f * (s32)(image[0])  +  frac_x * (s32)(image[1])))>>16;

    f = 0x10000 - frac_y;

    return ((f * r_1 + frac_y * r_2)>>16);
    
}


void pdp_resample_scale_bilin(s16 *src_image, s16 *dst_image, s32 src_w, s32 src_h, s32 dst_w, s32 dst_h)
{
    s32 i,j;
    s32 virt_x=0;
    s32 virt_y=0; /* virtual coordinates in 30 bit */
    s32 scale_x = 0x40000000 / dst_w;
    s32 scale_y = 0x40000000 / dst_h;

    for (j=0; j<dst_h; j++){
	for (i=0; i<dst_w; i++){
	    *dst_image++ = pdp_resample_bilin(src_image, src_w, src_h, virt_x>>14, virt_y>>14);
	    virt_x += scale_x;
	}
	virt_x = 0;
	virt_y += scale_y;
    }

}

void pdp_resample_scale_nn(s16 *src_image, s16 *dst_image, s32 src_w, s32 src_h, s32 dst_w, s32 dst_h)
{
    s32 i,j;
    s32 x=0;
    s32 y=0;
    s32 frac_x=0;
    s32 frac_y=0;
    s32 scale_x = (src_w << 20 ) / dst_w;
    s32 scale_y = (src_h << 20 ) / dst_h;

    for (j=0; j<dst_h; j++){
	for (i=0; i<dst_w; i++){
	    *dst_image++ = src_image[x+y];
	    frac_x += scale_x;
	    x = frac_x >> 20;
	}
	x = 0;
	frac_x = 0;
	frac_y += scale_y;
	y = (frac_y >> 20) * src_w;
    }

}

/* USE pdp_resample_affinemap
void pdp_resample_zoom_tiled_bilin(s16 *src_image, s16 *dst_image, s32 w, s32 h, 
				   float zoom_x, float zoom_y, float center_x_relative, float center_y_relative)
{
    float izx = 1.0f / zoom_x;
    float izy = 1.0f / zoom_y;
    s32 scale_x = (s32)((float)0x100000 * izx / (float)w);
    s32 scale_y = (s32)((float)0x100000 * izy / (float)h);

    s32 top_virt_x = (s32)((1.0f - izx) * (float)0x100000 * center_x_relative);
    s32 top_virt_y = (s32)((1.0f - izy) * (float)0x100000 * center_y_relative);

    s32 virt_x = top_virt_x;
    s32 virt_y = top_virt_y; 

    s32 i,j;

    for (j=0; j<h; j++){
	for (i=0; i<w; i++){
	    *dst_image++ = pdp_resample_bilin(src_image, w, h, virt_x>>4, virt_y>>4);
	    virt_x += scale_x;
	}
	virt_x = top_virt_x;
	virt_y += scale_y;
    }

}
*/

void pdp_resample_halve(s16 *src_image, s16 *dst_image, s32 src_w, s32 src_h)
{

    int dst_x,dst_y;
    int src_x = 0;
    int src_y = 0;
    int dst_w = src_w >> 1;
    int dst_h = src_h >> 1;
    s32 tmp1,tmp2,tmp3,tmp4;

    //post("%x %x %d %d\n", src_image, dst_image, src_w, src_h);

    for(dst_y = 0; dst_y < dst_h * dst_w; dst_y += dst_w){
	for (dst_x = 0; dst_x < dst_w; dst_x++){

	    tmp1 = (s32)src_image[src_y + src_x];
	    tmp2 = (s32)src_image[src_y + src_x + 1];
	    tmp3 = (s32)src_image[src_y + src_x + src_w];
	    tmp4 = (s32)src_image[src_y + src_x + src_w + 1];

	    tmp1 += tmp2;
	    tmp3 += tmp4;

	    src_x += 2;

	    dst_image[dst_x+dst_y] = (s16)((tmp1 + tmp3)>>2);
	}
	src_y += src_w << 1;
	src_x = 0;
    }
}

void pdp_resample_double(s16 *src_image, s16 *dst_image, s32 src_w, s32 src_h)
{
    int src_x = 0;
    int src_y = 0;
    int dst = 0;
    int dst_w = src_w << 1;

    s16 tmp;

    for(src_y = 0; src_y < src_h * src_w; src_y += src_w){
	for (src_x = 0; src_x < src_w; src_x++){

	    tmp = *src_image++;
      	    dst = (src_y << 2) + (src_x << 1);
	    dst_image[dst]   = tmp;
	    dst_image[dst+1] = tmp;
	    dst+=dst_w;
	    dst_image[dst]   = tmp;
	    dst_image[dst+1] = tmp;
	}
    }
}

/* $$$TODO: finish this */
void pdp_resample_padcrop(s16 *src_image, s16 *dst_image, s32 src_w, s32 src_h, s32 dst_w, s32 dst_h)
{

    int shift_x = (dst_w - src_w) / 2;
    int shift_y = (dst_h - src_h) / 2;
}


--- NEW FILE: pdp_imageproc_common.c ---
/*
 *   Pure Data Packet. common image processing routines.
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */


/*
  This file contains common code for (portable) low level image processing objects
  pdp_imageproc_* methods
  The rest is int pdp_imageproc_<platform>.c

  There are also highlevel dispatcher methods that operate on packets:
  pdp_imageproc_dispatch_* methods

*/

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "pdp_imageproc.h"
#include "pdp_image.h"
#include "pdp_mem.h"
#include "pdp_packet.h"

#define CLAMP16(x) (((x) > 0x7fff) ? 0x7fff : (((x) < -0x7fff) ? -0x7fff : (x)))

u32 pdp_imageproc_legalwidth(int i)
{
    if (i>1024) return 1024;
    if (i>0) return  ((((i-1)>>3)+1)<<3);
    return 8;
    
}

u32 pdp_imageproc_legalheight(int i)
{
    if (i>1024) return 1024;
    if (i>0) return  ((((i-1)>>3)+1)<<3);
    return 8;
}
u32 pdp_imageproc_legalwidth_round_down(int i)
{
    if (i>1024) return 1024;
    if (i>8) return  ((i>>3)<<3);
    return 8;
    
}

u32 pdp_imageproc_legalheight_round_down(int i)
{
    if (i>1024) return 1024;
    if (i>8) return  ((i>>3)<<3);
    return 8;
}


/* check if two packets are allocated and of the same type */
bool pdp_packet_compat(int packet0, int packet1)
{

    t_pdp *header0 = pdp_packet_header(packet0);
    t_pdp *header1 = pdp_packet_header(packet1);
    if (!(header1)) return 0;
    if (!(header0)) return 0;
    if (header0->type != header1->type)	return 0;
    return 1;
}

/* some operations */

/* logic operators */

void pdp_imageproc_xor_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    u32 *plane = (u32 *)image;
    u32 *plane2 = (u32 *)image2;
    int count = (width * height) >> 1;
    int i;

    for (i=0; i<count; i++){
	plane[i] ^= plane2[i];
    }
}

void pdp_imageproc_and_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    u32 *plane = (u32 *)image;
    u32 *plane2 = (u32 *)image2;
    int count = (width * height) >> 1;
    int i;

    for (i=0; i<count; i++){
	plane[i] &= plane2[i];
    }
}

void pdp_imageproc_or_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    u32 *plane = (u32 *)image;
    u32 *plane2 = (u32 *)image2;
    int count = (width * height) >> 1;
    int i;

    for (i=0; i<count; i++){
	plane[i] |= plane2[i];
    }
}

void pdp_imageproc_not_process(void *x, u32 width, u32 height, s16 *image)
{
    u32 *plane = (u32 *)image;
    int count = (width * height) >> 1;
    int i;

    for (i=0; i<count; i++){
	plane[i] ^= 0xffffffff;
    }
}

void pdp_imageproc_mask_process(void *x, u32 width, u32 height, s16 *image)
{
    u32 mask = (u32)x;
    u32 *plane = (u32 *)image;
    int count = (width * height) >> 1;
    int i;

    mask = (mask & 0xffff) | (mask << 16);

    for (i=0; i<count; i++){
	plane[i] &= mask;
    }
}

// produce a plasma image
// note: random number generator can be platform specific
// however, it should be seeded. (same seed produces the same result)

typedef struct
{
    u32 seed;
    s32 scale;
} t_plasma;

static inline s16 _rand_s16(void)
{
  return (s16)(random()<<0);
}

static inline s16 _new_color(s32 one, s32 two, s32 scale)
{
  return CLAMP16((one >> 1) + (two >> 1) + ((scale * _rand_s16()) >> 16));
  //return (one >> 1) + (two >> 1);
}

void *pdp_imageproc_plasma_new(void){return pdp_alloc(sizeof(t_plasma));}
void pdp_imageproc_plasma_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_plasma_setseed(void *x, float seed)
{
    *((float *)x) = seed;
}
void pdp_imageproc_plasma_setturbulence(void *x, float f)
{
    ((t_plasma *)x)->scale = CLAMP16(f * ((float)0x7fff));
}

static void _plasma_subdiv(u32 w, u32 h, u32 s, s16 *image, int calc_left, int calc_top, s32 scale)
{
  int w0 = ((w-1)>>1);  // width of left segments
  int h0 = ((h-1)>>1);  // heigth of top segments
  int w1 = w - w0;
  int h1 = h - h0;

  /* conditions: w0 <= w1, h0 <= h1 */

  /* original coordinates */
  int topleft = 0;
  int topright = w-1;
  int bottomleft = s * (h-1);
  int bottomright = bottomleft + topright;

  /* new subdivision coordinates */
  int top = w0;
  int left = s * h0;
  int bottom = bottomleft + w0;
  int right = topright + left;
  int center = left + top;

  if (w0 && h0){ /* left-right and top-bottom subdivide */

    /* calculate corner pixel colours */
    if (calc_top)  image[top]    = _new_color(image[topleft], image[topright], scale);
    if (calc_left) image[left]   = _new_color(image[topleft], image[bottomleft], scale);
    image[right]  = _new_color(image[topright], image[bottomright], scale);
    image[bottom] = _new_color(image[bottomleft], image[bottomright], scale);
    image[center] = (_new_color(image[top], image[bottom], scale) >> 1)
	+(_new_color(image[left], image[right], scale) >> 1);


    /* subdivide (with overlap) */
    _plasma_subdiv(w0+1, h0+1, s, &image[topleft], 1, 1, scale);
    _plasma_subdiv(w1, h0+1, s, &image[top], 0, 1, scale);
    _plasma_subdiv(w0+1, h1, s, &image[left], 1, 0, scale);
    _plasma_subdiv(w1, h1, s, &image[center], 0, 0, scale);
    
  }
  

  else if(h0) { /* top-bottom subdivide */

      //post("h:%d", h);

    /* calculate corner pixel colours */
    if(calc_left) image[left]   = _new_color(image[topleft], image[bottomleft], scale);
    image[right]  = _new_color(image[topright], image[bottomright], scale);

    /* subdivide (without overlap) */
    _plasma_subdiv(w, h0+1, s, &image[topleft], 1, 0, scale);
    _plasma_subdiv(w, h1, s, &image[left], 1, 0, scale);
    
  }
  
  else if (w0){ /* left-right subdivide */

    /* calculate corner pixel colours */
    if (calc_top) image[top]    = _new_color(image[topleft], image[topright], scale);
    image[bottom] = _new_color(image[bottomleft], image[bottomright],scale);

    /* subdivide with overlap */
    _plasma_subdiv(w0+1, h, s, &image[topleft], 0, 1, scale);
    _plasma_subdiv(w1, h, s, &image[top], 0, 1, scale);

  }

}

void pdp_imageproc_plasma_process(void *x, u32 width, u32 height, s16 *image)
{
    s32 scale = (((t_plasma *)x)->scale);
    srandom (((t_plasma *)x)->seed);

    /* set initial border colours */
    image[0]                  = _rand_s16();
    image[width-1]            = _rand_s16();
    image[width * (height-1)] = _rand_s16();
    image[width * height - 1] = _rand_s16();

    /* subdivide */
    _plasma_subdiv(width, height, width, image, 1, 1, scale);

    ((t_plasma *)x)->seed = random();

}
  


void pdp_imageproc_zero_process(void *x, u32 width, u32 height, s16 *image)
{
    int bytesize = (width * height) << 1;
    memset(image, 0, bytesize);
}

void pdp_imageproc_constant_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    u32 value = (u32)x;
    u32 *plane = (u32 *)image;
    int wordsize = (width * height) >> 1;
    value = (value & 0xffff) | (value << 16);
    for (i=0; i<wordsize; i++){
	plane[i] = value;
    }
}


/* other stateless operators */

/* some 2x16bit vector ops */

/* some bit shuffling to ensure 32 bit accesses 
   get the sign bit extended as a mask: - : 0xffff +: 0x0000 */
static inline u32 _sign(s32 invec)
{
    s32 mask_top = invec;
    s32 mask_bot = invec;

    mask_top &= 0x80000000; /* isolate top sign bit */
    mask_bot <<= 16;        /* shift bottom word to top word */
    mask_bot &= 0x80000000; /* isolate bottom sign bit */
    mask_top >>= 15;        /* shift sign bit into top word */
    mask_bot >>= 15;
    mask_bot = (s32)(((u32)mask_bot) >> 16);  /* shift top word into bottom word */
    return mask_top |mask_bot;    
}

/* clear the least significant bit of the top word
   to ensure a decoupled vector add */
static inline void _decouple(s32 *invec)
{
    *invec &= 0xfffeffff;
}

void pdp_imageproc_abs_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    s32 *wimage = (s32 *)image;
    int wsize = (width * height) >> 1;
    for (i=0; i<wsize; i++){
	/* this computes  c = (c >= 0) ? (c) : (~c) */
	/* not is used instead of neg to prevent overflow on 0x8000 */
	/* this maps both 0 and -1 to 0 */

	wimage[i] ^= _sign(wimage[i]);
	
    }
}

void pdp_imageproc_zthresh_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    s32 *wimage = (s32 *)image;
    int wsize = (width * height) >> 1;
    for (i=0; i<wsize; i++){
	/* this computes  c = (c >= 0) ? (c) : (0) */
	wimage[i] &= ~_sign(wimage[i]);
    }
}

/* hard thresholding: x contains a positive unsigned short int */
void pdp_imageproc_hardthresh_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    s32 thresh = (s32)x;
    s32 sign1, isign2, a;
    s32 *wimage = (s32 *)image;
    int wsize = (width * height) >> 1;
    thresh |= (thresh << 16);
    for (i=0; i<wsize; i++){
	a = wimage[i];
	sign1 = _sign(a);   
	a ^= sign1;           /* take abs */
	_decouple(&a);
	a -= thresh;          /* subtract threshold */
	isign2 = ~ _sign(a);
	a &= isign2;          /* zero thresh */
	_decouple(&a);
	a += thresh & isign2; /* add threshold (if not zero thresholded)*/
	a ^= sign1;
	wimage[i] = a;
    }
}

/* soft thresholding: x contains a positive unsigned short int */
void pdp_imageproc_softthresh_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    s32 thresh = (s32)x;
    s32 sign1, sign2, a;
    s32 *wimage = (s32 *)image;
    int wsize = (width * height) >> 1;
    thresh |= thresh << 16;
    for (i=0; i<wsize; i++){
	a = wimage[i];
	sign1 = _sign(a);   
	a ^= sign1;           /* take abs */
	_decouple(&a);
	a -= thresh;          /* subtract threshold */
	sign2 = _sign(a);
	a &= ~ sign2;         /* zero thresh */
	_decouple(&a);
	//a += thresh;        /* add threshold */
	a ^= sign1;
	wimage[i] = a;
	
    }

}


/* turns an image into a positive andmask */
void pdp_imageproc_ispositive_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    s32 *wimage = (s32 *)image;
    int wsize = (width * height) >> 1;
    for (i=0; i<wsize; i++){
	wimage[i] = ~_sign(wimage[i]);
    }

}

/* get sign */
void pdp_imageproc_sign_process(void *x, u32 width, u32 height, s16 *image)
{
    int i;
    s32 *wimage = (s32 *)image;
    int wsize = (width * height) >> 1;
    for (i=0; i<wsize; i++){
	wimage[i] = _sign(wimage[i]) ^ 0x7fff7fff;
    }

}

/* flip left <-> right */
void pdp_imageproc_flip_lr_process(void *dummy, u32 width, u32 height, s16 *image)
{
    u32 y;
    s16 tmp, *l, *r;
    for (y=0; y<height; y++){
	l = image;
	r = image + width - 1;
	while (l < r){
	    tmp = *l;
	    *l = *r;
	    *r = tmp;
	    l++;
	    r--;
	}
	image += width;
    }

}

void pdp_llconv_flip_top_bottom(s16 *data, int width, int height, int pixelsize);

void pdp_imageproc_flip_tb_process(void *dummy, u32 width, u32 height, s16 *image)
{
    pdp_llconv_flip_top_bottom(image, width, height, 2);
}


/* image processing dispatcher methods  */
/* if the first packet contains a nonzero channel mask, it will be used instead
   of the one supplied as argument to the dispatcher functions.
   the packet's channel mask will be reset to 0 */

void pdp_imageproc_dispatch_1buf(void (*process_routine)(void*, u32, u32, s16*), void *x, u32 chanmask, int packet0)
{
    t_pdp *header0;
    t_image *image0;
    s16  *idata0;
    unsigned int w,h,d,plane_size,mask;
 
    /* if packet is not a valid image return without doing anything */
    if (!(pdp_packet_image_isvalid(packet0))) return;

    header0 = pdp_packet_header(packet0);
    image0 = pdp_packet_image_info(packet0);
    idata0   = pdp_packet_data  (packet0);

    w = image0->width;
    h = image0->height;
    d = image0->depth;
    plane_size = w*h;

    if (image0->chanmask) chanmask = image0->chanmask;
    image0->chanmask = 0;


    switch(image0->encoding){
    case PDP_IMAGE_GREY:
	if (chanmask & 1) (*process_routine)(x, w, h, idata0);
	break;
    case PDP_IMAGE_YV12:
	if (chanmask & 1) (*process_routine)(x, w, h, idata0);
	idata0 += plane_size;
	plane_size >>= 2;
	w >>= 1;
	h >>= 1;
	if (chanmask & 2) (*process_routine)(x, w, h, idata0);
	idata0 += plane_size;
	if (chanmask & 4) (*process_routine)(x, w, h, idata0);
	break;
    case PDP_IMAGE_MCHP:
	mask = 1;
	while (d--){
	    if (chanmask & mask) (*process_routine)(x, w, h, idata0);
	    idata0 += plane_size;
	    mask <<= 1;
	}
	break;
    default:
	break;
    }
}


void pdp_imageproc_dispatch_2buf(void (*process_routine)(void*, u32, u32, s16*, s16 *), void *x, u32 chanmask, int packet0, int packet1)
{
    t_pdp *header0;
    t_image *image0;
    s16  *idata0, *idata1;
    unsigned int w,h,d,plane_size,mask;
 
    /* if packets are not compatible images, return without doing anything */
    if (!(pdp_packet_image_compat(packet0, packet1))) return;

    header0 = pdp_packet_header(packet0);
    image0 = pdp_packet_image_info(packet0);
    idata0   = pdp_packet_data  (packet0);
    idata1   = pdp_packet_data  (packet1);

    w = image0->width;
    h = image0->height;
    d = image0->depth;
    plane_size = w*h;

    if (image0->chanmask) chanmask = image0->chanmask;
    image0->chanmask = 0;

    switch(image0->encoding){
    case PDP_IMAGE_GREY:
	if (chanmask & 1) (*process_routine)(x, w, h, idata0, idata1);
	break;
    case PDP_IMAGE_YV12:
	if (chanmask & 1) (*process_routine)(x, w, h, idata0, idata1);
	idata0 += plane_size;
	idata1 += plane_size;
	plane_size >>= 2;
	w >>= 1;
	h >>= 1;
	if (chanmask & 2) (*process_routine)(x, w, h, idata0, idata1);
	idata0 += plane_size;
	idata1 += plane_size;
	if (chanmask & 4) (*process_routine)(x, w, h, idata0, idata1);
	break;
    case PDP_IMAGE_MCHP:
	mask = 1;
	while (d--){
	    if (chanmask & mask) (*process_routine)(x, w, h, idata0, idata1);
	    idata0 += plane_size;
	    idata1 += plane_size;
	    mask <<= 1;
	}
	break;
    default:
	break;
    }
}
void pdp_imageproc_dispatch_3buf(void (*process_routine)(void*, u32, u32, s16*, s16 *, s16 *), void *x, u32 chanmask, int packet0, int packet1, int packet2)
{
    t_pdp *header0;
    t_image *image0;
    s16  *idata0, *idata1, *idata2;
    unsigned int w,h,d,plane_size, mask;
 
    /* if packets are not compatible images, return without doing anything */
    if (!((pdp_packet_image_compat(packet0, packet1))
	  &&(pdp_packet_image_compat(packet0, packet1)))) return;

    header0 = pdp_packet_header(packet0);
    image0 = pdp_packet_image_info(packet0);
    idata0   = pdp_packet_data  (packet0);
    idata1   = pdp_packet_data  (packet1);
    idata2   = pdp_packet_data  (packet2);

    w = image0->width;
    h = image0->height;
    d = image0->depth;
    plane_size = w*h;

    if (image0->chanmask) chanmask = image0->chanmask;
    image0->chanmask = 0;

    switch(image0->encoding){
    case PDP_IMAGE_GREY:
	if (chanmask & 1)(*process_routine)(x, w, h, idata0, idata1, idata2);
	break;
    case PDP_IMAGE_YV12:
	if (chanmask & 1)(*process_routine)(x, w, h, idata0, idata1, idata2);
	idata0 += plane_size;
	idata1 += plane_size;
	idata2 += plane_size;
	plane_size >>= 2;
	w >>= 1;
	h >>= 1;
	if (chanmask & 2)(*process_routine)(x, w, h, idata0, idata1, idata2);
	idata0 += plane_size;
	idata1 += plane_size;
	idata2 += plane_size;
	if (chanmask & 4)(*process_routine)(x, w, h, idata0, idata1, idata2);
	break;
    case PDP_IMAGE_MCHP:
	mask = 1;
	while (d--){
	    if (chanmask & mask) (*process_routine)(x, w, h, idata0, idata1, idata2);
	    idata0 += plane_size;
	    idata1 += plane_size;
	    idata2 += plane_size;
	    mask <<= 1;
	}
	break;
    default:
	break;
    }
}

--- NEW FILE: pdp_imageproc_mmx.c ---
/*
 *   Pure Data Packet. c wrapper for mmx image processing routines.
 *   Copyright (c) by Tom Schouten <pdp at zzz.kotnet.org>
 *
 *   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., 675 Mass Ave, Cambridge, MA 02139, USA.
 *
 */


/* this is a c wrapper around platform specific (mmx) code */
#include <stdlib.h>
#include <math.h>
#include "pdp_mmx.h"
#include "pdp_imageproc.h"

/* pdp memory alloc/dealloc prototype */
void *pdp_alloc(int size);
void pdp_dealloc(void *);

// utility stuff
inline static s16 float2fixed(float f)
{
    if (f > 1) f = 1;
    if (f < -1) f = -1;
    f *= 0x7fff;
    return (s16)f;
}

inline static void setvec(s16 *v, float f)
{
    s16 a = float2fixed(f);
    v[0] = a;
    v[1] = a;
    v[2] = a;
    v[3] = a;
}



// add two images
void pdp_imageproc_add_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    unsigned int totalnbpixels = width * height;
    pixel_add_s16(image, image2, totalnbpixels>>2);
}

// mul two images
void pdp_imageproc_mul_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    unsigned int totalnbpixels = width * height;
    pixel_mul_s16(image, image2, totalnbpixels>>2);
}


// mix 2 images
#define MIX_SIZE 8*sizeof(s16)
void *pdp_imageproc_mix_new(void){return pdp_alloc(MIX_SIZE);}
int pdp_imageproc_mix_nb_stackwords(void){return PDP_IMAGEPROC_NB_STACKWORDS(MIX_SIZE);}
void pdp_imageproc_mix_delete(void *x) {pdp_dealloc (x);}
void pdp_imageproc_mix_setleftgain(void *x, float gain){setvec((s16 *)x, gain);}
void pdp_imageproc_mix_setrightgain(void *x, float gain){setvec((s16 *)x + 4, gain);}
void pdp_imageproc_mix_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    s16 *d = (s16 *)x;
    unsigned int totalnbpixels = width * height;
    pixel_mix_s16(image, image2, totalnbpixels>>2, d, d+4);
}


// random mix 2 images
#define RANDMIX_SIZE 8*sizeof(s16)
void *pdp_imageproc_randmix_new(void){return pdp_alloc(RANDMIX_SIZE);}
int  pdp_imageproc_randmix_nb_stackwords(void) {return  PDP_IMAGEPROC_NB_STACKWORDS(RANDMIX_SIZE);}
void pdp_imageproc_randmix_delete(void *x) {pdp_dealloc (x);}
void pdp_imageproc_randmix_setthreshold(void *x, float threshold){setvec((s16 *)x, 2*threshold-1);}
void pdp_imageproc_randmix_setseed(void *x, float seed)
{
    s16 *d = (s16 *)x;
    srandom((u32)seed);
    d[4] = (s16)random();
    d[5] = (s16)random();
    d[6] = (s16)random();
    d[7] = (s16)random();
    
}
void pdp_imageproc_randmix_process(void *x, u32 width, u32 height, s16 *image, s16 *image2)
{
    s16 *d = (s16 *)x;
    unsigned int totalnbpixels = width * height;
    pixel_randmix_s16(image, image2, totalnbpixels>>2, d+4, d);
}


// 3x1 or 1x3 in place convolution
// orientation
typedef struct
{
    s16 min1[4];
    s16 zero[4];
    s16 plus1[4];
    s16 border[4];
    u32 orientation;
    u32 nbpasses;
} t_conv;
void *pdp_imageproc_conv_new(void){return(pdp_alloc(sizeof(t_conv)));}
void pdp_imageproc_conv_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_conv_setmin1(void *x, float val){setvec(((t_conv *)x)->min1, val);}
void pdp_imageproc_conv_setzero(void *x, float val){setvec(((t_conv *)x)->zero, val);}
void pdp_imageproc_conv_setplus1(void *x, float val){setvec(((t_conv *)x)->plus1, val);}
void pdp_imageproc_conv_setbordercolor(void *x, float val){setvec(((t_conv *)x)->border, val);}
void pdp_imageproc_conv_setorientation(void *x, u32 val){((t_conv *)x)->orientation = val;}
void pdp_imageproc_conv_setnbpasses(void *x, u32 val){((t_conv *)x)->nbpasses = val;}
void pdp_imageproc_conv_process(void *x, u32 width, u32 height, s16 *image)

{
    t_conv *d = (t_conv *)x;
    
    u32 orientation = d->orientation;
    u32 nbp = d->nbpasses;
    u32 i,j;

    if (orientation == PDP_IMAGEPROC_CONV_HORIZONTAL)
    {
	for(i=0; i<width*height; i+=width)
	    for (j=0; j<nbp; j++)
		pixel_conv_hor_s16(image+i, width>>2, d->border, d->min1);
    }

    else
    {
	for (j=0; j<nbp; j++)
	    for(i=0; i<width; i +=4) pixel_conv_ver_s16(image+i,  height, width, d->border, d->min1);
    }

	
	
}

// apply a gain to an image
void *pdp_imageproc_gain_new(void){return(pdp_alloc(8*sizeof(s16)));}
void pdp_imageproc_gain_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_gain_setgain(void *x, float gain)
{
    /* convert float to s16 + shift */
    s16 *d = (s16 *)x;
    s16 g;
    int i;
    float sign;
    int shift = 0;
    
    sign = (gain < 0) ? -1 : 1;
    gain *= sign;

    /* max shift = 16 */
    for(i=0; i<=16; i++){
	if (gain < 0x4000){
	    gain *= 2;
	    shift++;
	}
	else break;
    }

    gain *= sign;
    g = (s16) gain;

    //g = 0x4000;
    //shift = 14;

    d[0]=g;
    d[1]=g;
    d[2]=g;
    d[3]=g;
    d[4]=(s16)shift;
    d[5]=0;
    d[6]=0;
    d[7]=0;
}
void pdp_imageproc_gain_process(void *x, u32 width, u32 height, s16 *image)
{
    s16 *d = (s16 *)x;
    pixel_gain_s16(image, (width*height)>>2, d, (u64 *)(d+4));
}

// colour rotation for 2 colour planes
void *pdp_imageproc_crot2d_new(void){return pdp_alloc(16*sizeof(s16));}
void pdp_imageproc_crot2d_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_crot2d_setmatrix(void *x, float *matrix)
{
    s16 *d = (s16 *)x;
    setvec(d, matrix[0]);
    setvec(d+4, matrix[1]);
    setvec(d+8, matrix[2]);
    setvec(d+12, matrix[3]);
}
void pdp_imageproc_crot2d_process(void *x, s16 *image, u32 width, u32 height)
{
    s16 *d = (s16 *)x;
    pixel_crot2d_s16(image, width*height >> 2, d);
}

// biquad and biquad time
typedef struct
{
    s16 ma1[4];
    s16 ma2[4];
    s16 b0[4];
    s16 b1[4];
    s16 b2[4];
    s16 u0[4];
    s16 u1[4];
    s16 u0_save[4];
    s16 u1_save[4];
    u32 nbpasses;
    u32 direction;
} t_bq;

void *pdp_imageproc_bq_new(void){return pdp_alloc(sizeof(t_bq));}
void pdp_imageproc_bq_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_bq_setcoef(void *x, float *coef) // a0,-a1,-a2,b0,b1,b2,u0,u1
{
    t_bq *d = (t_bq *)x;
    float ia0 = 1.0f / coef[0];

    /* all coefs are s1.14 fixed point */
    /* representing values -2 < x < 2  */
    /* so scale down before using the ordinary s0.15 float->fixed routine */

    ia0 *= 0.5f;

    // coef
    setvec(d->ma1, ia0*coef[1]);
    setvec(d->ma2, ia0*coef[2]);
    setvec(d->b0, ia0*coef[3]);
    setvec(d->b1, ia0*coef[4]);
    setvec(d->b2, ia0*coef[5]);

    // state to reset too
    setvec(d->u0_save, coef[6]);
    setvec(d->u1_save, coef[7]);

}
void pdp_imageproc_bq_setnbpasses(void *x, u32 nbpasses){((t_bq *)x)->nbpasses = nbpasses;}
void pdp_imageproc_bq_setdirection(void *x, u32 direction){((t_bq *)x)->direction = direction;}
void pdp_imageproc_bq_process(void *x, u32 width, u32 height, s16* image);


void pdp_imageproc_bqt_process(void *x, u32 width, u32 height, s16 *image, s16 *state0, s16 *state1)
{
    s16 *d = (s16 *)x;
    pixel_biquad_time_s16(image, state0, state1, d, (width*height)>>2);
}

void pdp_imageproc_bq_process(void *x, u32 width, u32 height, s16 *image)
{
    t_bq *d = (t_bq *)x;
    s16 *c = d->ma1; /* coefs */
    s16 *s = d->u0;  /* state */
    u32 direction = d->direction;
    u32 nbp = d->nbpasses;
    unsigned int i,j;



    /* VERTICAL */

    if ((direction & PDP_IMAGEPROC_BIQUAD_TOP2BOTTOM)
	&& (direction &  PDP_IMAGEPROC_BIQUAD_BOTTOM2TOP)){

	for(i=0; i<width; i +=4){
	    for (j=0; j<nbp; j++){
		pixel_biquad_vertb_s16(image+i,    height>>2, width, c, s);
		pixel_biquad_verbt_s16(image+i,    height>>2, width, c, s);
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_TOP2BOTTOM){
	for(i=0; i<width; i +=4){
	    for (j=0; j<nbp; j++){
		pixel_biquad_vertb_s16(image+i,    height>>2, width, c, s);
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_BOTTOM2TOP){
	for(i=0; i<width; i +=4){
	    for (j=0; j<nbp; j++){
		pixel_biquad_verbt_s16(image+i,    height>>2, width, c, s);
	    }
	}
    }

    /* HORIZONTAL */

    if ((direction & PDP_IMAGEPROC_BIQUAD_LEFT2RIGHT)
	&& (direction & PDP_IMAGEPROC_BIQUAD_RIGHT2LEFT)){

	for(i=0; i<(width*height); i +=(width<<2)){
	    for (j=0; j<nbp; j++){
		pixel_biquad_horlr_s16(image+i,    width>>2, width, c, s);
		pixel_biquad_horrl_s16(image+i,    width>>2, width, c, s);
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_LEFT2RIGHT){
	for(i=0; i<(width*height); i +=(width<<2)){
	    for (j=0; j<nbp; j++){
		pixel_biquad_horlr_s16(image+i,    width>>2, width, c, s);
	    }
	}
    }

    else if (direction & PDP_IMAGEPROC_BIQUAD_RIGHT2LEFT){
	for(i=0; i<(width*height); i +=(width<<2)){
	    for (j=0; j<nbp; j++){
		pixel_biquad_horrl_s16(image+i,    width>>2, width, c, s);
	    }
	}
    }

}

// produce a random image
// note: random number generator can be platform specific
// however, it should be seeded. (same seed produces the same result)
void *pdp_imageproc_random_new(void){return pdp_alloc(4*sizeof(s16));}
void pdp_imageproc_random_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_random_setseed(void *x, float seed)
{
    s16 *d = (s16 *)x;
    srandom((u32)seed);
    d[0] = (s16)random();
    d[1] = (s16)random();
    d[2] = (s16)random();
    d[3] = (s16)random();
    
}
void pdp_imageproc_random_process(void *x, u32 width, u32 height, s16 *image)
{
    s16 *d = (s16 *)x;
    unsigned int totalnbpixels = width * height;
    pixel_rand_s16(image, totalnbpixels>>2, d);
}


/* resampling stuff
   this is quite a zoo of data structures
   the major point is this: the resampler mmx code is shared for all resampling code
   it uses data specified in t_resample_cbrd (Cooked Bilinear Resampler Data)

   then the there are several feeder algorithms. one is the linear mapper. it's
   data is specified in t_resample_clrd (Cooked Linear Remapper Data)

   for each feeder algorithm, there are several high level algorithms. like zoom,
   rotate, ... 
*/

typedef struct
{
    u32 lineoffset;
    s16 *image;
    u32 width;
    u32 height;
    
} t_resample_id; // Image Data

/* initialize image meta data (dimensions + location) */
static void pdp_imageproc_resample_init_id(t_resample_id *x, u32 offset, s16* image, u32 w, u32 h)
{
    x->lineoffset = offset;
    x->image = image;
    x->width = w;
    x->height = h;
}

// mmx resampling source image resampling data + coefs
typedef struct
{
    // vector data for resampling routine (resampling computation)
    u8  reserved[0x60];  //internal data
    s16 *address[2];     //64 bit splatted offset address
    s16 twowidthm1[4];   //64 bit splatted 2*(width-1)
    s16 twoheightm1[4];  //64 bit splatted 2*(height-1)
    s16 lineoffset[4];   //64 bit splatted line offset in pixels

} t_resample_cid; // Cooked Image Data

/* convert image meta data into a cooked format used by the resampler routine */
static void pdp_imageproc_resample_init_cid(t_resample_cid *r, t_resample_id *i)
{
    u32 twowm1 = (i->width-1)<<1;
    u32 twohm1 = (i->height-1)<<1;
    r->address[0] = i->image;
    r->address[1] = i->image;
    r->twowidthm1[0] = twowm1;
    r->twowidthm1[1] = twowm1;
    r->twowidthm1[2] = twowm1;
    r->twowidthm1[3] = twowm1;
    r->twoheightm1[0] = twohm1;
    r->twoheightm1[1] = twohm1;
    r->twoheightm1[2] = twohm1;
    r->twoheightm1[3] = twohm1;
    r->lineoffset[0] = i->lineoffset;
    r->lineoffset[1] = i->lineoffset;
    r->lineoffset[2] = i->lineoffset;
    r->lineoffset[3] = i->lineoffset;
}

// linear mapping data struct (zoom, scale, rotate, shear, ...)
typedef struct
{
    s32 rowstatex[2]; // row state x coord
    s32 rowstatey[2]; // row state y coord
    s32 colstatex[2]; // column state x coord
    s32 colstatey[2]; // column state y coord
    s32 rowincx[2];   // row inc vector x coord
    s32 rowincy[2];   // row inc vector y coord
    s32 colincx[2];   // column inc vector x coord
    s32 colincy[2];   // column inc vector y coord
} t_resample_clmd; // Cooked Linear Mapping Data

/* convert incremental linear remapping vectors to internal cooked format */
static void pdp_imageproc_resample_cookedlinmap_init(t_resample_clmd *l, s32 sx, s32 sy, s32 rix, s32 riy, s32 cix, s32 ciy)
{
    l->colstatex[0] = l->rowstatex[0] = sx;
    l->colstatex[1] = l->rowstatex[1] = sx + rix;
    l->colstatey[0] = l->rowstatey[0] = sy;
    l->colstatey[1] = l->rowstatey[1] = sy + riy;
    l->rowincx[0] = rix << 1;
    l->rowincx[1] = rix << 1;
    l->rowincy[0] = riy << 1;
    l->rowincy[1] = riy << 1;
    l->colincx[0] = cix;
    l->colincx[1] = cix;
    l->colincy[0] = ciy;
    l->colincy[1] = ciy;
}


/* this struct contains all the data necessary for
   bilin interpolation from src -> dst image
   (src can be == dst) */
typedef struct
{
    t_resample_cid csrc;     //cooked src image meta data for bilinear interpolator
    t_resample_id src;       //src image meta
    t_resample_id dst;       //dst image meta
} t_resample_cbrd;            //Bilinear Resampler Data


/* this struct contains high level zoom parameters,
   all image relative */
typedef struct
{
    float centerx;
    float centery;
    float zoomx;
    float zoomy;
    float angle;
} t_resample_zrd;


/* convert floating point center and zoom data to incremental linear remapping vectors */
static void pdp_imageproc_resample_clmd_init_from_id_zrd(t_resample_clmd *l, t_resample_id *i, t_resample_zrd *z)
{
    double izx = 1.0f / (z->zoomx);
    double izy = 1.0f / (z->zoomy);
    double scale = (double)0xffffffff;
    double scalew = scale / ((double)(i->width - 1));
    double scaleh = scale / ((double)(i->height - 1));
    double cx = ((double)z->centerx) * ((double)(i->width - 1));
    double cy = ((double)z->centery) * ((double)(i->height - 1));
    double angle = z->angle * (-M_PI / 180.0);
    double c = cos(angle);
    double s = sin(angle);

    /* affine x, y mappings in screen coordinates */
    double mapx(double x, double y){return cx + izx * ( c * (x-cx) + s * (y-cy));}
    double mapy(double x, double y){return cy + izy * (-s * (x-cx) + c * (y-cy));}

    u32 tl_x = (u32)(scalew * mapx(0,0));
    u32 tl_y = (u32)(scaleh * mapy(0,0));


    u32 row_inc_x = (u32)(scalew * (mapx(1,0)-mapx(0,0)));
    u32 row_inc_y = (u32)(scaleh * (mapy(1,0)-mapy(0,0)));
    u32 col_inc_x = (u32)(scalew * (mapx(0,1)-mapx(0,0)));
    u32 col_inc_y = (u32)(scaleh * (mapy(0,1)-mapy(0,0)));


    pdp_imageproc_resample_cookedlinmap_init(l, tl_x, tl_y, row_inc_x, row_inc_y, col_inc_x, col_inc_y);
}

/* this struct contains all data for the zoom object */
typedef struct
{
    t_resample_cbrd cbrd;      // Bilinear Resampler Data
    t_resample_clmd clmd;      // Cooked Linear Mapping data
    t_resample_zrd   zrd;      // Zoom / Rotate Data
} t_resample_zoom_rotate;

// zoom + rotate
void *pdp_imageproc_resample_affinemap_new(void)
{
    t_resample_zoom_rotate *z = (t_resample_zoom_rotate *)pdp_alloc(sizeof(t_resample_zoom_rotate));
    z->zrd.centerx = 0.5;
    z->zrd.centery = 0.5;
    z->zrd.zoomx = 1.0;
    z->zrd.zoomy = 1.0;
    z->zrd.angle = 0.0f;
    return (void *)z;
}
void pdp_imageproc_resample_affinemap_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_resample_affinemap_setcenterx(void *x, float f){((t_resample_zoom_rotate *)x)->zrd.centerx = f;}
void pdp_imageproc_resample_affinemap_setcentery(void *x, float f){((t_resample_zoom_rotate *)x)->zrd.centery = f;}
void pdp_imageproc_resample_affinemap_setzoomx(void *x, float f){((t_resample_zoom_rotate *)x)->zrd.zoomx = f;}
void pdp_imageproc_resample_affinemap_setzoomy(void *x, float f){((t_resample_zoom_rotate *)x)->zrd.zoomy = f;}
void pdp_imageproc_resample_affinemap_setangle(void *x, float f){((t_resample_zoom_rotate *)x)->zrd.angle = f;}
void pdp_imageproc_resample_affinemap_process(void *x, u32 width, u32 height, s16 *srcimage, s16 *dstimage)
{
    t_resample_zoom_rotate *z = (t_resample_zoom_rotate *)x;

    /* setup resampler image meta data */
    pdp_imageproc_resample_init_id(&(z->cbrd.src), width, srcimage, width, height);
    pdp_imageproc_resample_init_id(&(z->cbrd.dst), width, dstimage, width, height);
    pdp_imageproc_resample_init_cid(&(z->cbrd.csrc),&(z->cbrd.src)); 

    /* setup linmap data from zoom_rotate parameters */
    pdp_imageproc_resample_clmd_init_from_id_zrd(&(z->clmd), &(z->cbrd.src), &(z->zrd));


    /* call assembler routine */
    pixel_resample_linmap_s16(z);   
}



// polynomials


typedef struct
{
    u32 order;
    u32 nbpasses;
    s16 coefs[0];
} t_cheby;

void *pdp_imageproc_cheby_new(int order)
{
    t_cheby *z;
    int i;
    if (order < 2) order = 2;
    z = (t_cheby *)pdp_alloc(sizeof(t_cheby) + (order + 1) * sizeof(s16[4]));
    z->order = order;
    setvec(z->coefs + 0*4, 0);
    setvec(z->coefs + 1*4, 0.25);
    for (i=2; i<=order; i++)  setvec(z->coefs + i*4, 0);

    return z;
}
void pdp_imageproc_cheby_delete(void *x){pdp_dealloc(x);}
void pdp_imageproc_cheby_setcoef(void *x, u32 n, float f)
{
    t_cheby *z = (t_cheby *)x;
    if (n <= z->order){
	setvec(z->coefs + n*4, f * 0.25); // coefs are in s2.13 format
    }
}
void pdp_imageproc_cheby_setnbpasses(void *x, u32 n){((t_cheby *)x)->nbpasses = n;}

void pdp_imageproc_cheby_process(void *x, u32 width, u32 height, s16 *image)
{
    t_cheby *z = (t_cheby *)x;
    u32 iterations = z->nbpasses;
    u32 i,j;
    for (j=0; j < (height*width); j += width)
	for (i=0; i<iterations; i++)
	    pixel_cheby_s16_3plus(image+j, width>>2, z->order+1, z->coefs);

    //pixel_cheby_s16_3plus(image, (width*height)>>2, z->order+1, z->coefs);
}





More information about the Pd-cvs mailing list