commit 05da8712f62879b20741a9475d4892870b2d24d7 Author: Jens Ahrensfeld Date: Sat Jul 19 07:44:42 2014 +0000 Initial import git-svn-id: http://moon:8086/svn/software/trunk/libsrc/_to_be_added@1 b431acfa-c32f-4a4a-93f1-934dc6c82436 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