根据硬件计算的架构，用C++建立了一个软件模型，证明硬件计算的可行性。

Shear之后的图像，C++结果：

 

http://s1/middle/679f93564bc33ebdb6710&690

 

MatLab Shear结果:

 

http://s1/middle/679f93564bc33ebdb9110&690

 

Warp结果，C++：

 

http://s10/middle/679f93564bc33ebea3429&690

 

MatLab Warp结果：

 

http://s14/middle/679f93564bc33ebf1533d&690

 

//shear-warp sw model

//zf2266#mail.ustc.edu.cn, Zhanxiang Zhao

 

#include <iostream>

#include <fstream>

#include "windows.h"

#include "math.h"

 

using std::endl;

using std::cout;

using std::cin;

 

#define LINE_LENGTH 81

#define LINE_SHEAR_LENGTH int(LINE_LENGTH*1.7321+2)

#define ALPHA_TABLE_SIZE 1000

#define DIM 4

#define IMAGE_SIZE 400

 

double volume[LINE_LENGTH][LINE_LENGTH][LINE_LENGTH];

double lineBuf[5][LINE_LENGTH];

double shearedBuf[LINE_SHEAR_LENGTH];

double interImgBuf[LINE_SHEAR_LENGTH*LINE_SHEAR_LENGTH];

double renderImage[IMAGE_SIZE*IMAGE_SIZE];

double alphaTable[ALPHA_TABLE_SIZE];

double Mview[DIM*DIM];

double Mshear[DIM*DIM];

double Mwarp[(DIM-1)*(DIM-1)];

double MshearInv[DIM*DIM];

double MwarpInv[(DIM-1)*(DIM-1)];

 

 

struct LUPMat

{

    double* L, *U, *P; //dim*dim

    int* pai; //dim

    int dim; //dimension of LUP, size is d*d

};

 

void DelLUPMat( LUPMat lm)

{

    delete [] lm.L;

    delete [] lm.U;

    delete [] lm.P;

    delete [] lm.pai;

}

 

void printMatrix(double* mat, int dim)

{

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

            printf("%2.5f ", mat[i*dim+j]);

        printf("\n");

    }

}

 

void printVector(double* vec, int dim)

{

    for(int i=0; i<dim; i++)

    {

        printf("%2.2f ", vec[i]);

    }

    printf("\n");

}

 

 

//LUP decomposition of dim*dim matrix

LUPMat LUPDecomp(double *A, int dim)

{

    double* mat = new double[dim*dim];

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

            mat[i*dim+j] = A[i*dim+j];

    double* l = new double[dim*dim];

    double* u = new double[dim*dim];

    double* p = new double[dim*dim];

    int* pai = new int[dim]; //pai[i]=k, then p[i][k]=1

    for(int i=0; i<dim; i++)

        pai[i] = i;

    double cmax; //max value of column

    int maxLineNo; //line of cmax

    for(int k=0; k<dim; k++)

    {

        //find max value of the k-th column

        cmax=0;

        for(int i=k; i<dim; i++)

        {

            if(abs(mat[i*dim+k])>abs(cmax))

            {

                cmax = mat[i*dim+k];

                maxLineNo = i;

            }

        }

        if(cmax==0)

        {

            printf("singular matrix!\n");

            exit(1);

        }

        int tmpk = pai[k];

        pai[k] = pai[maxLineNo];

        pai[maxLineNo] = tmpk;

        double tmpLn;

        //exchange line

        for(int li=0; li<dim; li++)

        {

            tmpLn = mat[k*dim+li];

            mat[k*dim+li] = mat[maxLineNo*dim+li];

            mat[maxLineNo*dim+li] = tmpLn;

        }

        //LU decomposition

        for(int i=k+1; i<dim; i++)

        {

            mat[i*dim+k] = mat[i*dim+k]/mat[k*dim+k];

            for(int j=k+1; j<dim; j++)

                mat[i*dim+j] = mat[i*dim+j] - mat[i*dim+k]*mat[k*dim+j];

        }

    }

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

        {

            if(i>j)

            {

                l[i*dim+j] = mat[i*dim+j];

                u[i*dim+j] = 0;

            }

            else if(i==j)

            {

                u[i*dim+j] = mat[i*dim+j];

                l[i*dim+j] = 1;

            }

            else

            {

                u[i*dim+j] = mat[i*dim+j];

                l[i*dim+j] = 0;

            }

        }

    }

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

            p[i*dim+j] = 0;

        p[i*dim+pai[i]] = 1;

    }

    LUPMat ret;

    ret.L = l;

    ret.U = u;

    ret.P = p;

    ret.pai = pai;

    ret.dim = dim;

    delete [] mat;

    return ret;

}

 

//testbench

void LUPDecomp_tb()

{

    double mat[DIM*DIM] = {

        2, 0, 2, 0.6,

        3, 3, 4, -2,

        5, 5, 4, 2,

        -1, -2, 3.4, -1};

        printMatrix(mat, DIM);

        LUPMat lm = LUPDecomp(mat, DIM);

 

        cout<<"P = "<<endl;

        printMatrix(lm.P, DIM);

 

        cout<<"L = "<<endl;

        printMatrix(lm.L, DIM);

 

        cout<<"U = "<<endl;

        printMatrix(lm.U, DIM);

        DelLUPMat(lm);

}

 

 

double *SolveLinearEq(double* A, double* b, int dim)

