From 05da8712f62879b20741a9475d4892870b2d24d7 Mon Sep 17 00:00:00 2001 From: Jens Ahrensfeld Date: Sat, 19 Jul 2014 07:44:42 +0000 Subject: [PATCH] Initial import git-svn-id: http://moon:8086/svn/software/trunk/libsrc/_to_be_added@1 b431acfa-c32f-4a4a-93f1-934dc6c82436 --- clob.c | 418 +++++++++++++++ clob.h | 64 +++ color_conv.c | 575 ++++++++++++++++++++ color_conv.h | 41 ++ imageio.c | 225 ++++++++ imageio.h | 30 ++ lapack_cwrap.c | 350 ++++++++++++ lapack_cwrap.h | 1389 ++++++++++++++++++++++++++++++++++++++++++++++++ matte.c | 977 ++++++++++++++++++++++++++++++++++ matte.h | 66 +++ mlio.c | 82 +++ mlio.h | 46 ++ mtypes.h | 34 ++ mutil.c | 220 ++++++++ mutil.h | 26 + pixel.c | 413 ++++++++++++++ pixel.h | 55 ++ tictoc.c | 29 + tictoc.h | 24 + 19 files changed, 5064 insertions(+) create mode 100755 clob.c create mode 100755 clob.h create mode 100755 color_conv.c create mode 100755 color_conv.h create mode 100755 imageio.c create mode 100755 imageio.h create mode 100755 lapack_cwrap.c create mode 100755 lapack_cwrap.h create mode 100755 matte.c create mode 100755 matte.h create mode 100755 mlio.c create mode 100755 mlio.h create mode 100755 mtypes.h create mode 100755 mutil.c create mode 100755 mutil.h create mode 100755 pixel.c create mode 100755 pixel.h create mode 100755 tictoc.c create mode 100755 tictoc.h diff --git a/clob.c b/clob.c new file mode 100755 index 0000000..b72ba04 --- /dev/null +++ b/clob.c @@ -0,0 +1,418 @@ +// ------------------------------------------------------------ +// clob.c +// +// A color quantizer using binary-split technique +// From: +// "Color Quantization of Images", Orchard, Bouman +// IEEE Trans. on Sig. Proc., vol. 39, no. 12, pp. 2677-2690, Dec. 1991. +// +// 10.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#include +#include +#include +#include +#include +#include "bintree.h" +#include "mutil.h" +#include "lapack_cwrap.h" +#include "clob.h" + +// ------------------------------------------------------------ +void CodebookAlloc(codebook_t *pObj, int N_max, int M_max) +{ + int i; + + pObj->pColor = (color_vector_t*)malloc(M_max*sizeof(color_vector_t)); + pObj->size_C = (int*)malloc(M_max*sizeof(int)); + memset(pObj->size_C, 0, M_max*sizeof(int)); + pObj->pC = (int*)malloc(N_max*sizeof(int)); + pObj->ppC = (int**)malloc(M_max*sizeof(int*)); + + pObj->ppC[0] = pObj->pC;; + for (i=1; i < M_max; i++) + pObj->ppC[i] = NULL; + + pObj->max_colors = M_max; + pObj->max_samples = N_max; + pObj->num_colors = 0; + +} + +// ------------------------------------------------------------ +void CodebookFree(codebook_t *pObj) +{ + int i; + + if(!pObj) + return; + + free(pObj->pC); + free(pObj->ppC); + free(pObj->pColor); + free(pObj->size_C); +} + +// ------------------------------------------------------------ +void ColorQuant_Init(colorquant_t *pCQ, int N_max, int M_max) +{ + M_max++; + pCQ->N_max = N_max; + pCQ->M_max = M_max; + + pCQ->pRootData = (node_data_t*)malloc(sizeof(node_data_t)); + memset(pCQ->pRootData, 0, sizeof(node_data_t)); + + pCQ->pData = (node_data_t*)malloc(2*M_max*sizeof(node_data_t)); + memset(pCQ->pData, 0, 2*M_max*sizeof(node_data_t)); + + pCQ->pTemp = (double*)malloc(N_max*sizeof(double)); + pCQ->pXs = (color_vector_t*)malloc(N_max*sizeof(color_vector_t)); + +} + +// ------------------------------------------------------------ +void ColorQuant_Free(colorquant_t *pCQ) +{ + int m; + + FreeNodeData(pCQ->pRootData); + for (m=0; m < 2*pCQ->M_max; m++) + { + FreeNodeData(&pCQ->pData[m]); + } + free(pCQ->pRootData); + pCQ->pRootData = NULL; + free(pCQ->pData); + pCQ->pData = NULL; + free(pCQ->pXs); + free(pCQ->pTemp); +} + +// ------------------------------------------------------------ +void FreeNodeData(node_data_t *pObj) +{ + if(pObj) + { + if(pObj->pC) + { + free(pObj->pC); + pObj->pC = NULL; + } + if(pObj->pExpr_left) + { + free(pObj->pExpr_left); + pObj->pExpr_left = NULL; + } + } +} + +// ------------------------------------------------------------ +void ColorNodeStat(colorquant_t *pCQ, node_data_t *pObj, color_vector_t *pXs) +{ + int i, j, s, info; + + if (pObj->N == 0) + { + fprintf(stderr, "ColorNodeStat(): N is zero!\n"); + return; + } + + // Calc node statistics + // R = sum_s(xs_s*xs_s') ; mit s = Element aus Menge C + memset(pObj->R, 0, ROW_DIM*ROW_DIM*sizeof(double)); + dsyrk('U', 'N', ROW_DIM, pObj->N, 1.0, (double*)pXs, ROW_DIM, 1.0, (double*)pObj->R, ROW_DIM); +// PrintMatrix("R=",(double*)pObj->R, ROW_DIM, ROW_DIM); + + // m = sum_s(xs_s) + memset(pObj->m, 0, ROW_DIM*sizeof(double)); + for(i=0; i < pObj->N; i++) + { + pObj->m[0] += pXs[i][0]; + pObj->m[1] += pXs[i][1]; + pObj->m[2] += pXs[i][2]; + + } +// PrintMatrix("m=",(double*)pObj->m, ROW_DIM, 1); + +} + +// ------------------------------------------------------------ +double ColorLambda(colorquant_t *pCQ, node_data_t *pObj, color_vector_t *pXs) +{ + + int i, s, info, lwork; + color_vector_t d = {0}; + double *pWork; + double *pDiff, Ls; + + if (pObj->N == 0) + { + fprintf(stderr, "ColorLambda(): N is zero!\n"); + pObj->lambda = 0; + return 0; + } + + pObj->pExpr_left = (double*)malloc(pObj->N*sizeof(double)); + + // q = m/N; + memcpy(pObj->q, pObj->m, sizeof(color_vector_t)); + dscal (ROW_DIM, 1.0/pObj->N, pObj->q, 1); +// PrintMatrix("q=",(double*)pObj->q, ROW_DIM, 1); + + // cov = R - 1/N*m*m' + memcpy(pObj->cov, pObj->R, ROW_DIM*sizeof(color_vector_t)); + dsyr('U', ROW_DIM, -1.0/pObj->N, pObj->m, 1, (double*)pObj->cov, ROW_DIM); +// PrintMatrix("Cov=",(double*)pObj->cov, ROW_DIM, ROW_DIM); + + // Query dsyev() for optimal work buffer size + pWork = pCQ->pTemp; + pWork[0] = 0; + info = dsyev('V', 'U', ROW_DIM, (double*)NULL, ROW_DIM, NULL, pWork, -1); + lwork = (int)pWork[0]; + + // d = eig(cov) + info = dsyev('V', 'U', ROW_DIM, (double*)pObj->cov, ROW_DIM, d, pWork, lwork); + if (info != 0) + { + fprintf (stderr, "Eigenvalue solver failure with error %d\n", info); + printf("N=%d\n", pObj->N); + PrintMatrix("Cov=",(double*)pObj->cov, ROW_DIM, ROW_DIM); + } + +// PrintMatrix("d=",(double*)d, ROW_DIM, 1); +// PrintMatrix("V=",(double*)pObj->cov, ROW_DIM, ROW_DIM); + + // Form unit-vector e in direction of max. element in d + dcopy (ROW_DIM, d, 1, pObj->e, 1); + GetPrincipleEigenVector((double*)pObj->cov, (double*)pObj->e, ROW_DIM); +// PrintMatrix("e=",(double*)pObj->e, ROW_DIM, 1); + + // ------------------------- + // Compute Lambda + + // Expr. left: dot(xs,e) = e'*xs_s = xs_s'*e + // Also used later in ColorSplit() + dgemv('T', ROW_DIM, pObj->N, 1.0, (double*)pXs, ROW_DIM, pObj->e, 1, 0.0, pObj->pExpr_left, 1); + + // Expr. right: dot(q,e) = e'*q = q'*e + // Also used later in ColorSplit() + pObj->expr_right = ddot(ROW_DIM, pObj->q, 1, pObj->e, 1); + + // lambda = sum_s((xs_s - q)'*e).^2); + pObj->lambda = 0; + for (i=0; i < pObj->N; i++) + { + // xs_s'*e - q'*e + Ls = pObj->pExpr_left[i] - pObj->expr_right; + pObj->lambda += Ls*Ls; + } + return pObj->lambda; +} + +// ------------------------------------------------------------ +void ColorNodeSplit(colorquant_t *pCQ, node_data_t *pObj, node_data_t *pData2n, node_data_t *pData2n1) +{ + int N2n, N2n1, s, i; + + memset(pData2n, 0, sizeof(node_data_t)); + memset(pData2n1, 0, sizeof(node_data_t)); + + pData2n->pC = (int*)malloc(pObj->N*sizeof(int)); + if(!pData2n->pC) + { + fprintf(stderr, "pData2n->pC = NULL for N=%d\n", pObj->N); + } + + pData2n1->pC = (int*)malloc(pObj->N*sizeof(int)); + if(!pData2n1->pC) + { + fprintf(stderr, "pData2n1->pC = NULL for N=%d\n", pObj->N); + } + + // Split into two new sets + N2n = 0; + N2n1 = 0; + for (i=0; i < pObj->N; i++) + { + s = pObj->pC[i]; + + if(pObj->pExpr_left[i] > pObj->expr_right) + { + pData2n1->pC[N2n1] = s; + N2n1++; + } + else + { + pData2n->pC[N2n] = s; + N2n++; + } + } + + // Compute statistics of two new nodes + // Node(2n) + pData2n->N = N2n; + for (i=0; i < N2n; i++) + { + s = pData2n->pC[i]; + pCQ->pXs[i][0] = pCQ->pX[s][0]; + pCQ->pXs[i][1] = pCQ->pX[s][1]; + pCQ->pXs[i][2] = pCQ->pX[s][2]; + } + + ColorNodeStat(pCQ, pData2n, pCQ->pXs); + ColorLambda(pCQ, pData2n, pCQ->pXs); + + // Node(2n+1) + pData2n1->N = N2n1; + for(i=0; i < ROW_DIM; i++) + { + pData2n1->R[i][0] = pObj->R[i][0] - pData2n->R[i][0]; + pData2n1->R[i][1] = pObj->R[i][1] - pData2n->R[i][1]; + pData2n1->R[i][2] = pObj->R[i][2] - pData2n->R[i][2]; + } + + for(i=0; i < ROW_DIM; i++) + pData2n1->m[i] = pObj->m[i] - pData2n->m[i]; + + for (i=0; i < N2n1; i++) + { + s = pData2n1->pC[i]; + pCQ->pXs[i][0] = pCQ->pX[s][0]; + pCQ->pXs[i][1] = pCQ->pX[s][1]; + pCQ->pXs[i][2] = pCQ->pX[s][2]; + } + ColorLambda(pCQ, pData2n1, pCQ->pXs); + +} + +// ------------------------------------------------------------ +int ColorQuant(colorquant_t *pCQ, codebook_t *pCB, color_vector_t *pX, int K, int N, int M, double lambda_min) +{ + int i, m, *pC, bst_state, s, sum_pixel, codebook_size, size, *ptr, *ptr_next; + node_data_t *pData, *pData2n, *pData2n1; + node_t *node_max, *node_temp; + bintree_t bst; + + if (N == 0) + { + fprintf(stderr, "ColorQuant(): N is zero!\n"); + return 0; + } + + pCQ->nRow = K; + pCQ->nCol = N; + + pCQ->pRootData->N = N; + pCQ->pX = (color_vector_t*)pX; + + pCQ->pRootData->pC = (int*)malloc(N*sizeof(int)); + + // Create root data set + for(i=0; i < N; i++) + { + pCQ->pRootData->pC[i] = i; // store indices of vector boundaries + } + + ColorNodeStat(pCQ, pCQ->pRootData, pCQ->pX); + ColorLambda(pCQ, pCQ->pRootData, pCQ->pX); + + // Create root-node + bintree_init(&bst, 2*N); + + bintree_insert(&bst, pCQ->pRootData->lambda, pCQ->pRootData); + + for (m=0; m < M-1; m++) + { + node_max = bintree_GetNodeMax(&bst); + if (!node_max) + { + fprintf(stderr,"node_max = NULL!\n"); + break; + } + + if (node_max->value < lambda_min) + break; + + pData2n = &pCQ->pData[2*m]; + pData2n1 = &pCQ->pData[2*m+1]; + ColorNodeSplit(pCQ, (node_data_t*)node_max->pData, pData2n, pData2n1); + if (pData2n->N > 0) + { + bintree_insert(&bst, pData2n->lambda, pData2n); + } + else + { + FreeNodeData(pData2n); + } + + if (pData2n1->lambda > 0) + { + bintree_insert(&bst, pData2n1->lambda, pData2n1); + } + else + { + FreeNodeData(pData2n1); + } + +// if ((pData2n->lambda == 0) || (pData2n1->lambda == 0)) +// { +// continue; +// } + + node_temp = bintree_node_remove(&bst, node_max); + if (node_temp == bst.pRoot) + { + break; + } + FreeNodeData((node_data_t*)node_max->pData); + + if (!bst.pRoot) + break; + + bst_state = bintree_isBST(bst.pRoot); + if (!bst_state) + { + fprintf(stderr,"Binary Tree is not BST!\n"); + } + } + +// bintree_print(pCQ->root); +// printf("\n"); + sum_pixel = 0; + codebook_size = 0; + pCB->ppC[0] = pCB->pC; + do + { + node_max = bintree_GetNodeMax(&bst); + if (!node_max) + { + fprintf(stderr,"node_max = NULL!\n"); + } +// printf("Lambda=%5.5g\n",node_max->value); + + size = ((node_data_t*)node_max->pData)->N; + pCB->size_C[codebook_size] = size; + ptr = pCB->ppC[codebook_size]; + ptr_next = ptr + size; + sum_pixel += size; + memcpy(ptr,((node_data_t*)node_max->pData)->pC, size*sizeof(int)); + memcpy(pCB->pColor[codebook_size],((node_data_t*)(node_max->pData))->q, sizeof(color_vector_t)); + if ((codebook_size+1) < pCB->max_colors) + pCB->ppC[codebook_size+1] = (int*)(ptr_next); + + FreeNodeData((node_data_t*)(node_max->pData)); + node_temp = bintree_node_remove(&bst, node_max); + codebook_size++; + + } while(node_temp); + pCB->num_colors = codebook_size; + +// printf("Num. pixels = %d\n",sum_pixel); +// printf("Codebook size = %d\n",codebook_size); +// PrintMatrix("Codebook=",(double*)pCB->pColor, ROW_DIM, codebook_size); + + FreeNodeData(pCQ->pRootData); + bintree_free(&bst); + return codebook_size; +} diff --git a/clob.h b/clob.h new file mode 100755 index 0000000..e23911b --- /dev/null +++ b/clob.h @@ -0,0 +1,64 @@ +// ------------------------------------------------------------ +// clob.h +// +// A color quantizer using binary-split technique +// From: +// "Color Quantization of Images", Orchard, Bouman +// IEEE Trans. on Sig. Proc., vol. 39, no. 12, pp. 2677-2690, Dec. 1991. +// +// 10.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef CLOB_H +#define CLOB_H + +#include "bintree.h" +#include "colortypes.h" + +// ------------------------------------------------------------ +#define ROW_DIM 3 + +// ------------------------------------------------------------ +typedef struct _node_data_t +{ + int N, *pC, len_C; + color_vector_t R[ROW_DIM]; + color_vector_t cov[ROW_DIM]; + color_vector_t m; + color_vector_t q; + color_vector_t e; + double lambda; + double *pExpr_left, expr_right; + +} node_data_t; + +typedef struct _colorquant_t +{ + int nRow, nCol, N_max, M_max, nPixel; + double *pTemp; + color_vector_t *pX, *pXs; + node_data_t *pRootData, *pData; + +} colorquant_t; + +typedef struct _codebook_t +{ + int num_colors, max_colors, max_samples; + + color_vector_t *pColor; + int **ppC, *pC, *size_C; + +} codebook_t; + +// ------------------------------------------------------------ +void ColorQuant_Init(colorquant_t *pCQ, int N_max, int M_max); +void ColorQuant_Free(colorquant_t *pCQ); +void FreeNodeData(node_data_t *pObj); +void ColorNodeStat(colorquant_t *pCQ, node_data_t *pObj, color_vector_t *pXs); +double ColorLambda(colorquant_t *pCQ, node_data_t *pObj, color_vector_t *pXs); +void ColorNodeSplit(colorquant_t *pCQ, node_data_t *pObj, node_data_t *pData2n, node_data_t *pData2n1); +int ColorQuant(colorquant_t *pCQ, codebook_t *pCB, color_vector_t *pX, int K, int N, int M, double lambda_min); +void CodebookAlloc(codebook_t *pObj, int N_max, int M_max); +void CodebookFree(codebook_t *pObj); + +// ------------------------------------------------------------ +#endif // CLOB_H diff --git a/color_conv.c b/color_conv.c new file mode 100755 index 0000000..373d627 --- /dev/null +++ b/color_conv.c @@ -0,0 +1,575 @@ +// ---------------------------------------------------------------------- +// color_conv.c +// 21.02.2005, Jens Ahrensfeld +// ---------------------------------------------------------------------- +#include +#include +#include "color_conv.h" +#include "colortypes.h" +#include "mutil.h" + +// ---------------------------------------------------------------------- +void rgb2hsv(color_vector_t *pRGB, color_vector_t *pHSV, int len) +{ + int i; + double *rgb, *hsv, v_min, v_max, diff_max, diff_max_r, diff_r, diff_g, diff_b; + + for (i=0; i < len; i++) + { + rgb = (double*)&pRGB[i]; + hsv = (double*)&pHSV[i]; + + v_min = smin(rgb[0], rgb[1]); + v_min = smin(v_min, rgb[2]); + + v_max = smax(rgb[0], rgb[1]); + v_max = smax(v_max, rgb[2]); + + diff_max = v_max - v_min; + + hsv[2] = v_max; + + if (diff_max == 0) + { + hsv[0] = diff_max; + hsv[1] = diff_max; + } + else + { + diff_max_r = 1.0/diff_max; + hsv[1] = diff_max/v_max; + + diff_r = (((v_max - rgb[0])*(1.0/6)) + (0.5*diff_max))*diff_max_r; + diff_g = (((v_max - rgb[1])*(1.0/6)) + (0.5*diff_max))*diff_max_r; + diff_b = (((v_max - rgb[2])*(1.0/6)) + (0.5*diff_max))*diff_max_r; + + if (rgb[0] == v_max) + { + hsv[0] = diff_b - diff_g; + } + else + { + if (rgb[1] == v_max) + { + hsv[0] = 1.0/3 + diff_r - diff_b; + } + else + { + if (rgb[2] == v_max) + { + hsv[0] = 2.0/3 + diff_g - diff_r; + } + } + } + + if (hsv[0] < 0) + hsv[0] += 1; + + if (hsv[0] > 1.0) + hsv[0] -= 1; + + } + } +} + +// ---------------------------------------------------------------------- +void hsv2rgb(color_vector_t *pHSV, color_vector_t *pRGB, int len) +{ + int i, j; + double h, t1, t2, t3, *rgb, *hsv; + + for (i=0; i < len; i++) + { + rgb = (double*)&pRGB[i]; + hsv = (double*)&pHSV[i]; + + if (hsv[1] == 0) + { + rgb[0] = hsv[2]; + rgb[1] = hsv[2]; + rgb[2] = hsv[2]; + + continue; + } + + h = hsv[0] * 6; + j = (int)h; + t1 = hsv[2]*(1 - hsv[1]); + t2 = hsv[2]*(1 - hsv[1]*(h-(double)j)); + t3 = hsv[2]*(1 - hsv[1]*(1-(h-(double)j))); + + switch(j) + { + case 0: + rgb[0] = hsv[2]; + rgb[1] = t3; + rgb[2] = t1; + break; + + case 1: + rgb[0] = t2; + rgb[1] = hsv[2]; + rgb[2] = t1; + break; + + case 2: + rgb[0] = t1; + rgb[1] = hsv[2]; + rgb[2] = t3; + break; + + case 3: + rgb[0] = t1; + rgb[1] = t2; + rgb[2] = hsv[2]; + break; + + case 4: + rgb[0] = t3; + rgb[1] = t1; + rgb[2] = hsv[2]; + break; + + default: + rgb[0] = hsv[2]; + rgb[1] = t1; + rgb[2] = t2; + break; + } + } +} + +// ---------------------------------------------------------------------- +void rgb2hsl(color_vector_t *pRGB, color_vector_t *pHSL, int len) +{ + int i; + double *rgb, *hsl, v_min, v_max, diff_max, diff_max_r, diff_r, diff_g, diff_b; + + for (i=0; i < len; i++) + { + rgb = (double*)&pRGB[i]; + hsl = (double*)&pHSL[i]; + + v_min = smin(rgb[0], rgb[1]); + v_min = smin(v_min, rgb[2]); + + v_max = smax(rgb[0], rgb[1]); + v_max = smax(v_max, rgb[2]); + + diff_max = v_max - v_min; + + hsl[2] = (v_max + v_min)/2; + + if (diff_max == 0) + { + hsl[0] = diff_max; + hsl[1] = diff_max; + } + else + { + diff_max_r = 1.0/diff_max; + + if(hsl[2] < 0.5) + { + hsl[1] = diff_max/(v_max + v_min); + } + else + { + hsl[1] = diff_max/(2 - v_max - v_min); + } + + diff_r = (((v_max - rgb[0])*(1.0/6)) + (0.5*diff_max))*diff_max_r; + diff_g = (((v_max - rgb[1])*(1.0/6)) + (0.5*diff_max))*diff_max_r; + diff_b = (((v_max - rgb[2])*(1.0/6)) + (0.5*diff_max))*diff_max_r; + + if (rgb[0] == v_max) + { + hsl[0] = diff_b - diff_g; + } + else + { + if (rgb[1] == v_max) + { + hsl[0] = 1.0/3 + diff_r - diff_b; + } + else + { + if (rgb[2] == v_max) + { + hsl[0] = 2.0/3 + diff_g - diff_r; + } + } + } + + if (hsl[0] < 0) + hsl[0] += 1; + + if (hsl[0] > 1.0) + hsl[0] -= 1; + + } + } +} + +// ---------------------------------------------------------------------- +inline double Hue_2_RGB(double v1, double v2, double vH) //Function Hue_2_RGB +{ + if (vH < 0) + vH += 1; + if (vH > 1) + vH -= 1; + + if ((6*vH) < 1) + return v1 + (v2 - v1)*6*vH; + + if ((2*vH) < 1) + return v2; + + if ((3*vH) < 2) + return v1 + (v2 - v1)*((2.0/3) - vH)*6; + + return v1; +} + +// ---------------------------------------------------------------------- +void hsl2rgb(color_vector_t *pHSL, color_vector_t *pRGB, int len) +{ + int i, j; + double h, t1, t2, t3, *rgb, *hsl; + + for (i=0; i < len; i++) + { + rgb = (double*)&pRGB[i]; + hsl = (double*)&pHSL[i]; + + if (hsl[1] == 0) + { + rgb[0] = hsl[2]; + rgb[1] = hsl[2]; + rgb[2] = hsl[2]; + + continue; + } + + if(hsl[2] < 0.5) + { + t2 = hsl[2]*(1 + hsl[1]); + } + else + { + t2 = (hsl[2] + hsl[1]) - hsl[2]*hsl[1]; + } + + t1 = 2*hsl[2] - t2; + + rgb[0] = Hue_2_RGB(t1, t2, hsl[0] + 1.0/3); + rgb[1] = Hue_2_RGB(t1, t2, hsl[0]); + rgb[2] = Hue_2_RGB(t1, t2, hsl[0] - 1.0/3); + } +} + +// ---------------------------------------------------------------------- +void rgb2xyz(double **ppRGB, double **ppXYZ, int len) +{ + int i, j; + double rgb[3]; + double *pRGB, *pXYZ; + + for (i=0; i < len; i++) + { + pRGB = ((double*)ppRGB[i]); + pXYZ = ((double*)ppXYZ[i]); + + for (j=0; j < 3; j++) + { + if (pRGB[j] > 0.04045) + { + rgb[j] = pow((pRGB[j]+0.055)/1.055,2.4); + } + else + { + rgb[j] = 1.0/12.92 * pRGB[j]; + } + + rgb[j] *= K_XYZ; + + } + + // Observer = 2°, Illuminant = D65 + pXYZ[0] = 0.412453*rgb[0] + 0.357580*rgb[1] + 0.180423*rgb[2]; + pXYZ[1] = 0.212671*rgb[0] + 0.715160*rgb[1] + 0.072169*rgb[2]; + pXYZ[2] = 0.019334*rgb[0] + 0.119193*rgb[1] + 0.950227*rgb[2]; + + } +} + +// ---------------------------------------------------------------------- +void xyz2rgb(double **ppXYZ, double **ppRGB, int len) +{ + int i, j; + double xyz[3]; + double *pRGB, *pXYZ; + + for (i=0; i < len; i++) + { + pRGB = ((double*)ppRGB[i]); + pXYZ = ((double*)ppXYZ[i]); + + for (j=0; j < 3; j++) + { + xyz[j] = 1.0/K_XYZ * pXYZ[j]; + } + + // Observer = 2°, Illuminant = D65 + pRGB[0] = 3.240479*xyz[0] - 1.537150*xyz[1] - 0.498535*xyz[2]; + pRGB[1] = -0.969256*xyz[0] + 1.875992*xyz[1] + 0.041556*xyz[2]; + pRGB[2] = 0.055648*xyz[0] - 0.204043*xyz[1] + 1.057311*xyz[2]; + + for (j=0; j < 3; j++) + { + if (pRGB[j] > 0.0031308) + { + pRGB[j] = 1.055 * pow(pRGB[j],1.0/2.4) - 0.055; + } + else + { + pRGB[j] = 12.92 * pRGB[j]; + } + } + } +} + +// ---------------------------------------------------------------------- +void xyz2lab(double **ppXYZ, double **ppLAB, double *pRef, int len) +{ + int i, j; + double xyz[3] = {0}; + double *pLAB, *pXYZ; + + for (i=0; i < len; i++) + { + pXYZ = ((double*)ppXYZ[i]); + pLAB = ((double*)ppLAB[i]); + + for (j=0; j < 3; j++) + { + xyz[j] = pXYZ[j] / pRef[j]; + } + + for (j=0; j < 3; j++) + { + if (xyz[j] > 0.008856) + { + xyz[j] = pow(xyz[j],1.0/3.0); + } + else + { + xyz[j] = 7.787*xyz[j] + 16.0/116.0; + } + + } + + pLAB[0] = 116.0*xyz[1] - 16; + pLAB[1] = 500.0*(xyz[0] - xyz[1]); + pLAB[2] = 200.0*(xyz[1] - xyz[2]); + + } +} + +// ---------------------------------------------------------------------- +void lab2xyz(double **ppLAB, double **ppXYZ, double *pRef, int len) +{ + int i, j; + double xyz[3], xyz3; + double *pLAB, *pXYZ; + + for (i=0; i < len; i++) + { + pXYZ = ((double*)ppXYZ[i]); + pLAB = ((double*)ppLAB[i]); + + xyz[1] = 1.0/116.0*(pLAB[0] + 16); + xyz[0] = 1.0/500.0*pLAB[1] + xyz[1]; + xyz[2] = xyz[1] - 1.0/200.0*pLAB[2]; + + for (j=0; j < 3; j++) + { + xyz3 = pow(xyz[j],3.0); + if (xyz3 > 0.008856) + { + xyz[j] = xyz3; + } + else + { + xyz[j] = 1.0/7.787*(xyz[j] - 16.0/116); + } + } + + for (j=0; j < 3; j++) + { + pXYZ[j] = xyz[j] *pRef[j]; + } + + } +} + +// ---------------------------------------------------------------------- +void tristimulus(double *pRef, unsigned illuminant_type, unsigned observer_type) +{ + double ref[3] = {K_XYZ*0.95047, K_XYZ*1.0, K_XYZ*1.08883}; + + if(observer_type == OBS_TYPE_2DEG) + { + + switch (illuminant_type) + { + case ILLU_TYPE_A: + case ILLU_TYPE_INCANDESCENT: + pRef[0] = K_XYZ*1.09850; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.35585; + break; + + case ILLU_TYPE_C: + pRef[0] = K_XYZ*0.98074; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.18232; + break; + + case ILLU_TYPE_D50: + pRef[0] = K_XYZ*0.96422; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.82521; + break; + + case ILLU_TYPE_D55: + pRef[0] = K_XYZ*0.95682; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.92149; + break; + + case ILLU_TYPE_D65: + case ILLU_TYPE_DAYLIGHT: + pRef[0] = K_XYZ*0.95047; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.08883; + break; + + case ILLU_TYPE_D75: + pRef[0] = K_XYZ*0.94972; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.22638; + break; + + case ILLU_TYPE_F2: + case ILLU_TYPE_FLUORESCENT: + pRef[0] = K_XYZ*0.99187; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.67395; + break; + + case ILLU_TYPE_F7: + pRef[0] = K_XYZ*0.95044; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.08755; + break; + + case ILLU_TYPE_F11: + pRef[0] = K_XYZ*1.00966; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.64370; + break; + + default: + memcpy(pRef, ref, 3*sizeof(double)); + break; + } + } + + if(observer_type == OBS_TYPE_10DEG) + { + + switch (illuminant_type) + { + case ILLU_TYPE_A: + case ILLU_TYPE_INCANDESCENT: + pRef[0] = K_XYZ*1.11144; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.35200; + break; + + case ILLU_TYPE_C: + pRef[0] = K_XYZ*0.97285; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.16145; + break; + + case ILLU_TYPE_D50: + pRef[0] = K_XYZ*0.96720; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.81427; + break; + + case ILLU_TYPE_D55: + pRef[0] = K_XYZ*0.95799; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.90926; + break; + + case ILLU_TYPE_D65: + case ILLU_TYPE_DAYLIGHT: + pRef[0] = K_XYZ*0.94811; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.07304; + break; + + case ILLU_TYPE_D75: + pRef[0] = K_XYZ*0.94416; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.20641; + break; + + case ILLU_TYPE_F2: + case ILLU_TYPE_FLUORESCENT: + pRef[0] = K_XYZ*1.03280; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.69026; + break; + + case ILLU_TYPE_F7: + pRef[0] = K_XYZ*0.95792; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*1.07687; + break; + + case ILLU_TYPE_F11: + pRef[0] = K_XYZ*1.03866; + pRef[1] = K_XYZ*1.0; + pRef[2] = K_XYZ*0.65627; + break; + + default: + memcpy(pRef, ref, 3*sizeof(double)); + break; + } + } +} + +void rgb2rgb(double **ppSrc, double **ppDst, int len) +{ + int i, j; + double *pSrc, *pDst; + + for (i=0; i < len; i++) + { + pSrc = ((double*)ppSrc[i]); + pDst = ((double*)ppDst[i]); + + for (j=0; j < 3; j++) + { + pDst[j] = pSrc[j]; + } + + } +} + diff --git a/color_conv.h b/color_conv.h new file mode 100755 index 0000000..acc01d3 --- /dev/null +++ b/color_conv.h @@ -0,0 +1,41 @@ +// ---------------------------------------------------------------------- +// color_conv.h +// 21.02.2005, Jens Ahrensfeld +// ---------------------------------------------------------------------- +#ifndef COLOR_CONV_H +#define COLOR_CONV_H + +#include "colortypes.h" + +// ---------------------------------------------------------------------- +#define K_XYZ 1.0 + +#define OBS_TYPE_2DEG 2 +#define OBS_TYPE_10DEG 10 +#define ILLU_TYPE_A 1 +#define ILLU_TYPE_C 2 +#define ILLU_TYPE_D50 3 +#define ILLU_TYPE_D55 4 +#define ILLU_TYPE_D65 5 +#define ILLU_TYPE_D75 6 +#define ILLU_TYPE_F2 7 +#define ILLU_TYPE_F7 8 +#define ILLU_TYPE_F11 9 +#define ILLU_TYPE_INCANDESCENT 10 +#define ILLU_TYPE_DAYLIGHT 11 +#define ILLU_TYPE_FLUORESCENT 12 + +// ---------------------------------------------------------------------- +void rgb2hsv(color_vector_t *pRGB, color_vector_t *pHSV, int len); +void hsv2rgb(color_vector_t *pHSV, color_vector_t *pRGB, int len); +void rgb2hsl(color_vector_t *pRGB, color_vector_t *pHSL, int len); +void hsl2rgb(color_vector_t *pHSL, color_vector_t *pRGB, int len); +void rgb2xyz(double **ppRGB, double **ppXYZ, int len); +void xyz2rgb(double **ppXYZ, double **ppRGB, int len); +void xyz2lab(double **ppXYZ, double **ppLAB, double *pRef, int len); +void lab2xyz(double **ppLAB, double **ppXYZ, double *pRef, int len); +void rgb2rgb(double **ppSrc, double **ppDst, int len); +void tristimulus(double *pRef, unsigned illuminant_type, unsigned observer_type); + +// ---------------------------------------------------------------------- +#endif // COLOR_CONV_H diff --git a/imageio.c b/imageio.c new file mode 100755 index 0000000..c92197d --- /dev/null +++ b/imageio.c @@ -0,0 +1,225 @@ +// ------------------------------------------------------------ +// imageio.c +// +// Reading writing tiff images +// +// 12.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#include +#include +#include +#include + +#ifndef WIN32 +#include +else +#include +#endif + +#include +#include "mtypes.h" +#include "imageio.h" +#include "mutil.h" + +// ------------------------------------------------------------ +void ImageInit(image_t *pObj, int num_row, int num_col, int num_ch) +{ + pObj->num_row = num_row; + pObj->num_col = num_col; + pObj->num_ch = num_ch; + pObj->num_pixel = num_row*num_col; + + pObj->pColor_data = (double*)malloc(pObj->num_pixel*num_ch*sizeof(double)); + memset(pObj->pColor_data, 0, pObj->num_pixel*num_ch*sizeof(double)); + +} + +// ------------------------------------------------------------ +void ImageFree(image_t *pObj) +{ + pObj->num_row = 0; + pObj->num_col = 0; + pObj->num_pixel = 0; + pObj->num_ch = 0; + + if (pObj->pColor_data) + + free(pObj->pColor_data); + +} + +// ------------------------------------------------------------ +void ImageCopy(image_t *pSrc, image_t *pDst) +{ + pDst->num_row = pSrc->num_row; + pDst->num_col = pSrc->num_col; + pDst->num_pixel = pSrc->num_pixel; + pDst->num_ch = pSrc->num_ch; + + if (!pDst->pColor_data) + pDst->pColor_data = (double*)malloc(pSrc->num_ch*pSrc->num_pixel*sizeof(double)); + + + memcpy(pDst->pColor_data, pSrc->pColor_data, pSrc->num_ch*pSrc->num_pixel*sizeof(double)); + +} + +// ------------------------------------------------------------ +int ImageLoadTiff(image_t *pObj, char *filename) +{ + + TIFF *pTifComp; + UINT16 config, nBits, nSamples, photometric; + UINT32 i, num_col, num_row, scan_size, size, num_color_ch; + double k_norm; + double c; + + UINT8 *pImage_comp8; + int status; + + pTifComp = TIFFOpen(filename,"r"); + + if (!pTifComp) + { + fprintf(stderr,"Unable to open input file %s\n",filename); + return -1; + } + + status = TIFFGetField(pTifComp, TIFFTAG_IMAGEWIDTH, &num_col); + status = TIFFGetField(pTifComp, TIFFTAG_IMAGELENGTH, &num_row); + status = TIFFGetField(pTifComp, TIFFTAG_BITSPERSAMPLE, &nBits); + status = TIFFGetField(pTifComp, TIFFTAG_PLANARCONFIG, &config); + status = TIFFGetField(pTifComp, TIFFTAG_PHOTOMETRIC, &photometric); + status = TIFFGetField(pTifComp, TIFFTAG_SAMPLESPERPIXEL, &nSamples); + + num_color_ch = (UINT32)nSamples; + pObj->num_row = num_row; + pObj->num_col = num_col; + pObj->num_pixel = num_row*num_col; + pObj->num_ch = num_color_ch; + + pObj->pColor_data = (double*)malloc(num_color_ch*pObj->num_pixel*sizeof(double)); + + scan_size = TIFFScanlineSize(pTifComp); + + pImage_comp8 = (UINT8*)_TIFFmalloc(scan_size*num_row*sizeof(UINT8)); + + for (i=0; i < num_row; i++) + status = TIFFReadScanline(pTifComp, &pImage_comp8[i*scan_size], i, 0); + + TIFFClose(pTifComp); + + k_norm = 1.0/(pow(2.0,nBits)-1); + + for (i=0; i < num_color_ch*num_col*num_row; i++) + { + c = k_norm*(double)pImage_comp8[i]; + pObj->pColor_data[i] = c; + } + + free(pImage_comp8); + +} + +// ------------------------------------------------------------ +int ImageSaveTiff(image_t *pObj, double *pAlpha, char *filename) +{ + + TIFF *pTifComp; + ttag_t config, nBits, nSamples, photometric, extra_sample_info; + UINT32 i, num_col, num_row, scan_size, size; + double k_norm; + double c, a; + + UINT8 *pImage_comp8; + int status; + + pTifComp = TIFFOpen(filename,"w"); + + if (!pTifComp) + { + fprintf(stderr,"Unable to open input file %s\n",filename); + return 1; + } + + nBits = 8; + num_row = pObj->num_row; + num_col = pObj->num_col; + config = PLANARCONFIG_CONTIG; + photometric = PHOTOMETRIC_RGB; + + nSamples = (ttag_t)pObj->num_ch; + + extra_sample_info = 0; + if (pAlpha) + { + nSamples++; + extra_sample_info = EXTRASAMPLE_ASSOCALPHA; + status = TIFFSetField(pTifComp, TIFFTAG_EXTRASAMPLES, 1, &extra_sample_info); + + } + + status = TIFFSetField(pTifComp, TIFFTAG_IMAGEWIDTH, num_col); + status = TIFFSetField(pTifComp, TIFFTAG_IMAGELENGTH, num_row); + status = TIFFSetField(pTifComp, TIFFTAG_ROWSPERSTRIP, num_row); + status = TIFFSetField(pTifComp, TIFFTAG_BITSPERSAMPLE, nBits); + status = TIFFSetField(pTifComp, TIFFTAG_PLANARCONFIG, config); + status = TIFFSetField(pTifComp, TIFFTAG_PHOTOMETRIC, photometric); + status = TIFFSetField(pTifComp, TIFFTAG_SAMPLESPERPIXEL, nSamples); + status = TIFFSetField(pTifComp, TIFFTAG_ARTIST, "JDI, J. Ahrensfeld"); + +//TIFFTAG_IMAGEDESCRIPTION +//TIFFTAG_SOFTWARE 305 /* name & release */ +//TIFFTAG_DATETIME 306 /* creation date and time */ + + scan_size = nSamples*num_col; + + pImage_comp8 = (UINT8*)_TIFFmalloc(scan_size*num_row*sizeof(UINT8)); + + k_norm = pow(2.0,nBits)-1; + + if (pAlpha) + { + for (i=0; i < num_col*num_row; i++) + { + a = smin(smax(pAlpha[i],0),1); + + c = pObj->pColor_data[(nSamples-1)*i+0]; + c = a*smin(smax(c,0),1); + pImage_comp8[nSamples*i+0] = (UINT8)(k_norm*c); + + c = pObj->pColor_data[(nSamples-1)*i+1]; + c = a*smin(smax(c,0),1); + pImage_comp8[nSamples*i+1] = (UINT8)(k_norm*c); + + c = pObj->pColor_data[(nSamples-1)*i+2]; + c = a*smin(smax(c,0),1); + pImage_comp8[nSamples*i+2] = (UINT8)(k_norm*c); + + pImage_comp8[nSamples*i+3] = (UINT8)(k_norm*a); + + } + } + else + { + for (i=0; i < nSamples*num_col*num_row; i++) + { + c = pObj->pColor_data[i]; + c = smin(smax(c,0),1); + pImage_comp8[i] = (UINT8)(k_norm*c); + } + } + + for (i=0; i < num_row; i++) + { + status = TIFFWriteScanline(pTifComp, &pImage_comp8[i*scan_size], i, 0); + if (status < 0) + TIFFError("ImageSaveTiff", "%s", "Hallo") ; + } + + TIFFClose(pTifComp); + free(pImage_comp8); + + return 0; + +} diff --git a/imageio.h b/imageio.h new file mode 100755 index 0000000..b7626ce --- /dev/null +++ b/imageio.h @@ -0,0 +1,30 @@ +// ------------------------------------------------------------ +// imageio.h +// +// Reading writing tiff images +// +// 12.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef IMAGEIO_H +#define IMAGEIO_H + +#include "mtypes.h" + +// ------------------------------------------------------------ +typedef struct _image_t +{ + UINT32 num_row, num_col, num_pixel, num_ch; + double *pColor_data; + +} image_t; + + +// ------------------------------------------------------------ +void ImageInit(image_t *pObj, int num_row, int num_col, int num_ch); +void ImageFree(image_t *pObj); +void ImageCopy(image_t *pSrc, image_t *pDst); +int ImageLoadTiff(image_t *pObj, char *filename); +int ImageSaveTiff(image_t *pObj, double *pAlpha, char *filename); + +// ------------------------------------------------------------ +#endif // IMAGEIO_H diff --git a/lapack_cwrap.c b/lapack_cwrap.c new file mode 100755 index 0000000..160f4b0 --- /dev/null +++ b/lapack_cwrap.c @@ -0,0 +1,350 @@ +// ------------------------------------------------------------ +// lapack_cwrap.c +// Wrapping LAPACK FORTRAN-function to C-Functions +// Libraries needed: +// lapack.a blas.a and libg2c.a (if available) +// +// 26.02.2005, J.Ahrensfeld +// ------------------------------------------------------------ + +// ------------------------------------------------------------ +// LAPACK functions +// ------------------------------------------------------------ +// SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) +extern +void dgesv_(int *N, int *NHRS, double *A, int *LDA, + int *IPIV, double *B, int *LDB, int *INFO); + +int dgesv(int N, int NHRS, double *A, int LDA, int *IPIV, + double *B, int LDB) +{ + int info; + + dgesv_(&N, &NHRS, A, &LDA, IPIV, B, &LDB, &info); + + return info; +} + +// ------------------------------------------------------------ +// SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) +extern +void dsyev_(char *JOBZ, char *UPLO, int *N, double *A, int *LDA, + double *W, double *WORK, int *LWORK, int *INFO); + +int dsyev(char JOBZ, char UPLO, int N, double *A, int LDA, + double *W, double *WORK, int LWORK) +{ + int info; + + dsyev_(&JOBZ, &UPLO, &N, A, &LDA, W, WORK, &LWORK, &info); + + return info; +} + +// ------------------------------------------------------------ +// Extern declaration of FORTRAN routine +extern +void dgtsv_(int *Np, int *NRHSp, double *DL, double *D, double *DU, + double *B, int *LDBp, int *INFOp); + +// C-Wrapper function +int dgtsv(int N, int NRHS, double *DL, double *D, double *DU, + double *B, int ldb) +{ + int info; + dgtsv_ (&N, &NRHS, DL, D, DU, B, &ldb, &info); + + return info; +} + + +// ------------------------------------------------------------ +// SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO ) +extern +void dgebrd_(int *M, int *N, double *A, int *LDA, double *D, double *E, + double *TAUQ, double *TAUP, double *WORK, int *LWORK, int *INFO); + +int dgebrd(int M, int N, double *A, int LDA, double *D, double *E, + double *TAUQ, double *TAUP, double *WORK, int LWORK) + +{ + int info; + dgebrd_(&M, &N, A, &LDA, D, E, TAUQ, TAUP, WORK, &LWORK, &info); + + return info; +} + +// ------------------------------------------------------------ +// SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, +// LDU, C, LDC, WORK, INFO ) +extern +void dbdsqr_(char *UPLO, int *N, int *NCVT, int *NRU, int *NCC, double *D, + double *E, double *VT, int *LDVT, double *U, int *LDU, + double *C, int *LDC, double *WORK, int *INFO); + +int dbdsqr(char UPLO, int N, int NCVT, int NRU, int NCC, double *D, + double *E, double *VT, int LDVT, double *U, int LDU, + double *C, int LDC, double *WORK) +{ + int info; + dbdsqr_(&UPLO, &N, &NCVT, &NRU, &NCC, D, E, VT, &LDVT, U, &LDU, + C, &LDC, WORK, &info); + + return info; +} + +int dgesvd(char JOBU, char JOBVT, int M, int N, double *A, int LDA, double *S, + double *U, int LDU, double *VT, int LDVT, double *WORK, int LWORK); + +// ------------------------------------------------------------ +// SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, +// WORK, LWORK, INFO ) +extern +void dgesvd_(char *JOBU, char *JOBVT, int *M, int *N, double *A, int *LDA, + double *S, double *U, int *LDU, double *VT, int *LDVT, + double *WORK, int *LWORK, int *INFO); + +int dgesvd(char JOBU, char JOBVT, int M, int N, double *A, int LDA, double *S, + double *U, int LDU, double *VT, int LDVT, double *WORK, int LWORK) +{ + int info; + dgesvd_(&JOBU, &JOBVT, &M, &N, A, &LDA, S, U, &LDU, VT, &LDVT, + WORK, &LWORK, &info); + + + return info; +} + +// ------------------------------------------------------------ +// SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) +extern +void dgetri_(int *N, double *A, int *LDA, int *IPIV, double *WORK, int *LWORK, int *INFO); + +int dgetri(int N, double *A, int LDA, int *IPIV, double *WORK, int LWORK) +{ + + int info; + dgetri_(&N, A, &LDA, IPIV, WORK, &LWORK, &info); + + return info; +} + +// ------------------------------------------------------------ +// SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO ) +extern +void dpotri_(char *UPLO, int *N, double *A, int *LDA, int *INFO); + +int dpotri(char UPLO, int N, double *A, int LDA) +{ + + int info; + dpotri_(&UPLO, &N, A, &LDA, &info); + + return info; +} + +// ------------------------------------------------------------ +// SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, INFO ) +extern +void dgecon_(char *NORM, int *N, double *A, int *LDA, double *ANORM, double *RCOND, + double *DWORK, int *IWORK, int *INFO); + +int dgecon(char NORM, int N, double *A, int LDA, double ANORM, double *RCOND, + double *DWORK, int *IWORK) +{ + int info; + dgecon_(&NORM, &N, A, &LDA, &ANORM, RCOND, DWORK, IWORK, &info); + + return info; +} + +// ------------------------------------------------------------ +// BLAS functions LEVEL 1 +// ------------------------------------------------------------ +// subroutine dswap(n,dx,incx,dy,incy) +// interchanges two vectors. +extern +void dswap_ (int *n, double *dx, int *incx, double *dy, int *incy); + +void dswap (int n, double *dx, int incx, double *dy, int incy) +{ + dswap_ (&n, dx, &incx, dy, &incy); +} + +// ------------------------------------------------------------ +// subroutine dcopy(n,dx,incx,dy,incy) +// copies a vector, x, to a vector, y. +extern +void dcopy_ (int *n, double *dx, int *incx, double *dy, int *incy); + +void dcopy (int n, double *dx, int incx, double *dy, int incy) +{ + dcopy_ (&n, dx, &incx, dy, &incy); +} + +// ------------------------------------------------------------ +// subroutine dscal(n,da,dx,incx) +// scales a vector by a constant. +extern +void dscal_ (int *n, double *da, double *dx, int *incx); + +void dscal (int n, double da, double *dx, int incx) +{ + dscal_ (&n, &da, dx, &incx); +} + +// ------------------------------------------------------------ +// subroutine daxpy(n,da,dx,incx,dy,incy) +// constant times a vector plus a vector. +extern +void daxpy_ (int *n, double *da, double *dx, int *incx, double *dy, int *incy); + +void daxpy (int n, double da, double *dx, int incx, double *dy, int incy) +{ + daxpy_ (&n, &da, dx, &incx, dy, &incy); +} + +// ------------------------------------------------------------ +// function ddot(n,dx,incx,dy,incy) +// forms the dot product of two vectors. +extern +double ddot_ (int *n, double *dx, int *incx, double *dy, int *incy); + +double ddot (int n, double *dx, int incx, double *dy, int incy) +{ + return ddot_ (&n, dx, &incx, dy, &incy); +} + +// ------------------------------------------------------------ +// function DNRM2 ( N, X, INCX ) +// DNRM2 returns the euclidean norm of a vector via the function +// DNRM2 := sqrt( x'*x ) +extern +double dnrm2_ (int *n, double *dx, int *incx); + +double dnrm2 (int n, double *dx, int incx) +{ + return dnrm2_ (&n, dx, &incx); +} + +// ------------------------------------------------------------ +// function dasum(n,dx,incx) +// takes the sum of the absolute values. +extern +double dasum_ (int *n, double *dx, int *incx); + +double dasum (int n, double *dx, int incx) +{ + return dasum_ (&n, dx, &incx); +} + +// ------------------------------------------------------------ +// function idamax(n,dx,incx) +// finds the index of element having max. absolute value. +extern +int idamax_ (int *n, double *dx, int *incx); + +int idamax (int n, double *dx, int incx) +{ + int i; + + i = idamax_ (&n, dx, &incx); + + if (i > 0) + return i - 1; + else + return i; + +} + +// ------------------------------------------------------------ +// BLAS functions LEVEL 2 +// ------------------------------------------------------------ + +// SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) +extern +void dgemv_(char *TRANS, int *M, int *N, double *ALPHA, double *A, int *LDA, + double *X, int *INCX, double *BETA, double *Y, int *INCY); + +void dgemv(char TRANS, int M, int N, double ALPHA, double *A, int LDA, + double *X, int INCX, double BETA, double *Y, int INCY) +{ + + dgemv_(&TRANS, &M, &N, &ALPHA, A, &LDA, X, &INCX, &BETA, Y, &INCY); +} + +// ------------------------------------------------------------ +// SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) +extern void dger_(int *m,int *n, double *alpha, double *x, int *incx, + double *y, int *incy, double *A, int *lda); + +void dger(int m,int n, double alpha, double *x, int incx, + double *y, int incy, double *A, int lda) +{ + dger_(&m, &n, &alpha, x, &incx, y, &incy, A, &lda); +} + +// ------------------------------------------------------------ +// SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) +extern void dsyr_(char *UPLO, int *n, double *alpha, double *x, int *incx, + double *A, int *lda); + +void dsyr(char UPLO, int n, double alpha, double *x, int incx, double *A, int lda) +{ + dsyr_(&UPLO, &n, &alpha, x, &incx, A, &lda); +} + +// ------------------------------------------------------------ +// BLAS functions LEVEL 3 +// ------------------------------------------------------------ + +// SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, +extern +void dsyrk_(char *UPLO, char *TRANS, int *N, int *K, double *ALPHA, + double *A, int *LDA, double *BETA, double *C, int *LDC); + +void dsyrk(char UPLO, char TRANS, int N, int K, double ALPHA, + double *A, int LDA, double BETA, double *C, int LDC) +{ + dsyrk_(&UPLO, &TRANS, &N, &K, &ALPHA, A, &LDA, &BETA, C, &LDC); +} + +// SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, +extern +void dgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, double *ALPHA, + double *A, int *LDA, double *B, int *LDB, double *BETA, double *C, int *LDC); + +void dgemm(char TRANSA, char TRANSB, int M, int N, int K, double ALPHA, + double *A, int LDA, double *B, int LDB, double BETA, double *C, int LDC) +{ + + dgemm_(&TRANSA, &TRANSB, &M, &N, &K, &ALPHA, A, &LDA, B, &LDB, &BETA, C, &LDC); + +} + +// SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, +extern +void dtrmm_(char *SIDE, char *UPLO, char *TRANSA, char *DIAG, int *M, int *N, + double *ALPHA, double *A, int *LDA, double *B, int *LDB); + +void dtrmm(char SIDE, char UPLO, char TRANSA, char DIAG, int M, int N, + double ALPHA, double *A, int LDA, double *B, int LDB) +{ + + dtrmm_(&SIDE, &UPLO, &TRANSA, &DIAG, &M, &N, &ALPHA, A, &LDA, B, &LDB); + +} + +// SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, +extern +void dtrsm_(char *SIDE, char *UPLO, char *TRANSA, char *DIAG, int *M, int *N, + double *ALPHA, double *A, int *LDA, double *B, int *LDB); + +void dtrsm(char SIDE, char UPLO, char TRANSA, char DIAG, int M, int N, + double ALPHA, double *A, int LDA, double *B, int LDB) +{ + + dtrsm_(&SIDE, &UPLO, &TRANSA, &DIAG, &M, &N, &ALPHA, A, &LDA, B, &LDB); + +} + diff --git a/lapack_cwrap.h b/lapack_cwrap.h new file mode 100755 index 0000000..1ead4a9 --- /dev/null +++ b/lapack_cwrap.h @@ -0,0 +1,1389 @@ +// ------------------------------------------------------------ +// lapack_cwrap.h +// Wrapping LAPACK FORTRAN-function to C-Functions +// +// 26.02.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef LAPACK_CWRAP_H +#define LAPACK_CWRAP_H + +// ------------------------------------------------------------ +// LAPACK functions +// ------------------------------------------------------------ +int dgtsv(int N, int NRHS, double *DL, double *D, double *DU, + double *B, int ldb); + +// ------------------------------------------------------------ +int dgesv(int N, int NHRS, double *A, int LDA, int *IPIV, + double *B, int LDB); + +// ------------------------------------------------------------ +// SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) +// +// -- LAPACK driver routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// March 31, 1993 +// +// .. Scalar Arguments .. +// INTEGER INFO, LDA, LDB, N, NRHS +// .. +// .. Array Arguments .. +// INTEGER IPIV( * ) +// DOUBLE PRECISION A( LDA, * ), B( LDB, * ) +// .. +// +// Purpose +// ======= +// +// DGESV computes the solution to a real system of linear equations +// A * X = B, +// where A is an N-by-N matrix and X and B are N-by-NRHS matrices. +// +// The LU decomposition with partial pivoting and row interchanges is +// used to factor A as +// A = P * L * U, +// where P is a permutation matrix, L is unit lower triangular, and U is +// upper triangular. The factored form of A is then used to solve the +// system of equations A * X = B. +// +// Arguments +// ========= +// +// N (input) INTEGER +// The number of linear equations, i.e., the order of the +// matrix A. N >= 0. +// +// NRHS (input) INTEGER +// The number of right hand sides, i.e., the number of columns +// of the matrix B. NRHS >= 0. +// +// A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +// On entry, the N-by-N coefficient matrix A. +// On exit, the factors L and U from the factorization +// A = P*L*U; the unit diagonal elements of L are not stored. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,N). +// +// IPIV (output) INTEGER array, dimension (N) +// The pivot indices that define the permutation matrix P; +// row i of the matrix was interchanged with row IPIV(i). +// +// B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) +// On entry, the N-by-NRHS matrix of right hand side matrix B. +// On exit, if INFO = 0, the N-by-NRHS solution matrix X. +// +// LDB (input) INTEGER +// The leading dimension of the array B. LDB >= max(1,N). +// +// INFO (output) INTEGER +// = 0: successful exit +// < 0: if INFO = -i, the i-th argument had an illegal value +// > 0: if INFO = i, U(i,i) is exactly zero. The factorization +// has been completed, but the factor U is exactly +// singular, so the solution could not be computed. + +// ------------------------------------------------------------ +int dsyev(char JOBZ, char UPLO, int N, double *A, int LDA, + double *W, double *WORK, int LWORK); +// ------------------------------------------------------------ +// SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) +// +// DSYEV computes all eigenvalues and, optionally, eigenvectors of a +// real symmetric matrix A. +// +// Arguments +// ========= +// +// JOBZ (input) CHARACTER*1 +// = 'N': Compute eigenvalues only; +// = 'V': Compute eigenvalues and eigenvectors. +// +// UPLO (input) CHARACTER*1 +// = 'U': Upper triangle of A is stored; +// = 'L': Lower triangle of A is stored. +// +// N (input) INTEGER +// The order of the matrix A. N >= 0. +// +// A (input/output) DOUBLE PRECISION array, dimension (LDA, N) +// On entry, the symmetric matrix A. If UPLO = 'U', the +// leading N-by-N upper triangular part of A contains the +// upper triangular part of the matrix A. If UPLO = 'L', +// the leading N-by-N lower triangular part of A contains +// the lower triangular part of the matrix A. +// On exit, if JOBZ = 'V', then if INFO = 0, A contains the +// orthonormal eigenvectors of the matrix A. +// If JOBZ = 'N', then on exit the lower triangle (if UPLO='L') +// or the upper triangle (if UPLO='U') of A, including the +// diagonal, is destroyed. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,N). +// +// W (output) DOUBLE PRECISION array, dimension (N) +// If INFO = 0, the eigenvalues in ascending order. +// +// WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) +// On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +// +// LWORK (input) INTEGER +// The length of the array WORK. LWORK >= max(1,3*N-1). +// For optimal efficiency, LWORK >= (NB+2)*N, +// where NB is the blocksize for DSYTRD returned by ILAENV. +// +// If LWORK = -1, then a workspace query is assumed; the routine +// only calculates the optimal size of the WORK array, returns +// this value as the first entry of the WORK array, and no error +// message related to LWORK is issued by XERBLA. +// +// INFO (output) INTEGER +// = 0: successful exit +// < 0: if INFO = -i, the i-th argument had an illegal value +// > 0: if INFO = i, the algorithm failed to converge; i +// off-diagonal elements of an intermediate tridiagonal +// form did not converge to zero. + +// ------------------------------------------------------------ +int dgebrd(int M, int N, double *A, int LDA, double *D, double *E, + double *TAUQ, double *TAUP, double *WORK, int LWORK); +// ------------------------------------------------------------ +// SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, +// $ INFO ) +// +// -- LAPACK routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// June 30, 1999 +// +// .. Scalar Arguments .. +// INTEGER INFO, LDA, LWORK, M, N +// .. +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ), +// $ TAUQ( * ), WORK( * ) +// .. +// +// Purpose +// ======= +// +// DGEBRD reduces a general real M-by-N matrix A to upper or lower +// bidiagonal form B by an orthogonal transformation: Q**T * A * P = B. +// +// If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal. +// +// Arguments +// ========= +// +// M (input) INTEGER +// The number of rows in the matrix A. M >= 0. +// +// N (input) INTEGER +// The number of columns in the matrix A. N >= 0. +// +// A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +// On entry, the M-by-N general matrix to be reduced. +// On exit, +// if m >= n, the diagonal and the first superdiagonal are +// overwritten with the upper bidiagonal matrix B; the +// elements below the diagonal, with the array TAUQ, represent +// the orthogonal matrix Q as a product of elementary +// reflectors, and the elements above the first superdiagonal, +// with the array TAUP, represent the orthogonal matrix P as +// a product of elementary reflectors; +// if m < n, the diagonal and the first subdiagonal are +// overwritten with the lower bidiagonal matrix B; the +// elements below the first subdiagonal, with the array TAUQ, +// represent the orthogonal matrix Q as a product of +// elementary reflectors, and the elements above the diagonal, +// with the array TAUP, represent the orthogonal matrix P as +// a product of elementary reflectors. +// See Further Details. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,M). +// +// D (output) DOUBLE PRECISION array, dimension (min(M,N)) +// The diagonal elements of the bidiagonal matrix B: +// D(i) = A(i,i). +// +// E (output) DOUBLE PRECISION array, dimension (min(M,N)-1) +// The off-diagonal elements of the bidiagonal matrix B: +// if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1; +// if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1. +// +// TAUQ (output) DOUBLE PRECISION array dimension (min(M,N)) +// The scalar factors of the elementary reflectors which +// represent the orthogonal matrix Q. See Further Details. +// +// TAUP (output) DOUBLE PRECISION array, dimension (min(M,N)) +// The scalar factors of the elementary reflectors which +// represent the orthogonal matrix P. See Further Details. +// +// WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) +// On exit, if INFO = 0, WORK(1) returns the optimal LWORK. +// +// LWORK (input) INTEGER +// The length of the array WORK. LWORK >= max(1,M,N). +// For optimum performance LWORK >= (M+N)*NB, where NB +// is the optimal blocksize. +// +// If LWORK = -1, then a workspace query is assumed; the routine +// only calculates the optimal size of the WORK array, returns +// this value as the first entry of the WORK array, and no error +// message related to LWORK is issued by XERBLA. +// +// INFO (output) INTEGER +// = 0: successful exit +// < 0: if INFO = -i, the i-th argument had an illegal value. +// + +// ------------------------------------------------------------ +int dbdsqr(char UPLO, int N, int NCVT, int NRU, int NCC, double *D, + double *E, double *VT, int LDVT, double *U, int LDU, + double *C, int LDC, double *WORK); +// ------------------------------------------------------------ +// SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, +// $ LDU, C, LDC, WORK, INFO ) +// +// -- LAPACK routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// October 31, 1999 +// +// .. Scalar Arguments .. +// CHARACTER UPLO +// INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU +// .. +// .. Array Arguments .. +// DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), +// $ VT( LDVT, * ), WORK( * ) +// .. +// +// Purpose +// ======= +// +// DBDSQR computes the singular value decomposition (SVD) of a real +// N-by-N (upper or lower) bidiagonal matrix B: B = Q * S * P' (P' +// denotes the transpose of P), where S is a diagonal matrix with +// non-negative diagonal elements (the singular values of B), and Q +// and P are orthogonal matrices. +// +// The routine computes S, and optionally computes U * Q, P' * VT, +// or Q' * C, for given real input matrices U, VT, and C. +// +// See "Computing Small Singular Values of Bidiagonal Matrices With +// Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, +// LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, +// no. 5, pp. 873-912, Sept 1990) and +// "Accurate singular values and differential qd algorithms," by +// B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics +// Department, University of California at Berkeley, July 1992 +// for a detailed description of the algorithm. +// +// Arguments +// ========= +// +// UPLO (input) CHARACTER*1 +// = 'U': B is upper bidiagonal; +// = 'L': B is lower bidiagonal. +// +// N (input) INTEGER +// The order of the matrix B. N >= 0. +// +// NCVT (input) INTEGER +// The number of columns of the matrix VT. NCVT >= 0. +// +// NRU (input) INTEGER +// The number of rows of the matrix U. NRU >= 0. +// +// NCC (input) INTEGER +// The number of columns of the matrix C. NCC >= 0. +// +// D (input/output) DOUBLE PRECISION array, dimension (N) +// On entry, the n diagonal elements of the bidiagonal matrix B. +// On exit, if INFO=0, the singular values of B in decreasing +// order. +// +// E (input/output) DOUBLE PRECISION array, dimension (N) +// On entry, the elements of E contain the +// offdiagonal elements of the bidiagonal matrix whose SVD +// is desired. On normal exit (INFO = 0), E is destroyed. +// If the algorithm does not converge (INFO > 0), D and E +// will contain the diagonal and superdiagonal elements of a +// bidiagonal matrix orthogonally equivalent to the one given +// as input. E(N) is used for workspace. +// +// VT (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT) +// On entry, an N-by-NCVT matrix VT. +// On exit, VT is overwritten by P' * VT. +// VT is not referenced if NCVT = 0. +// +// LDVT (input) INTEGER +// The leading dimension of the array VT. +// LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0. +// +// U (input/output) DOUBLE PRECISION array, dimension (LDU, N) +// On entry, an NRU-by-N matrix U. +// On exit, U is overwritten by U * Q. +// U is not referenced if NRU = 0. +// +// LDU (input) INTEGER +// The leading dimension of the array U. LDU >= max(1,NRU). +// +// C (input/output) DOUBLE PRECISION array, dimension (LDC, NCC) +// On entry, an N-by-NCC matrix C. +// On exit, C is overwritten by Q' * C. +// C is not referenced if NCC = 0. +// +// LDC (input) INTEGER +// The leading dimension of the array C. +// LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. +// +// WORK (workspace) DOUBLE PRECISION array, dimension (4*N) +// +// INFO (output) INTEGER +// = 0: successful exit +// < 0: If INFO = -i, the i-th argument had an illegal value +// > 0: the algorithm did not converge; D and E contain the +// elements of a bidiagonal matrix which is orthogonally +// similar to the input matrix B; if INFO = i, i +// elements of E have not converged to zero. +// + +// ------------------------------------------------------------ +int dgesvd(char JOBU, char JOBVT, int M, int N, double *A, int LDA, double *S, + double *U, int LDU, double *VT, int LDVT, double *WORK, int LWORK); + +// ------------------------------------------------------------ +// SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, +// $ WORK, LWORK, INFO ) +// +// -- LAPACK driver routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// October 31, 1999 +// +// .. Scalar Arguments .. +// CHARACTER JOBU, JOBVT +// INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N +// .. +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ), +// $ VT( LDVT, * ), WORK( * ) +// .. +// +// Purpose +// ======= +// +// DGESVD computes the singular value decomposition (SVD) of a real +// M-by-N matrix A, optionally computing the left and/or right singular +// vectors. The SVD is written +// +// A = U * SIGMA * transpose(V) +// +// where SIGMA is an M-by-N matrix which is zero except for its +// min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and +// V is an N-by-N orthogonal matrix. The diagonal elements of SIGMA +// are the singular values of A; they are real and non-negative, and +// are returned in descending order. The first min(m,n) columns of +// U and V are the left and right singular vectors of A. +// +// Note that the routine returns V**T, not V. +// +// Arguments +// ========= +// +// JOBU (input) CHARACTER*1 +// Specifies options for computing all or part of the matrix U: +// = 'A': all M columns of U are returned in array U: +// = 'S': the first min(m,n) columns of U (the left singular +// vectors) are returned in the array U; +// = 'O': the first min(m,n) columns of U (the left singular +// vectors) are overwritten on the array A; +// = 'N': no columns of U (no left singular vectors) are +// computed. +// +// JOBVT (input) CHARACTER*1 +// Specifies options for computing all or part of the matrix +// V**T: +// = 'A': all N rows of V**T are returned in the array VT; +// = 'S': the first min(m,n) rows of V**T (the right singular +// vectors) are returned in the array VT; +// = 'O': the first min(m,n) rows of V**T (the right singular +// vectors) are overwritten on the array A; +// = 'N': no rows of V**T (no right singular vectors) are +// computed. +// +// JOBVT and JOBU cannot both be 'O'. +// +// M (input) INTEGER +// The number of rows of the input matrix A. M >= 0. +// +// N (input) INTEGER +// The number of columns of the input matrix A. N >= 0. +// +// A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +// On entry, the M-by-N matrix A. +// On exit, +// if JOBU = 'O', A is overwritten with the first min(m,n) +// columns of U (the left singular vectors, +// stored columnwise); +// if JOBVT = 'O', A is overwritten with the first min(m,n) +// rows of V**T (the right singular vectors, +// stored rowwise); +// if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A +// are destroyed. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,M). +// +// S (output) DOUBLE PRECISION array, dimension (min(M,N)) +// The singular values of A, sorted so that S(i) >= S(i+1). +// +// U (output) DOUBLE PRECISION array, dimension (LDU,UCOL) +// (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'. +// If JOBU = 'A', U contains the M-by-M orthogonal matrix U; +// if JOBU = 'S', U contains the first min(m,n) columns of U +// (the left singular vectors, stored columnwise); +// if JOBU = 'N' or 'O', U is not referenced. +// +// LDU (input) INTEGER +// The leading dimension of the array U. LDU >= 1; if +// JOBU = 'S' or 'A', LDU >= M. +// +// VT (output) DOUBLE PRECISION array, dimension (LDVT,N) +// If JOBVT = 'A', VT contains the N-by-N orthogonal matrix +// V**T; +// if JOBVT = 'S', VT contains the first min(m,n) rows of +// V**T (the right singular vectors, stored rowwise); +// if JOBVT = 'N' or 'O', VT is not referenced. +// +// LDVT (input) INTEGER +// The leading dimension of the array VT. LDVT >= 1; if +// JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). +// +// WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) +// On exit, if INFO = 0, WORK(1) returns the optimal LWORK; +// if INFO > 0, WORK(2:MIN(M,N)) contains the unconverged +// superdiagonal elements of an upper bidiagonal matrix B +// whose diagonal is in S (not necessarily sorted). B +// satisfies A = U * B * VT, so it has the same singular values +// as A, and singular vectors related by U and VT. +// +// LWORK (input) INTEGER +// The dimension of the array WORK. LWORK >= 1. +// LWORK >= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). +// For good performance, LWORK should generally be larger. +// +// If LWORK = -1, then a workspace query is assumed; the routine +// only calculates the optimal size of the WORK array, returns +// this value as the first entry of the WORK array, and no error +// message related to LWORK is issued by XERBLA. +// +// INFO (output) INTEGER +// = 0: successful exit. +// < 0: if INFO = -i, the i-th argument had an illegal value. +// > 0: if DBDSQR did not converge, INFO specifies how many +// superdiagonals of an intermediate bidiagonal form B +// did not converge to zero. See the description of WORK +// above for details. +// + +int dgecon(char NORM, int N, double *A, int LDA, double ANORM, double *RCOND, + double *DWORK, int *IWORK); + +// SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK, +// $ INFO ) +// +// -- LAPACK routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// February 29, 1992 +// +// .. Scalar Arguments .. +// CHARACTER NORM +// INTEGER INFO, LDA, N +// DOUBLE PRECISION ANORM, RCOND +// .. +// .. Array Arguments .. +// INTEGER IWORK( * ) +// DOUBLE PRECISION A( LDA, * ), WORK( * ) +// .. +// +// Purpose +// ======= +// +// DGECON estimates the reciprocal of the condition number of a general +// real matrix A, in either the 1-norm or the infinity-norm, using +// the LU factorization computed by DGETRF. +// +// An estimate is obtained for norm(inv(A)), and the reciprocal of the +// condition number is computed as +// RCOND = 1 / ( norm(A) * norm(inv(A)) ). +// +// Arguments +// ========= +// +// NORM (input) CHARACTER*1 +// Specifies whether the 1-norm condition number or the +// infinity-norm condition number is required: +// = '1' or 'O': 1-norm; +// = 'I': Infinity-norm. +// +// N (input) INTEGER +// The order of the matrix A. N >= 0. +// +// A (input) DOUBLE PRECISION array, dimension (LDA,N) +// The factors L and U from the factorization A = P*L*U +// as computed by DGETRF. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,N). +// +// ANORM (input) DOUBLE PRECISION +// If NORM = '1' or 'O', the 1-norm of the original matrix A. +// If NORM = 'I', the infinity-norm of the original matrix A. +// +// RCOND (output) DOUBLE PRECISION +// The reciprocal of the condition number of the matrix A, +// computed as RCOND = 1/(norm(A) * norm(inv(A))). +// +// WORK (workspace) DOUBLE PRECISION array, dimension (4*N) +// +// IWORK (workspace) INTEGER array, dimension (N) +// +// INFO (output) INTEGER +// = 0: successful exit + +// ------------------------------------------------------------ +// BLAS functions LEVEL 1 +// ------------------------------------------------------------ +void dswap (int n, double *dx, int incx, double *dy, int incy); +void dcopy (int n, double *dx, int incx, double *dy, int incy); +void dscal (int N, double da, double *dx, int incx); +void daxpy (int n, double da, double *dx, int incx, double *dy, int incy); +double ddot (int n, double *dx, int incx, double *dy, int incy); +double dnrm2 (int n, double *dx, int incx); +double dasum (int n, double *dx, int incx); +int idamax (int n, double *dx, int incx); + +// ------------------------------------------------------------ +// BLAS functions LEVEL 2 +// ------------------------------------------------------------ + + +// ------------------------------------------------------------ +void dgemv(char TRANS, int M, int N, double ALPHA, double *A, int LDA, + double *X, int INCX, double BETA, double *Y, int INCY); +// ------------------------------------------------------------ +// SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY ) +// +// DGEMV performs one of the matrix-vector operations +// +// y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, +// +// where alpha and beta are scalars, x and y are vectors and A is an +// m by n matrix. +// +// Parameters +// ========== +// +// TRANS - CHARACTER*1. +// On entry, TRANS specifies the operation to be performed as +// follows: +// +// TRANS = 'N' or 'n' y := alpha*A*x + beta*y. +// +// TRANS = 'T' or 't' y := alpha*A'*x + beta*y. +// +// TRANS = 'C' or 'c' y := alpha*A'*x + beta*y. +// +// Unchanged on exit. +// +// M - INTEGER. +// On entry, M specifies the number of rows of the matrix A. +// M must be at least zero. +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the number of columns of the matrix A. +// N must be at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). +// Before entry, the leading m by n part of the array A must +// contain the matrix of coefficients. +// Unchanged on exit. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. LDA must be at least +// max( 1, m ). +// Unchanged on exit. +// +// X - DOUBLE PRECISION array of DIMENSION at least +// ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n' +// and at least +// ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. +// Before entry, the incremented array X must contain the +// vector x. +// Unchanged on exit. +// +// INCX - INTEGER. +// On entry, INCX specifies the increment for the elements of +// X. INCX must not be zero. +// Unchanged on exit. +// +// BETA - DOUBLE PRECISION. +// On entry, BETA specifies the scalar beta. When BETA is +// supplied as zero then Y need not be set on input. +// Unchanged on exit. +// +// Y - DOUBLE PRECISION array of DIMENSION at least +// ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n' +// and at least +// ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. +// Before entry with BETA non-zero, the incremented array Y +// must contain the vector y. On exit, Y is overwritten by the +// updated vector y. +// +// INCY - INTEGER. +// On entry, INCY specifies the increment for the elements of +// Y. INCY must not be zero. +// Unchanged on exit. + + +// ------------------------------------------------------------ +// BLAS functions LEVEL 2 +// ------------------------------------------------------------ + +void dger(int m,int n, double alpha, double *x, int incx, double *y, int incy, double *A, int lda); +// ------------------------------------------------------------ +// SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) +// DGER performs the rank 1 operation +// A := alpha*x*y' + A, +// where alpha is a scalar, x is an m element vector, y is an n element +// vector and A is an m by n matrix. +// +// Parameters +// ========== +// +// M - INTEGER. +// On entry, M specifies the number of rows of the matrix A. +// M must be at least zero. +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the number of columns of the matrix A. +// N must be at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. +// Unchanged on exit. +// +// X - DOUBLE PRECISION array of dimension at least +// ( 1 + ( m - 1 )//abs( INCX ) ). +// Before entry, the incremented array X must contain the m +// element vector x. +// Unchanged on exit. +// +// INCX - INTEGER. +// On entry, INCX specifies the increment for the elements of +// X. INCX must not be zero. +// Unchanged on exit. +// +// Y - DOUBLE PRECISION array of dimension at least +// ( 1 + ( n - 1 )//abs( INCY ) ). +// Before entry, the incremented array Y must contain the n +// element vector y. +// Unchanged on exit. +// +// INCY - INTEGER. +// On entry, INCY specifies the increment for the elements of +// Y. INCY must not be zero. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). +// Before entry, the leading m by n part of the array A must +// contain the matrix of coefficients. On exit, A is +// overwritten by the updated matrix. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. LDA must be at least +// max( 1, m ). +// Unchanged on exit. + +// ------------------------------------------------------------ +void dsyr(char UPLO, int n, double alpha, double *x, int incx, double *A, int lda); +// ------------------------------------------------------------ +// SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) +// .. Scalar Arguments .. +// DOUBLE PRECISION ALPHA +// INTEGER INCX, LDA, N +// CHARACTER*1 UPLO +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ), X( * ) +// .. +// +// Purpose +// ======= +// +// DSYR performs the symmetric rank 1 operation +// +// A := alpha*x*x' + A, +// +// where alpha is a real scalar, x is an n element vector and A is an +// n by n symmetric matrix. +// +// Parameters +// ========== +// +// UPLO - CHARACTER*1. +// On entry, UPLO specifies whether the upper or lower +// triangular part of the array A is to be referenced as +// follows: +// +// UPLO = 'U' or 'u' Only the upper triangular part of A +// is to be referenced. +// +// UPLO = 'L' or 'l' Only the lower triangular part of A +// is to be referenced. +// +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the order of the matrix A. +// N must be at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. +// Unchanged on exit. +// +// X - DOUBLE PRECISION array of dimension at least +// ( 1 + ( n - 1 )*abs( INCX ) ). +// Before entry, the incremented array X must contain the n +// element vector x. +// Unchanged on exit. +// +// INCX - INTEGER. +// On entry, INCX specifies the increment for the elements of +// X. INCX must not be zero. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, n ). +// Before entry with UPLO = 'U' or 'u', the leading n by n +// upper triangular part of the array A must contain the upper +// triangular part of the symmetric matrix and the strictly +// lower triangular part of A is not referenced. On exit, the +// upper triangular part of the array A is overwritten by the +// upper triangular part of the updated matrix. +// Before entry with UPLO = 'L' or 'l', the leading n by n +// lower triangular part of the array A must contain the lower +// triangular part of the symmetric matrix and the strictly +// upper triangular part of A is not referenced. On exit, the +// lower triangular part of the array A is overwritten by the +// lower triangular part of the updated matrix. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. LDA must be at least +// max( 1, n ). +// Unchanged on exit. + +// ------------------------------------------------------------ +// BLAS functions LEVEL 3 +// ------------------------------------------------------------ + +void dsyrk(char UPLO, char TRANS, int N, int K, double ALPHA, + double *A, int LDA, double BETA, double *C, int LDC); +// ------------------------------------------------------------ +// SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, +// $ BETA, C, LDC ) +// .. Scalar Arguments .. +// CHARACTER*1 UPLO, TRANS +// INTEGER N, K, LDA, LDC +// DOUBLE PRECISION ALPHA, BETA +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ), C( LDC, * ) +// .. +// +// Purpose +// ======= +// +// DSYRK performs one of the symmetric rank k operations +// +// C := alpha*A*A' + beta*C, +// +// or +// +// C := alpha*A'*A + beta*C, +// +// where alpha and beta are scalars, C is an n by n symmetric matrix +// and A is an n by k matrix in the first case and a k by n matrix +// in the second case. +// +// Parameters +// ========== +// +// UPLO - CHARACTER*1. +// On entry, UPLO specifies whether the upper or lower +// triangular part of the array C is to be referenced as +// follows: +// +// UPLO = 'U' or 'u' Only the upper triangular part of C +// is to be referenced. +// +// UPLO = 'L' or 'l' Only the lower triangular part of C +// is to be referenced. +// +// Unchanged on exit. +// +// TRANS - CHARACTER*1. +// On entry, TRANS specifies the operation to be performed as +// follows: +// +// TRANS = 'N' or 'n' C := alpha*A*A' + beta*C. +// +// TRANS = 'T' or 't' C := alpha*A'*A + beta*C. +// +// TRANS = 'C' or 'c' C := alpha*A'*A + beta*C. +// +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the order of the matrix C. N must be +// at least zero. +// Unchanged on exit. +// +// K - INTEGER. +// On entry with TRANS = 'N' or 'n', K specifies the number +// of columns of the matrix A, and on entry with +// TRANS = 'T' or 't' or 'C' or 'c', K specifies the number +// of rows of the matrix A. K must be at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is +// k when TRANS = 'N' or 'n', and is n otherwise. +// Before entry with TRANS = 'N' or 'n', the leading n by k +// part of the array A must contain the matrix A, otherwise +// the leading k by n part of the array A must contain the +// matrix A. +// Unchanged on exit. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. When TRANS = 'N' or 'n' +// then LDA must be at least max( 1, n ), otherwise LDA must +// be at least max( 1, k ). +// Unchanged on exit. +// +// BETA - DOUBLE PRECISION. +// On entry, BETA specifies the scalar beta. +// Unchanged on exit. +// +// C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). +// Before entry with UPLO = 'U' or 'u', the leading n by n +// upper triangular part of the array C must contain the upper +// triangular part of the symmetric matrix and the strictly +// lower triangular part of C is not referenced. On exit, the +// upper triangular part of the array C is overwritten by the +// upper triangular part of the updated matrix. +// Before entry with UPLO = 'L' or 'l', the leading n by n +// lower triangular part of the array C must contain the lower +// triangular part of the symmetric matrix and the strictly +// upper triangular part of C is not referenced. On exit, the +// lower triangular part of the array C is overwritten by the +// lower triangular part of the updated matrix. +// +// LDC - INTEGER. +// On entry, LDC specifies the first dimension of C as declared +// in the calling (sub) program. LDC must be at least +// max( 1, n ). +// Unchanged on exit. +// ------------------------------------------------------------ + +void dgemm(char TRANSA, char TRANSB, int M, int N, int K, double ALPHA, + double *A, int LDA, double *B, int LDB, double BETA, double *C, int LDC); +// ------------------------------------------------------------ +// SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC ) +// +// DGEMM performs one of the matrix-matrix operations +// +// C := alpha*op( A )*op( B ) + beta*C, +// +// where op( X ) is one of +// +// op( X ) = X or op( X ) = X', +// +// alpha and beta are scalars, and A, B and C are matrices, with op( A ) +// an m by k matrix, op( B ) a k by n matrix and C an m by n matrix. +// +// Parameters +// ========== +// +// TRANSA - CHARACTER*1. +// On entry, TRANSA specifies the form of op( A ) to be used in +// the matrix multiplication as follows: +// +// TRANSA = 'N' or 'n', op( A ) = A. +// +// TRANSA = 'T' or 't', op( A ) = A'. +// +// TRANSA = 'C' or 'c', op( A ) = A'. +// +// Unchanged on exit. +// +// TRANSB - CHARACTER*1. +// On entry, TRANSB specifies the form of op( B ) to be used in +// the matrix multiplication as follows: +// +// TRANSB = 'N' or 'n', op( B ) = B. +// +// TRANSB = 'T' or 't', op( B ) = B'. +// +// TRANSB = 'C' or 'c', op( B ) = B'. +// +// Unchanged on exit. +// +// M - INTEGER. +// On entry, M specifies the number of rows of the matrix +// op( A ) and of the matrix C. M must be at least zero. +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the number of columns of the matrix +// op( B ) and the number of columns of the matrix C. N must be +// at least zero. +// Unchanged on exit. +// +// K - INTEGER. +// On entry, K specifies the number of columns of the matrix +// op( A ) and the number of rows of the matrix op( B ). K must +// be at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is +// k when TRANSA = 'N' or 'n', and is m otherwise. +// Before entry with TRANSA = 'N' or 'n', the leading m by k +// part of the array A must contain the matrix A, otherwise +// the leading k by m part of the array A must contain the +// matrix A. +// Unchanged on exit. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. When TRANSA = 'N' or 'n' then +// LDA must be at least max( 1, m ), otherwise LDA must be at +// least max( 1, k ). +// Unchanged on exit. +// +// B - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is +// n when TRANSB = 'N' or 'n', and is k otherwise. +// Before entry with TRANSB = 'N' or 'n', the leading k by n +// part of the array B must contain the matrix B, otherwise +// the leading n by k part of the array B must contain the +// matrix B. +// Unchanged on exit. +// +// LDB - INTEGER. +// On entry, LDB specifies the first dimension of B as declared +// in the calling (sub) program. When TRANSB = 'N' or 'n' then +// LDB must be at least max( 1, k ), otherwise LDB must be at +// least max( 1, n ). +// Unchanged on exit. +// +// BETA - DOUBLE PRECISION. +// On entry, BETA specifies the scalar beta. When BETA is +// supplied as zero then C need not be set on input. +// Unchanged on exit. +// +// C - DOUBLE PRECISION array of DIMENSION ( LDC, n ). +// Before entry, the leading m by n part of the array C must +// contain the matrix C, except when beta is zero, in which +// case C need not be set on entry. +// On exit, the array C is overwritten by the m by n matrix +// ( alpha*op( A )*op( B ) + beta*C ). +// +// LDC - INTEGER. +// On entry, LDC specifies the first dimension of C as declared +// in the calling (sub) program. LDC must be at least +// max( 1, m ). +// Unchanged on exit. +// ------------------------------------------------------------ + +void dtrmm(char SIDE, char UPLO, char TRANSA, char DIAG, int M, int N, + double ALPHA, double *A, int LDA, double *B, int LDB); +// ------------------------------------------------------------ +// SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, +// B, LDB ) +// .. Scalar Arguments .. +// CHARACTER*1 SIDE, UPLO, TRANSA, DIAG +// INTEGER M, N, LDA, LDB +// DOUBLE PRECISION ALPHA +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ), B( LDB, * ) +// .. +// +// Purpose +// ======= +// +// DTRMM performs one of the matrix-matrix operations +// +// B := alpha*op( A )*B, or B := alpha*B*op( A ), +// +// where alpha is a scalar, B is an m by n matrix, A is a unit, or +// non-unit, upper or lower triangular matrix and op( A ) is one of +// +// op( A ) = A or op( A ) = A'. +// +// Parameters +// ========== +// +// SIDE - CHARACTER*1. +// On entry, SIDE specifies whether op( A ) multiplies B from +// the left or right as follows: +// +// SIDE = 'L' or 'l' B := alpha*op( A )*B. +// +// SIDE = 'R' or 'r' B := alpha*B*op( A ). +// +// Unchanged on exit. +// +// UPLO - CHARACTER*1. +// On entry, UPLO specifies whether the matrix A is an upper or +// lower triangular matrix as follows: +// +// UPLO = 'U' or 'u' A is an upper triangular matrix. +// +// UPLO = 'L' or 'l' A is a lower triangular matrix. +// +// Unchanged on exit. +// +// TRANSA - CHARACTER*1. +// On entry, TRANSA specifies the form of op( A ) to be used in +// the matrix multiplication as follows: +// +// TRANSA = 'N' or 'n' op( A ) = A. +// +// TRANSA = 'T' or 't' op( A ) = A'. +// +// TRANSA = 'C' or 'c' op( A ) = A'. +// +// Unchanged on exit. +// +// DIAG - CHARACTER*1. +// On entry, DIAG specifies whether or not A is unit triangular +// as follows: +// +// DIAG = 'U' or 'u' A is assumed to be unit triangular. +// +// DIAG = 'N' or 'n' A is not assumed to be unit +// triangular. +// +// Unchanged on exit. +// +// M - INTEGER. +// On entry, M specifies the number of rows of B. M must be at +// least zero. +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the number of columns of B. N must be +// at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. When alpha is +// zero then A is not referenced and B need not be set before +// entry. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m +// when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. +// Before entry with UPLO = 'U' or 'u', the leading k by k +// upper triangular part of the array A must contain the upper +// triangular matrix and the strictly lower triangular part of +// A is not referenced. +// Before entry with UPLO = 'L' or 'l', the leading k by k +// lower triangular part of the array A must contain the lower +// triangular matrix and the strictly upper triangular part of +// A is not referenced. +// Note that when DIAG = 'U' or 'u', the diagonal elements of +// A are not referenced either, but are assumed to be unity. +// Unchanged on exit. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. When SIDE = 'L' or 'l' then +// LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +// then LDA must be at least max( 1, n ). +// Unchanged on exit. +// +// B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). +// Before entry, the leading m by n part of the array B must +// contain the matrix B, and on exit is overwritten by the +// transformed matrix. +// +// LDB - INTEGER. +// On entry, LDB specifies the first dimension of B as declared +// in the calling (sub) program. LDB must be at least +// max( 1, m ). +// Unchanged on exit. +// ------------------------------------------------------------ + +void dtrsm(char SIDE, char UPLO, char TRANSA, char DIAG, int M, int N, + double ALPHA, double *A, int LDA, double *B, int LDB); +// ------------------------------------------------------------ +// SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, +// $ B, LDB ) +// .. Scalar Arguments .. +// CHARACTER*1 SIDE, UPLO, TRANSA, DIAG +// INTEGER M, N, LDA, LDB +// DOUBLE PRECISION ALPHA +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ), B( LDB, * ) +// .. +// +// Purpose +// ======= +// +// DTRSM solves one of the matrix equations +// +// op( A )*X = alpha*B, or X*op( A ) = alpha*B, +// +// where alpha is a scalar, X and B are m by n matrices, A is a unit, or +// non-unit, upper or lower triangular matrix and op( A ) is one of +// +// op( A ) = A or op( A ) = A'. +// +// The matrix X is overwritten on B. +// +// Parameters +// ========== +// +// SIDE - CHARACTER*1. +// On entry, SIDE specifies whether op( A ) appears on the left +// or right of X as follows: +// +// SIDE = 'L' or 'l' op( A )*X = alpha*B. +// +// SIDE = 'R' or 'r' X*op( A ) = alpha*B. +// +// Unchanged on exit. +// +// UPLO - CHARACTER*1. +// On entry, UPLO specifies whether the matrix A is an upper or +// lower triangular matrix as follows: +// +// UPLO = 'U' or 'u' A is an upper triangular matrix. +// +// UPLO = 'L' or 'l' A is a lower triangular matrix. +// +// Unchanged on exit. +// +// TRANSA - CHARACTER*1. +// On entry, TRANSA specifies the form of op( A ) to be used in +// the matrix multiplication as follows: +// +// TRANSA = 'N' or 'n' op( A ) = A. +// +// TRANSA = 'T' or 't' op( A ) = A'. +// +// TRANSA = 'C' or 'c' op( A ) = A'. +// +// Unchanged on exit. +// +// DIAG - CHARACTER*1. +// On entry, DIAG specifies whether or not A is unit triangular +// as follows: +// +// DIAG = 'U' or 'u' A is assumed to be unit triangular. +// +// DIAG = 'N' or 'n' A is not assumed to be unit +// triangular. +// +// Unchanged on exit. +// +// M - INTEGER. +// On entry, M specifies the number of rows of B. M must be at +// least zero. +// Unchanged on exit. +// +// N - INTEGER. +// On entry, N specifies the number of columns of B. N must be +// at least zero. +// Unchanged on exit. +// +// ALPHA - DOUBLE PRECISION. +// On entry, ALPHA specifies the scalar alpha. When alpha is +// zero then A is not referenced and B need not be set before +// entry. +// Unchanged on exit. +// +// A - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m +// when SIDE = 'L' or 'l' and is n when SIDE = 'R' or 'r'. +// Before entry with UPLO = 'U' or 'u', the leading k by k +// upper triangular part of the array A must contain the upper +// triangular matrix and the strictly lower triangular part of +// A is not referenced. +// Before entry with UPLO = 'L' or 'l', the leading k by k +// lower triangular part of the array A must contain the lower +// triangular matrix and the strictly upper triangular part of +// A is not referenced. +// Note that when DIAG = 'U' or 'u', the diagonal elements of +// A are not referenced either, but are assumed to be unity. +// Unchanged on exit. +// +// LDA - INTEGER. +// On entry, LDA specifies the first dimension of A as declared +// in the calling (sub) program. When SIDE = 'L' or 'l' then +// LDA must be at least max( 1, m ), when SIDE = 'R' or 'r' +// then LDA must be at least max( 1, n ). +// Unchanged on exit. +// +// B - DOUBLE PRECISION array of DIMENSION ( LDB, n ). +// Before entry, the leading m by n part of the array B must +// contain the right-hand side matrix B, and on exit is +// overwritten by the solution matrix X. +// +// LDB - INTEGER. +// On entry, LDB specifies the first dimension of B as declared +// in the calling (sub) program. LDB must be at least +// max( 1, m ). +// Unchanged on exit. + +// ------------------------------------------------------------ +int dgetri(int N, double *A, int LDA, int *IPIV, double *WORK, int LWORK); + +// ------------------------------------------------------------ +// SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) +// +// -- LAPACK routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// June 30, 1999 +// +// .. Scalar Arguments .. +// INTEGER INFO, LDA, LWORK, N +// .. +// .. Array Arguments .. +// INTEGER IPIV( * ) +// DOUBLE PRECISION A( LDA, * ), WORK( * ) +// .. +// +// Purpose +// ======= +// +// DGETRI computes the inverse of a matrix using the LU factorization +// computed by DGETRF. +// +// This method inverts U and then computes inv(A) by solving the system +// inv(A)*L = inv(U) for inv(A). +// +// Arguments +// ========= +// +// N (input) INTEGER +// The order of the matrix A. N >= 0. +// +// A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +// On entry, the factors L and U from the factorization +// A = P*L*U as computed by DGETRF. +// On exit, if INFO = 0, the inverse of the original matrix A. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,N). +// +// IPIV (input) INTEGER array, dimension (N) +// The pivot indices from DGETRF; for 1<=i<=N, row i of the +// matrix was interchanged with row IPIV(i). +// +// WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) +// On exit, if INFO=0, then WORK(1) returns the optimal LWORK. +// +// LWORK (input) INTEGER +// The dimension of the array WORK. LWORK >= max(1,N). +// For optimal performance LWORK >= N*NB, where NB is +// the optimal blocksize returned by ILAENV. +// +// If LWORK = -1, then a workspace query is assumed; the routine +// only calculates the optimal size of the WORK array, returns +// this value as the first entry of the WORK array, and no error +// message related to LWORK is issued by XERBLA. +// +// INFO (output) INTEGER +// = 0: successful exit +// < 0: if INFO = -i, the i-th argument had an illegal value +// > 0: if INFO = i, U(i,i) is exactly zero; the matrix is +// singular and its inverse could not be computed. + +// ------------------------------------------------------------ +int dpotri(char UPLO, int N, double *A, int LDA); + +// ------------------------------------------------------------ +// SUBROUTINE DPOTRI( UPLO, N, A, LDA, INFO ) +// +// -- LAPACK routine (version 3.0) -- +// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., +// Courant Institute, Argonne National Lab, and Rice University +// March 31, 1993 +// +// .. Scalar Arguments .. +// CHARACTER UPLO +// INTEGER INFO, LDA, N +// .. +// .. Array Arguments .. +// DOUBLE PRECISION A( LDA, * ) +// .. +// +// Purpose +// ======= +// +// DPOTRI computes the inverse of a real symmetric positive definite +// matrix A using the Cholesky factorization A = U**T*U or A = L*L**T +// computed by DPOTRF. +// +// Arguments +// ========= +// +// UPLO (input) CHARACTER*1 +// = 'U': Upper triangle of A is stored; +// = 'L': Lower triangle of A is stored. +// +// N (input) INTEGER +// The order of the matrix A. N >= 0. +// +// A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +// On entry, the triangular factor U or L from the Cholesky +// factorization A = U**T*U or A = L*L**T, as computed by +// DPOTRF. +// On exit, the upper or lower triangle of the (symmetric) +// inverse of A, overwriting the input factor U or L. +// +// LDA (input) INTEGER +// The leading dimension of the array A. LDA >= max(1,N). +// +// INFO (output) INTEGER +// = 0: successful exit +// < 0: if INFO = -i, the i-th argument had an illegal value +// > 0: if INFO = i, the (i,i) element of the factor U or L is +// zero, and the inverse could not be computed. + +#endif // LAPACK_CWRAP_H diff --git a/matte.c b/matte.c new file mode 100755 index 0000000..25cd3da --- /dev/null +++ b/matte.c @@ -0,0 +1,977 @@ +// ------------------------------------------------------------ +// matte.c +// +// Natural image matting using bayesian statistics +// From: +// +// +// +// 19.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#include +#include +#include +#include +#include +#include "lapack_cwrap.h" +#include "mlio.h" +#include "mutil.h" +#include "clob.h" +#include "mtypes.h" +#include "imageio.h" +#include "pixel.h" +#include "tictoc.h" +#include "matte.h" + +// ------------------------------------------------------------ +//#define MATTE_DEBUG +//#define WITH_COND_COMPUTING +#define USE_SLIDING_WINDOW +#define WITH_SETSIZE_BY_MAJOR +//#define WITH_ANIMATION + + +// ------------------------------------------------------------ +#ifdef USE_SLIDING_WINDOW +#undef WITH_SETSIZE_BY_MAJOR +#endif + +#define SIGMA_C 0.3 +#define SIGMA_GF 8 +#define SQRT_2PI 2.506628274631000502415765284811 +#define LAMBDA_THRESH 0.6 +#define MAX_COLORS 16 +#define RADIUS_MIN 4.0 +#define RADIUS_MAX 50 +#define WEIGHTS_MIN 0.001 +#define ALPHA_CONV_THRESH 0.001 +#define ALPHA_TAU 0.33 +#define SOLVER_MAX_ITER 30 +#define SOLVER_EPS_COND 1E-12 +#define REPORT_INTERVAL 1 +// ------------------------------------------------------------ +void AlphaInit(set_data_t *pObj, double *pAlpha, UINT32 type) +{ + int i; + pixel_t *pPixel = pObj->pPixel; + + if (type == pixel_state_foreground) + { + for (i=0; i < pObj->num_pixel; i++) + { + switch(pPixel[i].state) + { + case pixel_state_background: + pAlpha[i] = 0.0; + break; + + case pixel_state_foreground: + pAlpha[i] = 1.0; + break; + + case pixel_state_unknown: + pAlpha[i] = 0.5; + break; + } + } + } + else + { + + for (i=0; i < pObj->num_pixel; i++) + { + switch(pPixel[i].state) + { + case pixel_state_background: + pAlpha[i] = 1.0; + break; + + case pixel_state_foreground: + pAlpha[i] = 0.0; + break; + + case pixel_state_unknown: + pAlpha[i] = 0.5; + break; + } + } + } +} + +// ------------------------------------------------------------ +void MatteInit(matte_t *pObj, image_t *pComposite, image_t *pTrimap, char *name, matte_settings_t *pSettings) +{ + double k_norm[2]; + + int result, i, s; + tictoc_t tictoc; + char filename[1024]; + image_t im_dist_f, im_dist_b; + int num_row, num_col, num_unknown; + pxl_id_t t_f, t_b; + pixel_t *pPixel; + + color_vector_t *pColor_dist_f, *pColor_dist_b; + + num_row = pComposite->num_row; + num_col = pComposite->num_col; + + pObj->pWork = (double*)malloc(num_row*num_col*sizeof(double)); + + // Copy settings + memcpy(&pObj->settings, pSettings, sizeof(matte_settings_t)); + strcpy(pObj->name, name); + + SetInit(&pObj->comp_set, num_row, num_col); + + printf("Evaluating Trimap data..."); + Tic(&tictoc); + PixelStateInit(&pObj->comp_set, pTrimap); + ImageFree(pTrimap); + + // Init Alpha map + ImageInit(&pObj->im_alpha_f, num_row, num_col, 1); + AlphaInit(&pObj->comp_set, pObj->im_alpha_f.pColor_data, pixel_state_foreground); + + // Init 1-Alpha map + ImageInit(&pObj->im_alpha_b, num_row, num_col, 1); + AlphaInit(&pObj->comp_set, pObj->im_alpha_b.pColor_data, pixel_state_background); + + // Init probability map + ImageInit(&pObj->im_P, num_row, num_col, 1); + + // Get number of unknown pixels + num_unknown = GetSetByState(&pObj->comp_set, NULL, pixel_state_unknown); + pObj->pCalcSet = (pxl_id_t*)malloc(num_unknown*sizeof(pxl_id_t)); + + printf("done\n"); + Toc(&tictoc); + + // Create fore- and background images + ImageInit(&pObj->im_f, num_row, num_col, 3); + ImageInit(&pObj->im_b, num_row, num_col, 3); + ImageCopy(pComposite, &pObj->im_f); + ImageCopy(pComposite, &pObj->im_b); + + sprintf(filename, "%s_calcset.dat",name); + result = LoadCalcSet(filename, &pObj->comp_set, pObj->pCalcSet, num_unknown); + + // If loading of calcset fails - compute a new one + if (result < 0) + { + printf("Computing distance map..."); + Tic(&tictoc); + + // Get unknown pixels + num_unknown = GetSetByState(&pObj->comp_set, pObj->pCalcSet, pixel_state_unknown); + + // Compute distance map for unknown pixel + ComputeDistanceMap(&pObj->comp_set, pObj->pCalcSet, num_unknown); + + printf("done\n"); + Toc(&tictoc); + + printf("Computing calc set..."); + Tic(&tictoc); + + // Compute calcset base on distance map for unknown pixel + ComputeCalcSet(&pObj->comp_set, pObj->pCalcSet, num_unknown); + SaveCalcSet(filename, &pObj->comp_set, pObj->pCalcSet, num_unknown); + + printf("done\n"); + Toc(&tictoc); + + } + + ImageInit(&im_dist_f, num_row, num_col, 3); + ImageInit(&im_dist_b, num_row, num_col, 3); + + printf("Computing images from distance maps for %d unknown pixels...",num_unknown); + Tic(&tictoc); + + k_norm[0] = 1.0/num_col; + k_norm[1] = 1.0/num_row; + + pPixel = pObj->comp_set.pPixel; + pColor_dist_f = (color_vector_t*)im_dist_f.pColor_data; + pColor_dist_b = (color_vector_t*)im_dist_b.pColor_data; + + for (i=0; i < num_unknown; i++) + { + s = pObj->pCalcSet[i]; + t_f = pPixel[s].id_closest_f; + t_b = pPixel[s].id_closest_b; + + pColor_dist_f[s][0] = k_norm[0]*(double)pPixel[t_f].raster_coord[0]; + pColor_dist_f[s][2] = k_norm[1]*(double)pPixel[t_f].raster_coord[1]; + + pColor_dist_b[s][0] = k_norm[0]*(double)pPixel[t_b].raster_coord[0]; + pColor_dist_b[s][2] = k_norm[1]*(double)pPixel[t_b].raster_coord[1]; + } + + printf("done\n"); + Toc(&tictoc); + + sprintf(filename, "%s_dist_f.tif",name); + result = ImageSaveTiff(&im_dist_f, NULL, filename); + sprintf(filename, "%s_dist_b.tif",name); + result = ImageSaveTiff(&im_dist_b, NULL, filename); + + ImageFree(&im_dist_f); + ImageFree(&im_dist_b); + +} + +// ------------------------------------------------------------ +void MatteFree(matte_t *pObj) +{ + if (!pObj) + return; + + if (pObj->pWork) + free(pObj->pWork); + + if (pObj->pCalcSet) + free(pObj->pCalcSet); + + SetFree(&pObj->comp_set); +} + +// ------------------------------------------------------------ +int MatteGetSet(set_data_t *pObj, pxl_id_t *pPid, pxl_id_t pid, double *pAlpha, double *radius, UINT32 types) +{ + + int i, s, num_found; + pos_t center; + pixel_t *pPixel = pObj->pPixel; + double err, errlast, errslope, mean_weight; + double diff[2], a, d, s2_i, d_offset, sum_weight, w; + + s2_i = 1.0/(SIGMA_GF*SIGMA_GF); + + memcpy(center, pPixel[pid].raster_coord, sizeof(pos_t)); + + err = 0; + errlast = 0; + errslope = 1; + num_found = 0; + // A = r^2*pi + // r = sqrt(A/pi) + + do + { + err = 1 - num_found; + + if (err < 0) + *radius -= sqrt(-err/M_PI); + else + *radius += sqrt(err/M_PI); + + num_found = GetCircleSet(pObj, pPid, center, *radius, types); + + errslope = err - errlast; + errlast = err; + + } while (num_found == 0); + + sum_weight = 0; + do + { + + sum_weight = 0; + for (i=0; i < num_found; i++) + { + s = pPid[i]; + diff[0] = (double)(center[0] - pPixel[s].raster_coord[0]); + diff[1] = (double)(center[1] - pPixel[s].raster_coord[1]); + d = dnrm2(2, diff, 1); + a = pAlpha[s]; +// w = a*a/(d*d); + w = a*a*exp(-s2_i*d*d)*(1.0/SQRT_2PI); + sum_weight += w; + } + mean_weight = sum_weight/num_found; + + if (sum_weight < WEIGHTS_MIN) + { + *radius += 1; + num_found = GetCircleSet(pObj, pPid, center, *radius, types); + } + + } while ((sum_weight < WEIGHTS_MIN) && (*radius < RADIUS_MAX)); + +// printf("numfound = %d\n",num_found); + return num_found; +} + +// ------------------------------------------------------------ +int MatteGetColor(pxl_id_t *pPid, color_vector_t *pSrc, color_vector_t *pDst, int num_pixel) +{ + int i, s; + + for(i=0; i < num_pixel; i++) + { + s = pPid[i]; + pDst[i][0] = pSrc[s][0]; + pDst[i][1] = pSrc[s][1]; + pDst[i][2] = pSrc[s][2]; + } + + return i; +} + +// ------------------------------------------------------------ +int MatteGetData(pxl_id_t *pPid, double *pSrc, double *pDst, int num_pixel) +{ + int i, s; + + for(i=0; i < num_pixel; i++) + { + s = pPid[i]; + pDst[i] = pSrc[s]; + } + + return i; +} + +// ------------------------------------------------------------ +int MatteComputeWeights(set_data_t *pObj, pxl_id_t pid_u, pxl_id_t *pPid, double *pAlpha, double *pWeight, int num_pixel, double radius, UINT32 type) +{ + int i, s; + pos_t center; + double diff[2], a, d, s2_i, d_offset, sum_weight, alpha_mean; + + pixel_t *pPixel = pObj->pPixel; + + memcpy(center, pPixel[pid_u].raster_coord, sizeof(pos_t)); + +#ifndef USE_SLIDING_WINDOW + if (type == pixel_state_foreground) + { + d_offset = pPixel[pid_u].dist_closest_f; + } + else + { + d_offset = pPixel[pid_u].dist_closest_b; + } + + d_offset = smax(d_offset-radius, 0); +#else + d_offset = 0; +#endif + + s2_i = 1.0/(SIGMA_GF*SIGMA_GF); + +// d = exp(-s2_i*radius*radius)*(1.0/SQRT_2PI); +// printf("r=%f, min(d)=%f\n",radius,d); + + sum_weight = 0; + alpha_mean = 0; + for (i=0; i < num_pixel; i++) + { + s = pPid[i]; + diff[0] = (double)(center[0] - pPixel[s].raster_coord[0]); + diff[1] = (double)(center[1] - pPixel[s].raster_coord[1]); + d = dnrm2(2, diff, 1); + d = smax(d - d_offset, 0); + a = pAlpha[s]; + pWeight[i] = a*a*exp(-s2_i*d*d)*(1.0/SQRT_2PI); +// pWeight[i] = 1.0; + sum_weight += pWeight[i]; + alpha_mean += a; + } + + alpha_mean /= num_pixel; + if (alpha_mean < SOLVER_EPS_COND) + fprintf(stderr, "ComputeWeights(): alpha_mean = %f\n", alpha_mean); + + if (sum_weight < SOLVER_EPS_COND) + fprintf(stderr, "ComputeWeights(): Sum_i(w_i) = %f\n", sum_weight); + + return num_pixel; +} + +// ------------------------------------------------------------ +void ColorStatInit(color_stat_t *pObj, int max_pixel, int max_colors, UINT32 type) +{ + int info; + double lwork_svd, lwork_inv; + + pObj->type = type; + // Init SVD + // Query work size + info = dgesvd(/*JOBU*/'A', /*JOBVT*/'A', /*M*/3, /*N*/3, /*A*/NULL, /*LDA*/3, /*S*/NULL, + /*U*/NULL, /*LDU*/3, /*VT*/NULL, /*LDVT*/3, /*WORK*/&lwork_svd, /*LWORK*/-1); + + if (info != 0) fprintf (stderr, "Init SVD failure with error %d\n", info); + + // Matrix inversion + // Query work size + info = dgetri(3, NULL, 3, NULL, &lwork_inv, -1); + if (info != 0) fprintf (stderr, "Init dgetri failure with error %d\n", info); + + pObj->work_size = (int)smax(lwork_svd, lwork_inv); + pObj->work_size_svd = (int)lwork_svd; + pObj->work_size_inv = (int)lwork_inv; + pObj->pWork = (double*)malloc(pObj->work_size*sizeof(double)); + pObj->pColor_mean = (color_vector_t*)malloc(max_colors*sizeof(color_vector_t)); + pObj->pCov_inv = (color_matrix_t*)malloc(max_colors*sizeof(color_matrix_t)); +} + +// ------------------------------------------------------------ +void ColorStatFree(color_stat_t *pObj) +{ + + if (!pObj) + return; + + if (pObj->pWork) + free(pObj->pWork); + + if (pObj->pColor_mean) + free(pObj->pColor_mean); + + if (pObj->pCov_inv) + free(pObj->pCov_inv); + +} + +// ------------------------------------------------------------ +int ColorStatProcess(matte_t *pObj, color_stat_t *pStat, codebook_t *pCB, color_vector_t *pColor_set, int set_len, double *pWeight, double sigma_c) +{ + int i, j, num_colors, info, ipiv[3]={0}; + color_vector_t *pColor_data, *pColor_data_w, sv={0}; + color_vector_t color_mean={0}; + color_matrix_t cov={0}, cov1={0}, cov2={0}, cov_i={0}, U={0}, S={0}, VT={0}, T={0}; + double *pW, k_weight, sum_weight, s2, norm_cov, norm_cov_i, cond_cov; + + pW = (double*)malloc(set_len*sizeof(double)); + pColor_data = (color_vector_t*)malloc(set_len*sizeof(color_vector_t)); + + s2 = sigma_c*sigma_c; + + // Loop over number of colors found + pStat->num_colors = 0; + for (i=0; i < pCB->num_colors; i++) + { + num_colors = pCB->size_C[i]; + + // Get color values from Cluster i + MatteGetColor((pxl_id_t*)pCB->ppC[i], pColor_set, pColor_data, num_colors); + + // Get corresponding weights of pixels from Cluster i + MatteGetData((pxl_id_t*)pCB->ppC[i], pWeight, pW, num_colors); + + // k_weight = 1.0/sum_i(pW[i]) + sum_weight = dasum(num_colors, pW, 1); + if (sum_weight == 0) + { +// for(j=0; j < num_colors; j++) +// pW[j] = 1.0; + +// sum_weight = num_colors; + continue; + } + + k_weight = 1.0/sum_weight; + if (isinf(k_weight)) + continue; + + // color_mean := k_weight*[sum_i(pW[i]*pColorData[i])] + memset(color_mean, 0, sizeof(color_vector_t)); + dgemv('N', 3, num_colors, k_weight, (double*)pColor_data, 3, (double*)pW, 1, 0.0, (double*)color_mean, 1); + + // Remove mean from color data + for (j=0; j < num_colors; j++) + { + pColor_data[j][0] -= color_mean[0]; + pColor_data[j][1] -= color_mean[1]; + pColor_data[j][2] -= color_mean[2]; + } + + // cov = k_weight*pColor_data*pColor_data' + memset(cov, 0, sizeof(color_matrix_t)); + for (j=0; j < num_colors; j++) + { + dsyr('U', 3, pW[j], (double*)pColor_data[j], 1, (double*)cov, 3); + } + + cov[0][1] = cov[1][0]; + cov[0][2] = cov[2][0]; + cov[1][2] = cov[2][1]; + + // Scaling of cov + for (j=0; j < 3; j++) + { + cov[j][0] *= k_weight; + cov[j][1] *= k_weight; + cov[j][2] *= k_weight; + } + + // Do the Singular Value Decomposition + memcpy(cov1, cov, sizeof(color_matrix_t)); + info = dgesvd(/*JOBU*/'A', /*JOBVT*/'A', /*M*/3, /*N*/3, /*A*/(double*)cov1, /*LDA*/3, + /*S*/(double*)sv, /*U*/(double*)U, /*LDU*/3, + /*VT*/(double*)VT, /*LDVT*/3, + /*WORK*/(double*)pStat->pWork, /*LWORK*/pStat->work_size_svd); + + if (info != 0) fprintf (stderr, "SVD failure with error %d\n", info); + + // Add camera variance to singular value + sv[0] += s2; + sv[1] += s2; + sv[2] += s2; + + // Compose new covariance matrix + MakeDiag((double*)S, sv, 3); + memset(T, 0, 3*sizeof(color_vector_t)); + memset(cov2, 0, 3*sizeof(color_vector_t)); + + dgemm('N', 'N', 3, 3, 3, 1.0, (double*)U, 3, (double*)S, 3, 0.0, (double*)T, 3); + dgemm('N', 'N', 3, 3, 3, 1.0, (double*)T, 3, (double*)VT, 3, 0.0, (double*)cov2, 3); + +//#ifdef WITH_COND_COMPUTING + norm_cov = dnrm(3, 3, (double*)cov2, '1'); // 1-norm +//#endif + // Cov_inf = cov^-1; + Eye((double*)cov_i, 1.0, 3); + info = dgesv(3, 3, (double*)*cov2, 3, ipiv, (double*)cov_i, 3); + if (info != 0) + { + fprintf (stderr, "Matrix inversion failure with error %d\n", info); + } + +//#ifdef WITH_COND_COMPUTING + norm_cov_i = dnrm(3, 3, (double*)cov_i, '1'); // 1-norm + cond_cov = 1.0/(norm_cov*norm_cov_i); + + if (cond_cov < SOLVER_EPS_COND) + fprintf(stderr, "Condition(cov) = %f\n", cond_cov); +//#endif + + memcpy(pStat->pColor_mean[pStat->num_colors], color_mean, sizeof(color_vector_t)); + memcpy(pStat->pCov_inv[pStat->num_colors], cov_i, sizeof(color_matrix_t)); + pStat->num_colors++; + } + + free(pColor_data); + free(pW); + + return pStat->num_colors; +} + +// ------------------------------------------------------------ +void InsertMatrix(double *A, int lda, double *B, int ldb, int subid) +{ + int i, j, k, m; + + m = lda/ldb; + + k = ldb*(subid%m); + if (ldb*subid >= lda) + k += lda*ldb; + + for (i=0; i < ldb; i++) + for (j=0; j < ldb; j++) + A[j+i*lda+k] = B[i*ldb+j]; + +} + +// ------------------------------------------------------------ +int MatteEstimate(matte_sol_t *pObj, color_vector_t colorComp, color_stat_t *pStat_f, color_stat_t *pStat_b, double sigma_c) +{ + int i, j, k, iteration, info, ipiv[6], is_valid; + double *pW_f, *pW_b, alpha, b0, a0; + double a_s2, b_s2, ab_s2, s2_i, c[6], c_orig[6], c_last[6], x[6]; + color_matrix_t A[4] = {0}, A_orig[4], I, T, cov_f, cov_b; + color_vector_t colorMean_f, colorMean_b, t, c1, c2; + double a1, a2, delta_alpha, alpha_lt, L_f, L_b, L_sol; + double norm_A, rcond_A, dwork[24]; + int iwork[6]; + + // Find solution loop + s2_i = 1.0/(sigma_c*sigma_c); + + // ---------------------------- + // Loop begin + memset(pObj, 0, sizeof(matte_sol_t)); + pObj->L = -1.0/SOLVER_EPS_COND; + + if ((pStat_f->num_colors * pStat_b->num_colors) == 0) + fprintf(stderr, "MatteEstimate(): No solution possible!\n"); + + for (i=0; i < pStat_f->num_colors; i++) + { + memcpy(cov_f, pStat_f->pCov_inv[i], sizeof(color_matrix_t)); + memcpy(colorMean_f, pStat_f->pColor_mean[i], sizeof(color_vector_t)); + + for (j=0; j < pStat_b->num_colors; j++) + { + + memcpy(cov_b, pStat_b->pCov_inv[j], sizeof(color_matrix_t)); + memcpy(colorMean_b, pStat_b->pColor_mean[j], sizeof(color_vector_t)); + + delta_alpha = 1; + iteration = 0; + alpha = 0.5; + alpha_lt = 1.0; + while((delta_alpha > ALPHA_CONV_THRESH) && (iteration < SOLVER_MAX_ITER)) + { + memcpy(c_last, c, 6*sizeof(double)); + + is_valid = 0; + a0 = smin(smax(alpha,0),1); + b0 = 1.0-a0; + + // alpha/sigma_c^2 + a_s2 = a0*s2_i; + + // (1-alpha)/sigma_c^2 + b_s2 = b0*s2_i; + + // A = [[inv_C_F+(I*a0^2)/sigma_c^2 (I*a0*(1-a0))/sigma_c^2] + // ; [(I*a0*(1-a0))/sigma_c^2 inv_C_B+(I*(1-a0)^2)/sigma_c^2]]; + + memcpy(T, cov_f, sizeof(color_matrix_t)); + for (k=0; k < 3; k++) + { + T[k][k] += (a0*a_s2); + } + InsertMatrix((double*)A, 6, (double*)T, 3, 0); + + Eye((double*)T, a0*b_s2, 3); + InsertMatrix((double*)A, 6, (double*)T, 3, 1); + InsertMatrix((double*)A, 6, (double*)T, 3, 2); + + memcpy(T, cov_b, sizeof(color_matrix_t)); + for (k=0; k < 3; k++) + { + T[k][k] += (b0*b_s2); + } + InsertMatrix((double*)A, 6, (double*)T, 3, 3); + +#ifdef WITH_COND_COMPUTING + // Checking condition of A + norm_A = dnrm(6, 6, (double*)A, '1'); // 1-norm + info = dgecon('1', 6, (double*)A, 6, norm_A, &rcond_A, dwork, iwork); + if (info != 0) + { + fprintf (stderr, "Condition estimator failure with error %d\n", info); + break; + } + + if (rcond_A < SOLVER_EPS_COND) + { + fprintf(stderr, "Cond(A) = %g\n", rcond_A); + break; + } +#endif + // c = [[inv_C_F*color_mean_F+(color_C*a0)/sigma_c^2] + // ; [inv_C_B*color_mean_B+(color_C*(1-a0))/sigma_c^2]]; + memcpy(t, colorComp, sizeof(color_vector_t)); + dgemv('N', 3, 3, 1.0, (double*)cov_f, 3, (double*)colorMean_f, 1, a_s2, t, 1); + memcpy(&c[0], t, sizeof(color_vector_t)); + + memcpy(t, colorComp, sizeof(color_vector_t)); + dgemv('N', 3, 3, 1.0, (double*)cov_b, 3, (double*)colorMean_b, 1, b_s2, t, 1); + memcpy(&c[3], t, sizeof(color_vector_t)); + + // Solve + // x = A/c + memcpy(A_orig, A, 4*sizeof(color_matrix_t)); + memcpy(c_orig, c, 6*sizeof(double)); + info = dgesv(6, 1, (double*)A, 6, ipiv, (double*)c, 6); + if (info != 0) + { + fprintf (stderr, "Solver failure with error %d\n", info); +#ifdef MATTE_DEBUG + PrintMatrix("A=", (double*)A_orig, 6, 6); + PrintMatrix("c=", (double*)c_orig, 1, 6); + PrintMatrix("Cov_f=", (double*)cov_f, 3, 3); + PrintMatrix("Cov_b=", (double*)cov_b, 3, 3); + fprintf(stderr, "a0 = %f\n", a0); +#endif + // If solver fails just copy meanF into F and meanB into B + memcpy(&c[0], colorMean_f, sizeof(color_vector_t)); + memcpy(&c[3], colorMean_b, sizeof(color_vector_t)); + } + + // C - B + c1[0] = colorComp[0] - c[3]; + c1[1] = colorComp[1] - c[4]; + c1[2] = colorComp[2] - c[5]; + + // F - B + c2[0] = c[0] - c[3]; + c2[1] = c[1] - c[4]; + c2[2] = c[2] - c[5]; + + a1 = ddot(3, (double*)c1, 1, (double*)c2, 1); + a2 = dnrm2(3, c2, 1); + alpha = a1/(a2*a2); + alpha_lt = ALPHA_TAU*alpha + (1-ALPHA_TAU)*alpha_lt; + delta_alpha = fabs(1 - alpha/alpha_lt); + alpha = alpha_lt; + is_valid = 1; + iteration++; + } + + if (!is_valid && !iteration) + continue; + + // If we have a solution from previous iteration + if (!is_valid && iteration) + memcpy(c, c_last, 6*sizeof(double)); + + c1[0] = c[0] - colorMean_f[0]; + c1[1] = c[1] - colorMean_f[1]; + c1[2] = c[2] - colorMean_f[2]; + c2[0] = c[3] - colorMean_b[0]; + c2[1] = c[4] - colorMean_b[1]; + c2[2] = c[5] - colorMean_b[2]; + + // L_F = -(color_F-color_mean_F)'*inv_C_F*(color_F-color_mean_F)/2; + dgemv('N', 3, 3, 0.5, (double*)cov_f, 3, (double*)c1, 1, 0.0, t, 1); + L_f = -ddot(3, (double*)c1, 1, (double*)t, 1); + + // L_B = -(color_B-color_mean_B)'*inv_C_B*(color_B-color_mean_B)/2; + dgemv('N', 3, 3, 0.5, (double*)cov_b, 3, (double*)c2, 1, 0.0, t, 1); + L_b = -ddot(3, (double*)c2, 1, (double*)t, 1); + + c1[0] = colorComp[0] - (alpha*c[0] + (1-alpha)*c[3]); + c1[1] = colorComp[1] - (alpha*c[1] + (1-alpha)*c[4]); + c1[2] = colorComp[2] - (alpha*c[2] + (1-alpha)*c[5]); + + // Likelihood(k) = -norm(color_C-a*color_F-(1-a)*color_B)^2/sigma_c^2 + L_F + L_B; + L_sol = dnrm2(3, c1, 1); + L_sol = -(L_sol*L_sol)*s2_i + L_f + L_b; + + if (L_sol > pObj->L) + { + pObj->L = L_sol; + pObj->P = pow(10,L_sol); + pObj->alpha = alpha; + memcpy(pObj->color_f, &c[0], sizeof(color_vector_t)); + memcpy(pObj->color_b, &c[3], sizeof(color_vector_t)); + + pObj->alpha_c = smin(smax(alpha,0),1); + + pObj->color_f_c[0] = smin(smax(pObj->color_f[0],0),1); + pObj->color_f_c[1] = smin(smax(pObj->color_f[1],0),1); + pObj->color_f_c[2] = smin(smax(pObj->color_f[2],0),1); + + pObj->color_b_c[0] = smin(smax(pObj->color_b[0],0),1); + pObj->color_b_c[1] = smin(smax(pObj->color_b[1],0),1); + pObj->color_b_c[2] = smin(smax(pObj->color_b[2],0),1); + + pObj->solver_rounds = iteration; + } + +#ifdef MATTE_DEBUG + PrintMatrix("Mean_f=", (double*)colorMean_f, 1, 3); + PrintMatrix("color_f=", (double*)&c[0], 1, 3); + PrintMatrix("Mean_b=", (double*)colorMean_b, 1, 3); + PrintMatrix("color_b=", (double*)&c[3], 1, 3); + printf("Alpha=%g\n", alpha); + printf("P=%g\n", pow(10,L_sol)); + printf("\n"); +#endif + + } + } + // Loop end + // ---------------------------- + if (pObj->solver_rounds == 0) + fprintf(stderr, "MatteEstimate(): No solution found!\n"); +} + +// ------------------------------------------------------------ +int MatteProcess(matte_t *pObj) +{ + + int i, num_fore, num_back, num_colors_f, num_colors_b; + int max_npixel, num_computed; + pxl_id_t pid_u, pid_closest_f, pid_closest_b, *pPid_f, *pPid_b; + pxl_id_t *pCalcSet = pObj->pCalcSet; + set_data_t *pCompSet = &pObj->comp_set; + pixel_t *pPixel = pCompSet->pPixel; + color_vector_t *pColor_set; + double *pW_f, *pW_b, *pAlpha_f, *pAlpha_b, *pLikely, radius_f, radius_b; + color_vector_t colorComp, *pColorComp, *pColor_f, *pColor_b; + matte_sol_t sol; + char filename[1024]; + double progress, last_progress; + + last_progress = -REPORT_INTERVAL; + + if(pCompSet->num_fore > pCompSet->num_back) + + max_npixel = pCompSet->num_fore + pCompSet->num_unknown; + else + max_npixel = pCompSet->num_back + pCompSet->num_unknown; + + + colorquant_t clob_f, clob_b; + codebook_t cb_f, cb_b; + color_stat_t colorStat_f, colorStat_b; + + ColorQuant_Init(&clob_f, max_npixel, MAX_COLORS); + ColorQuant_Init(&clob_b, max_npixel, MAX_COLORS); + CodebookAlloc(&cb_f, max_npixel, MAX_COLORS); + CodebookAlloc(&cb_b, max_npixel, MAX_COLORS); + + pPid_f = (pxl_id_t*)malloc(max_npixel*sizeof(pxl_id_t)); + pPid_b = (pxl_id_t*)malloc(max_npixel*sizeof(pxl_id_t)); + + pW_f = (double*)malloc(max_npixel*sizeof(double)); + pW_b = (double*)malloc(max_npixel*sizeof(double)); + + pColor_set = (color_vector_t*)malloc(max_npixel*sizeof(color_vector_t)); + + ColorStatInit(&colorStat_f, max_npixel, MAX_COLORS, pixel_state_foreground); + ColorStatInit(&colorStat_b, max_npixel, MAX_COLORS, pixel_state_background); + + pColor_f = (color_vector_t*)pObj->im_f.pColor_data; + pColor_b = (color_vector_t*)pObj->im_b.pColor_data; + pColorComp = pColor_f; + pAlpha_f = (double*)pObj->im_alpha_f.pColor_data; + pAlpha_b = (double*)pObj->im_alpha_b.pColor_data; + pLikely = (double*)pObj->im_P.pColor_data; + + num_computed = 0; + for(i=0; i < pCompSet->num_unknown; i++) + { + + pid_u = pCalcSet[i]; + + if(pixel_state_unknown != pObj->comp_set.pPixel[pid_u].state) + { + fprintf(stderr,"Pixel #%d is not unknown\n", i); + } + + memcpy(colorComp, pColorComp[pid_u], sizeof(color_vector_t)); + + pid_closest_f = pPixel[pid_u].id_closest_f; + pid_closest_b = pPixel[pid_u].id_closest_b; + + // Get foreground set and quantize + radius_f = RADIUS_MIN; + +#ifdef USE_SLIDING_WINDOW + num_fore = MatteGetSet(pCompSet, pPid_f, pid_u, pAlpha_f, &radius_f, pixel_state_foreground | pixel_state_computed); +#else + num_fore = MatteGetSet(pCompSet, pPid_f, pid_closest_f, pAlpha_f, &radius_f, pixel_state_foreground | pixel_state_computed); +#endif + MatteGetColor(pPid_f, pColor_f, pColor_set, num_fore); + + num_colors_f = ColorQuant(&clob_f, &cb_f, pColor_set, 3, num_fore, MAX_COLORS, LAMBDA_THRESH); + + MatteComputeWeights(pCompSet, pid_u, pPid_f, pAlpha_f, pW_f, num_fore, radius_f, pixel_state_foreground); + ColorStatProcess(pObj, &colorStat_f, &cb_f, pColor_set, num_fore, pW_f, SIGMA_C); + + // Get background set and quantize + radius_b = RADIUS_MIN; + +#ifdef USE_SLIDING_WINDOW + num_back = MatteGetSet(pCompSet, pPid_b, pid_u, pAlpha_b, &radius_b, pixel_state_background | pixel_state_computed); +#else + num_back = MatteGetSet(pCompSet, pPid_b, pid_closest_b, pAlpha_b, &radius_b, pixel_state_background | pixel_state_computed); +#endif + MatteGetColor(pPid_b, pColor_b, pColor_set, num_back); + + num_colors_b = ColorQuant(&clob_b, &cb_b, pColor_set, 3, num_back, MAX_COLORS, LAMBDA_THRESH); + + MatteComputeWeights(pCompSet, pid_u, pPid_b, pAlpha_b, pW_b, num_back, radius_b, pixel_state_background); + ColorStatProcess(pObj, &colorStat_b, &cb_b, pColor_set, num_back, pW_b, SIGMA_C); + + // Find solution + MatteEstimate(&sol, colorComp, &colorStat_f, &colorStat_b, SIGMA_C); + +#ifdef MATTE_DEBUG + printf("\n------------------------------------------\n"); + printf("Unknown Pixel #%d\n", i); + PrintMatrix("color_c=", (double*)colorComp, 3, 1); + PrintMatrix("color_f=", (double*)sol.color_f, 3, 1); + PrintMatrix("color_b=", (double*)sol.color_b, 3, 1); + printf("Alpha= %f\n", sol.alpha); + printf("P_sol=%f\n", sol.P); + printf("Number of solver round=%d\n",sol.solver_rounds); + printf("Number of colors int foreground statistic=%d\n",colorStat_f.num_colors); + printf("Number of colors int background statistic=%d\n",colorStat_b.num_colors); +#endif + memcpy(pColor_f[pid_u], sol.color_f_c, sizeof(color_vector_t)); + memcpy(pColor_b[pid_u], sol.color_b_c, sizeof(color_vector_t)); + pAlpha_f[pid_u] = sol.alpha_c; + pAlpha_b[pid_u] = 1-sol.alpha_c; + pLikely[pid_u] = sol.P; + + // Label unknown pixel as computed + if (sol.solver_rounds > 0) + { + pObj->comp_set.pPixel[pid_u].state = pixel_state_computed; + num_computed++; + } + else + { + fprintf(stderr, "Unknown Pixel #%d(ID=%d) not computed.\n", i, pid_u); + } + + +#ifdef WITH_ANIMATION + sprintf(filename, "ani/%s_alpha_%5.5d.tif",pObj->name,i); + ImageSaveTiff(&pObj->im_alpha_f, NULL, filename); + sprintf(filename, "ani/%s_est_fore_%5.5d.tif",pObj->name, i); + ImageSaveTiff(&pObj->im_f, NULL, filename); + sprintf(filename, "ani/%s_est_back_%5.5d.tif",pObj->name, i); + ImageSaveTiff(&pObj->im_b, NULL, filename); +#endif + // Output progress + progress = 100*(double)num_computed/(pCompSet->num_unknown-1); + + if ((progress-last_progress) >= REPORT_INTERVAL) + { + printf("\r%d %% computed...", (int)progress); + last_progress = progress; + } + } + + + printf("\r%d %% computed. \n", (int)progress); + printf("%d of %d pixels computed.\n\n", num_computed, pCompSet->num_unknown); + + sprintf(filename, "%s_alpha.tif", pObj->name); + ImageSaveTiff(&pObj->im_alpha_f, NULL, filename); + sprintf(filename, "%s_fore.tif",pObj->name); + ImageSaveTiff(&pObj->im_f, pAlpha_f, filename); + sprintf(filename, "%s_back.tif",pObj->name); + ImageSaveTiff(&pObj->im_b, pAlpha_b, filename); + sprintf(filename, "%s_est_fore.tif",pObj->name); + ImageSaveTiff(&pObj->im_f, NULL, filename); + sprintf(filename, "%s_est_back.tif",pObj->name); + ImageSaveTiff(&pObj->im_b, NULL, filename); + sprintf(filename, "%s_p.tif",pObj->name); + ImageSaveTiff(&pObj->im_P, NULL, filename); + + ImageFree(&pObj->im_f); + ImageFree(&pObj->im_b); + ImageFree(&pObj->im_alpha_f); + ImageFree(&pObj->im_alpha_b); + ImageFree(&pObj->im_P); + + CodebookFree(&cb_f); + CodebookFree(&cb_b); + + ColorQuant_Free(&clob_f); + ColorQuant_Free(&clob_b); + + ColorStatFree(&colorStat_f); + ColorStatFree(&colorStat_b); + + free(pPid_f); + free(pPid_b); + free(pW_f); + free(pW_b); + free(pColor_set); + +} +// ------------------------------------------------------------ diff --git a/matte.h b/matte.h new file mode 100755 index 0000000..95c4a07 --- /dev/null +++ b/matte.h @@ -0,0 +1,66 @@ +// ------------------------------------------------------------ +// matte.h +// +// Reading writing tiff images +// +// 12.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef MATTE_H +#define MATTE_H + +#include "mtypes.h" +#include "imageio.h" +#include "pixel.h" + +// ------------------------------------------------------------ +typedef struct _matte_settings_t +{ + int debug; + int use_lab; + int use_computed_pixel; + int use_sliding_window_computing; // instead of clustering fore and back separately + int use_constant_color_background; + double sigma_c, sigma_gf; + int min_set_size; + +} matte_settings_t; + +typedef struct _matte_t +{ + char name[256]; + image_t im_f, im_b, im_alpha_f, im_alpha_b, im_P; + double *pL, *pWork; + set_data_t comp_set; + pxl_id_t *pCalcSet; + matte_settings_t settings; + +} matte_t; + +typedef struct _color_stat_t +{ + color_vector_t *pColor_mean; + color_matrix_t *pCov_inv; + double *pWork; + int work_size, work_size_svd, work_size_inv; + int num_colors; + UINT32 type; + +} color_stat_t; + +typedef struct _matte_sol_t +{ + double L, P; + double alpha, alpha_c; + color_vector_t color_f, color_b, color_f_c, color_b_c; + int solver_rounds; + +} matte_sol_t; + +// ------------------------------------------------------------ +void MatteInit(matte_t *pObj, image_t *pComposite, image_t *pTrimap, char *name, matte_settings_t *pSettings); +int MatteGetSet(set_data_t *pObj, pxl_id_t *pPid, pxl_id_t pid, double *pAlpha, double *radius, UINT32 types); +int MatteGetColor(pxl_id_t *pPid, color_vector_t *pSrc, color_vector_t *pDst, int num_pixel); +int MatteGetData(pxl_id_t *pPid, double *pSrc, double *pDst, int num_pixel); + +// ------------------------------------------------------------ +#endif // MATTE_H diff --git a/mlio.c b/mlio.c new file mode 100755 index 0000000..1c9fce3 --- /dev/null +++ b/mlio.c @@ -0,0 +1,82 @@ +// ----------------------------------------------------------------------------------------- +// mlio.c +// Input/output for Matlab matrices +// ----------------------------------------------------------------------------------------- +#include +#include +#include +#include "mlio.h" + +// ----------------------------------------------------------------------------------------- +int mlio_open(mlio_t *pObj, char *pName, unsigned mode) +{ + + pObj->io_mode = mode; + + if (mode == MLIO_READ) + { + pObj->pFile = fopen(pName, "rb"); + + // Read input file + if (!pObj->pFile) + { + fprintf(stderr,"Unable to open input file %s\n", pName); + return -1; + } + + // Header + fread(&pObj->hdr_file.version, 1, sizeof(unsigned), pObj->pFile); + fscanf(pObj->pFile, "%s", pObj->hdr_file.prec_typename); + fseek(pObj->pFile, 1, SEEK_CUR); + fread(&pObj->hdr_data, 1, sizeof(data_hdr_t), pObj->pFile); + printf("Rows=%d, Columns=%d, File size=%d bytes\n", pObj->hdr_data.num_row, pObj->hdr_data.num_col, pObj->hdr_data.size); + + } + else + { + pObj->pFile = fopen(pName, "wb"); + + // Write output file + if (!pObj->pFile) + { + fprintf(stderr,"Unable to open output file %s\n", pName); + return -1; + } + + } + return 0; + +} + +int mlio_close(mlio_t *pObj) +{ + return fclose(pObj->pFile); +} + +int mlio_read(mlio_t *pObj, double *pData) +{ + // Read data + fread(pData, pObj->hdr_data.num_row*pObj->hdr_data.num_col, sizeof(double), pObj->pFile); + + return 0; + +} + +int mlio_write(mlio_t *pObj, double *pData, unsigned num_rows, unsigned num_cols) +{ + + pObj->hdr_data.num_row = num_rows; + pObj->hdr_data.num_col = num_cols; + + // Header + fwrite(&pObj->hdr_file.version, 1, sizeof(unsigned), pObj->pFile); + fprintf(pObj->pFile, "%s\r", pObj->hdr_file.prec_typename); + fwrite(&pObj->hdr_data, 1, sizeof(data_hdr_t), pObj->pFile); + + // Write data + fwrite(pData, pObj->hdr_data.num_row*pObj->hdr_data.num_col, sizeof(double), pObj->pFile); + + return 0; +} + +// ----------------------------------------------------------------------------------------- diff --git a/mlio.h b/mlio.h new file mode 100755 index 0000000..2a41194 --- /dev/null +++ b/mlio.h @@ -0,0 +1,46 @@ +// ----------------------------------------------------------------------------------------- +// mlio.h +// Input/output for Matlab matrices +// 26.02.2005, J. Ahrensfeld +// +// ----------------------------------------------------------------------------------------- +#ifndef MLIO_H +#define MLIO_H + +// ----------------------------------------------------------------------------------------- +#define MLIO_READ 1 +#define MLIO_WRITE 2 + +// ----------------------------------------------------------------------------------------- +typedef struct _data_hdr_t +{ + unsigned num_row, num_col; + unsigned size; + +} data_hdr_t; + +typedef struct _file_hdr_t +{ + unsigned version; + char prec_typename[32]; + +} file_hdr_t; + +typedef struct _mlio_t +{ + FILE *pFile; + unsigned io_mode; + + file_hdr_t hdr_file; + data_hdr_t hdr_data; +} mlio_t; + +// ----------------------------------------------------------------------------------------- +int mlio_open(mlio_t *pObj, char *pName, unsigned mode); +int mlio_close(mlio_t *pObj); +int mlio_read(mlio_t *pObj, double *pData); +int mlio_write(mlio_t *pObj, double *pData, unsigned num_rows, unsigned num_cols); + +// ----------------------------------------------------------------------------------------- +#endif //MLIO_H + diff --git a/mtypes.h b/mtypes.h new file mode 100755 index 0000000..faec5f8 --- /dev/null +++ b/mtypes.h @@ -0,0 +1,34 @@ +// ------------------------------------------------------------ +// imageio.h +// +// Reading writing tiff images +// +// 12.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef MTYPES_H +#define MTYPES_H + +// ------------------------------------------------------------ +#define COMP_RED 0 +#define COMP_GREEN 1 +#define COMP_BLUE 2 +#define COMP_L 0 +#define COMP_A 1 +#define COMP_B 2 +#define COORD_X 0 +#define COORD_Y 1 + +typedef unsigned int UINT32; +typedef signed int INT32; +typedef unsigned short UINT16; +typedef signed short INT16; +typedef unsigned char UINT8; +typedef signed char INT8; + +typedef UINT32 pxl_id_t; + +typedef int pos_t[2]; + +// ------------------------------------------------------------ + +#endif // MTYPES_H diff --git a/mutil.c b/mutil.c new file mode 100755 index 0000000..6f9381a --- /dev/null +++ b/mutil.c @@ -0,0 +1,220 @@ +// ------------------------------------------------------------ +// mutil.c +// Matrix utilities +// 5.3.2005, J.Ahrensfeld +// +// ------------------------------------------------------------ +#include +#include +#include + +#include "lapack_cwrap.h" + +// ------------------------------------------------------------ +double lp_work_buf[1024] = {0}; + +// ------------------------------------------------------------ +void Eye(double *pI, double value, int N) +{ + int len; + + memset(pI, 0, N*N*sizeof(double)); + + len = 0; + + while(len < N*N) + { + pI[len] = value; + len += (N+1); + } + +} + +// ------------------------------------------------------------ +void RandMat(double *pA, int N) +{ + int i; + + for(i=0; i < N; i++) + pA[i] = 2.0*(0.5 - (double)rand()/RAND_MAX); +} + +// ------------------------------------------------------------ +void IncMat(double *pA, int nRow, int nCol, char order) +{ + int i, j, k; + + if (order == 'R') + { + k = 1; + for(i=0; i < nCol; i++) + for(j=0; j < nRow; j++) + pA[i*nRow+j] = (double)k++; + } + else + { + k = 1; + for(i=0; i < nRow; i++) + for(j=0; j < nCol; j++) + pA[j*nRow+i] = (double)k++; + } +} + +// ------------------------------------------------------------ +void CopyTrans(double *pA, double *pAT, int nRow, int nCol) +{ + int i, j; + + for(i=0; i < nRow; i++) + for(j=0; j < nCol; j++) + pAT[j+i*nCol] = pA[j*nRow+i]; + +} + +// ------------------------------------------------------------ +void GetDiag(double *pA, double *pY, int N) +{ + int i; + for (i=0; i < N; i++) + pY[i] = pA[i*N+i]; + +} + +// ------------------------------------------------------------ +void MakeDiag(double *pD, double *pX, int N) +{ + int len; + + memset(pD, 0, N*N*sizeof(double)); + + len = 0; + while(len < N*N) + { + pD[len] = *(pX++); + len += (N+1); + } + +} + +// ------------------------------------------------------------ +void PrintMatrix(char *pStr, double *A, int num_row, int num_col) +{ + int i, j; + double v; + + printf("%s\n",pStr); + for (i=0; i v2) + return v1; + + return v2; +} + +// ------------------------------------------------------------ +void vsum(double *pX, double *pY, int lda, int N) +{ + int i, j; + + for (j=0; j < lda; j++) + pY[j] = *(pX++); + + for(i=1; i < N; i++) + for (j=0; j < lda; j++) + pY[j] += *(pX++); + +} + +// ------------------------------------------------------------ +void GetPrincipleEigenVector(double *pA, double *pV, int lda) +{ + double max_x; + int max_i, i; + + max_x = pV[0]; + max_i = 0; + + for (i=1; i < lda; i++) + { + if(pV[i] > max_x) + { + max_x = pV[0]; + max_i = i; + } + } + + for (i=0; i < lda; i++) + pV[i] = pA[max_i*lda + i]; +} + +// ------------------------------------------------------------ +double dnrm(int N, int M, double *A, char mode) +{ + int i; + double maxval, val; + + maxval = 0; + if (mode == '1') + { + for (i=0; i < N*M; i += N) + { + val = dasum(N, &A[i], 1); + if (val > maxval) + maxval = val; + } + } + else + { + for (i=0; i < N; i++) + { + val = dasum(M, &A[i], N); + if (val > maxval) + maxval = val; + } + } + + return maxval; +} diff --git a/mutil.h b/mutil.h new file mode 100755 index 0000000..7e134a7 --- /dev/null +++ b/mutil.h @@ -0,0 +1,26 @@ +// ------------------------------------------------------------ +// mutil.h +// Matrix utilities +// 5.3.2005, J.Ahrensfeld +// +// ------------------------------------------------------------ +#ifndef MUTIL_H +#define MUTIL_H + +// ------------------------------------------------------------ +// ------------------------------------------------------------ +void RandMat(double *pA, int N); +void IncMat(double *pA, int nRow, int nCol, char order); +void CopyTrans(double *pA, double *pAT, int nRow, int nCol); +void Eye(double *pI, double value, int N); +void Diag(double *pA, double *pY, int N); +void MakeDiag(double *pD, double *pX, int N); +void PrintMatrix(char *pStr, double *A, int num_row, int num_col); +double smin(double v1, double v2); +double smax(double v1, double v2); +void vsum(double *pX, double *pY, int lda, int N); +void GetPrincipleEigenVector(double *pA, double *pV, int lda); +double dnrm(int N, int M, double *A, char mode); + +// ------------------------------------------------------------ +#endif //MUTIL_H diff --git a/pixel.c b/pixel.c new file mode 100755 index 0000000..575d2f5 --- /dev/null +++ b/pixel.c @@ -0,0 +1,413 @@ +// ------------------------------------------------------------ +// pixel.c +// +// Retrieving pixels in image and maintaining pixel attributes +// +// 13.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#include +#include +#include +#include +#include +#include "bintree.h" +#include "mutil.h" +#include "lapack_cwrap.h" +#include "clob.h" +#include "pixel.h" + +// ------------------------------------------------------------ +#define CALCSET_RADIUS_INITIAL 5 +#define CALCSET_RADIUS_INCREMENT 1.0 +#define CALCSET_DISTMIN_INITIAL 1E12 + +// ------------------------------------------------------------ +void SetInit(set_data_t *pObj, int num_row, int num_col) +{ + int i, x, y; + + pObj->num_row = num_row; + pObj->num_col = num_col; + pObj->num_pixel = num_row*num_col; + + pObj->pPixel = (pixel_t*)malloc(pObj->num_pixel*sizeof(pixel_t)); + + x = y = 0; + for (i=0; i < pObj->num_pixel; i++) + { + pObj->pPixel[i].id_this = i; + pObj->pPixel[i].state = 0; + pObj->pPixel[i].raster_coord[0] = x++; + pObj->pPixel[i].raster_coord[1] = y; + if (x >= num_col) + { + x = 0; + y++; + } + } +} + +// ------------------------------------------------------------ +void SetFree(set_data_t *pObj) +{ + if (!pObj) + return; + + if (pObj->pPixel) + free(pObj->pPixel); +} + +// ------------------------------------------------------------ +void PixelStateInit(set_data_t *pObj, image_t *pTrimap_image) +{ + int i, k; + double *pTrimap = (double*)pTrimap_image->pColor_data; + + k = pTrimap_image->num_ch; + + pObj->num_fore = 0; + pObj->num_back = 0; + pObj->num_unknown = 0; + + for (i=0; i < pObj->num_pixel; i++) + { + pObj->pPixel[i].state = pixel_state_background; + if(pTrimap[k*i] > 0.7) + { + pObj->pPixel[i].state = pixel_state_foreground; + pObj->num_fore++; + } + else + { + if(pTrimap[k*i] > 0.2) + { + pObj->pPixel[i].state = pixel_state_unknown; + pObj->num_unknown++; + } + } + } + pObj->num_back = pObj->num_pixel - (pObj->num_fore + pObj->num_unknown); +} + +// ------------------------------------------------------------ +int GetCropSet(set_data_t *pObj, int *pCrop, pos_t center, pos_t rect_size, UINT32 select_state) +{ + + int i, j, x, y, min_x, min_y, max_x, max_y; + + if (!(rect_size[COORD_X]%2)) + rect_size[COORD_X]++; + + if (!(rect_size[COORD_Y]%2)) + rect_size[COORD_Y]++; + + min_x = center[COORD_X]-(rect_size[COORD_X]-1)/2; + if (min_x < 0) + min_x = 0; + + max_x = center[COORD_X]+(rect_size[COORD_X]-1)/2; + if (max_x >= pObj->num_col) + max_x = pObj->num_col-1; + + min_y = center[COORD_Y]-(rect_size[COORD_Y]-1)/2; + if (min_y < 0) + min_y = 0; + + max_y = center[COORD_Y]+(rect_size[COORD_Y]-1)/2; + if (max_y >= pObj->num_row) + max_y = pObj->num_row-1; + + + i = 0; + if (select_state == pixel_state_all) + { + for (y=min_y; y <= max_y; y++) + { + for (x=min_x; x <= max_x; x++) + { + pCrop[i] = y*pObj->num_col + x; + i++; + } + } + } + else + { + for (y=min_y; y <= max_y; y++) + { + for (x=min_x; x <= max_x; x++) + { + j = y*pObj->num_col + x; + if (pObj->pPixel[j].state & select_state) + { + pCrop[i] = j; + i++; + } + } + } + } + + return i; +} + +// ------------------------------------------------------------ +int GetCircleSet(set_data_t *pObj, int *pCrop, pos_t center, double radius, UINT32 select_state) +{ + int i, j, num_found; + + double distance, diff[2]; + pos_t rect_size; + + rect_size[0] = 2*(int)radius+1; + rect_size[1] = 2*(int)radius+1; + + num_found = GetCropSet(pObj, pCrop, center, rect_size, select_state); + + j=0; + for (i=0; i < num_found; i++) + { + diff[0] = (double)(center[0] - pObj->pPixel[pCrop[i]].raster_coord[0]); + diff[1] = (double)(center[1] - pObj->pPixel[pCrop[i]].raster_coord[1]); + distance = dnrm2(2, diff, 1); + + if(distance < radius) + { + pCrop[j++] = pCrop[i]; + } + } + + return j; +} + +// ------------------------------------------------------------ +int GetSetByState(set_data_t *pObj, pxl_id_t *pSet, UINT32 select_state) +{ + + UINT32 i, j; + + j=0; + if (pSet) + { + for (i=0; i < pObj->num_pixel; i++) + { + + if (pObj->pPixel[i].state & select_state) + { + pSet[j] = pObj->pPixel[i].id_this; + j++; + } + } + } + else // Query only + { + for (i=0; i < pObj->num_pixel; i++) + { + + if (pObj->pPixel[i].state & select_state) + { + j++; + } + } + } + return j; +} + +// ------------------------------------------------------------ +// Calcset set stuff +// Creates optimal computing order for unknown pixels +// ------------------------------------------------------------ +int ComputeDistanceMap(set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown) +{ + int i, j, s, t, num_fore, num_back; + pxl_id_t *pSet_f, *pSet_b, id_min; + double radius, distance, distance_min, diff[2]; + pos_t curr_pos; + node_t *node; + bintree_t bst; + + void *data; + + pSet_f = (pxl_id_t*)malloc(pObj->num_fore*sizeof(pxl_id_t)); + pSet_b = (pxl_id_t*)malloc(pObj->num_back*sizeof(pxl_id_t)); + + for (i=0; i < num_unknown; i++) + { + s = pSet_unknown[i]; + curr_pos[0] = pObj->pPixel[s].raster_coord[0]; + curr_pos[1] = pObj->pPixel[s].raster_coord[1]; + + radius = CALCSET_RADIUS_INITIAL; + do + { + num_fore = GetCircleSet(pObj, pSet_f, curr_pos, radius, pixel_state_foreground); + radius += CALCSET_RADIUS_INCREMENT; + + } while (num_fore == 0); + + radius = CALCSET_RADIUS_INITIAL; + do + { + num_back = GetCircleSet(pObj, pSet_b, curr_pos, radius, pixel_state_background); + radius += CALCSET_RADIUS_INCREMENT; + + } while (num_back == 0); + + bintree_init(&bst, num_fore); + distance_min = CALCSET_DISTMIN_INITIAL; + for (j=0; j < num_fore; j++) + { + t = pSet_f[j]; + diff[0] = (double)(curr_pos[0] - pObj->pPixel[t].raster_coord[0]); + diff[1] = (double)(curr_pos[1] - pObj->pPixel[t].raster_coord[1]); + + distance = diff[0]*diff[0] + diff[1]*diff[1]; + + data = (void*)t; + bintree_insert(&bst, distance, data); + + } + node = bintree_GetNodeMin(&bst); + id_min = (int)node->pData; + distance_min = node->value; + pObj->pPixel[s].id_closest_f = id_min; + pObj->pPixel[s].dist_closest_f = sqrt(distance_min); + bintree_free(&bst); + + bintree_init(&bst, num_back); + distance_min = CALCSET_DISTMIN_INITIAL; + for (j=0; j < num_back; j++) + { + t = pSet_b[j]; + diff[0] = (double)(curr_pos[0] - pObj->pPixel[t].raster_coord[0]); + diff[1] = (double)(curr_pos[1] - pObj->pPixel[t].raster_coord[1]); + + distance = diff[0]*diff[0] + diff[1]*diff[1]; + + data = (void*)t; + bintree_insert(&bst, distance, data); + + } + node = bintree_GetNodeMin(&bst); + + id_min = (int)node->pData; + distance_min = node->value; + pObj->pPixel[s].id_closest_b = id_min; + pObj->pPixel[s].dist_closest_b = sqrt(distance_min); + bintree_free(&bst); + } + + free(pSet_f); + free(pSet_b); + + return num_unknown; +} + +// ------------------------------------------------------------ +int ComputeCalcSet(set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown) +{ + int i, s; + void *data; + + double v[2], distance; + bintree_t bst; + node_t *node; + + bintree_init(&bst, num_unknown); + + for (i=0; i < num_unknown; i++) + { + s = pSet_unknown[i]; + v[0] = pObj->pPixel[s].dist_closest_f; + v[1] = pObj->pPixel[s].dist_closest_b; + + if (v[0] < v[1]) + distance = v[0]; + else + distance = v[1]; + + data = (void*)s; + bintree_insert(&bst, distance, data); + } + + i=0; + do + { + node = bintree_GetNodeMin(&bst); + s = (int)node->pData; + pSet_unknown[i] = s; + i++; + + } while(bintree_node_remove(&bst, node)); + + bintree_free(&bst); + + return i; +} +// ------------------------------------------------------------ +int SaveCalcSet(char *filename, set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown) +{ + int i; + pxl_id_t s; + FILE *pFile; + + pFile = fopen(filename, "wb"); + if (!pFile) + { + fprintf(stderr, "Unable to open %s\n",filename); + return -1; + } + + fwrite(&num_unknown, sizeof(int), 1, pFile); + fwrite(pSet_unknown, sizeof(pxl_id_t), num_unknown, pFile); + + for(i=0; i < num_unknown; i++) + { + s = pSet_unknown[i]; + fwrite(&pObj->pPixel[s], sizeof(pixel_t), 1, pFile); + } + + + fclose(pFile); + + return 0; +} + +// ------------------------------------------------------------ +int LoadCalcSet(char *filename, set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown) +{ + int i, N; + pxl_id_t s; + FILE *pFile; + + pFile = fopen(filename, "rb"); + if (!pFile) + { + fprintf(stderr, "Unable to open %s\n",filename); + return -1; + } + + fread(&N, sizeof(int), 1, pFile); + if (N != num_unknown) + { + fclose(pFile); + return -1; + } + + fread(pSet_unknown, sizeof(pxl_id_t), num_unknown, pFile); + + for(i=0; i < num_unknown; i++) + { + s = pSet_unknown[i]; + fread(&pObj->pPixel[s], sizeof(pixel_t), 1, pFile); + } + + fclose(pFile); + + for(i=0; i < num_unknown; i++) + { + s = pSet_unknown[i]; + if(pObj->pPixel[s].state != pixel_state_unknown) + return -1; + } + + return 0; +} diff --git a/pixel.h b/pixel.h new file mode 100755 index 0000000..933acf7 --- /dev/null +++ b/pixel.h @@ -0,0 +1,55 @@ +// ------------------------------------------------------------ +// pixel.h +// +// Reading writing tiff images +// +// 12.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef PIXEL_H +#define PIXEL_H + +#include "mtypes.h" +#include "imageio.h" + +// ------------------------------------------------------------ +#define MAJOR_STATE_MASK 0x03 + +enum PixelState +{ + pixel_state_all = 0, + pixel_state_foreground = 1, + pixel_state_background = 2, + pixel_state_unknown = 4, + pixel_state_computed = 8 +}; + +typedef struct _pixel_t +{ + UINT32 state; + + pxl_id_t id_this, id_closest_f, id_closest_b; + pos_t raster_coord; + double dist_closest_f, dist_closest_b; + +} pixel_t; + +typedef struct _set_data_t +{ + UINT32 num_row, num_col, num_pixel, num_unknown, num_fore, num_back; + pixel_t *pPixel; + +} set_data_t; + +// ------------------------------------------------------------ +void SetInit(set_data_t *pObj, int num_row, int num_col); +void PixelStateInit(set_data_t *pObj, image_t *pTrimap_image); +int GetCropSet(set_data_t *pObj, int *pCrop, pos_t center, pos_t rect_size, UINT32 select_state); +int GetCircleSet(set_data_t *pObj, int *pCrop, pos_t center, double radius, UINT32 select_state); +int GetSetByState(set_data_t *pObj, pxl_id_t *pSet, UINT32 select_state); +int ComputeDistanceMap(set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown); +int ComputeCalcSet(set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown); +int SaveCalcSet(char *filename, set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown); +int LoadCalcSet(char *filename, set_data_t *pObj, pxl_id_t *pSet_unknown, int num_unknown); + +// ------------------------------------------------------------ +#endif // PIXEL_H diff --git a/tictoc.c b/tictoc.c new file mode 100755 index 0000000..5eda45a --- /dev/null +++ b/tictoc.c @@ -0,0 +1,29 @@ +// ------------------------------------------------------------ +// tictoc.c +// +// Utility to measure timing +// +// 10.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#include +#include +#include +#include +#include +#include "tictoc.h" + +// ------------------------------------------------------------ +void Tic(struct _tictoc_t *pObj) +{ + gettimeofday(&pObj->s, &pObj->tz); +} + +// ------------------------------------------------------------ +void Toc(struct _tictoc_t *pObj) +{ + gettimeofday(&pObj->e, &pObj->tz); + printf("Operation takes %f milliseconds.\n", (double)((pObj->e.tv_sec - pObj->s.tv_sec)*1000) + (double)(pObj->e.tv_usec - pObj->s.tv_usec)/1000); + +} + +// ------------------------------------------------------------ diff --git a/tictoc.h b/tictoc.h new file mode 100755 index 0000000..d8a382d --- /dev/null +++ b/tictoc.h @@ -0,0 +1,24 @@ +// ------------------------------------------------------------ +// tictoc.h +// +// Utility to measure timing +// +// 10.03.2005, J.Ahrensfeld +// ------------------------------------------------------------ +#ifndef TICTOC_H +#define TICTOC_H + +// ------------------------------------------------------------ +typedef struct _tictoc_t +{ + struct timeval s, e; + struct timezone tz; + +} tictoc_t; + +// ------------------------------------------------------------ +void Tic(struct _tictoc_t *pObj); +void Toc(struct _tictoc_t *pObj); + +// ------------------------------------------------------------ +#endif // TICTOC_H