src/Objects/DataManip/mrcImage/src/lmrcImageMasking.c

   1 /*
   2 # lmrcImageMasking.c  1.1
   3 # The latest update : 01/23/97 at 10:38:49
   4 #
   5 #@(#) lmrcImageMasking ver 1.1
   6 #@(#) Created by
   7 #@(#)
   8 #@(#) Usage : lmrcImageMasking
   9 #@(#) Attention
  10 #@(#)
  11 */
  12 static char __sccs_id[] = "@(#)lmrcImageMasking ver1.1; Date:97/01/23 @(#)";
  13
  14 #include <stdio.h>
  15 #include <stdlib.h>
  16
  17 #undef DEBUG
  18 #include "genUtil.h"
  19 #include "mrcImage.h"
  20 #include "lmrcImageMasking.h"
  21
  22 void
  23 lmrcImageMasking(mrcImage* out, mrcImage* in, lmrcImageMaskingInfo* linfo, int mode)
  24 {
  25         mrcImageParaTypeInteger x, y, z;
  26         double data;
  27         double sum;
  28         int count;
  29         double avg;
  30         mrcImageParaTypeReal sx, sy, sz;
  31         floatVector nx, ny, nz, cx0, cx1, cy0, cy1, cz0, cz1;
  32         floatVector tmp;
  33         Matrix3D mat;
  34         Matrix3D matRotAnti;
  35         Matrix3D matRot;
  36         float x0, x1, y0, y1, z0, z1;
  37     int flagIn = 0;
  38     int flagAvg = 0;
  39     int orgData;
  40
  41         out->Header =   in->Header;
  42         mrcInit (out, NULL);
  43
  44         if(linfo->flagRotation) {
  45                 floatVectorInit(&nx, 4);
  46                 floatVectorInit(&ny, 4);
  47                 floatVectorInit(&nz, 4);
  48                 floatVectorInit(&cx0, 4);
  49                 floatVectorInit(&cx1, 4);
  50                 floatVectorInit(&cy0, 4);
  51                 floatVectorInit(&cy1, 4);
  52                 floatVectorInit(&cz0, 4);
  53                 floatVectorInit(&cz1, 4);
  54                 floatVectorInit(&tmp, 4);
  55
  56                 nx.data[0] = 1; nx.data[1] = 0; nx.data[2] = 0; nx.data[3] = 1;
  57                 ny.data[0] = 0; ny.data[1] = 1; ny.data[2] = 0; ny.data[3] = 1;
  58         nz.data[0] = 0; nz.data[1] = 0; nz.data[2] = 1; nz.data[3] = 1;
  59
  60                 cx0.data[0] = -linfo->n.x/2.0;   cx0.data[1] = 0; cx0.data[2] = 0; cx0.data[3] = 1;
  61                 cx1.data[0] = +linfo->n.x/2.0;   cx1.data[1] = 0; cx1.data[2] = 0; cx1.data[3] = 1;
  62                 cy0.data[0] = 0;   cy0.data[1] = -linfo->n.y/2.0; cy0.data[2] = 0; cy0.data[3] = 1;
  63                 cy1.data[0] = 0;   cy1.data[1] = +linfo->n.y/2.0; cy1.data[2] = 0; cy1.data[3] = 1;
  64                 cz0.data[0] = 0;   cz0.data[1] = 0; cz0.data[2] = -linfo->n.z/2.0; cz0.data[3] = 1;
  65                 cz1.data[0] = 0;   cz1.data[1] = 0; cz1.data[2] = +linfo->n.z/2.0; cz1.data[3] = 1;
  66
  67                 tmp.data[0] = 0; tmp.data[1] = 0; tmp.data[2] = 0; tmp.data[3] = 1;
  68
  69                 matrix3DRotationAntiSetFollowingEulerAngle(matRotAnti, linfo->Euler, linfo->Rot1, linfo->Rot2, linfo->Rot3,  MATRIX_3D_MODE_INITIALIZE);
  70                 matrix3DRotationSetFollowingEulerAngle    (matRot,     linfo->Euler, linfo->Rot1, linfo->Rot2, linfo->Rot3,  MATRIX_3D_MODE_INITIALIZE);
  71
  72                 matrix3DRotationSetFollowingEulerAngle(mat, linfo->Euler, linfo->Rot1, linfo->Rot2, linfo->Rot3,  MATRIX_3D_MODE_INITIALIZE);
  73                 matrix3DTranslationSet(mat, linfo->c.x , linfo->c.y, linfo->c.z, MATRIX_3D_MODE_NOT_INITIALIZE);
  74
  75                 matrix3DMultiplyVector(&nx, matRot); DEBUGPRINT3("nx %f %f %f\n", nx.data[0], nx.data[1], nx.data[2]);
  76                 matrix3DMultiplyVector(&ny, matRot); DEBUGPRINT3("ny %f %f %f\n", ny.data[0], ny.data[1], ny.data[2]);
  77                 matrix3DMultiplyVector(&nz, matRot); DEBUGPRINT3("nz %f %f %f\n", nz.data[0], nz.data[1], nz.data[2]);
  78                 matrix3DMultiplyVector(&cx0, mat);   DEBUGPRINT3("cx0 %f %f %f\n", cx0.data[0], cx0.data[1], cx0.data[2]);
  79                 matrix3DMultiplyVector(&cy0, mat);   DEBUGPRINT3("cy0 %f %f %f\n", cy0.data[0], cy0.data[1], cy0.data[2]);
  80                 matrix3DMultiplyVector(&cz0, mat);   DEBUGPRINT3("cz0 %f %f %f\n", cz0.data[0], cz0.data[1], cz0.data[2]);
  81                 matrix3DMultiplyVector(&cx1, mat);   DEBUGPRINT3("cx1 %f %f %f\n", cx1.data[0], cx1.data[1], cx1.data[2]);
  82                 matrix3DMultiplyVector(&cy1, mat);   DEBUGPRINT3("cy1 %f %f %f\n", cy1.data[0], cy1.data[1], cy1.data[2]);
  83                 matrix3DMultiplyVector(&cz1, mat);   DEBUGPRINT3("cz1 %f %f %f\n", cz1.data[0], cz1.data[1], cz1.data[2]);
  84
  85         }
  86
  87     if(linfo->mode==1) {
  88             orgData = 1;
  89     } else {
  90         orgData = 0;
  91     }
  92
  93         DEBUGPRINT3("shape %d mode %d rotation %d\n", linfo->shape, linfo->mode, linfo->flagRotation);
  94         switch(linfo->shape) {
  95                 case lmrcImageMaskingInfoShapeSphere:
  96                         DEBUGPRINT("lmrcImageMaskingInfoShapeSphere\n");
  97                         sum      = 0.0; count = 0; avg = 0;
  98                         if (linfo->mode==0){ // Calculate Average density outside a sphere
  99                                 for(z=0; z<out->HeaderN.z; z++) {
 100                                 for(y=0; y<out->HeaderN.y; y++) {
 101                                 for(x=0; x<out->HeaderN.x; x++) {
 102                                         if(!linfo->flagRotation) {
 103                                                 if(SQR((x - linfo->c.x )/(linfo->n.x/2.0))
 104                                                  + SQR((y - linfo->c.y )/(linfo->n.y/2.0))
 105                                                  + SQR((z - linfo->c.z )/(linfo->n.z/2.0)) >= 1) { // Outside
 106                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 107                                                         sum += data;
 108                                                         count ++;
 109                                                 }
 110                                         } else {
 111                                                 tmp.data[0] = x - linfo->c.x;
 112                                                 tmp.data[1] = y - linfo->c.y;
 113                                                 tmp.data[2] = z - linfo->c.z;
 114                                                 matrix3DMultiplyVector(&tmp, matRotAnti);
 115                                                 if(SQR((tmp.data[0])/(linfo->n.x/2.0))
 116                                                  + SQR((tmp.data[1])/(linfo->n.y/2.0))
 117                                                  + SQR((tmp.data[2])/(linfo->n.z/2.0)) >= 1) {    // Outside
 118                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 119                                                         sum += data;
 120                                                         count ++;
 121                                                 }
 122                                         }
 123                                 }
 124                                 }
 125                                 }
 126                                 if(0<count) {
 127                                         avg = sum/count;
 128                                 }
 129                         }
 130                         DEBUGPRINT2("average %f %d\n", avg, count);
 131
 132                         for(z=0; z<out->HeaderN.z; z++) {
 133                         for(y=0; y<out->HeaderN.y; y++) {
 134                         for(x=0; x<out->HeaderN.x; x++) {
 135                                 data = orgData;
 136                 flagIn = 0;
 137                                 if(!linfo->flagRotation) {
 138                                         if(SQR((x - linfo->c.x )/(linfo->n.x/2.0))
 139                                          + SQR((y - linfo->c.y )/(linfo->n.y/2.0))
 140                                          + SQR((z - linfo->c.z )/(linfo->n.z/2.0)) <= 1) { // Inside
 141                                                 mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 142                         flagIn = 1;
 143                                         }
 144                                 } else {
 145                                         tmp.data[0] = x -linfo->c.x;
 146                                         tmp.data[1] = y -linfo->c.y;
 147                                         tmp.data[2] = z -linfo->c.z;
 148                                         matrix3DMultiplyVector(&tmp, matRotAnti);
 149                                         if(SQR((tmp.data[0])/(linfo->n.x/2.0))
 150                                          + SQR((tmp.data[1])/(linfo->n.y/2.0))
 151                                          + SQR((tmp.data[2])/(linfo->n.z/2.0)) <= 1) {     // Inside
 152                                                 mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 153                         flagIn=1;
 154                                         }
 155                                 }
 156                 if(flagIn==1) {
 157                                         if(linfo->mode==0) {
 158                                                 data -= avg;
 159                                         } else if(linfo->mode==1)  {
 160                                                 data = 0;
 161                                         } else if(linfo->mode==2) {
 162                         data = 1;
 163                     }
 164                 }
 165                                 mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 166                         }
 167                         }
 168                         }
 169                         break;
 170
 171                 case lmrcImageMaskingInfoSquare:
 172                         sum      = 0.0; count = 0; avg = 0;
 173                         if (linfo->mode==0){
 174                                 for(z=0; z<out->HeaderN.z; z++) {
 175                                 for(y=0; y<out->HeaderN.y; y++) {
 176                                 for(x=0; x<out->HeaderN.x; x++) {
 177                                         if(!linfo->flagRotation) {
 178                                                 if(     fabs(x - linfo->c.x ) > (linfo->n.x/2.0)
 179                                                  || fabs(y - linfo->c.y ) > (linfo->n.y/2.0)
 180                                                  || fabs(z - linfo->c.z ) > (linfo->n.z/2.0) ) { // Outside
 181                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 182                                                         sum += data;
 183                                                         count++;
 184                                                 }
 185                                         } else {
 186                                                 x0 = (nx.data[0]*(x-cx0.data[0])+nx.data[1]*(y-cx0.data[1])+nx.data[2]*(z-cx0.data[2]));
 187                                                 x1 = (nx.data[0]*(x-cx1.data[0])+nx.data[1]*(y-cx1.data[1])+nx.data[2]*(z-cx1.data[2]));
 188                                                 y0 = (ny.data[0]*(x-cy0.data[0])+ny.data[1]*(y-cy0.data[1])+ny.data[2]*(z-cy0.data[2]));
 189                                                 y1 = (ny.data[0]*(x-cy1.data[0])+ny.data[1]*(y-cy1.data[1])+ny.data[2]*(z-cy1.data[2]));
 190                                                 z0 = (nz.data[0]*(x-cz0.data[0])+nz.data[1]*(y-cz0.data[1])+nz.data[2]*(z-cz0.data[2]));
 191                                                 z1 = (nz.data[0]*(x-cz1.data[0])+nz.data[1]*(y-cz1.data[1])+nz.data[2]*(z-cz1.data[2]));
 192                                                 if(x0*x1 > 0 || y0*y1 > 0 || z0*z1 > 0) { //
 193                                                         //DEBUGPRINT6("x0-z1: %f %f %f %f %f %f\n", x0, x1, y0, y1, z0, z1);
 194                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 195                                                         sum += data;
 196                                                         count++;
 197                                                 }
 198                                         }
 199                                 }
 200                                 }
 201                                 }
 202                                 if(0<count) {
 203                                         avg = sum/count;
 204                                 }
 205                         }
 206                         DEBUGPRINT2("average %f %d\n", avg, count);
 207
 208                         for(z=0; z<out->HeaderN.z; z++) {
 209                         for(y=0; y<out->HeaderN.y; y++) {
 210                         for(x=0; x<out->HeaderN.x; x++) {
 211                                 data = orgData;
 212                 flagIn = 0;
 213                                 if(!linfo->flagRotation) {
 214                                         //DEBUGPRINT("Not Rotated\n")
 215                                         if(     fabs(x - linfo->c.x ) <= (linfo->n.x/2.0)
 216                                          && fabs(y - linfo->c.y ) <= (linfo->n.y/2.0)
 217                                          && fabs(z - linfo->c.z ) <= (linfo->n.z/2.0) ) { // Inside
 218                                                 mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 219                                                 DEBUGPRINT3("(%d, %d, %d) inside square\n", x, y, z);
 220                         flagIn = 1;
 221                                         }
 222                                 } else {
 223                                         x0 = (nx.data[0]*(x-cx0.data[0])+nx.data[1]*(y-cx0.data[1])+nx.data[2]*(z-cx0.data[2]));
 224                                         x1 = (nx.data[0]*(x-cx1.data[0])+nx.data[1]*(y-cx1.data[1])+nx.data[2]*(z-cx1.data[2]));
 225                                         y0 = (ny.data[0]*(x-cy0.data[0])+ny.data[1]*(y-cy0.data[1])+ny.data[2]*(z-cy0.data[2]));
 226                                         y1 = (ny.data[0]*(x-cy1.data[0])+ny.data[1]*(y-cy1.data[1])+ny.data[2]*(z-cy1.data[2]));
 227                                         z0 = (nz.data[0]*(x-cz0.data[0])+nz.data[1]*(y-cz0.data[1])+nz.data[2]*(z-cz0.data[2]));
 228                                         z1 = (nz.data[0]*(x-cz1.data[0])+nz.data[1]*(y-cz1.data[1])+nz.data[2]*(z-cz1.data[2]));
 229                     if(x0*x1 <= 0 && y0*y1 <= 0 && z0*z1 <= 0) {
 230                                             mrcPixelDataGet(in, x, y, z,        &data, mrcPixelRePart, mrcPixelHowNearest);
 231                         flagIn = 1;
 232                     }
 233                                 }
 234                 if(flagIn==1) {
 235                                     if(linfo->mode==0) {
 236                                             data -= avg;
 237                                     } else if(linfo->mode==1) {
 238                                         data = 0.0;
 239                                     } else if(linfo->mode==2) {
 240                         data = 1.0;
 241                     }
 242                 }
 243                                 mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 244                         }
 245                         }
 246                         }
 247                         break;
 248
 249                 case lmrcImageMaskingInfoCylinder:
 250                         sum      = 0.0; count = 0; avg = 0;
 251                         if (linfo->mode==0){
 252                                 for(z=0; z<out->HeaderN.z; z++) {
 253                                 for(y=0; y<out->HeaderN.y; y++) {
 254                                 for(x=0; x<out->HeaderN.x; x++) {
 255                                         if(!linfo->flagRotation) {
 256                                                 if(     SQR((x - linfo->c.x )/(linfo->n.x/2.0))
 257                                                   + SQR((y - linfo->c.y )/(linfo->n.y/2.0)) > 1
 258                                                  || fabs(z - linfo->c.z ) > (linfo->n.z/2.0) ) { // Outside
 259                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 260                                                         sum += data;
 261                                                         count++;
 262                                                 }
 263                                         } else {
 264                                                 tmp.data[0] = x -linfo->c.x;
 265                                                 tmp.data[1] = y -linfo->c.y;
 266                                                 tmp.data[2] = z -linfo->c.z;
 267                                                 matrix3DMultiplyVector(&tmp, matRotAnti);
 268                                                 z0 = (nz.data[0]*(x-cz0.data[0])+nz.data[1]*(y-cz0.data[1])+nz.data[2]*(z-cz0.data[2]));
 269                                                 z1 = (nz.data[0]*(x-cz1.data[0])+nz.data[1]*(y-cz1.data[1])+nz.data[2]*(z-cz1.data[2]));
 270                                                 if(SQR((tmp.data[0])/(linfo->n.x/2.0))
 271                                                  + SQR((tmp.data[1])/(linfo->n.y/2.0))  > 1
 272                                                 || z0*z1 > 0) {                                  // Outside
 273                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 274                                                         sum += data;
 275                                                         count++;
 276                                                 }
 277                                         }
 278                                 }
 279                                 }
 280                                 }
 281                                 if(0<count) {
 282                                         avg = sum/count;
 283                                 }
 284                         }
 285                         DEBUGPRINT2("average %f %d\n", avg, count);
 286
 287                         for(z=0; z<out->HeaderN.z; z++) {
 288                         for(y=0; y<out->HeaderN.y; y++) {
 289                         for(x=0; x<out->HeaderN.x; x++) {
 290                                 data = orgData;
 291                 flagIn = 0;
 292                                 if(!linfo->flagRotation) {
 293                                         //DEBUGPRINT("Not Rotated\n")
 294                                         if(SQR((x - linfo->c.x )/(linfo->n.x/2.0))
 295                                          + SQR((y - linfo->c.y )/(linfo->n.y/2.0)) <= 1
 296                                         && fabs(z - linfo->c.z ) <= (linfo->n.z/2.0) ) { // Inside
 297                                                 mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 298                         flagIn = 1;
 299                                                 DEBUGPRINT3("(%d, %d, %d) inside square\n", x, y, z);
 300                                         }
 301                                 } else {
 302                                         tmp.data[0] = x -linfo->c.x;
 303                                         tmp.data[1] = y -linfo->c.y;
 304                                         tmp.data[2] = z -linfo->c.z;
 305                                         matrix3DMultiplyVector(&tmp, matRotAnti);
 306                                         z0 = (nz.data[0]*(x-cz0.data[0])+nz.data[1]*(y-cz0.data[1])+nz.data[2]*(z-cz0.data[2]));
 307                                         z1 = (nz.data[0]*(x-cz1.data[0])+nz.data[1]*(y-cz1.data[1])+nz.data[2]*(z-cz1.data[2]));
 308                                         if(SQR((tmp.data[0])/(linfo->n.x/2.0))
 309                                          + SQR((tmp.data[1])/(linfo->n.y/2.0)) <= 1
 310                                         && z0*z1 <= 0) {   // Inside
 311                                                 mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 312                         flagIn = 1;
 313                                         }
 314                                 }
 315                 if(flagIn==1) {
 316                                     if(linfo->mode==0) {
 317                                             data -= avg;
 318                                     } else if(linfo->mode==1) {
 319                                         data = 0.0;
 320                                     } else if(linfo->mode==2) {
 321                         data = 1.0;
 322                     }
 323                 }
 324                                 mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 325                         }
 326                         }
 327                         }
 328                         break;
 329
 330                 case lmrcImageMaskingInfoRoundedSquare:
 331                         sum      = 0.0; count = 0; avg = 0;
 332                         if (linfo->mode==0){
 333                                 for(z=0; z<out->HeaderN.z; z++) {
 334                                 for(y=0; y<out->HeaderN.y; y++) {
 335                                 for(x=0; x<out->HeaderN.x; x++) {
 336                                         if(!linfo->flagRotation) {
 337                                                 if(     fabs(x - linfo->c.x ) > (linfo->n.x/2.0)
 338                                                  || fabs(y - linfo->c.y ) > (linfo->n.y/2.0)
 339                                                  || fabs(z - linfo->c.z ) > (linfo->n.z/2.0) ) { // Outside
 340                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 341                                                         sum += data;
 342                                                         count++;
 343                                                 }
 344                                         } else {
 345                                                 x0 = (nx.data[0]*(x-cx0.data[0])+nx.data[1]*(y-cx0.data[1])+nx.data[2]*(z-cx0.data[2]));
 346                                                 x1 = (nx.data[0]*(x-cx1.data[0])+nx.data[1]*(y-cx1.data[1])+nx.data[2]*(z-cx1.data[2]));
 347                                                 y0 = (ny.data[0]*(x-cy0.data[0])+ny.data[1]*(y-cy0.data[1])+ny.data[2]*(z-cy0.data[2]));
 348                                                 y1 = (ny.data[0]*(x-cy1.data[0])+ny.data[1]*(y-cy1.data[1])+ny.data[2]*(z-cy1.data[2]));
 349                                                 z0 = (nz.data[0]*(x-cz0.data[0])+nz.data[1]*(y-cz0.data[1])+nz.data[2]*(z-cz0.data[2]));
 350                                                 z1 = (nz.data[0]*(x-cz1.data[0])+nz.data[1]*(y-cz1.data[1])+nz.data[2]*(z-cz1.data[2]));
 351                                                 if(x0*x1 > 0 || y0*y1 > 0 || z0*z1 > 0) { //
 352                                                         //DEBUGPRINT6("x0-z1: %f %f %f %f %f %f\n", x0, x1, y0, y1, z0, z1);
 353                                                         mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 354                                                         sum += data;
 355                                                         count++;
 356                                                 }
 357                                         }
 358                                 }
 359                                 }
 360                                 }
 361                                 if(0<count) {
 362                                         avg = sum/count;
 363                                 }
 364                         }
 365                         DEBUGPRINT2("average %f %d\n", avg, count);
 366
 367                         for(z=0; z<out->HeaderN.z; z++) {
 368                         for(y=0; y<out->HeaderN.y; y++) {
 369                         for(x=0; x<out->HeaderN.x; x++) {
 370                                 data = orgData;
 371                 flagIn = 0;
 372                                 if(!linfo->flagRotation) {
 373                                         //DEBUGPRINT("Not Rotated\n")
 374                                         if(     fabs(x - linfo->c.x ) <= (linfo->n.x/2.0)
 375                                          && fabs(y - linfo->c.y ) <= (linfo->n.y/2.0)
 376                                          && fabs(z - linfo->c.z ) <= (linfo->n.z/2.0) ) { // Inside
 377                                                 mrcPixelDataGet(in, x, y, z,    &data, mrcPixelRePart, mrcPixelHowNearest);
 378                                                 DEBUGPRINT3("(%d, %d, %d) inside square\n", x, y, z);
 379                         flagIn = 1;
 380                                         }
 381                                 } else {
 382                                         x0 = (nx.data[0]*(x-cx0.data[0])+nx.data[1]*(y-cx0.data[1])+nx.data[2]*(z-cx0.data[2]));
 383                                         x1 = (nx.data[0]*(x-cx1.data[0])+nx.data[1]*(y-cx1.data[1])+nx.data[2]*(z-cx1.data[2]));
 384                                         y0 = (ny.data[0]*(x-cy0.data[0])+ny.data[1]*(y-cy0.data[1])+ny.data[2]*(z-cy0.data[2]));
 385                                         y1 = (ny.data[0]*(x-cy1.data[0])+ny.data[1]*(y-cy1.data[1])+ny.data[2]*(z-cy1.data[2]));
 386                                         z0 = (nz.data[0]*(x-cz0.data[0])+nz.data[1]*(y-cz0.data[1])+nz.data[2]*(z-cz0.data[2]));
 387                                         z1 = (nz.data[0]*(x-cz1.data[0])+nz.data[1]*(y-cz1.data[1])+nz.data[2]*(z-cz1.data[2]));
 388                     if(x0*x1 <= 0 && y0*y1 <= 0 && z0*z1 <= 0) {
 389                                             mrcPixelDataGet(in, x, y, z,        &data, mrcPixelRePart, mrcPixelHowNearest);
 390                         flagIn = 1;
 391                     }
 392                                 }
 393                 if(flagIn==1) {
 394                                     if(linfo->mode==0) {
 395                                             data -= avg;
 396                                     } else if(linfo->mode==1) {
 397                                         data = 0.0;
 398                                     } else if(linfo->mode==2) {
 399                         data = 1.0;
 400                     }
 401                 }
 402                                 mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 403                         }
 404                         }
 405                         }
 406                         break;
 407
 408                 default: {
 409                         fprintf(stderr, "Not supported shape: %d\n", linfo->shape);
 410                         exit(EXIT_FAILURE);
 411                         break;
 412                 }
 413         }
 414 }