{

    LUPMat lm = LUPDecomp(A, dim);

    double* x = new double[dim];

    double* y = new double[dim];

    for(int i=0; i<dim; i++)

    {

        y[i] = b[lm.pai[i]];

        for(int j=0; j<i; j++)

            y[i] -= y[j]*lm.L[i*dim+j];

    }

    //cout<<"y:"<<endl;

    //    printVector(y, dim);

    for(int i=dim-1; i>=0; i--)

    {

        for(int j=dim-1; j>i; j--)

            y[i] -= lm.U[i*dim+j]*x[j];

        x[i] = y[i]/lm.U[i*dim+i];

    }    

    //cout<<"y:"<<endl;

    //    printVector(y, dim);

    delete [] y;

    DelLUPMat(lm);

    return x;

}

 

void SolveLinearEq_tb()

{

    const int dim=3;

    double A[dim*dim] =

    {1, 2, 0,

    3, 4, 4,

    5, 6, 3};

    cout<<"A: "<<endl;

    printMatrix(A, dim);

    double b[dim] ={3, 7, 8};

    double* x = SolveLinearEq(A, b, dim);

    cout<<"x: "<<endl;

    printVector(x, dim);

    delete [] x;

}

 

double* InverseMatrix(double *A, int dim)

{

    double *invA = new double[dim*dim];

    double *e = new double[dim];

    double *x;

    for(int i=0; i<dim; i++)

        e[i] = 0;

 

    for(int i=0; i<dim; i++)

    {

        e[i] = 1;

        if(i>0) e[i-1]=0;

        x = SolveLinearEq(A, e, dim);

        //    cout<<"No. "<<i<<" x: ";

        //    printVector(x, dim);

        //    cout<<"e: ";

        //    printVector(e, dim);

        for(int j=0; j<dim; j++)

            invA[j*dim+i] = x[j];

    }

    delete [] x;

    return invA;

}

 

//matrix multiply

double* MatMult(double* A, double* invA, int dim)

{

    double* aij = new double[dim*dim];

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

        {

            aij[i*dim+j]=0;

            for(int k=0; k<dim; k++)

                aij[i*dim+j] += A[i*dim+k]*invA[k*dim+j];

        }

        return aij;

}

 

void InverseMatrix_tb()

{

    const int dim=3;

    double A[dim*dim] =

    {1, 2, 0,

    3, 4, 4,

    5, 6, 3};

    cout<<"A: "<<endl;

    printMatrix(A, dim);

    double* invA = InverseMatrix(A, dim);

    cout<<"inverse of A: "<<endl;

    printMatrix(invA, dim);

    double* aij = MatMult(A, invA, dim);

    cout<<"A*invA:"<<endl;

    printMatrix(aij, dim);

    delete [] invA;

    delete [] aij;

}

 

void copyMatrix2D(double* src, double* des, int dim)

{

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

            des[i*dim+j] = src[i*dim+j];

}

 

//compute shear warp matrix

//return principle viewing axis: c

int MakeShearWarpMat()

{

    //Find the principal viewing axis

    double Vo[DIM-1] ={

        Mview[0*DIM+1]*Mview[1*DIM+2] - Mview[1*DIM+1]*Mview[0*DIM+2],

        Mview[1*DIM+0]*Mview[0*DIM+2] - Mview[0*DIM+0]*Mview[1*DIM+2],

        Mview[0*DIM+0]*Mview[1*DIM+1] - Mview[1*DIM+0]*Mview[0*DIM+1]

    };

 

    double maxv=Vo[0];

    int c=0;

    for(int i=1; i<DIM-1; i++)

    {

        if(Vo[i]>maxv)

        {    

            maxv=Vo[i];

            c=i;

        }

    }

    // Choose the corresponding Permutation matrix P

    double pyzx[DIM*DIM] = {0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1};

    double pzxy[DIM*DIM] = {0,0,1,0, 1,0,0,0, 0,1,0,0, 0,0,0,1};

    double pxyz[DIM*DIM] = {1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1};

    double* p;

    double P[DIM*DIM];

 

    if(c==0)

        p = pyzx;

    else if(c==1)

        p = pzxy;

    else    //2

        p = pxyz;

    copyMatrix2D(p, P, DIM);

 

    //Compute the permuted view matrix from Mview and P    

    double* invP = InverseMatrix(P, DIM);

    //xiangge

    cout<<"inv P:"<<endl;

    printMatrix(invP, DIM);

    double* Mview_p = MatMult(Mview, invP, DIM);

    //xiangge

    cout<<"inv(Mview_p):"<<endl;

    printMatrix(Mview_p, DIM);

    delete [] invP;

    //180 degrees rotate detection

    if(Mview_p[2*DIM+2]<0) c=c+3;

 

    //Compute the shear coeficients from the permuted view matrix

    double Si = (Mview_p[1*DIM+1]* Mview_p[0*DIM+2] - Mview_p[0*DIM+1]* Mview_p[1*DIM+2]) / (Mview_p[0*DIM+0]* Mview_p[1*DIM+1] - Mview_p[1*DIM+0]* Mview_p[0*DIM+1]);

    double Sj = (Mview_p[0*DIM+0]* Mview_p[1*DIM+2] - Mview_p[1*DIM+0]* Mview_p[0*DIM+2]) / (Mview_p[0*DIM+0]* Mview_p[1*DIM+1] - Mview_p[1*DIM+0]* Mview_p[0*DIM+1]);

    double Ti, Tj;

    // Compute the translation between the orgins of standard object coordinates

    // and intermdiate image coordinates

    double kmax = LINE_LENGTH-1;

    if ((Si>=0)&&(Sj>=0)){ Ti = 0; Tj = 0; }

    if ((Si>=0)&&(Sj<0)){ Ti = 0; Tj = -Sj*kmax; }

    if ((Si<0)&&(Sj>=0)){ Ti = -Si*kmax; Tj = 0; }

    if ((Si<0)&&(Sj<0)){ Ti = -Si*kmax; Tj = -Sj*kmax; }

 

    // Compute the shear matrix

    double mshear[DIM*DIM]=

    {

        1, 0, Si, Ti,

        0, 1, Sj, Tj,

        0, 0, 1 ,0 ,

        0, 0, 0 ,1

    };

    copyMatrix2D(mshear, Mshear, DIM);

 

    // Compute the 2Dwarp matrix

    double mwarp[(DIM-1)*(DIM-1)]={

        Mview_p[0*DIM+0], Mview_p[0*DIM+1], (Mview_p[0*DIM+3]-Ti*Mview_p[0*DIM+0]-Tj*Mview_p[0*DIM+1]),

        Mview_p[1*DIM+0], Mview_p[1*DIM+1], (Mview_p[1*DIM+3]-Ti*Mview_p[1*DIM+0]-Tj*Mview_p[1*DIM+1]),

        0, 0, 1

    };

    copyMatrix2D(mwarp, Mwarp, DIM-1);

    return c;

}

 

