|
Отправлено: 15.02.13 22:44. Заголовок: The code is shown be..
The code is shown below: It's working perfect, What I'm trying to do is code what I attached in the excel sheet into this one below: // Upscaling-FMM.cpp : Defines the entry point for the console application. // 28jan-2012.cpp : Defines the entry point for the console application. // 27jan2012.cpp : Defines the entry point for the console application. // // 11-dec-noon.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include <iostream> #include <vector> #include <math.h> #include <time.h> #include <ctime> #include <fstream> #include <iomanip> //#include <vector> //#include <algorithm> #include <stdio.h> #include <stdlib.h> using namespace std; //========================================================================================================= double *dxc; double *dyc; double *dzc; double *dx; double *dy; double *dz; //========================================================================================================= //#define eps 2.2204460492503131e-16 #define eps 3.e-16 #define doublemax 1e50 #define INF 2e50 //#define INF 111 #define listINF 2.345e50 //#define listINF 111 #ifndef min #define min(a,b) ((a) < (b) ? (a): (b)) #endif #ifndef max #define max(a,b) ((a) > (b) ? (a): (b)) #endif /* Find minimum value of an array and return its index */ int minarray(double *A, int l) { int i; int minind=0; for (i=0; i<l; i++) { if(A<A[minind]){ minind=i; } } return minind; } double pow2(double val) { return val*val; } int iszero(double a){ return a*a<eps; } int isnotzero(double a){ return a*a>eps; } void roots(double* Coeff, double* ans) { double a=Coeff[0]; double b=Coeff[1]; double c=Coeff[2]; double r1, r2; double d; d=max(pow2(b)-4.0*a*c,0); if(isnotzero(a)) { ans[0]= (-b - sqrt(d)) / (2.0*a); ans[1]= (-b + sqrt(d)) / (2.0*a); } else { r1=(-b - sqrt(d)); r2=(-b + sqrt(d)); if(isnotzero(r1)) { if(isnotzero(r2)) { ans[0]= (2.0*c)/r1; ans[1]= (2.0*c)/r2; } else { ans[0]= (2.0*c)/r1; ans[1]= (2.0*c)/r1; } } else if(isnotzero(r2)) { ans[0]= (2.0*c)/r2; ans[1]= (2.0*c)/r2; } else { ans[0]=0; ans[1]=0; } } } int maxarray(double *A, int l) { int i; int maxind=0; for (i=0; i<l; i++) { if(A>A[maxind]){ maxind=i; } } return maxind; } __inline int mindex3(int x, int y, int z, int sizx, int sizy) { return x+y*sizx+z*sizx*sizy; } __inline bool IsFinite(double x) { return (x <= doublemax && x >= -doublemax ); } __inline bool IsInf(double x) { return (x >= doublemax ); } __inline bool IsListInf(double x){ return (x == listINF ); } __inline bool isntfrozen3d(int i, int j, int k, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(k>=0)&&(i<dims[0])&&(j<dims[1])&&(k<dims[2])&&(Frozen[mindex3(i, j, k, dims[0], dims[1])]==0); } __inline bool isfrozen3d(int i, int j, int k, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(k>=0)&&(i<dims[0])&&(j<dims[1])&&(k<dims[2])&&(Frozen[mindex3(i, j, k, dims[0], dims[1])]==1); } int p2x(int x) /* 2^x */ { /* return pow(2,x); */ int y=1; int p2x[16] ={1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768}; while(x>15) { x=x-15; y=y*32768; } return y*p2x[x]; } void show_list(double **listval, int *listprop) { int z, k; double mm; for(z=0;z<listprop[1]; z++) { for(k=0;k<p2x(z+1); k++) { mm=k; if((z>0)&&(listval[z-1][(int)floor(mm/2)]>listval[z][k])) { printf("*%6.2f", listval[z][k]); } else { printf(" %6.2f", listval[z][k]); } } printf("\n"); } } void initialize_list(double ** listval, int *listprop) { /* Loop variables */ int i; /* Current Length, Orde and Current Max Length */ listprop[0]=0; listprop[1]=1; listprop[2]=2; /* Make first orde storage of 2 values */ listval[0]=new double[2]; /* Initialize on infinite */ for(i=0;i<2;i++) { listval[0]=listINF; } } void destroy_list(double ** listval, int *listprop) { /* Loop variables */ int i, list_orde; /* Get list orde */ list_orde=listprop[1]; /* Free memory */ for(i=0;i<list_orde;i++) { free(listval); } free(listval); free(listprop); } void list_add(double ** listval, int *listprop, double val) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i, j; /* Temporary list location */ int listp; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* If list is full expand list */ if(list_length==list_lengthmax) { list_lengthmax=list_lengthmax*2; for (i=list_orde; i>0; i--) { listval=listval[i-1]; listval = (double *)realloc(listval, p2x(i+1)*sizeof(double)); for(j=p2x(i); j<p2x(i+1); j++) { listval[j]=listINF; } } listval[0]=new double[2]; listval[0][0]=min(listval[1][0], listval[1][1]); listval[0][1]=listINF; list_orde++; } /* Add a value to the list */ listp=list_length; list_length++; listval[list_orde-1][listp]=val; /* Update the links minimum */ for (i=list_orde-1; i>0; i--) { listp=(int)floor(((double)listp)/2); if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; } } /* Set list parameters */ listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax; } int list_minimum(double ** listval, int *listprop) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i; /* Temporary list location */ int listp; /* Index of Minimum */ int minindex; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Follow the minimum through the binary tree */ listp=0; for(i=0;i<(list_orde-1);i++) { /* <= ????? */ if(listval[listp]<=listval[listp+1]) { listp=listp*2; } else { listp=(listp+1)*2; } } i=list_orde-1; if(listval[listp]<=listval[listp+1]){minindex=listp; } else { minindex=listp+1; } return minindex; } void list_remove(double ** listval, int *listprop, int index) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i; /* Temp index */ int index2; double val; double valmin; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Temporary store current value */ val=listval[list_orde-1][index]; valmin=listINF; /* Replace value by infinite */ listval[list_orde-1][index]=listINF; /* Follow the binary tree to replace value by minimum values from */ /** the other values in the binary tree. */ i=list_orde-1; while(true) { if((index%2)==0) { index2=index+1; } else { index2=index-1; } if(val<listval[index2]) { index=(int)floor(((double)index2)/2.0); if(listval[index2]<valmin) { valmin=listval[index2]; } if(i==0) { break; } listval[i-1][index]=valmin; i--; if(i==0) { break; } } else { break; } } } void list_remove_replace(double ** listval, int *listprop, int index) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* Loop variable */ int i, listp; /* Temporary store value */ double val; int templ; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Remove the value the minimum number and goup in tree by comparision its adjacent number to the parents */ list_remove(listval, listprop, index); /* Replace the removed value by the last in the list. (if it was */ /* not already the last value) */ if(index<(list_length-1)) { /* Temporary store last value in the list */ val=listval[list_orde-1][list_length-1]; /* Remove last value in the list */ list_remove(listval, listprop, list_length-1); /* Add a value to the list */ listp=index; listval[list_orde-1][index]=val; /* Update the links minimum */ for (i=(list_orde-1); i>0; i--) { listp=(int)floor(((double)listp)/2); if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; } } } /* List is now shorter */ list_length--; /* Remove trunk of binary tree / Free memory if list becomes shorter */ if(list_orde>2&&IsListInf(listval[0][1])) { //cout<<"truncate-start"<<endl; //show_list(listval, listprop); list_orde--; list_lengthmax=list_lengthmax/2; /* Remove trunk array */ free(listval[0]); /* Move the other arrays one up */ templ=2; for (i=0; i<list_orde; i++) { listval=listval[i+1]; /* Resize arrays to their new shorter size */ listval = (double *)realloc(listval, templ*sizeof(double)); templ*=2; } } /* Set list parameters */ listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax; //show_list(listval, listprop); //cout<<"truncate-end"<<endl; } void listupdate(double **listval, int *listprop, int index, double val) { /* List parameters */ int list_length, list_orde, list_lengthmax; /* loop variable */ int i, listp; /* Get list parameters */ list_length=listprop[0]; list_orde=listprop[1]; list_lengthmax=listprop[2]; /* Add a value to the list */ listp=index; listval[list_orde-1][index]=val; /* Update the links minimum */ for (i=(list_orde-1); i>0; i--) { listp=(int)floor(((double)listp)/2); if(val<listval[i-1][listp]) { listval[i-1][listp]=val; } else { break; } } /* Set list parameters */ listprop[0]=list_length; listprop[1]=list_orde; listprop[2]=list_lengthmax; } __inline int mindex2(int x, int y, int sizx) { return x+y*sizx; } __inline bool isntfrozen2d(int i, int j, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(i<dims[0])&&(j<dims[1])&&(Frozen[i+j*dims[0]]==0); } __inline bool isfrozen2d(int i, int j, int *dims, bool *Frozen) { return (i>=0)&&(j>=0)&&(i<dims[0])&&(j<dims[1])&&(Frozen[i+j*dims[0]]==1); } //========================================================================================================= double second_derivative(double Txm1, double Txm2, double Txp1, double Txp2) { bool ch1, ch2; double Tm; Tm=INF; ch1=(Txm2<Txm1)&&IsFinite(Txm1); ch2=(Txp2<Txp1)&&IsFinite(Txp1); if(ch1&&ch2) { Tm =min( (4.0*Txm1-Txm2)/3.0 , (4.0*Txp1-Txp2)/3.0);} else if(ch1) { Tm =(4.0*Txm1-Txm2)/3.0; } else if(ch2) { Tm =(4.0*Txp1-Txp2)/3.0; } return Tm; } double CalculateDistance(double *T, double Fijk, int *dims, int i, int j, int k, bool *Frozen) { /* Loop variables */ int q, t; /* Current location */ int in, jn, kn; /* Derivatives */ double Tm[3]={0, 0, 0}; double Coeff[3]; /* local derivatives in distance image */ double Txm1, Txm2, Txp1, Txp2; double Tym1, Tym2, Typ1, Typ2; double Tzm1, Tzm2, Tzp1, Tzp2; double Tt, Tt2,xm1,xp1,ym1,yp1,zm1,zp1,x1,x2,y1,y2,z1,z2; /* Return values root of polynomial */ double ansroot[2]={0, 0}; /* Order derivatives in a certain direction */ int Order[3]={0, 0, 0}; /* Neighbours 4x2 */ int ne[18]={-1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1}; /* Stencil constants */ double G1[3]={1, 1, 1}; x1 = dx[mindex3(i, j, k, dims[0], dims[1])]; y1 = dy[mindex3(i, j, k, dims[0], dims[1])]; z1 = dz[mindex3(i, j, k, dims[0], dims[1])]; /*Get First order derivatives (only use frozen pixel) */ in=i-1; jn=j+0; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Txm1=T[mindex3(in, jn, kn, dims[0], dims[1])];xm1=dx[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Txm1=INF; xm1=0; } in=i+1; jn=j+0; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Txp1=T[mindex3(in, jn, kn, dims[0], dims[1])];xp1=dx[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Txp1=INF; xp1=0; } in=i+0; jn=j-1; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tym1=T[mindex3(in, jn, kn, dims[0], dims[1])];ym1=dy[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tym1=INF; ym1=0;} in=i+0; jn=j+1; kn=k+0; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Typ1=T[mindex3(in, jn, kn, dims[0], dims[1])];yp1=dy[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Typ1=INF; yp1=0;} in=i+0; jn=j+0; kn=k-1; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tzm1=T[mindex3(in, jn, kn, dims[0], dims[1])];zm1=dz[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tzm1=INF; zm1=0;} in=i+0; jn=j+0; kn=k+1; if(isfrozen3d(in, jn, kn, dims, Frozen)) { Tzp1=T[mindex3(in, jn, kn, dims[0], dims[1])];zp1=dz[mindex3(in, jn, kn, dims[0], dims[1])]; } else { Tzp1=INF; zp1=0;} /*Make 1e order derivatives in x and y direction */ Tm[0] = min( Txm1 , Txp1); if(IsFinite(Tm[0])){ Order[0]=1; } else { Order[0]=0; } Tm[1] = min( Tym1 , Typ1); if(IsFinite(Tm[1])){ Order[1]=1; } else { Order[1]=0; } Tm[2] = min( Tzm1 , Tzp1); if(IsFinite(Tm[2])){ Order[2]=1; } else { Order[2]=0; } // X2=XM1 IF CONDITIONAL EXPRESSION IS TRUE AND XP1 IF IT IS FALSE x2 = Tm[0] == Txm1 ? xm1 : xp1; y2 = Tm[1] == Tym1 ? ym1 : yp1; z2 = Tm[2] == Tzm1 ? zm1 : zp1; G1[0] = abs(x2 + x1)/2; G1[1] = abs(y2 + y1)/2; G1[2] = abs(z2 + z1)/2; /*Calculate the distance using x and y direction */ Coeff[0]=0; Coeff[1]=0; Coeff[2]=-1/(max(pow2(Fijk),eps)); for (t=0; t<3; t++) { switch(Order[t]) { case 1: Coeff[0]+=G1[t]; Coeff[1]+=-2.0*Tm[t]*G1[t]; Coeff[2]+=pow2(Tm[t])*G1[t]; break; } } roots(Coeff, ansroot); Tt=max(ansroot[0], ansroot[1]); return Tt; } //starting your modificatiob //============================================================================== int main() { clock_t startClock,finishClock,currentClock; vector<int>::iterator pos; double timeCount,temp_time; startClock = clock(); //-----your sort function goes here here-------- // Declare size of two dimensional array and initialize. int nx,ny,nz,nt,i,j,k,initial_zero,index,inf; cout << endl << "enter nx == " ; cin>> nx; cout << endl << "enter ny == " ; cin>> ny; cout << endl << "enter nz == " ; cin>> nz; nt=nx*ny*nz; startClock = clock(); //ofstream timefile; cout<<"t1= "; cout<<startClock/1000<<endl; int s,x,y,z,XYZ_index,w,IJK_index,itt,npixels,q; double Tt; double *F, *SourcePoints; int ne[18]={-1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0, 0, -1, 0, 0, 1}; /* Neighbour list */ int neg_free; int neg_pos; double *neg_listv; double *neg_listx; double *neg_listy; double *neg_listz; double *neg_listo; /* Matrix containing the Frozen Pixels" */ bool *Frozen; /* Size of input image */ int *dims_c; int dims[3]; /* Size of SourcePoints array */ int num_source; cout << endl << "enter num_source == " ; cin>> num_source; // assign the dimensions dims[0]=nx; dims[1]=ny; dims[2]=nz; npixels=dims[0]*dims[1]*dims[2]; //double* dx = NULL; dx = new double[nt]; //double* dy = NULL; dy = new double[nt]; //double* dz = NULL; dz = new double[nt]; for (int ii=0;ii<npixels;ii++){ dx[ii]=1; dy[ii]=1; dz[ii]=1; } /* The output distance image */ double *T; int *listprop; double **listval; double *Final_T; double *perm; double *permx; double *permx_orig; double *time_orig; double *dtc; double *porosity; double *porosity_orig; double *source; permx_orig=new double[npixels]; time_orig=new double[npixels]; dtc=new double[npixels]; porosity_orig=new double[npixels]; Final_T=new double[npixels]; perm=new double[npixels]; permx=new double[npixels]; porosity=new double[npixels]; source=new double[num_source*num_source*3]; int kk; kk=-1; for(int i=0;i<num_source;i++){ for(int j=0;j<num_source;j++){ k=0; kk=kk+1; source[3*kk]=nx/2+1; source[3*kk+1]=ny/2+1; source[3*kk+2]=nz/2+1; cout<<source[3*(kk)]<<"-"<<source[3*kk+1]<<"-"<<source[3*kk+2]<<endl; } } for (int ii=0;ii<npixels;ii++){ permx[ii]=1; permx_orig[ii]=1; porosity[ii]=1; porosity_orig[ii]=1; } int ii; T= Final_T; F=permx; SourcePoints=source; /* Pixels which are processed and have a final distance are frozen */ Frozen = new bool[npixels]; for(q=0;q<npixels;q++){Frozen[q]=0; T[q]=-1;} /*Free memory to store neighbours of the (segmented) region */ neg_free = 1000000; neg_pos=0; neg_listx = new double [neg_free]; neg_listy = new double [neg_free]; neg_listz = new double [neg_free]; /* List parameters array */ listprop=new int [3]; /* Make jagged list to store a maximum of 2^64 values */ //listval= (double **)malloc( 64* sizeof(double *) );/// num of rows=64? listval= new double* [64]; /* Initialize parameter list */ initialize_list(listval, listprop); neg_listv=listval[listprop[1]-1]; for (s=0; s<num_source*num_source; s++) { /*starting point */ x= (int)SourcePoints[0+s*3]-1; y= (int)SourcePoints[1+s*3]-1; z= (int)SourcePoints[2+s*3]-1; XYZ_index=mindex3(x, y, z, dims[0], dims[1]); Frozen[XYZ_index]=1; T[XYZ_index]=0; } for (s=0; s<num_source*num_source; s++) { /*starting point */ x= (int)SourcePoints[0+s*3]-1; y= (int)SourcePoints[1+s*3]-1; z= (int)SourcePoints[2+s*3]-1; XYZ_index=mindex3(x, y, z, dims[0], dims[1]); for (w=0; w<6; w++) { /*Location of neighbour */ i=x+ne[w]; j=y+ne[w+6]; k=z+ne[w+12]; IJK_index=mindex3(i, j, k, dims[0], dims[1]); /*Check if current neighbour is not yet frozen and inside the */ /*picture */ if(isntfrozen3d(i, j, k, dims, Frozen)) { Tt=(1/(max(F[IJK_index],eps))); /*Update distance in neigbour list or add to neigbour list */ if(T[IJK_index]>0) { if(neg_listv[(int)T[IJK_index]]>Tt) { listupdate(listval, listprop, (int)T[IJK_index], Tt); } } else { /*If running out of memory at a new block */ //if(neg_pos>=neg_free) { // neg_free+=100000; // neg_listx = (double *)realloc(neg_listx, neg_free*sizeof(double) ); // neg_listy = (double *)realloc(neg_listy, neg_free*sizeof(double) ); // neg_listz = (double *)realloc(neg_listz, neg_free*sizeof(double) ); //} list_add(listval, listprop, Tt); neg_listv=listval[listprop[1]-1]; neg_listx[neg_pos]=i; neg_listy[neg_pos]=j; neg_listz[neg_pos]=k; T[IJK_index]=neg_pos; neg_pos++; } } } } /*Loop through all pixels of the image */ for (itt=0; itt<(npixels); itt++) /* */ { /*Get the pixel from narrow list (boundary list) with smallest */ /*distance value and set it to current pixel location */ index=list_minimum(listval, listprop); neg_listv=listval[listprop[1]-1]; /* Stop if pixel distance is infinite (all pixels are processed) */ if(IsInf(neg_listv[index])) { break; } /*index=minarray(neg_listv, neg_pos); */ x=(int)neg_listx[index]; y=(int)neg_listy[index]; z=(int)neg_listz[index]; XYZ_index=mindex3(x, y, z, dims[0], dims[1]); Frozen[XYZ_index]=1; T[XYZ_index]=neg_listv[index]; //cout<< T[XYZ_index]<<endl; //show_list(listval, listprop); /*Remove min value by replacing it with the last value in the array */ //1- remove the min value and go up in tree by comparing it with its adjacent bothers and parent //2- remove the last value in the list and again go up in tree //3- replace the minimum by the last value and again go up in the tree //- if we have extra tree delete them list_remove_replace(listval, listprop, index) ; //show_list(listval, listprop); neg_listv=listval[listprop[1]-1]; if(index<(neg_pos-1)) { neg_listx[index]=neg_listx[neg_pos-1]; neg_listy[index]=neg_listy[neg_pos-1]; neg_listz[index]=neg_listz[neg_pos-1]; T[(int)mindex3((int)neg_listx[index], (int)neg_listy[index], (int)neg_listz[index], dims[0], dims[1])]=index; } neg_pos =neg_pos-1; /*Loop through all 6 neighbours of current pixel */ for (w=0;w<6;w++) { /*Location of neighbour */ i=x+ne[w]; j=y+ne[w+6]; k=z+ne[w+12]; IJK_index=mindex3(i, j, k, dims[0], dims[1]); /*Check if current neighbour is not yet frozen and inside the */ /*picture */ if(isntfrozen3d(i, j, k, dims, Frozen)) { //Tt=CalculateDistance(T, F[IJK_index], dims, i, j, k, usesecond, usecross, Frozen); Tt=CalculateDistance(T, F[IJK_index], dims, i, j, k, Frozen); /*Update distance in neigbour list or add to neigbour list */ IJK_index=mindex3(i, j, k, dims[0], dims[1]); if((T[IJK_index]>-1)&&T[IJK_index]<=listprop[0]) { if(neg_listv[(int)T[IJK_index]]>Tt) { listupdate(listval, listprop, (int)T[IJK_index], Tt); } } else { list_add(listval, listprop, Tt); //show_list(listval, listprop); neg_listv=listval[listprop[1]-1]; neg_listx[neg_pos]=i; neg_listy[neg_pos]=j; neg_listz[neg_pos]=k; T[IJK_index]=neg_pos; neg_pos++; } } } } /* Free memory */ /* Destroy parameter list */ destroy_list(listval, listprop); free(neg_listx); free(neg_listy); free(neg_listz); free(Frozen); finishClock=clock(); timeCount = finishClock - startClock ; cout<<"run time = "; cout<<timeCount/1000/60<<endl; //================================================================================== ofstream myfile1; myfile1.open ("fine-time.txt"); //myfile << "Writing this to a file.\n"; double temp=0; myfile1 << setw(16) ; for (int k=1; k<nz+1;k++){ for (int j=1; j<ny+1;j++){ for (int i=1; i<nx+1;i++){ index= (k-1) * nx * ny +(j-1) * (nx) + (i-1) ; /*if(Final_T[index]==inf){ Final_T[index]=1.; }*/ myfile1 <<Final_T[index]<< setw(16) ; // save fine scale time time_orig[index]= Final_T[index]; //myfile1 <<F[index]<< setw(16) ; temp=temp+Final_T[index]; } myfile1 << endl ; } } myfile1.close(); cout<<"Summation-time = "; cout<<temp/1000; //===================================================================================== }
|