src/Tools/mrcImage/mrcImageAbnormalValueRemove/src/mrcImageAbnormalValueRemove.c

   1 /*
   2 # %M% %Y% %I%
   3 # The latest update : %G% at %U%
   4 #
   5 #%Z% mrcImageAbnormalValueRemove ver %I%
   6 #%Z% Created by
   7 #%Z%
   8 #%Z% Usage : mrcImageAbnormalValueRemove
   9 #%Z% Attention
  10 #%Z%
  11 */
  12 static char __sccs_id[] = "%Z%mrcImageAbnormalValueRemove ver%I%; Date:%D% %Z%";
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <string.h>
  16 #include <math.h>
  17 #define GLOBAL_DECLARATION
  18 #include "../inc/config.h"
  19
  20 #define DEBUG
  21 #include "genUtil.h"
  22 #include "Random.h"
  23 #include "mrcImage.h"
  24
  25 typedef struct lmrcImageAbnormalValueRemoveInfo {
  26         float cutValueOverSigma;
  27         int   flagUnsigned;
  28         float UnsignedMax;
  29         float UnsignedHalf;
  30     double Threshold;
  31     mrcImageInformation ImageInfo;
  32     FILE* paramOut;
  33 } lmrcImageAbnormalValueRemoveInfo;
  34
  35 typedef enum lmrcImageAbnormalValueRemoveMode {
  36         lmrcImageAbnormalValueModeDefault = 0,
  37         lmrcImageAbnormalValueModeNear    = 1,
  38         lmrcImageAbnormalValueModeRandom  = 2,
  39         lmrcImageAbnormalValueModeMean    = 3
  40 } lmrcImageAbnormalValueRemoveMode;
  41
  42 typedef enum lmrcImageAbnormalValueRemoveAreaPattern {
  43         lmrcImageAbnormalValueModeAllArea       = 0,
  44         lmrcImageAbnormalValueModeHighValueArea = 1,
  45         lmrcImageAbnormalValueModeLowValueArea  = 2
  46 } lmrcImageAbnormalValueRemoveAreaPattern;
  47
  48 extern void
  49 lmrcImageAbnormalValueRemove(mrcImage* out,
  50                                                          mrcImage* in,
  51                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
  52                                                          lmrcImageAbnormalValueRemoveAreaPattern pattern,
  53                                                          lmrcImageAbnormalValueRemoveMode mode);
  54
  55 extern void
  56 lmrcImageAbnormalValueRemoveAllArea(mrcImage* out,
  57                                                          mrcImage* in,
  58                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
  59                                                          lmrcImageAbnormalValueRemoveMode mode);
  60
  61 extern void
  62 lmrcImageAbnormalValueRemoveHighValueArea(mrcImage* out,
  63                                                          mrcImage* in,
  64                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
  65                                                          lmrcImageAbnormalValueRemoveMode mode);
  66
  67 extern void
  68 lmrcImageAbnormalValueRemoveLowValueArea(mrcImage* out,
  69                                                          mrcImage* in,
  70                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
  71                                                          lmrcImageAbnormalValueRemoveMode mode);
  72
  73 int
  74 main(int argc, char* argv[])
  75 {
  76         mrcImage in;
  77         mrcImage out;
  78         mrcImageAbnormalValueRemoveInfo info;
  79         lmrcImageAbnormalValueRemoveInfo linfo;
  80
  81         init0(&info);
  82     argCheck(&info, argc, argv);
  83     init1(&info);
  84
  85         linfo.cutValueOverSigma = info.cutValue;
  86         linfo.flagUnsigned = info.flagUnsignedMax;
  87         linfo.UnsignedMax = info.UnsignedMax;
  88         linfo.UnsignedHalf= info.UnsignedHalf;
  89         linfo.Threshold = info.Threshold;
  90
  91     linfo.paramOut = info.fptOutParams;
  92
  93         DEBUGPRINT("Program Start\n");
  94         mrcFileRead(&in, info.In, "in main", 0);
  95
  96         out.Header = in.Header;
  97         out.HeaderMode = mrcFloatImage;
  98         mrcInit(&out, NULL);
  99
 100         lmrcImageAbnormalValueRemove(&out, &in, &linfo, info.Pattern, info.mode);
 101
 102         mrcFileWrite(&out, info.Out, "in main", 0);
 103
 104         exit(EXIT_SUCCESS);
 105 }
 106
 107 void
 108 additionalUsage()
 109 {
 110         fprintf(stderr, "----- Additional Usage -----\n");
 111         fprintf(stderr, "----- Mode: \n");
 112         fprintf(stderr, "    %d: Random Value \n", lmrcImageAbnormalValueModeDefault);
 113         fprintf(stderr, "    %d: Near Value \n", lmrcImageAbnormalValueModeNear);
 114         fprintf(stderr, "    %d: Random Value using SD and mean within Normal Area\n", lmrcImageAbnormalValueModeRandom);
 115         fprintf(stderr, "    %d: Mean Value using SD and mean within Normal Area\n", lmrcImageAbnormalValueModeMean);
 116     fprintf(stderr, "----- Area Pattern \n");
 117         fprintf(stderr, "    %d: All area\n", lmrcImageAbnormalValueModeAllArea);
 118         fprintf(stderr, "    %d: HIgh value area\n", lmrcImageAbnormalValueModeHighValueArea);
 119         fprintf(stderr, "    %d: Low  value area\n", lmrcImageAbnormalValueModeLowValueArea);
 120 }
 121
 122 void
 123 lmrcImageAbnormalValueRemove(mrcImage* out,
 124                                                          mrcImage* in,
 125                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
 126                                                          lmrcImageAbnormalValueRemoveAreaPattern pattern,
 127                                                          lmrcImageAbnormalValueRemoveMode mode)
 128 {
 129     mrcImage tmp;
 130     mrcImageParaTypeReal x, y, z;
 131     double data, dst, tmpD;
 132
 133         if(linfo->flagUnsigned){
 134                 tmp.Header = in->Header;
 135                 tmp.HeaderMode = mrcFloatImage;
 136                 mrcInit(&tmp, NULL);
 137                 for(z=0; z<in->HeaderN.z; z++) {
 138                 for(y=0; y<in->HeaderN.y; y++) {
 139                 for(x=0; x<in->HeaderN.x; x++) {
 140                         mrcPixelDataGet(in,   x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 141                         if(linfo->UnsignedHalf <= data) {
 142                                 data -= linfo->UnsignedMax;
 143                         }
 144                 }
 145                 }
 146                 }
 147                 mrcPixelDataSet(&tmp, x, y, z,  data, mrcPixelRePart);
 148                 in = &tmp;
 149         }
 150
 151     switch(pattern) {
 152         case lmrcImageAbnormalValueModeAllArea: {
 153             linfo->ImageInfo.mode = meanOfAll;
 154             lmrcImageInformation(&(linfo->ImageInfo), in);
 155             fprintf(linfo->paramOut, "AbnormalMode: AllArea\n");
 156             fprintf(linfo->paramOut, "ImageInformation(mean,sd):    %g %g\n", linfo->ImageInfo.mean, linfo->ImageInfo.sd);
 157             fprintf(linfo->paramOut, "CutValue:                     %g\n", linfo->cutValueOverSigma);
 158             fprintf(linfo->paramOut, "NormalArea(mean-sd,mean+sd):  %g %g\n", linfo->ImageInfo.mean-linfo->ImageInfo.sd, linfo->ImageInfo.mean+ linfo->ImageInfo.sd);
 159             lmrcImageAbnormalValueRemoveAllArea(out, in, linfo, mode);
 160             break;
 161         }
 162         case lmrcImageAbnormalValueModeHighValueArea: {
 163             linfo->ImageInfo.mode = meanOfHighValueArea;
 164             linfo->ImageInfo.thresHigh = linfo->Threshold;
 165             lmrcImageInformation(&(linfo->ImageInfo), in);
 166             fprintf(linfo->paramOut, "AbnormalMode: HighValueArea\n");
 167             fprintf(linfo->paramOut, "ImageInformation(mean,sd):    %g %g\n", linfo->ImageInfo.meanOfHighValueArea, linfo->ImageInfo.sdOfHighValueArea);
 168             fprintf(linfo->paramOut, "CutValue:                     %g\n", linfo->cutValueOverSigma);
 169             fprintf(linfo->paramOut, "NormalArea(thres,mean+sd):  %g %g\n", linfo->ImageInfo.thresOfHighValueArea,
 170                                                                             linfo->ImageInfo.meanOfHighValueArea+linfo->cutValueOverSigma*linfo->ImageInfo.sdOfHighValueArea);
 171             lmrcImageAbnormalValueRemoveHighValueArea(out, in, linfo, mode);
 172             break;
 173         }
 174         case lmrcImageAbnormalValueModeLowValueArea: {
 175             linfo->ImageInfo.mode = meanOfLowValueArea;
 176             linfo->ImageInfo.thresLow = linfo->Threshold;
 177             lmrcImageInformation(&(linfo->ImageInfo), in);
 178             fprintf(linfo->paramOut, "AbnormalMode: LowValueArea\n");
 179             fprintf(linfo->paramOut, "ImageInformation(mean,sd):    %g %g\n", linfo->ImageInfo.mean, linfo->ImageInfo.sd);
 180             fprintf(linfo->paramOut, "CutValue:                     %g\n", linfo->cutValueOverSigma);
 181             fprintf(linfo->paramOut, "NormalArea(mean-cutValue*sd,thres):  %g %g\n", linfo->ImageInfo.meanOfLowValueArea-linfo->cutValueOverSigma*linfo->ImageInfo.sdOfLowValueArea,
 182                                                                                      linfo->ImageInfo.thresOfLowValueArea);
 183             lmrcImageAbnormalValueRemoveLowValueArea(out, in, linfo, mode);
 184             break;
 185         }
 186         default: {
 187             fprintf(stderr, "Not supported pattern: %d\n", pattern);
 188             break;
 189         }
 190     }
 191 }
 192
 193 void
 194 lmrcImageAbnormalValueRemoveAllArea(mrcImage* out,
 195                                                          mrcImage* in,
 196                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
 197                                                          lmrcImageAbnormalValueRemoveMode mode)
 198 {
 199         mrcImageInformation inInfo;
 200         mrcImage tmp;
 201         mrcImageParaTypeReal x, y, z;
 202         double data, dst, tmpD;
 203
 204     inInfo = linfo->ImageInfo;
 205         switch(mode) {
 206                 case lmrcImageAbnormalValueModeDefault: {
 207                         for(z=0; z<in->HeaderN.z; z++) {
 208                         for(y=0; y<in->HeaderN.y; y++) {
 209                         for(x=0; x<in->HeaderN.x; x++) {
 210                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 211                                 if(data < inInfo.mean - linfo->cutValueOverSigma*inInfo.sd
 212                                   ||      inInfo.mean + linfo->cutValueOverSigma*inInfo.sd < data) {
 213                                         dst = inInfo.mean + inInfo.sd*randomNormalGet(2);
 214                                         switch(out->HeaderMode) {
 215                                                 case mrcCharImage:
 216                                                 case mrcShortImage:
 217                                                 case mrcLongImage: {
 218                                                         if(dst<0) {
 219                                                     DEBUGPRINT1("data is negative: %lf\n", data);
 220                                                                 dst = 0;
 221                                                         }
 222                                                         break;
 223                                                 }
 224                                         }
 225                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 226                                 }
 227                         }
 228                         }
 229                         }
 230                         break;
 231                 }
 232                 case lmrcImageAbnormalValueModeRandom: {
 233                         double sum, mean, SD;
 234                         int n;
 235                         sum = n = 0;
 236                         for(z=0; z<in->HeaderN.z; z++) {
 237                         for(y=0; y<in->HeaderN.y; y++) {
 238                         for(x=0; x<in->HeaderN.x; x++) {
 239                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 240                                 if(inInfo.mean - linfo->cutValueOverSigma*inInfo.sd < data
 241                                  && data < inInfo.mean + linfo->cutValueOverSigma*inInfo.sd ) {
 242                                         sum += data;
 243                                         n++;
 244                                 }
 245                         }
 246                         }
 247                         }
 248                         if(0<n) {
 249                                 mean = sum/n;
 250                         } else {
 251                                 fprintf(stderr, "No data within proper area\n");
 252                         }
 253                         for(z=0; z<in->HeaderN.z; z++) {
 254                         for(y=0; y<in->HeaderN.y; y++) {
 255                         for(x=0; x<in->HeaderN.x; x++) {
 256                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 257                                 if(inInfo.mean - linfo->cutValueOverSigma*inInfo.sd < data
 258                                  && data < inInfo.mean + linfo->cutValueOverSigma*inInfo.sd ) {
 259                                         sum += SQR(data - mean);
 260                                 }
 261                         }
 262                         }
 263                         }
 264                         SD = sqrt(sum/n);
 265
 266                         for(z=0; z<in->HeaderN.z; z++) {
 267                         for(y=0; y<in->HeaderN.y; y++) {
 268                         for(x=0; x<in->HeaderN.x; x++) {
 269                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 270                                 if(data < mean - linfo->cutValueOverSigma*inInfo.sd
 271                                   ||      mean + linfo->cutValueOverSigma*inInfo.sd < data) {
 272                                         dst = mean + SD*randomNormalGet(2);
 273                                         switch(out->HeaderMode) {
 274                                                 case mrcCharImage:
 275                                                 case mrcShortImage:
 276                                                 case mrcLongImage: {
 277                                                         if(dst<0) {
 278                                                     DEBUGPRINT1("data is negative: %lf\n", data);
 279                                                                 dst = 0;
 280                                                         }
 281                                                         break;
 282                                                 }
 283                                         }
 284                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 285                                 }
 286                         }
 287                         }
 288                         }
 289                         break;
 290                 }
 291                 case lmrcImageAbnormalValueModeMean: {
 292                         double sum, mean, SD;
 293                         int n;
 294                         sum = n = 0;
 295                         for(z=0; z<in->HeaderN.z; z++) {
 296                         for(y=0; y<in->HeaderN.y; y++) {
 297                         for(x=0; x<in->HeaderN.x; x++) {
 298                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 299                                 if(inInfo.mean - linfo->cutValueOverSigma*inInfo.sd < data
 300                                  && data < inInfo.mean + linfo->cutValueOverSigma*inInfo.sd ) {
 301                                         sum += data;
 302                                         n++;
 303                                 }
 304                         }
 305                         }
 306                         }
 307                         if(0<n) {
 308                                 mean = sum/n;
 309                         } else {
 310                                 fprintf(stderr, "No data within proper area\n");
 311                         }
 312
 313                         for(z=0; z<in->HeaderN.z; z++) {
 314                         for(y=0; y<in->HeaderN.y; y++) {
 315                         for(x=0; x<in->HeaderN.x; x++) {
 316                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 317                                 if(data < mean - linfo->cutValueOverSigma*inInfo.sd
 318                                   ||      mean + linfo->cutValueOverSigma*inInfo.sd < data) {
 319                                         dst = mean;
 320                                         switch(out->HeaderMode) {
 321                                                 case mrcCharImage:
 322                                                 case mrcShortImage:
 323                                                 case mrcLongImage: {
 324                                                         if(data<0) {
 325                                                     DEBUGPRINT1("data is negative: %lf\n", data);
 326                                                                 dst = 0;
 327                                                         }
 328                                                         break;
 329                                                 }
 330                                         }
 331                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 332                                 }
 333                         }
 334                         }
 335                         }
 336                         break;
 337                 }
 338                 case lmrcImageAbnormalValueModeNear: {
 339                         mrcImageParaTypeReal ox, oy, oz;
 340                         mrcImageParaTypeReal dx, dy, dz;
 341                         mrcImageParaTypeReal min, max;
 342                         mrcImageParaTypeReal minHalf, maxHalf;
 343                         double sum;
 344                         mrcImageParaTypeInteger count;
 345                         mrcImageParaTypeInteger flag;
 346
 347                         min = inInfo.mean - linfo->cutValueOverSigma*inInfo.sd;
 348                         minHalf = inInfo.mean - linfo->cutValueOverSigma*inInfo.sd/2.0;
 349                         max = inInfo.mean + linfo->cutValueOverSigma*inInfo.sd;
 350                         maxHalf = inInfo.mean + linfo->cutValueOverSigma*inInfo.sd/2.0;
 351
 352                         for(z=0; z<in->HeaderN.z; z++) {
 353                         for(y=0; y<in->HeaderN.y; y++) {
 354                         for(x=0; x<in->HeaderN.x; x++) {
 355                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 356
 357                                 if(data < min  || max < data) {
 358                                         DEBUGPRINT6("data %f out of range %f <-> %f at (%f, %f, %f) \n", data, min, max, x, y, z);
 359                                         flag  = 0;
 360                                         dx = dy = dz = 1;
 361                                         while(flag<4) {
 362                                                 count = 0;
 363                                                 sum = 0;
 364                                                 for(oz=z-dz; oz<=z+dz; oz++) {
 365                                                 for(oy=y-dy; oy<=y+dy; oy++) {
 366                                                 for(ox=x-dx; ox<=x+dx; ox++) {
 367                                                         if(SQR((ox-x)/dx)+SQR((oy-y)/dy)+SQR((oz-z)/dz)<=1) {
 368                                                         mrcPixelDataGet(in, ox, oy, oz, &tmpD, mrcPixelRePart, mrcPixelHowNearest);
 369                                                         if( min < tmpD && tmpD < max) {
 370                                                                         sum += tmpD;
 371                                                                 count++;
 372                                                                 }
 373                                                         }
 374                                                 }
 375                                                 }
 376                                                 }
 377                                                 if(4<=count) {
 378                                                         flag ++;
 379                                                         dst = sum/count;
 380                                                         DEBUGPRINT5("data %f count %d at (%f, %f, %f) \n", dst, count, x, y, z);
 381                                                         if(dst < minHalf || maxHalf < dst ) {
 382                                                                 DEBUGPRINT("wider range +1 because of wide abnormal value\n");
 383                                                                 dx++;
 384                                                                 dy++;
 385                                                                 dz++;
 386                                                         } else {
 387                                                                 flag = 4;
 388                                                         }
 389                                                 } else {
 390                                                         DEBUGPRINT("wider range +1\n");
 391                                                         flag=0;
 392                                                         dx++;
 393                                                         dy++;
 394                                                         dz++;
 395                                                 }
 396                                         }
 397                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 398                                 }
 399                         }
 400                         }
 401                         }
 402                         break;
 403                 }
 404                 default: {
 405                         fprintf(stderr, "Not supported mode: %d\n", mode);
 406                         exit(EXIT_FAILURE);
 407                         break;
 408                 }
 409         }
 410         mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 411         return;
 412 }
 413
 414
 415 void
 416 lmrcImageAbnormalValueRemoveHighValueArea(mrcImage* out,
 417                                                          mrcImage* in,
 418                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
 419                                                          lmrcImageAbnormalValueRemoveMode mode)
 420 {
 421         mrcImageInformation inInfo;
 422         mrcImage tmp;
 423         mrcImageParaTypeReal x, y, z;
 424         double data, dst, tmpD;
 425
 426     inInfo = linfo->ImageInfo;
 427         switch(mode) {
 428                 case lmrcImageAbnormalValueModeDefault: {
 429                         for(z=0; z<in->HeaderN.z; z++) {
 430                         for(y=0; y<in->HeaderN.y; y++) {
 431                         for(x=0; x<in->HeaderN.x; x++) {
 432                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 433                                 if(inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea < data) {
 434                                         dst = inInfo.meanOfHighValueArea + inInfo.sdOfHighValueArea*randomNormalGet(2);
 435                                         switch(out->HeaderMode) {
 436                                                 case mrcCharImage:
 437                                                 case mrcShortImage:
 438                                                 case mrcLongImage: {
 439                                                         if(dst<0) {
 440                                                                 dst= 0;
 441                                                         }
 442                                                         break;
 443                                                 }
 444                                         }
 445                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 446                                 }
 447                         }
 448                         }
 449                         }
 450                         break;
 451                 }
 452                 case lmrcImageAbnormalValueModeMean: {
 453             double sum, mean;
 454             int n;
 455             sum = 0;
 456             n = 0;
 457                         for(z=0; z<in->HeaderN.z; z++) {
 458                         for(y=0; y<in->HeaderN.y; y++) {
 459                         for(x=0; x<in->HeaderN.x; x++) {
 460                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 461                                 if( inInfo.thresOfHighValueArea < data &&
 462                     data < inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea ) {
 463                                         sum += data;
 464                                 }
 465             }
 466             }
 467             }
 468             mean = sum/n;
 469                         for(z=0; z<in->HeaderN.z; z++) {
 470                         for(y=0; y<in->HeaderN.y; y++) {
 471                         for(x=0; x<in->HeaderN.x; x++) {
 472                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 473                                 if(inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea < data) {
 474                                         dst = mean;
 475                                         switch(out->HeaderMode) {
 476                                                 case mrcCharImage:
 477                                                 case mrcShortImage:
 478                                                 case mrcLongImage: {
 479                                                         if(dst <0) {
 480                                                                 dst = 0;
 481                                                         }
 482                                                         break;
 483                                                 }
 484                                         }
 485                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 486                                 }
 487             }
 488             }
 489             }
 490             break;
 491         }
 492                 case lmrcImageAbnormalValueModeNear: {
 493                         mrcImageParaTypeReal ox, oy, oz;
 494                         mrcImageParaTypeReal dx, dy, dz;
 495                         mrcImageParaTypeReal min, max;
 496                         mrcImageParaTypeReal minHalf, maxHalf;
 497                         double sum;
 498                         mrcImageParaTypeInteger count;
 499                         mrcImageParaTypeInteger flag;
 500
 501                         min = inInfo.thresOfHighValueArea;
 502                         minHalf = (inInfo.meanOfHighValueArea + min)/2.0;
 503                         max = inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea;
 504                         maxHalf = inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea/2.0;
 505
 506                         for(z=0; z<in->HeaderN.z; z++) {
 507                         for(y=0; y<in->HeaderN.y; y++) {
 508                         for(x=0; x<in->HeaderN.x; x++) {
 509                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 510                                 if(max < data) {
 511                                         DEBUGPRINT6("data %f out of range %f <-> %f at (%f, %f, %f) \n", data, min, max, x, y, z);
 512                                         flag  = 0;
 513                                         dx = dy = dz = 1;
 514                                         while(flag<4) {
 515                                                 count = 0;
 516                                                 sum = 0;
 517                                                 for(oz=z-dz; oz<=z+dz; oz++) {
 518                                                 for(oy=y-dy; oy<=y+dy; oy++) {
 519                                                 for(ox=x-dx; ox<=x+dx; ox++) {
 520                                                         if(SQR((ox-x)/dx)+SQR((oy-y)/dy)+SQR((oz-z)/dz)<=1) {
 521                                                         mrcPixelDataGet(in, ox, oy, oz, &tmpD, mrcPixelRePart, mrcPixelHowNearest);
 522                                                         if(tmpD < max) {
 523                                                                         sum += tmpD;
 524                                                                 count++;
 525                                                                 }
 526                                                         }
 527                                                 }
 528                                                 }
 529                                                 }
 530                                                 if(4<=count) {
 531                                                         flag ++;
 532                                                         dst = sum/count;
 533                                                         DEBUGPRINT5("data %f count %d at (%f, %f, %f) \n", data, count, x, y, z);
 534                                                         if( maxHalf < data ) {
 535                                                                 DEBUGPRINT("wider range +1 because of wide abnormal value\n");
 536                                                                 dx++;
 537                                                                 dy++;
 538                                                                 dz++;
 539                                                         } else {
 540                                                                 flag = 4;
 541                                                         }
 542                                                 } else {
 543                                                         DEBUGPRINT("wider range +1\n");
 544                                                         flag=0;
 545                                                         dx++;
 546                                                         dy++;
 547                                                         dz++;
 548                                                 }
 549                                         }
 550                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 551                                 }
 552                         }
 553                         }
 554                         }
 555             break;
 556         }
 557                 case lmrcImageAbnormalValueModeRandom: {
 558             double sum, mean, SD;
 559             int n;
 560             n=0;
 561             sum = 0;
 562                         for(z=0; z<in->HeaderN.z; z++) {
 563                         for(y=0; y<in->HeaderN.y; y++) {
 564                         for(x=0; x<in->HeaderN.x; x++) {
 565                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 566                                 if(inInfo.thresOfHighValueArea < data
 567                  &&data < inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea) {
 568                     sum += data;
 569                     n++;
 570                                 }
 571             }
 572             }
 573             }
 574             mean = sum/n;
 575             sum = 0;
 576                         for(z=0; z<in->HeaderN.z; z++) {
 577                         for(y=0; y<in->HeaderN.y; y++) {
 578                         for(x=0; x<in->HeaderN.x; x++) {
 579                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 580                                 if(inInfo.thresOfHighValueArea < data
 581                  &&data < inInfo.meanOfHighValueArea + linfo->cutValueOverSigma*inInfo.sdOfHighValueArea) {
 582                     sum += SQR(data-mean);
 583                                 }
 584             }
 585             }
 586             }
 587             SD = sqrt(sum/n);
 588
 589                         for(z=0; z<in->HeaderN.z; z++) {
 590                         for(y=0; y<in->HeaderN.y; y++) {
 591                         for(x=0; x<in->HeaderN.x; x++) {
 592                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 593                                 if(mean + linfo->cutValueOverSigma*SD < data) {
 594                                         dst = mean + SD*randomNormalGet(2);
 595                                         switch(out->HeaderMode) {
 596                                                 case mrcCharImage:
 597                                                 case mrcShortImage:
 598                                                 case mrcLongImage: {
 599                                                         if(dst<0) {
 600                                                     DEBUGPRINT1("data is negative: %lf\n", data);
 601                                                                 dst = 0;
 602                                                         }
 603                                                         break;
 604                                                 }
 605                                         }
 606                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 607                                 }
 608                         }
 609                         }
 610                         }
 611                         break;
 612                 }
 613                 default: {
 614                         fprintf(stderr, "Not supported mode: %d\n", mode);
 615                         exit(EXIT_FAILURE);
 616                         break;
 617                 }
 618         }
 619         mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 620         return;
 621 }
 622
 623
 624 void
 625 lmrcImageAbnormalValueRemoveLowValueArea(mrcImage* out,
 626                                                          mrcImage* in,
 627                                                          lmrcImageAbnormalValueRemoveInfo* linfo,
 628                                                          lmrcImageAbnormalValueRemoveMode mode)
 629 {
 630         mrcImageInformation inInfo;
 631         mrcImage tmp;
 632         mrcImageParaTypeReal x, y, z;
 633         double data, dst, tmpD;
 634
 635     inInfo = linfo->ImageInfo;
 636         switch(mode) {
 637                 case lmrcImageAbnormalValueModeDefault: {
 638                         for(z=0; z<in->HeaderN.z; z++) {
 639                         for(y=0; y<in->HeaderN.y; y++) {
 640                         for(x=0; x<in->HeaderN.x; x++) {
 641                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 642                                 if(data < inInfo.meanOfLowValueArea - linfo->cutValueOverSigma*inInfo.sdOfLowValueArea) {
 643                                         dst = inInfo.meanOfLowValueArea + inInfo.sdOfLowValueArea*randomNormalGet(2);
 644                                         switch(out->HeaderMode) {
 645                                                 case mrcCharImage:
 646                                                 case mrcShortImage:
 647                                                 case mrcLongImage: {
 648                                                         if(dst<0) {
 649                                                                 dst = 0;
 650                                                         }
 651                                                         break;
 652                                                 }
 653                                         }
 654                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 655                                 }
 656                         }
 657                         }
 658                         }
 659                         break;
 660                 }
 661                 case lmrcImageAbnormalValueModeMean: {
 662             double sum, mean;
 663             int n;
 664             n = 0;
 665             sum = 0;
 666                         for(z=0; z<in->HeaderN.z; z++) {
 667                         for(y=0; y<in->HeaderN.y; y++) {
 668                         for(x=0; x<in->HeaderN.x; x++) {
 669                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 670                                 if(inInfo.meanOfLowValueArea - linfo->cutValueOverSigma*inInfo.sdOfLowValueArea < data
 671                    && data < inInfo.thresOfLowValueArea ) {
 672                     sum += data;
 673                     n++;
 674                                 }
 675                         }
 676                         }
 677                         }
 678             mean = sum/n;
 679                         for(z=0; z<in->HeaderN.z; z++) {
 680                         for(y=0; y<in->HeaderN.y; y++) {
 681                         for(x=0; x<in->HeaderN.x; x++) {
 682                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 683                                 if(data < inInfo.meanOfLowValueArea - linfo->cutValueOverSigma*inInfo.sdOfLowValueArea) {
 684                                         dst = mean;
 685                                         switch(out->HeaderMode) {
 686                                                 case mrcCharImage:
 687                                                 case mrcShortImage:
 688                                                 case mrcLongImage: {
 689                                                         if(dst<0) {
 690                                                                 dst = 0;
 691                                                         }
 692                                                         break;
 693                                                 }
 694                                         }
 695                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 696                                 }
 697                         }
 698                         }
 699                         }
 700             break;
 701         }
 702                 case lmrcImageAbnormalValueModeNear: {
 703                         mrcImageParaTypeReal ox, oy, oz;
 704                         mrcImageParaTypeReal dx, dy, dz;
 705                         mrcImageParaTypeReal min, max;
 706                         mrcImageParaTypeReal minHalf, maxHalf;
 707                         double sum;
 708                         mrcImageParaTypeInteger count;
 709                         mrcImageParaTypeInteger flag;
 710
 711                         min = inInfo.meanOfLowValueArea - linfo->cutValueOverSigma*inInfo.sdOfLowValueArea;
 712                         minHalf = (inInfo.meanOfLowValueArea + min)/2.0;
 713                         max = inInfo.thresOfLowValueArea;
 714                         maxHalf = (inInfo.meanOfLowValueArea + max)/2.0;
 715
 716                         for(z=0; z<in->HeaderN.z; z++) {
 717                         for(y=0; y<in->HeaderN.y; y++) {
 718                         for(x=0; x<in->HeaderN.x; x++) {
 719                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 720                                 if(data < min) {
 721                                         DEBUGPRINT6("data %f out of range %f <-> %f at (%f, %f, %f) \n", data, min, max, x, y, z);
 722                                         flag  = 0;
 723                                         dx = dy = dz = 1;
 724                                         while(flag<4) {
 725                                                 count = 0;
 726                                                 sum = 0;
 727                                                 for(oz=z-dz; oz<=z+dz; oz++) {
 728                                                 for(oy=y-dy; oy<=y+dy; oy++) {
 729                                                 for(ox=x-dx; ox<=x+dx; ox++) {
 730                                                         if(SQR((ox-x)/dx)+SQR((oy-y)/dy)+SQR((oz-z)/dz)<=1) {
 731                                                         mrcPixelDataGet(in, ox, oy, oz, &data, mrcPixelRePart, mrcPixelHowNearest);
 732                                                         if( min < data && data < max) {
 733                                                                         sum += data;
 734                                                                 count++;
 735                                                                 }
 736                                                         }
 737                                                 }
 738                                                 }
 739                                                 }
 740                                                 if(4<=count) {
 741                                                         flag ++;
 742                                                         data = sum/count;
 743                                                         DEBUGPRINT5("data %f count %d at (%f, %f, %f) \n", data, count, x, y, z);
 744                                                         if(data<minHalf) {
 745                                                                 DEBUGPRINT("wider range +1 because of wide abnormal value\n");
 746                                                                 dx++;
 747                                                                 dy++;
 748                                                                 dz++;
 749                                                         } else {
 750                                                                 flag = 4;
 751                                                         }
 752
 753                                                 } else {
 754                                                         DEBUGPRINT("wider range +1\n");
 755                                                         flag=0;
 756                                                         dx++;
 757                                                         dy++;
 758                                                         dz++;
 759                                                 }
 760                                         }
 761                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 762                                 }
 763                         }
 764                         }
 765                         }
 766             break;
 767         }
 768                 case lmrcImageAbnormalValueModeRandom: {
 769             double sum, mean, SD;
 770             int n;
 771             n=0;
 772             sum = 0;
 773                         for(z=0; z<in->HeaderN.z; z++) {
 774                         for(y=0; y<in->HeaderN.y; y++) {
 775                         for(x=0; x<in->HeaderN.x; x++) {
 776                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 777                                 if(data < inInfo.thresOfLowValueArea
 778                  &&inInfo.meanOfLowValueArea - linfo->cutValueOverSigma*inInfo.sdOfLowValueArea < data) {
 779                     sum += data;
 780                     n++;
 781                                 }
 782             }
 783             }
 784             }
 785             mean = sum/n;
 786             sum = 0;
 787                         for(z=0; z<in->HeaderN.z; z++) {
 788                         for(y=0; y<in->HeaderN.y; y++) {
 789                         for(x=0; x<in->HeaderN.x; x++) {
 790                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 791                                 if(data < inInfo.thresOfLowValueArea
 792                  &&inInfo.meanOfLowValueArea - linfo->cutValueOverSigma*inInfo.sdOfLowValueArea  <data) {
 793                     sum += SQR(data-mean);
 794                                 }
 795             }
 796             }
 797             }
 798             SD = sqrt(sum/n);
 799
 800                         for(z=0; z<in->HeaderN.z; z++) {
 801                         for(y=0; y<in->HeaderN.y; y++) {
 802                         for(x=0; x<in->HeaderN.x; x++) {
 803                                 mrcPixelDataGet(in, x, y, z, &data, mrcPixelRePart, mrcPixelHowNearest);
 804                                 if(data < mean - linfo->cutValueOverSigma*SD) {
 805                                         dst = mean + SD*randomNormalGet(2);
 806                                         switch(out->HeaderMode) {
 807                                                 case mrcCharImage:
 808                                                 case mrcShortImage:
 809                                                 case mrcLongImage: {
 810                                                         if(dst<0) {
 811                                                     DEBUGPRINT1("data is negative: %lf\n", data);
 812                                                                 dst = 0;
 813                                                         }
 814                                                         break;
 815                                                 }
 816                                         }
 817                     fprintf(linfo->paramOut,"AbnormalValue: %f -> %f at %f %f %f\n", data, dst, x, y, z);
 818                                 }
 819                         }
 820                         }
 821                         }
 822             break;
 823         }
 824                 default: {
 825                         fprintf(stderr, "Not supported mode: %d\n", mode);
 826                         exit(EXIT_FAILURE);
 827                         break;
 828                 }
 829         }
 830         mrcPixelDataSet(out, x, y, z, data, mrcPixelRePart);
 831         return;
 832 }
 833
 834