//bilinear interpolation of 4 points

double Interp(double* intensity, double* perc)

{

    double result=0.0f;

    for(int i=0; i<4; i++)

        result += intensity[i] * perc[i];

    //Attention: should not use: *intensity++ * *perc++,

    //for system can delete array normally without the change of pointer

    return result;

}

 

//

void Shear()

{

    cout<<"Start Shear"<<endl;

    double* invshear = InverseMatrix(Mshear, DIM);

    copyMatrix2D(invshear, MshearInv, DIM);

    double* invwarp = InverseMatrix(Mwarp, DIM-1);

    copyMatrix2D(invwarp, MwarpInv, DIM-1);

    cout<<"MshearInv:"<<endl;

    printMatrix(MshearInv, DIM);

    cout<<"MwarpInv:"<<endl;

    printMatrix(MwarpInv, DIM-1);

    delete [] invshear;

    delete [] invwarp;

    //c=0

    int py[LINE_LENGTH-1];

    int px[LINE_LENGTH-1];

    for(int z=0; z<(LINE_LENGTH-1); z++)

    {

        cout<<"Shear Slice "<<z+1<<endl;

        //Offset calculation

        double xd=(0-(double)LINE_SHEAR_LENGTH/2)+MshearInv[0*DIM+2]*((double)z-(double)LINE_LENGTH/2)+(double)LINE_LENGTH/2;

        double yd=(0-(double)LINE_SHEAR_LENGTH/2)+MshearInv[1*DIM+2]*((double)z-(double)LINE_LENGTH/2)+(double)LINE_LENGTH/2;

 

        double xdfloor=floor(xd);

        double ydfloor=floor(yd);

 

        //calculate the coordinates on which a image slice starts and

        //ends in the temporary shear image (buffer)

        double pystart=-ydfloor;

        if(pystart<0) pystart=0;

        double pyend = LINE_LENGTH-ydfloor;

        if(pyend>LINE_SHEAR_LENGTH) pyend=LINE_SHEAR_LENGTH;

 

        double pxstart=-xdfloor;

        if(pxstart<0) pxstart=0;

        double pxend = LINE_LENGTH-xdfloor;

        if(pxend>LINE_SHEAR_LENGTH) pxend=LINE_SHEAR_LENGTH;

 

        int npy = pyend-pystart-1;

        int npx = pxend-pxstart-1;

        py[0] = pystart;

        for(int i=1; i<npy; i++)

            py[i] = py[i-1]+1;

        px[0] = pxstart;

        for(int i=1; i<npx; i++)

            px[i] = px[i-1]+1;

        //Linear interpolation constants (percentages)

        double xCom=xd-floor(xd);

        double yCom=yd-floor(yd);

        double perc[4]={(1-xCom)*(1-yCom), (1-xCom)*yCom, xCom*(1-yCom), xCom*yCom};

        double intensity_loc;

        double intensityInterp[4];

        double alphaimage;

        double alphaimage_inv;

        double shearedimage;

        int indexAlpha;

 

        //Determine x and y coordinates of pixel(s) which will be come current pixel

        //bilinear

        for(int iy=0; iy<npy-1; iy++)

        {

            for(int ix=0; ix<npx-1; ix++)

            {

                intensityInterp[0] = volume[z][int(px[ix])+int(xdfloor)][int(py[iy])+int(ydfloor)];

                intensityInterp[1] = volume[z][int(px[ix])+int(xdfloor)][int(py[iy])+int(ydfloor)+1];

                intensityInterp[2] = volume[z][int(px[ix])+int(xdfloor)+1][int(py[iy])+int(ydfloor)];

                intensityInterp[3] = volume[z][int(px[ix])+int(xdfloor)+1][int(py[iy])+int(ydfloor)+1];

                intensity_loc = Interp(intensityInterp, perc);

                //calculate current alphaimage

                indexAlpha = floor(intensity_loc*(ALPHA_TABLE_SIZE-1)+0.5)+1;

                alphaimage=alphaTable[indexAlpha];

                alphaimage_inv=(1-alphaimage);

                //Update the current pixel in the shear image buffer

                interImgBuf[int(px[ix])*LINE_SHEAR_LENGTH+int(py[iy])]=alphaimage_inv*interImgBuf[int(px[ix])*LINE_SHEAR_LENGTH+int(py[iy])]+alphaimage*intensity_loc;

            }

        }

        cout<<"sheared intermediate image: "<<endl;

        std::ofstream ofShear("shearImg.dat");

        for(int i=0; i<LINE_SHEAR_LENGTH; i++)

        {

            for(int j=0; j<LINE_SHEAR_LENGTH; j++)

                ofShear<<interImgBuf[i*LINE_SHEAR_LENGTH+j]<<" ";

            ofShear<<endl;

        }

        ofShear.close();

    }

    cout<<"exit shear"<<endl;

}

 

