[PD-cvs] externals/pdp/system/type Makefile, 1.2, 1.3 pdp_bitmap.c, 1.2, 1.3 pdp_image.c, 1.2, 1.3 pdp_matrix.c, 1.2, 1.3

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


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

Added Files:
	Makefile pdp_bitmap.c pdp_image.c pdp_matrix.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_image.c ---
/*
 *   Pure Data Packet system implementation. : 16 bit image packet implementation
 *   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 methods for the image packets
  pdp_packet_new_* methods are several image packet constructors
  pdp_type_* are image type checkers & converters

*/

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

#include "pdp_internals.h"
#include "pdp_packet.h"
#include "pdp_imageproc.h"
#include "pdp_resample.h"
#include "pdp_list.h"
#include "pdp_post.h"
#include "pdp_type.h"


/* the class object */
static t_pdp_class* image_class;



/* check dimensions */
static void _checkdim(u32 width, u32 height){
    if 	((pdp_imageproc_legalwidth(width) != width) ||
	 (pdp_imageproc_legalheight(height) != height)){
	pdp_post("WARNING: request to create image packet with illegal dimensions %d x %d", width, height);
    }
}	    


/* image packet constructors */
int pdp_packet_new_image_YCrCb(u32 w, u32 h)
{
    t_pdp *header;
    t_image *image;
    int packet;


    u32 size = w*h;
    u32 totalnbpixels = size + (size >> 1);
    u32 packet_size = totalnbpixels << 1;

    _checkdim(w,h);

    packet = pdp_packet_new(PDP_IMAGE, packet_size);
    header = pdp_packet_header(packet);
    image = pdp_packet_image_info(packet);
    if (!header) return -1;

    image->encoding = PDP_IMAGE_YV12;
    image->width = w;
    image->height = h;
    header->desc = pdp_packet_image_get_description(packet);
    header->theclass = image_class;

    return packet;
}

int pdp_packet_new_image_grey(u32 w, u32 h)
{
    t_pdp *header;
    t_image *image;
    int packet;


    u32 size = w*h;
    u32 totalnbpixels = size;
    u32 packet_size = totalnbpixels << 1;

    _checkdim(w,h);

    packet = pdp_packet_new(PDP_IMAGE, packet_size);
    header = pdp_packet_header(packet);
    image = pdp_packet_image_info(packet);
    if (!header) return -1;

    image->encoding = PDP_IMAGE_GREY;
    image->width = w;
    image->height = h;
    header->desc = pdp_packet_image_get_description(packet);
    header->theclass = image_class;

    return packet;
}

int pdp_packet_new_image_mchp(u32 w, u32 h, u32 d)
{
    t_pdp *header;
    t_image *image;
    int packet;


    u32 size = w*h*d;
    u32 totalnbpixels = size;
    u32 packet_size = totalnbpixels << 1;

    _checkdim(w,h);

    packet = pdp_packet_new(PDP_IMAGE, packet_size);
    header = pdp_packet_header(packet);
    image = pdp_packet_image_info(packet);
    if (!header) return -1;


    image->encoding = PDP_IMAGE_MCHP;
    image->width = w;
    image->height = h;
    image->depth = d;
    header->desc = pdp_packet_image_get_description(packet);
    header->theclass = image_class;

    return packet;
}


int pdp_packet_new_image(u32 type, u32 w, u32 h)
{
    switch (type){
    case PDP_IMAGE_YV12:
	return pdp_packet_new_image_YCrCb(w,h);
    case PDP_IMAGE_GREY:
	return pdp_packet_new_image_grey(w,h);
    default:
	return -1;
    }
}


/****************** packet type checking and conversion methods ********************/



/* check if two image packets are allocated and of the same type */
int pdp_packet_image_compat(int packet0, int packet1)
{
    t_pdp *header0 = pdp_packet_header(packet0);
    t_pdp *header1 = pdp_packet_header(packet1);
    t_image *image0 = pdp_packet_image_info(packet0);
    t_image *image1 = pdp_packet_image_info(packet1);


    if (!(pdp_packet_compat(packet0, packet1))) return 0;
    if (header0->type != PDP_IMAGE){
	//pdp_post("pdp_type_compat_image: not a PDP_IMAGE");
	return 0;
    }
    if (image0->encoding != image1->encoding){
	//pdp_post("pdp_type_compat_image: encodings differ");
	return 0;
    }
    if (image0->width != image1->width){
	//pdp_post("pdp_type_compat_image: image withs differ");
	return 0;
    }
    if (image0->height != image1->height){
	//pdp_post("pdp_type_compat_image: image heights differ");
	return 0;
    }
    return 1;
}

/* check if packet is a valid image packet */
int pdp_packet_image_isvalid(int packet)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    if (!header) return 0;
    if (PDP_IMAGE != header->type) return 0;
    if ((PDP_IMAGE_YV12 != image->encoding)
	&& (PDP_IMAGE_GREY != image->encoding)
	&& (PDP_IMAGE_MCHP != image->encoding)) return 0;

    return 1;

}

/* set the channel mask for the image */
void pdp_packet_image_set_chanmask(int packet, unsigned int chanmask)
{
    if (pdp_packet_image_isvalid(packet)) pdp_packet_image_info(packet)->chanmask = chanmask;
    
}


t_image *pdp_packet_image_info(int packet)
{
    return (t_image *)pdp_packet_subheader(packet);
}


t_pdp_symbol *pdp_packet_image_get_description(int packet)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    char description[1024];
    char *c = description;
    int encoding;

    if (!header) return pdp_gensym("invalid");
    else if (!header->desc){
	/* if description is not defined, try to reconstruct it (for backwards compat) */
	if (header->type == PDP_IMAGE){
	    c += sprintf(c, "image");
	    encoding = image->encoding;
	    switch(encoding){
	    case PDP_IMAGE_YV12: c += sprintf(c, "/YCrCb"); break;
	    case PDP_IMAGE_GREY: c += sprintf(c, "/grey"); break;
	    case PDP_IMAGE_MCHP: c += sprintf(c, "/multi"); break;
	    default:
		c += sprintf(c, "/unknown"); goto exit;
	    }
	    if (encoding == PDP_IMAGE_MCHP){
		c += sprintf(c, "/%dx%dx%d", 
			     image->width, 
			     image->height,
			     image->depth);
	    }
	    else {
		c += sprintf(c, "/%dx%d", 
			     image->width, 
			     image->height);
	    }

	exit:
	    return pdp_gensym(description);
	} 
	else return pdp_gensym("unknown");
    }
    else return header->desc;
}





/* IMAGE PACKAGE CONVERSION ROUTINES */

/* note: these are internal: no extra checking is done
   it is assumed the packets are of correct type (the type template associated with the conversion program) */

// image/YCrCb/* -> image/grey/*
// image/multi/* -> image/grey/*  (only first channel)
static int _pdp_packet_image_convert_YCrCb_to_grey(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    int w = image->width;
    int h = image->height;
    int p = pdp_packet_new_image_grey(w,h);
    if (p == -1) return p;
    memcpy(pdp_packet_data(p), pdp_packet_data(packet), 2*w*h);
    return p;
}

// image/grey/* -> image/YCrCb/*
static int _pdp_packet_image_convert_grey_to_YCrCb(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    int w = image->width;
    int h = image->height;
    int p = pdp_packet_new_image_YCrCb(w,h);
    int y_bytes = 2*w*h;
    void *data;
    if (p == -1) return p;
    data = pdp_packet_data(p);
    memcpy(data, pdp_packet_data(packet), y_bytes);
    memset(data+y_bytes, 0, y_bytes >> 1);
    return p;
}

