From ebb2c2f0e0e6871f0c03bb859327cdee2d7dedd6 Mon Sep 17 00:00:00 2001 From: Takuo Yasunaga Date: Sun, 8 Mar 2015 00:56:13 +0900 Subject: [PATCH] ticket: #34564 and #34965 new file: util/X86MAC64/cuda x86Mac64.inc ticket: #34964 mrcImageShapeSearch Instable version --- src/Config/x86Mac64.inc | 2 +- .../mrcImageShapeSearch/src/mrcImageShapeSearch.c | 225 ++++++++++++++------- .../src/mrcImageShapeSearch.html | 2 +- .../mrcImage/mrcImageShapeSearch/src/test/Makefile | 2 +- 4 files changed, 157 insertions(+), 74 deletions(-) diff --git a/src/Config/x86Mac64.inc b/src/Config/x86Mac64.inc index d21be5bb0a..51218d10b7 100755 --- a/src/Config/x86Mac64.inc +++ b/src/Config/x86Mac64.inc @@ -274,7 +274,7 @@ INCCUDA = -I$(EOS_HOME)/util/$(OSTYPE)/cuda/inc -I$(EOS_HOME)/util/$(OSTYPE)/cud #LIBCUDA = -L$(EOS_HOME)/util/$(OSTYPE)/cuda/NVIDIA_CUDA_SDK/common/lib -L$(EOS_HOME)/util/$(OSTYPE)/cuda/NVIDIA_CUDA_SDK/lib -L$(EOS_HOME)/util/$(OSTYPE)/cuda/lib -lcudart -lcutil -lcufft #LIBCUDA = -L$(EOS_HOME)/util/$(OSTYPE)/cuda/lib64 -L$(EOS_HOME)/util/$(OSTYPE)/cuda/lib64 -lcudart -lcufft -LIBCUDA = -L$(EOS_HOME)/util/$(OSTYPE)/cuda/lib64 -L$(EOS_HOME)/util/$(OSTYPE)/cuda/lib -lcudart -lcufft +LIBCUDA = -L$(EOS_HOME)/util/$(OSTYPE)/cuda/lib -lcudart -lcufft #HOSTDEPENDENTLIB += $(LIBCUDA) diff --git a/src/Tools/mrcImage/mrcImageShapeSearch/src/mrcImageShapeSearch.c b/src/Tools/mrcImage/mrcImageShapeSearch/src/mrcImageShapeSearch.c index dcf748f674..13712b8f23 100755 --- a/src/Tools/mrcImage/mrcImageShapeSearch/src/mrcImageShapeSearch.c +++ b/src/Tools/mrcImage/mrcImageShapeSearch/src/mrcImageShapeSearch.c @@ -23,7 +23,7 @@ #include "mrcImage.h" typedef struct lmrcImageShapeSearchInfo { - float radius; // Sylinder, plane + float radius; // Sylinder, half disk float minRadius; float maxRadius; float delRadius; @@ -33,11 +33,6 @@ typedef struct lmrcImageShapeSearchInfo { float maxLength; float delLength; - float thickness; // Plane - float minThickness; - float maxThickness; - float delThickness; - float minTheta; float maxTheta; float delTheta; @@ -50,6 +45,11 @@ typedef struct lmrcImageShapeSearchInfo { float maxPsi; float delPsi; + // temporary + float x; + float y; + float z; + int interpMode; int thresZscore; @@ -66,12 +66,13 @@ typedef struct lmrcImageShapeSearchInfo { typedef enum lmrcImageShapeSearchMode { lmrcImageShapeSearchModeSylinder=0, - lmrcImageShapeSearchModePlane=1 + lmrcImageShapeSearchModeHalfDisk=1 } lmrcImageShapeSearchMode; extern void lmrcImageShapeSearch(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode); -extern void lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode); -extern void lmrcImageShapeSearchPlane(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode); +extern void lmrcImageShapeSearchCalc0(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode); +extern void lmrcImageShapeSearchSylinder(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode); +extern void lmrcImageShapeSearchHalfDisk(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode); extern void lmrcImageShapeSearchModePrint(FILE* fpt); int @@ -145,7 +146,7 @@ void lmrcImageShapeSearchModePrint(FILE* fpt) { fprintf(fpt, "%d: Cylinder with radius and length\n", lmrcImageShapeSearchModeSylinder); - fprintf(fpt, "%d: Plane with radius and thickness\n", lmrcImageShapeSearchModePlane); + fprintf(fpt, "%d: Half disk with radius and length(thickness)\n", lmrcImageShapeSearchModeHalfDisk); } @@ -154,11 +155,11 @@ lmrcImageShapeSearch(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linf { switch(mode) { case lmrcImageShapeSearchModeSylinder: { - lmrcImageShapeSearchSylinder(out, in, linfo, mode); + lmrcImageShapeSearchCalc0(out, in, linfo, mode); break; } - case lmrcImageShapeSearchModePlane: { - lmrcImageShapeSearchPlane(out, in, linfo, mode); + case lmrcImageShapeSearchModeHalfDisk: { + lmrcImageShapeSearchCalc0(out, in, linfo, mode); break; } default: { @@ -170,7 +171,7 @@ lmrcImageShapeSearch(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linf void -lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode) +lmrcImageShapeSearchCalc0(mrcImage* out, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode) { float x, y, z; float srcx, srcy, srcz; @@ -178,7 +179,7 @@ lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchIn float theta, phi, psi, radius, length; int iPhi, iTheta, iPsi, iRadius, iLength; int nPhi, nTheta, nPsi, nRadius, nLength; - float max, ave, sum, sd, zscore; + float max, ave, sum, sd, zscore, min; double* data; double d, cp, sp; double score; @@ -187,7 +188,7 @@ lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchIn Matrix3D mat; floatVector v; int i, imax; - int k, k0; + int k, k0, nVoxel; nPhi = (int)((linfo->maxPhi - linfo->minPhi) /linfo->delPhi) +1; nTheta = (int)((linfo->maxTheta - linfo->minTheta)/linfo->delTheta)+1; @@ -220,9 +221,13 @@ lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchIn for(z=0; zHeaderN.z; z++) { DEBUGPRINT1("z: %f\n", z); + linfo->z = z; for(y=0; yHeaderN.y; y++) { //DEBUGPRINT1("y: %f\n", y); + linfo->y = y; for(x=0; xHeaderN.x; x++) { // Each voxel + linfo->x = x; + k0 = (int)(x + y*in->HeaderN.x + z*in->HeaderN.x*in->HeaderN.y); for(iPhi=0; iPhiminPhi + iPhi*linfo->delPhi; @@ -246,54 +251,28 @@ lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchIn +iRadius*nPhi*nTheta*nPsi +iLength*nPhi*nTheta*nPsi*nRadius; - score = 0; - for(p=0; p<2*M_PI; p+=linfo->delPsi) { - cp = cos(p); - sp = sin(p); - for(l=0; l<=length; l+=linfo->delLength) { - v.data[2] = l-length/2; - for(r=0; r<=radius; r+=linfo->delRadius) { - v.data[0] = r*cp; - v.data[1] = r*sp; - - //matrix3DMultiplyVector(&v, mat); - srcx = mat[0][0]*v.data[0] + mat[1][0]*v.data[1] + mat[2][0]*v.data[2]; - srcy = mat[0][1]*v.data[0] + mat[1][1]*v.data[1] + mat[2][1]*v.data[2]; - srcz = mat[0][2]*v.data[0] + mat[1][2]*v.data[1] + mat[2][2]*v.data[2]; - - dstx = x + srcx; - dsty = y + srcy; - dstz = z + srcz; - //DEBUGPRINT3("dst: %f %f %f\n", dstx, dsty, dstz); - //mrcPixelDataGet(in, dstx, dsty, dstz, &d, mrcPixelRePart, linfo->interpMode); - if(-0.5HeaderN.x-0.5 - &&-0.5HeaderN.y-0.5 - &&-0.5HeaderN.z-0.5) { - k = (int)(dstx+0.5) + (int)(dsty+0.5)*in->HeaderN.x + (int)(dstz+0.5)*in->HeaderN.x*in->HeaderN.y; - d = in->ImageFloatImage[k]; - } else { - d = 0; - } - score += d; - } - - // Edge Check - v.data[0] = r*cos(p); - v.data[1] = r*sin(p); - matrix3DMultiplyVector(&v, mat); - dstx = x + v.data[0]; - dsty = y + v.data[1]; - dstz = z + v.data[2]; - mrcPixelDataGet(in, dstx, dsty, dstz, &d, mrcPixelRePart, linfo->interpMode); - score -= d; - } - } - data[index] = score; + linfo->radius = radius; + linfo->length = length; + switch(mode) { + case lmrcImageShapeSearchModeSylinder: { + lmrcImageShapeSearchSylinder(&(data[index]), in, mat, linfo, mode); #ifdef DEBUG2 - if(0Max, x, y, z, max, mrcPixelRePart); linfo->Max.ImageFloatImage[k0] = max; - zscore = ((sd>0)?(max-ave)/sd:0); - //mrcPixelDataSet(&linfo->Zscore, x, y, z, zscore, mrcPixelRePart); - linfo->Zscore.ImageFloatImage[k0] = zscore; - - d = ((zscore>linfo->thresZscore)?1:0); - //mrcPixelDataSet(out, x, y, z, d, mrcPixelRePart); - out->ImageFloatImage[k0] = d; - i = imax%nPhi; //mrcPixelDataSet(&linfo->shapeInfo[0], x, y, z, (i*linfo->delPhi + linfo->minPhi)*DEGREE, mrcPixelRePart); linfo->shapeInfo[0].ImageFloatImage[k0] = (i*linfo->delPhi + linfo->minPhi)*DEGREE; @@ -366,10 +337,122 @@ lmrcImageShapeSearchSylinder(mrcImage* out, mrcImage* in, lmrcImageShapeSearchIn } } } + + nVoxel = in->HeaderN.x*in->HeaderN.y*in->HeaderN.z; + max = -1e30; + min = 1e30; + ave = 0; + sd = 0; + for(k=0; kMax.ImageFloatImage[k]; + ave += d; + if(sdSD.ImageFloatImage[k]) { + sd = linfo->SD.ImageFloatImage[k]; + } + if(maxMax.ImageFloatImage[k]; + zscore = (d-linfo->average.ImageFloatImage[k])/linfo->SD.ImageFloatImage[k]; + //zscore = (d-ave)/sd; + linfo->Zscore.ImageFloatImage[k] = zscore*in->ImageFloatImage[k]; + + d = (d-min)*sum; + out->ImageFloatImage[k] = d*in->ImageFloatImage[k]; + + //DEBUGPRINT3("%d z: %f out: %f\n", k, zscore, d); + } +} + +void +lmrcImageShapeSearchSylinder(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode) +{ + double p, l, r; + double x, y, z; + double srcx, srcy, srcz; + double dstx, dsty, dstz; + double scorePos, scoreNeg; + double cp, sp, d; + int k, k0; + int countPos, countNeg; + + //DEBUGPRINT("lmrcImageShapeSearchSylinder start\n"); + scorePos = 0; + scoreNeg = 0; + countPos = 0; + countNeg = 0; + for(p=0; p<2*M_PI; p+=linfo->delPsi) { + cp = cos(p); + sp = sin(p); + for(l=0; l<=linfo->length; l+=linfo->delLength) { + z = l-linfo->length/2; + for(r=0; r<=linfo->radius; r+=linfo->delRadius) { + x = r*cp; + y = r*sp; + + srcx = mat[0][0]*x + mat[1][0]*y + mat[2][0]*z; + srcy = mat[0][1]*x + mat[1][1]*y + mat[2][1]*z; + srcz = mat[0][2]*x + mat[1][2]*y + mat[2][2]*z; + + dstx = linfo->x + srcx; + dsty = linfo->y + srcy; + dstz = linfo->z + srcz; + if(-0.5HeaderN.x-0.5 + &&-0.5HeaderN.y-0.5 + &&-0.5HeaderN.z-0.5) { + k = (int)(dstx+0.5) + (int)(dsty+0.5)*in->HeaderN.x + (int)(dstz+0.5)*in->HeaderN.x*in->HeaderN.y; + d = in->ImageFloatImage[k]; + } else { + d = 0; + } + + countPos++; + scorePos += d; + } + // Edge Check + x = r*cos(p); + y = r*sin(p); + srcx = mat[0][0]*x + mat[1][0]*y + mat[2][0]*z; + srcy = mat[0][1]*x + mat[1][1]*y + mat[2][1]*z; + srcz = mat[0][2]*x + mat[1][2]*y + mat[2][2]*z; + + dstx = linfo->x + srcx; + dsty = linfo->y + srcy; + dstz = linfo->z + srcz; + if(-0.5HeaderN.x-0.5 + &&-0.5HeaderN.y-0.5 + &&-0.5HeaderN.z-0.5) { + k = (int)(dstx+0.5) + (int)(dsty+0.5)*in->HeaderN.x + (int)(dstz+0.5)*in->HeaderN.x*in->HeaderN.y; + d = in->ImageFloatImage[k]; + } else { + d = 0; + } + scoreNeg += d; + countNeg++; + } + } + k0 = (int)(linfo->x+0.5) + (int)(linfo->y+0.5)*in->HeaderN.x + (int)(linfo->z+0.5)*in->HeaderN.x*in->HeaderN.y; + *data = (0.8*scorePos/countPos - 0.2*scoreNeg/countNeg)*in->ImageFloatImage[k0]; +#ifdef DEBUG2 + //DEBUGPRINT1("data: %f\n", score); + if(0 diff --git a/src/Tools/mrcImage/mrcImageShapeSearch/src/test/Makefile b/src/Tools/mrcImage/mrcImageShapeSearch/src/test/Makefile index 021860d06b..2658be269c 100755 --- a/src/Tools/mrcImage/mrcImageShapeSearch/src/test/Makefile +++ b/src/Tools/mrcImage/mrcImageShapeSearch/src/test/Makefile @@ -11,7 +11,7 @@ help: exec: @echo "----- Execution Check -----" - ../$(OSTYPE)/$(OBJECTNAME) -i data/test-sylinder.bin -o data/test-sylinder.syl -r 3 5 1 -l 5 5 1 -Phi 0 359 30 -Theta - 0 189 30 -Psi 0 0 90 -m 0 -thres 2.5 + time ../$(OSTYPE)/$(OBJECTNAME) -i data/test-sylinder.bin -o data/test-sylinder.syl -r 3 5 1 -l 5 5 1 -Phi 0 359 45 -Theta 0 179 45 -Psi 0 0 90 -m 0 -thres 0.5 @echo "----- Calc check -----" init: -- 2.11.0