//set x to boundary

void setBoundary(double &x, double floor, double ceil)

{

    if(x<floor) x=floor;

    if(x>ceil-1) x=ceil-1;

}

 

//

void Warp()

{

    //format the output image

    for(int i=0; i<IMAGE_SIZE; i++)

    {

        for(int j=0; j<IMAGE_SIZE; j++)

        {

            renderImage[i*IMAGE_SIZE+j] = 0;

        }

    }

    cout<<"start Warp"<<endl;

    double *Tlocalx = new double[IMAGE_SIZE*IMAGE_SIZE];

    double *Tlocaly = new double[IMAGE_SIZE*IMAGE_SIZE];

    double *xbas = new double[IMAGE_SIZE*IMAGE_SIZE];

    double *ybas = new double[IMAGE_SIZE*IMAGE_SIZE];

    double *tx = new double[IMAGE_SIZE*IMAGE_SIZE];

    double *ty = new double[IMAGE_SIZE*IMAGE_SIZE];

    //center of output image

    double mean_out = IMAGE_SIZE/2;

    //center of input image

    double mean_in = LINE_SHEAR_LENGTH/2;

 

    double percInterp[4];

    double imageInterp[4];    

    bool check_xbas0, check_ybas0, check_xbas1, check_ybas1;

    for(int i=0; i<IMAGE_SIZE; i++)

    {

        for(int j=0; j<IMAGE_SIZE; j++)

        {

            //center of output image coordinates is (0,0)

            double xd = i-mean_out;

            double yd = j-mean_out;

 

            //Calculate the Transformed coordinates

            Tlocalx[i*IMAGE_SIZE+j] = mean_in + MwarpInv[0*(DIM-1)+0]*xd + MwarpInv[0*(DIM-1)+1]*yd + (MwarpInv[0*(DIM-1)+2]+MshearInv[0*DIM+3]) * 1;

            Tlocaly[i*IMAGE_SIZE+j] = mean_in + MwarpInv[1*(DIM-1)+0]*xd + MwarpInv[1*(DIM-1)+1]*yd + (MwarpInv[1*(DIM-1)+2]+MshearInv[1*DIM+3]) * 1;

 

            //Interp coordinates

            xbas[i*IMAGE_SIZE+j] = floor(Tlocalx[i*IMAGE_SIZE+j]);

            ybas[i*IMAGE_SIZE+j] = floor(Tlocaly[i*IMAGE_SIZE+j]);

            tx[i*IMAGE_SIZE+j] = Tlocalx[i*IMAGE_SIZE+j] - xbas[i*IMAGE_SIZE+j];

            ty[i*IMAGE_SIZE+j] = Tlocaly[i*IMAGE_SIZE+j] - ybas[i*IMAGE_SIZE+j];

        }

    }

    for(int i=0; i<IMAGE_SIZE-1; i++)

    {

        for(int j=0; j<IMAGE_SIZE-1; j++)

        {

            percInterp[0] = (1-tx[i*IMAGE_SIZE+j])*(1-ty[i*IMAGE_SIZE+j]);

            percInterp[1] = (1-tx[i*IMAGE_SIZE+j])*ty[i*IMAGE_SIZE+j];

            percInterp[2] = tx[i*IMAGE_SIZE+j]*(1-ty[i*IMAGE_SIZE+j]);

            percInterp[3] = tx[i*IMAGE_SIZE+j]*ty[i*IMAGE_SIZE+j];

            

            //limit indices to boundaries

            check_xbas0 = xbas[i*IMAGE_SIZE+j]<0 || xbas[i*IMAGE_SIZE+j]>(LINE_SHEAR_LENGTH-1);

            check_ybas0 = ybas[i*IMAGE_SIZE+j]<0 || ybas[i*IMAGE_SIZE+j]>(LINE_SHEAR_LENGTH-1);

            check_xbas1 = (xbas[i*IMAGE_SIZE+j]+1)<0 || (xbas[i*IMAGE_SIZE+j]+1)>(LINE_SHEAR_LENGTH-1);

            check_ybas1 = (ybas[i*IMAGE_SIZE+j]+1)<0 || (ybas[i*IMAGE_SIZE+j]+1)>(LINE_SHEAR_LENGTH-1);

 

            //set to boundary

            setBoundary(xbas[i*IMAGE_SIZE+j], 0, LINE_SHEAR_LENGTH);

            setBoundary(xbas[i*IMAGE_SIZE+j], -1, LINE_SHEAR_LENGTH-1); //xbas[i]+1

            setBoundary(ybas[i*IMAGE_SIZE+j], 0, LINE_SHEAR_LENGTH);

            setBoundary(ybas[i*IMAGE_SIZE+j], -1, LINE_SHEAR_LENGTH-1);

 

            imageInterp[0] = interImgBuf[int(xbas[i*IMAGE_SIZE+j])*LINE_SHEAR_LENGTH+int(ybas[i*IMAGE_SIZE+j])];

            imageInterp[1] = interImgBuf[int(xbas[i*IMAGE_SIZE+j])*LINE_SHEAR_LENGTH+int(ybas[i*IMAGE_SIZE+j]+1)];

            imageInterp[2] = interImgBuf[int(xbas[i*IMAGE_SIZE+j]+1)*LINE_SHEAR_LENGTH+int(ybas[i*IMAGE_SIZE+j])];

            imageInterp[3] = interImgBuf[int(xbas[i*IMAGE_SIZE+j]+1)*LINE_SHEAR_LENGTH+int(ybas[i*IMAGE_SIZE+j]+1)];

            

            imageInterp[0] = (check_xbas0 || check_ybas0) ? 0 : imageInterp[0];

            imageInterp[1] = (check_xbas0 || check_ybas1) ? 0 : imageInterp[1];

            imageInterp[2] = (check_xbas1 || check_ybas0) ? 0 : imageInterp[2];

            imageInterp[3] = (check_xbas0 || check_ybas1) ? 0 : imageInterp[3];

 

            renderImage[i*IMAGE_SIZE+j] = Interp(imageInterp, percInterp);

        }

    }

    delete [] Tlocalx ;

    delete [] Tlocaly ;

    delete [] xbas ;

    delete [] ybas ;

    delete [] tx ;

    delete [] ty ;

    cout<<"exit Warp"<<endl;

}

 

