Please enter the commit message for your changes. Lines starting

 with '#' will be ignored, and an empty message aborts the commit.

 On branch master
 Your branch is up to date with 'origin/master'.

 Changes to be committed:
	modified:   src/Objects/DataManip/mrcImage/src/lmrcImageCrystalCreate.c

diff --git a/src/Objects/DataManip/mrcImage/src/lmrcImageCrystalCreate.c b/src/Objects/DataManip/mrcImage/src/lmrcImageCrystalCreate.c
index eb0129d..8674713 100755 (executable)
--- a/src/Objects/DataManip/mrcImage/src/lmrcImageCrystalCreate.c
+++ b/src/Objects/DataManip/mrcImage/src/lmrcImageCrystalCreate.c
@@ -17,12 +17,12 @@
 void 
 lmrcImageCrystalCreate(mrcImage* out, mrcImage* in, lmrcImageCrystalCreateInfo* info)
 {
-       mrcImageParaTypeReal x, y, z;
+       mrcImageParaTypeReal xout, yout, zout;
        mrcImageParaTypeReal orgx, orgy, orgz;
-       mrcImageParaTypeReal realx, realy, realz;
+       mrcImageParaTypeReal realxout, realyout, realzout;
        mrcImageParaTypeReal ix, iy, iz;
-       mrcImageParaTypeReal da, db, dc;
-       mrcImageParaTypeReal orgda, orgdb, orgdc;
+       mrcImageParaTypeReal na, nb, nc;
+       mrcImageParaTypeReal orgna, orgnb, orgnc;
        mrcImageParaTypeReal La, Lb, Lc;
   double det, Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz;
 
@@ -64,6 +64,8 @@ lmrcImageCrystalCreate(mrcImage* out, mrcImage* in, lmrcImageCrystalCreateInfo*
   } else {
     DEBUGPRINT1("Det of A, B, C: %f\n", det);
   }
+
+  // Inverse of (a b c)
   Ax =  (info->by*info->cz - info->cy*info->bz)/det;
   Ay = -(info->bx*info->cz - info->cx*info->bz)/det;
   Az =  (info->bx*info->cy - info->cx*info->by)/det;
@@ -79,19 +81,37 @@ lmrcImageCrystalCreate(mrcImage* out, mrcImage* in, lmrcImageCrystalCreateInfo*
   DEBUGPRINT3("C-1: %f %f %f\n", Cx, Cy, Cz);  
   DEBUGPRINT3("L: %f %f %f\n", La, Lb, Lc);
 
-       for(z=0; z<out->HeaderN.z; z++) {
-       for(y=0; y<out->HeaderN.y; y++) {
-       for(x=0; x<out->HeaderN.x; x++) {
-    realx = (x+out->HeaderStartN.x)*in->HeaderLength.x - info->sx;
-    realy = (y+out->HeaderStartN.y)*in->HeaderLength.y - info->sy;
-    realz = (z+out->HeaderStartN.z)*in->HeaderLength.z - info->sz;
+  // out->StartN
+  out->HeaderStartN.x = round(info->sx/out->HeaderLength.x); 
+  out->HeaderStartN.y = round(info->sy/out->HeaderLength.y); 
+  out->HeaderStartN.z = round(info->sz/out->HeaderLength.z); 
+
+       for(zout=0; zout<out->HeaderN.z; zout++) {
+       for(yout=0; yout<out->HeaderN.y; yout++) {
+       for(xout=0; xout<out->HeaderN.x; xout++) {
+  /*
+    xout      ax      bx      cx
+    yout = na ay + nb by + nc cy 
+    zout      az      bz      cz
+
+        ax bx cx    xout      na
+    A = ay by cy    yout  = A nb 
+        az bz cz ,  zout      nc
 
+    na        xout
+    nb = A^-1 yout
+    nc        zout
+
+  */
+    realxout = xout*in->HeaderLength.x;
+    realyout = yout*in->HeaderLength.y;
+    realzout = zout*in->HeaderLength.z;
     
-    da = Ax*realx + Ay*realy + Az*realz; 
-    db = Bx*realx + By*realy + Bz*realz;
-    dc = Cx*realx + Cy*realy + Cz*realz;
-    //DEBUGPRINT3("d: %f %f %f\n", da, db, dc);
-    DEBUGPRINT6("d: %f %f %f at %f %f %f\n", da, db, dc, realx, realy, realz);
+    na = Ax*realxout + Ay*realyout + Az*realzout; 
+    nb = Bx*realxout + By*realyout + Bz*realzout;
+    nc = Cx*realxout + Cy*realyout + Cz*realzout;
+    //DEBUGPRINT3("d: %f %f %f\n", na, nb, nc);
+    DEBUGPRINT6("n: %f %f %f at %f %f %f\n", na, nb, nc, realxout, realyout, realzout);
 
     /*
     if(0<=da && da<info->nx 
@@ -105,17 +125,20 @@ lmrcImageCrystalCreate(mrcImage* out, mrcImage* in, lmrcImageCrystalCreateInfo*
     flagInCrystal = 1;
 
     if(flagInCrystal) { // in cyrstal
-        orgda = da - floor(da);
-        orgdb = db - floor(db);
-        orgdc = dc - floor(dc);
-
-        orgx = (orgda*info->ax + orgdb*info->bx + orgdc*info->cx)/in->HeaderLength.x - in->HeaderStartN.x;  
-        orgy = (orgda*info->ay + orgdb*info->by + orgdc*info->cy)/in->HeaderLength.y - in->HeaderStartN.y;  
-        orgz = (orgda*info->az + orgdb*info->bz + orgdc*info->cz)/in->HeaderLength.z - in->HeaderStartN.z;  
-        DEBUGPRINT6("org: %f %f %f at %f %f %f\n", orgda, orgdb, orgdc, orgx, orgy, orgz);
+        // Position of In Unit Cell
+        orgna = na - floor(na);
+        orgnb = nb - floor(nb);
+        orgnc = nc - floor(nc);
+        
+        // pixel                                       Angstrom  + start              -> pixel    -> Start              
+        orgx = (orgna*info->ax + orgnb*info->bx + orgnc*info->cx + info->sx)/in->HeaderLength.x - in->HeaderStartN.x;  
+        orgy = (orgna*info->ay + orgnb*info->by + orgnc*info->cy + info->sy)/in->HeaderLength.y - in->HeaderStartN.y;  
+        orgz = (orgna*info->az + orgnb*info->bz + orgnc*info->cz + info->sz)/in->HeaderLength.z - in->HeaderStartN.z;  
+        
+        DEBUGPRINT6("n(unit): %f %f %f coord(unit):  %f %f %f\n", orgna, orgnb, orgnc, orgx, orgy, orgz);
 
                                mrcPixelDataGet(in,  orgx, orgy, orgz, &data, mrcPixelRePart, info->pixelMode);
-                               mrcPixelDataSet(out,    x,    y,    z,  data, mrcPixelRePart);
+                               mrcPixelDataSet(out, xout, yout, zout,  data, mrcPixelRePart);
         /*
                                for(ix=0; ix<info->nx; ix++) {
                                for(iy=0; iy<info->ny; iy++) {