// image/grey/* -> image/multi/*
static int _pdp_packet_image_convert_grey_to_multi(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    int w = image->width;
    int h = image->height;
    int p = pdp_packet_new_image_mchp(w,h,1);
    int y_bytes = 2*w*h;
    void *data;
    if (p == -1) return p;
    data = pdp_packet_data(p);
    memcpy(data, pdp_packet_data(packet), y_bytes);
    return p;
}

// image/multi/* -> image/YCrCb/*  (only first 3 channels)
static int _pdp_packet_image_convert_multi_to_YCrCb(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    s16 *data, *newdata;

    /* get info */
    int w = image->width;
    int h = image->height;
    int d = image->depth;
    int plane_size = w*h;

    /* create new packet */
    int newpacket = pdp_packet_new_image_YCrCb(w, h);
    if (-1 == newpacket) return -1;

    data = pdp_packet_data(packet);
    newdata = pdp_packet_data(newpacket);

    /* copy channel 0 */
    memcpy(newdata, data, plane_size<<1);
    newdata += plane_size;
    data += plane_size;

    /* copy channel 1 */
    if (d >= 1) pdp_resample_halve(data, newdata, w, h);
    else        memset(newdata, 0, plane_size >> 1);
    data += plane_size;
    newdata += (plane_size >> 2);


    /* copy channel 2 */
    if (d >= 2) pdp_resample_halve(data, newdata, w, h);
    else        memset(newdata, 0, plane_size >> 1);

    return newpacket;
    
}

// image/YCrCb/* -> image/multi/*
static int _pdp_packet_image_convert_YCrCb_to_multi(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    s16 *data, *newdata;

    /* get info */
    int w = image->width;
    int h = image->height;
    int plane_size = w*h;

    /* create new packet */
    int newpacket = pdp_packet_new_image_mchp(w, h, 3);
    if (-1 == newpacket) return -1;

    data = pdp_packet_data(packet);
    newdata = pdp_packet_data(newpacket);

    /* copy channel 0 */
    memcpy(newdata, data, plane_size<<1);
    newdata += plane_size;
    data += plane_size;
    w >>= 1;
    h >>= 1;

    /* copy channel 1 */
    pdp_resample_double(data, newdata, w, h);
    data += (plane_size >> 2);
    newdata += plane_size;

    /* copy channel 2 */
    pdp_resample_double(data, newdata, w, h);

    return newpacket;
    
}

static void _pdp_description_get_dims(t_pdp_symbol *template, int *w, int *h, int *d)
{
    char *c = template->s_name;
    // get requested dimensions
    *w = 0;
    *h = 0;
    *d = 0;
    while (*c++ != '/');
    while (*c++ != '/');
    sscanf(c, "%dx%dx%d", w, h, d);

}

// resample image/YCrCb/*
static int _pdp_packet_image_convert_resample_YCrCb(int packet, t_pdp_symbol *dest_template)
{
    int quality = 1;
    int dst_w, dst_h, dummy;
    int new_packet;
    unsigned int src_size, src_voffset, src_uoffset;
    unsigned int dst_size, dst_voffset, dst_uoffset;
    t_pdp *header0 = pdp_packet_header(packet);
    t_image *image0 = pdp_packet_image_info(packet);
    unsigned int src_w = image0->width;
    unsigned int src_h = image0->height;
    short int *src_image = (short int *)pdp_packet_data(packet);
    short int *dst_image;
    _pdp_description_get_dims(dest_template, &dst_w, &dst_h, &dummy);
    new_packet = pdp_packet_new_image_YCrCb(dst_w, dst_h);
    if (-1 == new_packet) return -1;
    dst_image = (short int*)pdp_packet_data(new_packet);
    src_size = src_w*src_h;
    src_voffset = src_size;
    src_uoffset = src_size + (src_size>>2);
    dst_size = dst_w*dst_h;
    dst_voffset = dst_size;
    dst_uoffset = dst_size + (dst_size>>2);
    if (quality){
	pdp_resample_scale_bilin(src_image, dst_image, src_w, src_h, dst_w, dst_h);
	pdp_resample_scale_bilin(src_image+src_voffset, dst_image+dst_voffset, src_w>>1, src_h>>1, dst_w>>1, dst_h>>1);
	pdp_resample_scale_bilin(src_image+src_uoffset, dst_image+dst_uoffset, src_w>>1, src_h>>1, dst_w>>1, dst_h>>1);
    }
    else{
	pdp_resample_scale_nn(src_image, dst_image, src_w, src_h, dst_w, dst_h);
	pdp_resample_scale_nn(src_image+src_voffset, dst_image+dst_voffset, src_w>>1, src_h>>1, dst_w>>1, dst_h>>1);
	pdp_resample_scale_nn(src_image+src_uoffset, dst_image+dst_uoffset, src_w>>1, src_h>>1, dst_w>>1, dst_h>>1);
    }
    return new_packet;
}

// resample image/grey/* and image/multi/*
static int _pdp_packet_image_convert_resample_multi(int packet, t_pdp_symbol *dest_template)
{
    int quality = 1;
    int dst_w, dst_h, depth;
    int new_packet;
    unsigned int src_size;
    unsigned int dst_size;
    t_pdp *header0 = pdp_packet_header(packet);
    t_image *image0 = pdp_packet_image_info(packet);
    unsigned int src_w = image0->width;
    unsigned int src_h = image0->height;
    short int *src_image = (short int *)pdp_packet_data(packet);
    short int *dst_image;
    _pdp_description_get_dims(dest_template, &dst_w, &dst_h, &depth);

    if (depth == 0){
	depth = 1;
	new_packet = pdp_packet_new_image_grey(dst_w, dst_h);
    }
    else {
	new_packet = pdp_packet_new_image_mchp(dst_w, dst_h, depth);
    }
    if (-1 == new_packet) return -1;

    dst_image = (short int*)pdp_packet_data(new_packet);

    src_size = src_w*src_h;
    dst_size = dst_w*dst_h;
    while (depth--){
	if (quality){
	    pdp_resample_scale_bilin(src_image, dst_image, src_w, src_h, dst_w, dst_h);
	}
	else{
	    pdp_resample_scale_nn(src_image, dst_image, src_w, src_h, dst_w, dst_h);
	}
	src_image += src_size;
	dst_image += dst_size;
    }

    return new_packet;
}

static int _pdp_packet_image_convert_fallback(int packet, t_pdp_symbol *dest_template)
{
    pdp_post("can't convert image type %s to %s",
	 pdp_packet_get_description(packet)->s_name, dest_template->s_name);

    return -1;
}



/* the expensive factory method */
static int pdp_image_factory(t_pdp_symbol *type)
{
    t_pdp_list *l;
    t_pdp_symbol *s;
    int t;
    int w = 0;
    int h = 0;
    int d = 0;
    int p = -1;
    int n = 0;
    int m = 0;
    char *garbage = 0;

    //pdp_post("creating:");
    //pdp_post("%s", type->s_name);

    l = pdp_type_to_list(type);
    s = pdp_list_pop(l).w_symbol; // first element is "image"
    s = pdp_list_pop(l).w_symbol;

    /* get image type */
    if (s == pdp_gensym("grey")) t = PDP_IMAGE_GREY;
    else if (s == pdp_gensym("YCrCb")) t = PDP_IMAGE_YV12;
    else if (s == pdp_gensym("multi")) t = PDP_IMAGE_MCHP;
    else goto exit;
    
    /* get image dimensions and create image */
    s = pdp_list_pop(l).w_symbol;
    switch (t){
    case PDP_IMAGE_MCHP:
	m  = sscanf(s->s_name, "%dx%dx%d", &w, &h, &d);
	p = pdp_packet_new_image_mchp(w,h,d);
	break;
    default:
	sscanf(s->s_name, "%dx%d", &w, &h);
	p = pdp_packet_new_image(t,w,h);
	break;
    }
    if (p != -1){
	t_pdp *h = pdp_packet_header(p);
	/* if type is not exact, delete the packet */
	if (type != h->desc) {
	    pdp_packet_delete(p);
	    p = -1;
	}
    }
 exit:
    pdp_list_free(l);
    return p;
}