//

bool render(const char* volumefile, const char* alphafile)

{

    for(int i=0; i<LINE_SHEAR_LENGTH; i++)

    {

        for(int j=0; j<LINE_SHEAR_LENGTH; j++)

            interImgBuf[i*LINE_SHEAR_LENGTH+j]=0;

    }

    for(int i=0; i<LINE_LENGTH; i++)

    {

        for(int j=0; j<LINE_LENGTH; j++)

            renderImage[i*LINE_LENGTH+j]=0;

    }

    std::ifstream ifVolume(volumefile);

    if(ifVolume.fail())

    {

        cout<<"unable to open volume file!"<<endl;

        return false;

    }

    //read volume

    for(int i=0; i<LINE_LENGTH; i++)

        for(int j=0; j<LINE_LENGTH; j++)

            for(int k=0; k<LINE_LENGTH; k++)

            {

                ifVolume>>volume[i][j][k];

            }

    ifVolume.close();

    std::ifstream ifAlpha(alphafile);

    if(ifAlpha.fail())

    {

        cout<<"unable to open alpha table file!"<<endl;

        return false;

    }

    //read alpha table

    for(int i=0; i<ALPHA_TABLE_SIZE; i++)

        ifAlpha>>alphaTable[i];

    ifAlpha.close();

    int c=MakeShearWarpMat();

    cout<<"Mshear:"<<endl;

    printMatrix(Mshear, DIM);

    cout<<"Mwarp:"<<endl;

    printMatrix(Mwarp, DIM-1);

    Shear();

    Warp();

    std::ofstream ofImage("render.dat");

    for(int i=0; i<IMAGE_SIZE; i++)

    {

        for(int j=0; j<IMAGE_SIZE; j++)

            ofImage<<renderImage[i*IMAGE_SIZE+j]<<" ";

        ofImage<<endl;

    }

    ofImage.close();

}

 

int main()

{

    double mview[DIM*DIM]=

{     0.6783, 1.6190, -0.2468, 0,

0.9310, -0.1614, 1.4998, 0,

1.3473, -0.7035, -0.9121, 0,

0, 0, 0, 1.0000

    };

    copyMatrix2D(mview, Mview, DIM);

    const char* volumefile = "volume.dat";

    const char* alphafile = "alphatable.dat";

    render(volumefile, alphafile);

    system("pause");

}

在xilinx SDK基于microblaze测试pcie示例程序一直不能run，简单的下面代码

printf("x=%8.8x\n", x);
xil_printf("x=%8.8x\n", x);

 

第一个就不能运行，第二个就可以。很奇怪。

 

http://www.eefocus.com/walkie/blog/09-02/165549_bf04c.html

这篇文章讲了一些区别，但没解释我看到的现象。

按照算法导论上矩阵求逆的原理，写了个C++代码。

有A*X=I，则X为inv（A）。

1. 用高斯消元法对原矩阵LUP分解，使得PA=LU。P为置换矩阵，为了使每次消元时取有列最大值的行。L下三角阵，U上三角阵。
2. 根据分解结果求出线性方程组A*xi=ei的xi向量，ei是单位阵I的列向量i，则xi为X的对应列向量。
3. 把xi组合成X，即为逆矩阵。

    

struct LUPMat

{

    float* L, *U, *P; //dim*dim

    int* pai; //dim

    int dim; //dimension of LUP, size is d*d

};

 

void DelLUPMat( LUPMat lm)

{

    delete [] lm.L;

    delete [] lm.U;

    delete [] lm.P;

    delete [] lm.pai;

}

 

void printMatrix(float* mat, int dim)

{

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

            printf("%2.2f ", mat[i*dim+j]);

        printf("\n");

    }

}

 

void printVector(float* vec, int dim)

{

    for(int i=0; i<dim; i++)

    {

        printf("%2.2f ", vec[i]);

    }

        printf("\n");

}

 

 

//LUP decomposition of dim*dim matrix

LUPMat LUPDecomp(float *A, int dim)

