00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "processor.h"
00037 #include "sparx/processor_sparx.h"
00038 #include "ctf.h"
00039 #include "xydata.h"
00040 #include "emdata.h"
00041 #include "emassert.h"
00042 #include "randnum.h"
00043 #include "symmetry.h"
00044
00045 #include <gsl/gsl_randist.h>
00046 #include <gsl/gsl_statistics.h>
00047 #include <gsl/gsl_wavelet.h>
00048 #include <gsl/gsl_wavelet2d.h>
00049 #include <algorithm>
00050 #include <ctime>
00051
00052 #ifdef EMAN2_USING_CUDA
00053 #include "cuda/cuda_util.h"
00054 #include "cuda/cuda_processor.h"
00055 #endif // EMAN2_USING_CUDA
00056
00057 using namespace EMAN;
00058 using std::reverse;
00059
00060 template <> Factory < Processor >::Factory()
00061 {
00062 force_add(&LowpassSharpCutoffProcessor::NEW);
00063 force_add(&HighpassSharpCutoffProcessor::NEW);
00064 force_add(&LowpassGaussProcessor::NEW);
00065 force_add(&HighpassGaussProcessor::NEW);
00066 force_add(&HighpassAutoPeakProcessor::NEW);
00067 force_add(&LinearRampFourierProcessor::NEW);
00068
00069 force_add(&LowpassTanhProcessor::NEW);
00070 force_add(&HighpassTanhProcessor::NEW);
00071 force_add(&HighpassButterworthProcessor::NEW);
00072 force_add(&AmpweightFourierProcessor::NEW);
00073 force_add(&Wiener2DFourierProcessor::NEW);
00074
00075 force_add(&LinearPyramidProcessor::NEW);
00076 force_add(&LinearRampProcessor::NEW);
00077 force_add(&AbsoluateValueProcessor::NEW);
00078 force_add(&BooleanProcessor::NEW);
00079 force_add(&ValuePowProcessor::NEW);
00080 force_add(&ValueSquaredProcessor::NEW);
00081 force_add(&ValueSqrtProcessor::NEW);
00082 force_add(&Rotate180Processor::NEW);
00083 force_add(&TransformProcessor::NEW);
00084 force_add(&ScaleTransformProcessor::NEW);
00085 force_add(&IntTranslateProcessor::NEW);
00086 force_add(&InvertCarefullyProcessor::NEW);
00087
00088 force_add(&ClampingProcessor::NEW);
00089 force_add(&NSigmaClampingProcessor::NEW);
00090
00091 force_add(&ToZeroProcessor::NEW);
00092 force_add(&ToMinvalProcessor::NEW);
00093 force_add(&CutToZeroProcessor::NEW);
00094 force_add(&BinarizeProcessor::NEW);
00095 force_add(&BinarizeFourierProcessor::NEW);
00096 force_add(&CollapseProcessor::NEW);
00097 force_add(&LinearXformProcessor::NEW);
00098
00099 force_add(&ExpProcessor::NEW);
00100 force_add(&RangeThresholdProcessor::NEW);
00101 force_add(&SigmaProcessor::NEW);
00102 force_add(&LogProcessor::NEW);
00103 force_add(&FiniteProcessor::NEW);
00104
00105 force_add(&BinaryOperateProcessor<MaxPixelOperator>::NEW);
00106 force_add(&BinaryOperateProcessor<MinPixelOperator>::NEW);
00107
00108 force_add(&PaintProcessor::NEW);
00109 force_add(&WatershedProcessor::NEW);
00110 force_add(&MaskSharpProcessor::NEW);
00111 force_add(&MaskEdgeMeanProcessor::NEW);
00112 force_add(&MaskNoiseProcessor::NEW);
00113 force_add(&MaskGaussProcessor::NEW);
00114 force_add(&MaskGaussNonuniformProcessor::NEW);
00115 force_add(&MaskGaussInvProcessor::NEW);
00116
00117 force_add(&MaxShrinkProcessor::NEW);
00118 force_add(&MinShrinkProcessor::NEW);
00119 force_add(&MeanShrinkProcessor::NEW);
00120 force_add(&MedianShrinkProcessor::NEW);
00121 force_add(&FFTResampleProcessor::NEW);
00122
00123 force_add(&MakeRadiusSquaredProcessor::NEW);
00124 force_add(&MakeRadiusProcessor::NEW);
00125
00126 force_add(&ComplexNormPixel::NEW);
00127
00128 force_add(&LaplacianProcessor::NEW);
00129 force_add(&ZeroConstantProcessor::NEW);
00130
00131 force_add(&BoxMedianProcessor::NEW);
00132 force_add(&BoxSigmaProcessor::NEW);
00133 force_add(&BoxMaxProcessor::NEW);
00134
00135 force_add(&MinusPeakProcessor::NEW);
00136 force_add(&PeakOnlyProcessor::NEW);
00137 force_add(&DiffBlockProcessor::NEW);
00138
00139 force_add(&CutoffBlockProcessor::NEW);
00140 force_add(&GradientRemoverProcessor::NEW);
00141 force_add(&GradientPlaneRemoverProcessor::NEW);
00142 force_add(&FlattenBackgroundProcessor::NEW);
00143 force_add(&VerticalStripeProcessor::NEW);
00144 force_add(&RealToFFTProcessor::NEW);
00145 force_add(&SigmaZeroEdgeProcessor::NEW);
00146 force_add(&RampProcessor::NEW);
00147
00148 force_add(&BeamstopProcessor::NEW);
00149 force_add(&MeanZeroEdgeProcessor::NEW);
00150 force_add(&AverageXProcessor::NEW);
00151 force_add(&DecayEdgeProcessor::NEW);
00152 force_add(&ZeroEdgeRowProcessor::NEW);
00153 force_add(&ZeroEdgePlaneProcessor::NEW);
00154
00155 force_add(&BilateralProcessor::NEW);
00156
00157 force_add(&ConvolutionProcessor::NEW);
00158
00159 force_add(&NormalizeStdProcessor::NEW);
00160 force_add(&NormalizeUnitProcessor::NEW);
00161 force_add(&NormalizeUnitSumProcessor::NEW);
00162 force_add(&NormalizeMaskProcessor::NEW);
00163 force_add(&NormalizeEdgeMeanProcessor::NEW);
00164 force_add(&NormalizeCircleMeanProcessor::NEW);
00165 force_add(&NormalizeLREdgeMeanProcessor::NEW);
00166 force_add(&NormalizeMaxMinProcessor::NEW);
00167 force_add(&NormalizeByMassProcessor::NEW);
00168 force_add(&NormalizeRowProcessor::NEW);
00169 force_add(&NormalizeRampNormVar::NEW);
00170
00171 force_add(&HistogramBin::NEW);
00172
00173 force_add(&NormalizeToLeastSquareProcessor::NEW);
00174
00175 force_add(&RotationalAverageProcessor::NEW);
00176 force_add(&RotationalSubstractProcessor::NEW);
00177 force_add(&FlipProcessor::NEW);
00178 force_add(&TransposeProcessor::NEW);
00179 force_add(&MirrorProcessor::NEW);
00180
00181 force_add(&AddNoiseProcessor::NEW);
00182 force_add(&AddSigmaNoiseProcessor::NEW);
00183 force_add(&AddRandomNoiseProcessor::NEW);
00184
00185 force_add(&PhaseToCenterProcessor::NEW);
00186 force_add(&PhaseToCornerProcessor::NEW);
00187 force_add(&FourierToCenterProcessor::NEW);
00188 force_add(&FourierToCornerProcessor::NEW);
00189 force_add(&AutoMask2DProcessor::NEW);
00190 force_add(&AutoMask3DProcessor::NEW);
00191 force_add(&AutoMask3D2Processor::NEW);
00192 force_add(&AddMaskShellProcessor::NEW);
00193 force_add(&AutoMaskAsymUnit::NEW);
00194
00195 force_add(&CTFSNRWeightProcessor::NEW);
00196
00197 force_add(&ToMassCenterProcessor::NEW);
00198 force_add(&PhaseToMassCenterProcessor::NEW);
00199 force_add(&ACFCenterProcessor::NEW);
00200 force_add(&SNRProcessor::NEW);
00201
00202 force_add(&XGradientProcessor::NEW);
00203 force_add(&YGradientProcessor::NEW);
00204 force_add(&ZGradientProcessor::NEW);
00205
00206 force_add(&FileFourierProcessor::NEW);
00207
00208 force_add(&SymSearchProcessor::NEW);
00209 force_add(&LocalNormProcessor::NEW);
00210
00211 force_add(&IndexMaskFileProcessor::NEW);
00212 force_add(&CoordinateMaskFileProcessor::NEW);
00213 force_add(&SetSFProcessor::NEW);
00214 force_add(&MatchSFProcessor::NEW);
00215
00216 force_add(&SmartMaskProcessor::NEW);
00217 force_add(&IterBinMaskProcessor::NEW);
00218
00219 force_add(&TestImageGaussian::NEW);
00220 force_add(&TestImagePureGaussian::NEW);
00221 force_add(&TestImageSinewave::NEW);
00222 force_add(&TestImageSphericalWave::NEW);
00223 force_add(&TestImageSinewaveCircular::NEW);
00224 force_add(&TestImageSquarecube::NEW);
00225 force_add(&TestImageCirclesphere::NEW);
00226 force_add(&TestImageAxes::NEW);
00227 force_add(&TestImageNoiseUniformRand::NEW);
00228 force_add(&TestImageNoiseGauss::NEW);
00229 force_add(&TestImageScurve::NEW);
00230 force_add(&TestImageCylinder::NEW);
00231 force_add(&TestImageGradient::NEW);
00232 force_add(&TestTomoImage::NEW);
00233 force_add(&TestImageLineWave::NEW);
00234 force_add(&TestImageEllipse::NEW);
00235 force_add(&TestImageHollowEllipse::NEW);
00236 force_add(&TestImageFourierNoiseGaussian::NEW);
00237 force_add(&TestImageFourierNoiseProfile::NEW);
00238
00239 force_add(&TomoTiltEdgeMaskProcessor::NEW);
00240 force_add(&TomoTiltAngleWeightProcessor::NEW);
00241
00242 force_add(&NewLowpassTopHatProcessor::NEW);
00243 force_add(&NewHighpassTopHatProcessor::NEW);
00244 force_add(&NewBandpassTopHatProcessor::NEW);
00245 force_add(&NewHomomorphicTopHatProcessor::NEW);
00246 force_add(&NewLowpassGaussProcessor::NEW);
00247 force_add(&NewHighpassGaussProcessor::NEW);
00248 force_add(&NewBandpassGaussProcessor::NEW);
00249 force_add(&NewHomomorphicGaussProcessor::NEW);
00250 force_add(&NewInverseGaussProcessor::NEW);
00251 force_add(&NewLowpassButterworthProcessor::NEW);
00252 force_add(&NewHighpassButterworthProcessor::NEW);
00253 force_add(&NewHomomorphicButterworthProcessor::NEW);
00254 force_add(&NewLowpassTanhProcessor::NEW);
00255 force_add(&NewHighpassTanhProcessor::NEW);
00256 force_add(&NewBandpassTanhProcessor::NEW);
00257 force_add(&NewHomomorphicTanhProcessor::NEW);
00258 force_add(&NewRadialTableProcessor::NEW);
00259 force_add(&InverseKaiserI0Processor::NEW);
00260 force_add(&InverseKaiserSinhProcessor::NEW);
00261 force_add(&CCDNormProcessor::NEW);
00262 force_add(&CTF_Processor::NEW);
00263 force_add(&SHIFTProcessor::NEW);
00264
00265 force_add(&WaveletProcessor::NEW);
00266 force_add(&FFTProcessor::NEW);
00267 force_add(&RadialProcessor::NEW);
00268
00269 force_add(&DirectionalSumProcessor::NEW);
00270
00271
00272 force_add(&ModelEMCylinderProcessor::NEW);
00273 force_add(&ApplyPolynomialProfileToHelix::NEW);
00274 force_add(&BinarySkeletonizerProcessor::NEW);
00275
00276 #ifdef EMAN2_USING_CUDA
00277 force_add(&CudaMultProcessor::NEW);
00278 force_add(&CudaCorrelationProcessor::NEW);
00279 #endif // EMAN2_USING_CUDA
00280 }
00281
00282 void FiniteProcessor::process_pixel(float *x) const
00283 {
00284 if ( !Util::goodf(x) ) {
00285 *x = to;
00286 }
00287 }
00288
00289
00290 EMData* Processor::process(const EMData * const image)
00291 {
00292 EMData * result = image->copy();
00293 process_inplace(result);
00294 return result;
00295 }
00296
00297 void ImageProcessor::process_inplace(EMData * image)
00298 {
00299 if (!image) {
00300 LOGWARN("NULL image");
00301 return;
00302 }
00303
00304 EMData *processor_image = create_processor_image();
00305
00306 if (image->is_complex()) {
00307 (*image) *= *processor_image;
00308 }
00309 else {
00310 EMData *fft = image->do_fft();
00311 (*fft) *= (*processor_image);
00312 EMData *ift = fft->do_ift();
00313
00314 float *data = image->get_data();
00315 float *t = data;
00316 float *ift_data = ift->get_data();
00317
00318 data = ift_data;
00319 ift_data = t;
00320
00321 ift->update();
00322
00323 if( fft )
00324 {
00325 delete fft;
00326 fft = 0;
00327 }
00328
00329 if( ift )
00330 {
00331 delete ift;
00332 ift = 0;
00333 }
00334 }
00335
00336 image->update();
00337 }
00338
00339 #define FFTRADIALOVERSAMPLE 4
00340 void FourierProcessor::process_inplace(EMData * image)
00341 {
00342 if (!image) {
00343 LOGWARN("NULL Image");
00344 return;
00345 }
00346
00347 preprocess(image);
00348
00349 int array_size = FFTRADIALOVERSAMPLE * image->get_ysize();
00350 float step=0.5f/array_size;
00351
00352 vector < float >yarray(array_size);
00353
00354 create_radial_func(yarray);
00355
00356 if (image->is_complex()) {
00357 image->apply_radial_func(0, step, yarray);
00358 }
00359 else {
00360 EMData *fft = image->do_fft();
00361 fft->apply_radial_func(0, step, yarray);
00362 EMData *ift = fft->do_ift();
00363
00364 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));
00365
00366
00367
00368 if( fft )
00369 {
00370 delete fft;
00371 fft = 0;
00372 }
00373
00374 if( ift )
00375 {
00376 delete ift;
00377 ift = 0;
00378 }
00379 }
00380
00381 image->update();
00382 }
00383
00384 void FourierAnlProcessor::process_inplace(EMData * image)
00385 {
00386 if (!image) {
00387 LOGWARN("NULL Image");
00388 return;
00389 }
00390
00391 preprocess(image);
00392
00393
00394
00395
00396
00397
00398
00399 if (image->is_complex()) {
00400 vector <float>yarray = image->calc_radial_dist(image->get_ysize()/2,0,1.0,1);
00401 create_radial_func(yarray,image);
00402 image->apply_radial_func(0, 0.5f/yarray.size(), yarray);
00403 }
00404 else {
00405 EMData *fft = image->do_fft();
00406 vector <float>yarray = fft->calc_radial_dist(fft->get_ysize()/2,0,1.0,1);
00407 create_radial_func(yarray,image);
00408 fft->apply_radial_func(0, 0.5f/yarray.size(), yarray);
00409 EMData *ift = fft->do_ift();
00410
00411 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));
00412
00413
00414
00415 delete fft;
00416 delete ift;
00417
00418 }
00419
00420 image->update();
00421 }
00422
00423 void LowpassFourierProcessor::preprocess(EMData * image)
00424 {
00425 if(params.has_key("apix")) {
00426 image->set_attr("apix_x", (float)params["apix"]);
00427 image->set_attr("apix_y", (float)params["apix"]);
00428 image->set_attr("apix_z", (float)params["apix"]);
00429 }
00430
00431 const Dict dict = image->get_attr_dict();
00432
00433 if( params.has_key("cutoff_abs") ) {
00434 lowpass = params["cutoff_abs"];
00435 }
00436 else if( params.has_key("cutoff_freq") ) {
00437 lowpass = (float)params["cutoff_freq"] * (float)dict["apix_x"] * (float)dict["nx"] / 2.0f;
00438 }
00439 else if( params.has_key("cutoff_pixels") ) {
00440 lowpass = (float)params["cutoff_pixels"] / (float)dict["nx"];
00441 }
00442 }
00443
00444 void HighpassFourierProcessor::preprocess(EMData * image)
00445 {
00446 if(params.has_key("apix")) {
00447 image->set_attr("apix_x", (float)params["apix"]);
00448 image->set_attr("apix_y", (float)params["apix"]);
00449 image->set_attr("apix_z", (float)params["apix"]);
00450 }
00451
00452 const Dict dict = image->get_attr_dict();
00453
00454 if( params.has_key("cutoff_abs") ) {
00455 highpass = params["cutoff_abs"];
00456 }
00457 else if( params.has_key("cutoff_freq") ) {
00458 highpass = (float)params["cutoff_freq"] * (float)dict["apix_x"] * (float)dict["nx"] / 2.0f;
00459 }
00460 else if( params.has_key("cutoff_pixels") ) {
00461 highpass = (float)params["cutoff_pixels"] / (float)dict["nx"];
00462 }
00463 }
00464
00465 void SNREvalProcessor::process_inplace(EMData * image)
00466 {
00467 int ys=image->get_ysize();
00468
00469 EMData *mask1=new EMData(ys,ys,1);
00470 mask1->process_inplace("mask.gaussian",Dict("outer_radius", ys/2.0));
00471 EMData *mask2=mask1->copy();
00472 mask2->mult(-1.0f);
00473 mask2->add(1.0);
00474 mask2->process_inplace("mask.decayedge2d",Dict("width",4));
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495 }
00496
00497 void AmpweightFourierProcessor::process_inplace(EMData * image)
00498 {
00499 EMData *fft;
00500 float *fftd;
00501 int i,f=0;
00502
00503
00504
00505 if (!image) {
00506 LOGWARN("NULL Image");
00507 return;
00508 }
00509
00510 if (!image->is_complex()) {
00511 fft = image->do_fft();
00512 fftd = fft->get_data();
00513 f=1;
00514 }
00515 else {
00516 fft=image;
00517 fftd=image->get_data();
00518 }
00519 float *sumd = NULL;
00520 if (sum) sumd=sum->get_data();
00521
00522 int n = fft->get_xsize()*fft->get_ysize()*fft->get_zsize();
00523 for (i=0; i<n; i+=2) {
00524 float c;
00525 if (dosqrt) c=pow(fftd[i]*fftd[i]+fftd[i+1]*fftd[i+1],0.25f);
00526 #ifdef _WIN32
00527 else c = static_cast<float>(_hypot(fftd[i],fftd[i+1]));
00528 #else
00529 else c = static_cast<float>(hypot(fftd[i],fftd[i+1]));
00530 #endif //_WIN32
00531 if (c==0) c=1.0e-30f;
00532 fftd[i]*=c;
00533 fftd[i+1]*=c;
00534 if (sumd) { sumd[i]+=c; sumd[i+1]+=0; }
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547 }
00548
00549 if (f) {
00550 fft->update();
00551 EMData *ift=fft->do_ift();
00552 memcpy(image->get_data(),ift->get_data(),n*sizeof(float));
00553 delete fft;
00554 delete ift;
00555 }
00556
00557 sum->update();
00558 image->update();
00559
00560 }
00561
00562 void LinearPyramidProcessor::process_inplace(EMData *image) {
00563
00564 if (image->get_zsize()!=1) { throw ImageDimensionException("Only 2-D images supported"); }
00565
00566 float *d=image->get_data();
00567 int nx=image->get_xsize();
00568 int ny=image->get_ysize();
00569
00570 for (int y=0; y<ny; y++) {
00571 for (int x=0; x<nx; x++) {
00572 int l=x+y*nx;
00573 d[l]*=1.0f-abs(x-nx/2)*abs(y-ny/2)*4.0f/(nx*ny);
00574 }
00575 }
00576 image->update();
00577 }
00578
00579 EMData * Wiener2DAutoAreaProcessor::process(const EMData * image)
00580 {
00581
00582 EMData *ret = 0;
00583 const EMData *fft;
00584 float *fftd;
00585 int f=0;
00586
00587 if (!image) {
00588 LOGWARN("NULL Image");
00589 return ret;
00590 }
00591 throw NullPointerException("Processor not yet implemented");
00592
00593 if (!image->is_complex()) {
00594 fft = image->do_fft();
00595 fftd = fft->get_data();
00596 f=1;
00597 }
00598 else {
00599 fft=image;
00600 fftd=image->get_data();
00601 }
00602
00603 return ret;
00604 }
00605
00606 void Wiener2DAutoAreaProcessor::process_inplace(EMData *image) {
00607 EMData *tmp=process(image);
00608 memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
00609 delete tmp;
00610 image->update();
00611 return;
00612 }
00613
00614
00615 EMData * Wiener2DFourierProcessor::process(const EMData *in)
00616 {
00617 const EMData *in2 = 0;
00618 if (in->is_complex()) in2=in;
00619 else in=in->do_fft();
00620
00621 EMData *filt = in->copy_head();
00622 Ctf *ictf = ctf;
00623
00624 if (!ictf) ctf=(Ctf *)in->get_attr("ctf");
00625
00626 ictf->compute_2d_complex(filt,Ctf::CTF_WIENER_FILTER);
00627 filt->mult(*in2);
00628 EMData *ret=filt->do_ift();
00629
00630 delete filt;
00631 if (!in->is_complex()) delete in2;
00632
00633 if(!ictf && ctf) {delete ctf; ctf=0;}
00634 return(ret);
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663 return ret;
00664
00665 }
00666
00667 void Wiener2DFourierProcessor::process_inplace(EMData *image) {
00668 EMData *tmp=process(image);
00669 memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
00670 delete tmp;
00671 image->update();
00672 return;
00673 }
00674
00675 void LinearRampFourierProcessor::create_radial_func(vector < float >&radial_mask) const
00676 {
00677 Assert(radial_mask.size() > 0);
00678 for (size_t i = 0; i < radial_mask.size(); i++) {
00679 radial_mask[i] = (float)i;
00680 }
00681 }
00682
00683 void LowpassSharpCutoffProcessor::create_radial_func(vector < float >&radial_mask) const
00684 {
00685 Assert(radial_mask.size() > 0);
00686 float x = 0.0f , step = 0.5f/radial_mask.size();
00687 for (size_t i = 0; i < radial_mask.size(); i++) {
00688 if (x <= lowpass) {
00689 radial_mask[i] = 1.0f;
00690 }
00691 else {
00692 radial_mask[i] = 0;
00693 }
00694 x += step;
00695 }
00696 }
00697
00698
00699 void HighpassSharpCutoffProcessor::create_radial_func(vector < float >&radial_mask) const
00700 {
00701 Assert(radial_mask.size() > 0);
00702 float x = 0.0f , step = 0.5f/radial_mask.size();
00703 for (size_t i = 0; i < radial_mask.size(); i++) {
00704 if (x >= highpass) {
00705 radial_mask[i] = 1.0f;
00706 }
00707 else {
00708 radial_mask[i] = 0;
00709 }
00710 x += step;
00711 }
00712 }
00713
00714
00715 void LowpassGaussProcessor::create_radial_func(vector < float >&radial_mask) const
00716 {
00717
00718
00719 float x = 0.0f , step = 0.5f/radial_mask.size();
00720 float sig = 1;
00721 if (lowpass > 0) {
00722 sig = -1;
00723 }
00724
00725 for (size_t i = 0; i < radial_mask.size(); i++) {
00726 radial_mask[i] = exp(sig * x * x / (lowpass * lowpass));
00727 x += step;
00728 }
00729
00730 }
00731
00732 void HighpassGaussProcessor::create_radial_func(vector < float >&radial_mask) const
00733 {
00734 Assert(radial_mask.size() > 0);
00735 float x = 0.0f , step = 0.5f/radial_mask.size();
00736 for (size_t i = 0; i < radial_mask.size(); i++) {
00737 radial_mask[i] = 1.0f - exp(-x * x / (highpass * highpass));
00738 x += step;
00739 }
00740 }
00741
00742 void HighpassAutoPeakProcessor::preprocess(EMData * image)
00743 {
00744 if(params.has_key("apix")) {
00745 image->set_attr("apix_x", (float)params["apix"]);
00746 image->set_attr("apix_y", (float)params["apix"]);
00747 image->set_attr("apix_z", (float)params["apix"]);
00748 }
00749
00750 const Dict dict = image->get_attr_dict();
00751
00752 if( params.has_key("cutoff_abs") ) {
00753 highpass = params["cutoff_abs"];
00754 }
00755 else if( params.has_key("cutoff_freq") ) {
00756 highpass = (float)params["cutoff_freq"] * (float)dict["apix_x"] * (float)dict["nx"] / 2.0f;
00757 }
00758 else if( params.has_key("cutoff_pixels") ) {
00759 highpass = (float)params["cutoff_pixels"] / (float)dict["nx"];
00760 }
00761 }
00762
00763 void HighpassAutoPeakProcessor::create_radial_func(vector < float >&radial_mask, EMData *image) const
00764 {
00765 unsigned int c;
00766
00767
00768 for (c=2; c<radial_mask.size(); c++) if (radial_mask[c-1]<=radial_mask[c]) break;
00769 if (c>highpass) c=highpass;
00770
00771 radial_mask[0]=0.0;
00772
00773 for (unsigned int i=1; i<radial_mask.size(); i++) radial_mask[i]=(i<=c?0.0:1.0);
00774
00775 printf("%f %d\n",highpass,c);
00776
00777
00778 }
00779
00780
00781 void LowpassTanhProcessor::create_radial_func(vector < float >&radial_mask) const
00782 {
00783 Assert(radial_mask.size() > 0);
00784 float x = 0.0f , step = 0.5f/radial_mask.size();
00785 for (size_t i = 0; i < radial_mask.size(); i++) {
00786 radial_mask[i] = tanh((lowpass - x)*60.0f) / 2.0f + 0.5f;
00787 x += step;
00788 }
00789 }
00790
00791 void HighpassTanhProcessor::create_radial_func(vector < float >&radial_mask) const
00792 {
00793 Assert(radial_mask.size() > 0);
00794 float x = 0.0f , step = 0.5f/radial_mask.size();
00795 for (size_t i = 0; i < radial_mask.size(); i++) {
00796 radial_mask[i] = tanh(x - highpass) / 2.0f + 0.5f;
00797 x += step;
00798 }
00799 }
00800
00801
00802 void HighpassButterworthProcessor::create_radial_func(vector < float >&radial_mask) const
00803 {
00804 Assert(radial_mask.size() > 0);
00805 float x = 0.0f , step = 0.5f/radial_mask.size();
00806 for (size_t i = 0; i < radial_mask.size(); i++) {
00807 float t = highpass / 1.5f / (x + 0.001f);
00808 radial_mask[i] = 1.0f / (1.0f + t * t);
00809 x += step;
00810 }
00811 }
00812
00813
00814 void LinearRampProcessor::create_radial_func(vector < float >&radial_mask) const
00815 {
00816 Assert(radial_mask.size() > 0);
00817 float x = 0.0f , step = 0.5f/radial_mask.size();
00818 size_t size=radial_mask.size();
00819 for (size_t i = 0; i < size; i++) {
00820 radial_mask[i] = intercept + ((slope - intercept) * i) / (size - 1.0f);
00821 x += step;
00822 }
00823 }
00824
00825
00826 void RealPixelProcessor::process_inplace(EMData * image)
00827 {
00828 if (!image) {
00829 LOGWARN("NULL Image");
00830 return;
00831 }
00832
00833 maxval = image->get_attr("maximum");
00834 mean = image->get_attr("mean");
00835 sigma = image->get_attr("sigma");
00836
00837 calc_locals(image);
00838
00839 size_t size = (size_t)image->get_xsize() *
00840 (size_t)image->get_ysize() *
00841 (size_t)image->get_zsize();
00842 float *data = image->get_data();
00843
00844 for (size_t i = 0; i < size; i++) {
00845 process_pixel(&data[i]);
00846 }
00847 image->update();
00848 }
00849
00850 void CoordinateProcessor::process_inplace(EMData * image)
00851 {
00852 if (!image) {
00853 LOGWARN("NULL Image");
00854 return;
00855 }
00856
00857 maxval = image->get_attr("maximum");
00858 mean = image->get_attr("mean");
00859 sigma = image->get_attr("sigma");
00860 nx = image->get_xsize();
00861 ny = image->get_ysize();
00862 nz = image->get_zsize();
00863 is_complex = image->is_complex();
00864
00865 calc_locals(image);
00866
00867
00868 if (!is_valid()) {
00869 return;
00870 }
00871
00872 float *data = image->get_data();
00873 size_t i = 0;
00874
00875 for (int z = 0; z < nz; z++) {
00876 for (int y = 0; y < ny; y++) {
00877 for (int x = 0; x < nx; x++) {
00878 process_pixel(&data[i], x, y, z);
00879 ++i;
00880 }
00881 }
00882 }
00883 image->update();
00884 }
00885
00886 void PaintProcessor::process_inplace(EMData *image) {
00887 int nx=image->get_xsize();
00888 int ny=image->get_ysize();
00889 int nz=image->get_zsize();
00890
00891 if (nz==1) {
00892 float r;
00893 for (int j=(y<r2?0:y-r2); j<(y+r2>ny?ny:y+r2); j++) {
00894 for (int i=(x<r2?0:x-r2); i<(x+r2>nx?nx:x+r2); i++) {
00895 r=sqrt(float(Util::square(i-x)+Util::square(j-y)));
00896 if (r>r2 && r>r1) continue;
00897 if (r>r1) image->set_value_at(i,j,0,v2*(r-r1)/(r2-r1)+v1*(r2-r)/(r2-r1));
00898 else image->set_value_at(i,j,0,v1);
00899 }
00900 }
00901 }
00902 else {
00903 float r;
00904 for (int k=(z<r2?0:z-r2); k<(z+r2>nz?nz:z+r2); k++) {
00905 for (int j=(y<r2?0:y-r2); j<(y+r2>ny?ny:y+r2); j++) {
00906 for (int i=(x<r2?0:x-r2); i<(x+r2>nx?nx:x+r2); i++) {
00907 r=sqrt(float(Util::square(i-x)+Util::square(j-y)+Util::square(k-z)));
00908 if (r>r2 && r>r1) continue;
00909 if (r>r1) image->set_value_at(i,j,k,v2*(r-r1)/(r2-r1)+v1*(r2-r)/(r2-r1));
00910 else image->set_value_at(i,j,k,v1);
00911 }
00912 }
00913 }
00914 }
00915 image->update();
00916 }
00917
00918 void WatershedProcessor::process_inplace(EMData * image) {
00919 vector<float> xpoints = params["xpoints"];
00920 vector<float> ypoints = params["ypoints"];
00921 vector<float> zpoints = params["zpoints"];
00922
00923 vector<int> x(xpoints.begin(),xpoints.end());
00924 vector<int> y(ypoints.begin(),ypoints.end());
00925 vector<int> z(zpoints.begin(),zpoints.end());
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939 float minval = params["minval"];
00940
00941 EMData* mask = new EMData(*image);
00942 mask->to_zero();
00943
00944
00945 for(unsigned int i = 0; i < xpoints.size(); ++i) {
00946 try {
00947 mask->set_value_at(x[i],y[i],z[i], (float)(i+1));
00948 } catch (...) {
00949 continue;
00950 }
00951 }
00952 mask->write_image("seeds2.mrc");
00953
00954
00955
00956
00957
00958
00959
00960 while( true ) {
00961 bool cont= false;
00962 for(unsigned int j = 0; j < xpoints.size(); ++j)
00963 {
00964
00965 Vec3i coord(x[j],y[j],z[j]);
00966 vector<Vec3i> region;
00967 region.push_back(coord);
00968 vector<Vec3i> find_region_input = region;
00969 while (true) {
00970 vector<Vec3i> v = find_region(mask,find_region_input, j+1, region);
00971 if (v.size() == 0 ) break;
00972 else find_region_input = v;
00973 }
00974
00975 vector<Vec3i> tmp(region.begin(),region.end());
00976 region.clear();
00977 for(vector<Vec3i>::const_iterator it = tmp.begin(); it != tmp.end(); ++it ) {
00978 vector<Vec3i> tmp2 = watershed(mask, image, minval, *it, j+1);
00979 copy(tmp2.begin(),tmp2.end(),back_inserter(region));
00980 }
00981 if (region.size() != 0) cont = true;
00982 }
00983
00984 if (!cont) break;
00985 }
00986
00987
00988 memcpy(image->get_data(),mask->get_data(),sizeof(float)*image->get_size());
00989 image->update();
00990 }
00991
00992
00993 vector<Vec3i > WatershedProcessor::find_region(EMData* mask,const vector<Vec3i >& coords, const int mask_value, vector<Vec3i >& region)
00994 {
00995 static vector<Vec3i> two_six_connected;
00996 if (two_six_connected.size() == 0) {
00997 for(int i = -1; i <= 1; ++i) {
00998 for(int j = -1; j <= 1; ++j) {
00999 for(int k = -1; k <= 1; ++k) {
01000 if ( j != 0 || i != 0 || k != 0) {
01001 two_six_connected.push_back(Vec3i(i,j,k));
01002 }
01003 }
01004 }
01005 }
01006 }
01007
01008 vector<Vec3i> ret;
01009 for(vector<Vec3i>::const_iterator it = two_six_connected.begin(); it != two_six_connected.end(); ++it ) {
01010 for(vector<Vec3i>::const_iterator it2 = coords.begin(); it2 != coords.end(); ++it2 ) {
01011 if (mask->get_value_at((*it2)[0],(*it2)[1],(*it2)[2]) != mask_value) throw;
01012 Vec3i c = (*it)+(*it2);
01013
01014 if ( c[0] < 0 || c[0] >= mask->get_xsize()) continue;
01015 if ( c[1] < 0 || c[1] >= mask->get_ysize()) continue;
01016 if ( c[2] < 0 || c[2] >= mask->get_zsize()) continue;
01017
01018 if( mask->get_value_at(c[0],c[1],c[2]) == mask_value ) {
01019 if (find(ret.begin(),ret.end(),c) == ret.end()) {
01020 if (find(region.begin(),region.end(),c) == region.end()) {
01021 region.push_back(c);
01022 ret.push_back(c);
01023 }
01024 }
01025 }
01026 }
01027 }
01028 return ret;
01029 }
01030
01031 vector<Vec3i > WatershedProcessor::watershed(EMData* mask, EMData* image, const float& threshold, const Vec3i& coordinate, const int mask_value)
01032 {
01033 static vector<Vec3i> two_six_connected;
01034 if (two_six_connected.size() == 0) {
01035 for(int i = -1; i <= 1; ++i) {
01036 for(int j = -1; j <= 1; ++j) {
01037 for(int k = -1; k <= 1; ++k) {
01038 if ( j != 0 || i != 0 || k != 0) {
01039 two_six_connected.push_back(Vec3i(i,j,k));
01040 }
01041 }
01042 }
01043 }
01044 }
01045
01046 if (mask->get_value_at(coordinate[0],coordinate[1],coordinate[2]) != mask_value) throw;
01047
01048 vector<Vec3i> ret;
01049 for(vector<Vec3i>::const_iterator it = two_six_connected.begin(); it != two_six_connected.end(); ++it ) {
01050 Vec3i c = (*it)+coordinate;
01051
01052 if ( c[0] < 0 || c[0] >= image->get_xsize()) continue;
01053 if ( c[1] < 0 || c[1] >= image->get_ysize()) continue;
01054 if ( c[2] < 0 || c[2] >= image->get_zsize()) continue;
01055
01056
01057 if( image->get_value_at(c[0],c[1],c[2]) != 0 && (mask->get_value_at(c[0],c[1],c[2]) == 0 )) {
01058
01059 mask->set_value_at(c[0],c[1],c[2], (float)mask_value);
01060 ret.push_back(c);
01061 }
01062 }
01063 return ret;
01064 }
01065
01066
01067 void CircularMaskProcessor::calc_locals(EMData *)
01068 {
01069 xc = nx / 2.0f + dx;
01070 yc = ny / 2.0f + dy;
01071 zc = nz / 2.0f + dz;
01072
01073
01074 if (outer_radius <= 0) {
01075 outer_radius = nx / 2;
01076 outer_radius_square = outer_radius * outer_radius;
01077 }
01078
01079 if (inner_radius <= 0) {
01080 inner_radius_square = 0;
01081 }
01082 }
01083
01084
01085 void MaskEdgeMeanProcessor::calc_locals(EMData * image)
01086 {
01087 if (!image) {
01088 throw NullPointerException("NULL image");
01089 }
01090 int nitems = 0;
01091 float sum = 0;
01092 float *data = image->get_data();
01093 size_t i = 0;
01094
01095 for (int z = 0; z < nz; ++z) {
01096 for (int y = 0; y < ny; ++y) {
01097 for (int x = 0; x < nx; ++x) {
01098 float x1 = sqrt((x - xc) * (x - xc) + (y - yc) * (y - yc) + (z - zc) * (z - zc));
01099 if (x1 <= outer_radius + ring_width && x1 >= outer_radius - ring_width) {
01100 sum += data[i];
01101 ++nitems;
01102 }
01103 ++i;
01104 }
01105 }
01106 }
01107
01108 ring_avg = sum / nitems;
01109 }
01110
01111 void ComplexPixelProcessor::process_inplace(EMData * image)
01112 {
01113 if (!image) {
01114 LOGWARN("NULL image");
01115 return;
01116 }
01117 if (!image->is_complex()) {
01118 LOGWARN("cannot apply complex processor on a real image. Nothing is done.");
01119 return;
01120 }
01121
01122 size_t size = (size_t)image->get_xsize() *
01123 (size_t)image->get_ysize() *
01124 (size_t)image->get_zsize();
01125 float *data = image->get_data();
01126
01127 image->ri2ap();
01128
01129 for (size_t i = 0; i < size; i += 2) {
01130 process_pixel(data);
01131 data += 2;
01132 }
01133
01134 image->update();
01135 image->ap2ri();
01136 }
01137
01138
01139
01140 void AreaProcessor::process_inplace(EMData * image)
01141 {
01142 if (!image) {
01143 LOGWARN("NULL Image");
01144 return;
01145 }
01146
01147 float *data = image->get_data();
01148
01149 nx = image->get_xsize();
01150 ny = image->get_ysize();
01151 nz = image->get_zsize();
01152
01153 int n = (areasize - 1) / 2;
01154 matrix_size = areasize * areasize;
01155
01156 if (nz > 1) {
01157 matrix_size *= areasize;
01158 }
01159
01160 float *matrix = new float[matrix_size];
01161 kernel = new float[matrix_size];
01162
01163 int cpysize = areasize * (int) sizeof(float);
01164 int start = (nx * ny + nx + 1) * n;
01165
01166 int xend = nx - n;
01167 int yend = ny - n;
01168
01169 int zstart = n;
01170 int zend = nz - n;
01171
01172 int zbox_start = 0;
01173 int zbox_end = areasize;
01174
01175 if (nz == 1) {
01176 zstart = 0;
01177 zend = 1;
01178 zbox_end = 1;
01179 }
01180
01181 size_t nsec = (size_t)nx * (size_t)ny;
01182 int box_nsec = areasize * areasize;
01183
01184 create_kernel();
01185
01186 size_t total_size = (size_t)nx * (size_t)ny * (size_t)nz;
01187 float *data2 = new float[total_size];
01188 memcpy(data2, data, total_size * sizeof(float));
01189
01190 size_t k;
01191 for (int z = zstart; z < zend; z++) {
01192 for (int y = n; y < yend; y++) {
01193 for (int x = n; x < xend; x++) {
01194
01195 k = z * nsec + y * nx + x;
01196
01197 for (int bz = zbox_start; bz < zbox_end; bz++) {
01198 for (int by = 0; by < areasize; by++) {
01199 memcpy(&matrix[bz * box_nsec + by * areasize],
01200 &data2[k - start + bz * nsec + by * nx], cpysize);
01201 }
01202 }
01203
01204 process_pixel(&data[k], (float) x, (float) y, (float) z, matrix);
01205 }
01206 }
01207 }
01208
01209 if( matrix )
01210 {
01211 delete[]matrix;
01212 matrix = 0;
01213 }
01214
01215 if( kernel )
01216 {
01217 delete[]kernel;
01218 kernel = 0;
01219 }
01220 image->update();
01221 }
01222
01223
01224 void LaplacianProcessor::create_kernel() const
01225 {
01226 if (nz == 1) {
01227 memset(kernel, 0, areasize * areasize);
01228 kernel[1] = -0.25f;
01229 kernel[3] = -0.25f;
01230 kernel[5] = -0.25f;
01231 kernel[7] = -0.25f;
01232 kernel[4] = 1;
01233 }
01234 else {
01235 memset(kernel, 0, areasize * areasize * areasize);
01236 kernel[4] = -1.0f / 6.0f;
01237 kernel[10] = -1.0f / 6.0f;
01238 kernel[12] = -1.0f / 6.0f;
01239 kernel[14] = -1.0f / 6.0f;
01240 kernel[16] = -1.0f / 6.0f;
01241 kernel[22] = -1.0f / 6.0f;
01242 kernel[13] = 1;
01243 }
01244 }
01245
01246 void BoxStatProcessor::process_inplace(EMData * image)
01247 {
01248 if (!image) {
01249 LOGWARN("NULL Image");
01250 return;
01251 }
01252
01253 int nx = image->get_xsize();
01254 int ny = image->get_ysize();
01255 int nz = image->get_zsize();
01256
01257 int n = params.set_default("radius",1);
01258 int areasize = 2 * n + 1;
01259
01260 int matrix_size = areasize * areasize;
01261 if (nz > 1) {
01262 matrix_size *= areasize;
01263 }
01264
01265 float *array = new float[matrix_size];
01266
01267
01268 float *data = image->get_data();
01269 size_t total_size = (size_t)nx * (size_t)ny * (size_t)nz;
01270 float *data2 = new float[total_size];
01271 memcpy(data2, data, total_size * sizeof(float));
01272
01273 int z_begin = 0;
01274 int z_end = 1;
01275 int nzz=0;
01276 if (nz > 1) {
01277 z_begin = n;
01278 z_end = nz - n;
01279 nzz=n;
01280 }
01281
01282 int nxy = nx * ny;
01283
01284 for (int k = z_begin; k < z_end; k++) {
01285 size_t knxy = k * nxy;
01286
01287 for (int j = n; j < ny - n; j++) {
01288 int jnx = j * nx;
01289
01290 for (int i = n; i < nx - n; i++) {
01291 size_t s = 0;
01292
01293 for (int i2 = i - n; i2 <= i + n; i2++) {
01294 for (int j2 = j - n; j2 <= j + n; j2++) {
01295 for (int k2 = k - nzz; k2 <= k + nzz; k2++) {
01296 array[s] = data2[i2 + j2 * nx + k2 * nxy];
01297 ++s;
01298 }
01299 }
01300 }
01301
01302 process_pixel(&data[i + jnx + knxy], array, matrix_size);
01303 }
01304 }
01305 }
01306
01307 image->update();
01308
01309 if( data2 )
01310 {
01311 delete[]data2;
01312 data2 = 0;
01313 }
01314 }
01315
01316 void DiffBlockProcessor::process_inplace(EMData * image)
01317 {
01318 if (!image) {
01319 LOGWARN("NULL Image");
01320 return;
01321 }
01322
01323 int nz = image->get_zsize();
01324
01325 if (nz > 1) {
01326 LOGERR("%s Processor doesn't support 3D", get_name().c_str());
01327 throw ImageDimensionException("3D model not supported");
01328 }
01329
01330 int nx = image->get_xsize();
01331 int ny = image->get_ysize();
01332
01333 int v1 = params["cal_half_width"];
01334 int v2 = params["fill_half_width"];
01335
01336 int v0 = v1 > v2 ? v1 : v2;
01337
01338 if (v2 <= 0) {
01339 v2 = v1;
01340 }
01341
01342 float *data = image->get_data();
01343
01344 for (int y = v0; y <= ny - v0 - 1; y += v2) {
01345 for (int x = v0; x <= nx - v0 - 1; x += v2) {
01346
01347 float sum = 0;
01348 for (int y1 = y - v1; y1 <= y + v1; y1++) {
01349 for (int x1 = x - v1; x1 <= x + v1; x1++) {
01350 sum += data[x1 + y1 * nx];
01351 }
01352 }
01353 float mean = sum / ((v1 * 2 + 1) * (v1 * 2 + 1));
01354
01355 for (int j = y - v2; j <= y + v2; j++) {
01356 for (int i = x - v2; i <= x + v2; i++) {
01357 data[i + j * nx] = mean;
01358 }
01359 }
01360 }
01361 }
01362
01363 image->update();
01364 }
01365
01366
01367 void CutoffBlockProcessor::process_inplace(EMData * image)
01368 {
01369 if (!image) {
01370 LOGWARN("NULL Image");
01371 return;
01372 }
01373 int nz = image->get_zsize();
01374
01375 if (nz > 1) {
01376 LOGERR("%s Processor doesn't support 3D", get_name().c_str());
01377 throw ImageDimensionException("3D model not supported");
01378 }
01379
01380 int nx = image->get_xsize();
01381 int ny = image->get_ysize();
01382
01383 float value1 = params["value1"];
01384 float value2 = params["value2"];
01385
01386 int v1 = (int) value1;
01387 int v2 = (int) value2;
01388 if (v2 > v1 / 2) {
01389 LOGERR("invalid value2 '%f' in CutoffBlockProcessor", value2);
01390 return;
01391 }
01392
01393 if (v2 <= 0) {
01394 v2 = v1;
01395 }
01396
01397 float *data = image->get_data();
01398 int y = 0, x = 0;
01399 for (y = 0; y <= ny - v1; y += v1) {
01400 for (x = 0; x <= nx - v1; x += v1) {
01401
01402 EMData *clip = image->get_clip(Region(x, y, v1, v1));
01403 EMData *fft = clip->do_fft();
01404
01405 float *fft_data = fft->get_data();
01406 float sum = 0;
01407 int nitems = 0;
01408
01409 for (int i = -v2; i < v2; i++) {
01410 for (int j = 0; j < v2; j++) {
01411 if (j == 0 && i == 0) {
01412 continue;
01413 }
01414
01415 #ifdef _WIN32
01416 if (_hypot(j, i) < value2) {
01417 #else
01418 if (hypot(j, i) < value2) {
01419 #endif
01420 int t = j * 2 + (i + v1 / 2) * (v1 + 2);
01421 sum += (fft_data[t] * fft_data[t] + fft_data[t + 1] * fft_data[t + 1]);
01422 nitems++;
01423 }
01424 }
01425 }
01426
01427 if( clip )
01428 {
01429 delete clip;
01430 clip = 0;
01431 }
01432
01433 float mean = sum / nitems;
01434
01435 for (int i = y; i < y + v1; i++) {
01436 for (int j = x; j < x + v1; j++) {
01437 data[i * nx + j] = mean;
01438 }
01439 }
01440 }
01441 }
01442
01443 memset(&data[y * nx], 0, (ny - y) * nx * sizeof(float));
01444
01445 for (int i = 0; i < ny; i++) {
01446 memset(&data[i * nx + x], 0, (nx - x) * sizeof(float));
01447 }
01448
01449 image->update();
01450 }
01451
01452 void MedianShrinkProcessor::process_inplace(EMData * image)
01453 {
01454 if (image->is_complex()) throw ImageFormatException("Error, the median shrink processor does not work on complex images");
01455
01456 int shrink_factor = params.set_default("n",0);
01457 if (shrink_factor <= 1) {
01458 throw InvalidValueException(shrink_factor,
01459 "median shrink: shrink factor must > 1");
01460 }
01461
01462 int nx = image->get_xsize();
01463 int ny = image->get_ysize();
01464 int nz = image->get_zsize();
01465
01466
01467
01468
01469
01470
01471 int shrunken_nx = nx / shrink_factor;
01472 int shrunken_ny = ny / shrink_factor;
01473 int shrunken_nz = 1;
01474 if (nz > 1) shrunken_nz = nz / shrink_factor;
01475
01476 EMData* copy = image->copy();
01477 image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
01478 accrue_median(image,copy,shrink_factor);
01479 image->update();
01480 if( copy )
01481 {
01482 delete copy;
01483 copy = 0;
01484 }
01485 }
01486
01487
01488 EMData* MedianShrinkProcessor::process(const EMData *const image)
01489 {
01490 if (image->is_complex()) throw ImageFormatException("Error, the median shrink processor does not work on complex images");
01491
01492 int shrink_factor = params.set_default("n",0);
01493 if (shrink_factor <= 1) {
01494 throw InvalidValueException(shrink_factor,
01495 "median shrink: shrink factor must > 1");
01496 }
01497 int nx = image->get_xsize();
01498 int ny = image->get_ysize();
01499 int nz = image->get_zsize();
01500
01501
01502
01503
01504
01505
01506
01507 int shrunken_nx = nx / shrink_factor;
01508 int shrunken_ny = ny / shrink_factor;
01509 int shrunken_nz = 1;
01510 if (nz > 1) shrunken_nz = nz / shrink_factor;
01511
01512
01513 EMData *ret = image->copy_head();
01514 ret->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
01515
01516 accrue_median(ret,image,shrink_factor);
01517 ret->update();
01518 return ret;
01519 }
01520
01521 void MedianShrinkProcessor::accrue_median(EMData* to, const EMData* const from,const int shrink_factor)
01522 {
01523
01524 int nx_old = from->get_xsize();
01525 int ny_old = from->get_ysize();
01526
01527 int threed_shrink_factor = shrink_factor * shrink_factor;
01528 int z_shrink_factor = 1;
01529 if (from->get_zsize() > 1) {
01530 threed_shrink_factor *= shrink_factor;
01531 z_shrink_factor = shrink_factor;
01532 }
01533
01534 float *mbuf = new float[threed_shrink_factor];
01535
01536
01537 int nxy_old = nx_old * ny_old;
01538
01539 int nx = to->get_xsize();
01540 int ny = to->get_ysize();
01541 int nz = to->get_zsize();
01542 int nxy_new = nx * ny;
01543
01544 float * rdata = to->get_data();
01545 const float *const data_copy = from->get_const_data();
01546
01547 for (int l = 0; l < nz; l++) {
01548 int l_min = l * shrink_factor;
01549 int l_max = l * shrink_factor + z_shrink_factor;
01550 int cur_l = l * nxy_new;
01551
01552 for (int j = 0; j < ny; j++) {
01553 int j_min = j * shrink_factor;
01554 int j_max = (j + 1) * shrink_factor;
01555 int cur_j = j * nx + cur_l;
01556
01557 for (int i = 0; i < nx; i++) {
01558 int i_min = i * shrink_factor;
01559 int i_max = (i + 1) * shrink_factor;
01560
01561 int k = 0;
01562 for (int l2 = l_min; l2 < l_max; l2++) {
01563 int cur_l2 = l2 * nxy_old;
01564
01565 for (int j2 = j_min; j2 < j_max; j2++) {
01566 int cur_j2 = j2 * nx_old + cur_l2;
01567
01568 for (int i2 = i_min; i2 < i_max; i2++) {
01569 mbuf[k] = data_copy[i2 + cur_j2];
01570 ++k;
01571 }
01572 }
01573 }
01574
01575 for (k = 0; k < threed_shrink_factor / 2 + 1; k++) {
01576 for (int i2 = k + 1; i2 < threed_shrink_factor; i2++) {
01577 if (mbuf[i2] < mbuf[k]) {
01578 float f = mbuf[i2];
01579 mbuf[i2] = mbuf[k];
01580 mbuf[k] = f;
01581 }
01582 }
01583 }
01584
01585 rdata[i + cur_j] = mbuf[threed_shrink_factor / 2];
01586 }
01587 }
01588 }
01589
01590 if( mbuf )
01591 {
01592 delete[]mbuf;
01593 mbuf = 0;
01594 }
01595
01596 to->scale_pixel((float)shrink_factor);
01597 }
01598
01599 EMData* FFTResampleProcessor::process(const EMData *const image)
01600 {
01601 float sample_rate = params.set_default("n",0.0f);
01602 if (sample_rate <= 0.0F ) {
01603 throw InvalidValueException(sample_rate, "sample rate must be >0 ");
01604 }
01605
01606 EMData* result;
01607 if (image->is_complex()) result = image->copy();
01608 else result = image->do_fft();
01609 fft_resample(result,image,sample_rate);
01610
01611 result->update();
01612 result->scale_pixel(sample_rate);
01613 return result;
01614 }
01615
01616 void FFTResampleProcessor::process_inplace(EMData * image)
01617 {
01618 if (image->is_complex()) throw ImageFormatException("Error, the fft resampling processor does not work on complex images");
01619
01620
01621 float sample_rate = params.set_default("n",0.0f);
01622 if (sample_rate <= 0.0F ) {
01623 throw InvalidValueException(sample_rate, "sample rate (n) must be >0 ");
01624 }
01625
01626 fft_resample(image,image,sample_rate);
01627
01628 image->scale_pixel(sample_rate);
01629 image->update();
01630
01631
01632 }
01633
01634 void FFTResampleProcessor::fft_resample(EMData* to, const EMData *const from, const float& sample_rate) {
01635 int nx = from->get_xsize();
01636 int ny = from->get_ysize();
01637 int nz = from->get_zsize();
01638
01639 int new_nx = static_cast<int>( static_cast<float> (nx) / sample_rate);
01640 int new_ny = static_cast<int>( static_cast<float> (ny) / sample_rate);
01641 int new_nz = static_cast<int>( static_cast<float> (nz) / sample_rate);
01642
01643 if (new_nx == 0) throw UnexpectedBehaviorException("The resample rate causes the pixel dimensions in the x direction to go to zero");
01644 if (new_ny == 0) new_ny = 1;
01645 if (new_nz == 0) new_nz = 1;
01646
01647 int ndim = from->get_ndim();
01648 if ( ndim < 3 ) {
01649 new_nz = 1;
01650 }
01651 if ( ndim < 2 ) {
01652 new_ny = 1;
01653 }
01654
01655 int fft_x_correction = 1;
01656 if (new_nx % 2 == 0) fft_x_correction = 2;
01657
01658 int fft_y_correction = 0;
01659 if (ny != 1 && new_ny % 2 == 0 && ny % 2 == 1) fft_y_correction = 1;
01660 else if (ny != 1 && new_ny % 2 == 1 && ny % 2 == 0) fft_y_correction = -1;
01661
01662 int fft_z_correction = 0;
01663 if (nz != 1 && new_nz % 2 == 0 && nz % 2 == 1) fft_z_correction = 1;
01664 else if (nz != 1 && new_nz % 2 == 1 && nz % 2 == 0) fft_z_correction = -1;
01665
01666 if ( ! to->is_complex()) to->do_fft_inplace();
01667
01668 if (ndim != 1) to->process_inplace("xform.fourierorigin.tocenter");
01669
01670 Region clip(0,(ny-new_ny)/2-fft_y_correction,(nz-new_nz)/2-fft_z_correction,new_nx+fft_x_correction,new_ny,new_nz);
01671 to->clip_inplace(clip);
01672
01673 if (fft_x_correction == 1) to->set_fftodd(true);
01674 else to->set_fftodd(false);
01675
01676 if (ndim != 1) to->process_inplace("xform.fourierorigin.tocorner");
01677
01678 to->do_ift_inplace();
01679 to->depad_corner();
01680
01681 }
01682
01683
01684 EMData* MeanShrinkProcessor::process(const EMData *const image)
01685 {
01686 if (image->is_complex()) throw ImageFormatException("Error, the mean shrink processor does not work on complex images");
01687
01688 if (image->get_ndim() == 1) { throw ImageDimensionException("Error, mean shrink works only for 2D & 3D images"); }
01689
01690 float shrink_factor0 = params.set_default("n",0.0f);
01691 int shrink_factor = int(shrink_factor0);
01692 if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {
01693 throw InvalidValueException(shrink_factor0,
01694 "mean shrink: shrink factor must be >1 integer or 1.5");
01695 }
01696
01697 int nx = image->get_xsize();
01698 int ny = image->get_ysize();
01699 int nz = image->get_zsize();
01700
01701
01702
01703 if (shrink_factor0==1.5 ) {
01704 if (nz > 1 ) throw InvalidValueException(shrink_factor0, "mean shrink: only support 2D images for shrink factor = 1.5");
01705
01706 int shrunken_nx = (int(nx / 1.5)+1)/2*2;
01707 int shrunken_ny = (int(ny / 1.5)+1)/2*2;
01708 EMData* result = new EMData(shrunken_nx,shrunken_ny,1);
01709
01710 accrue_mean_one_p_five(result,image);
01711 result->update();
01712
01713 return result;
01714 }
01715
01716 int shrunken_nx = nx / shrink_factor;
01717 int shrunken_ny = ny / shrink_factor;
01718 int shrunken_nz = 1;
01719
01720 if (nz > 1) {
01721 shrunken_nz = nz / shrink_factor;
01722 }
01723
01724
01725 EMData* result = image->copy_head();
01726 result->set_size(shrunken_nx,shrunken_ny,shrunken_nz);
01727 accrue_mean(result,image,shrink_factor);
01728
01729 result->update();
01730
01731 return result;
01732 }
01733
01734 void MeanShrinkProcessor::process_inplace(EMData * image)
01735 {
01736 if (image->is_complex()) throw ImageFormatException("Error, the mean shrink processor does not work on complex images");
01737
01738 if (image->get_ndim() == 1) { throw ImageDimensionException("Error, mean shrink works only for 2D & 3D images"); }
01739
01740 float shrink_factor0 = params.set_default("n",0.0f);
01741 int shrink_factor = int(shrink_factor0);
01742 if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {
01743 throw InvalidValueException(shrink_factor0,
01744 "mean shrink: shrink factor must be >1 integer or 1.5");
01745 }
01746
01747
01748
01749
01750
01751
01752
01753 int nx = image->get_xsize();
01754 int ny = image->get_ysize();
01755 int nz = image->get_zsize();
01756
01757 if (shrink_factor0==1.5 ) {
01758 if (nz > 1 ) throw InvalidValueException(shrink_factor0, "mean shrink: only support 2D images for shrink factor = 1.5");
01759
01760 int shrunken_nx = (int(nx / 1.5)+1)/2*2;
01761 int shrunken_ny = (int(ny / 1.5)+1)/2*2;
01762
01763 EMData* orig = image->copy();
01764 image->set_size(shrunken_nx, shrunken_ny, 1);
01765 image->to_zero();
01766
01767 accrue_mean_one_p_five(image,orig);
01768
01769 if( orig ) {
01770 delete orig;
01771 orig = 0;
01772 }
01773 image->update();
01774
01775 return;
01776 }
01777
01778 accrue_mean(image,image,shrink_factor);
01779
01780 int shrunken_nx = nx / shrink_factor;
01781 int shrunken_ny = ny / shrink_factor;
01782 int shrunken_nz = 1;
01783 if (nz > 1) shrunken_nz = nz / shrink_factor;
01784
01785 image->update();
01786 image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
01787 }
01788
01789 void MeanShrinkProcessor::accrue_mean(EMData* to, const EMData* const from,const int shrink_factor)
01790 {
01791 const float * const data = from->get_const_data();
01792 float* rdata = to->get_data();
01793
01794 int nx = from->get_xsize();
01795 int ny = from->get_ysize();
01796 int nz = from->get_zsize();
01797 int nxy = nx*ny;
01798
01799
01800 int shrunken_nx = nx / shrink_factor;
01801 int shrunken_ny = ny / shrink_factor;
01802 int shrunken_nz = 1;
01803 int shrunken_nxy = shrunken_nx * shrunken_ny;
01804
01805 int normalize_shrink_factor = shrink_factor * shrink_factor;
01806 int z_shrink_factor = 1;
01807
01808 if (nz > 1) {
01809 shrunken_nz = nz / shrink_factor;
01810 normalize_shrink_factor *= shrink_factor;
01811 z_shrink_factor = shrink_factor;
01812 }
01813
01814 float invnormfactor = 1.0f/(float)normalize_shrink_factor;
01815
01816 for (int k = 0; k < shrunken_nz; k++) {
01817 int k_min = k * shrink_factor;
01818 int k_max = k * shrink_factor + z_shrink_factor;
01819 size_t cur_k = k * shrunken_nxy;
01820
01821 for (int j = 0; j < shrunken_ny; j++) {
01822 int j_min = j * shrink_factor;
01823 int j_max = j * shrink_factor + shrink_factor;
01824 size_t cur_j = j * shrunken_nx + cur_k;
01825
01826 for (int i = 0; i < shrunken_nx; i++) {
01827 int i_min = i * shrink_factor;
01828 int i_max = i * shrink_factor + shrink_factor;
01829
01830 float sum = 0;
01831 for (int kk = k_min; kk < k_max; kk++) {
01832 size_t cur_kk = kk * nxy;
01833
01834 for (int jj = j_min; jj < j_max; jj++) {
01835 size_t cur_jj = jj * nx + cur_kk;
01836 for (int ii = i_min; ii < i_max; ii++) {
01837 sum += data[ii + cur_jj];
01838 }
01839 }
01840 }
01841 rdata[i + cur_j] = sum * invnormfactor;
01842 }
01843 }
01844 }
01845 to->scale_pixel((float)shrink_factor);
01846 }
01847
01848
01849 void MeanShrinkProcessor::accrue_mean_one_p_five(EMData* to, const EMData * const from)
01850 {
01851 int nx0 = from->get_xsize(), ny0 = from->get_ysize();
01852
01853 int nx = to->get_xsize(), ny = to->get_ysize();
01854
01855 float *data = to->get_data();
01856 const float * const data0 = from->get_const_data();
01857
01858 for (int j = 0; j < ny; j++) {
01859 int jj = int(j * 1.5);
01860 float jw0 = 1.0F, jw1 = 0.5F;
01861 if ( j%2 ) {
01862 jw0 = 0.5F;
01863 jw1 = 1.0F;
01864 }
01865 for (int i = 0; i < nx; i++) {
01866 int ii = int(i * 1.5);
01867 float iw0 = 1.0F, iw1 = 0.5F;
01868 float w = 0.0F;
01869
01870 if ( i%2 ) {
01871 iw0 = 0.5F;
01872 iw1 = 1.0F;
01873 }
01874 if ( jj < ny0 ) {
01875 if ( ii < nx0 ) {
01876 data[j * nx + i] = data0[ jj * nx0 + ii ] * jw0 * iw0 ;
01877 w += jw0 * iw0 ;
01878 if ( ii+1 < nx0 ) {
01879 data[j * nx + i] += data0[ jj * nx0 + ii + 1] * jw0 * iw1;
01880 w += jw0 * iw1;
01881 }
01882 }
01883 if ( jj +1 < ny0 ) {
01884 if ( ii < nx0 ) {
01885 data[j * nx + i] += data0[ (jj+1) * nx0 + ii ] * jw1 * iw0;
01886 w += jw1 * iw0;
01887 if ( ii+1 < nx0 ) {
01888 data[j * nx + i] += data0[ (jj+1) * nx0 + ii + 1] * jw1 * iw1;
01889 w += jw1 * iw1;
01890 }
01891 }
01892 }
01893 }
01894 if ( w>0 ) data[j * nx + i] /= w;
01895 }
01896 }
01897
01898 to->update();
01899 to->scale_pixel((float)1.5);
01900 }
01901
01902
01903 template<class LogicOp>
01904 EMData* BooleanShrinkProcessor::process(const EMData *const image, Dict& params)
01905 {
01906
01907
01908
01909 if (!image) throw NullPointerException("Attempt to max shrink a null image");
01910
01911 if (image->is_complex() ) throw ImageFormatException("Can not max shrink a complex image");
01912
01913
01914 int shrink = params.set_default("n",2);
01915 int search = params.set_default("search",2);
01916
01917 if ( shrink < 0 ) throw InvalidValueException(shrink, "Can not shrink by a value less than 0");
01918
01919
01920 int nz = image->get_zsize();
01921 int ny = image->get_ysize();
01922 int nx = image->get_xsize();
01923
01924 if (nx == 1 && ny == 1 && nz == 1 ) return image->copy();
01925
01926 LogicOp op;
01927 EMData* return_image = new EMData();
01928
01929 int shrinkx = shrink;
01930 int shrinky = shrink;
01931 int shrinkz = shrink;
01932
01933 int searchx = search;
01934 int searchy = search;
01935 int searchz = search;
01936
01937
01938
01939 if ( shrinkx > nx ) shrinkx = nx;
01940 if ( shrinky > ny ) shrinky = ny;
01941 if ( shrinkz > nz ) shrinkz = nz;
01942
01943 if ( nz == 1 && ny == 1 )
01944 {
01945 return_image->set_size(nx/shrinkx);
01946 for(int i = 0; i < nx/shrinkx; ++i)
01947 {
01948 float tmp = op.get_start_val();
01949 for(int s=0; s < searchx; ++s)
01950 {
01951 int idx = shrinkx*i+s;
01952
01953 if ( idx > nx ) break;
01954 else
01955 {
01956 float val = image->get_value_at(idx);
01957 if ( op( val,tmp) ) tmp = val;
01958 }
01959 }
01960 return_image->set_value_at(i,tmp);
01961 }
01962 }
01963 else if ( nz == 1 )
01964 {
01965 int ty = ny/shrinky;
01966 int tx = nx/shrinkx;
01967 return_image->set_size(tx,ty);
01968 for(int y = 0; y < ty; ++y) {
01969 for(int x = 0; x < tx; ++x) {
01970 float tmp = op.get_start_val();
01971 for(int sy=0; sy < searchy; ++sy) {
01972 int yidx = shrinky*y+sy;
01973 if ( yidx >= ny) break;
01974 for(int sx=0; sx < searchx; ++sx) {
01975 int xidx = shrinkx*x+sx;
01976 if ( xidx >= nx) break;
01977
01978 float val = image->get_value_at(xidx,yidx);
01979 if ( op( val,tmp) ) tmp = val;
01980 }
01981 }
01982 return_image->set_value_at(x,y,tmp);
01983 }
01984 }
01985 }
01986 else
01987 {
01988 int tz = nz/shrinkz;
01989 int ty = ny/shrinky;
01990 int tx = nx/shrinkx;
01991
01992 return_image->set_size(tx,ty,tz);
01993 for(int z = 0; z < tz; ++z) {
01994 for(int y = 0; y < ty; ++y) {
01995 for(int x = 0; x < tx; ++x) {
01996 float tmp = op.get_start_val();
01997
01998 for(int sz=0; sz < searchz; ++sz) {
01999 int zidx = shrinkz*z+sz;
02000 if ( zidx >= nz) break;
02001
02002 for(int sy=0; sy < searchy; ++sy) {
02003 int yidx = shrinky*y+sy;
02004 if ( yidx >= ny) break;
02005
02006 for(int sx=0; sx < searchx; ++sx) {
02007 int xidx = shrinkx*x+sx;
02008 if ( xidx >= nx) break;
02009 float val = image->get_value_at(xidx,yidx,zidx);
02010 if ( op( val,tmp) ) tmp = val;
02011 }
02012 }
02013 }
02014 return_image->set_value_at(x,y,z,tmp);
02015 }
02016 }
02017 }
02018 }
02019 return_image->update();
02020
02021 return return_image;
02022 }
02023
02024 template<class LogicOp>
02025 void BooleanShrinkProcessor::process_inplace(EMData * image, Dict& params)
02026 {
02027
02028
02029 if (!image) throw NullPointerException("Attempt to max shrink a null image");
02030
02031 if (image->is_complex() ) throw ImageFormatException("Can not max shrink a complex image");
02032
02033
02034 int shrink = params.set_default("shrink",2);
02035 int search = params.set_default("search",2);
02036
02037 if ( shrink < 0 ) throw InvalidValueException(shrink, "Can not shrink by a value less than 0");
02038
02039
02040 int nz = image->get_zsize();
02041 int ny = image->get_ysize();
02042 int nx = image->get_xsize();
02043
02044 LogicOp op;
02045
02046 int shrinkx = shrink;
02047 int shrinky = shrink;
02048 int shrinkz = shrink;
02049
02050 int searchx = search;
02051 int searchy = search;
02052 int searchz = search;
02053
02054
02055
02056 if ( shrinkx > nx ) shrinkx = nx;
02057 if ( shrinky > ny ) shrinkx = ny;
02058 if ( shrinkz > nz ) shrinkx = nz;
02059
02060 if (nx == 1 && ny == 1 && nz == 1 ) return;
02061
02062 if ( nz == 1 && ny == 1 )
02063 {
02064 for(int i = 0; i < nx/shrink; ++i)
02065 {
02066 float tmp = op.get_start_val();
02067 for(int s=0; s < searchx; ++s)
02068 {
02069 int idx = shrinkx*i+s;
02070 if ( idx > nx ) break;
02071 else
02072 {
02073 float val = image->get_value_at(idx);
02074 if ( op( val,tmp) ) tmp = val;
02075 }
02076 }
02077 image->set_value_at(i,tmp);
02078 }
02079
02080 image->set_size(nx/shrinkx);
02081 }
02082 else if ( nz == 1 )
02083 {
02084 int ty = ny/shrinky;
02085 int tx = nx/shrinkx;
02086 for(int y = 0; y < ty; ++y) {
02087 for(int x = 0; x < tx; ++x) {
02088 float tmp = op.get_start_val();
02089 for(int sy=0; sy < searchy; ++sy) {
02090 int yidx = shrinky*y+sy;
02091 if ( yidx >= ny) break;
02092 for(int sx=0; sx < searchx; ++sx) {
02093 int xidx = shrinkx*x+sx;
02094 if ( xidx >= nx) break;
02095
02096 float val = image->get_value_at(xidx,yidx);
02097 if ( op( val,tmp) ) tmp = val;
02098 }
02099 }
02100 (*image)(x+tx*y) = tmp;
02101 }
02102 }
02103 image->set_size(tx,ty);
02104 }
02105 else
02106 {
02107 int tnxy = nx/shrinkx*ny/shrinky;
02108 int tz = nz/shrinkz;
02109 int ty = ny/shrinky;
02110 int tx = nx/shrinkx;
02111
02112 for(int z = 0; z < tz; ++z) {
02113 for(int y = 0; y < ty; ++y) {
02114 for(int x = 0; x < tx; ++x) {
02115 float tmp = op.get_start_val();
02116 for(int sz=0; sz < searchz; ++sz) {
02117 int zidx = shrinkz*z+sz;
02118 if ( zidx >= nz) break;
02119 for(int sy=0; sy < searchy; ++sy) {
02120 int yidx = shrinky*y+sy;
02121 if ( yidx >= ny) break;
02122 for(int sx=0; sx < shrinkx; ++sx) {
02123 int xidx = shrinkx*x+sx;
02124 if ( xidx >= nx) break;
02125
02126 float val = image->get_value_at(xidx,yidx,zidx);
02127 if ( op( val,tmp) ) tmp = val;
02128 }
02129 }
02130 }
02131 (*image)(x+tx*y+tnxy*z) = tmp;
02132 }
02133 }
02134 }
02135 image->set_size(tx,ty,tz);
02136 }
02137 image->update();
02138 }
02139
02140
02141
02142 void GradientRemoverProcessor::process_inplace(EMData * image)
02143 {
02144 if (!image) {
02145 LOGWARN("NULL Image");
02146 return;
02147 }
02148
02149 int nz = image->get_zsize();
02150 if (nz > 1) {
02151 LOGERR("%s Processor doesn't support 3D model", get_name().c_str());
02152 throw ImageDimensionException("3D model not supported");
02153 }
02154
02155 int nx = image->get_xsize();
02156 int ny = image->get_ysize();
02157 float *dy = new float[ny];
02158 float m = 0;
02159 float b = 0;
02160 float sum_y = 0;
02161 float *data = image->get_data();
02162
02163 for (int i = 0; i < ny; i++) {
02164 Util::calc_least_square_fit(nx, 0, data + i * nx, &m, &b, false);
02165 dy[i] = b;
02166 sum_y += m;
02167 }
02168
02169 float mean_y = sum_y / ny;
02170 float sum_x = 0;
02171 Util::calc_least_square_fit(ny, 0, dy, &sum_x, &b, false);
02172
02173 for (int j = 0; j < ny; j++) {
02174 for (int i = 0; i < nx; i++) {
02175 data[i + j * nx] -= i * sum_x + j * mean_y + b;
02176 }
02177 }
02178
02179 image->update();
02180 }
02181
02182 void FlattenBackgroundProcessor::process_inplace(EMData * image)
02183 {
02184
02185 EMData* mask = params.set_default("mask",(EMData*)0);
02186 int radius = params.set_default("radius",0);
02187
02188 if (radius != 0 && mask != 0) throw InvalidParameterException("Error - the mask and radius parameters are mutually exclusive.");
02189
02190 if (mask == 0 && radius == 0) throw InvalidParameterException("Error - you must specify either the mask or the radius parameter.");
02191
02192
02193 bool deletemask = false;
02194 if (radius != 0) {
02195 mask = new EMData;
02196 int n = image->get_ndim();
02197 if (n==1){
02198 mask->set_size(2*radius+1);
02199 } else if (n==2) {
02200 mask->set_size(2*radius+1,2*radius+1);
02201 }
02202 else {
02203 mask->set_size(2*radius+1,2*radius+1,2*radius+1);
02204 }
02205
02206 mask->process_inplace("testimage.circlesphere");
02207 }
02208
02209
02210 int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();
02211 int nx = image->get_xsize(); int ny = image->get_ysize(); int nz = image->get_zsize();
02212 int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;
02213 if (nz == 1) nzc = 1;
02214 if (ny == 1) nyc = 1;
02215
02216 if ( mnx > nx || mny > ny || mnz > nz)
02217 throw ImageDimensionException("Can not flatten using a mask that is larger than the image.");
02218
02219
02220 float normfac = 0.0;
02221 for (int i=0; i<mask->get_xsize()*mask->get_ysize()*mask->get_zsize(); ++i){
02222 normfac += mask->get_value_at(i);
02223 }
02224
02225
02226
02227 if (normfac == 0) throw InvalidParameterException("Error - the pixels in the mask sum to zero. This breaks the flattening procedure");
02228 normfac = 1.0f/normfac;
02229
02230
02231
02232
02233 Region r;
02234 if (ny == 1) r = Region((mnx-nxc)/2,nxc);
02235 else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
02236 else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
02237 mask->clip_inplace(r,0);
02238
02239
02240
02241
02242
02243
02244
02245 Region r2;
02246 if (ny == 1) r2 = Region((nx-nxc)/2,nxc);
02247 else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
02248 else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
02249 image->clip_inplace(r2,image->get_edge_mean());
02250
02251 EMData* m = image->convolute(mask);
02252
02253 m->mult(normfac);
02254
02255 m->process_inplace("xform.phaseorigin.tocenter");
02256
02257 image->sub(*m);
02258 delete m;
02259
02260 if (deletemask) {
02261 delete mask;
02262 } else {
02263 Region r;
02264 if (ny == 1) r = Region((nxc-mnx)/2,mnx);
02265 else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
02266 else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
02267 mask->clip_inplace(r);
02268 }
02269
02270 Region r3;
02271 if (ny == 1) r3 = Region((nxc-nx)/2,nx);
02272 else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
02273 else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
02274 image->clip_inplace(r3);
02275
02276
02277
02278
02279
02280 }
02281
02282 #include <gsl/gsl_linalg.h>
02283 void GradientPlaneRemoverProcessor::process_inplace(EMData * image)
02284 {
02285 if (!image) {
02286 LOGWARN("NULL Image");
02287 return;
02288 }
02289
02290 int nz = image->get_zsize();
02291 if (nz > 1) {
02292 LOGERR("%s Processor doesn't support 3D model", get_name().c_str());
02293 throw ImageDimensionException("3D map not supported");
02294 }
02295
02296 int nx = image->get_xsize();
02297 int ny = image->get_ysize();
02298 float *d = image->get_data();
02299 EMData *mask = 0;
02300 float *dm = 0;
02301 if (params.has_key("mask")) {
02302 mask = params["mask"];
02303 if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) {
02304 LOGERR("%s Processor requires same size mask image", get_name().c_str());
02305 throw ImageDimensionException("wrong size mask image");
02306 }
02307 dm = mask->get_data();
02308 }
02309 int count = 0;
02310 if (dm) {
02311 for(int i=0; i<nx*ny; i++) {
02312 if(dm[i]) count++;
02313 }
02314 }
02315 else {
02316 count = nx * ny;
02317 }
02318 if(count<3) {
02319 LOGERR("%s Processor requires at least 3 pixels to fit a plane", get_name().c_str());
02320 throw ImageDimensionException("too few usable pixels to fit a plane");
02321 }
02322
02323 gsl_vector *S=gsl_vector_calloc(3);
02324 gsl_matrix *A=gsl_matrix_calloc(count,3);
02325 gsl_matrix *V=gsl_matrix_calloc(3,3);
02326
02327 double m[3] = {0, 0, 0};
02328 int index=0;
02329 if (dm) {
02330 for(int j=0; j<ny; j++){
02331 for(int i=0; i<nx; i++){
02332 int ij=j*nx+i;
02333 if(dm[ij]) {
02334 m[0]+=i;
02335 m[1]+=j;
02336 m[2]+=d[ij];
02337
02338
02339 index++;
02340 }
02341 }
02342 }
02343 }
02344 else {
02345 for(int j=0; j<ny; j++){
02346 for(int i=0; i<nx; i++){
02347 int ij=j*nx+i;
02348 m[0]+=i;
02349 m[1]+=j;
02350 m[2]+=d[ij];
02351
02352
02353 index++;
02354 }
02355 }
02356 }
02357
02358 for(int i=0; i<3; i++) m[i]/=count;
02359
02360 index=0;
02361 if (dm) {
02362 for(int j=0; j<ny; j++){
02363 for(int i=0; i<nx; i++){
02364 int ij=j*nx+i;
02365 if(dm[ij]) {
02366
02367 gsl_matrix_set(A,index,0,i-m[0]);
02368 gsl_matrix_set(A,index,1,j-m[1]);
02369 gsl_matrix_set(A,index,2,d[ij]-m[2]);
02370 index++;
02371 }
02372 }
02373 }
02374 mask->update();
02375 }
02376 else {
02377 for(int j=0; j<ny; j++){
02378 for(int i=0; i<nx; i++){
02379 int ij=j*nx+i;
02380
02381 gsl_matrix_set(A,index,0,i-m[0]);
02382 gsl_matrix_set(A,index,1,j-m[1]);
02383 gsl_matrix_set(A,index,2,d[ij]-m[2]);
02384 index++;
02385 }
02386 }
02387 }
02388
02389
02390 gsl_linalg_SV_decomp_jacobi(A, V, S);
02391
02392 double n[3];
02393 for(int i=0; i<3; i++) n[i] = gsl_matrix_get(V, i, 2);
02394
02395 #ifdef DEBUG
02396 printf("S=%g,%g,%g\n",gsl_vector_get(S,0), gsl_vector_get(S,1), gsl_vector_get(S,2));
02397 printf("V[0,:]=%g,%g,%g\n",gsl_matrix_get(V,0,0), gsl_matrix_get(V,0,1),gsl_matrix_get(V,0,2));
02398 printf("V[1,:]=%g,%g,%g\n",gsl_matrix_get(V,1,0), gsl_matrix_get(V,1,1),gsl_matrix_get(V,1,2));
02399 printf("V[2,:]=%g,%g,%g\n",gsl_matrix_get(V,2,0), gsl_matrix_get(V,2,1),gsl_matrix_get(V,2,2));
02400 printf("Fitted plane: p0=%g,%g,%g\tn=%g,%g,%g\n",m[0],m[1],m[2],n[0],n[1],n[2]);
02401 #endif
02402
02403 int changeZero = 0;
02404 if (params.has_key("changeZero")) changeZero = params["changeZero"];
02405 if (changeZero) {
02406 for(int j=0; j<nx; j++){
02407 for(int i=0; i<ny; i++){
02408 int ij = j*nx+i;
02409 d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);
02410 }
02411 }
02412 }
02413 else {
02414 for(int j=0; j<nx; j++){
02415 for(int i=0; i<ny; i++){
02416 int ij = j*nx+i;
02417 if(d[ij]) d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);
02418 }
02419 }
02420 }
02421 image->update();
02422
02423 vector< float > planeParam;
02424 planeParam.resize(6);
02425 for(int i=0; i<3; i++) planeParam[i] = static_cast<float>(n[i]);
02426 for(int i=0; i<3; i++) planeParam[i+3] = static_cast<float>(m[i]);
02427 params["planeParam"]=EMObject(planeParam);
02428 }
02429
02430 void VerticalStripeProcessor::process_inplace(EMData * image)
02431 {
02432 if (!image) {
02433 LOGWARN("NULL Image");
02434 return;
02435 }
02436
02437 int nx = image->get_xsize();
02438 int ny = image->get_ysize();
02439 int nz = image->get_zsize();
02440
02441 float *data = image->get_data();
02442 float sigma = image->get_attr("sigma");
02443
02444 for (int k = 0; k < nz; k++) {
02445 for (int i = 0; i < nx; i++) {
02446 double sum = 0;
02447 for (int j = ny / 4; j < 3 * ny / 4; j++) {
02448 sum += data[i + j * nx];
02449 }
02450
02451 float mean = (float)sum / (ny / 2);
02452 for (int j = 0; j < ny; j++) {
02453 data[i + j * nx] = (data[i + j * nx] - mean) / sigma;
02454 }
02455 }
02456 }
02457
02458 image->update();
02459 }
02460
02461 void RealToFFTProcessor::process_inplace(EMData *image)
02462 {
02463 if (!image) {
02464 LOGWARN("NULL Image");
02465 return;
02466 }
02467
02468
02469 if(image->is_complex()) {
02470 LOGERR("%s Processor only operates on real images", get_name().c_str());
02471 throw ImageFormatException("apply to real image only");
02472 }
02473
02474
02475 int nz = image->get_zsize();
02476 if (nz > 1) {
02477 LOGERR("%s Processor doesn't support 3D models", get_name().c_str());
02478 throw ImageDimensionException("3D model not supported");
02479 }
02480
02481 EMData *ff=image->do_fft();
02482 ff->ri2ap();
02483
02484 int nx=image->get_xsize();
02485 int ny=image->get_ysize();
02486
02487 int x,y;
02488 float norm=static_cast<float>(nx*ny);
02489
02490 for (y=0; y<ny; y++) image->set_value_at(0,y,0);
02491
02492 for (x=1; x<nx/2; x++) {
02493 for (y=0; y<ny; y++) {
02494 int y2;
02495 if (y<ny/2) y2=y+ny/2;
02496 else if (y==ny/2) y2=ny;
02497 else y2=y-ny/2;
02498 image->set_value_at(x,y,ff->get_value_at(nx-x*2,ny-y2)/norm);
02499 }
02500 }
02501
02502 for (x=nx/2; x<nx; x++) {
02503 for (y=0; y<ny; y++) {
02504 int y2;
02505 if (y<ny/2) y2=y+ny/2;
02506 else y2=y-ny/2;
02507 image->set_value_at(x,y,ff->get_value_at(x*2-nx,y2)/norm);
02508 }
02509 }
02510
02511 image->update();
02512 if( ff )
02513 {
02514 delete ff;
02515 ff = 0;
02516 }
02517 }
02518
02519 void SigmaZeroEdgeProcessor::process_inplace(EMData * image)
02520 {
02521 if (!image) {
02522 LOGWARN("NULL Image");
02523 return;
02524 }
02525
02526 if (image->get_zsize() > 1) {
02527 LOGERR("%s Processor doesn't support 3D model", get_name().c_str());
02528 throw ImageDimensionException("3D model not supported");
02529 }
02530 float *d = image->get_data();
02531 int i = 0;
02532 int j = 0;
02533
02534 int nx = image->get_xsize();
02535 int ny = image->get_ysize();
02536
02537 for (j = 0; j < ny; j++) {
02538 for (i = 0; i < nx - 1; i++) {
02539 if (d[i + j * nx] != 0) {
02540 break;
02541 }
02542 }
02543
02544 float v = d[i + j * nx];
02545 while (i >= 0) {
02546 d[i + j * nx] = v;
02547 i--;
02548 }
02549
02550 for (i = nx - 1; i > 0; i--) {
02551 if (d[i + j * nx] != 0)
02552 break;
02553 }
02554 v = d[i + j * nx];
02555 while (i < nx) {
02556 d[i + j * nx] = v;
02557 i++;
02558 }
02559 }
02560
02561 for (i = 0; i < nx; i++) {
02562 for (j = 0; j < ny; j++) {
02563 if (d[i + j * nx] != 0)
02564 break;
02565 }
02566
02567 float v = d[i + j * nx];
02568 while (j >= 0) {
02569 d[i + j * nx] = v;
02570 j--;
02571 }
02572
02573 for (j = ny - 1; j > 0; j--) {
02574 if (d[i + j * nx] != 0)
02575 break;
02576 }
02577 v = d[i + j * nx];
02578 while (j < ny) {
02579 d[i + j * nx] = v;
02580 j++;
02581 }
02582 }
02583
02584
02585 image->update();
02586 }
02587
02588
02589
02590 void BeamstopProcessor::process_inplace(EMData * image)
02591 {
02592 if (!image) {
02593 LOGWARN("NULL Image");
02594 return;
02595 }
02596 if (image->get_zsize() > 1) {
02597 LOGERR("BeamstopProcessor doesn't support 3D model");
02598 throw ImageDimensionException("3D model not supported");
02599 }
02600
02601 float value1 = params["value1"];
02602 float value2 = params["value2"];
02603 float value3 = params["value3"];
02604
02605 float thr = fabs(value1);
02606 float *data = image->get_data();
02607 int cenx = (int) value2;
02608 int ceny = (int) value3;
02609
02610 int nx = image->get_xsize();
02611 int ny = image->get_ysize();
02612
02613 if (cenx <= 0) {
02614 cenx = nx / 2;
02615 }
02616
02617 if (ceny <= 0) {
02618 ceny = ny / 2;
02619 }
02620
02621 int mxr = (int) floor(sqrt(2.0f) * nx / 2);
02622
02623 float *mean_values = new float[mxr];
02624 float *sigma_values = new float[mxr];
02625 double sum = 0;
02626 int count = 0;
02627 double square_sum = 0;
02628
02629 for (int i = 0; i < mxr; i++) {
02630 sum = 0;
02631 count = 0;
02632 square_sum = 0;
02633 int nitems = 6 * i + 2;
02634
02635 for (int j = 0; j < nitems; j++) {
02636 float ang = j * 2 * M_PI / nitems;
02637 int x0 = (int) floor(cos(ang) * i + cenx);
02638 int y0 = (int) floor(sin(ang) * i + ceny);
02639
02640 if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny) {
02641 continue;
02642 }
02643
02644 float f = data[x0 + y0 * nx];
02645 sum += f;
02646 square_sum += f * f;
02647 count++;
02648 }
02649
02650 mean_values[i] = (float)sum / count;
02651 sigma_values[i] = (float) sqrt(square_sum / count - mean_values[i] * mean_values[i]);
02652 }
02653
02654
02655 for (int k = 0; k < 5; k++) {
02656 for (int i = 0; i < mxr; i++) {
02657 sum = 0;
02658 count = 0;
02659 square_sum = 0;
02660 int nitems = 6 * i + 2;
02661 double thr1 = mean_values[i] - sigma_values[i] * thr;
02662 double thr2 = mean_values[i] + sigma_values[i];
02663
02664 for (int j = 0; j < nitems; j++) {
02665 float ang = j * 2 * M_PI / nitems;
02666 int x0 = (int) floor(cos(ang) * i + cenx);
02667 int y0 = (int) floor(sin(ang) * i + ceny);
02668
02669 if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny ||
02670 data[x0 + y0 * nx] < thr1 || data[x0 + y0 * nx] > thr2) {
02671 continue;
02672 }
02673
02674 sum += data[x0 + y0 * nx];
02675 square_sum += data[x0 + y0 * nx] * data[x0 + y0 * nx];
02676 count++;
02677 }
02678
02679 mean_values[i] = (float) sum / count;
02680 sigma_values[i] = (float) sqrt(square_sum / count - mean_values[i] * mean_values[i]);
02681 }
02682 }
02683
02684 for (int i = 0; i < nx; i++) {
02685 for (int j = 0; j < ny; j++) {
02686
02687 #ifdef _WIN32
02688 int r = Util::round(_hypot((float) i - cenx, (float) j - ceny));
02689 #else
02690 int r = Util::round(hypot((float) i - cenx, (float) j - ceny));
02691 #endif //_WIN32
02692
02693 if (value1 < 0) {
02694 if (data[i + j * nx] < (mean_values[r] - sigma_values[r] * thr)) {
02695 data[i + j * nx] = 0;
02696 }
02697 else {
02698 data[i + j * nx] -= mean_values[r];
02699 }
02700 continue;
02701 }
02702 if (data[i + j * nx] > (mean_values[r] - sigma_values[r] * thr)) {
02703 continue;
02704 }
02705 data[i + j * nx] = mean_values[r];
02706 }
02707 }
02708
02709 if( mean_values )
02710 {
02711 delete[]mean_values;
02712 mean_values = 0;
02713 }
02714
02715 if( sigma_values )
02716 {
02717 delete[]sigma_values;
02718 sigma_values = 0;
02719 }
02720
02721 image->update();
02722 }
02723
02724
02725
02726 void MeanZeroEdgeProcessor::process_inplace(EMData * image)
02727 {
02728 if (!image) {
02729 LOGWARN("NULL Image");
02730 return;
02731 }
02732 if (image->get_zsize() > 1) {
02733 LOGERR("MeanZeroEdgeProcessor doesn't support 3D model");
02734 throw ImageDimensionException("3D model not supported");
02735 }
02736
02737 int nx = image->get_xsize();
02738 int ny = image->get_ysize();
02739 Dict dict = image->get_attr_dict();
02740 float mean_nonzero = dict.get("mean_nonzero");
02741
02742 float *d = image->get_data();
02743 int i = 0;
02744 int j = 0;
02745
02746 for (j = 0; j < ny; j++) {
02747 for (i = 0; i < nx - 1; i++) {
02748 if (d[i + j * nx] != 0) {
02749 break;
02750 }
02751 }
02752
02753 if (i == nx - 1) {
02754 i = -1;
02755 }
02756
02757 float v = d[i + j * nx] - mean_nonzero;
02758
02759 while (i >= 0) {
02760 v *= 0.9f;
02761 d[i + j * nx] = v + mean_nonzero;
02762 i--;
02763 }
02764
02765
02766 for (i = nx - 1; i > 0; i--) {
02767 if (d[i + j * nx] != 0) {
02768 break;
02769 }
02770 }
02771
02772 if (i == 0) {
02773 i = nx;
02774 }
02775
02776 v = d[i + j * nx] - mean_nonzero;
02777
02778 while (i < nx) {
02779 v *= .9f;
02780 d[i + j * nx] = v + mean_nonzero;
02781 i++;
02782 }
02783 }
02784
02785
02786 for (i = 0; i < nx; i++) {
02787 for (j = 0; j < ny; j++) {
02788 if (d[i + j * nx] != 0)
02789 break;
02790 }
02791
02792 float v = d[i + j * nx] - mean_nonzero;
02793
02794 while (j >= 0) {
02795 v *= .9f;
02796 d[i + j * nx] = v + mean_nonzero;
02797 j--;
02798 }
02799
02800 for (j = ny - 1; j > 0; j--) {
02801 if (d[i + j * nx] != 0)
02802 break;
02803 }
02804
02805 v = d[i + j * nx] - mean_nonzero;
02806
02807 while (j < ny) {
02808 v *= .9f;
02809 d[i + j * nx] = v + mean_nonzero;
02810 j++;
02811 }
02812 }
02813
02814 image->update();
02815 }
02816
02817
02818
02819 void AverageXProcessor::process_inplace(EMData * image)
02820 {
02821 if (!image) {
02822 LOGWARN("NULL Image");
02823 return;
02824 }
02825
02826 float *data = image->get_data();
02827 int nx = image->get_xsize();
02828 int ny = image->get_ysize();
02829 int nz = image->get_zsize();
02830 size_t nxy = (size_t)nx * ny;
02831
02832 size_t idx;
02833 for (int z = 0; z < nz; z++) {
02834 for (int x = 0; x < nx; x++) {
02835 double sum = 0;
02836 for (int y = 0; y < ny; y++) {
02837 idx = x + y * nx + z * nxy;
02838 sum += data[idx];
02839 }
02840 float mean = (float) sum / ny;
02841
02842 for (int y = 0; y < ny; y++) {
02843 idx = x + y * nx + z * nxy;
02844 data[idx] = mean;
02845 }
02846 }
02847 }
02848
02849 image->update();
02850 }
02851
02852 void DecayEdgeProcessor::process_inplace(EMData * image)
02853 {
02854 if (!image) {
02855 LOGWARN("NULL Image");
02856 return;
02857 }
02858
02859 if (image->get_zsize() > 1) throw ImageDimensionException("3D model not supported");
02860
02861 int nx = image->get_xsize();
02862 int ny = image->get_ysize();
02863
02864 float *d = image->get_data();
02865 int width = params["width"];
02866
02867 for (int i=0; i<width; i++) {
02868 float frac=i/(float)width;
02869 for (int j=0; j<nx; j++) {
02870 d[j+i*nx]*=frac;
02871 d[nx*ny-j-i*nx-1]*=frac;
02872 }
02873 for (int j=0; j<ny; j++) {
02874 d[j*nx+i]*=frac;
02875 d[nx*ny-j*nx-i-1]*=frac;
02876 }
02877 }
02878
02879 image->update();
02880 }
02881
02882 void ZeroEdgeRowProcessor::process_inplace(EMData * image)
02883 {
02884 if (!image) {
02885 LOGWARN("NULL Image");
02886 return;
02887 }
02888
02889 if (image->get_zsize() > 1) {
02890 LOGERR("ZeroEdgeRowProcessor is not supported in 3D models");
02891 throw ImageDimensionException("3D model not supported");
02892 }
02893
02894 int nx = image->get_xsize();
02895 int ny = image->get_ysize();
02896
02897 float *d = image->get_data();
02898 int top_nrows = params["y0"];
02899 int bottom_nrows = params["y1"];
02900
02901 int left_ncols = params["x0"];
02902 int right_ncols = params["x1"];
02903
02904 size_t row_size = nx * sizeof(float);
02905
02906 memset(d, 0, top_nrows * row_size);
02907 memset(d + (ny - bottom_nrows) * nx, 0, bottom_nrows * row_size);
02908
02909 for (int i = top_nrows; i < ny - bottom_nrows; i++) {
02910 memset(d + i * nx, 0, left_ncols * sizeof(float));
02911 memset(d + i * nx + nx - right_ncols, 0, right_ncols * sizeof(float));
02912 }
02913 image->update();
02914 }
02915
02916 void ZeroEdgePlaneProcessor::process_inplace(EMData * image)
02917 {
02918 if (!image) {
02919 LOGWARN("NULL Image");
02920 return;
02921 }
02922
02923 if (image->get_zsize() <= 1) {
02924 LOGERR("ZeroEdgePlaneProcessor only support 3D models");
02925 throw ImageDimensionException("3D model only");
02926 }
02927
02928 int nx = image->get_xsize();
02929 int ny = image->get_ysize();
02930 int nz = image->get_zsize();
02931
02932 float *d = image->get_data();
02933
02934 int x0=params["x0"];
02935 int x1=params["x1"];
02936 int y0=params["y0"];
02937 int y1=params["y1"];
02938 int z0=params["z0"];
02939 int z1=params["z1"];
02940
02941 size_t row_size = nx * sizeof(float);
02942 size_t nxy = nx * ny;
02943 size_t sec_size = nxy * sizeof(float);
02944 size_t y0row = y0 * row_size;
02945 size_t y1row = y1 * row_size;
02946 int max_y = ny-y1;
02947 size_t x0size = x0*sizeof(float);
02948 size_t x1size = x1*sizeof(float);
02949
02950 memset(d,0,z0*sec_size);
02951 memset(d+(nxy*(nz-z1)),0,sec_size*z1);
02952
02953 for (int z=z0; z<nz-z1; z++) {
02954 memset(d+z*nxy,0,y0row);
02955 memset(d+z*nxy+(ny-y1)*nx,0,y1row);
02956
02957 int znxy = z * nxy;
02958 int znxy2 = znxy + nx - x1;
02959
02960 for (int y=y0; y<max_y; y++) {
02961 memset(d+znxy+y*nx,0,x0size);
02962 memset(d+znxy2+y*nx,0,x1size);
02963 }
02964 }
02965
02966 image->update();
02967 }
02968
02969
02970 float NormalizeProcessor::calc_sigma(EMData * image) const
02971 {
02972 return image->get_attr("sigma");
02973 }
02974
02975 void NormalizeProcessor::process_inplace(EMData * image)
02976 {
02977 if (!image) {
02978 LOGWARN("cannot do normalization on NULL image");
02979 return;
02980 }
02981
02982 if (image->is_complex()) {
02983 LOGWARN("cannot do normalization on complex image");
02984 return;
02985 }
02986
02987 float sigma = calc_sigma(image);
02988 if (sigma == 0 || !Util::goodf(&sigma)) {
02989 LOGWARN("cannot do normalization on image with sigma = 0");
02990 return;
02991 }
02992
02993 float mean = calc_mean(image);
02994
02995 size_t size = (size_t)image->get_xsize() * (size_t)image->get_ysize() *
02996 (size_t)image->get_zsize();
02997 float *data = image->get_data();
02998
02999 for (size_t i = 0; i < size; i++) {
03000 data[i] = (data[i] - mean) / sigma;
03001 }
03002
03003 image->update();
03004 }
03005
03006 float NormalizeUnitProcessor::calc_sigma(EMData * image) const
03007 {
03008 if (!image) {
03009 LOGWARN("NULL Image");
03010 return 0;
03011 }
03012 float ret=sqrt((float)image->get_attr("square_sum"));
03013 return ret==0.0f?1.0f:ret;
03014 }
03015
03016 float NormalizeUnitSumProcessor::calc_sigma(EMData * image) const
03017 {
03018 if (!image) {
03019 LOGWARN("NULL Image");
03020 return 0;
03021 }
03022 float ret=(float)image->get_attr("mean")*image->get_xsize()*image->get_ysize()*image->get_zsize();
03023 return ret==0.0f?1.0f:ret;
03024 }
03025
03026 float NormalizeMaskProcessor::calc_sigma(EMData * image) const
03027 {
03028 if (!image) {
03029 LOGWARN("NULL Image");
03030 return 0;
03031 }
03032 EMData *mask = params["mask"];
03033 int no_sigma = params["no_sigma"];
03034
03035 if(no_sigma == 0) {
03036 return 1;
03037 }
03038 else {
03039 if (!EMUtil::is_same_size(mask, image)) {
03040 LOGERR("normalize.maskProcessor: mask and image must be the same size");
03041 throw ImageDimensionException("mask and image must be the same size");
03042 }
03043
03044 float *data = image->get_data();
03045 float *mask_data = mask->get_data();
03046 size_t size = image->get_xsize() * image->get_ysize() * image->get_zsize();
03047 double sum = 0;
03048 double sq2 = 0;
03049 size_t n_norm = 0;
03050
03051 for (size_t i = 0; i < size; i++) {
03052 if (mask_data[i] > 0.5f) {
03053 sum += data[i];
03054 sq2 += data[i]*double (data[i]);
03055 n_norm++;
03056 }
03057 }
03058 return sqrt(static_cast<float>((sq2 - sum * sum /n_norm)/(n_norm -1))) ;
03059 }
03060 }
03061
03062 float NormalizeMaskProcessor::calc_mean(EMData * image) const
03063 {
03064 if (!image) {
03065 LOGWARN("NULL Image");
03066 return 0;
03067 }
03068 EMData *mask = params["mask"];
03069
03070 if (!EMUtil::is_same_size(mask, image)) {
03071 LOGERR("normalize.maskProcessor: mask and image must be the same size");
03072 throw ImageDimensionException("mask and image must be the same size");
03073 }
03074
03075 float *data = image->get_data();
03076 float *mask_data = mask->get_data();
03077 size_t size = image->get_xsize() * image->get_ysize() * image->get_zsize();
03078 double sum = 0;
03079 size_t n_norm = 0;
03080
03081 for (size_t i = 0; i < size; i++) {
03082 if (mask_data[i] > 0.5f) {
03083 sum += data[i];
03084 n_norm++;
03085 }
03086 }
03087
03088 float mean = 0;
03089 if (n_norm == 0) {
03090 mean = image->get_edge_mean();
03091 }
03092 else {
03093 mean = (float) sum / n_norm;
03094 }
03095
03096 return mean;
03097 }
03098
03099 void NormalizeRampNormVar::process_inplace(EMData * image)
03100 {
03101 if (!image) {
03102 LOGWARN("cannot do normalization on NULL image");
03103 return;
03104 }
03105
03106 if (image->is_complex()) {
03107 LOGWARN("cannot do normalization on complex image");
03108 return;
03109 }
03110
03111 image->process_inplace( "filter.ramp" );
03112 int nx = image->get_xsize();
03113 EMData mask(nx,nx);
03114 mask.process_inplace("testimage.circlesphere", Dict("radius",nx/2-2,"fill",1));
03115
03116 vector<float> rstls = Util::infomask( image, &mask, false);
03117 image->add((float)-rstls[0]);
03118 image->mult((float)1.0/rstls[1]);
03119 image->update();
03120 }
03121
03122 void NormalizeByMassProcessor::process_inplace(EMData * image)
03123 {
03124 float mass = params.set_default("mass",-1.0f);
03125
03126 if (mass <= 0) throw InvalidParameterException("You must specify a positive non zero mass");
03127
03128 float thr = params.set_default("thr",(float)image->get_attr("mean")+(float)image->get_attr("sigma"));
03129
03130 float apix = params.set_default("apix",-1.123456789f);
03131 if (apix == -1.123456789 ) {
03132 if (image->has_attr("apix_x")) {
03133 apix = image->get_attr("apix_x");
03134 }
03135 }
03136
03137 if (apix <= 0) throw InvalidParameterException("You must specify a positive non zero apix");
03138
03139 float step = ((float)image->get_attr("sigma"))/2.0f;
03140
03141 int count=0;
03142 size_t n = image->get_size();
03143 float* d = image->get_data();
03144
03145 for (size_t i=0; i<n; ++i) {
03146 if (d[i]>=thr) ++count;
03147 }
03148
03149 float max = image->get_attr("maximum");
03150 float min = image->get_attr("minimum");
03151 for (int j=0; j<4; j++) {
03152 while (thr<max && count*apix*apix*apix*.81/1000.0>mass) {
03153 thr+=step;
03154 count=0;
03155 for (size_t i=0; i<n; ++i) {
03156 if (d[i]>=thr) ++count;
03157 }
03158 }
03159
03160 step/=4.0;
03161
03162 while (thr>min && count*apix*apix*apix*.81/1000.0<mass) {
03163 thr-=step;
03164 count=0;
03165 for (size_t i=0; i<n; ++i) {
03166 if (d[i]>=thr) ++count;
03167 }
03168 }
03169
03170 step/=4.0;
03171 }
03172
03173 image->mult((float)1.0/thr);
03174 image->update();
03175 }
03176
03177 float NormalizeEdgeMeanProcessor::calc_mean(EMData * image) const
03178 {
03179 if (!image) {
03180 LOGWARN("NULL Image");
03181 return 0;
03182 }
03183 return image->get_edge_mean();
03184 }
03185
03186 float NormalizeCircleMeanProcessor::calc_mean(EMData * image) const
03187 {
03188 if (!image) {
03189 LOGWARN("NULL Image");
03190 return 0;
03191 }
03192 return image->get_circle_mean();
03193 }
03194
03195
03196 float NormalizeMaxMinProcessor::calc_sigma(EMData * image) const
03197 {
03198 if (!image) {
03199 LOGWARN("NULL Image");
03200 return 0;
03201 }
03202 float maxval = image->get_attr("maximum");
03203 float minval = image->get_attr("minimum");
03204 return (maxval + minval) / 2;
03205 }
03206
03207 float NormalizeMaxMinProcessor::calc_mean(EMData * image) const
03208 {
03209 if (!image) {
03210 LOGWARN("NULL Image");
03211 return 0;
03212 }
03213 float maxval = image->get_attr("maximum");
03214 float minval = image->get_attr("minimum");
03215 return (maxval - minval) / 2;
03216 }
03217
03218 float NormalizeLREdgeMeanProcessor::calc_mean(EMData * image) const
03219 {
03220 if (!image) {
03221 LOGWARN("NULL Image");
03222 return 0;
03223 }
03224 double sum = 0;
03225 int nx = image->get_xsize();
03226 int ny = image->get_ysize();
03227 int nz = image->get_zsize();
03228 float *d = image->get_data();
03229 int nyz = ny * nz;
03230
03231 for (int i = 0; i < nyz; i++) {
03232 size_t l = i * nx;
03233 size_t r = l + nx - 2;
03234 sum += d[l] + d[l + 1] + d[r] + d[r + 1];
03235 }
03236 float mean = (float) sum / (4 * nyz);
03237 return mean;
03238 }
03239
03240 void NormalizeRowProcessor::process_inplace(EMData * image)
03241 {
03242 if (!image) {
03243 LOGWARN("NULL Image");
03244 return;
03245 }
03246
03247 if (image->get_zsize() > 1) {
03248 LOGERR("row normalize only works for 2D image");
03249 return;
03250 }
03251
03252 float *rdata = image->get_data();
03253 int nx = image->get_xsize();
03254 int ny = image->get_ysize();
03255
03256 for (int y = 0; y < ny; y++) {
03257 double row_sum = 0;
03258 for (int x = 0; x < nx; x++) {
03259 row_sum += rdata[x + y * nx];
03260 }
03261
03262 double row_mean = row_sum / nx;
03263 if (row_mean <= 0) {
03264 row_mean = 1;
03265 }
03266
03267 for (int x = 0; x < nx; x++) {
03268 rdata[x + y * nx] /= (float)row_mean;
03269 }
03270 }
03271
03272 image->update();
03273 }
03274
03275 float NormalizeStdProcessor::calc_mean(EMData * image) const
03276 {
03277 if (!image) {
03278 LOGWARN("NULL Image");
03279 return 0;
03280 }
03281 return image->get_attr("mean");
03282 }
03283
03284
03285
03286 void NormalizeToLeastSquareProcessor::process_inplace(EMData * image)
03287 {
03288 if (!image) {
03289 LOGWARN("NULL Image");
03290 return;
03291 }
03292
03293 EMData *to = params["to"];
03294
03295 float low_threshold = FLT_MIN;
03296 string low_thr_name = "low_threshold";
03297 if (params.has_key(low_thr_name)) {
03298 low_threshold = params[low_thr_name];
03299 }
03300
03301 float high_threshold = FLT_MAX;
03302 string high_thr_name = "high_threshold";
03303 if (params.has_key(high_thr_name)) {
03304 high_threshold = params[high_thr_name];
03305 }
03306
03307 float *rawp = image->get_data();
03308 float *refp = to->get_data();
03309
03310 int nx = image->get_xsize();
03311 int ny = image->get_ysize();
03312 int nz = image->get_zsize();
03313 size_t size = nx * ny * nz;
03314
03315 float sum_x = 0;
03316 float sum_y = 0;
03317 int count = 0;
03318
03319 for (size_t i = 0; i < size; i++) {
03320 if (refp[i] >= low_threshold && refp[i] <= high_threshold && refp[i] != 0.0 && rawp[i] != 0.0) {
03321 count++;
03322 sum_x += refp[i];
03323 sum_y += rawp[i];
03324 }
03325 }
03326
03327 float sum_x_mean = sum_x / count;
03328 float sum_tt = 0;
03329 float b = 0;
03330
03331 float t;
03332 for (size_t i = 0; i < size; i++) {
03333 if (refp[i] >= low_threshold && refp[i] <= high_threshold && refp[i] != 0.0 && rawp[i] != 0.0) {
03334 t = refp[i] - sum_x_mean;
03335 sum_tt += t * t;
03336 b += t * rawp[i];
03337 }
03338 }
03339
03340 b /= sum_tt;
03341
03342 float a = (sum_y - sum_x * b) / count;
03343 float scale = 1 / b;
03344 float shift = -a / b;
03345
03346 for (size_t i = 0; i < size; i++) {
03347 rawp[i] = (rawp[i] - a) / b;
03348 }
03349
03350 image->update();
03351
03352 params["scale"] = scale;
03353 params["shift"] = shift;
03354
03355 image->set_attr("norm_mult",scale);
03356 image->set_attr("norm_add",shift);
03357
03358 }
03359
03360
03361 void BinarizeFourierProcessor::process_inplace(EMData* image) {
03362 ENTERFUNC;
03363 if (!image->is_complex()) throw ImageFormatException("Fourier binary thresholding processor only works for complex images");
03364
03365 float threshold = params.set_default("value",-1.0f);
03366 if (threshold < 0) throw InvalidParameterException("For fourier amplitude-based thresholding, the threshold must be greater than or equal to 0.");
03367
03368 image->ri2ap();
03369
03370 #ifdef EMAN2_USING_CUDA
03371 if (image->gpu_operation_preferred()) {
03372 EMDataForCuda tmp = image->get_data_struct_for_cuda();
03373 binarize_fourier_amp_processor(&tmp,threshold);
03374 image->set_ri(true);
03375 image->gpu_update();
03376 EXITFUNC;
03377 return;
03378 }
03379 #endif
03380
03381 float* d = image->get_data();
03382 for( size_t i = 0; i < image->get_size()/2; ++i, d+=2) {
03383 float v = *d;
03384 if ( v >= threshold ) {
03385 *d = 1;
03386 *(d+1) = 0;
03387 } else {
03388 *d = 0;
03389 *(d+1) = 0;
03390 }
03391 }
03392
03393
03394 image->set_ri(true);
03395 image->update();
03396 EXITFUNC;
03397 }
03398
03399
03400 void BilateralProcessor::process_inplace(EMData * image)
03401 {
03402 if (!image) {
03403 LOGWARN("NULL Image");
03404 return;
03405 }
03406
03407 float distance_sigma = params["distance_sigma"];
03408 float value_sigma = params["value_sigma"];
03409 int max_iter = params["niter"];
03410 int half_width = params["half_width"];
03411
03412 if (half_width < distance_sigma) {
03413 LOGWARN("localwidth(=%d) should be larger than distance_sigma=(%f)\n",
03414 half_width, distance_sigma);
03415 }
03416
03417 distance_sigma *= distance_sigma;
03418
03419 float image_sigma = image->get_attr("sigma");
03420 if (image_sigma > value_sigma) {
03421 LOGWARN("image sigma(=%f) should be smaller than value_sigma=(%f)\n",
03422 image_sigma, value_sigma);
03423 }
03424 value_sigma *= value_sigma;
03425
03426 int nx = image->get_xsize();
03427 int ny = image->get_ysize();
03428 int nz = image->get_zsize();
03429
03430 if(nz==1) {
03431 int width=nx, height=ny;
03432
03433 int i,j,m,n;
03434
03435 float tempfloat1,tempfloat2,tempfloat3;
03436 int index1,index2,index;
03437 int Iter;
03438 int tempint1,tempint3;
03439
03440 tempint1=width;
03441 tempint3=width+2*half_width;
03442
03443 float* mask=(float*)calloc((2*half_width+1)*(2*half_width+1),sizeof(float));
03444 float* OrgImg=(float*)calloc((2*half_width+width)*(2*half_width+height),sizeof(float));
03445 float* NewImg=image->get_data();
03446
03447 for(m=-(half_width);m<=half_width;m++)
03448 for(n=-(half_width);n<=half_width;n++) {
03449 index=(m+half_width)*(2*half_width+1)+(n+half_width);
03450 mask[index]=exp((float)(-(m*m+n*n)/distance_sigma/2.0));
03451 }
03452
03453
03454
03455 Iter=0;
03456 while(Iter<max_iter) {
03457 for(i=0;i<height;i++)
03458 for(j=0;j<width;j++) {
03459 index1=(i+half_width)*tempint3+(j+half_width);
03460 index2=i*tempint1+j;
03461 OrgImg[index1]=NewImg[index2];
03462 }
03463
03464
03465 for(i=0;i<height;i++){
03466 for(j=0;j<half_width;j++) OrgImg[(i+half_width)*tempint3+(j)]=OrgImg[(i+half_width)*tempint3+(2*half_width-j)];
03467 for(j=0;j<half_width;j++) OrgImg[(i+half_width)*tempint3+(j+width+half_width)]=OrgImg[(i+half_width)*tempint3+(width+half_width-j-2)];
03468 }
03469 for(i=0;i<half_width;i++){
03470 for(j=0;j<(width+2*half_width);j++) OrgImg[i*tempint3+j]=OrgImg[(2*half_width-i)*tempint3+j];
03471 for(j=0;j<(width+2*half_width);j++) OrgImg[(i+height+half_width)*tempint3+j]=OrgImg[(height+half_width-2-i)*tempint3+j];
03472 }
03473
03474
03475
03476
03477 for(i=0;i<height;i++){
03478
03479 for(j=0;j<width;j++){
03480 tempfloat1=0.0; tempfloat2=0.0;
03481 for(m=-(half_width);m<=half_width;m++)
03482 for(n=-(half_width);n<=half_width;n++){
03483 index =(m+half_width)*(2*half_width+1)+(n+half_width);
03484 index1=(i+half_width)*tempint3+(j+half_width);
03485 index2=(i+half_width+m)*tempint3+(j+half_width+n);
03486 tempfloat3=(OrgImg[index1]-OrgImg[index2])*(OrgImg[index1]-OrgImg[index2]);
03487
03488 tempfloat3=mask[index]*(1.0f/(1+tempfloat3/value_sigma));
03489
03490 tempfloat1+=tempfloat3;
03491
03492 tempfloat2+=tempfloat3*OrgImg[(i+half_width+m)*tempint3+(j+half_width+n)];
03493 }
03494 NewImg[i*width+j]=tempfloat2/tempfloat1;
03495 }
03496 }
03497 Iter++;
03498 }
03499
03500
03501
03502 free(mask);
03503 free(OrgImg);
03504
03505
03506 }
03507 else {
03508 int width = nx;
03509 int height = ny;
03510 int slicenum = nz;
03511
03512 int slice_size = width * height;
03513 int new_width = width + 2 * half_width;
03514 int new_slice_size = (width + 2 * half_width) * (height + 2 * half_width);
03515
03516 int width1 = 2 * half_width + 1;
03517 int mask_size = width1 * width1;
03518 int old_img_size = (2 * half_width + width) * (2 * half_width + height);
03519
03520 int zstart = -half_width;
03521 int zend = -half_width;
03522 int is_3d = 0;
03523 if (nz > 1) {
03524 mask_size *= width1;
03525 old_img_size *= (2 * half_width + slicenum);
03526 zend = half_width;
03527 is_3d = 1;
03528 }
03529
03530 float *mask = (float *) calloc(mask_size, sizeof(float));
03531 float *old_img = (float *) calloc(old_img_size, sizeof(float));
03532
03533 float *new_img = image->get_data();
03534
03535 for (int p = zstart; p <= zend; p++) {
03536 int cur_p = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1);
03537
03538 for (int m = -half_width; m <= half_width; m++) {
03539 int cur_m = (m + half_width) * (2 * half_width + 1) + half_width;
03540
03541 for (int n = -half_width; n <= half_width; n++) {
03542 int l = cur_p + cur_m + n;
03543 mask[l] = exp((float) (-(m * m + n * n + p * p * is_3d) / distance_sigma / 2.0f));
03544 }
03545 }
03546 }
03547
03548 int iter = 0;
03549 while (iter < max_iter) {
03550 for (int k = 0; k < slicenum; k++) {
03551 int cur_k1 = (k + half_width) * new_slice_size * is_3d;
03552 int cur_k2 = k * slice_size;
03553
03554 for (int i = 0; i < height; i++) {
03555 int cur_i1 = (i + half_width) * new_width;
03556 int cur_i2 = i * width;
03557
03558 for (int j = 0; j < width; j++) {
03559 int k1 = cur_k1 + cur_i1 + (j + half_width);
03560 int k2 = cur_k2 + cur_i2 + j;
03561 old_img[k1] = new_img[k2];
03562 }
03563 }
03564 }
03565
03566 for (int k = 0; k < slicenum; k++) {
03567 int cur_k = (k + half_width) * new_slice_size * is_3d;
03568
03569 for (int i = 0; i < height; i++) {
03570 int cur_i = (i + half_width) * new_width;
03571
03572 for (int j = 0; j < half_width; j++) {
03573 int k1 = cur_k + cur_i + j;
03574 int k2 = cur_k + cur_i + (2 * half_width - j);
03575 old_img[k1] = old_img[k2];
03576 }
03577
03578 for (int j = 0; j < half_width; j++) {
03579 int k1 = cur_k + cur_i + (width + half_width + j);
03580 int k2 = cur_k + cur_i + (width + half_width - j - 2);
03581 old_img[k1] = old_img[k2];
03582 }
03583 }
03584
03585
03586 for (int i = 0; i < half_width; i++) {
03587 int i2 = i * new_width;
03588 int i3 = (2 * half_width - i) * new_width;
03589 for (int j = 0; j < (width + 2 * half_width); j++) {
03590 int k1 = cur_k + i2 + j;
03591 int k2 = cur_k + i3 + j;
03592 old_img[k1] = old_img[k2];
03593 }
03594
03595 i2 = (height + half_width + i) * new_width;
03596 i3 = (height + half_width - 2 - i) * new_width;
03597 for (int j = 0; j < (width + 2 * half_width); j++) {
03598 int k1 = cur_k + i2 + j;
03599 int k2 = cur_k + i3 + j;
03600 old_img[k1] = old_img[k2];
03601 }
03602 }
03603 }
03604
03605 size_t idx;
03606 for (int k = 0; k < slicenum; k++) {
03607 int cur_k = (k + half_width) * new_slice_size;
03608
03609 for (int i = 0; i < height; i++) {
03610 int cur_i = (i + half_width) * new_width;
03611
03612 for (int j = 0; j < width; j++) {
03613 float f1 = 0;
03614 float f2 = 0;
03615 int k1 = cur_k + cur_i + (j + half_width);
03616
03617 for (int p = zstart; p <= zend; p++) {
03618 int cur_p1 = (p + half_width) * (2 * half_width + 1) * (2 * half_width + 1);
03619 int cur_p2 = (k + half_width + p) * new_slice_size;
03620
03621 for (int m = -half_width; m <= half_width; m++) {
03622 int cur_m1 = (m + half_width) * (2 * half_width + 1);
03623 int cur_m2 = cur_p2 + cur_i + m * new_width + j + half_width;
03624
03625 for (int n = -half_width; n <= half_width; n++) {
03626 int k = cur_p1 + cur_m1 + (n + half_width);
03627 int k2 = cur_m2 + n;
03628 float f3 = Util::square(old_img[k1] - old_img[k2]);
03629
03630 f3 = mask[k] * (1.0f / (1 + f3 / value_sigma));
03631 f1 += f3;
03632 int l1 = cur_m2 + n;
03633 f2 += f3 * old_img[l1];
03634 }
03635
03636 idx = k * height * width + i * width + j;
03637 new_img[idx] = f2 / f1;
03638 }
03639 }
03640 }
03641 }
03642 }
03643 iter++;
03644 }
03645 if( mask ) {
03646 free(mask);
03647 mask = 0;
03648 }
03649
03650 if( old_img ) {
03651 free(old_img);
03652 old_img = 0;
03653 }
03654 }
03655
03656 image->update();
03657 }
03658
03659 void RotationalAverageProcessor::process_inplace(EMData * image)
03660 {
03661 if (!image || image->is_complex()) {
03662 LOGWARN("only works on real image. do nothing.");
03663 return;
03664 }
03665
03666 if (image->get_ndim() <= 0 || image->get_ndim() > 3) throw ImageDimensionException("radial average processor only works for 2D and 3D images");
03667
03668 float *rdata = image->get_data();
03669 int nx = image->get_xsize();
03670 int ny = image->get_ysize();
03671
03672 vector < float >dist = image->calc_radial_dist(nx / 2, 0, 1,0);
03673
03674 float midx = (float)((int)nx/2);
03675 float midy = (float)((int)ny/2);
03676
03677 size_t c = 0;
03678 if (image->get_ndim() == 2) {
03679 for (int y = 0; y < ny; y++) {
03680 for (int x = 0; x < nx; x++, c++) {
03681 #ifdef _WIN32
03682 float r = (float) _hypot(x - midx, y - midy);
03683 #else
03684 float r = (float) hypot(x - midx, y - midy);
03685 #endif //_WIN32
03686
03687
03688 int i = (int) floor(r);
03689 r -= i;
03690 if (i >= 0 && i < nx / 2 - 1) {
03691 rdata[c] = dist[i] * (1.0f - r) + dist[i + 1] * r;
03692 }
03693 else if (i < 0) {
03694 rdata[c] = dist[0];
03695 }
03696 else {
03697 rdata[c] = 0;
03698 }
03699 }
03700 }
03701 }
03702 else if (image->get_ndim() == 3) {
03703 int nz = image->get_zsize();
03704 float midz = (float)((int)nz/2);
03705 float r;
03706 int i;
03707 for (int z = 0; z < nz; ++z) {
03708 for (int y = 0; y < ny; ++y) {
03709 for (int x = 0; x < nx; ++x, ++c) {
03710
03711 r = (float) Util::hypot3(x - midx, y - midy, z - midz);
03712
03713 i = Util::fast_floor(r);
03714 r -= i;
03715 if (i >= 0 && i < nx / 2 - 1) {
03716 rdata[c] = dist[i] * (1.0f - r) + dist[i + 1] * r;
03717 }
03718 else if (i < 0) {
03719 rdata[c] = dist[0];
03720 }
03721 else {
03722 rdata[c] = 0;
03723 }
03724 }
03725 }
03726 }
03727 }
03728
03729 image->update();
03730 }
03731
03732
03733
03734 void RotationalSubstractProcessor::process_inplace(EMData * image)
03735 {
03736 if (!image || image->is_complex()) {
03737 LOGWARN("only works on real image. do nothing.");
03738 return;
03739 }
03740
03741 if (image->get_ndim() != 2) throw ImageDimensionException("This processor works only for 2D images");
03742
03743 float *rdata = image->get_data();
03744 int nx = image->get_xsize();
03745 int ny = image->get_ysize();
03746
03747 vector < float >dist = image->calc_radial_dist(nx / 2, 0, 1,0);
03748
03749 int c = 0;
03750 for (int y = 0; y < ny; y++) {
03751 for (int x = 0; x < nx; x++, c++) {
03752 #ifdef _WIN32
03753 float r = (float) _hypot(x - nx / 2, y - ny / 2);
03754 #else
03755 float r = (float) hypot(x - nx / 2, y - ny / 2);
03756 #endif
03757 int i = (int) floor(r);
03758 r -= i;
03759 if (i >= 0 && i < nx / 2 - 1) {
03760 rdata[c] -= dist[i] * (1.0f - r) + dist[i + 1] * r;
03761 }
03762 else {
03763 rdata[c] = 0;
03764 }
03765 }
03766 }
03767
03768 image->update();
03769 }
03770
03771
03772 EMData* TransposeProcessor::process(const EMData* const image) {
03773 if (image->get_ndim() != 2) throw UnexpectedBehaviorException("Transpose processor only works with 2D images");
03774 if (image->is_complex()) throw UnexpectedBehaviorException("Transpose processor only works with real images");
03775
03776 EMData* ret = new EMData(image->get_ysize(),image->get_xsize(),1);
03777
03778 for(int j = 0; j< image->get_ysize();++j) {
03779 for(int i = 0; i< image->get_xsize();++i) {
03780 ret->set_value_at(j,i,image->get_value_at(i,j));
03781 }
03782 }
03783
03784 return ret;
03785
03786 }
03787
03788 void TransposeProcessor::process_inplace(EMData* image) {
03789 if (image->get_ndim() != 2) throw UnexpectedBehaviorException("Transpose processor only works with 2D images");
03790 if (image->is_complex()) throw UnexpectedBehaviorException("Transpose processor only works with real images");
03791
03792 float* data = (float*)malloc(image->get_ysize()*image->get_xsize()*sizeof(float));
03793
03794 int nx = image->get_ysize();
03795 for(int j = 0; j< image->get_ysize();++j) {
03796 for(int i = 0; i< image->get_xsize();++i) {
03797 data[i*nx+j] = image->get_value_at(i,j);
03798 }
03799 }
03800
03801 image->set_data(data,image->get_ysize(),image->get_xsize(),1);
03802
03803 }
03804
03805 void FlipProcessor::process_inplace(EMData * image)
03806 {
03807 ENTERFUNC;
03808 if (!image) {
03809 LOGWARN("NULL Image");
03810 return;
03811 }
03812 string axis = (const char*)params["axis"];
03813 #ifdef EMAN2_USING_CUDA
03814 if (image->gpu_operation_preferred()) {
03815 float array[12] = {1.0, 0.0, 0.0, 0.0,
03816 0.0, 1.0, 0.0, 0.0,
03817 0.0, 0.0, 1.0, 0.0};
03818 if (axis == "x" || axis == "X") {
03819 array[0] = -1.0;
03820 }else if (axis == "y" || axis == "Y") {
03821 array[5] = -1.0;
03822 }
03823 else if (axis == "z" || axis == "Z") {
03824 array[10] = -1.0;
03825 }
03826 Transform t(array);
03827 Dict params("transform",(Transform*)&t);
03828 image->process_inplace("math.transform",params);
03829 EXITFUNC;
03830 return;
03831 }
03832 #endif
03833
03834
03835 float *d = image->get_data();
03836 int nx = image->get_xsize();
03837 int ny = image->get_ysize();
03838 int nz = image->get_zsize();
03839
03840 size_t nxy = nx * ny;
03841
03842
03843
03844
03845
03846
03847 if (axis == "x" || axis == "X") {
03848 int offset = (nx%2 == 0);
03849 size_t idx1, idx2;
03850 for(int z = 0; z < nz; ++z) {
03851 for(int y = 0; y < ny; ++y) {
03852 if (offset != 0 ) {
03853 idx1 = z*nxy + y*nx;
03854 d[idx1] = 0;
03855 }
03856 for(int x = offset; x < nx / 2; ++x) {
03857 idx1 = z*nxy + y*nx + x;
03858 idx2 = z*nxy + y*nx + (nx-x-1+offset);
03859 std::swap(d[idx1], d[idx2]);
03860 }
03861
03862 }
03863 }
03864 }
03865
03866 else if (axis == "y" || axis == "Y") {
03867 int offset = (ny%2 == 0);
03868 for(int z=0; z<nz; ++z) {
03869 if (offset != 0) {
03870 std::fill(d+z*nxy,d+z*nxy+nx,0);
03871 }
03872 for(int y=offset; y<ny/2; ++y) {
03873 for(int x=0; x<nx; ++x) {
03874 std::swap(d[z*nxy + y*nx +x], d[z*nxy + (ny -y -1+offset)*nx +x]);
03875 }
03876 }
03877 }
03878 }
03879 else if (axis == "z" || axis == "Z") {
03880 int offset = (nz%2 == 0);
03881 if (offset != 0) {
03882 std::fill(d,d+nxy,0);
03883 }
03884 size_t idx1, idx2;
03885 for(int z=offset; z<nz/2; ++z) {
03886 for(int y=0; y<ny; ++y) {
03887 for(int x=0; x<nx; ++x) {
03888 idx1 = z*nxy + y*nx + x;
03889 idx2 = (nz-z-1+offset)*nxy + y*nx + x;
03890 std::swap(d[idx1], d[idx2]);
03891 }
03892 }
03893 }
03894 }
03895
03896 image->update();
03897 EXITFUNC;
03898 }
03899
03900 void AddNoiseProcessor::process_inplace(EMData * image)
03901 {
03902 if (!image) {
03903 LOGWARN("NULL Image");
03904 return;
03905 }
03906
03907 Randnum * randnum = Randnum::Instance();
03908 if(params.has_key("seed")) {
03909 randnum->set_seed((int)params["seed"]);
03910 }
03911
03912 float addnoise = params["noise"];
03913 addnoise *= get_sigma(image);
03914 float *dat = image->get_data();
03915
03916 for (size_t j = 0; j < image->get_size(); ++j) {
03917 dat[j] += randnum->get_gauss_rand(addnoise, addnoise / 2);
03918 }
03919
03920 image->update();
03921 }
03922
03923 float AddSigmaNoiseProcessor::get_sigma(EMData * image)
03924 {
03925 if (!image) {
03926 LOGWARN("NULL Image");
03927 return 0;
03928 }
03929 return image->get_attr("sigma");
03930 }
03931
03932 void FourierToCornerProcessor::process_inplace(EMData * image)
03933 {
03934 if ( !image->is_complex() ) throw ImageFormatException("Can not Fourier origin shift an image that is not complex");
03935
03936 int nx=image->get_xsize();
03937 int ny=image->get_ysize();
03938 int nz=image->get_zsize();
03939
03940 int nxy = nx*ny;
03941
03942 if ( ny == 1 && nz == 1 ){
03943 cout << "Warning- attempted Fourier origin shift a 1D image - no action taken" << endl;
03944 return;
03945 }
03946 int yodd = (ny%2==1);
03947 int zodd = (nz%2==1);
03948
03949 float* rdata = image->get_data();
03950
03951 float tmp[2];
03952 float* p1;
03953 float* p2;
03954
03955 if (yodd){
03956
03957
03958
03959
03960 float prev[2];
03961 size_t idx;
03962 for( int s = 0; s < nz; s++ ) {
03963 for( int c =0; c < nx; c += 2 ) {
03964 idx = s*nxy+ny/2*nx+c;
03965 prev[0] = rdata[idx];
03966 prev[1] = rdata[idx+1];
03967 for( int r = 0; r <= ny/2; ++r ) {
03968 idx = s*nxy+r*nx+c;
03969 float* p1 = &rdata[idx];
03970 tmp[0] = p1[0];
03971 tmp[1] = p1[1];
03972
03973 p1[0] = prev[0];
03974 p1[1] = prev[1];
03975
03976 prev[0] = tmp[0];
03977 prev[1] = tmp[1];
03978 }
03979 }
03980 }
03981 }
03982
03983
03984 size_t idx1, idx2;
03985 for( int s = 0; s < nz; ++s ) {
03986 for( int r = 0 + yodd; r < ny/2+yodd; ++r ) {
03987 for( int c =0; c < nx; c += 2 ) {
03988 idx1 = s*nxy+r*nx+c;
03989 idx2 = s*nxy+(r+ny/2)*nx+c;
03990 p1 = &rdata[idx1];
03991 p2 = &rdata[idx2];
03992
03993 tmp[0] = p1[0];
03994 tmp[1] = p1[1];
03995
03996 p1[0] = p2[0];
03997 p1[1] = p2[1];
03998
03999 p2[0] = tmp[0];
04000 p2[1] = tmp[1];
04001 }
04002 }
04003 }
04004
04005 if ( nz != 1 )
04006 {
04007
04008 if (zodd){
04009
04010
04011
04012 float prev[2];
04013 size_t idx;
04014 for( int r = 0; r < ny; ++r ) {
04015 for( int c =0; c < nx; c += 2 ) {
04016 idx = nz/2*nxy+r*nx+c;
04017 prev[0] = rdata[idx];
04018 prev[1] = rdata[idx+1];
04019 for( int s = 0; s <= nz/2; ++s ) {
04020 idx = s*nxy+r*nx+c;
04021 float* p1 = &rdata[idx];
04022 tmp[0] = p1[0];
04023 tmp[1] = p1[1];
04024
04025 p1[0] = prev[0];
04026 p1[1] = prev[1];
04027
04028 prev[0] = tmp[0];
04029 prev[1] = tmp[1];
04030 }
04031 }
04032 }
04033 }
04034
04035
04036 size_t idx1, idx2;
04037 for( int s = 0+zodd; s < nz/2 + zodd; ++s ) {
04038 for( int r = 0; r < ny; ++r ) {
04039 for( int c =0; c < nx; c += 2 ) {
04040 idx1 = s*nxy+r*nx+c;
04041 idx2 = (s+nz/2)*nxy+r*nx+c;
04042 p1 = &rdata[idx1];
04043 p2 = &rdata[idx2];
04044
04045 tmp[0] = p1[0];
04046 tmp[1] = p1[1];
04047
04048 p1[0] = p2[0];
04049 p1[1] = p2[1];
04050
04051 p2[0] = tmp[0];
04052 p2[1] = tmp[1];
04053 }
04054 }
04