void pdp_image_words_setup(t_pdp_class *c);

void pdp_image_setup(void)
{
    t_pdp_conversion_program *program;

    /* setup the class object */
    image_class = pdp_class_new(pdp_gensym("image/*/*"), pdp_image_factory);
    

    /* setup conversion programs */
    program = pdp_conversion_program_new(_pdp_packet_image_convert_YCrCb_to_grey, 0);
    pdp_type_register_conversion(pdp_gensym("image/YCrCb/*"), pdp_gensym("image/grey/*"), program);
    pdp_type_register_conversion(pdp_gensym("image/multi/*"), pdp_gensym("image/grey/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_image_convert_grey_to_YCrCb, 0);
    pdp_type_register_conversion(pdp_gensym("image/grey/*"), pdp_gensym("image/YCrCb/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_image_convert_grey_to_multi, 0);
    pdp_type_register_conversion(pdp_gensym("image/grey/*"), pdp_gensym("image/multi/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_image_convert_multi_to_YCrCb, 0);
    pdp_type_register_conversion(pdp_gensym("image/multi/*"), pdp_gensym("image/YCrCb/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_image_convert_YCrCb_to_multi, 0);
    pdp_type_register_conversion(pdp_gensym("image/YCrCb/*"), pdp_gensym("image/multi/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_image_convert_resample_YCrCb, 0);
    pdp_type_register_conversion(pdp_gensym("image/YCrCb/*"), pdp_gensym("image/YCrCb/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_image_convert_resample_multi, 0);
    pdp_type_register_conversion(pdp_gensym("image/multi/*"), pdp_gensym("image/multi/*"), program);
    pdp_type_register_conversion(pdp_gensym("image/grey/*"), pdp_gensym("image/grey/*"), program);

    /* catch-all fallback */
    program = pdp_conversion_program_new(_pdp_packet_image_convert_fallback, 0);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("image/*/*"), program);
							     
}

--- NEW FILE: pdp_bitmap.c ---
/*
 *   Pure Data Packet system implementation. : 8 bit image packet implementation
 *   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 <stdio.h>
#include <string.h>
#include "pdp_packet.h"
#include "pdp_type.h"
#include "pdp_llconv.h"
#include "pdp_internals.h"
#include "pdp_post.h"


/* the class object */
static t_pdp_class* bitmap_class;


t_bitmap *pdp_packet_bitmap_info(int packet)
{
    return (t_bitmap *)pdp_packet_subheader(packet);
}

/* check if packet is a valid image packet */
int pdp_packet_bitmap_isvalid(int packet)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *image = pdp_packet_bitmap_info(packet);
    if (!header) return 0;
    if (PDP_BITMAP != header->type) return 0;

    return 1;

}



/* bitmap constructors */
t_pdp_symbol *pdp_packet_bitmap_get_description(int packet)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    char description[1024];
    char *c = description;
    int encoding;

    if (!header) return pdp_gensym("invalid");
    else if (!header->desc){
	/* if description is not defined, try to reconstruct it (for backwards compat) */
	if (header->type == PDP_BITMAP){
	    c += sprintf(c, "bitmap");
	    encoding = bitmap->encoding;
	    switch(encoding){
	    case PDP_BITMAP_RGB: c += sprintf(c, "/rgb"); break;
	    case PDP_BITMAP_RGBA: c += sprintf(c, "/rgba"); break;
	    case PDP_BITMAP_GREY: c += sprintf(c, "/grey"); break;
	    case PDP_BITMAP_YV12: c += sprintf(c, "/yv12"); break;
	    default:
		c += sprintf(c, "/unknown"); goto exit;
	    }
	    c += sprintf(c, "/%dx%d", 
			 bitmap->width, 
			 bitmap->height);
	exit:
	    return pdp_gensym(description);
	} 
	else return pdp_gensym("unknown");
    }
    else return header->desc;
}

int pdp_packet_new_bitmap_yv12(u32 w, u32 h)
{
    t_pdp *header;
    t_bitmap *bitmap;
    int packet;


    u32 size = w*h;
    u32 totalnbpixels = size;
    u32 packet_size = size + (size >> 1);

    packet = pdp_packet_new(PDP_BITMAP, packet_size);
    header = pdp_packet_header(packet);
    bitmap = pdp_packet_bitmap_info(packet);
    if (!header) return -1;

    bitmap->encoding = PDP_BITMAP_YV12;
    bitmap->width = w;
    bitmap->height = h;
    header->desc = pdp_packet_bitmap_get_description(packet);
    header->theclass = bitmap_class;

    return packet;
}

int pdp_packet_new_bitmap_grey(u32 w, u32 h)
{
    t_pdp *header;
    t_bitmap *bitmap;
    int packet;

    u32 size = w*h;
    u32 totalnbpixels = size;
    u32 packet_size = totalnbpixels;

    packet = pdp_packet_new(PDP_BITMAP, packet_size);
    header = pdp_packet_header(packet);
    bitmap = pdp_packet_bitmap_info(packet);
    if (!header) return -1;

    bitmap->encoding = PDP_BITMAP_GREY;
    bitmap->width = w;
    bitmap->height = h;
    header->desc = pdp_packet_bitmap_get_description(packet);
    header->theclass = bitmap_class;

    return packet;
}

int pdp_packet_new_bitmap_rgb(u32 w, u32 h)
{
    t_pdp *header;
    t_bitmap *bitmap;
    int packet;


    u32 size = w*h;
    u32 totalnbpixels = size;
    u32 packet_size = totalnbpixels * 3;

    packet = pdp_packet_new(PDP_BITMAP, packet_size);
    header = pdp_packet_header(packet);
    bitmap = pdp_packet_bitmap_info(packet);
    if (!header) return -1;

    bitmap->encoding = PDP_BITMAP_RGB;
    bitmap->width = w;
    bitmap->height = h;
    header->desc = pdp_packet_bitmap_get_description(packet);
    header->theclass = bitmap_class;

    return packet;
}

int pdp_packet_new_bitmap_rgba(u32 w, u32 h)
{
    t_pdp *header;
    t_bitmap *bitmap;
    int packet;


    u32 size = w*h;
    u32 totalnbpixels = size;
    u32 packet_size = totalnbpixels * 4;

    packet = pdp_packet_new(PDP_BITMAP, packet_size);
    header = pdp_packet_header(packet);
    bitmap = pdp_packet_bitmap_info(packet);
    if (!header) return -1;

    bitmap->encoding = PDP_BITMAP_RGBA;
    bitmap->width = w;
    bitmap->height = h;
    header->desc = pdp_packet_bitmap_get_description(packet);
    header->theclass = bitmap_class;

    return packet;
}

int pdp_packet_new_bitmap(int type, u32 w, u32 h)
{
    switch(type){
    case PDP_BITMAP_GREY: return pdp_packet_new_bitmap_grey(w,h);
    case PDP_BITMAP_YV12: return pdp_packet_new_bitmap_yv12(w,h);
    case PDP_BITMAP_RGB:  return pdp_packet_new_bitmap_rgb(w,h);
    case PDP_BITMAP_RGBA: return pdp_packet_new_bitmap_rgba(w,h);
    default: return -1;
    }
}


/* some utility methods */
void pdp_llconv_flip_top_bottom(char *data, int width, int height, int pixelsize);

/* flip top & bottom */
void pdp_packet_bitmap_flip_top_bottom(int packet)
{
    t_bitmap *b = (t_bitmap *)pdp_packet_subheader(packet);
    char *d = (char *)pdp_packet_data(packet);
    int w,h;
    if (!pdp_packet_bitmap_isvalid(packet)) return;
    if (!b) return;
    w = b->width;
    h = b->height;

    switch(b->encoding){
    case PDP_BITMAP_GREY: pdp_llconv_flip_top_bottom(d,w,h,1); break;
    case PDP_BITMAP_RGB:  pdp_llconv_flip_top_bottom(d,w,h,3); break;
    case PDP_BITMAP_RGBA: pdp_llconv_flip_top_bottom(d,w,h,4); break;
    default: break;
    }

}



/* conversion methods */
// image/grey/* -> bitmap/grey/*
static int _pdp_packet_bitmap_convert_grey_to_grey8(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    s16 *data = (s16 *)pdp_packet_data(packet);
    u8 *new_data;
    u32 w,h;
    int new_p;

    if (!pdp_packet_image_isvalid(packet)) return -1;
    w = image->width;
    h = image->height;

    if (!((image->encoding == PDP_IMAGE_GREY) ||
	(image->encoding == PDP_IMAGE_YV12))) return -1;

    new_p = pdp_packet_new_bitmap_grey(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    /* convert image to greyscale 8 bit */
    pdp_llconv(data,RIF_GREY______S16, new_data, RIF_GREY______U8, w, h);

    return new_p;
}

static int _pdp_packet_bitmap_convert_YCrCb_to_rgb8(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    s16 *data = (s16 *)pdp_packet_data(packet);
    u8 *new_data;
    u32 w,h;
    int new_p;

    if (!pdp_packet_image_isvalid(packet)) return -1;
    w = image->width;
    h = image->height;

    if (!((image->encoding == PDP_IMAGE_YV12))) return -1;

    new_p = pdp_packet_new_bitmap_rgb(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    /* convert image to greyscale 8 bit */
    pdp_llconv(data,RIF_YVU__P411_S16, new_data, RIF_RGB__P____U8, w, h);

    return new_p;
}

static int _pdp_packet_bitmap_convert_image_to_yv12(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    s16 *data = (s16 *)pdp_packet_data(packet);
    u8 *new_data;
    u32 w,h, nbpixels;
    int new_p;
    int encoding = image->encoding;

    if (!pdp_packet_image_isvalid(packet)) return -1;
    w = image->width;
    h = image->height;
    nbpixels = w*h;

    new_p = pdp_packet_new_bitmap_yv12(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    switch (encoding){
    case PDP_IMAGE_YV12:
	pdp_llconv(data, RIF_YVU__P411_S16, new_data, RIF_YVU__P411_U8, w, h);
	break;
    case PDP_IMAGE_GREY:
	pdp_llconv(data, RIF_GREY______S16, new_data, RIF_GREY______U8, w, h);
	memset(new_data + nbpixels, 0x80, nbpixels>>1);
	break;
    default:
	/* not supported, $$$TODO add more */
	pdp_packet_mark_unused(new_p);
	new_p = -1;
	break;
    }

    return new_p;

}

static int _pdp_packet_bitmap_convert_rgb8_to_YCrCb(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    u8 *data = (u8 *)pdp_packet_data(packet);
    s16 *new_data;
    u32 w,h;
    int new_p;

    w = bitmap->width;
    h = bitmap->height;
    new_p = pdp_packet_new_image_YCrCb(w,h);
    if (-1 == new_p) return -1;
    new_data = (s16 *)pdp_packet_data(new_p);

    /* convert image to greyscale 8 bit */
    pdp_llconv(data, RIF_RGB__P____U8, new_data, RIF_YVU__P411_S16, w, h);

    return new_p;
}

static int _pdp_packet_bitmap_convert_rgb8_to_bitmap_YCrCb(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    u8 *data = (u8 *)pdp_packet_data(packet);
    u8 *new_data;
    u32 w,h;
    int new_p;

    w = bitmap->width;
    h = bitmap->height;
    new_p = pdp_packet_new_bitmap_yv12(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    /* convert image to greyscale 8 bit */
    pdp_llconv(data, RIF_RGB__P____U8, new_data, RIF_YVU__P411_U8, w, h);

    return new_p;
}

static int _pdp_packet_bitmap_convert_grey8_to_grey(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    s16 *data = (s16 *)pdp_packet_data(packet);
    u8 *new_data;
    u32 w,h;
    int new_p;

    w = bitmap->width;
    h = bitmap->height;
    new_p = pdp_packet_new_image_grey(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    /* convert image to greyscale 8 bit */
    pdp_llconv(data, RIF_GREY______U8, new_data, RIF_GREY______S16, w, h);

    return new_p;
}
static int _pdp_packet_bitmap_convert_rgb8_to_rgba8(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    u8 *data = (u8 *)pdp_packet_data(packet);
    u8 *new_data;
    int w,h, new_p, i;

    w = bitmap->width;
    h = bitmap->height;
    new_p = pdp_packet_new_bitmap_rgba(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    for(i=0; i<w*h; i++){
	new_data[4*i+0]    = data[3*i + 0];
	new_data[4*i+1]    = data[3*i + 1];
	new_data[4*i+2]    = data[3*i + 2];
	new_data[4*i+3]    = 0;
    }

    return new_p;
}
static int _pdp_packet_bitmap_convert_rgba8_to_rgb8(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    u8 *data = (u8 *)pdp_packet_data(packet);
    u8 *new_data;
    int w,h, new_p, i;

    w = bitmap->width;
    h = bitmap->height;
    new_p = pdp_packet_new_bitmap_rgb(w,h);
    if (-1 == new_p) return -1;
    new_data = (u8 *)pdp_packet_data(new_p);

    for(i=0; i<w*h; i++){
	new_data[3*i+0]    = data[4*i + 0];
	new_data[3*i+1]    = data[4*i + 1];
	new_data[3*i+2]    = data[4*i + 2];
    }

    return new_p;
}
static int _pdp_packet_bitmap_convert_rgb8_to_mchp(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    u8 *data = (u8 *)pdp_packet_data(packet);
    s16 *new_data;
    int w,h;
    int new_p, i, plane;

    w = bitmap->width;
    h = bitmap->height;
    plane = w*h;
    new_p = pdp_packet_new_image_multi(w,h,3);
    if (-1 == new_p) return -1;
    new_data = (s16 *)pdp_packet_data(new_p);

    for(i=0; i<w*h; i++){
	new_data[i]          = ((u32)data[3*i + 0]) << 7;
	new_data[i+plane]    = ((u32)data[3*i + 1]) << 7;
	new_data[i+plane*2]  = ((u32)data[3*i + 2]) << 7;
    }

    return new_p;
}

static int _pdp_packet_bitmap_convert_yv12_to_image(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_bitmap *bitmap = pdp_packet_bitmap_info(packet);
    u8 *data = (u8 *)pdp_packet_data(packet);
    s16 *new_data;
    int w,h;
    int new_p;

    w = bitmap->width;
    h = bitmap->height;
    new_p = pdp_packet_new_image_YCrCb(w,h);
    new_data = pdp_packet_data(new_p);
    if (-1 == new_p || !new_data) return -1;
    pdp_llconv(data, RIF_YVU__P411_U8, new_data, RIF_YVU__P411_S16, w, h);

    return new_p;
}

static int _pdp_packet_bitmap_convert_mchp_to_rgb8(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = (t_image *)pdp_packet_subheader(packet);
    s16 *data = (s16 *)pdp_packet_data(packet);
    u8 *new_data;
    int w = image->width;
    int h = image->height;
    int plane = w*h;
    int nb_channels = image->depth;
    int new_p, i;

    static inline u8 _map(s32 pixel){
	s32 mask = ~(pixel>>16);
	return ((pixel >> 7) & mask);
    }

    switch(nb_channels){
    default: return -1;
    case 1:
	if (-1 == (new_p = pdp_packet_new_bitmap_grey(w,h))) return -1;
	new_data = (u8*)pdp_packet_data(new_p);
	for(i=0; i<plane; i++) new_data[i] = _map(data[i]); 
	break;
    case 3:
	if (-1 == (new_p = pdp_packet_new_bitmap_rgb(w,h)))  return -1;
	new_data = (u8*)pdp_packet_data(new_p);
	for(i=0; i<plane; i++){
	    new_data[3*i+0] = _map(data[i]); 
	    new_data[3*i+1] = _map(data[i+plane]); 
	    new_data[3*i+2] = _map(data[i+plane*2]);
	}
	break;
    case 4:
	if (-1 == (new_p = pdp_packet_new_bitmap_rgba(w,h))) return -1;
	new_data = (u8*)pdp_packet_data(new_p);
	for(i=0; i<plane; i++){
	    new_data[4*i+0] = _map(data[i]); 
	    new_data[4*i+1] = _map(data[i+plane]); 
	    new_data[4*i+2] = _map(data[i+plane*2]);
	    new_data[4*i+3] = _map(data[i+plane*3]);
	}
	break;
	
	
    }

    return new_p;

}

static int _pdp_packet_bitmap_convert_fallback(int packet, t_pdp_symbol *dest_template)
{
    pdp_post("can't convert image type %s to %s",
	 pdp_packet_get_description(packet)->s_name, dest_template->s_name);

    return -1;
}

static int pdp_bitmap_factory(t_pdp_symbol *type)
{
    t_pdp_list *l;
    t_pdp_symbol *s;
    int t;
    int w = 0;
    int h = 0;
    int d = 0;
    int p = -1;
    int n = 0;
    int m = 0;
    char *garbage = 0;

    //pdp_post("creating:");
    //pdp_post("%s", type->s_name);

    l = pdp_type_to_list(type);
    s = pdp_list_pop(l).w_symbol; // first element is "bitmap"
    s = pdp_list_pop(l).w_symbol;

    /* get image type */
    if (s == pdp_gensym("grey")) t = PDP_BITMAP_GREY;
    else if (s == pdp_gensym("yv12")) t = PDP_BITMAP_YV12;
    else if (s == pdp_gensym("rgb")) t = PDP_BITMAP_RGB;
    else goto exit;
    
    /* get image dimensions and create image */
    s = pdp_list_pop(l).w_symbol;
    sscanf(s->s_name, "%dx%d", &w, &h);
    p = pdp_packet_new_bitmap(t,w,h);

    if (p != -1){
	t_pdp *h = pdp_packet_header(p);
	/* if type is not exact, delete the packet */
	if (type != h->desc) {
	    pdp_packet_delete(p);
	    p = -1;
	}
    }
 exit:
    pdp_list_free(l);
    return p;
}

void pdp_bitmap_setup(void)
{
    t_pdp_conversion_program *program;

    /* setup class object */
    bitmap_class = pdp_class_new(pdp_gensym("bitmap/*/*"), pdp_bitmap_factory);
				 

    /* setup conversion programs */
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_grey_to_grey8, 0);
    pdp_type_register_conversion(pdp_gensym("image/grey/*"), pdp_gensym("bitmap/grey/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_grey8_to_grey, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/grey/*"), pdp_gensym("image/grey/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_YCrCb_to_rgb8, 0);
    pdp_type_register_conversion(pdp_gensym("image/YCrCb/*"), pdp_gensym("bitmap/rgb/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_image_to_yv12, 0);
    pdp_type_register_conversion(pdp_gensym("image/YCrCb/*"), pdp_gensym("bitmap/yv12/*"), program);
    pdp_type_register_conversion(pdp_gensym("image/grey/*"), pdp_gensym("bitmap/yv12/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_yv12_to_image, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/yv12/*"), pdp_gensym("image/YCrCb/*"), program);

    /* rgb->YCrCb converts the colour space */
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_rgb8_to_YCrCb, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/rgb/*"), pdp_gensym("image/YCrCb/*"), program);


    /* rgb <-> multi does not convert the colour space */
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_rgb8_to_mchp, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/rgb/*"), pdp_gensym("image/multi/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_mchp_to_rgb8, 0);
    pdp_type_register_conversion(pdp_gensym("image/multi/*"), pdp_gensym("bitmap/*/*"), program);


    /* rgb <-> rgba */
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_rgb8_to_rgba8, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/rgb/*"), pdp_gensym("bitmap/rgba/*"), program);
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_rgba8_to_rgb8, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/rgba/*"), pdp_gensym("bitmap/rgb/*"), program);


    /* fallback rgb convertor */
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_rgb8_to_YCrCb, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/rgb/*"), pdp_gensym("image/*/*"), program);

    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_rgb8_to_bitmap_YCrCb, 0);
    pdp_type_register_conversion(pdp_gensym("bitmap/rgb/*"), pdp_gensym("bitmap/yv12/*"), program);


    /* fallbacks */
    program = pdp_conversion_program_new(_pdp_packet_bitmap_convert_fallback, 0);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("bitmap/*/*"), program);
    pdp_type_register_conversion(pdp_gensym("bitmap/*/*"), pdp_gensym("image/*/*"), program);
    pdp_type_register_conversion(pdp_gensym("bitmap/*/*"), pdp_gensym("bitmap/*/*"), program);

}

--- NEW FILE: Makefile ---

OBJECTS = pdp_bitmap.o pdp_image.o $(PDP_OPTTYPES)


include ../../Makefile.config

all: $(OBJECTS)

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

--- NEW FILE: pdp_matrix.c ---

#include <string.h>
#include "pdp_matrix.h"
#include "pdp_packet.h"
#include "pdp_symbol.h"
#include "pdp_internals.h" // for pdp_packet_new, which is not a proper constructor
#include "pdp_post.h"
#include "pdp_type.h"

/* the class object */
static t_pdp_class *matrix_class;



/* declare header and subheader variables and exit with errval if invalid */
#define VALID_MATRIX_HEADER(packet, header, subheader, mtype, errval)   \
t_pdp *    header    = pdp_packet_header( packet );			\
t_matrix * subheader = (t_matrix *)pdp_packet_subheader( packet );	\
if (! header ) return errval;						\
if (PDP_MATRIX != header->type) return errval;				\
if (mtype) {if (subheader->type != mtype) return errval;}


int pdp_packet_matrix_isvalid(int p)
{
    VALID_MATRIX_HEADER(p, h, m, 0, 0);
    return 1;
}


void *pdp_packet_matrix_get_gsl_matrix(int p, u32 type)
{
    VALID_MATRIX_HEADER(p, h, m, type, 0);
    return &m->matrix;
}

void *pdp_packet_matrix_get_gsl_vector(int p, u32 type)
{
    VALID_MATRIX_HEADER(p, h, m, type, 0);
    return &m->vector;
}

int pdp_packet_matrix_get_type(int p)
{
    VALID_MATRIX_HEADER(p, h, m, 0, 0);
    return m->type;
}

int pdp_packet_matrix_isvector(int p)
{
    VALID_MATRIX_HEADER(p, h, m, 0, 0);
    return ((m->rows == 1) || (m->columns == 1));
}

int pdp_packet_matrix_ismatrix(int p)
{
    VALID_MATRIX_HEADER(p, h, m, 0, 0);
    return ((m->rows != 1) && (m->columns != 1));
}




/* gsl blas matrix/vector multiplication:

vector.vector

 (const gsl_vector * x, const gsl_vector * y, double * result) 

gsl_blas_sdot
gsl_blas_ddot
gsl_blas_cdot
gsl_blas_zdot
gsl_blas_cdotu
gsl_blas_zdotu

matrix.vector
 (
   CBLAS_TRANSPOSE_t 
   TransA, 
   double alpha, 
   const gsl_matrix * A, 
   const gsl_vector * x, 
   double beta, 
   gsl_vector * y
 )

gsl_blas_sgemv
gsl_blas_dgemv
gsl_blas_cgemv
gsl_blas_zgemv

matrix.matrix
 (
   CBLAS_TRANSPOSE_t TransA, 
   CBLAS_TRANSPOSE_t TransB, 
   double alpha, 
   const gsl_matrix * A, 
   const gsl_matrix * B, 
   double beta, 
   gsl_matrix * C
 )

gsl_blas_sgemm
gsl_blas_dgemm
gsl_blas_cgemm
gsl_blas_zgemm

*/

/* compute the matrix inverse using the LU decomposition */
/* it only works for double real/complex */
int pdp_packet_matrix_LU_to_inverse(int p)
{
    int new_p;
    u32 n;
    int type = pdp_packet_matrix_get_type(p);
    t_matrix *sheader = (t_matrix *)pdp_packet_subheader(p);
    gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p, type);
    gsl_matrix *im;
    if (!m) return -1;
    n = m->gsl_rows; 
    if (n != m->gsl_columns) return -1;
    new_p = pdp_packet_new_matrix(n, n, type);
    if (-1 == new_p) return -1;
    im = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(new_p, type);

    switch(type){
    case PDP_MATRIX_TYPE_RDOUBLE:
	gsl_linalg_LU_invert (m, &sheader->perm, im); break;
    case PDP_MATRIX_TYPE_CDOUBLE:
	gsl_linalg_complex_LU_invert ((gsl_matrix_complex *)m, &sheader->perm, (gsl_matrix_complex *)im); break;
    default:
	pdp_packet_mark_unused(new_p);
	new_p = -1;
    }
    return new_p;
}
/* compute the LU decomposition of a square matrix */
/* it only works for double real/complex */
int pdp_packet_matrix_LU(int p)
{
    int p_LU, bytes;
    u32 n;
    int type = pdp_packet_matrix_get_type(p);
    t_matrix *sh_m_LU;
    t_matrix *sh_m = (t_matrix *)pdp_packet_subheader(p);
    gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p, type);
    gsl_matrix *m_LU;
    if (!m) return -1;
    n = m->gsl_rows; 
    if (n != m->gsl_columns) return -1;
    p_LU = pdp_packet_new_matrix(n, n, type);
    if (-1 == p_LU) return -1;
    sh_m_LU  = (t_matrix *)pdp_packet_subheader(p_LU);
    m_LU = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p_LU, type);
    /* copy matrix data: move this to copy method */
    memcpy(pdp_packet_data(p_LU), pdp_packet_data(p), sh_m->block.size); 

    switch(type){
    case PDP_MATRIX_TYPE_RDOUBLE:
	gsl_linalg_LU_decomp (m_LU, &sh_m_LU->perm, &sh_m_LU->signum); break; 
    case PDP_MATRIX_TYPE_CDOUBLE:
	gsl_linalg_complex_LU_decomp ((gsl_matrix_complex *)m_LU, &sh_m_LU->perm, &sh_m_LU->signum); break; 
    default:
	pdp_packet_mark_unused(p_LU);
	p_LU = -1;
    }
    return p_LU;
}

int pdp_packet_matrix_LU_solve(int p_matrix, int p_vector)
{
    int type = pdp_packet_matrix_get_type(p_matrix);
    int p_result_vector = pdp_packet_clone_rw(p_vector);
    gsl_matrix *m = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(p_matrix, type);
    gsl_vector *v = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p_vector, type);
    gsl_vector *rv = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(p_result_vector, type);
    t_matrix *sh_m = (t_matrix *)pdp_packet_subheader(p_matrix);

    if (!(m && v && rv)) goto error;

    switch(type){
    case PDP_MATRIX_TYPE_RDOUBLE:
	if (gsl_linalg_LU_solve (m, &sh_m->perm, v, rv)) goto error; break;
    case PDP_MATRIX_TYPE_CDOUBLE:
	if(gsl_linalg_complex_LU_solve ((gsl_matrix_complex*)m, &sh_m->perm, 
		(gsl_vector_complex *)v, (gsl_vector_complex *)rv)) goto error; break;
    default:
	goto error;
    }
    return p_result_vector;

 error:
    pdp_packet_mark_unused(p_result_vector);
    //post("error");
    return -1;
}

/* matrix matrix mul: C is defining type 
   returns 0 on success */
int pdp_packet_matrix_blas_mm
(
 CBLAS_TRANSPOSE_t TransA,
 CBLAS_TRANSPOSE_t TransB,
 int pA,
 int pB,
 int pC,
 float scale_r,
 float scale_i
)
{
    gsl_complex_float cf_scale = {{scale_r, scale_i}};
    gsl_complex cd_scale = {{(double)scale_r, (double)scale_i}};
    gsl_complex_float cf_one = {{1.0f, 0.0f}};
    gsl_complex cd_one = {{1.0, 0.0}};
    gsl_matrix *mA, *mB, *mC;
    int type;
    type = pdp_packet_matrix_get_type(pC);
    mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA,type);
    mB = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pB,type);
    mC = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pC,type);

    if (!(mA && mB)) return 1;


    switch(type){
    case PDP_MATRIX_TYPE_RFLOAT:
	return gsl_blas_sgemm(TransA, TransB, scale_r, (gsl_matrix_float *)mA, 
			      (gsl_matrix_float *)mB, 1.0f, (gsl_matrix_float *)mC); 
    case PDP_MATRIX_TYPE_RDOUBLE:
	return gsl_blas_dgemm(TransA, TransB, (double)scale_r, (gsl_matrix *)mA, 
			      (gsl_matrix *)mB, 1.0, (gsl_matrix *)mC);
    case PDP_MATRIX_TYPE_CFLOAT:
	return gsl_blas_cgemm(TransA, TransB, cf_scale, (gsl_matrix_complex_float *)mA, 
			      (gsl_matrix_complex_float *)mB, cf_one, (gsl_matrix_complex_float *)mC); 
    case PDP_MATRIX_TYPE_CDOUBLE:
	return gsl_blas_zgemm(TransA, TransB, cd_scale, (gsl_matrix_complex *)mA, 
			      (gsl_matrix_complex *)mB, cd_one, (gsl_matrix_complex *)mC); 
    default:
	return 0;
    }
}

/* matrix vector mul: C is defining type 
   returns 0 on success */
int pdp_packet_matrix_blas_mv
(
 CBLAS_TRANSPOSE_t TransA,
 int pA,
 int pb,
 int pc,
 float scale_r,
 float scale_i
)
{
    gsl_complex_float cf_scale = {{scale_r, scale_i}};
    gsl_complex cd_scale = {{(double)scale_r, (double)scale_i}};
    gsl_complex_float cf_one = {{1.0f, 0.0f}};
    gsl_complex cd_one = {{1.0, 0.0}};
    gsl_matrix *mA;
    gsl_vector *vb, *vc;
    int type;
    type = pdp_packet_matrix_get_type(pA);
    mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA,type);
    vb = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(pb,type);
    vc = (gsl_vector *)pdp_packet_matrix_get_gsl_vector(pc,type);

    if (!(vb && vc)) return 1;


    switch(type){
    case PDP_MATRIX_TYPE_RFLOAT:
	return gsl_blas_sgemv(TransA, scale_r, (gsl_matrix_float *)mA, 
			      (gsl_vector_float *)vb, 1.0f, (gsl_vector_float *)vc); 
    case PDP_MATRIX_TYPE_RDOUBLE:
	return gsl_blas_dgemv(TransA, (double)scale_r, (gsl_matrix *)mA, 
			      (gsl_vector *)vb, 1.0, (gsl_vector *)vc);
    case PDP_MATRIX_TYPE_CFLOAT:
	return gsl_blas_cgemv(TransA, cf_scale, (gsl_matrix_complex_float *)mA, 
			      (gsl_vector_complex_float *)vb, cf_one, (gsl_vector_complex_float *)vc); 
    case PDP_MATRIX_TYPE_CDOUBLE:
	return gsl_blas_zgemv(TransA, cd_scale, (gsl_matrix_complex *)mA, 
			      (gsl_vector_complex *)vb, cd_one, (gsl_vector_complex *)vc); 
    default:
	return 0;
    }
}



t_pdp_symbol *_pdp_matrix_get_description(t_pdp *header)
{
    t_matrix *m = (t_matrix *)(&header->info.raw);
    char description[100];
    char *c = description;
    int encoding;

    if (!header) return pdp_gensym("invalid");
    else if (!header->desc){
	/* if description is not defined, try to reconstruct it (for backwards compat) */
	if (header->type == PDP_MATRIX){

	    c += sprintf(c, "matrix");

	    switch(m->type){
	    case PDP_MATRIX_TYPE_RFLOAT: c += sprintf(c, "/float/real"); break;
	    case PDP_MATRIX_TYPE_CFLOAT: c += sprintf(c, "/float/complex"); break;
	    case PDP_MATRIX_TYPE_RDOUBLE: c += sprintf(c, "/double/real"); break;
	    case PDP_MATRIX_TYPE_CDOUBLE: c += sprintf(c, "/double/complex"); break;
	    default:
		c += sprintf(c, "/unknown"); goto exit;
	    }

	    c += sprintf(c, "/%dx%d", (int)m->rows, (int)m->columns);
	    
	exit:
	    return pdp_gensym(description);
	} 
	else return pdp_gensym("unknown");
    }
    else return header->desc;
}



static void _pdp_matrix_copy(t_pdp *dst, t_pdp *src);
static void _pdp_matrix_clone(t_pdp *dst, t_pdp *src);
static void _pdp_matrix_reinit(t_pdp *dst);
static void _pdp_matrix_cleanup(t_pdp *dst);

static size_t _pdp_matrix_mdata_byte_size(u32 rows, u32 columns, u32 type)
{
    size_t dsize = rows * columns;
    switch (type){
    case PDP_MATRIX_TYPE_RFLOAT:   dsize *= sizeof(float);    break;
    case PDP_MATRIX_TYPE_CFLOAT:   dsize *= 2*sizeof(float);  break;
    case PDP_MATRIX_TYPE_RDOUBLE:  dsize *= sizeof(double);   break;
    case PDP_MATRIX_TYPE_CDOUBLE:  dsize *= 2*sizeof(double); break;
    default: dsize = 0;
    }
    return dsize;
}


static size_t _pdp_matrix_pdata_vector_size(u32 rows, u32 columns)
{
    return (rows>columns ? rows : columns);
}


static void _pdp_matrix_init(t_pdp *dst, u32 rows, u32 columns, u32 type)
{
    int i;
    t_matrix *m = (t_matrix *)(&dst->info.raw);
    void *d = ((void *)dst) + PDP_HEADER_SIZE;
    int matrix_bytesize = _pdp_matrix_mdata_byte_size(rows, columns, type);
    int permsize =  _pdp_matrix_pdata_vector_size(rows, columns);

    /* set meta data */
    m->type = type;
    m->rows = rows;
    m->columns = columns;

    /* set the block data */
    m->block.size = matrix_bytesize;
    m->block.data = (double *)d;

    /* set the vector view */
    m->vector.size = (1==columns) ? rows : columns;
    m->vector.stride = 1;
    m->vector.data = (double *)d;
    m->vector.block = &m->block;
    m->vector.owner = 0;

    /* set the matrix view */
    m->matrix.gsl_rows = rows;
    m->matrix.gsl_columns = columns;
    m->matrix.tda = columns;
    m->matrix.data = (double *)d;
    m->matrix.block = &m->block;
    m->matrix.owner = 0;

    /* set the permutation object & init */
    m->perm.size = permsize;
    m->perm.data = (size_t *)(d + matrix_bytesize);
    for(i=0; i<permsize; i++) m->perm.data[i] = i;
    m->signum = -1;
    
    /* init packet header */
    dst->theclass = matrix_class;
    dst->desc = _pdp_matrix_get_description(dst);

}

static void _pdp_matrix_clone(t_pdp *dst, t_pdp *src)
{
    t_matrix *m = (t_matrix *)(&src->info.raw);
    _pdp_matrix_init(dst, m->rows, m->columns, m->type);

}
static void _pdp_matrix_copy(t_pdp *dst, t_pdp *src)
{
    _pdp_matrix_clone(dst, src);
    memcpy(dst + PDP_HEADER_SIZE, src + PDP_HEADER_SIZE, src->size - PDP_HEADER_SIZE);
    //post("matrix copy successful");
}
static void _pdp_matrix_reinit(t_pdp *dst)
{
    /* nothing to do, assuming data is correct */
}
static void _pdp_matrix_cleanup(t_pdp *dst)
{
    /* no extra memory to free */
}

int pdp_packet_new_matrix(u32 rows, u32 columns, u32 type)
{
    t_pdp *header;
    int dsize, p;

    /* compute the blocksize for the matrix data */
    /* if 0, something went wrong -> return invalid packet */
    if (!(dsize = _pdp_matrix_mdata_byte_size(rows, columns, type))) return -1;
    dsize += sizeof(size_t) * _pdp_matrix_pdata_vector_size(rows, columns);

    p = pdp_packet_new(PDP_MATRIX, dsize);
    if (-1 == p) return -1;
    header = pdp_packet_header(p);
    
    _pdp_matrix_init(header, rows, columns, type);

    return p;
}

int pdp_packet_new_matrix_product_result(CBLAS_TRANSPOSE_t TransA, CBLAS_TRANSPOSE_t TransB, int pA, int pB)
{
    gsl_matrix *mA, *mB;
    u32 type = pdp_packet_matrix_get_type(pA);

    u32 colA, colB, rowA, rowB;
   
    /* check if A is a matrix */
    if (!type) return -1;

    mA = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pA, type);
    mB = (gsl_matrix *)pdp_packet_matrix_get_gsl_matrix(pB, type);

    /* check if A and B are same type */
    if (!mB) return -1;

    /* get dims A */
    if (TransA == CblasNoTrans){
	rowA = mA->gsl_rows;
	colA = mA->gsl_columns;
    }
    else {
	rowA = mA->gsl_columns;
	colA = mA->gsl_rows;
    }

    /* get dims B */
    if (TransB == CblasNoTrans){
	rowB = mB->gsl_rows;
	colB = mB->gsl_columns;
    }
    else {
	rowB = mB->gsl_columns;
	colB = mB->gsl_rows;
    }

    /* check if sizes are compatible */
    if (colA != rowB) return -1;

    /* create new packet */
    return pdp_packet_new_matrix(rowA, colB, type);
}

void pdp_packet_matrix_setzero(int p)
{
    t_matrix *m;
    if (!pdp_packet_matrix_isvalid(p)) return;
    m = (t_matrix *) pdp_packet_subheader(p);
    memset(m->block.data, 0, m->block.size);
}



/* type conversion programs */
static int _pdp_packet_convert_matrix_to_greyimage(int packet, t_pdp_symbol *dest_template)
{
    t_pdp *header = pdp_packet_header(packet);
    t_matrix *matrix = (t_matrix *)pdp_packet_subheader(packet);
    void *data = pdp_packet_data(packet);
    s16 *new_data;
    u32 c,r, nbelements;
    int new_p;
    u32 i;

    c = matrix->columns;
    r = matrix->rows;
    nbelements = c*r;

    new_p = pdp_packet_new_image_grey(c,r);
    if (-1 == new_p) return -1;
    new_data = (s16 *)pdp_packet_data(new_p);

    /* convert first channel */
    switch(matrix->type){
    case PDP_MATRIX_TYPE_RFLOAT:
	for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((float *)data)[i] * (float)0x8000); break;
    case PDP_MATRIX_TYPE_CFLOAT: //only copy real channel
	for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((float *)data)[i<<1] * (float)0x8000); break;
    case PDP_MATRIX_TYPE_RDOUBLE:
	for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((double *)data)[i] * (float)0x8000); break;
    case PDP_MATRIX_TYPE_CDOUBLE: //only copy real channel
	for (i=0; i<nbelements; i++) (new_data)[i] = (s16)(((double *)data)[i<<1] * (float)0x8000); break;
    default:
	pdp_packet_mark_unused(new_p);
	new_p = -1;
    }	
    return new_p;
}

static int _pdp_packet_convert_image_to_matrix(int packet, t_pdp_symbol *dest_template, int type)
{
    t_pdp *header = pdp_packet_header(packet);
    t_image *image = pdp_packet_image_info(packet);
    s16 *data = (s16 *)pdp_packet_data(packet);
    void *new_data;
    u32 w,h, nbelements;
    int new_p;
    int encoding = image->encoding;
    u32 i;

    if (!pdp_packet_image_isvalid(packet)) return -1;
    w = image->width;
    h = image->height;
    nbelements = w*h;

    new_p = pdp_packet_new_matrix(h,w, type);
    if (-1 == new_p) return -1;
    new_data = pdp_packet_data(new_p);

    switch (encoding){
    case PDP_IMAGE_YV12:
    case PDP_IMAGE_GREY:
    case PDP_IMAGE_MCHP:
	/* convert first channel */
	switch(type){
	case PDP_MATRIX_TYPE_RFLOAT:
	    for (i=0; i<nbelements; i++)
		((float *)new_data)[i] = ((float)data[i]) * (1.0f / (float)0x8000); break;
	case PDP_MATRIX_TYPE_RDOUBLE:
	    for (i=0; i<nbelements; i++)
		((double *)new_data)[i] = ((double)data[i]) * (1.0f / (double)0x8000); break;
	case PDP_MATRIX_TYPE_CFLOAT:
	    for (i=0; i<nbelements; i++){
		((float *)new_data)[i*2] = ((float)data[i]) * (1.0f / (float)0x8000); 
		((float *)new_data)[i*2+1] = 0.0f;
	    }
	    break;
	case PDP_MATRIX_TYPE_CDOUBLE:
	    for (i=0; i<nbelements; i++){
		((double *)new_data)[i*2] = ((double)data[i]) * (1.0f / (double)0x8000); 
		((double *)new_data)[i*2+1] = 0.0;
	    }
	    break;
	default:
	    pdp_post("_pdp_packet_convert_image_to_matrix: INTERNAL ERROR");
	}
	break;
    default:
	pdp_packet_mark_unused(new_p);
	new_p = -1;
	break;
    }

    return new_p;

}

static int _pdp_packet_convert_image_to_floatmatrix(int packet, t_pdp_symbol *dest_template){
    return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_RFLOAT);}
static int _pdp_packet_convert_image_to_doublematrix(int packet, t_pdp_symbol *dest_template){
    return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_RDOUBLE);}
static int _pdp_packet_convert_image_to_complexfloatmatrix(int packet, t_pdp_symbol *dest_template){
    return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_CFLOAT);}
static int _pdp_packet_convert_image_to_complexdoublematrix(int packet, t_pdp_symbol *dest_template){
    return _pdp_packet_convert_image_to_matrix(packet, dest_template, PDP_MATRIX_TYPE_CDOUBLE);}


static int _pdp_packet_matrix_convert_fallback(int packet, t_pdp_symbol *dest_template)
{
    pdp_post("can't convert image type %s to %s",
	 pdp_packet_get_description(packet)->s_name, dest_template->s_name);

    return -1;
}


/* this seems like a nice spot to place a dummy gsl signal handler */
static void _gsl_error_handler (const char * reason, 
              const char * file, 
              int line, 
              int gsl_errno)
{
    //pdp_post("gsl error:\nREASON: %s\nFILE:%s\nLINE:%d\nERRNO:%d", reason, file, line, gsl_errno);
}

static int pdp_matrix_factory(t_pdp_symbol *type)
{
    return -1;
}
void pdp_matrix_setup(void)
{
    t_pdp_conversion_program *program;


    /* setup the class */
    matrix_class = pdp_class_new(pdp_gensym("matrix/*/*/*"), pdp_matrix_factory);
    matrix_class->copy = _pdp_matrix_copy;
    //matrix_class->clone = _pdp_matrix_clone; // is now solved through (symbol based) constructor
    matrix_class->wakeup = _pdp_matrix_reinit;
    matrix_class->cleanup = _pdp_matrix_cleanup;
    

    /* image -> matrix: default = double/real (most ops available) */
    program = pdp_conversion_program_new(_pdp_packet_convert_image_to_doublematrix, 0);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/double/real/*"), program);
    program = pdp_conversion_program_new(_pdp_packet_convert_image_to_floatmatrix, 0);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/float/real/*"), program);
    program = pdp_conversion_program_new(_pdp_packet_convert_image_to_complexdoublematrix, 0);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/double/complex/*"), program);
    program = pdp_conversion_program_new(_pdp_packet_convert_image_to_complexfloatmatrix, 0);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/float/complex/*"), program);

    /* matrix -> image */
    program = pdp_conversion_program_new(_pdp_packet_convert_matrix_to_greyimage, 0);
    pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("image/grey/*"), program);

    /* fallbacks */
    program = pdp_conversion_program_new(_pdp_packet_matrix_convert_fallback, 0);
    pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("image/*/*"), program);
    pdp_type_register_conversion(pdp_gensym("matrix/*/*/*"), pdp_gensym("bitmap/*/*"), program);
    pdp_type_register_conversion(pdp_gensym("image/*/*"), pdp_gensym("matrix/*/*/*/*"), program);
    pdp_type_register_conversion(pdp_gensym("bitmap/*/*"), pdp_gensym("matrix/*/*/*/*"), program);

    /* setup gsl handler */
    if(gsl_set_error_handler(_gsl_error_handler)){
	pdp_post("pdp_matrix_setup: WARNING: overriding gsl error handler.");
    }    

}





More information about the Pd-cvs mailing list