{

    float* mat = new float[dim*dim];

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

            mat[i*dim+j] = A[i*dim+j];

    float* l = new float[dim*dim];

    float* u = new float[dim*dim];

    float* p = new float[dim*dim];

    int* pai = new int[dim]; //pai[i]=k, then p[i][k]=1

    for(int i=0; i<dim; i++)

        pai[i] = i;

    float cmax; //max value of column

    int maxLineNo; //line of cmax

    for(int k=0; k<dim; k++)

    {

        //find max value of the k-th column

        cmax=0;

        for(int i=k; i<dim; i++)

        {

            if(abs(mat[i*dim+k])>abs(cmax))

            {

                cmax = mat[i*dim+k];

                maxLineNo = i;

            }

        }

        if(cmax==0)

        {

            printf("singular matrix!\n");

            exit(1);

        }

        int tmpk = pai[k];

        pai[k] = pai[maxLineNo];

        pai[maxLineNo] = tmpk;

        float tmpLn;

        //exchange line

        for(int li=0; li<dim; li++)

        {

            tmpLn = mat[k*dim+li];

            mat[k*dim+li] = mat[maxLineNo*dim+li];

            mat[maxLineNo*dim+li] = tmpLn;

        }

        //LU decomposition

        for(int i=k+1; i<dim; i++)

        {

            mat[i*dim+k] = mat[i*dim+k]/mat[k*dim+k];

            for(int j=k+1; j<dim; j++)

                mat[i*dim+j] = mat[i*dim+j] - mat[i*dim+k]*mat[k*dim+j];

        }

    }

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

        {

            if(i>j)

            {

                l[i*dim+j] = mat[i*dim+j];

                u[i*dim+j] = 0;

            }

            else if(i==j)

            {

                u[i*dim+j] = mat[i*dim+j];

                l[i*dim+j] = 1;

            }

            else

            {

                u[i*dim+j] = mat[i*dim+j];

                l[i*dim+j] = 0;

            }

        }

    }

    for(int i=0; i<dim; i++)

    {

        for(int j=0; j<dim; j++)

            p[i*dim+j] = 0;

        p[i*dim+pai[i]] = 1;

    }

    LUPMat ret;

    ret.L = l;

    ret.U = u;

    ret.P = p;

    ret.pai = pai;

    ret.dim = dim;

    delete [] mat;

    return ret;

}

 

//testbench

void LUPDecomp_tb()

{

    float mat[DIM*DIM] = {

        2, 0, 2, 0.6,

        3, 3, 4, -2,

        5, 5, 4, 2,

        -1, -2, 3.4, -1};

    printMatrix(mat, DIM);

    LUPMat lm = LUPDecomp(mat, DIM);

 

    cout<<"P = "<<endl;

    printMatrix(lm.P, DIM);

 

    cout<<"L = "<<endl;

    printMatrix(lm.L, DIM);

 

    cout<<"U = "<<endl;

    printMatrix(lm.U, DIM);

    DelLUPMat(lm);

}

 

 

float *SolveLinearEq(float* A, float* b, int dim)

{

    LUPMat lm = LUPDecomp(A, dim);

    float* x = new float[dim];

    float* y = new float[dim];

    for(int i=0; i<dim; i++)

    {

        y[i] = b[lm.pai[i]];

        for(int j=0; j<i; j++)

            y[i] -= y[j]*lm.L[i*dim+j];

    }

    //cout<<"y:"<<endl;

    printVector(y, dim);

    for(int i=dim-1; i>=0; i--)

    {

        for(int j=dim-1; j>i; j--)

            y[i] -= lm.U[i*dim+j]*x[j];

        x[i] = y[i]/lm.U[i*dim+i];

    }    

    //cout<<"y:"<<endl;

    printVector(y, dim);

    delete [] y;

    DelLUPMat(lm);

    return x;

}

 

void SolveLinearEq_tb()

{

    const int dim=3;

    float A[dim*dim] =

    {1, 2, 0,

     3, 4, 4,

     5, 6, 3};

    cout<<"A: "<<endl;

    printMatrix(A, dim);

    float b[dim] ={3, 7, 8};

    float* x = SolveLinearEq(A, b, dim);

    cout<<"x: "<<endl;

    printVector(x, dim);

    delete [] x;

}

 

float* InverseMatrix(float *A, int dim)

{

    float *invA = new float[dim*dim];

    float *e = new float[dim];

    float *x;

    for(int i=0; i<dim; i++)

        e[i] = 0;

 

    for(int i=0; i<dim; i++)

    {

        e[i] = 1;

        if(i>0) e[i-1]=0;

        x = SolveLinearEq(A, e, dim);

    //    cout<<"No. "<<i<<" x: ";

        printVector(x, dim);

    //    cout<<"e: ";

        printVector(e, dim);

        for(int j=0; j<dim; j++)

            invA[j*dim+i] = x[j];

    }

    delete [] x;

    return invA;

}

 

float* isInverse(float* A, float* invA, int dim)

{

    float* aij = new float[dim*dim];

    for(int i=0; i<dim; i++)

        for(int j=0; j<dim; j++)

        {

            aij[i*dim+j]=0;

            for(int k=0; k<dim; k++)

                aij[i*dim+j] += A[i*dim+k]*invA[k*dim+j];

        }

    return aij;

}

 

void InverseMatrix_tb()

{

    const int dim=3;

    float A[dim*dim] =

    {1, 2, 0,

     3, 4, 4,

     5, 6, 3};

    cout<<"A: "<<endl;

    printMatrix(A, dim);

    float* invA = InverseMatrix(A, dim);

    cout<<"inverse of A: "<<endl;

    printMatrix(invA, dim);

    float* aij = isInverse(A, invA, dim);

    cout<<"A*invA:"<<endl;

    printMatrix(aij, dim);

    delete [] invA;

    delete [] aij;

}

 

