"-i","-i[nput]","Input","Essential","1","1","In","inFile::mrcImage","NULL"
"-o","-o[utput]","Output","Essential","1","1","Out","outFile::mrcImage","NULL"
"-Shape","-Shape","Output","Optional","1","1","Shape","outFile::mrcImage","NULL"
+"-n","-n[Cube]","Mode","Optional","1","1","nCube","Integer","0"
"-r","-r[adius]","Radius[pixel]","Optional","3","1","minR","Integer","3","2","maxR","Integer","5","3","delR","Integer","1"
"-l","-l[ength]","Length[pixel]","Optional","3","1","minL","Integer","3","2","maxL","Integer","5","3","delL","Integer","1"
"-Phi","-Phi","around z-axis[degree]","Optional","3","1","minPhi","Real","0.0","2","maxPhi","Real","359.0","3","delPhi","Real","1.0"
float maxPsi;
float delPsi;
+ int nCube;
+
// temporary
float x;
float y;
int interpMode;
int thresZscore;
- mrcImage shape; // Template Structure
-
// Output
mrcImage average; // Average for all orientations
mrcImage SD; // SD for all
mrcImage Max; // Max for all
mrcImage Zscore; // Z-score;
mrcImage PCA; // PCA
+ mrcImage PCAaccuracy;
int nShapeInfo;
mrcImage* shapeInfo; // Shape, Orientation, ...
+ mrcImage shape; // Template Structure
+
} lmrcImageShapeSearchInfo;
typedef enum lmrcImageShapeSearchMode {
linfo.thresZscore = info.thresZscore;
linfo.interpMode = info.interpMode;
+ //
+ linfo.nCube = info.nCube;
+
mrcFileRead(&in, info.In, "in main", 0);
lmrcImageShapeSearch(&out, &in, &linfo, info.mode);
mrcFileWrite(&linfo.Zscore, filename, "in main", 0);
sprintf(filename, "%s.max", info.Out);
mrcFileWrite(&linfo.Max, filename, "in main", 0);
+
sprintf(filename, "%s.PCA", info.Out);
mrcFileWrite(&linfo.PCA, filename, "in main", 0);
+ sprintf(filename, "%s.PCAaccuracy", info.Out);
+ mrcFileWrite(&linfo.PCAaccuracy, filename, "in main", 0);
+
for(i=0; i<linfo.nShapeInfo; i++) {
sprintf(filename, "%s.shapeinfo-%02d", info.Out, i);
mrcFileWrite(&linfo.shapeInfo[i], filename, "in main", 0);
}
- //mrcFileWrite(&linfo.shape, info.Shape, "in main", 0);
+ if(info.flagShape) {
+ mrcFileWrite(&linfo.shape, info.Shape, "in main", 0);
+ }
exit(EXIT_SUCCESS);
}
Array X;
Array Ave;
mrcImageParaTypeInteger nCube;
- double W;
+ double W, Wsylinder, Wdisk, Wsphere, Waccuracy;
nPhi = (int)((linfo->maxPhi - linfo->minPhi) /linfo->delPhi) +1;
nTheta = (int)((linfo->maxTheta - linfo->minTheta)/linfo->delTheta)+1;
mrcInit(&linfo->Zscore, NULL);
linfo->Max.Header = in->Header;
mrcInit(&linfo->Max, NULL);
+
linfo->PCA.Header = in->Header;
mrcInit(&linfo->PCA, NULL);
+ linfo->shape.Header = in->Header;
+ mrcInit(&linfo->shape, NULL);
+ linfo->PCAaccuracy.Header = in->Header;
+ mrcInit(&linfo->PCAaccuracy, NULL);
+
linfo->nShapeInfo = 5;
linfo->shapeInfo = (mrcImage*)memoryAllocate(sizeof(mrcImage)*linfo->nShapeInfo, "in MFPCalc");
for(i=0; i<linfo->nShapeInfo; i++) {
X.dim = 2;
X.n[0] = 3;
- nCube = ((int)(MAX(linfo->maxLength, 2*linfo->maxRadius+1)+3)/2)*2;
+ if(0<linfo->nCube) {
+ nCube = linfo->nCube;
+ } else {
+ nCube = ((int)(MAX(linfo->maxLength, 2*linfo->maxRadius+1)+3)/2)*2;
+ }
X.n[1] = nCube*nCube*nCube;
arrayInit(&X, "in");
DEBUGPRINT2("nCube: %d %lu\n", nCube, X.n[1]);
linfo->x = x;
k0 = (int)(x + y*in->HeaderN.x + z*in->HeaderN.x*in->HeaderN.y);
+ linfo->shape.ImageFloatImage[k0] = -1;
//mrcPixelDataGet(in, x, y, z, &d, mrcPixelRePart, mrcPixelHowNearest);
k = x + y*in->HeaderN.x + z*in->HeaderN.x*in->HeaderN.y;
d = in->ImageFloatImage[k];
- if(0<d) {
+ if(0<d) {
numPoint = 0;
- for(srcz=-nCube/2; srcz<nCube/2; srcz++) {
- for(srcy=-nCube/2; srcy<nCube/2; srcy++) {
- for(srcx=-nCube/2; srcx<nCube/2; srcx++) {
- mrcPixelDataGet(in, x+srcx, y+srcy, z+srcz, &d, mrcPixelRePart, mrcPixelHowNearest);
- //srcx0 = x + srcx;
- //srcy0 = y + srcy;
- //srcz0 = z + srcz;
- //k = srcx0 + srcy0*in->HeaderN.x + srcz0*in->HeaderN.x*in->HeaderN.y;
+ for(srcz=-nCube/2; srcz<=nCube/2; srcz++) {
+ for(srcy=-nCube/2; srcy<=nCube/2; srcy++) {
+ for(srcx=-nCube/2; srcx<=nCube/2; srcx++) {
+ //mrcPixelDataGet(in, x+srcx, y+srcy, z+srcz, &d, mrcPixelRePart, mrcPixelHowNearest);
+ srcx0 = ((int)(x + srcx)+in->HeaderN.x)%in->HeaderN.x;
+ srcy0 = ((int)(y + srcy)+in->HeaderN.y)%in->HeaderN.y;
+ srcz0 = ((int)(z + srcz)+in->HeaderN.z)%in->HeaderN.z;
+ k = srcx0 + srcy0*in->HeaderN.x + srcz0*in->HeaderN.x*in->HeaderN.y;
d = in->ImageFloatImage[k];
if(0<d) {
arrayDataSet2(X, 0, numPoint, srcx0);
}
}
DEBUGPRINT4("%f %f %f numPoint: %d\n", x, y, z, numPoint);
- if(16<numPoint) {
+ if(5<numPoint) {
X.n[1] = numPoint;
arrayPCA(&U, &C, &Lambda, &X, &Ave, PCAMode);
if(PCAMode==0) PCAMode=1;
L1 = arrayDataGet1(Lambda, 1);
L2 = arrayDataGet1(Lambda, 2);
- DEBUGPRINT4("W: %f %f %f %f\n", W, L0, L1, L2);
-
+ DEBUGPRINT3("Lambda: %f %f %f\n", L0, L1, L2);
L = L0+L1+L2;
L0 /= L;
L1 /= L;
L2 /= L;
+ DEBUGPRINT3("Lambda: %f %f %f\n", L0, L1, L2);
+
+ Wsphere = 100*exp(-((L1/L0) + (L2/L1) + 0.5*(L0 - L2)/L0));
+ Wdisk = 100*exp(-((L2/L0) + (L2/L1) + 0.5*(L0 - L1)/L0));
+ Wsylinder = 100*exp(-((L1/L0) + (L2/L0) + 0.5*(L1 - L2)/L1));
+ DEBUGPRINT3("W: sphere %f disk %f sylinder %f\n", Wsphere, Wdisk, Wsylinder);
+
+ if(Wsphere<Wdisk) {
+ if(Wdisk<Wsylinder) { // sphere -> disk -> sylinder
+ linfo->shape.ImageFloatImage[k0] = lmrcImageShapeSearchModeSylinder;
+ Waccuracy = Wsylinder - Wdisk;
+ } else { // sphere/sylinder -> disk
+ linfo->shape.ImageFloatImage[k0] = lmrcImageShapeSearchModeDisk;
+ if(Wsphere<Wsylinder) { // sphere -> sylinder -> disk
+ Waccuracy = Wdisk - Wsylinder;
+ } else { // sphere -> sylinder -> disk
+ Waccuracy = Wdisk - Wsphere;
+ }
+ }
+ } else {
+ if(Wsphere<Wsylinder) { // Wdisk -> Wsphere -> Wsylinder
+ linfo->shape.ImageFloatImage[k0] = lmrcImageShapeSearchModeSylinder;
+ Waccuracy = Wsylinder - Wsphere;
+ } else { // Wdisk/Wsylinder -> Wsphere
+ linfo->shape.ImageFloatImage[k0] = lmrcImageShapeSearchModeSphere;
+ if(Wdisk<Wsylinder) { // Wdisk -> Wsylinder -> Wsphere
+ Waccuracy = Wsphere - Wsylinder;
+ } else { // Wsylinder -> Wdisk -> Wsphere
+ Waccuracy = Wsphere - Wdisk;
+ }
+ }
+ }
switch(mode) {
case lmrcImageShapeSearchModeSylinder: {
- //W = ((L1!=0)?(L0/L1):1)*(L2!=0?(L0/L2):1);
- W = 1.0/(SQR(1-L0)+SQR(L1)+SQR(L2));
+ W = Wsylinder;
DEBUGPRINT4("W: %f %f %f %f\n", W, L0, L1, L2);
-#ifdef DEBUG
- if(1<W) {
+#ifdef DEBUG2
+ if(0.2<W) {
DEBUGPRINT3("%f %f %f\n", x, y, z);
DEBUGPRINT4("W: %f %f %f %f\n", W, L0, L1, L2);
}
break;
}
case lmrcImageShapeSearchModeDisk: {
- W = ((L2!=0)?((L0/L2)*(L1/L2)):1);
-#ifdef DEBUG
- if(1000<W) {
+ W = Wdisk;
+#ifdef DEBUG2
+ if(0.2<W) {
DEBUGPRINT3("%f %f %f\n", x, y, z);
DEBUGPRINT4("W: %f %f %f %f\n", W, L0, L1, L2);
}
break;
}
case lmrcImageShapeSearchModeSphere: {
- W = L0*L1+L1*L2+L2*L0/SQR(L0+L1+L2);
+ W = Wsphere;
break;
}
default: {
}
} else {
W = 1;
- }
+ } // numPoint
+
} else {
- W = 1;
- }
+ W = 0;
+ } // 0 < d
//mrcPixelDataSet(&linfo->PCA, x, y, z, W, mrcPixelRePart);
- linfo->PCA.ImageFloatImage[k] = W;
- if(1<W) {
+ linfo->PCA.ImageFloatImage[k0] = W;
+ if(linfo->shape.ImageFloatImage[k0]==mode) {
+ linfo->PCAaccuracy.ImageFloatImage[k0] = W*Waccuracy/100/100;
+ } else {
+ linfo->PCAaccuracy.ImageFloatImage[k0] = 0;
+ }
+
+ if(1<W) {
for(iPhi=0; iPhi<nPhi; iPhi++) { // Three Rotation
phi = linfo->minPhi + iPhi*linfo->delPhi;
for(iTheta=0; iTheta<nTheta; iTheta++) {
DEBUGPRINT3("%f %f %f \n", max, ave, sd);
}
#endif
- } else {
+ } else {
ave = 0;
sd = 0;
max = 0;
imax = 0;
- }
+ }
//mrcPixelDataSet(&linfo->average, x, y, z, ave, mrcPixelRePart);
linfo->average.ImageFloatImage[k0] = ave;
sum = 100./(max - min);
for(k=0; k<nVoxel; k++) {
d = linfo->Max.ImageFloatImage[k];
- zscore = (d-linfo->average.ImageFloatImage[k])/linfo->SD.ImageFloatImage[k];
+
+ zscore = (d-linfo->average.ImageFloatImage[k])/linfo->SD.ImageFloatImage[k];
//zscore = (d-ave)/sd;
linfo->Zscore.ImageFloatImage[k] = zscore*in->ImageFloatImage[k];
OSDN Git Service
