extern void arrayReadFromFile(FILE* fpt, Array* a, char* message);
extern void arrayWriteToFile(FILE* fpt, Array* a, char* message);
+// mode: 0
+// memoryAllocate u, C, lambda, ave;
+// 0<
+// Not Allocate
extern void arrayPCA(Array* u, Array* C, Array* lambda, Array* X, Array* ave, int mode);
extern void arrayPCATransform(Array* XX, Array* X, Array* ave, Array* u, int mode);
extern void arrayPCAPrint(FILE* fpt, Array* u, Array* lambda, int mode);
<H2>Information from source codes</H2>
</A>
<PRE>
-../src/lAllCommonLineCalculation.c:
-
-../src/lAngularWeightCalculation.c:
-
-../src/lCommonLineCalculation.c:
-
-../src/lCommonLineDataSetInPlaneRotation.c:
-
-../src/lCommonLinesDataSet.c:
-
-../src/lDescartesIntoPolar.c:
-
-../src/lEvaluateCorrelationMapwithCommonLine.c:
-
-../src/lFETOrientationSearchByAnnealing.c:
-
-../src/lFETOrientationSearchByFeatureAlignment.c:
-
-../src/lFETOrientationSearchByFeatureAlignmentInitSet.c:
-
-../src/lFETOrientationSearchByFeatureAlignmentSphere.c:
-
-../src/lFETmapOrientationSearchBySimultaneousFitting.c:
-
-../src/lFETsmallMapSetCreate_forSimultaneousMinimization.c:
-
-../src/lInitialDataFileRead.c:
-
-../src/lJcalulation.c:
-
-../src/lLcalculation.c:
-
-../src/lPlaneRotationAngleSet.c:
-
-../src/lllDataSeparationInfoPrint.c:
-
-../src/lllDataSeparationInfoRead.c:
-
-../src/lllDataSeparationInfoSet.c:
-
-../src/lllDatarMaxLimitForSeparation.c:
-
-../src/lllExtract.c:
-
-../src/lllExtractCtfinfFileCreate.c:
-
-../src/lllExtractWithSeparation.c:
-
-../src/lllExtractdYEstimate.c:
-
-../src/lmrc2Dto3D.c:
-
-../src/lmrc2Dto3DSIRT.c:
-
-../src/lmrc3Dto2D.c:
-
-../src/lmrcImageFileListAndEulerAngleDataRead.c:
-
-../src/lmrcImageOneLineGet.c:
-
-../src/lmrcImageSinogramCorrelationAdditionalWeight.c:
</PRE>
<HR>
<A NAME="include">
lmrc2Dto3DSimpleBackProjectionForEach(mrcImage* out, mrcImage* prj, lmrc2Dto3DInfo* linfo, Matrix3D Matrix, long mode)
{
mrcImageParaTypeReal x, y, z;
+ mrcImageParaTypeReal tmpx, tmpy, tmpz;
mrcImageParaTypeReal prjx, prjy, prjz;
double projdata, data;
mrcImageParaTypeReal gx, gy, gz;
mrcImageParaTypeReal g3x, g3y, g3z;
floatVector v;
+ int index;
floatVectorInit(&v, 4);
v.data[3] = 1.0;
gz = 0.0;
for(z=0; z<out->HeaderN.z; z++) {
+ v.data[2] = z-g3z;
for(y=0; y<out->HeaderN.y; y++) {
+ v.data[1] = y-g3y;
for(x=0; x<out->HeaderN.x; x++) {
+ index = x + y*out->HeaderN.x + z*out->HeaderN.x*out->HeaderN.y;
v.data[0] = x-g3x;
- v.data[1] = y-g3y;
- v.data[2] = z-g3z;
- matrix3DMultiplyVector(&v, Matrix);
- prjx = v.data[0] + gx;
- prjy = v.data[1] + gy;
- prjz = v.data[2] + g3z;
+ //matrix3DMultiplyVector(&v, Matrix);
+
+ tmpx = Matrix[0][0]*v.data[0]+ Matrix[1][0]*v.data[1] + Matrix[2][0]*v.data[2];
+ tmpy = Matrix[0][1]*v.data[0]+ Matrix[1][1]*v.data[1] + Matrix[2][1]*v.data[2];
+ tmpz = Matrix[0][2]*v.data[0]+ Matrix[1][2]*v.data[1] + Matrix[2][2]*v.data[2];
+
+ prjx = tmpx + gx;
+ prjy = tmpy + gy;
+ prjz = tmpz + g3z;
+ //prjx = v.data[0] + gx;
+ //prjy = v.data[1] + gy;
+ //prjz = v.data[2] + g3z;
+ //mrcPixelDataGet(out, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
- mrcPixelDataGet(out, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
if( -0.5<=prjx && prjx<prj->HeaderN.x-0.5
&&-0.5<=prjy && prjy<prj->HeaderN.y-0.5) {
+ data = out->ImageFloatImage[index];
mrcPixelDataGet(prj, prjx, prjy, 0.0, &projdata, mrcPixelRePart, linfo->InterpolationMode);
- mrcPixelDataSet(out, x, y, z, data + projdata, mrcPixelRePart);
+ //mrcPixelDataSet(out, x, y, z, data + projdata, mrcPixelRePart);
+ out->ImageFloatImage[index] = data + projdata;
}
}
}
lmrc2Dto3DFilteredBackProjectionForEach(mrcImage* out, mrcImage* prj, lmrc2Dto3DInfo* linfo, Matrix3D Matrix, long mode)
{
mrcImageParaTypeReal x, y, z;
+ mrcImageParaTypeReal tmpx, tmpy, tmpz;
+ mrcImageParaTypeReal tmpxx, tmpyx, tmpzx;
+ mrcImageParaTypeReal tmpxy, tmpyy, tmpzy;
+ mrcImageParaTypeReal tmpxz, tmpyz, tmpzz;
mrcImageParaTypeReal prjx, prjy, prjz;
double projdata, data, normz;
mrcImageParaTypeReal gx, gy, gz;
mrcImageParaTypeReal g3x, g3y, g3z;
floatVector v;
+ int index;
floatVectorInit(&v, 4);
v.data[3] = 1.0;
gz = 0.0;
for(z=0; z<out->HeaderN.z; z++) {
+ v.data[2] = z-g3z;
+ tmpxz = Matrix[2][0]*v.data[2];
+ tmpyz = Matrix[2][1]*v.data[2];
+ tmpzz = Matrix[2][2]*v.data[2];
+
for(y=0; y<out->HeaderN.y; y++) {
+ v.data[1] = y-g3y;
+ tmpxy = Matrix[1][0]*v.data[1];
+ tmpyy = Matrix[1][1]*v.data[1];
+ tmpzy = Matrix[1][2]*v.data[1];
+
for(x=0; x<out->HeaderN.x; x++) {
+ index = x + y*out->HeaderN.x + z*out->HeaderN.x*out->HeaderN.y;
v.data[0] = x-g3x;
- v.data[1] = y-g3y;
- v.data[2] = z-g3z;
- matrix3DMultiplyVector(&v, Matrix);
- prjx = v.data[0] + gx;
- prjy = v.data[1] + gy;
- prjz = v.data[2] + g3z;
+ tmpxx = Matrix[0][0]*v.data[0];
+ tmpyx = Matrix[0][1]*v.data[0];
+ tmpzx = Matrix[0][2]*v.data[0];
+ //matrix3DMultiplyVector(&v, Matrix);
+ tmpx = tmpxx + tmpxy + tmpxz;
+ tmpy = tmpyx + tmpyy + tmpyz;
+ tmpz = tmpzx + tmpzy + tmpzz;
+
+ //prjx = v.data[0] + gx;
+ //prjy = v.data[1] + gy;
+ //prjz = v.data[2] + g3z;
+ prjx = tmpx + gx;
+ prjy = tmpy + gy;
+ prjz = tmpz + g3z;
if(-0.5<=prjz && prjz < out->HeaderN.z - 0.5) {
- mrcPixelDataGet(out, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
+ //mrcPixelDataGet(out, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
if( -0.5<=prjx && prjx<prj->HeaderN.x-0.5
&&-0.5<=prjy && prjy<prj->HeaderN.y-0.5) {
+ data = out->ImageFloatImage[index];
mrcPixelDataGet(prj, prjx, prjy, 0.0, &projdata, mrcPixelRePart, linfo->InterpolationMode);
if(linfo->flagDouble && linfo->rhoInfo.flagThicknessWeight) {
- normz = fabs(v.data[2]/g3z);
+ normz = fabs(tmpz/g3z);
if(normz < 0.9) {
} else if(normz < 1) {
projdata = 0;
}
}
- mrcPixelDataSet(out, x, y, z, data + projdata, mrcPixelRePart);
+ // mrcPixelDataSet(out, x, y, z, data + projdata, mrcPixelRePart);
+ out->ImageFloatImage[index] = data + projdata;
}
}
}
extern void arrayReadFromFile(FILE* fpt, Array* a, char* message);
extern void arrayWriteToFile(FILE* fpt, Array* a, char* message);
+// mode: 0
+// memoryAllocate u, C, lambda, ave;
+// 0<
+// Not Allocate
extern void arrayPCA(Array* u, Array* C, Array* lambda, Array* X, Array* ave, int mode);
extern void arrayPCATransform(Array* XX, Array* X, Array* ave, Array* u, int mode);
extern void arrayPCAPrint(FILE* fpt, Array* u, Array* lambda, int mode);
#%Z%
*/
static char __sccs_id[] = "%Z%arrayPCA ver%I%; Date:%D% %Z%";
+#undef DEBUG
+#include "genUtil.h"
#include "../inc/Array.h"
/*
TmpX.n[1] = X->n[1];
arrayInit(&TmpX, "TmpX in arrayPCA");
- lambda->dim = 1;
- lambda->n[0] = X->n[0];
- arrayInit(lambda, "lambda in arrayPCA");
-
- average->dim = 1;
- average->n[0] = X->n[0];
- arrayInit(average, "average in arrayPCA");
+ if(0==mode) {
+ lambda->dim = 1;
+ lambda->n[0] = X->n[0];
+ arrayInit(lambda, "lambda in arrayPCA");
- C->dim = 2;
- C->n[0] = X->n[0];
- C->n[1] = X->n[0];
- arrayInit(C, "C in arrayPCA");
+ average->dim = 1;
+ average->n[0] = X->n[0];
+ arrayInit(average, "average in arrayPCA");
- u->dim = 2;
- u->n[0] = X->n[0];
- u->n[1] = X->n[0];
- arrayInit(u, "u in arrayPCA");
+ C->dim = 2;
+ C->n[0] = X->n[0];
+ C->n[1] = X->n[0];
+ arrayInit(C, "C in arrayPCA");
+ u->dim = 2;
+ u->n[0] = X->n[0];
+ u->n[1] = X->n[0];
+ arrayInit(u, "u in arrayPCA");
+ }
/* Calc TmpX : subtraction of average */
for(i=0; i<X->n[0]; i++) {
ave = 0;
}
}
jacobi(a, n, d, v, &nrot); //
- fprintf(stderr, "nrot for PCA: %d\n", nrot);
+#ifdef DEBUB
+ DEBUGPRINT1("nrot for PCA: %d\n", nrot);
+#endif
eigsrt(d, v, n); //
for(i=0; i<n; i++) {
arrayDataSet2(*u, i, j, v[i+1][j+1]);
}
}
+
+ arrayFree(&TmpX, "TmpX in arrayPCA");
}
clusterLogOneRecord * cluster;
cluster = in->current;
- return clusterLogWriteOneRecord(cluster, fpt, mode);
+ return clusterLogWriteOneRecordBinary2(cluster, fpt, mode);
}
clusterLogOneRecord*
<H2>Information from source codes</H2>
</A>
<PRE>
-../src/eosFunc.c:
</PRE>
<HR>
<A NAME="include">
#define DEBUG
#include "genUtil.h"
+#include "File.h"
#include "Cluster.h"
/*
{
clusterLog* cluster=NULL;
clusterLogASCII2BinaryInfo info;
+ FILE* fpt;
init0(&info);
argCheck(&info, argc, argv);
cluster = clusterLogReadAll(cluster, info.fptIn, 0);
clusterLogWriteAllBinary(cluster, info.Out, 0);
+
+ fpt = fileOpen(info.Out, "w");
+ clusterLogWriteAllBinary2(cluster, fpt, 0);
exit(EXIT_SUCCESS);
}
exec:
@echo "----- Execution Check -----"
- ../$(OSTYPE)/$(OBJECTNAME) -i test.log -o test.logbin
+ ../$(OSTYPE)/$(OBJECTNAME) -i data/test.log -o data/test.logbin
@echo "----- Calc check -----"
clean:
WORLDNAME=Tools
WORLDNAME=Tools
WORLDNAME=Tools
+WORLDNAME=Tools
+WORLDNAME=Tools
+WORLDNAME=Tools
+++ /dev/null
-include ../Config/Define.inc
-include Config/Define.inc
-include .Source
-
-include Config/Target.inc
--- /dev/null
+/Users/tacyas/Eos/src/Config/Template/ToolsHomeTemplate.Dir/Makefile
\ No newline at end of file
"-Ob","-O[utput]b[lack]","Output: black","Optional","1","1","OutBlack","outFile::mrcImage","NULL"
"-Owe","-O[utput]w[hite]e[dge]","Output: whiteEdge","Optional","1","1","OutWhiteEdge","outFile::mrcImage","NULL"
"-Obe","-O[utput]b[lack]e[dge]","Output: blackEdge","Optional","1","1","OutBlackEdge","outFile::mrcImage","NULL"
+"-lt","-l[ine]t[hickness]","[pixel]","Optional","2","1","LineThicknessX","Real","1","2","LineThicknessY","Real","1"
"-c","-c[onfig]","ConfigurationFile","Optional","1","1","configFile","inFile","NULL"
"-m","-m[ode]","Mode","Optional","1","1","mode","Integer","0"
char* OutBlackEdge;
FILE* fptOutBlackEdge;
+ long flagLineThicknessX;
+ float LineThicknessX;
+
+ long flagLineThicknessY;
+ float LineThicknessY;
+
long flagconfigFile;
char* configFile;
FILE* fptconfigFile;
}
SBREAK;
}
+ SCASE("lt") {
+ if(i+2<argc) {
+ info->LineThicknessX = stringGetNthRealData(argv[i+1], 1, " ,");
+ i++;
+ info->flagLineThicknessX++;
+ info->LineThicknessY = stringGetNthRealData(argv[i+1], 1, " ,");
+ i++;
+ info->flagLineThicknessY++;
+ } else {
+ usage(argv[0]);
+ exit(EXIT_FAILURE);
+ }
+ SBREAK;
+ }
SCASE("c") {
if(i+1<argc) {
info->configFile = stringGetNthWord(argv[i+1], 1, " ,");
info->fptOutBlack = NULL; info->flagOutBlack = 0;
info->fptOutWhiteEdge = NULL; info->flagOutWhiteEdge = 0;
info->fptOutBlackEdge = NULL; info->flagOutBlackEdge = 0;
+ info->LineThicknessX = 1; info->flagLineThicknessX = 0;
+ info->LineThicknessY = 1; info->flagLineThicknessY = 0;
info->fptconfigFile = NULL; info->flagconfigFile = 0;
info->mode = 0; info->flagmode = 0;
}
info->fptOutBlackEdge = fileOpen(info->OutBlackEdge, "w");
}
+ if(info->flagLineThicknessX) {
+ }
+
+ if(info->flagLineThicknessY) {
+ }
+
if(info->flagconfigFile) {
info->fptconfigFile = fileOpen(info->configFile, "r");
}
// Control: Input
mrcImageParaTypeRealCoord Length;
+ mrcImageParaTypeRealCoord LineThickness;
} lmrcImageMeanFreePathCalcInfo;
mrcFileRead(&in, info.In, "in main", 0);
linfo.Length = in.HeaderLength;
+ linfo.LineThickness.x = info.LineThicknessX;
+ linfo.LineThickness.y = info.LineThicknessY;
+ linfo.LineThickness.z = 1;
lmrcImageMeanFreePathCalc(&in, &linfo, info.mode);
int* countBlackEdge;
int* countWhiteEdge;
int x, y, z;
+ int srcx, srcy;
+ int minx, miny;
+ int maxx, maxy;
int flagEdge;
double data;
+ double data0;
double prevData;
int lengthWhite;
int lengthBlack;
linfo->N = in->HeaderN.z;
for(y=0; y<in->HeaderN.y; y++) {
+ miny = floor(y - linfo->LineThickness.y/2.0 + 0.5);
+ maxy = floor(y + linfo->LineThickness.y/2.0 + 0.5);
for(x=0; x<in->HeaderN.x; x++) {
+ minx = floor(x - linfo->LineThickness.x/2.0 + 0.5);
+ maxx = floor(x + linfo->LineThickness.x/2.0 + 0.5);
+
flagEdge = 1;
lengthWhite = 0;
lengthBlack = 0;
-
+
// Start
- mrcPixelDataGet(in, x, y, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
+ data = 1;
+ for(srcy=miny; srcy<maxy; srcy++) {
+ for(srcx=minx; srcx<maxx; srcx++) {
+ mrcPixelDataGet(in, srcx, srcy, 0, &data0, mrcPixelRePart, mrcPixelHowNearest);
+ data *= data0;
+ }
+ }
+ //mrcPixelDataGet(in, x, y, 0, &data, mrcPixelRePart, mrcPixelHowNearest);
if(0<data) {
countWhiteEdge[0]++;
lengthWhite++;
// Interval
for(z=1; z<in->HeaderN.z-1; z++) {
- mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
+ data = 1;
+ for(srcy=miny; srcy<maxy; srcy++) {
+ for(srcx=minx; srcx<maxx; srcx++) {
+ mrcPixelDataGet(in, srcx, srcy, z, &data0, mrcPixelRePart, mrcPixelHowNearest);
+ data *= data0;
+ }
+ }
+ //mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
+
if(0<prevData) {
if(0<data) { // White continue
if(flagEdge) {
// End
z = in->HeaderN.z-1;
- mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
+ data = 1;
+ for(srcy=miny; srcy<maxy; srcy++) {
+ for(srcx=minx; srcx<maxx; srcx++) {
+ mrcPixelDataGet(in, srcx, srcy, z, &data0, mrcPixelRePart, mrcPixelHowNearest);
+ data *= data0;
+ }
+ }
+ //mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
if(0<prevData) {
if(0<data) { // White continue
[-O[utput]b[lack] OutBlack (NULL ).as(outFile::mrcImage ) ] :Optional :Output: black
[-O[utput]w[hite]e[dge]OutWhiteEdge (NULL ).as(outFile::mrcImage ) ] :Optional :Output: whiteEdge
[-O[utput]b[lack]e[dge]OutBlackEdge (NULL ).as(outFile::mrcImage ) ] :Optional :Output: blackEdge
+ [-l[ine]t[hickness] LineThicknessX (1 ).as(Real )
+ LineThicknessY (1 ).as(Real ) ] :Optional :[pixel]
[-c[onfig] configFile (NULL ).as(inFile ) ] :Optional :ConfigurationFile
[-m[ode] mode (0 ).as(Integer ) ] :Optional :Mode
----- Additional Usage -----
-O 1 0 1 0 0 1 -1x1+1+6.000000 ' ' 'OutBlack' 'Output: black' Ob
-O 1 0 1 0 0 1 -1x1+1+7.500000 ' ' 'OutWhiteEdge' 'Output: whiteEdge' Owe
-O 1 0 1 0 0 1 -1x1+1+9.000000 ' ' 'OutBlackEdge' 'Output: blackEdge' Obe
- -I 1 0 1 0 0 1 -1x1+1+10.500000 ' ' 'configFile' 'ConfigurationFile' c
- -i 1 0 1 0 0 -1x1+1+12.000000 0 0 0 0 0 'mode' 'Mode' m
+ -I 1 0 0 1 0 1 -1x1+1+10 ' ' '1' 'LineThicknessX' [pixel]
+ -I 1 0 1 0 0 1 -1x1+1+12.000000 ' ' 'configFile' 'ConfigurationFile' c
+ -i 1 0 1 0 0 -1x1+1+13.500000 0 0 0 0 0 'mode' 'Mode' m
-E
-E
-E
../$(OSTYPE)/$(OBJECTNAME) -i data/test.bin -o data/test.mfp16 -m 16 -Ow data/test.white16 -Owe data/test.whiteEdge16 -Ob data/test.black16 -Obe data/test.blackEdge16
@echo "----- Calc check -----"
+exec2:
+ @echo "----- Execution Check -----"
+ ../$(OSTYPE)/$(OBJECTNAME) -i data/test.bin -o data/test.mfp2 -m 0 -Ow data/test.white2 -Owe data/test.whiteEdge2 -Ob data/test.black2 -Obe data/test.blackEdge2 -lt 2 2
+ ../$(OSTYPE)/$(OBJECTNAME) -i data/test.bin -o data/test.mfp216 -m 16 -Ow data/test.white216 -Owe data/test.whiteEdge216 -Ob data/test.black216 -Obe data/test.blackEdge216 -lt 2 2
+ @echo "----- Calc check -----"
+
init:
pdb2mrc -i data/121p.pdb2 -o data/test.mrc -nx 32 -ny 32 -nz 32 -dx 2 -dy 2 -dz 2 -Sx -32 -Sy -32 -Sz -32 -sig 1.6 -w 1.0
mrcImageBinalization -i data/test.mrc -o data/test.bin -m 32
fprintf(stderr, " [-O[utput]b[lack] OutBlack (NULL ).as(outFile::mrcImage ) ] :Optional :Output: black\n");
fprintf(stderr, " [-O[utput]w[hite]e[dge]OutWhiteEdge (NULL ).as(outFile::mrcImage ) ] :Optional :Output: whiteEdge\n");
fprintf(stderr, " [-O[utput]b[lack]e[dge]OutBlackEdge (NULL ).as(outFile::mrcImage ) ] :Optional :Output: blackEdge\n");
+ fprintf(stderr, " [-l[ine]t[hickness] LineThicknessX (1 ).as(Real ) \n LineThicknessY (1 ).as(Real ) ] :Optional :[pixel]\n");
fprintf(stderr, " [-c[onfig] configFile (NULL ).as(inFile ) ] :Optional :ConfigurationFile\n");
fprintf(stderr, " [-m[ode] mode (0 ).as(Integer ) ] :Optional :Mode\n");
additionalUsage();
#include "mrcImage.h"
typedef struct lmrcImageShapeSearchInfo {
- float radius; // Sylinder, half disk
+ float radius; // Sylinder, half disk, sphere
float minRadius;
float maxRadius;
float delRadius;
mrcImage SD; // SD for all
mrcImage Max; // Max for all
mrcImage Zscore; // Z-score;
+ mrcImage PCA; // PCA
int nShapeInfo;
mrcImage* shapeInfo; // Shape, Orientation, ...
} lmrcImageShapeSearchInfo;
typedef enum lmrcImageShapeSearchMode {
lmrcImageShapeSearchModeSylinder=0,
- lmrcImageShapeSearchModeHalfDisk=1
+ lmrcImageShapeSearchModeDisk=1,
+ lmrcImageShapeSearchModeSphere=2
} lmrcImageShapeSearchMode;
extern void lmrcImageShapeSearch(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 lmrcImageShapeSearchDisk(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode);
+extern void lmrcImageShapeSearchSphere(double* data, mrcImage* in, lmrcImageShapeSearchInfo* linfo, int mode);
extern void lmrcImageShapeSearchModePrint(FILE* fpt);
int
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);
for(i=0; i<linfo.nShapeInfo; i++) {
sprintf(filename, "%s.shapeinfo-%02d", info.Out, i);
mrcFileWrite(&linfo.shapeInfo[i], filename, "in main", 0);
lmrcImageShapeSearchModePrint(FILE* fpt)
{
fprintf(fpt, "%d: Cylinder with radius and length\n", lmrcImageShapeSearchModeSylinder);
- fprintf(fpt, "%d: Half disk with radius and length(thickness)\n", lmrcImageShapeSearchModeHalfDisk);
+ fprintf(fpt, "%d: Disk with radius and length(thickness)\n", lmrcImageShapeSearchModeDisk);
+ fprintf(fpt, "%d: Sphere with radius \n", lmrcImageShapeSearchModeSphere);
}
lmrcImageShapeSearchCalc0(out, in, linfo, mode);
break;
}
- case lmrcImageShapeSearchModeHalfDisk: {
+ case lmrcImageShapeSearchModeDisk: {
+ lmrcImageShapeSearchCalc0(out, in, linfo, mode);
+ break;
+ }
+ case lmrcImageShapeSearchModeSphere: {
lmrcImageShapeSearchCalc0(out, in, linfo, mode);
break;
}
{
float x, y, z;
float srcx, srcy, srcz;
+ float srcx0, srcy0, srcz0;
float dstx, dsty, dstz;
float theta, phi, psi, radius, length;
int iPhi, iTheta, iPsi, iRadius, iLength;
floatVector v;
int i, imax;
int k, k0, nVoxel;
+ int PCAMode=0;
+ int numPoint;
+ Array U;
+ Array C;
+ Array Lambda;
+ double L0, L1, L2;
+ Array X;
+ Array Ave;
+ mrcImageParaTypeInteger nCube;
+ double W;
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->nShapeInfo = 5;
linfo->shapeInfo = (mrcImage*)memoryAllocate(sizeof(mrcImage)*linfo->nShapeInfo, "in MFPCalc");
for(i=0; i<linfo->nShapeInfo; i++) {
mrcInit(&linfo->shapeInfo[i], NULL);
}
+ X.dim = 2;
+ X.n[0] = 3;
+ nCube = ((int)(MAX(linfo->maxLength, 2*linfo->maxRadius+1)+3)/2)*2;
+ X.n[1] = nCube*nCube*nCube;
+ arrayInit(&X, "in");
+
+ DEBUGPRINT1("nCube: %d\n", nCube);
+
for(z=0; z<in->HeaderN.z; z++) {
DEBUGPRINT1("z: %f\n", z);
linfo->z = z;
linfo->x = x;
k0 = (int)(x + y*in->HeaderN.x + z*in->HeaderN.x*in->HeaderN.y);
+
+ //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) {
+ 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;
+ d = in->ImageFloatImage[k];
+ if(0<d) {
+ arrayDataSet2(X, 0, numPoint, srcx0);
+ arrayDataSet2(X, 1, numPoint, srcy0);
+ arrayDataSet2(X, 2, numPoint, srcz0);
+ numPoint++;
+ }
+ }
+ }
+ }
+ //DEBUGPRINT4("%f %f %f numPoint: %d\n", x, y, z, numPoint);
+ if(16<numPoint) {
+ X.n[1] = numPoint;
+ arrayPCA(&U, &C, &Lambda, &X, &Ave, PCAMode);
+ if(PCAMode==0) PCAMode=1;
+ L0 = arrayDataGet1(Lambda, 0);
+ L1 = arrayDataGet1(Lambda, 1);
+ L2 = arrayDataGet1(Lambda, 2);
+
+ switch(mode) {
+ case lmrcImageShapeSearchModeSylinder: {
+ W = ((L1!=0)?(L0/L1):1)*(L2!=0?(L0/L2):1);
+#ifdef DEBUG
+ if(4000<W) {
+ DEBUGPRINT3("%f %f %f\n", x, y, z);
+ DEBUGPRINT4("W: %f %f %f %f\n", W, L0, L1, L2);
+ }
+#endif
+ break;
+ }
+ case lmrcImageShapeSearchModeDisk: {
+ W = ((L2!=0)?((L0/L2)*(L1/L2)):1);
+#ifdef DEBUG
+ if(1000<W) {
+ DEBUGPRINT3("%f %f %f\n", x, y, z);
+ DEBUGPRINT4("W: %f %f %f %f\n", W, L0, L1, L2);
+ }
+#endif
+ break;
+ }
+ case lmrcImageShapeSearchModeSphere: {
+ W = L0*L1+L1*L2+L2*L0/SQR(L0+L1+L2);
+ break;
+ }
+ default: {
+ W = 1;
+ break;
+ }
+ }
+ } else {
+ W = 1;
+ }
+ } else {
+ W = 1;
+ }
+ //mrcPixelDataSet(&linfo->PCA, x, y, z, W, mrcPixelRePart);
+ linfo->PCA.ImageFloatImage[k] = W;
+ if(1<W) {
for(iPhi=0; iPhi<nPhi; iPhi++) { // Three Rotation
phi = linfo->minPhi + iPhi*linfo->delPhi;
for(iTheta=0; iTheta<nTheta; iTheta++) {
#endif
break;
}
- case lmrcImageShapeSearchModeHalfDisk: {
- lmrcImageShapeSearchHalfDisk(&data[index], in, mat, linfo, mode);
+ case lmrcImageShapeSearchModeDisk: {
+ lmrcImageShapeSearchDisk(&data[index], in, mat, linfo, mode);
+ break;
+ }
+ case lmrcImageShapeSearchModeSphere: {
+ lmrcImageShapeSearchSphere(&data[index], in, mat, linfo, mode);
break;
}
default: {
DEBUGPRINT3("%f %f %f \n", max, ave, sd);
}
#endif
+ } else {
+ ave = 0;
+ sd = 0;
+ max = 0;
+ imax = 0;
+ }
+
//mrcPixelDataSet(&linfo->average, x, y, z, ave, mrcPixelRePart);
linfo->average.ImageFloatImage[k0] = ave;
linfo->shapeInfo[3].ImageFloatImage[k0] = i*linfo->delRadius + linfo->minRadius;
i = imax/(nPhi*nTheta*nPsi*nRadius);
- mrcPixelDataSet(&linfo->shapeInfo[4], x, y, z, i*linfo->delLength + linfo->minLength, mrcPixelRePart);
+ //mrcPixelDataSet(&linfo->shapeInfo[4], x, y, z, i*linfo->delLength + linfo->minLength, mrcPixelRePart);
linfo->shapeInfo[4].ImageFloatImage[k0] = i*linfo->delLength + linfo->minLength;
}
}
}
void
-lmrcImageShapeSearchHalfDisk(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode)
+lmrcImageShapeSearchDisk(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode)
{
- DEBUGPRINT("Not yet supported\n");
+
+ double p, l, r;
+ double x, y, z;
+ double zz;
+ 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(r=0; r<=linfo->radius; r+=linfo->delRadius) {
+ x = r*cp;
+ y = r*sp;
+ for(l=0; l<=linfo->length; l+=linfo->delLength) {
+ z = l-linfo->length/2;
+
+ 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.5<dstx && dstx < in->HeaderN.x-0.5
+ &&-0.5<dsty && dsty < in->HeaderN.y-0.5
+ &&-0.5<dstz && dstz < in->HeaderN.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
+ for(zz = -1; zz <= 1; zz+=2) {
+ x = r*cos(p);
+ y = r*sin(p);
+ z = l - (linfo->length/2 + linfo->delLength)*zz;
+ 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.5<dstx && dstx < in->HeaderN.x-0.5
+ &&-0.5<dsty && dsty < in->HeaderN.y-0.5
+ &&-0.5<dstz && dstz < in->HeaderN.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);
+#ifdef DEBUG2
+ //DEBUGPRINT1("data: %f\n", score);
+ if(0<score) {
+ DEBUGPRINT1("data: %f\n", score);
+ }
+#endif
+}
+
+
+void
+lmrcImageShapeSearchSphere(double* data, mrcImage* in, Matrix3D mat, lmrcImageShapeSearchInfo* linfo, int mode)
+{
+
+ double phi, psi, r;
+ double x, y, z;
+ double zz;
+ 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(phi=0; phi<2*M_PI; phi+=linfo->delPhi) {
+ for(psi=0; psi<2*M_PI; psi+=linfo->delPsi) {
+ for(r=0; r<=linfo->radius; r+=linfo->delRadius) {
+
+ countPos++;
+ scorePos += d;
+ }
+ // Edge Check
+ 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);
}
+
[-m[ode] mode (0 ).as(Integer ) ] :Optional :Mode
----- Additional Usage -----
0: Cylinder with radius and length
-1: Half disk with radius and length(thickness)
+1: Disk with radius and length(thickness)
</PRE>
</BODY>
</HTML>
include ../../../../Config/Define.inc
include ../../../../../Config/Define.inc
-all: help exec
+all: help exec exec2
help:
@echo "----- Help Message Check -----"
exec:
@echo "----- Execution Check -----"
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
+
+exec2:
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test-disk.bin -o data/test-disk.disk -r 6 6 1 -l 3 5 1 -Phi 0 359 45 -Theta 0 179 45 -Psi 0 0 90 -m 1 -thres 0.5
+ @echo "----- Calc check -----"
+
+exec3:
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test-disk.bin -o data/test-disk.disk -r 6 6 1 -l 3 5 1 -Phi 0 359 45 -Theta 0 179 45 -Psi 0 0 90 -m 2 -thres 0.5
@echo "----- Calc check -----"
init:
mrcImageModelCreate -o data/test-sylinder.mrc -nx 32 -ny 32 -nz 32 -r 3 -l 20 -m 8
mrcImageBinalization -i data/test-sylinder.mrc -o data/test-sylinder.bin -m 32
+ mrcImageModelCreate -o data/test-disk.mrc -nx 32 -ny 32 -nz 32 -r 10 -l 3 -m 8
+ mrcImageBinalization -i data/test-disk.mrc -o data/test-disk.bin -m 32
clean:
include ../../../../Config/Define.inc
include ../../../../../Config/Define.inc
-all: help exec
+all: help exec exec2
help:
@echo "----- Help Message Check -----"
time ../$(OSTYPE)/$(OBJECTNAME) -i data/test.mrc2d3d -o data/test.mrc2d3d.mrc3d -Double -CounterThreshold 1 -WeightMode 3 -m 1 -InterpolationMode 2 -DoubleCounter data/test.mrc2d3d.mrc3d.counter 2> /dev/null
time ~/Eos/bin/$(OSTYPE)/$(OBJECTNAME) -i data/test.mrc2d3d -o data/test.mrc2d3d.mrc3d -Double -CounterThreshold 1 -WeightMode 3 -m 1 -InterpolationMode 2 -DoubleCounter data/test.mrc2d3d.mrc3d.counter 2> /dev/null
+exec2:
+ @echo "----- Execution Check 2-----"
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test2.mrc2d3d -o data/test2.mrc2d3d.mrc3d -Double -CounterThreshold 1 -WeightMode 3 -m 1 -InterpolationMode 2 -DoubleCounter data/test2.mrc2d3d.mrc3d.counter 2> /dev/null
+ time ~/Eos/bin/$(OSTYPE)/$(OBJECTNAME) -i data/test2.mrc2d3d -o data/test2.mrc2d3d.mrc3d -Double -CounterThreshold 1 -WeightMode 3 -m 1 -InterpolationMode 2 -DoubleCounter data/test2.mrc2d3d.mrc3d.counter 2> /dev/null
+
exec-pthread:
@echo "----- Execution Check -----"
time ../$(OSTYPE)/$(OBJECTNAME) -i data/test.mrc2d3d -o data/test.mrc2d3d.mrc3d-pthread -Double -CounterThreshold 1 -WeightMode 3 -m 1 -InterpolationMode 2 -DoubleCounter data/test.mrc2d3d.mrc3d.counter-pthread -pthread 8
pdb2mrc -i data/1vomTrans.pdb -o data/test2.mrc3d -nx 81 -ny 81 -nz 81 -dx 2.5 -dy 2.5 -dz 2.5 -Sx -100 -Sy -100 -Sz -100 -sig 2.5
pdb2mrc2d -i data/1vomTrans.pdb -o data/test2.mrc2d -d 2.5 2.5 -s -100 -100 -Rot 30 30 -n 81 81 -O data/test2.mrc2d3d -m 1 -sig 2.5
+prepare-all:
+ # GC: 6.727136e+01 3.631518e+01 -2.983928e+01
+ pdb2mrc -i data/1vomTrans.pdb -o data/test00.mrc3d -nx 10 -ny 10 -nz 10 -dx 20 -dy 20 -dz 20 -Sx -100 -Sy -100 -Sz -100 -sig 2.5
+ pdb2mrc2d -i data/1vomTrans.pdb -o data/test00.mrc2d -d 20 20 -s -100 -100 -Rot 10 10 -n 10 10 -O data/test00.mrc2d3d -m 1 -sig 2.5
+ pdb2mrc -i data/1vomTrans.pdb -o data/test01.mrc3d -nx 20 -ny 20 -nz 20 -dx 10 -dy 10 -dz 20 -Sx -100 -Sy -100 -Sz -100 -sig 2.5
+ pdb2mrc2d -i data/1vomTrans.pdb -o data/test01.mrc2d -d 10 10 -s -100 -100 -Rot 10 10 -n 20 20 -O data/test01.mrc2d3d -m 1 -sig 2.5
+ pdb2mrc -i data/1vomTrans.pdb -o data/test02.mrc3d -nx 40 -ny 40 -nz 40 -dx 5 -dy 5 -dz 5 -Sx -100 -Sy -100 -Sz -100 -sig 2.5
+ pdb2mrc2d -i data/1vomTrans.pdb -o data/test02.mrc2d -d 5 5 -s -100 -100 -Rot 10 10 -n 40 40 -O data/test02.mrc2d3d -m 1 -sig 2.5
+ pdb2mrc -i data/1vomTrans.pdb -o data/test03.mrc3d -nx 80 -ny 80 -nz 80 -dx 2.5 -dy 2.5 -dz 2.5 -Sx -100 -Sy -100 -Sz -100 -sig 2.5
+ pdb2mrc2d -i data/1vomTrans.pdb -o data/test03.mrc2d -d 2.5 2.5 -s -100 -100 -Rot 10 10 -n 80 80 -O data/test03.mrc2d3d -m 1 -sig 2.5
+ pdb2mrc -i data/1vomTrans.pdb -o data/test04.mrc3d -nx 160 -ny 160 -nz 160 -dx 1.25 -dy 1.25 -dz 1.25 -Sx -100 -Sy -100 -Sz -100 -sig 2.5
+ pdb2mrc2d -i data/1vomTrans.pdb -o data/test04.mrc2d -d 1.25 1.25 -s -100 -100 -Rot 10 10 -n 160 160 -O data/test04.mrc2d3d -m 1 -sig 2.5
+
+exec-all:
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test00.mrc2d3d -o data/test00.mrc2d3d.mrc3d -m 1 -InterpolationMode 2 -single 1 2> /dev/null
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test01.mrc2d3d -o data/test01.mrc2d3d.mrc3d -m 1 -InterpolationMode 2 -single 1 2> /dev/null
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test02.mrc2d3d -o data/test02.mrc2d3d.mrc3d -m 1 -InterpolationMode 2 -single 1 2> /dev/null
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test03.mrc2d3d -o data/test03.mrc2d3d.mrc3d -m 1 -InterpolationMode 2 -single 1 2> /dev/null
+ time ../$(OSTYPE)/$(OBJECTNAME) -i data/test04.mrc2d3d -o data/test04.mrc2d3d.mrc3d -m 1 -InterpolationMode 2 -single 1 2> /dev/null
clean:
OSDN Git Service