测试程序的结果：

 

http://s7/middle/679f93564bb3523c25716&amp;690

作为新型的二维超声换能器阵列，CMUT阵列相比原来的二维压电陶瓷换能器阵列有什么优势呢，我们来看一篇论文：引文1，斯坦福大学Omer Oralkan等人通过CMUT换能器与压电陶瓷换能器的实验比较，测试了CMUT的性能。里面比较了CMUT与二维压电陶瓷换能器。

 

二维压电陶瓷换能器缺点：

1. 要制造很多阵元困难，连线复杂；
2. 栅瓣。为了减少栅瓣，阵元间距必须小于声波波长，（声速1540m/s，频率5MHz时，波长308 um；10MHz，154um；频率高，探测精度高）这种小尺寸对与压电陶瓷材料限制了发射能量并降低了接收灵敏度，很难达到较大的动态范围。

 

而CMUT优势：

1. 更好带宽；
2. 灵敏度高；
3. 易于制造；
4. 尺寸小；
5. 自身噪声低；
6. 还可以与CMOS IC集成封装。

论文中用直径36um CMUT和400us*400um压电换能器做实验。下图为压电陶瓷PZT与CMUT灵敏度比较，相比CMUT 200% Cp（寄生电容）模型，压电陶瓷只在中心频率3MHz附近有优势。通过把接收放大器集成到CMUT芯片上，可以降低寄生电容。

http://s16/middle/679f93564baf73bbd0e0f&amp;690

 

二者信噪比比较，压电陶瓷只在中心频率3MHz附近有优势。

http://s13/middle/679f93564baf73bd2878c&amp;690

 

Reference

[1] O. Oralkan, X. Jin, F. L.Degertekin, and B. T.Khuri-Yakub, "Simulation and experimental characterization of a 2-D capacitive micromachined ultrasonic transducer array element," IEEE Trans. Ultrason. Ferroelect. Freq. Contr., vol. 46, pp. 1337–1340, Nov. 1999.

昨天开组会大家讨论问题，有几个疑问被提出来了：

1. Field II主要建模算法用C代码实现。

    

2. 二维压电陶瓷换能器阵列和CMUT二维换能器阵列有什么区别？

    

   本文只讨论第一个问题，第二个另文细述。

    

Field II代码

 

和大二小师弟讨论Field II建模CMUT换能器时，发现Field II的Matlab代码只是来做参数配置和调用函数的，这些函数编译为可执行文件，主要的建模算法是用C代码实现，而这个代码并没有开源。我给Dr. Jensen发了封邮件要C代码，但愿能回复我。

 

建模方法

 

使用空间冲击响应（spatial impulse response）的原理，依托于超声场的线性系统理论。发射时，当换能器被狄拉克delta函数激发后，空间中每个点的声场是关于时间的函数。声场可以通过空间冲激响应和激发函数卷积得到。接收时，把换能器激发函数和发射孔径的空间冲击响应、接受孔径的空间冲击响应卷积得到接收响应，通过换能器的机械电子传递函数得到接收电压曲线。

 

软件使用

 

如果只是用Field II仿真声场，可以直接调用Matlab函数，配置参数，比如探头的物理参数（阵元类型，数量，尺寸，探测物，接收孔径等），然后得到接收信号，可以拿来波束合成。我们用一段示例程序来产生血流的接收波束信号如图。

http://s13/middle/679f93564baf697eb57fc&amp;690

 

% This examples shows how the procedures can be used for making flow data from a number of scatters in a tube.

% Example of use of the new Field II program running under Matlab

%

% This example shows how flow can simulated

%

% This script assumes that the field_init procedure has been called

%

% Example by Joergen Arendt Jensen, March 22, 2011.

% Generate the transducer apertures for send and receive

f0=3e6; % Transducer center frequency [Hz]

fs=100e6; % Sampling frequency [Hz]

c=1540; % Speed of sound [m/s]，人体内声音平均速度

lambda=c/f0; % Wavelength

element_height=5/1000; % Height of element [m]

kerf=0.1/1000; % Kerf [m]

focus=[0 0 70]/1000; % Fixed focal point [m]

% Generate aperture

aperture = xdc_linear_array (128, lambda/2, element_height, kerf, 1, 1,focus);

% Set the impulse response and excitation of the emit aperture

impulse_response=sin(2*pi*f0*(0:1/fs:2/f0));

impulse_response=impulse_response.*hanning(max(size(impulse_response)))';

xdc_impulse (aperture, impulse_response);

excitation=sin(2*pi*f0*(0:1/fs:8/f0));

xdc_excitation (aperture, excitation);

% Set the seed of the random number generator

randn('seed',sum(100*clock))

% Initialize the ranges for the scatterers，血管里面的物质

% Notice that the coordinates are in meters

x_range=0.015; % x range for the scatterers [m]

y_range=0.015; % y range for the scatterers [m]

z_range=0.015; % z range for the scatterers [m]

z_offset=0.70; % Offset of the mid-point of the scatterers [m]

R=0.005; % Radius of blood vessel [m]

% Set the number of scatterers. It should be roughly

% 10 scatterers per resolution cell

c=1540; % Ultrasound propagation velocity [m/s]

f0=3e6; % Center frequency of transducer [Hz]

lambda=c/f0;

N=round(10*x_range/(5*lambda)*y_range/(5*lambda)*z_range/(lambda*2));

disp([num2str(N),' Scatterers'])

% Generate the coordinates and amplitude

% Coordinates are rectangular within the range.

% The amplitude has a Gaussian distribution.

x=x_range*(rand(1,N)-0.5);

y=y_range*(rand(1,N)-0.5);

z=z_range*(rand(1,N)-0.5);

% Find which scatterers that lie within the blood vessel

r=(y.^2+z.^2).^0.5;

within_vessel= (r < R)';

% Assign an amplitude and a velocity for each scatterer

v0=0.5; % Largest velocity of scatterers [m/s]

velocity=v0*(1-(r/R).^2).*within_vessel';

blood_to_stationary= 0.1; % Ratio between amplitude of blood to stationary tissue

amp=randn(N,1).*((1-within_vessel) + within_vessel*blood_to_stationary);

% Calculate a suitable Tprf

theta=45/180*pi;

f_max=2*v0*cos(theta)/c*f0;

fprf=3*f_max;

Tprf=1/fprf; % Time between pulse emissions [s]

Nshoots=128; % Number of shoots

% Find the response by calling field

for i=1:Nshoots

i

% Generate the rotated and offset block of sample

xnew=x*cos(theta)+z*sin(theta);

znew=z*cos(theta)-x*sin(theta) + z_offset;

scatterers=[xnew; y; znew;]' ;

% Calculate the received response

[v, t1]=calc_scat(aperture, aperture, scatterers, amp);

% Store the result

image_data(1:max(size(v)),i)=v;

times(i) = t1;

% Propagate the scatterers and alias them

% to lie within the correct range

x=x + velocity*Tprf;

outside_range= (x > x_range/2);

x=x - x_range*outside_range;

end

% Here the display of the data

% Here the display of the data is inserted

plot(image_data)

Volume Rendering

常见的三维成像有表面绘制和体绘制，surface rendering, volume rendering,后者需要庞大的计算量，为了利用硬件加速计算的的优势，我只关注volume rendering。

 

主要包括以下方法：

( 1 ) 基于图像空间扫描的绘制， 如光线投射( Ray Casting ) 算法；该方法是对射入观察者眼睛的所有光线，沿着光线路径对所有体元按照颜色和不透明度积分得到最后射出的颜色。Ray tracing：逆向跟踪射入的光线，遇到交点时类似积分。

( 2 ) 基于物体空间扫描的绘制， 如足迹表法( 又叫是splatting算法) 、 错切变形( Shear-Warp ) 算法、 体元投射( Cell Projection ) 法等。

 

Shear-Warp原理

其中Shear-Warp算法把三维投影变为三维数据场错切shear和二维图像变形warp两步实现，是目前最快的Volume Rendering算法。

我们来看看斯坦福大学Philippe Lacroute和 Marc Levoy的论文Fast Volume Rendering Using a Shear-Warp Factorization of the Viewing Transformation，Shear-Warp算法分为三个步骤：

1. 物体空间->错切图像空间。如图一，横线为一个个slice，就是三维物体的二维切片，从上看就是线。Viewing rays是平行投影光线。经过shear变换后，在图像空间，切片和光线垂直。对于透视的情况，如图二，切片还要收缩。

http://s9/middle/679f93564bab577d54e98&amp;690

http://s3/middle/679f93564bab577dc4602&amp;690

1. 按照从前往后的顺序组合切片。得到错切图像空间的二维中间图像。
2. 通过Warp，把中间图像变形为最终投影图像。

    

Shear-Warp详解

Shear-Warp相当于对原坐标系做了变换，变换矩阵为：

http://s11/middle/679f93564bab577e41dea&amp;690

P说明以那个方向为主轴来切片，有xyz三种可能。S从物体空间变到错切空间。M(warp)再变回图像空间。以z为主轴时，平行投影S有

http://s15/middle/679f93562f4b16b0aeb3e&amp;690

透视投影S有，

http://s1/middle/679f93564bab5781e5860&amp;690

其中切片z0被拉伸了http://s5/middle/679f93564bab5781ddd04&amp;690。矩阵元素从M(view)得到。

最后，

http://s16/middle/679f93564bab57821840f&amp;690

图3为前述三个步骤详解。

1. 沿着扫描线变换并重采样，甚至拉伸。
2. 把切片的对应扫描线投影并叠加，alpha处理透明度。
3. 根据M(warp)变换，需要再次重采样。

http://s2/middle/679f93564bab5782e6001&amp;690

Shear-warp算法

对于平行投影，首先对于原数据游程编码RLE（run-length encode），

1. 1. 根据透明度阈值在体素（物体的像素别称）数据结构中作记录，跳过透明的体素，在rendering之前进行，对于三种矩阵P形式，及三种切片方式都要计算。2. 对于中间图像，在像素数据结构中保存此不透明像素到可透光像素位移，在rendering同时进行。最后只处理不完全透明和可视的体素。

http://s1/middle/679f93564bab57830b870&amp;690

1. 开始重采样和组合。采用双线性插值滤波器（bilinear interp），每两个相邻扫描线计算得到一个中间图像扫描线。背光用聚合卷积，正面光用散射卷积。用查找表的方法实现透明度处理。

   http://s5/middle/679f93560791225a31624&amp;690

   http://s1/middle/679f93564bab5786da080&amp;690

2. 再次用双线性滤波器把大的中间图像变为最终小图像。（general-purpose affine image warper with a bilinear filter）

这篇文章还提出了一种alpha传递函数，通过递归的办法实现透明度可变时的行程编码。先得到一个立方体八个顶点确定的八叉树，查表计算立方体是否透明，不透明则继续递归计算更小立方体，直到体积小到一定值。

http://s12/middle/679f93564bab578794fdb&amp;690