33#include "sparx/processor_sparx.h"
47#include <gsl/gsl_randist.h>
48#include <gsl/gsl_statistics.h>
49#include <gsl/gsl_wavelet.h>
50#include <gsl/gsl_wavelet2d.h>
51#include <gsl/gsl_multimin.h>
53#include <gsl/gsl_fit.h>
57 typedef unsigned int uint;
61 typedef unsigned int uint;
64#ifdef EMAN2_USING_CUDA
72const string SNREvalProcessor::NAME =
"eval.maskedsnr";
73const string AmpweightFourierProcessor::NAME =
"filter.ampweight";
74const string Axis0FourierProcessor::NAME =
"filter.xyaxes0";
75const string ConvolutionProcessor::NAME =
"math.convolution";
76const string BispecSliceProcessor::NAME =
"math.bispectrum.slice";
77const string HarmonicProcessor::NAME =
"math.harmonic";
78const string XGradientProcessor::NAME =
"math.edge.xgradient";
79const string YGradientProcessor::NAME =
"math.edge.ygradient";
80const string ZGradientProcessor::NAME =
"math.edge.zgradient";
81const string ImageDivergenceProcessor::NAME =
"math.divergence";
82const string GradientMagnitudeProcessor::NAME =
"math.gradient.magnitude";
83const string GradientDirectionProcessor::NAME =
"math.gradient.direction";
84const string SDGDProcessor::NAME =
"math.laplacian.sdgd";
85const string LaplacianProcessor::NAME =
"math.laplacian";
86const string LaplacianDirectionProcessor::NAME =
"math.laplacian.direction";
87const string LaplacianMagnitudeProcessor::NAME =
"math.laplacian.magnitude";
88const string Wiener2DAutoAreaProcessor::NAME =
"filter.wiener2dauto";
89const string Wiener2DFourierProcessor::NAME =
"filter.wiener2d";
90const string CtfSimProcessor::NAME =
"math.simulatectf";
91const string LinearRampFourierProcessor::NAME =
"filter.linearfourier";
92const string LoGFourierProcessor::NAME =
"filter.LoG";
93const string DoGFourierProcessor::NAME =
"filter.DoG";
94const string AzSharpProcessor::NAME =
"filter.azimuthal.contrast";
95const string HighpassAutoPeakProcessor::NAME =
"filter.highpass.autopeak";
96const string LinearRampProcessor::NAME =
"eman1.filter.ramp";
97const string AbsoluteValueProcessor::NAME =
"math.absvalue";
98const string CCCSNRProcessor::NAME =
"math.ccc_snr_wiener";
99const string FloorValueProcessor::NAME =
"math.floor";
100const string BooleanProcessor::NAME =
"threshold.notzero";
101const string KmeansSegmentProcessor::NAME =
"segment.kmeans";
102const string GaussSegmentProcessor::NAME =
"segment.gauss";
103const string DistanceSegmentProcessor::NAME =
"segment.distance";
104const string WatershedProcessor::NAME =
"segment.watershed";
105const string SegmentSubunitProcessor::NAME =
"segment.subunit";
106const string RecipCarefullyProcessor::NAME =
"math.reciprocal";
107const string SubtractOptProcessor::NAME =
"math.sub.optimal";
108const string ValuePowProcessor::NAME =
"math.pow";
109const string ValueSquaredProcessor::NAME =
"math.squared";
110const string ValueSqrtProcessor::NAME =
"math.sqrt";
111const string DiscritizeProcessor::NAME =
"threshold.discritize.sigma";
112const string ToZeroProcessor::NAME =
"threshold.belowtozero";
113const string AboveToZeroProcessor::NAME=
"threshold.abovetozero";
114const string OutlierProcessor::NAME=
"threshold.outlier.localmean";
115const string Rotate180Processor::NAME =
"math.rotate.180";
116const string TransformProcessor::NAME =
"xform";
117const string IntTranslateProcessor::NAME =
"xform.translate.int";
118const string ScaleTransformProcessor::NAME =
"xform.scale";
119const string ApplySymProcessor::NAME =
"xform.applysym";
120const string ClampingProcessor::NAME =
"threshold.clampminmax";
121const string RangeZeroProcessor::NAME =
"threshold.rangetozero";
122const string NSigmaClampingProcessor::NAME =
"threshold.clampminmax.nsigma";
123const string ToMinvalProcessor::NAME =
"threshold.belowtominval";
124const string CutToZeroProcessor::NAME =
"threshold.belowtozero_cut";
125const string BinarizeProcessor::NAME =
"threshold.binary";
127const string BinarizeFourierProcessor::NAME =
"threshold.binary.fourier";
128const string CollapseProcessor::NAME =
"threshold.compress";
129const string LinearXformProcessor::NAME =
"math.linear";
130const string ExpProcessor::NAME =
"math.exp";
131const string FiniteProcessor::NAME =
"math.finite";
132const string RangeThresholdProcessor::NAME =
"threshold.binaryrange";
133const string SigmaProcessor::NAME =
"math.sigma";
134const string LogProcessor::NAME =
"math.log";
135const string MaskSharpProcessor::NAME =
"mask.sharp";
136const string MaskSoftProcessor::NAME =
"mask.soft";
137const string MaskEdgeMeanProcessor::NAME =
"mask.ringmean";
138const string MaskNoiseProcessor::NAME =
"mask.noise";
139const string MaskGaussProcessor::NAME =
"mask.gaussian";
140const string MaskAzProcessor::NAME =
"mask.cylinder";
141const string MaskGaussNonuniformProcessor::NAME =
"mask.gaussian.nonuniform";
142const string MaskGaussInvProcessor::NAME =
"math.gausskernelfix";
143const string GridKernelFixProcessor::NAME =
"math.gridkernelfix";
144const string LinearPyramidProcessor::NAME =
"math.linearpyramid";
145const string SetBitsProcessor::NAME =
"math.setbits";
146const string MakeRadiusSquaredProcessor::NAME =
"math.toradiussqr";
147const string MakeRadiusProcessor::NAME =
"math.toradius";
148const string ComplexNormPixel::NAME =
"complex.normpixels";
149const string ZeroConstantProcessor::NAME =
"mask.contract";
150const string BoxMedianProcessor::NAME =
"eman1.filter.median";
151const string BoxSigmaProcessor::NAME =
"math.localsigma";
152const string NonConvexProcessor::NAME =
"math.nonconvex";
153const string BoxMaxProcessor::NAME =
"math.localmax";
154const string LocalMinAbsProcessor::NAME =
"math.localminabs";
155const string MinusPeakProcessor::NAME =
"math.submax";
156const string PeakOnlyProcessor::NAME =
"mask.onlypeaks";
157const string DiffBlockProcessor::NAME =
"eman1.filter.blockrange";
158const string CutoffBlockProcessor::NAME =
"eman1.filter.blockcutoff";
159const string MaxShrinkProcessor::NAME =
"math.maxshrink";
160const string MinShrinkProcessor::NAME =
"math.minshrink";
161const string MeanShrinkProcessor::NAME =
"math.meanshrink";
162const string MedianShrinkProcessor::NAME =
"math.medianshrink";
163const string FFTResampleProcessor::NAME =
"math.fft.resample";
164const string GradientRemoverProcessor::NAME =
"math.lineargradientfix";
165const string GradientPlaneRemoverProcessor::NAME =
"filter.gradientPlaneRemover";
166const string FlattenBackgroundProcessor::NAME =
"filter.flattenbackground";
167const string RampProcessor::NAME =
"filter.ramp";
168const string VerticalStripeProcessor::NAME =
"math.verticalstripefix";
169const string RealToFFTProcessor::NAME =
"math.realtofft";
170const string SigmaZeroEdgeProcessor::NAME =
"mask.zeroedgefill";
171const string WedgeFillProcessor::NAME =
"mask.wedgefill";
172const string FFTPeakProcessor::NAME =
"mask.fft.peak";
173const string FFTConeProcessor::NAME =
"mask.fft.cone";
174const string FFTWedgeProcessor::NAME =
"mask.fft.wedge";
175const string BeamstopProcessor::NAME =
"mask.beamstop";
176const string MeanZeroEdgeProcessor::NAME =
"mask.dampedzeroedgefill";
177const string AverageXProcessor::NAME =
"math.averageovery";
178const string DecayEdgeProcessor::NAME =
"mask.decayedge2d";
179const string ZeroEdgeRowProcessor::NAME =
"mask.zeroedge2d";
180const string ZeroEdgePlaneProcessor::NAME =
"mask.zeroedge3d";
181const string BilateralProcessor::NAME =
"filter.bilateral";
182const string NormalizeUnitProcessor::NAME =
"normalize.unitlen";
183const string NormalizeUnitSumProcessor::NAME =
"normalize.unitsum";
184const string NormalizeStdProcessor::NAME =
"normalize";
185const string NormalizeMaskProcessor::NAME =
"normalize.mask";
186const string NormalizeRampNormVar::NAME =
"normalize.ramp.normvar";
187const string NormalizeByMassProcessor::NAME =
"normalize.bymass";
188const string NormalizeEdgeMeanProcessor::NAME =
"normalize.edgemean";
189const string NormalizeCircleMeanProcessor::NAME =
"normalize.circlemean";
190const string NormalizeLREdgeMeanProcessor::NAME =
"normalize.lredge";
191const string NormalizeHistPeakProcessor::NAME =
"normalize.histpeak";
192const string NormalizeMaxMinProcessor::NAME =
"normalize.maxmin";
193const string NormalizeRowProcessor::NAME =
"normalize.rows";
194const string NormalizeToLeastSquareProcessor::NAME =
"normalize.toimage";
195const string RotationalAverageProcessor::NAME =
"math.rotationalaverage";
196const string RotationalSubstractProcessor::NAME =
"math.rotationalsubtract";
197const string TransposeProcessor::NAME =
"xform.transpose";
198const string FlipProcessor::NAME =
"xform.flip";
199const string AddNoiseProcessor::NAME =
"math.addnoise";
200const string AddSigmaNoiseProcessor::NAME =
"math.addsignoise";
201const string AddRandomNoiseProcessor::NAME =
"math.addspectralnoise";
202const string FourierToCornerProcessor::NAME =
"xform.fourierorigin.tocorner";
203const string FourierToCenterProcessor::NAME =
"xform.fourierorigin.tocenter";
204const string PhaseToCenterProcessor::NAME =
"xform.phaseorigin.tocenter";
205const string PhaseToCornerProcessor::NAME =
"xform.phaseorigin.tocorner";
206const string AutoMask2DProcessor::NAME =
"mask.auto2d";
207const string AutoMaskAsymUnit::NAME =
"mask.asymunit";
208const string AutoMask3DProcessor::NAME =
"mask.auto3d.thresh";
209const string AutoMask3D2Processor::NAME =
"mask.auto3d";
210const string AutoMaskDustProcessor::NAME =
"mask.dust3d";
211const string AddMaskShellProcessor::NAME =
"mask.addshells";
212const string IterMultiMaskProcessor::NAME =
"mask.addshells.multilevel";
213const string IterBinMaskProcessor::NAME =
"mask.addshells.gauss";
214const string PhaseToMassCenterProcessor::NAME =
"xform.phasecenterofmass";
215const string ToMassCenterProcessor::NAME =
"xform.centerofmass";
216const string ToCenterProcessor::NAME =
"xform.center";
217const string ACFCenterProcessor::NAME =
"xform.centeracf";
218const string CTFCorrProcessor::NAME =
"filter.ctfcorr.simple";
219const string SNRProcessor::NAME =
"eman1.filter.snr";
220const string FileFourierProcessor::NAME =
"eman1.filter.byfile";
221const string FSCFourierProcessor::NAME =
"filter.wiener.byfsc";
222const string SymSearchProcessor::NAME =
"misc.symsearch";
223const string MaskPackProcessor::NAME =
"misc.mask.pack";
224const string LocalNormProcessor::NAME =
"normalize.local";
225const string StripeXYProcessor::NAME =
"math.xystripefix";
226const string BadLineXYProcessor::NAME =
"math.xybadlines";
227const string IndexMaskFileProcessor::NAME =
"mask.fromfile";
229const string PaintProcessor::NAME =
"mask.paint";
230const string DirectionalSumProcessor::NAME =
"misc.directional_sum";
233const string MaxPixelOperator::NAME =
"math.max";
234const string MinPixelOperator::NAME =
"math.min";
235const string MatchSFProcessor::NAME =
"filter.matchto";
236const string SetSFProcessor::NAME =
"filter.setstrucfac";
237const string SetIsoPowProcessor::NAME =
"filter.setisotropicpow";
238const string SmartMaskProcessor::NAME =
"mask.smart";
239const string TestImagePureGaussian::NAME =
"testimage.puregaussian";
240const string TestImageFourierGaussianBand::NAME =
"testimage.fourier.gaussianband";
241const string TestImageFourierNoiseGaussian::NAME =
"testimage.noise.fourier.gaussian";
242const string TestImageFourierNoiseProfile::NAME =
"testimage.noise.fourier.profile";
243const string CTFSNRWeightProcessor::NAME =
"ctf.snr.weight";
244const string TestImageLineWave::NAME =
"testimage.linewave";
245const string TestTomoImage::NAME =
"testimage.tomo.objects";
246const string TestImageGradient::NAME =
"testimage.gradient";
247const string TestImageAxes::NAME =
"testimage.axes";
248const string TestImageGaussian::NAME =
"testimage.gaussian";
249const string TestImageScurve::NAME =
"testimage.scurve";
250const string TestImageSphericalWave::NAME =
"testimage.sphericalwave";
251const string TestImageSinewave::NAME =
"testimage.sinewave";
252const string TestImageSinewaveCircular::NAME =
"testimage.sinewave.circular";
253const string TestImageSquarecube::NAME =
"testimage.squarecube";
254const string TestImageEllipse::NAME =
"testimage.ellipsoid";
255const string TestImageHollowEllipse::NAME =
"testimage.ellipsoid.hollow";
256const string TestImageCirclesphere::NAME =
"testimage.circlesphere";
257const string TestImageNoiseUniformRand::NAME =
"testimage.noise.uniform.rand";
258const string TestImageNoiseGauss::NAME =
"testimage.noise.gauss";
259const string TestImageCylinder::NAME =
"testimage.cylinder";
260const string TestImageDisc::NAME =
"testimage.disc";
261const string CCDNormProcessor::NAME =
"filter.ccdnorm";
262const string WaveletProcessor::NAME =
"basis.wavelet";
263const string EnhanceProcessor::NAME =
"filter.enhance";
264const string TomoTiltEdgeMaskProcessor::NAME =
"tomo.tiltedgemask";
265const string TomoTiltAngleWeightProcessor::NAME =
"tomo.tiltangleweight";
266const string FFTProcessor::NAME =
"basis.fft";
267const string RadialProcessor::NAME =
"mask.radialprofile";
268const string HistogramBin::NAME =
"histogram.bin";
269const string ModelEMCylinderProcessor::NAME =
"math.model_em_cylinder";
270const string ApplyPolynomialProfileToHelix::NAME =
"math.poly_radial_profile";
271const string BinarySkeletonizerProcessor::NAME=
"gorgon.binary_skel";
272const string MirrorProcessor::NAME =
"xform.mirror";
273const string ReverseProcessor::NAME =
"xform.reverse";
274const string NewLowpassTopHatProcessor::NAME =
"filter.lowpass.tophat";
275const string NewHighpassTopHatProcessor::NAME =
"filter.highpass.tophat";
276const string NewBandpassTopHatProcessor::NAME =
"filter.bandpass.tophat";
277const string NewHomomorphicTopHatProcessor::NAME =
"filter.homomorphic.tophat";
278const string NewLowpassGaussProcessor::NAME =
"filter.lowpass.gauss";
279const string LowpassAutoBProcessor::NAME=
"filter.lowpass.autob";
280const string NewHighpassGaussProcessor::NAME =
"filter.highpass.gauss";
281const string NewBandpassGaussProcessor::NAME =
"filter.bandpass.gauss";
282const string NewHomomorphicGaussProcessor::NAME =
"filter.homomorphic.gauss";
283const string NewInverseGaussProcessor::NAME =
"filter.gaussinverse";
284const string GaussZFourierProcessor::NAME =
"filter.lowpass.gaussz";
285const string SHIFTProcessor::NAME =
"filter.shift";
286const string InverseKaiserI0Processor::NAME =
"filter.kaiser_io_inverse";
287const string InverseKaiserSinhProcessor::NAME =
"filter.kaisersinhinverse";
288const string NewRadialTableProcessor::NAME =
"filter.radialtable";
289const string LowpassRandomPhaseProcessor::NAME =
"filter.lowpass.randomphase";
290const string NewLowpassButterworthProcessor::NAME =
"filter.lowpass.butterworth";
291const string NewHighpassButterworthProcessor::NAME =
"filter.highpass.butterworth";
292const string NewHomomorphicButterworthProcessor::NAME =
"filter.homomorphic.butterworth";
293const string NewLowpassTanhProcessor::NAME =
"filter.lowpass.tanh";
294const string NewHighpassTanhProcessor::NAME =
"filter.highpass.tanh";
295const string NewHomomorphicTanhProcessor::NAME =
"filter.homomorphic.tanh";
296const string NewBandpassTanhProcessor::NAME =
"filter.bandpass.tanh";
297const string CTF_Processor::NAME =
"filter.CTF_";
298const string ConvolutionKernelProcessor::NAME =
"filter.convolution.kernel";
299const string RotateInFSProcessor::NAME =
"rotateinfs";
300const string CircularAverageBinarizeProcessor::NAME =
"threshold.binary.circularmean";
301const string ObjDensityProcessor::NAME =
"morph.object.density";
302const string ObjLabelProcessor::NAME =
"morph.object.label";
303const string BwThinningProcessor::NAME =
"morph.thin";
304const string BwMajorityProcessor::NAME =
"morph.majority";
305const string PruneSkeletonProcessor::NAME =
"morph.prune";
306const string ManhattanDistanceProcessor::NAME =
"math.distance.manhattan";
307const string BinaryDilationProcessor::NAME =
"morph.dilate.binary";
308const string BinaryErosionProcessor::NAME =
"morph.erode.binary";
309const string BinaryOpeningProcessor::NAME =
"morph.open.binary";
310const string BinaryClosingProcessor::NAME =
"morph.close.binary";
311const string BinaryMorphGradientProcessor::NAME =
"morph.gradient.binary";
312const string BinaryExternalGradientProcessor::NAME =
"morph.ext_grad.binary";
313const string BinaryInternalGradientProcessor::NAME =
"morph.int_grad.binary";
314const string BinaryTopHatProcessor::NAME =
"morph.tophat.binary";
315const string BinaryBlackHatProcessor::NAME =
"morph.blackhat.binary";
316const string GrowSkeletonProcessor::NAME =
"morph.grow";
317const string FixSignProcessor::NAME =
"math.fixmode";
318const string ZThicknessProcessor::NAME =
"misc.zthick";
319const string ReplaceValuefromListProcessor::NAME =
"misc.colorlabel";
320const string PolyMaskProcessor::NAME =
"mask.poly";
321const string AmpMultProcessor::NAME =
"math.multamplitude";
335 force_add<HighpassAutoPeakProcessor>();
336 force_add<LinearRampFourierProcessor>();
337 force_add<LoGFourierProcessor>();
338 force_add<DoGFourierProcessor>();
340 force_add<FixSignProcessor>();
342 force_add<AmpweightFourierProcessor>();
343 force_add<Axis0FourierProcessor>();
344 force_add<Wiener2DFourierProcessor>();
345 force_add<LowpassAutoBProcessor>();
346 force_add<CtfSimProcessor>();
348 force_add<LinearPyramidProcessor>();
349 force_add<LinearRampProcessor>();
350 force_add<AbsoluteValueProcessor>();
351 force_add<FloorValueProcessor>();
352 force_add<BooleanProcessor>();
353 force_add<KmeansSegmentProcessor>();
354 force_add<GaussSegmentProcessor>();
355 force_add<SegmentSubunitProcessor>();
356 force_add<DistanceSegmentProcessor>();
357 force_add<ValuePowProcessor>();
358 force_add<ValueSquaredProcessor>();
359 force_add<ValueSqrtProcessor>();
360 force_add<DiscritizeProcessor>();
361 force_add<Rotate180Processor>();
362 force_add<TransformProcessor>();
363 force_add<ScaleTransformProcessor>();
364 force_add<ApplySymProcessor>();
365 force_add<IntTranslateProcessor>();
366 force_add<RecipCarefullyProcessor>();
367 force_add<SubtractOptProcessor>();
369 force_add<ClampingProcessor>();
370 force_add<NSigmaClampingProcessor>();
372 force_add<ToZeroProcessor>();
373 force_add<RangeZeroProcessor>();
374 force_add<AboveToZeroProcessor>();
375 force_add<OutlierProcessor>();
376 force_add<ToMinvalProcessor>();
377 force_add<CutToZeroProcessor>();
378 force_add<BinarizeProcessor>();
380 force_add<BinarizeFourierProcessor>();
381 force_add<CollapseProcessor>();
382 force_add<LinearXformProcessor>();
383 force_add<SetBitsProcessor>();
385 force_add<CCCSNRProcessor>();
386 force_add<ExpProcessor>();
387 force_add<RangeThresholdProcessor>();
388 force_add<SigmaProcessor>();
389 force_add<LogProcessor>();
390 force_add<FiniteProcessor>();
392 force_add< BinaryOperateProcessor<MaxPixelOperator> >();
393 force_add< BinaryOperateProcessor<MinPixelOperator> >();
395 force_add<PaintProcessor>();
396 force_add<WatershedProcessor>();
397 force_add<MaskSharpProcessor>();
398 force_add<MaskSoftProcessor>();
399 force_add<MaskEdgeMeanProcessor>();
400 force_add<MaskNoiseProcessor>();
401 force_add<MaskGaussProcessor>();
402 force_add<MaskGaussNonuniformProcessor>();
403 force_add<MaskGaussInvProcessor>();
404 force_add<GridKernelFixProcessor>();
405 force_add<MaskAzProcessor>();
407 force_add<MaxShrinkProcessor>();
408 force_add<MinShrinkProcessor>();
409 force_add<MeanShrinkProcessor>();
410 force_add<MedianShrinkProcessor>();
411 force_add<FFTResampleProcessor>();
412 force_add<NonConvexProcessor>();
414 force_add<MakeRadiusSquaredProcessor>();
415 force_add<MakeRadiusProcessor>();
417 force_add<ComplexNormPixel>();
419 force_add<LaplacianProcessor>();
422 force_add<BoxMedianProcessor>();
423 force_add<BoxSigmaProcessor>();
424 force_add<BoxMaxProcessor>();
425 force_add<LocalMinAbsProcessor>();
427 force_add<MinusPeakProcessor>();
428 force_add<PeakOnlyProcessor>();
429 force_add<DiffBlockProcessor>();
431 force_add<CutoffBlockProcessor>();
433 force_add<GradientPlaneRemoverProcessor>();
434 force_add<FlattenBackgroundProcessor>();
435 force_add<VerticalStripeProcessor>();
436 force_add<RealToFFTProcessor>();
437 force_add<SigmaZeroEdgeProcessor>();
438 force_add<WedgeFillProcessor>();
439 force_add<FFTPeakProcessor>();
440 force_add<FFTConeProcessor>();
441 force_add<FFTWedgeProcessor>();
442 force_add<RampProcessor>();
444 force_add<BeamstopProcessor>();
445 force_add<MeanZeroEdgeProcessor>();
446 force_add<AverageXProcessor>();
447 force_add<DecayEdgeProcessor>();
448 force_add<ZeroEdgeRowProcessor>();
449 force_add<ZeroEdgePlaneProcessor>();
451 force_add<BilateralProcessor>();
453 force_add<ConvolutionProcessor>();
454 force_add<BispecSliceProcessor>();
455 force_add<HarmonicProcessor>();
457 force_add<NormalizeStdProcessor>();
458 force_add<NormalizeHistPeakProcessor>();
459 force_add<NormalizeUnitProcessor>();
460 force_add<NormalizeUnitSumProcessor>();
461 force_add<NormalizeMaskProcessor>();
462 force_add<NormalizeEdgeMeanProcessor>();
463 force_add<NormalizeCircleMeanProcessor>();
464 force_add<NormalizeLREdgeMeanProcessor>();
465 force_add<NormalizeMaxMinProcessor>();
466 force_add<NormalizeByMassProcessor>();
467 force_add<NormalizeRowProcessor>();
468 force_add<NormalizeRampNormVar>();
470 force_add<HistogramBin>();
472 force_add<NormalizeToLeastSquareProcessor>();
474 force_add<RotationalAverageProcessor>();
475 force_add<RotationalSubstractProcessor>();
476 force_add<FlipProcessor>();
477 force_add<TransposeProcessor>();
478 force_add<MirrorProcessor>();
479 force_add<ReverseProcessor>();
481 force_add<AddNoiseProcessor>();
482 force_add<AddSigmaNoiseProcessor>();
483 force_add<AddRandomNoiseProcessor>();
485 force_add<PhaseToCenterProcessor>();
486 force_add<PhaseToCornerProcessor>();
487 force_add<FourierToCenterProcessor>();
488 force_add<FourierToCornerProcessor>();
489 force_add<AutoMask2DProcessor>();
490 force_add<AutoMask3DProcessor>();
491 force_add<AutoMask3D2Processor>();
492 force_add<AutoMaskDustProcessor>();
493 force_add<AddMaskShellProcessor>();
494 force_add<IterMultiMaskProcessor>();
495 force_add<IterBinMaskProcessor>();
496 force_add<AutoMaskAsymUnit>();
498 force_add<CTFSNRWeightProcessor>();
500 force_add<ToMassCenterProcessor>();
501 force_add<ToCenterProcessor>();
502 force_add<PhaseToMassCenterProcessor>();
503 force_add<ACFCenterProcessor>();
505 force_add<CTFCorrProcessor>();
506 force_add<FSCFourierProcessor>();
508 force_add<XGradientProcessor>();
509 force_add<YGradientProcessor>();
510 force_add<ZGradientProcessor>();
512 force_add<ImageDivergenceProcessor>();
513 force_add<GradientMagnitudeProcessor>();
514 force_add<GradientDirectionProcessor>();
515 force_add<LaplacianMagnitudeProcessor>();
516 force_add<LaplacianDirectionProcessor>();
517 force_add<SDGDProcessor>();
521 force_add<SymSearchProcessor>();
522 force_add<MaskPackProcessor>();
523 force_add<StripeXYProcessor>();
524 force_add<BadLineXYProcessor>();
525 force_add<LocalNormProcessor>();
527 force_add<IndexMaskFileProcessor>();
529 force_add<SetIsoPowProcessor>();
530 force_add<SetSFProcessor>();
531 force_add<MatchSFProcessor>();
533 force_add<SmartMaskProcessor>();
535 force_add<TestImageGaussian>();
536 force_add<TestImagePureGaussian>();
537 force_add<TestImageSinewave>();
538 force_add<TestImageSphericalWave>();
539 force_add<TestImageSinewaveCircular>();
540 force_add<TestImageSquarecube>();
541 force_add<TestImageCirclesphere>();
542 force_add<TestImageAxes>();
543 force_add<TestImageNoiseUniformRand>();
544 force_add<TestImageNoiseGauss>();
545 force_add<TestImageScurve>();
546 force_add<TestImageCylinder>();
547 force_add<TestImageDisc>();
548 force_add<TestImageGradient>();
549 force_add<TestTomoImage>();
550 force_add<TestImageLineWave>();
551 force_add<TestImageEllipse>();
552 force_add<TestImageHollowEllipse>();
553 force_add<TestImageFourierGaussianBand>();
554 force_add<TestImageFourierNoiseGaussian>();
555 force_add<TestImageFourierNoiseProfile>();
557 force_add<TomoTiltEdgeMaskProcessor>();
558 force_add<TomoTiltAngleWeightProcessor>();
560 force_add<NewLowpassTopHatProcessor>();
561 force_add<NewHighpassTopHatProcessor>();
562 force_add<NewBandpassTopHatProcessor>();
563 force_add<NewHomomorphicTopHatProcessor>();
564 force_add<NewLowpassGaussProcessor>();
565 force_add<NewHighpassGaussProcessor>();
566 force_add<NewBandpassGaussProcessor>();
567 force_add<NewHomomorphicGaussProcessor>();
568 force_add<NewInverseGaussProcessor>();
569 force_add<GaussZFourierProcessor>();
570 force_add<LowpassRandomPhaseProcessor>();
571 force_add<NewLowpassButterworthProcessor>();
572 force_add<NewHighpassButterworthProcessor>();
573 force_add<NewHomomorphicButterworthProcessor>();
574 force_add<NewLowpassTanhProcessor>();
575 force_add<NewHighpassTanhProcessor>();
576 force_add<NewBandpassTanhProcessor>();
577 force_add<NewHomomorphicTanhProcessor>();
578 force_add<NewRadialTableProcessor>();
579 force_add<InverseKaiserI0Processor>();
580 force_add<InverseKaiserSinhProcessor>();
581 force_add<CCDNormProcessor>();
582 force_add<CTF_Processor>();
583 force_add<SHIFTProcessor>();
586 force_add<EnhanceProcessor>();
587 force_add<FFTProcessor>();
588 force_add<RadialProcessor>();
590 force_add<DirectionalSumProcessor>();
591 force_add<ConvolutionKernelProcessor>();
594 force_add<ModelEMCylinderProcessor>();
595 force_add<ApplyPolynomialProfileToHelix>();
596 force_add<BinarySkeletonizerProcessor>();
597 force_add<RotateInFSProcessor>();
598 force_add<CircularAverageBinarizeProcessor>();
599 force_add<ObjDensityProcessor>();
600 force_add<ObjLabelProcessor>();
601 force_add<BwThinningProcessor>();
602 force_add<BwMajorityProcessor>();
603 force_add<PruneSkeletonProcessor>();
604 force_add<GrowSkeletonProcessor>();
606 force_add<ManhattanDistanceProcessor>();
607 force_add<BinaryDilationProcessor>();
608 force_add<BinaryErosionProcessor>();
609 force_add<BinaryClosingProcessor>();
610 force_add<BinaryOpeningProcessor>();
611 force_add<BinaryMorphGradientProcessor>();
612 force_add<BinaryExternalGradientProcessor>();
613 force_add<BinaryInternalGradientProcessor>();
614 force_add<BinaryTopHatProcessor>();
615 force_add<BinaryBlackHatProcessor>();
616 force_add<ZThicknessProcessor>();
617 force_add<ReplaceValuefromListProcessor>();
618 force_add<PolyMaskProcessor>();
619 force_add<AmpMultProcessor>();
629void FiniteProcessor::process_pixel(
float *
x)
const
639 EMData * result = image->copy();
654 if (image->is_complex()) {
655 (*image) *= *processor_image;
658 EMData *fft = image->do_fft();
659 (*fft) *= (*processor_image);
660 EMData *ift = fft->do_ift();
662 float *data = image->get_data();
664 float *ift_data = ift->get_data();
687#define FFTRADIALOVERSAMPLE 4
698 float step=0.5f/array_size;
701 vector < float >yarray(array_size);
705 if (image->is_complex()) {
709 EMData *fft = image->do_fft();
711 EMData *ift = fft->do_ift();
713 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*
sizeof(
float));
729 if (return_radial) image->set_attr(
"filter_curve",yarray);
752 if (image->get_zsize()>1) cornerscale=
sqrt(3.0);
753 else cornerscale=
sqrt(2.0);
755 if (image->is_complex()) {
756 vector <float>yarray = image->
calc_radial_dist(floor(image->get_ysize()*cornerscale/2),0,1.0,1);
759 if (return_radial) image->set_attr(
"filter_curve",yarray);
762 EMData *fft = image->do_fft();
763 vector <float>yarray = fft->
calc_radial_dist((
int)floor(fft->get_ysize()*cornerscale/2.0),0,1.0,1);
768 EMData *ift = fft->do_ift();
770 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*
sizeof(
float));
771 if (return_radial) image->set_attr(
"filter_curve",yarray);
789 if ((
float)mask->get_attr(
"mean_nonzero")!=1.0f)
throw InvalidParameterException(
"MaskPackProcessor requires a binary mask");
791 int n_nz=(int)mask->get_attr(
"square_sum");
795 ret =
new EMData(mask->get_xsize(),mask->get_ysize(),mask->get_zsize());
799 for (
size_t i=0; i<ret->get_size(); i++) {
800 if (mask->get_value_at_index(i)) ret->set_value_at_index(i,image->get_value_at_index(ix++));
804 ret =
new EMData(n_nz,1,1);
808 for (
size_t i=0; i<image->get_size(); i++) {
809 if (mask->get_value_at_index(i)) ret->set_value_at_index(ix++,image->get_value_at_index(i));
825 string mode=(string)
params[
"mode"];
827 int nx = image->get_xsize();
828 int ny = image->get_ysize();
829 int nz = image->get_zsize();
835 if (nx&1) fft->set_fftpad(1);
839 if (mode==
"gridding_5") {
840 for (
int z=-2; z<3; z++) {
841 for (
int y=-2;
y<3;
y++) {
842 for (
int x=0;
x<3;
x++) {
848 else if (mode==
"gridding_7") {
849 for (
int z=-3; z<4; z++) {
850 for (
int y=-3;
y<4;
y++) {
851 for (
int x=0;
x<4;
x++) {
857 else throw InvalidParameterException(
"Gridding kernel correction of unknown mode, only gridding_5 or gridding_7 allowed");
860 real->process_inplace(
"xform.phaseorigin.tocenter");
861 real->mult(2.0f/(
float)
real->get_attr(
"maximum"));
862 real->process_inplace(
"math.reciprocal");
863 real->process_inplace(
"threshold.clampminmax",
Dict(
"minval",-4.0f,
"maxval",4.0f));
883 if (image->is_complex()) fft=image;
884 else EMData *fft = image->do_fft();
886 int nx=fft->get_xsize();
887 int ny=fft->get_ysize();
888 int nz=fft->get_zsize();
894 if (image->is_complex()) {
895 EMData *ift = fft->do_ift();
897 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*
sizeof(
float));
911 if (!image->is_complex()) fft = image->do_fft();
914 int nx=fft->get_xsize();
915 int ny=fft->get_ysize();
916 int nz=fft->get_zsize();
924 for (
int i=0; i<nx/2; i++) thr[i]=thresh_sigma*
sqrt(thr[i]-amp[i]*amp[i])+amp[i];
927 for (
int z=-nz/2; z<nz/2; z++) {
928 for (
int y=-ny/2;
y<ny/2;
y++) {
929 for (
int x=0;
x<nx/2;
x++) {
932 if (r>=nx/2)
continue;
933 complex<float> v=fft->get_complex_at(
x,
y,z);
934 float va=std::abs(v);
939 fft->set_complex_at(
x,
y,z,v);
941 else fft->set_complex_at(
x,
y,z,0.0);
948 fft->set_complex_at(
x,
y,z,v);
950 else fft->set_complex_at(
x,
y,z,0.0);
958 for (
int y=-ny/2;
y<ny/2;
y++) {
959 for (
int x=0;
x<nx/2;
x++) {
961 if (r>=nx/2)
continue;
962 complex<float> v=fft->get_complex_at(
x,
y);
963 float va=std::abs(v);
968 fft->set_complex_at(
x,
y,v);
970 else fft->set_complex_at(
x,
y,0.0);
977 fft->set_complex_at(
x,
y,v);
979 else fft->set_complex_at(
x,
y,0.0);
987 EMData *ift=fft->do_ift();
988 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*
sizeof(
float));
1002 if (!image->is_complex()) fft = image->do_fft();
1006 int nx=fft->get_xsize();
1007 int ny=fft->get_ysize();
1008 int nz=fft->get_zsize();
1014 for (
int z=-nz/2; z<nz/2; z++) {
1015 for (
int y=-ny/2;
y<ny/2;
y++) {
1016 for (
int x=0;
x<nx/2;
x++) {
1018 if (
x!=0 ||
y!=0) ang=90.0-atan(fabs(
float(z)/nz)/hypot(
float(
x)/nx,
float(
y)/ny))*180.0/M_PI;
1020 fft->set_complex_at(
x,
y,z,0.0f);
1026 EMData *ift=fft->do_ift();
1027 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*
sizeof(
float));
1041 if (!image->is_complex()) fft = image->do_fft();
1045 int nx=fft->get_xsize();
1046 int ny=fft->get_ysize();
1047 int nz=fft->get_zsize();
1054 for (
int z=-nz/2; z<nz/2; z++) {
1055 for (
int y=-ny/2;
y<ny/2;
y++) {
1056 for (
int x=0;
x<nx/2;
x++) {
1058 if (z!=0) ang=atan((
float(
y)/ny)/fabs(
float(z)/nz))*180.0/M_PI;
1059 if (ang<anglemin||ang>anglemax ||
Util::hypot3(
x,
y,z)<rmin)
continue;
1060 fft->set_complex_at(
x,
y,z,0.0f);
1066 EMData *ift=fft->do_ift();
1067 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*
sizeof(
float));
1080 if (!image->is_complex())
throw ImageFormatException(
"WedgeFillProcessor: target image must be complex");
1085 if (source && !source->is_complex())
throw ImageFormatException(
"WedgeFillProcessor: fillsource must be complex");
1086 if (image->get_xsize()!=source->get_xsize()||image->get_ysize()!=source->get_ysize()||image->get_zsize()!=source->get_zsize())
throw ImageFormatException(
"WedgeFillProcessor: image/fill size mismatch");
1088 int nx=image->get_xsize();
1089 int ny=image->get_ysize();
1090 int nz=image->get_zsize();
1092 bool dosigma = 1 ? thresh_sigma>0.0 : 0;
1095 bool dotilt = 1 ? maxtilt <90.0 : 0;
1096 maxtilt*=M_PI/180.0;
1098 vector<float> sigmaimg;
1101 for (
int i=0; i<nx/2; i++) sigmaimg[i]*=sigmaimg[i]*thresh_sigma;
1104 vector<int> realpixel(nx/2);
1105 for (
int i=0; i<nx/2; i++) realpixel[i]=0;
1107 for (
int z=0; z<nz; z++) {
1108 for (
int y=0;
y<ny;
y++) {
1109 for (
int x=0;
x<nx;
x+=2) {
1115 if (dotilt) tilt=atan2((
float)(z<nz/2?z:nz-z),(
float)(
x/2));
1117 float v1r=image->get_value_at(
x,
y,z);
1118 float v1i=image->get_value_at(
x+1,
y,z);
1121 if ((!dosigma || v1>sigmaimg[r]) && r<nx/2 && (!dotilt || fabs(tilt)<maxtilt)){
1127 image->set_value_at_fast(
x,
y,z,0);
1128 image->set_value_at_fast(
x+1,
y,z,0);
1131 image->set_value_at_fast(
x,
y,z,source->get_value_at(
x,
y,z));
1132 image->set_value_at_fast(
x+1,
y,z,source->get_value_at(
x+1,
y,z));
1138 image->set_attr(
"real_pixels", realpixel);
1144 float apix=(float)image->get_attr(
"apix_x");
1145 int verbose=(int)
params[
"verbose"];
1149 float ds=1.0f/(apix*image->get_ysize());
1150 unsigned int start=(int)floor(1.0/(15.0*ds));
1151 unsigned int end=radial_mask.size()-2;
1152 if (noisecutoff>0) end=(int)floor(noisecutoff/ds);
1153 if (end>radial_mask.size()-2) {
1154 if (verbose) printf(
"WARNING: specified noisecutoff too close to Nyquist, reset !");
1155 end=radial_mask.size()-2;
1158 printf(
"WARNING: noise cutoff too close to 15 A ! Results will not be good...");
1164 printf(
"Autob -> %d - %d ds=%g apix=%g rdlmsk=%d\n",start,end,ds,apix,
int(radial_mask.size()));
1165 out=fopen(
"fitplot.txt",
"w");
1167 int N=(radial_mask.size()-start-2);
1168 float *
x=(
float *)malloc(N*
sizeof(
float));
1169 float *
y=(
float *)malloc(N*
sizeof(
float));
1170 float *dy=(
float *)malloc(N*
sizeof(
float));
1171 for (
unsigned int i=start; i<radial_mask.size()-2; i++ ) {
1172 x[i-start]=ds*ds*i*i;
1173 if (radial_mask[i]>0)
y[i-start]=
log(radial_mask[i]);
1174 else if (i>start)
y[i-start]=
y[i-start-1];
1175 else y[i-start]=0.0;
1176 if (i>start &&i<radial_mask.size()-3) dy[i-start]=
y[i-start]-
y[i-start-1];
1177 if (out) fprintf(out,
"%f\t%f\n",
x[i-start],
y[i-start]);
1179 if (out) fclose(out);
1181 float slope=0,intercept=0;
1184 if (verbose) printf(
"slope=%f intercept=%f\n",slope,intercept);
1186 float B=(float)
params[
"bfactor"]+slope*4.0f/2.0f;
1187 float B2=(float)
params[
"bfactor"];
1189 if (verbose) printf(
"User B = %1.2f Corrective B = %1.2f Total B=%1.3f\n",(
float)
params[
"bfactor"],slope*2.0,B);
1191 float cutval=exp(-B*pow(end*ds,2.0f)/4.0f)/exp(-B2*pow(end*ds,2.0f)/4.0f);
1192 for (
unsigned int i=0; i<radial_mask.size(); i++) {
1193 if (i<=end) radial_mask[i]=exp(-B*pow(i*ds,2.0f)/4.0f);
1194 else radial_mask[i]=cutval*exp(-B2*pow(i*ds,2.0f)/4.0f);
1204 size_t nxyz = image->get_size();
1208 for (
size_t i=0; i<nxyz; i++) {
1209 float v=image->get_value_at_index(i);
1210 float snr=(v>=.9999)?10000.0f:scale*v/(1.0f-v);
1211 if (snr<0) snr=0.0f;
1212 if (mode) image->set_value_at_index(i,snr/(1.0f+snr));
1213 else image->set_value_at_index(i,snr);
1220 int ys=image->get_ysize();
1223 mask1->process_inplace(
"mask.gaussian",
Dict(
"outer_radius", ys/2.0));
1224 EMData *mask2=mask1->copy();
1227 mask2->process_inplace(
"mask.decayedge2d",
Dict(
"width",4));
1263 if (!image->is_complex()) fft = image->do_fft();
1269 int nx=fft->get_xsize()/2;
1270 int ny=fft->get_ysize();
1273 for (
int x=1;
x<nx;
x++) fft->set_complex_at(
x,0,(fft->get_complex_at(
x,1)+fft->get_complex_at(
x,-1))/neighbornorm);
1276 for (
int x=1;
x<nx;
x++) fft->set_complex_at(
x,0,0.0f);
1281 for (
int y=1;
y<ny;
y++) fft->set_complex_at(0,
y,(fft->get_complex_at(-1,
y)+fft->get_complex_at(1,
y))/neighbornorm);
1284 for (
int y=1;
y<ny;
y++) fft->set_complex_at(0,
y,0.0f);
1290 EMData *ift=fft->do_ift();
1291 memcpy(image->get_data(),ift->get_data(),(nx*2-2)*ny*
sizeof(
float));
1303 image->set_attr(
"apix_x", (
float)
params[
"apix"]);
1304 image->set_attr(
"apix_y", (
float)
params[
"apix"]);
1305 image->set_attr(
"apix_z", (
float)
params[
"apix"]);
1310 const Dict dict = image->get_attr_dict();
1312 float val = (float)
params[
"cutoff_freq"] * (
float)dict[
"apix_x"];
1313 params[
"cutoff_abs"] = val;
1316 float val = ((float)
params[
"cutoff_pixels"] / (
float)dict[
"nx"]);
1317 params[
"cutoff_abs"] = val;
1320 float omega =
params[
"cutoff_abs"];
1324 omega = (omega<0?-1.0:1.0)*0.5f/omega/omega;
1334 if (!image->is_complex()) {
1335 fft = image->do_fft();
1342 int nx=fft->get_xsize();
1343 int ny=fft->get_ysize();
1344 int nz=fft->get_ysize();
1346 zcenter=zcenter*(float)dict[
"apix_x"]*nz;
1349 for (
int x=0;
x<nx/2;
x++) {
1350 for (
int z=(
x==0?0:-nz/2); z<nz/2; z++) {
1351 for (
int y=(
x==0&&z==0?0:-ny/2);
y<ny/2;
y++) {
1352 std::complex <float> v=fft->get_complex_at(
x,
y,z);
1354 fft->set_complex_at(
x,
y,z,v*exp(-omega*(r-zcenter)*(r-zcenter))*(r>hppix?1.0f:
Util::hypot3(
x,
y,z)/hppix));
1360 for (
int z=-nz/2; z<nz/2; z++) {
1361 for (
int y=-ny/2;
y<ny/2;
y++) {
1362 for (
int x=0;
x<nx/2;
x++) {
1363 std::complex <float> v=fft->get_complex_at(
x,
y,z);
1364 fft->set_complex_at(
x,
y,z,v*exp(-omega*(abs(z)-zcenter)*(abs(z)-zcenter))*(fabs(z)>hppix?1.0f:
Util::hypot3(
x,
y,z)/hppix));
1371 EMData *ift=fft->do_ift();
1372 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*
sizeof(
float));
1393 if (!image->is_complex()) {
1394 fft = image->do_fft();
1395 fftd = fft->get_data();
1400 fftd=image->get_data();
1403 if (
sum) sumd=
sum->get_data();
1405 size_t n = (size_t)fft->get_xsize()*fft->get_ysize()*fft->get_zsize();
1406 for (
size_t i=0; i<n; i+=2) {
1408 if (
dosqrt) c=pow(fftd[i]*fftd[i]+fftd[i+1]*fftd[i+1],0.25f);
1409 else c =
static_cast<float>(hypot(fftd[i],fftd[i+1]));
1410 if (c==0) c=1.0e-30f;
1413 if (sumd) { sumd[i]+=c; sumd[i+1]+=0; }
1430 EMData *ift=fft->do_ift();
1431 memcpy(image->get_data(),ift->get_data(),n*
sizeof(
float));
1443 printf(
"Process inplace not implemented. Please use process.\n");
1457 float apix=image->get_attr(
"apix_x");
1458 int nx=image->get_xsize();
1459 int ny=image->get_ysize();
1460 int nz=image->get_zsize();
1462 EMData * result = image->process(
"threshold.belowtozero",
Dict(
"minval",0.0f));
1467 result->mult(*mask);
1468 result->process_inplace(
"normalize.local",
Dict(
"radius",nx*apix/(3.0f*width),
"threshold",(
float)result->get_attr(
"sigma_nonzero")*1.2));
1474 result->process_inplace(
"filter.setstrucfac",
Dict(
"apix",apix,
"strucfac",&gsf));
1477 int gs=2*width/(float)image->get_attr(
"apix_x");
1480 gsub.process_inplace(
"mask.gaussian",
Dict(
"outer_radius",
int(width/(2.0*apix))));
1483 result->set_attr(
"render_bits",12);
1484 result->set_attr(
"render_min",(
float)result->get_attr(
"minimum"));
1485 result->set_attr(
"render_max",(
float)result->get_attr(
"maximum"));
1490 if (skipseg==2) cache=result->copy();
1515 vector<float> centers;
1517 if (maxnseg==0) maxnseg=2000000000;
1519 float thr=(float)result->get_attr(
"maximum")*minratio;
1520 while (amps.size()<maxnseg) {
1521 if ((
float)result->get_attr(
"maximum")<=thr)
break;
1523 vector<Pixel> pixels=result->calc_highest_locations(thr);
1524 if (pixels.size()==0)
break;
1525 if (verbose>1) printf(
"%d %f %f %f %d\n",pixels.size(),pixels[0].get_value(),pixels[1].get_value(),pixels[2].get_value(),amps.size());
1527 for (
int i=0; i<pixels.size(); i++) {
1528 IntPoint p = pixels[i].get_point();
1531 float amp=result->get_value_at(p[0],p[1],p[2]);
1532 if (amp!=pixels[i].get_value())
continue;
1535 if (amp<thr)
continue;
1537 amps.push_back(amp);
1538 centers.push_back(p[0]);
1539 centers.push_back(p[1]);
1540 centers.push_back(p[2]);
1542 result->insert_scaled_sum(&gsub,fp,1.0,-amp);
1546 if (verbose) printf(
"Found %d centers\n",amps.size());
1549 result->set_attr(
"render_bits",12);
1550 result->set_attr(
"render_min",(
float)result->get_attr(
"minimum"));
1551 result->set_attr(
"render_max",(
float)result->get_attr(
"maximum"));
1556 for (
int z=0; z<nz; z++) {
1557 for (
int y=0;
y<ny;
y++) {
1558 for (
int x=0;
x<nz;
x++) {
1564 float bdist=(float)(nx+ny+nz);
1565 for (
unsigned int c=0; c<centers.size()/3; c++) {
1566 float d=
Util::hypot3(
x-centers[c*3],
y-centers[c*3+1],z-centers[c*3+2]);
1567 if (d<bdist) { bdist=d; bcls=c; }
1569 result->set_value_at(
x,
y,z,(
float)bcls);
1573 if (verbose) printf(
"segmented\n");
1582 result->set_attr(
"segment_centers",centers);
1583 result->set_attr(
"segment_amps",amps);
1590 printf(
"Process inplace not implemented. Please use process.\n");
1597 EMData * result = image->copy();
1604 vector<Pixel> pixels=image->calc_highest_locations(thr);
1606 vector<float> centers(3);
1607 int nx=image->get_xsize();
1608 int ny=image->get_ysize();
1609 int nz=image->get_zsize();
1613 centers[0]=(float)pixels[0].
x;
1614 centers[1]=(float)pixels[0].
y;
1615 centers[2]=(float)pixels[0].z;
1616 pixels.erase(pixels.begin());
1620 while (pixels.size()>0) {
1624 for (
unsigned int i=0; i<pixels.size(); i++) {
1628 for (
unsigned int j=0; j<centers.size(); j+=3) {
1629 float d=
Util::hypot3(centers[j]-p.
x,centers[j+1]-p.
y,centers[j+2]-p.
z);
1631 pixels.erase(pixels.begin()+i);
1639 for (
unsigned int i=0; i<pixels.size() && found==0; i++) {
1643 for (
unsigned int j=centers.size()-3; j>0; j-=3) {
1644 float d=
Util::hypot3(centers[j]-p.
x,centers[j+1]-p.
y,centers[j+2]-p.
z);
1646 centers.push_back((
float)p.
x);
1647 centers.push_back((
float)p.
y);
1648 centers.push_back((
float)p.
z);
1649 pixels.erase(pixels.begin()+i);
1657 if (!found && pixels.size()) {
1658 if (verbose) printf(
"New chain\n");
1659 centers.push_back((
float)pixels[0].
x);
1660 centers.push_back((
float)pixels[0].
y);
1661 centers.push_back((
float)pixels[0].z);
1662 pixels.erase(pixels.begin());
1665 if (verbose) printf(
"%d points found\n",(
int)(centers.size()/3));
1669 for (
int z=0; z<nz; z++) {
1670 for (
int y=0;
y<ny;
y++) {
1671 for (
int x=0;
x<nz;
x++) {
1672 if (image->get_value_at(
x,
y,z)<thr) {
1673 result->set_value_at(
x,
y,z,-1.0);
1677 float bdist=(float)(nx+ny+nz);
1678 for (
unsigned int c=0; c<centers.size()/3; c++) {
1679 float d=
Util::hypot3(
x-centers[c*3],
y-centers[c*3+1],z-centers[c*3+2]);
1680 if (d<bdist) { bdist=d; bcls=c; }
1682 result->set_value_at(
x,
y,z,(
float)bcls);
1687 result->set_attr(
"segment_centers",centers);
1695 if (s.length()<2)
return image->copy();
1696 int n=atoi(s.c_str()+1);
1697 if ((s[0]==
'c' || s[0]==
'C') && n==1)
return image->copy();
1701 if (image->get_zsize()==1) {
1702 if (s[0]!=
'c' && s[0]!=
'C')
throw ImageDimensionException(
"xform.applysym: Cn symmetry required for 2-D symmetrization");
1705 for (
int i=0; i<n; i++) {
1707 EMData* transformed = image->process(
"xform",
Dict(
"transform",&t));
1717 vector<Transform> transforms = sym->
get_syms();
1719 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
1721 EMData* transformed = image->process(
"xform",
Dict(
"transform",&t));
1733 memcpy(image->get_data(),tmp->get_data(),(
size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*
sizeof(
float));
1741 EMData * result = image->copy();
1754 int nx=image->get_xsize();
1755 int ny=image->get_ysize();
1756 int nz=image->get_zsize();
1760 thr=float(image->get_attr(
"mean_nonzero"))+ 1.0 *float(image->get_attr(
"sigma_nonzero"));
1761 printf(
"Estimated map threshold: %4f\n", thr);
1764 vector<float> centers(nseg*3);
1765 vector<float> count(nseg);
1768 float ax=image->get_attr(
"apix_x");
1770 if (verbose) printf(
"Seeding .....\n");
1771 int sx=int(nx/sep)+1,sy=int(ny/sep)+1,sz=int(nz/sep)+1;
1774 for (
int i=0; i<nx; i++){
1775 for (
int j=0; j<ny; j++){
1776 for (
int k=0; k<nz; k++){
1777 int ni=(i/sep),nj=(j/sep),nk=(k/sep);
1778 float v=image->get_value_at(i,j,k);
1780 m.set_value_at(ni,nj,nk,(m.get_value_at(ni,nj,nk)+v));
1781 mcount.set_value_at(ni,nj,nk,(mcount.get_value_at(ni,nj,nk)+1));
1786 m.div((nx/sx)*(ny/sy)*(nz/sz));
1788 float l=image->get_attr(
"minimum"),r=image->get_attr(
"maximum"),th=0;
1789 while (abs(nsum-nseg)>0){
1792 for (
int i=0; i<sx; i++){
1793 for (
int j=0; j<sy; j++){
1794 for (
int k=0; k<sz; k++){
1795 if (m.get_value_at(i,j,k)>th) nsum+=1;
1800 if (nsum>nseg) l=th;
1801 if (nsum<nseg) r=th;
1802 if ((r-l)<.0001)
break;
1805 if (verbose) printf(
"%3d pseudoatoms seeded at threshold value %3f\n", nsum, th);
1807 for (
int i=0; i<sx; i++){
1808 for (
int j=0; j<sy; j++){
1809 for (
int k=0; k<sz; k++){
1810 if (m.get_value_at(i,j,k)>th){
1812 centers[q]= float(i+.5)*sep;
1813 centers[q+1]=float(j+.5)*sep;
1814 centers[q+2]=float(k+.5)*sep;
1825 for (
int i=0; i<nseg*3; i+=3) {
1832 for (
int iter=0; iter<maxiter; iter++) {
1835 for (
int z=0; z<nz; z++) {
1836 for (
int y=0;
y<ny;
y++) {
1837 for (
int x=0;
x<nx;
x++) {
1838 if (image->get_value_at(
x,
y,z)<thr) {
1839 result->set_value_at(
x,
y,z,-1.0);
1843 float bdist=(float)(nx+ny+nz);
1844 for (
int c=0; c<nseg; c++) {
1845 float d=
Util::hypot3(
x-centers[c*3],
y-centers[c*3+1],z-centers[c*3+2]);
1846 if (d<bdist) { bdist=d; bcls=c; }
1848 if ((
int)result->get_value_at(
x,
y,z)!=bcls) pixmov++;
1849 if (bdist>maxsegsize) result->set_value_at(
x,
y,z,-1);
1850 else result->set_value_at(
x,
y,z,(
float)bcls);
1856 for (
int i=0; i<nseg*3; i++) centers[i]=0;
1857 for (
int i=0; i<nseg; i++) count[i]=0;
1859 for (
int z=0; z<nz; z++) {
1860 for (
int y=0;
y<ny;
y++) {
1861 for (
int x=0;
x<nx;
x++) {
1862 int cls = (int)result->get_value_at(
x,
y,z);
1863 if (cls==-1)
continue;
1865 if (ampweight) w=image->get_value_at(
x,
y,z);
1867 centers[cls*3]+=
x*w;
1868 centers[cls*3+1]+=
y*w;
1869 centers[cls*3+2]+=z*w;
1876 for (
int c=0; c<nseg; c++) {
1884 }
while (image->get_value_at((
int)centers[c*3],(
int)centers[c*3+1],(
int)centers[c*3+2])<thr);
1888 centers[c*3]/=count[c];
1889 centers[c*3+1]/=count[c];
1890 centers[c*3+2]/=count[c];
1896 for (
int c1=0; c1<nseg-1; c1++) {
1897 for (
int c2=c1+1; c2<nseg; c2++) {
1898 if (
Util::hypot3(centers[c1*3]-centers[c2*3],centers[c1*3+1]-centers[c2*3+1],centers[c1*3+2]-centers[c2*3+2])<=minsegsep) {
1904 }
while (image->get_value_at((
int)centers[c1*3],(
int)centers[c1*3+1],(
int)centers[c1*3+2])<thr);
1911 if (verbose) printf(
"Iteration %3d: %6ld voxels moved, %3d classes reseeded\n",iter,pixmov,nreseed);
1912 if (nreseed==0 && pixmov<(
size_t)maxvoxmove)
break;
1915 result->set_attr(
"segment_centers",centers);
1922 printf(
"Process inplace not implemented. Please use process.\n");
1931 int nx=image->get_xsize();
1932 int ny=image->get_ysize();
1933 int nz=image->get_zsize();
1941 for (
int z=0; z<nz; z++) {
1942 for (
int y=0;
y<ny;
y++) {
1943 for (
int x=0;
x<nx;
x++) {
1944 float xf=1.0-fabs(
x-x0)*2.0/xwidth;
1945 float yf=1.0-fabs(
y-y0)*2.0/ywidth;
1946 float zf=1.0-fabs(z-z0)*2.0/zwidth;
1948 if (xf>0&&yf>0&&zf>0) val=xf*yf*zf;
1949 image->mult_value_at_fast(
x,
y,z,val);
1984 memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*
sizeof(
float));
1994 if (in->is_complex()) in2=in;
1995 else in=in->do_fft();
1997 EMData *filt = in->copy_head();
2000 if (!ictf)
ctf=(
Ctf *)in->get_attr(
"ctf");
2004 EMData *ret=filt->do_ift();
2007 if (!in->is_complex())
delete in2;
2045 memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*
sizeof(
float));
2053 Assert(radial_mask.size() > 0);
2054 for (
size_t i = 0; i < radial_mask.size(); i++) {
2055 radial_mask[i] = (float)i;
2066 cutoff=(float)
params[
"cutoff_abs"];
2069 printf(
"A cutoff_* parameter is required by filter.lowpass.randomphase\n");
2074 if (image->get_zsize()==1) {
2076 if (!image->is_complex()) { image->do_fft_inplace(); flag=1; }
2078 int nx=image->get_xsize();
2079 int ny=image->get_ysize();
2082 float *data=image->get_data();
2083 for (
int y=-ny/2;
y<ny/2;
y++) {
2084 for (
int x=0;
x<nx/2+1;
x++) {
2085 if (hypot(
x/
float(nx),
y/
float(ny))>=cutoff) {
2086 size_t idx=image->get_complex_index_fast(
x,
y,z);
2095 image->do_ift_inplace();
2101 if (!image->is_complex()) { image->do_fft_inplace(); flag=1; }
2103 int nx=image->get_xsize();
2104 int ny=image->get_ysize();
2105 int nz=image->get_zsize();
2107 float *data=image->get_data();
2108 for (
int z=-nz/2; z<nz/2; z++) {
2109 for (
int y=-ny/2;
y<ny/2;
y++) {
2110 for (
int x=0;
x<nx/2;
x++) {
2112 size_t idx=image->get_complex_index_fast(
x,
y,z);
2121 image->do_ift_inplace();
2130 image->set_attr(
"apix_x", (
float)
params[
"apix"]);
2131 image->set_attr(
"apix_y", (
float)
params[
"apix"]);
2132 image->set_attr(
"apix_z", (
float)
params[
"apix"]);
2135 const Dict dict = image->get_attr_dict();
2141 highpass = (float)
params[
"cutoff_freq"] * (
float)dict[
"apix_x"] * (float)dict[
"ny"] / 2.0f;
2153 for (c=2; c<radial_mask.size(); c++)
if (radial_mask[c-1]<=radial_mask[c])
break;
2158 for (
unsigned int i=1; i<radial_mask.size(); i++) radial_mask[i]=(i<=c?0.0f:1.0f);
2167 Assert(radial_mask.size() > 0);
2168 float x = 0.0f , step = 0.5f/radial_mask.size();
2169 size_t size=radial_mask.size();
2170 for (
size_t i = 0; i < size; i++) {
2179 Assert(radial_mask.size() > 0);
2180 float x = 0.0f , nqstep = 0.5f/radial_mask.size();
2181 size_t size=radial_mask.size();
2183 for (
size_t i = 0; i < size; i++) {
2184 radial_mask[i] = ((
x*
x - var)/var*var)*exp(-
x*
x/2*var);
2192 Assert(radial_mask.size() > 0);
2193 float x = 0.0f , nqstep = 0.5f/radial_mask.size();
2194 size_t size=radial_mask.size();
2195 float norm = 1.0f/
sqrt(2*M_PI);
2196 for (
size_t i = 0; i < size; i++) {
2211 if (gauss_width<=0) {
2212 size_t size = (size_t)image->get_xsize() *
2213 (size_t)image->get_ysize() *
2214 (size_t)image->get_zsize();
2215 float *data = image->get_data();
2217 for (
size_t i = 0; i < size; ++i) {
2223 int nx = image->get_xsize();
2224 int ny = image->get_ysize();
2225 int nz = image->get_zsize();
2226 float wid=gauss_width/(ny*ny);
2228 for (
int z=0; z<nz; z++) {
2229 for (
int y=0;
y<ny;
y++) {
2230 for (
int x=0;
x<nx;
x++) {
2232 if (image->get_value_at(
x,
y,z)>=
minval*cor && image->get_value_at(
x,
y,z)<=
maxval*cor) image->set_value_at(
x,
y,z,0.0f);
2250 maxval = image->get_attr(
"maximum");
2251 mean = image->get_attr(
"mean");
2252 sigma = image->get_attr(
"sigma");
2256 size_t size = (size_t)image->get_xsize() *
2257 (size_t)image->get_ysize() *
2258 (size_t)image->get_zsize();
2259 float *data = image->get_data();
2261 for (
size_t i = 0; i < size; ++i) {
2275 maxval = image->get_attr(
"maximum");
2276 mean = image->get_attr(
"mean");
2277 sigma = image->get_attr(
"sigma");
2278 nx = image->get_xsize();
2279 ny = image->get_ysize();
2280 nz = image->get_zsize();
2290 float *data = image->get_data();
2293 for (
int z = 0; z <
nz; z++) {
2294 for (
int y = 0;
y <
ny;
y++) {
2295 for (
int x = 0;
x <
nx;
x++) {
2305 int nx=image->get_xsize();
2306 int ny=image->get_ysize();
2307 int nz=image->get_zsize();
2314 if (r>
r2 && r>
r1)
continue;
2316 else image->set_value_at(i,j,0,
v1);
2326 if (r>
r2 && r>
r1)
continue;
2328 else image->set_value_at(i,j,k,
v1);
2337 if (image->is_complex())
throw ImageFormatException(
"MaskAzProcessor: target image must be real");
2339 int nx=image->get_xsize();
2340 int ny=image->get_ysize();
2341 int nz=image->get_zsize();
2356 for (
int x=0;
x<nx;
x++) {
2357 for (
int y=0;
y<ny;
y++) {
2358 float az=atan2(
y-cy,
x-cx);
2359 float az2=az-M_PI*2.0f;
2360 float az3=az+M_PI*2.0f;
2361 float r=hypot(
y-cy,
x-cx);
2363 if (r>inner_radius&&r<=outer_radius) {
2365 if (as<phirange) val=1.0f;
2366 else if (!phitriangle || as>phirange+phitrirange) val=0.0f;
2367 else val=1.0-(as-phirange)/
float(phitrirange);
2369 if (r==0 && inner_radius<=0) val=1.0;
2371 for (
int z=0; z<nz; z++) {
2372 if (z<zmin-ztri || z>zmax+ztri) image->mult_value_at_fast(
x,
y,z,0);
2373 else if (z>=zmin+ztri && z<=zmax-ztri) image->mult_value_at_fast(
x,
y,z,val);
2374 else if (z>=zmin-ztri && z<=zmin+ztri) image->mult_value_at_fast(
x,
y,z,val*((z-zmin)/(2.0f*ztri)+0.5));
2375 else image->mult_value_at_fast(
x,
y,z,val*((zmax-z)/(2.0f*ztri)+0.5));
2406 float *data = image->get_data();
2413 for (
int z = 0; z <
nz; ++z) {
2414 for (
int y = 0;
y <
ny; ++
y) {
2415 for (
int x = 0;
x <
nx; ++
x) {
2439 size_t size = (size_t)image->get_xsize() *
2440 (size_t)image->get_ysize() *
2441 (size_t)image->get_zsize();
2442 float *data = image->get_data();
2446 for (
size_t i = 0; i < size; ++i) {
2447 if (data[i]<minval) data[i]=newval;
2459 if (!image->is_complex()) {
2460 LOGWARN(
"cannot apply complex processor on a real image. Nothing is done.");
2464 size_t size = (size_t)image->get_xsize() *
2465 (size_t)image->get_ysize() *
2466 (size_t)image->get_zsize();
2467 float *data = image->get_data();
2471 for (
size_t i = 0; i < size; i += 2) {
2489 float *data = image->get_data();
2491 nx = image->get_xsize();
2492 ny = image->get_ysize();
2493 nz = image->get_zsize();
2505 size_t cpysize =
areasize *
sizeof(float);
2506 size_t start = (
nx *
ny +
nx + 1) * n;
2523 size_t nsec = (size_t)
nx * (
size_t)
ny;
2528 size_t total_size = (size_t)
nx * (
size_t)
ny * (size_t)
nz;
2529 float *data2 =
new float[total_size];
2530 memcpy(data2, data, total_size *
sizeof(
float));
2533 for (
int z = zstart; z < zend; z++) {
2534 for (
int y = n;
y < yend;
y++) {
2535 for (
int x = n;
x < xend;
x++) {
2537 k = (size_t)z * nsec +
y *
nx +
x;
2539 for (
int bz = zbox_start; bz < zbox_end; bz++) {
2540 for (
int by = 0; by <
areasize; by++) {
2541 memcpy(&matrix[(
size_t)bz * box_nsec + by *
areasize],
2542 &data2[k - start + bz * nsec + by *
nx], cpysize);
2568 d = image->process(
"math.edge.xgradient");
2569 d->process_inplace(
"math.edge.xgradient");
2570 image->process_inplace(
"math.edge.ygradient");
2571 image->process_inplace(
"math.edge.xgradient");
2573 image->process_inplace(
"normalize");
2589 kernel[4] = -1.0f / 6.0f;
2590 kernel[10] = -1.0f / 6.0f;
2591 kernel[12] = -1.0f / 6.0f;
2592 kernel[14] = -1.0f / 6.0f;
2593 kernel[16] = -1.0f / 6.0f;
2594 kernel[22] = -1.0f / 6.0f;
2606 float sigma = image->get_attr(
"sigma");
2607 float mean = image->get_attr(
"mean");
2608 float bitmax=pow(2.,bits);
2610 int nx = image->get_xsize();
2611 int ny = image->get_ysize();
2612 int nz = image->get_zsize();
2613 size_t total_size = (size_t)nx * (
size_t)ny * (size_t)nz;
2616 float min=image->get_attr(
"minimum");
2617 float max=image->get_attr(
"maximum");
2622 if (mean-nsigma*sigma>min) min=mean-nsigma*sigma;
2623 if (mean+nsigma*sigma<max) max=mean+nsigma*sigma;
2626 for (
size_t i=0; i<total_size; i++) {
2627 float newval=floor((image->get_value_at_index(i)-min)*bitmax/(max-min)
2630 if (newval<0) newval=0;
2631 if (newval>=bitmax) newval=bitmax-1.0f;
2632 image->set_value_at_index(i,newval);
2638 unsigned int bitmask= 0xffffffff << (23-floatbits);
2639 for (
size_t i=0; i<total_size; i++) {
2640 float newval=image->get_value_at_index(i);
2641 unsigned int *newvali = (
unsigned int*)&newval;
2643 image->set_value_at_index(i,newval);
2651 memcpy(image->get_data(),cpy->get_data(),image->get_size()*
sizeof(
float));
2656inline int MIRE(
int x,
int nx) {
return x<0?-
x:(
x>=nx?nx-(
x-nx+1):
x); }
2666 int nx = image->get_xsize();
2667 int ny = image->get_ysize();
2668 int nz = image->get_zsize();
2677 int matrix_size = (2*dx+1)*(2*dy+1)*(2*dz+1);
2679 vector<float> array(matrix_size);
2682 EMData *ret = image->copy_head();
2687 for (
int k=0; k<nz; k++) {
2688 for (
int j=0; j<ny; j++) {
2689 for (
int i=0; i<nx; i++) {
2691 for (
int kk=k-dz,s=0; kk<=k+dz; kk++) {
2692 for (
int jj=j-dy; jj<=j+dy; jj++) {
2693 for (
int ii=i-dx; ii<=i+dx; ii++,s++) {
2694 array[s]=image->get_value_at(
MIRE(ii,nx),
MIRE(jj,ny),
MIRE(kk,nz));
2698 float newv=image->get_value_at(i,j,k);
2700 ret->set_value_at(i,j,k,newv);
2717 int nz = image->get_zsize();
2724 int nx = image->get_xsize();
2725 int ny = image->get_ysize();
2727 int v1 =
params[
"cal_half_width"];
2728 int v2 =
params[
"fill_half_width"];
2730 int v0 = v1 > v2 ? v1 : v2;
2736 float *data = image->get_data();
2738 for (
int y =
v0;
y <= ny -
v0 - 1;
y += v2) {
2739 for (
int x =
v0;
x <= nx -
v0 - 1;
x += v2) {
2742 for (
int y1 =
y - v1; y1 <=
y + v1; y1++) {
2743 for (
int x1 =
x - v1; x1 <=
x + v1; x1++) {
2744 sum += data[x1 + y1 * nx];
2747 float mean = sum / ((v1 * 2 + 1) * (v1 * 2 + 1));
2749 for (
int j =
y - v2; j <=
y + v2; j++) {
2750 for (
int i =
x - v2; i <=
x + v2; i++) {
2751 data[i + j * nx] = mean;
2767 int nz = image->get_zsize();
2774 int nx = image->get_xsize();
2775 int ny = image->get_ysize();
2777 float value1 =
params[
"value1"];
2778 float value2 =
params[
"value2"];
2780 int v1 = (int) value1;
2781 int v2 = (int) value2;
2783 LOGERR(
"invalid value2 '%f' in CutoffBlockProcessor", value2);
2791 float *data = image->get_data();
2793 for (
y = 0;
y <= ny - v1;
y += v1) {
2794 for (
x = 0;
x <= nx - v1;
x += v1) {
2797 EMData *fft = clip->do_fft();
2799 float *fft_data = fft->get_data();
2803 for (
int i = -v2; i < v2; i++) {
2804 for (
int j = 0; j < v2; j++) {
2805 if (j == 0 && i == 0) {
2809 if (hypot(j, i) < value2) {
2810 int t = j * 2 + (i + v1 / 2) * (v1 + 2);
2811 sum += (fft_data[t] * fft_data[t] + fft_data[t + 1] * fft_data[t + 1]);
2823 float mean = sum / nitems;
2825 for (
int i =
y; i <
y + v1; i++) {
2826 for (
int j =
x; j <
x + v1; j++) {
2827 data[i * nx + j] = mean;
2833 memset(&data[
y * nx], 0, (ny -
y) * nx *
sizeof(
float));
2835 for (
int i = 0; i < ny; i++) {
2836 memset(&data[i * nx +
x], 0, (nx -
x) *
sizeof(
float));
2844 if (image->is_complex())
throw ImageFormatException(
"Error, the median shrink processor does not work on complex images");
2847 if (shrink_factor <= 1) {
2849 "median shrink: shrink factor must > 1");
2852 int nx = image->get_xsize();
2853 int ny = image->get_ysize();
2854 int nz = image->get_zsize();
2861 int shrunken_nx = nx / shrink_factor;
2862 int shrunken_ny = ny / shrink_factor;
2863 int shrunken_nz = 1;
2864 if (nz > 1) shrunken_nz = nz / shrink_factor;
2867 image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
2880 if (image->is_complex())
throw ImageFormatException(
"Error, the median shrink processor does not work on complex images");
2883 if (shrink_factor <= 1) {
2885 "median shrink: shrink factor must > 1");
2887 int nx = image->get_xsize();
2888 int ny = image->get_ysize();
2889 int nz = image->get_zsize();
2897 int shrunken_nx = nx / shrink_factor;
2898 int shrunken_ny = ny / shrink_factor;
2899 int shrunken_nz = 1;
2900 if (nz > 1) shrunken_nz = nz / shrink_factor;
2903 EMData *ret = image->copy_head();
2904 ret->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
2914 int nx_old = from->get_xsize();
2915 int ny_old = from->get_ysize();
2917 int threed_shrink_factor = shrink_factor * shrink_factor;
2918 int z_shrink_factor = 1;
2919 if (from->get_zsize() > 1) {
2920 threed_shrink_factor *= shrink_factor;
2921 z_shrink_factor = shrink_factor;
2924 float *mbuf =
new float[threed_shrink_factor];
2927 int nxy_old = nx_old * ny_old;
2929 int nx = to->get_xsize();
2930 int ny = to->get_ysize();
2931 int nz = to->get_zsize();
2932 int nxy_new = nx * ny;
2934 float *
rdata = to->get_data();
2935 const float *
const data_copy = from->get_const_data();
2937 for (
int l = 0; l < nz; l++) {
2938 int l_min = l * shrink_factor;
2939 int l_max = l * shrink_factor + z_shrink_factor;
2940 size_t cur_l = (size_t)l * nxy_new;
2942 for (
int j = 0; j < ny; j++) {
2943 int j_min = j * shrink_factor;
2944 int j_max = (j + 1) * shrink_factor;
2945 size_t cur_j = j * nx + cur_l;
2947 for (
int i = 0; i < nx; i++) {
2948 int i_min = i * shrink_factor;
2949 int i_max = (i + 1) * shrink_factor;
2952 for (
int l2 = l_min; l2 < l_max; l2++) {
2953 size_t cur_l2 = l2 * nxy_old;
2955 for (
int j2 = j_min; j2 < j_max; j2++) {
2956 size_t cur_j2 = j2 * nx_old + cur_l2;
2958 for (
int i2 = i_min; i2 < i_max; i2++) {
2959 mbuf[k] = data_copy[i2 + cur_j2];
2965 for (k = 0; k < size_t(threed_shrink_factor / 2 + 1); k++) {
2966 for (
int i2 = k + 1; i2 < threed_shrink_factor; i2++) {
2967 if (mbuf[i2] < mbuf[k]) {
2975 rdata[i + cur_j] = mbuf[threed_shrink_factor / 2];
2986 to->scale_pixel((
float)shrink_factor);
2991 if (image->is_complex())
throw ImageFormatException(
"FFTResampleProcessor does not support complex images.");
2994 if (sample_rate <= 0.0F ) {
2998 int nx=image->get_xsize();
2999 int ny=image->get_ysize();
3000 int nz=image->get_zsize();
3002 int nnx=(int)floor(nx/sample_rate+0.5f);
3003 int nny=(ny==1?1:(int)floor(ny/sample_rate+0.5f));
3004 int nnz=(nz==1?1:(int)floor(nz/sample_rate+0.5f));
3006 if (nnx==nx && nny==ny && nnz==nz)
return image->copy();
3010 if (nnx<nx||nny<ny||nnz<nz) {
3012 result=((
EMData *)image)->FourTruncate(nnx, nny, nnz, 1, 0);
3016 result=((
EMData *)image)->FourInterpol(nnx, nny, nnz, 1, 0);
3018 result->scale_pixel((
float)nx/(float)nnx);
3036 if (sample_rate <= 0.0F ) {
3041 int nx=image->get_xsize();
3042 int ny=image->get_ysize();
3043 int nz=image->get_zsize();
3045 int nnx=(int)floor(nx/sample_rate+0.5f);
3046 int nny=(ny==1?1:(int)floor(ny/sample_rate+0.5f));
3047 int nnz=(nz==1?1:(int)floor(nz/sample_rate+0.5f));
3049 if (nnx==nx && nny==ny && nnz==nz)
return;
3053 if (nnx<nx||nny<ny||nnz<nz) {
3055 result=((
EMData *)image)->FourTruncate(nnx, nny, nnz, 1, 0);
3059 result=((
EMData *)image)->FourInterpol(nnx, nny, nnz, 1, 0);
3062 image->set_size(nnx,nny,nnz);
3063 memcpy(image->get_data(),result->get_data(),nnx*nny*nnz*
sizeof(
float));
3064 image->scale_pixel((
float)nx/(
float)nnx);
3075 int nx = from->get_xsize();
3076 int ny = from->get_ysize();
3077 int nz = from->get_zsize();
3079 int new_nx =
static_cast<int>(
static_cast<float> (nx) / sample_rate);
3080 int new_ny =
static_cast<int>(
static_cast<float> (ny) / sample_rate);
3081 int new_nz =
static_cast<int>(
static_cast<float> (nz) / sample_rate);
3083 if (new_nx == 0)
throw UnexpectedBehaviorException(
"The resample rate causes the pixel dimensions in the x direction to go to zero");
3084 if (new_ny == 0) new_ny = 1;
3085 if (new_nz == 0) new_nz = 1;
3087 int ndim = from->get_ndim();
3095 int fft_x_correction = 1;
3096 if (new_nx % 2 == 0) fft_x_correction = 2;
3098 int fft_y_correction = 0;
3099 if (ny != 1 && new_ny % 2 == 0 && ny % 2 == 1) fft_y_correction = 1;
3100 else if (ny != 1 && new_ny % 2 == 1 && ny % 2 == 0) fft_y_correction = -1;
3102 int fft_z_correction = 0;
3103 if (nz != 1 && new_nz % 2 == 0 && nz % 2 == 1) fft_z_correction = 1;
3104 else if (nz != 1 && new_nz % 2 == 1 && nz % 2 == 0) fft_z_correction = -1;
3106 if ( ! to->is_complex()) to->do_fft_inplace();
3108 if (ndim != 1) to->process_inplace(
"xform.fourierorigin.tocenter");
3110 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);
3113 if (fft_x_correction == 1) to->set_fftodd(
true);
3114 else to->set_fftodd(
false);
3116 if (ndim != 1) to->process_inplace(
"xform.fourierorigin.tocorner");
3118 to->do_ift_inplace();
3126 if (image->is_complex())
throw ImageFormatException(
"Error, the mean shrink processor does not work on complex images");
3128 if (image->get_ndim() == 1) {
throw ImageDimensionException(
"Error, mean shrink works only for 2D & 3D images"); }
3131 int shrink_factor = int(shrink_factor0);
3132 if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {
3134 "mean shrink: shrink factor must be >1 integer or 1.5");
3137 int nx = image->get_xsize();
3138 int ny = image->get_ysize();
3139 int nz = image->get_zsize();
3143 if (shrink_factor0==1.5 ) {
3144 if (nz > 1 )
throw InvalidValueException(shrink_factor0,
"mean shrink: only support 2D images for shrink factor = 1.5");
3146 int shrunken_nx = (int(nx / 1.5)+1)/2*2;
3147 int shrunken_ny = (int(ny / 1.5)+1)/2*2;
3148 EMData* result =
new EMData(shrunken_nx,shrunken_ny,1);
3156 int shrunken_nx = nx / shrink_factor;
3157 int shrunken_ny = ny / shrink_factor;
3158 int shrunken_nz = 1;
3161 shrunken_nz = nz / shrink_factor;
3165 EMData* result = image->copy_head();
3166 result->set_size(shrunken_nx,shrunken_ny,shrunken_nz);
3176 if (image->is_complex())
throw ImageFormatException(
"Error, the mean shrink processor does not work on complex images");
3178 if (image->get_ndim() == 1) {
throw ImageDimensionException(
"Error, mean shrink works only for 2D & 3D images"); }
3181 int shrink_factor = int(shrink_factor0);
3182 if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {
3184 "mean shrink: shrink factor must be >1 integer or 1.5");
3193 int nx = image->get_xsize();
3194 int ny = image->get_ysize();
3195 int nz = image->get_zsize();
3197 if (shrink_factor0==1.5 ) {
3198 if (nz > 1 )
throw InvalidValueException(shrink_factor0,
"mean shrink: only support 2D images for shrink factor = 1.5");
3200 int shrunken_nx = (int(nx / 1.5)+1)/2*2;
3201 int shrunken_ny = (int(ny / 1.5)+1)/2*2;
3203 EMData* orig = image->copy();
3204 image->set_size(shrunken_nx, shrunken_ny, 1);
3220 int shrunken_nx = nx / shrink_factor;
3221 int shrunken_ny = ny / shrink_factor;
3222 int shrunken_nz = 1;
3223 if (nz > 1) shrunken_nz = nz / shrink_factor;
3226 image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
3231 const float *
const data = from->get_const_data();
3232 float*
rdata = to->get_data();
3234 size_t nx = from->get_xsize();
3235 size_t ny = from->get_ysize();
3236 size_t nz = from->get_zsize();
3240 size_t shrunken_nx = nx / shrink_factor;
3241 size_t shrunken_ny = ny / shrink_factor;
3242 size_t shrunken_nz = 1;
3243 size_t shrunken_nxy = shrunken_nx * shrunken_ny;
3245 int normalize_shrink_factor = shrink_factor * shrink_factor;
3246 int z_shrink_factor = 1;
3249 shrunken_nz = nz / shrink_factor;
3250 normalize_shrink_factor *= shrink_factor;
3251 z_shrink_factor = shrink_factor;
3254 float invnormfactor = 1.0f/(float)normalize_shrink_factor;
3256 for (
size_t k = 0; k < shrunken_nz; k++) {
3257 size_t k_min = k * shrink_factor;
3258 size_t k_max = k * shrink_factor + z_shrink_factor;
3259 size_t cur_k = k * shrunken_nxy;
3261 for (
size_t j = 0; j < shrunken_ny; j++) {
3262 size_t j_min = j * shrink_factor;
3263 size_t j_max = j * shrink_factor + shrink_factor;
3264 size_t cur_j = j * shrunken_nx + cur_k;
3266 for (
size_t i = 0; i < shrunken_nx; i++) {
3267 size_t i_min = i * shrink_factor;
3268 size_t i_max = i * shrink_factor + shrink_factor;
3271 for (
size_t kk = k_min; kk < k_max; kk++) {
3272 size_t cur_kk = kk * nxy;
3274 for (
size_t jj = j_min; jj < j_max; jj++) {
3275 size_t cur_jj = jj * nx + cur_kk;
3276 for (
size_t ii = i_min; ii < i_max; ii++) {
3277 sum += data[ii + cur_jj];
3281 rdata[i + cur_j] = sum * invnormfactor;
3285 to->scale_pixel((
float)shrink_factor);
3291 int nx0 = from->get_xsize(), ny0 = from->get_ysize();
3293 int nx = to->get_xsize(), ny = to->get_ysize();
3295 float *data = to->get_data();
3296 const float *
const data0 = from->get_const_data();
3298 for (
int j = 0; j < ny; j++) {
3299 int jj = int(j * 1.5);
3300 float jw0 = 1.0F, jw1 = 0.5F;
3305 for (
int i = 0; i < nx; i++) {
3306 int ii = int(i * 1.5);
3307 float iw0 = 1.0F, iw1 = 0.5F;
3316 data[j * nx + i] = data0[ jj * nx0 + ii ] * jw0 * iw0 ;
3319 data[j * nx + i] += data0[ jj * nx0 + ii + 1] * jw0 * iw1;
3323 if ( jj +1 < ny0 ) {
3325 data[j * nx + i] += data0[ (jj+1) * nx0 + ii ] * jw1 * iw0;
3328 data[j * nx + i] += data0[ (jj+1) * nx0 + ii + 1] * jw1 * iw1;
3334 if ( w>0 ) data[j * nx + i] /= w;
3339 to->scale_pixel((
float)1.5);
3343template<
class LogicOp>
3360 int nz = image->get_zsize();
3361 int ny = image->get_ysize();
3362 int nx = image->get_xsize();
3364 if (nx == 1 && ny == 1 && nz == 1 )
return image->copy();
3369 int shrinkx = shrink;
3370 int shrinky = shrink;
3371 int shrinkz = shrink;
3373 int searchx = search;
3374 int searchy = search;
3375 int searchz = search;
3379 if ( shrinkx > nx ) shrinkx = nx;
3380 if ( shrinky > ny ) shrinky = ny;
3381 if ( shrinkz > nz ) shrinkz = nz;
3383 if ( nz == 1 && ny == 1 )
3385 return_image->set_size(nx/shrinkx);
3386 for(
int i = 0; i < nx/shrinkx; ++i)
3388 float tmp = op.get_start_val();
3389 for(
int s=0; s < searchx; ++s)
3391 int idx = shrinkx*i+s;
3393 if ( idx > nx )
break;
3396 float val = image->get_value_at(idx);
3397 if ( op( val,tmp) ) tmp = val;
3400 return_image->set_value_at(i,tmp);
3405 int ty = ny/shrinky;
3406 int tx = nx/shrinkx;
3407 return_image->set_size(tx,ty);
3408 for(
int y = 0;
y < ty; ++
y) {
3409 for(
int x = 0;
x < tx; ++
x) {
3410 float tmp = op.get_start_val();
3411 for(
int sy=0; sy < searchy; ++sy) {
3412 int yidx = shrinky*
y+sy;
3413 if ( yidx >= ny)
break;
3414 for(
int sx=0; sx < searchx; ++sx) {
3415 int xidx = shrinkx*
x+sx;
3416 if ( xidx >= nx)
break;
3418 float val = image->get_value_at(xidx,yidx);
3419 if ( op( val,tmp) ) tmp = val;
3422 return_image->set_value_at(
x,
y,tmp);
3428 int tz = nz/shrinkz;
3429 int ty = ny/shrinky;
3430 int tx = nx/shrinkx;
3432 return_image->set_size(tx,ty,tz);
3433 for(
int z = 0; z < tz; ++z) {
3434 for(
int y = 0;
y < ty; ++
y) {
3435 for(
int x = 0;
x < tx; ++
x) {
3436 float tmp = op.get_start_val();
3438 for(
int sz=0; sz < searchz; ++sz) {
3439 int zidx = shrinkz*z+sz;
3440 if ( zidx >= nz)
break;
3442 for(
int sy=0; sy < searchy; ++sy) {
3443 int yidx = shrinky*
y+sy;
3444 if ( yidx >= ny)
break;
3446 for(
int sx=0; sx < searchx; ++sx) {
3447 int xidx = shrinkx*
x+sx;
3448 if ( xidx >= nx)
break;
3449 float val = image->get_value_at(xidx,yidx,zidx);
3450 if ( op( val,tmp) ) tmp = val;
3454 return_image->set_value_at(
x,
y,z,tmp);
3459 return_image->update();
3461 return return_image;
3464template<
class LogicOp>
3480 int nz = image->get_zsize();
3481 int ny = image->get_ysize();
3482 int nx = image->get_xsize();
3486 int shrinkx = shrink;
3487 int shrinky = shrink;
3488 int shrinkz = shrink;
3490 int searchx = search;
3491 int searchy = search;
3492 int searchz = search;
3496 if ( shrinkx > nx ) shrinkx = nx;
3497 if ( shrinky > ny ) shrinky = ny;
3498 if ( shrinkz > nz ) shrinkz = nz;
3500 if (nx == 1 && ny == 1 && nz == 1 )
return;
3502 if ( nz == 1 && ny == 1 )
3504 for(
int i = 0; i < nx/shrink; ++i)
3506 float tmp = op.get_start_val();
3507 for(
int s=0; s < searchx; ++s)
3509 int idx = shrinkx*i+s;
3510 if ( idx > nx )
break;
3513 float val = image->get_value_at(idx);
3514 if ( op( val,tmp) ) tmp = val;
3517 image->set_value_at(i,tmp);
3520 image->set_size(nx/shrinkx);
3524 int ty = ny/shrinky;
3525 int tx = nx/shrinkx;
3526 for(
int y = 0;
y < ty; ++
y) {
3527 for(
int x = 0;
x < tx; ++
x) {
3528 float tmp = op.get_start_val();
3529 for(
int sy=0; sy < searchy; ++sy) {
3530 int yidx = shrinky*
y+sy;
3531 if ( yidx >= ny)
break;
3532 for(
int sx=0; sx < searchx; ++sx) {
3533 int xidx = shrinkx*
x+sx;
3534 if ( xidx >= nx)
break;
3536 float val = image->get_value_at(xidx,yidx);
3537 if ( op( val,tmp) ) tmp = val;
3540 (*image)(
x+tx*
y) = tmp;
3543 image->set_size(tx,ty);
3547 int tnxy = nx/shrinkx*ny/shrinky;
3548 int tz = nz/shrinkz;
3549 int ty = ny/shrinky;
3550 int tx = nx/shrinkx;
3552 for(
int z = 0; z < tz; ++z) {
3553 for(
int y = 0;
y < ty; ++
y) {
3554 for(
int x = 0;
x < tx; ++
x) {
3555 float tmp = op.get_start_val();
3556 for(
int sz=0; sz < searchz; ++sz) {
3557 int zidx = shrinkz*z+sz;
3558 if ( zidx >= nz)
break;
3559 for(
int sy=0; sy < searchy; ++sy) {
3560 int yidx = shrinky*
y+sy;
3561 if ( yidx >= ny)
break;
3562 for(
int sx=0; sx < shrinkx; ++sx) {
3563 int xidx = shrinkx*
x+sx;
3564 if ( xidx >= nx)
break;
3566 float val = image->get_value_at(xidx,yidx,zidx);
3567 if ( op( val,tmp) ) tmp = val;
3571 (*image)(
x+tx*
y+tnxy*z) = tmp;
3575 image->set_size(tx,ty,tz);
3589 int nz = image->get_zsize();
3591 LOGERR(
"%s Processor doesn't support 3D model",
get_name().c_str());
3595 int nx = image->get_xsize();
3596 int ny = image->get_ysize();
3597 float *dy =
new float[ny];
3601 float *data = image->get_data();
3603 for (
int i = 0; i < ny; i++) {
3609 float mean_y = sum_y / ny;
3613 for (
int j = 0; j < ny; j++) {
3614 for (
int i = 0; i < nx; i++) {
3615 data[i + j * nx] -= i * sum_x + j * mean_y + b;
3628 if (radius != 0 && mask != 0)
throw InvalidParameterException(
"Error - the mask and radius parameters are mutually exclusive.");
3630 if (mask == 0 && radius == 0)
throw InvalidParameterException(
"Error - you must specify either the mask or the radius parameter.");
3633 bool deletemask =
false;
3636 int n = image->get_ndim();
3638 mask->set_size(2*radius+1);
3640 mask->set_size(2*radius+1,2*radius+1);
3643 mask->set_size(2*radius+1,2*radius+1,2*radius+1);
3646 mask->process_inplace(
"testimage.circlesphere");
3650 int mnx = mask->get_xsize();
int mny = mask->get_ysize();
int mnz = mask->get_zsize();
3651 int nx = image->get_xsize();
int ny = image->get_ysize();
int nz = image->get_zsize();
3652 int nxc = nx+mnx;
int nyc = ny+mny;
int nzc = nz+mnz;
3653 if (nz == 1) nzc = 1;
3654 if (ny == 1) nyc = 1;
3656 if ( mnx > nx || mny > ny || mnz > nz)
3660 float normfac = 0.0;
3661 for (
int i=0; i<mask->get_xsize()*mask->get_ysize()*mask->get_zsize(); ++i){
3662 normfac += mask->get_value_at(i);
3667 if (normfac == 0)
throw InvalidParameterException(
"Error - the pixels in the mask sum to zero. This breaks the flattening procedure");
3668 normfac = 1.0f/normfac;
3674 if (ny == 1) r =
Region((mnx-nxc)/2,nxc);
3675 else if (nz == 1) r =
Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);
3676 else r =
Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);
3686 if (ny == 1) r2 =
Region((nx-nxc)/2,nxc);
3687 else if (nz == 1) r2 =
Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);
3688 else r2 =
Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);
3695 m->process_inplace(
"xform.phaseorigin.tocenter");
3706 if (ny == 1) r =
Region((nxc-mnx)/2,mnx);
3707 else if (nz == 1) r =
Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);
3708 else r =
Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);
3713 if (ny == 1) r3 =
Region((nxc-nx)/2,nx);
3714 else if (nz == 1) r3 =
Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);
3715 else r3 =
Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);
3725 if (!image) {
LOGWARN(
"NULL IMAGE");
return; }
3729 if (image->get_ysize()==1) {
3733 else if (image->get_zsize()==1) {
3735 int nx2=image->get_xsize()/2;
3736 int ny2=image->get_ysize()/2;
3739 for (
int i=1; i<nx2; i++) {
3740 if (rdist[i]>rdist[i-1]) rdist[i]=rdist[i-1];
3743 image->process_inplace(
"xform.fourierorigin.tocenter");
3744 EMData* binary=image->copy();
3747 for (
int x=0;
x<image->get_xsize();
x+=2) {
3748 for (
int y=1;
y<image->get_ysize()-1;
y++) {
3749 int r=(int)hypot((
float)(
x/2),(
float)(
y-ny2));
3750 float cen=(*binary)(
x,
y);
3751 if (
x==0 ||
x==nx2*2-2 || (cen>(*binary)(
x+2,
y) || cen>(*binary)(
x-2,
y) || cen>(*binary)(
x,
y+1) || cen >(*binary)(
x,
y-1) || (*binary)(
x,
y)>rdist[r])) {
3752 image->set_value_at_fast(
x/2+nx2,
y,0.0);
3753 image->set_value_at_fast(nx2-
x/2,ny2*2-
y-1,0.0);
3756 image->set_value_at_fast(
x/2+nx2,
y,cen);
3757 image->set_value_at_fast(nx2-
x/2,ny2*2-
y-1,cen);
3761 image->set_value_at_fast(nx2+1,ny2,(*binary)(2,ny2));
3762 image->set_value_at_fast(nx2-1,ny2,(*binary)(2,ny2));
3763 image->set_value_at_fast(nx2,ny2+1,(*binary)(0,ny2+1));
3764 image->set_value_at_fast(nx2,ny2-1,(*binary)(0,ny2-1));
3765 for (
int y=0;
y<ny2*2;
y++) image->set_value_at_fast(0,
y,0.0f);
3768 float *idat=image->get_data();
3769 float *bdat=binary->get_data();
3770 int nxy=(nx2*ny2*4);
3771 for (
int i=0; i<nxy; i++) {
3772 bdat[i]=idat[i]==0?0:1.0f;
3777 image->set_complex(
false);
3778 binary->set_complex(
false);
3781 image->set_fftpad(
false);
3782 binary->set_fftpad(
false);
3785 image->process_inplace(
"filter.lowpass.gauss",
Dict(
"cutoff_abs",0.04f));
3786 binary->process_inplace(
"filter.lowpass.gauss",
Dict(
"cutoff_abs",0.04f));
3791 for (
int x=0;
x<image->get_xsize();
x+=2) {
3792 for (
int y=0;
y<image->get_ysize();
y++) {
3793 float bv=binary->get_value_at(
x/2+nx2,
y);
3794 image->set_value_at_fast(
x,
y,image->get_value_at(
x/2+nx2,
y)/(bv<=0?1.0f:bv));
3795 image->set_value_at_fast(
x+1,
y,0.0);
3798 image->set_complex(
true);
3799 image->set_fftpad(
true);
3800 image->process_inplace(
"xform.fourierorigin.tocorner");
3808#include <gsl/gsl_linalg.h>
3816 int nz = image->get_zsize();
3818 LOGERR(
"%s Processor doesn't support 3D model",
get_name().c_str());
3822 int nx = image->get_xsize();
3823 int ny = image->get_ysize();
3824 float *d = image->get_data();
3829 if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) {
3830 LOGERR(
"%s Processor requires same size mask image",
get_name().c_str());
3833 dm = mask->get_data();
3837 for(
int i=0; i<nx*ny; i++) {
3845 LOGERR(
"%s Processor requires at least 3 pixels to fit a plane",
get_name().c_str());
3849 gsl_vector *S=gsl_vector_calloc(3);
3850 gsl_matrix *A=gsl_matrix_calloc(count,3);
3851 gsl_matrix *
V=gsl_matrix_calloc(3,3);
3853 double m[3] = {0, 0, 0};
3856 for(
int j=0; j<ny; j++){
3857 for(
int i=0; i<nx; i++){
3871 for(
int j=0; j<ny; j++){
3872 for(
int i=0; i<nx; i++){
3884 for(
int i=0; i<3; i++) m[i]/=count;
3888 for(
int j=0; j<ny; j++){
3889 for(
int i=0; i<nx; i++){
3893 gsl_matrix_set(A,index,0,i-m[0]);
3894 gsl_matrix_set(A,index,1,j-m[1]);
3895 gsl_matrix_set(A,index,2,d[ij]-m[2]);
3903 for(
int j=0; j<ny; j++){
3904 for(
int i=0; i<nx; i++){
3907 gsl_matrix_set(A,index,0,i-m[0]);
3908 gsl_matrix_set(A,index,1,j-m[1]);
3909 gsl_matrix_set(A,index,2,d[ij]-m[2]);
3916 gsl_linalg_SV_decomp_jacobi(A,
V, S);
3919 for(
int i=0; i<3; i++) n[i] = gsl_matrix_get(
V, i, 2);
3922 printf(
"S=%g,%g,%g\n",gsl_vector_get(S,0), gsl_vector_get(S,1), gsl_vector_get(S,2));
3923 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));
3924 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));
3925 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));
3926 printf(
"Fitted plane: p0=%g,%g,%g\tn=%g,%g,%g\n",m[0],m[1],m[2],n[0],n[1],n[2]);
3932 for(
int j=0; j<nx; j++){
3933 for(
int i=0; i<ny; i++){
3935 d[ij]-=
static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);
3940 for(
int j=0; j<nx; j++){
3941 for(
int i=0; i<ny; i++){
3943 if(d[ij]) d[ij]-=
static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);
3949 vector< float > planeParam;
3950 planeParam.resize(6);
3951 for(
int i=0; i<3; i++) planeParam[i] = static_cast<float>(n[i]);
3952 for(
int i=0; i<3; i++) planeParam[i+3] = static_cast<float>(m[i]);
3963 int nx = image->get_xsize();
3964 int ny = image->get_ysize();
3965 int nz = image->get_zsize();
3967 float *data = image->get_data();
3968 float sigma = image->get_attr(
"sigma");
3970 for (
int k = 0; k < nz; k++) {
3971 for (
int i = 0; i < nx; i++) {
3973 for (
int j = ny / 4; j < 3 * ny / 4; j++) {
3974 sum += data[i + j * nx];
3977 float mean = (float)sum / (ny / 2);
3978 for (
int j = 0; j < ny; j++) {
3979 data[i + j * nx] = (data[i + j * nx] - mean) / sigma;
3995 if(image->is_complex()) {
3996 LOGERR(
"%s Processor only operates on real images",
get_name().c_str());
4001 int nz = image->get_zsize();
4003 LOGERR(
"%s Processor doesn't support 3D models",
get_name().c_str());
4007 EMData *ff=image->do_fft();
4010 int nx=image->get_xsize();
4011 int ny=image->get_ysize();
4014 float norm=
static_cast<float>(nx*ny);
4016 for (
y=0;
y<ny;
y++) image->set_value_at(0,
y,0);
4018 for (
x=1;
x<nx/2;
x++) {
4019 for (
y=0;
y<ny;
y++) {
4021 if (
y<ny/2) y2=
y+ny/2;
4022 else if (
y==ny/2) y2=ny;
4024 image->set_value_at(
x,
y,ff->get_value_at(nx-
x*2,ny-y2)/norm);
4028 for (
x=nx/2;
x<nx;
x++) {
4029 for (
y=0;
y<ny;
y++) {
4031 if (
y<ny/2) y2=
y+ny/2;
4033 image->set_value_at(
x,
y,ff->get_value_at(
x*2-nx,y2)/norm);
4052 if (image->get_zsize() > 1) {
4053 LOGERR(
"%s Processor doesn't support 3D model",
get_name().c_str());
4059 float *d = image->get_data();
4063 int nx = image->get_xsize();
4064 int ny = image->get_ysize();
4066 float zval=9.99e23f;
4069 size_t corn=nx*ny-1;
4072 for (
x=1;
x<nx;
x++) {
if (d[
x]!=d[0])
break;}
4073 if (
x==nx) zval=d[0];
4075 for (
y=1;
y<ny;
y++) {
if (d[
y*nx]!=d[0])
break; }
4076 if (
y==ny) zval=d[0];
4078 for (
x=1;
x<nx;
x++) {
if (d[corn-
x]!=d[corn])
break;}
4079 if (
x==nx) zval=d[corn];
4081 for (
y=1;
y<ny;
y++) {
if (d[corn-
y*nx]!=d[corn])
break; }
4082 if (
y==ny) zval=d[corn];
4084 if (zval!=9.99e23f) {
4085 image->set_attr(
"hadzeroedge",1);
4089 else image->set_attr(
"hadzeroedge",0);
4092 for (
x=nx/2-5;
x<nx/2+5;
x++) {
4093 if (d[
x]!=d[
x+nx] || d[
x]!=d[
x+nx*2] )
break;
4095 if (
x==nx/2+5) image->set_attr(
"hadzeroedge",2);
4097 for (
x=nx/2-5;
x<nx/2+5;
x++) {
4098 if (d[corn-
x]!=d[corn-
x-nx] || d[corn-
x]!=d[corn-
x-nx*2])
break;
4100 if (
x==nx/2+5) image->set_attr(
"hadzeroedge",2);
4102 for (
y=ny/2-5;
y<ny/2+5;
y++) {
4103 if (d[
y*nx]!=d[
y*nx+1] || d[
y*nx]!=d[
y*nx+2] )
break;
4105 if (
y==ny/2+5) image->set_attr(
"hadzeroedge",2);
4107 for (
y=ny/2-5;
y<ny/2+5;
y++) {
4108 if (d[corn-
y*nx]!=d[corn-
y*nx-1] || d[corn-
y*nx]!=d[corn-
y*nx-2])
break;
4110 if (
y==ny/2+5) image->set_attr(
"hadzeroedge",2);
4113 if (zval==9.99e23f) zval=0;
4115 for (j = 0; j < ny; j++) {
4116 for (i = 0; i < nx - 1; i++) {
4117 if (d[i + j * nx] != zval) {
4122 float v = d[i + j * nx];
4128 for (i = nx - 1; i > 0; i--) {
4129 if (d[i + j * nx] != zval)
4139 for (i = 0; i < nx; i++) {
4140 for (j = 0; j < ny; j++) {
4141 if (d[i + j * nx] != zval)
4145 float v = d[i + j * nx];
4151 for (j = ny - 1; j > 0; j--) {
4152 if (d[i + j * nx] != zval)
4176 if (sigmamult<=0.0)
throw InvalidValueException(sigmamult,
"threshold.outlier.localmean: sigma must be >0");
4178 float hithr=(float)image->get_attr(
"mean")+(float)(image->get_attr(
"sigma"))*sigmamult;
4179 float lothr=(float)image->get_attr(
"mean")-(float)(image->get_attr(
"sigma"))*sigmamult;
4181 int nx=image->get_xsize();
4182 int ny=image->get_ysize();
4183 int nz=image->get_zsize();
4188 im[1]=image->copy_head();
4193 memcpy(im[1]->
get_data(),im[0]->
get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*
sizeof(
float));
4195 for (
int y=0;
y<ny;
y++) {
4196 for (
int x=0;
x<nx;
x++) {
4198 float pix=im[1]->get_value_at(
x,
y);
4199 if (pix>hithr || pix<lothr || (pix==0 && fix_zero)) {
4201 int y1=
y>=ny-1?ny-1:
y+1;
4203 int x1=
x>=nx-1?nx-1:
x+1;
4204 float c=0.0f,nc=0.0f;
4205 for (
int yy=y0; yy<=y1; yy++) {
4206 for (
int xx=x0; xx<=x1; xx++) {
4207 float lpix=im[1]->get_value_at(xx,yy);
4208 if (lpix>hithr || lpix<lothr || (lpix==0 && fix_zero))
continue;
4213 if (nc!=0) im[0]->set_value_at(
x,
y,c/nc);
4233 if (image->get_zsize() > 1) {
4234 LOGERR(
"BeamstopProcessor doesn't support 3D model");
4238 float value1 =
params[
"value1"];
4239 float value2 =
params[
"value2"];
4240 float value3 =
params[
"value3"];
4242 float thr = fabs(value1);
4243 float *data = image->get_data();
4244 int cenx = (int) value2;
4245 int ceny = (int) value3;
4247 int nx = image->get_xsize();
4248 int ny = image->get_ysize();
4258 int mxr = (int) floor(
sqrt(2.0f) * nx / 2);
4260 float *mean_values =
new float[mxr];
4261 float *sigma_values =
new float[mxr];
4264 double square_sum = 0;
4266 for (
int i = 0; i < mxr; i++) {
4270 int nitems = 6 * i + 2;
4272 for (
int j = 0; j < nitems; j++) {
4273 float ang = j * 2 * M_PI / nitems;
4274 int x0 = (int) floor(cos(ang) * i + cenx);
4275 int y0 = (int) floor(sin(ang) * i + ceny);
4277 if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny) {
4281 float f = data[x0 + y0 * nx];
4283 square_sum += f * f;
4287 mean_values[i] = (float)sum / count;
4288 sigma_values[i] = (float)
sqrt(square_sum / count - mean_values[i] * mean_values[i]);
4292 for (
int k = 0; k < 5; k++) {
4293 for (
int i = 0; i < mxr; i++) {
4297 int nitems = 6 * i + 2;
4298 double thr1 = mean_values[i] - sigma_values[i] * thr;
4299 double thr2 = mean_values[i] + sigma_values[i];
4301 for (
int j = 0; j < nitems; j++) {
4302 float ang = j * 2 * M_PI / nitems;
4303 int x0 = (int) floor(cos(ang) * i + cenx);
4304 int y0 = (int) floor(sin(ang) * i + ceny);
4306 if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny ||
4307 data[x0 + y0 * nx] < thr1 || data[x0 + y0 * nx] > thr2) {
4311 sum += data[x0 + y0 * nx];
4312 square_sum += data[x0 + y0 * nx] * data[x0 + y0 * nx];
4316 mean_values[i] = (float) sum / count;
4317 sigma_values[i] = (float)
sqrt(square_sum / count - mean_values[i] * mean_values[i]);
4321 for (
int i = 0; i < nx; i++) {
4322 for (
int j = 0; j < ny; j++) {
4324 int r =
Util::round(hypot((
float) i - cenx, (
float) j - ceny));
4327 if (data[i + j * nx] < (mean_values[r] - sigma_values[r] * thr)) {
4328 data[i + j * nx] = 0;
4331 data[i + j * nx] -= mean_values[r];
4335 if (data[i + j * nx] > (mean_values[r] - sigma_values[r] * thr)) {
4338 data[i + j * nx] = mean_values[r];
4344 delete[]mean_values;
4350 delete[]sigma_values;
4365 if (image->get_zsize() > 1) {
4366 LOGERR(
"MeanZeroEdgeProcessor doesn't support 3D model");
4370 int nx = image->get_xsize();
4371 int ny = image->get_ysize();
4372 Dict dict = image->get_attr_dict();
4373 float mean_nonzero = dict.
get(
"mean_nonzero");
4375 float *d = image->get_data();
4379 for (j = 0; j < ny; j++) {
4380 for (i = 0; i < nx - 1; i++) {
4381 if (d[i + j * nx] != 0) {
4390 float v = d[i + j * nx] - mean_nonzero;
4394 d[i + j * nx] = v + mean_nonzero;
4399 for (i = nx - 1; i > 0; i--) {
4400 if (d[i + j * nx] != 0) {
4409 v = d[i + j * nx] - mean_nonzero;
4413 d[i + j * nx] = v + mean_nonzero;
4419 for (i = 0; i < nx; i++) {
4420 for (j = 0; j < ny; j++) {
4421 if (d[i + j * nx] != 0)
4425 float v = d[i + j * nx] - mean_nonzero;
4429 d[i + j * nx] = v + mean_nonzero;
4433 for (j = ny - 1; j > 0; j--) {
4434 if (d[i + j * nx] != 0)
4438 v = d[i + j * nx] - mean_nonzero;
4442 d[i + j * nx] = v + mean_nonzero;
4459 float *data = image->get_data();
4460 int nx = image->get_xsize();
4461 int ny = image->get_ysize();
4462 int nz = image->get_zsize();
4463 size_t nxy = (size_t)nx * ny;
4466 for (
int z = 0; z < nz; z++) {
4467 for (
int x = 0;
x < nx;
x++) {
4469 for (
int y = 0;
y < ny;
y++) {
4470 idx =
x +
y * nx + z * nxy;
4473 float mean = (float) sum / ny;
4475 for (
int y = 0;
y < ny;
y++) {
4476 idx =
x +
y * nx + z * nxy;
4492 int nx = image->get_xsize();
4493 int ny = image->get_ysize();
4495 float *d = image->get_data();
4496 int width =
params[
"width"];
4498 if (width > min(nx,ny)/2.){
4499 LOGERR(
"width parameter cannot be greater than min(nx,ny)/2");
4503 if (image->get_zsize() > 1){
4504 for (
int k=0; k<image->get_zsize(); k++){
4506 for (
int i=0; i<width; i++) {
4507 float frac=i/(float)width;
4508 for (
int j=0; j<nx; j++) {
4509 d[zidx+j+i*nx]*=frac;
4510 d[zidx+nx*ny-j-i*nx-1]*=frac;
4512 for (
int j=0; j<ny; j++) {
4513 d[zidx+j*nx+i]*=frac;
4514 d[zidx+nx*ny-j*nx-i-1]*=frac;
4521 for (
int i=0; i<width; i++) {
4522 float frac=i/(float)width;
4523 for (
int j=0; j<nx; j++) {
4525 d[nx*ny-j-i*nx-1]*=frac;
4527 for (
int j=0; j<ny; j++) {
4529 d[nx*ny-j*nx-i-1]*=frac;
4544 if (image->get_zsize() > 1) {
4545 LOGERR(
"ZeroEdgeRowProcessor is not supported in 3D models");
4549 int nx = image->get_xsize();
4550 int ny = image->get_ysize();
4552 float *d = image->get_data();
4564 size_t row_size = nx *
sizeof(float);
4566 memset(d, 0, top_nrows * row_size);
4567 memset(d + (ny - bottom_nrows) * nx, 0, bottom_nrows * row_size);
4569 for (
int i = top_nrows; i < ny - bottom_nrows; i++) {
4570 memset(d + i * nx, 0, left_ncols *
sizeof(
float));
4571 memset(d + i * nx + nx - right_ncols, 0, right_ncols *
sizeof(
float));
4583 if (image->get_zsize() <= 1) {
4584 LOGERR(
"ZeroEdgePlaneProcessor only support 3D models");
4588 int nx = image->get_xsize();
4589 int ny = image->get_ysize();
4590 int nz = image->get_zsize();
4592 float *d = image->get_data();
4601 size_t row_size = nx *
sizeof(float);
4602 size_t nxy = nx * ny;
4603 size_t sec_size = nxy *
sizeof(float);
4604 size_t y0row = y0 * row_size;
4605 size_t y1row = y1 * row_size;
4607 size_t x0size = x0*
sizeof(float);
4608 size_t x1size = x1*
sizeof(float);
4610 memset(d,0,z0*sec_size);
4611 memset(d+(nxy*(nz-z1)),0,sec_size*z1);
4613 for (
int z=z0; z<nz-z1; z++) {
4614 memset(d+z*nxy,0,y0row);
4615 memset(d+z*nxy+(ny-y1)*nx,0,y1row);
4618 int znxy2 = znxy + nx - x1;
4620 for (
int y=y0;
y<max_y;
y++) {
4621 memset(d+znxy+
y*nx,0,x0size);
4622 memset(d+znxy2+
y*nx,0,x1size);
4632 return image->get_attr(
"sigma");
4638 LOGWARN(
"cannot do normalization on NULL image");
4642 if (image->is_complex()) {
4643 LOGWARN(
"cannot do normalization on complex image");
4649 LOGWARN(
"cannot do normalization on image with sigma = 0");
4655 size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
4656 float *data = image->get_data();
4658 for (
size_t i = 0; i < size; ++i) {
4659 data[i] = (data[i] - mean) / sigma;
4671 float ret=
sqrt((
float)image->get_attr(
"square_sum"));
4672 return ret==0.0f?1.0f:ret;
4681 float ret=(float)image->get_attr(
"mean")*image->get_xsize()*image->get_ysize()*image->get_zsize();
4682 return ret==0.0f?1.0f:ret;
4695 int nx=image->get_xsize();
4696 int ny=image->get_ysize();
4697 int nz=image->get_zsize();
4703 for (
int z=0; z<nz; z++) {
4704 for (
int y=0;
y<ny;
y++) {
4705 for (
int x=0;
x<nx;
x++) {
4706 double m=mask->get_value_at(
x,
y,z);
4708 if (!apply_mask) m=1.0;
4709 double v=(double)image->get_value_at(
x,
y,z)*m;
4717 float mean=(float)(sum/nnz);
4718 float sigma=
sqrt((
float)(sumsq-sum*sum/nnz)/nnz);
4719 if (no_sigma) sigma=1.0f;
4721 for (
int z=0; z<nz; z++) {
4722 for (
int y=0;
y<ny;
y++) {
4723 for (
int x=0;
x<nx;
x++) {
4724 float m=mask->get_value_at(
x,
y,z);
4725 if (!apply_mask) m=1.0;
4726 if (m==0) image->set_value_at(
x,
y,z,0);
4727 else image->set_value_at(
x,
y,z,(image->get_value_at(
x,
y,z)*m-mean)/sigma);
4736 LOGWARN(
"cannot do normalization on NULL image");
4740 if (image->is_complex()) {
4741 LOGWARN(
"cannot do normalization on complex image");
4745 image->process_inplace(
"filter.ramp" );
4746 int nx = image->get_xsize();
4748 mask.process_inplace(
"testimage.circlesphere",
Dict(
"radius",nx/2-2,
"fill",1));
4750 vector<float> rstls = Util::infomask( image, &mask,
false);
4751 image->add((
float)-rstls[0]);
4752 image->mult((
float)1.0/rstls[1]);
4763 float tthr =
params.
set_default(
"thr",(
float)image->get_attr(
"mean")+(
float)image->get_attr(
"sigma"));
4765 float apix = image->get_attr_default(
"apix_x",1.0f);
4770 float step = ((float)image->get_attr(
"sigma"))/5.0f;
4775 size_t n = image->get_size();
4776 float* d = image->get_data();
4779 float thr=(float)image->get_attr(
"mean")+(float)image->get_attr(
"sigma")/2.0;
4781 for (
size_t i=0; i<n; ++i) {
4782 if (d[i]>=thr) ++count;
4784 if (verbose) printf(
"apix=%1.3f\tmass=%1.1f\tthr=%1.2f\tstep=%1.3g\n",apix,mass,thr,step);
4786 float max = image->get_attr(
"maximum");
4787 float min = image->get_attr(
"minimum");
4788 for (
int j=0; j<4; j++) {
4790 while (thr<max && count*apix*apix*apix*.81/1000.0>mass) {
4793 for (
size_t i=0; i<n; ++i) {
4794 if (d[i]>=thr) ++count;
4798 if (verbose>1) printf(
"%d\t%d\t%1.3f\t%1.2f\n",err,count,thr,count*apix*apix*apix*.81/1000.0);
4803 while (thr>min && count*apix*apix*apix*.81/1000.0<mass) {
4806 for (
size_t i=0; i<n; ++i) {
4807 if (d[i]>=thr) ++count;
4811 if (verbose>1) printf(
"%d\t%d\t%1.3f\t%1.2f\n",err,count,thr,count*apix*apix*apix*.81/1000.0);
4818 if ((
float)tthr/thr>0) image->mult((
float)tthr/thr);
4819 else printf(
"WARNING: could not normalize map to specified mass.");
4829 return image->get_edge_mean();
4839 int nx=image->get_xsize();
4840 int ny=image->get_ysize();
4841 int nz=image->get_zsize();
4845 if (radius<0) radius=ny/2+radius;
4847 static bool busy =
false;
4849 static int oldradius=radius;
4850 static int oldwidth=width;
4858 mask->set_size(nx, ny, nz);
4861 mask->process_inplace(
"mask.sharp",
Dict(
"inner_radius", radius,
4862 "outer_radius", radius + width));
4867 float *d = mask->get_data();
4868 float * data = image->get_data();
4869 size_t size = (size_t)nx*ny*nz;
4870 for (
size_t i = 0; i < size; ++i) {
4871 if (d[i]!=0.0f) { n+=1.0; s+=data[i]; }
4874 float result = (float)(s/n);
4888 float maxval = image->get_attr(
"maximum");
4889 float minval = image->get_attr(
"minimum");
4890 return (maxval + minval) / 2;
4899 float maxval = image->get_attr(
"maximum");
4900 float minval = image->get_attr(
"minimum");
4901 return (maxval - minval) / 2;
4911 int nx = image->get_xsize();
4912 int ny = image->get_ysize();
4913 int nz = image->get_zsize();
4914 float *d = image->get_data();
4915 size_t nyz = ny * nz;
4917 for (
size_t i = 0; i < nyz; i++) {
4919 size_t r = l + nx - 2;
4920 sum += d[l] + d[l + 1] + d[r] + d[r + 1];
4922 float mean = (float) sum / (4 * nyz);
4936 float *
rdata = image->get_data();
4937 int nx = image->get_xsize();
4938 int ny = image->get_ysize();
4939 int nz = image->get_zsize();
4944 for (
int z = 0; z < nz; z++) {
4947 for (
int y = 0;
y < ny;
y++) {
4948 for (
int x = 0;
x < nx;
x++) {
4949 float d =
rdata[
x +
y * nx + z*nx*ny];
4955 double row_mean = row_sum / nx / ny;
4956 double row_std = row_sqr / nx / ny;
4957 row_std =
sqrt(row_std - row_mean*row_mean);
4958 for (
int y = 0;
y < ny;
y++) {
4959 for (
int x = 0;
x < nx;
x++) {
4960 rdata[
x +
y * nx+ z*nx*ny] -= (float)row_mean;
4961 rdata[
x +
y * nx+ z*nx*ny] /= (float)row_std;
4974 for (
int y = 0;
y < ny;
y++) {
4976 for (
int x = 0;
x < nx;
x++) {
4977 row_len += pow((
double)
rdata[
x +
y * nx],2.0);
4979 row_len=
sqrt(row_len);
4980 if (row_len==0) row_len=1.0;
4982 for (
int x = 0;
x < nx;
x++) {
4983 rdata[
x +
y * nx] /= (float)row_len;
4990 for (
int y = 0;
y < ny;
y++) {
4992 for (
int x = 0;
x < nx;
x++) {
4996 double row_mean = row_sum / nx;
4997 if (row_mean <= 0) {
5001 for (
int x = 0;
x < nx;
x++) {
5002 rdata[
x +
y * nx] /= (float)row_mean;
5015 return image->get_attr(
"mean");
5024 float mean = image->get_attr(
"mean");
5025 float sig = image->get_attr(
"sigma");
5026 size_t n=image->get_size();
5029 int hist_size = image->get_size()/100;
5030 if (hist_size>1000) hist_size=1000;
5031 if (hist_size<5) hist_size=5;
5032 vector<float> hist(hist_size,0.0f);
5034 float histmin=mean-2.0f*sig;
5035 float histmax=mean+2.0f*sig;
5036 float step=(histmax-histmin)/(hist_size-1);
5037 for (
size_t i=0; i<n; i++) {
5038 int bin = (image->get_value_at_index(i)-histmin)/step;
5039 bin=bin<0?0:(bin>=hist_size?hist_size-1:bin);
5044 int maxv=hist[1],maxn=1;
5045 for (
int i=2; i<hist_size-1; i++) {
5046 if (hist[i]>maxv) { maxv=hist[i]; maxn=i; }
5050 float v1=hist[maxn-1],v2=hist[maxn],v3=hist[maxn+1];
5051 float ret= (v1*(maxn-1)+v2*maxn+v3*(maxn+1))/(v1+v2+v3)*step+histmin;
5052 printf(
"%f %f %f %f %f %f\n",v1,v2,v3,mean,sig,ret);
5074 int si0=(int)floor(
params.
set_default(
"low_cutoff_frequency",0.0f)*image->get_ysize());
5075 int si1=(int)ceil(
params.
set_default(
"high_cutoff_frequency",0.7071f)*image->get_ysize());
5079 if (image->is_complex()) imf=image->copy();
5080 else imf=image->do_fft();
5083 EMData *ctfi=imf->copy_head();
5086 ctf=(
Ctf *)(image->get_attr(
"ctf"));
5094 if (refr->is_complex()) ref=refr;
5095 else ref=refr->do_fft();
5098 if (actual==NULL) actf=ref;
5100 if (ctfweight)
throw InvalidCallException(
"math.sub.optimal: Sorry, cannot use ctfweight in combination with actual");
5101 if (actual->is_complex()) actf=actual;
5102 else actf=actual->do_fft();
5105 int ny2=(int)(image->get_ysize()*
sqrt(2.0)/2);
5106 vector <double>rad(ny2+1);
5107 vector <double>norm(ny2+1);
5112 for (
int y=-ny2;
y<ny2;
y++) {
5113 for (
int x=0;
x<ny2;
x++) {
5115 if (r>ny2)
continue;
5116 std::complex<float> v1=imf->get_complex_at(
x,
y);
5117 std::complex<float> v2=ref->get_complex_at(
x,
y);
5118 rad[r]+=(double)(v1.real()*v2.real()+v1.imag()*v2.imag());
5120 norm[r]+=(double)(v2.real()*v2.real()+v2.imag()*v2.imag());
5123 for (
int i=1; i<ny2; i++) rad[i]/=norm[i];
5132 if (return_presigma) {
5133 for (
int y=-ny2;
y<ny2;
y++) {
5134 for (
int x=0;
x<imf->get_xsize()/2;
x++) {
5136 if (r>=ny2 || r>=si1 || r<si0) {
5137 imf->set_complex_at(
x,
y,0);
5140 std::complex<float> v1=imf->get_complex_at(
x,
y);
5141 imf->set_complex_at(
x,
y,v1);
5144 EMData *tmp=imf->do_ift();
5145 oldsig=(float)tmp->get_attr(
"sigma");
5149 for (
int y=-ny2;
y<ny2;
y++) {
5150 for (
int x=0;
x<imf->get_xsize()/2;
x++) {
5152 if (r>=ny2 || r>=si1 || r<si0) {
5153 imf->set_complex_at(
x,
y,0);
5156 std::complex<float> v1=imf->get_complex_at(
x,
y);
5157 std::complex<float> v2=actf->get_complex_at(
x,
y);
5159 if (return_subim) imf->set_complex_at(
x,
y,v2);
5160 else imf->set_complex_at(
x,
y,v1-v2);
5164 if (!refr->is_complex())
delete ref;
5165 if (actual!=NULL && !actual->is_complex())
delete actf;
5168 if (return_radial) {
5170 for (
int i=0; i<ny2; i++) radf[i]=(
float)rad[i];
5174 EMData *ret=imf->do_ift();
5176 if (return_radial) ret->set_attr(
"filter_curve",radf);
5177 if (return_presigma) {
5178 ret->set_attr(
"sigma_presub",oldsig);
5183 if (return_radial) imf->set_attr(
"filter_curve",radf);
5184 if (return_presigma) imf->set_attr(
"sigma_presub",oldsig);
5196 memcpy(image->get_data(),tmp->get_data(),(
size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*
sizeof(
float));
5219 int nx = image->get_xsize();
5220 int ny = image->get_ysize();
5221 int nz = image->get_zsize();
5222 size_t size = (size_t)nx * ny * nz;
5231 float meani=(float)image->get_attr(
"mean");
5232 float meant=(float)to->get_attr(
"mean");
5233 float sigi=(float)image->get_attr(
"sigma")*ignore_lowsig;
5234 float sigt=(float)to->get_attr(
"sigma")*ignore_lowsig;
5237 fim=image->do_fft();
5241 dimage=fim->get_data();
5242 dto=fto->get_data();
5244 size2=(size_t)(nx+2) * ny * nz;
5249 for (
size_t i=0; i<size2; i+=2) { sigi+=pow(dimage[i],2.0f); sigt+=pow(dto[i],2.0f); }
5250 sigi=ignore_lowsig*
sqrt(sigi/(size2/2));
5251 sigt=ignore_lowsig*
sqrt(sigt/(size2/2));
5254 dimage = image->get_data();
5255 dto = to->get_data();
5262 if (ignore_lowsig<0) ignore_lowsig=0;
5265 for (
size_t i = 0; i < size2; i+=step) {
5266 if (dto[i] >= low_threshold && dto[i] <= high_threshold
5267 && (dto[i]>=meant+sigt || dto[i]<=meant-sigt)
5268 && (dimage[i]>=meani+sigi || dimage[i]<=meani-sigi)
5269 && (!ignore_zero ||(dto[i] != 0.0f && dimage[i] != 0.0f))) {
5275 double *
x=(
double *)malloc(count*
sizeof(
double));
5276 double *
y=(
double *)malloc(count*
sizeof(
double));
5278 for (
size_t i = 0; i < size2; i+=step) {
5279 if (dto[i] >= low_threshold && dto[i] <= high_threshold
5280 && (dto[i]>=meant+sigt || dto[i]<=meant-sigt)
5281 && (dimage[i]>=meani+sigi || dimage[i]<=meani-sigi)
5282 && (!ignore_zero ||(dto[i] != 0.0f && dimage[i] != 0.0f))) {
5289 double cov00,cov01,cov11,sumsq;
5290 gsl_fit_linear (
x, 1,
y, 1, count, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
5293 FILE*out=fopen(
"debug.txt",
"w");
5294 for (
size_t i = 0; i < count; i++) {
5295 fprintf(out,
"%lf\t%lf\n",
x[i],
y[i]);
5298 printf(
"add %lf\tmul %lf\t%lf\t%lf\t%lf\t%lf\n",c0,c1,cov00,cov01,cov11,sumsq);
5303 if (fim!=NULL)
delete fim;
5304 if (fto!=NULL)
delete fto;
5306 if (fourieramp) c0=0;
5307 dimage = image->get_data();
5308 for (
size_t i = 0; i < size; ++i) dimage[i]=dimage[i]*c1+c0;
5310 image->set_attr(
"norm_mult",c1);
5311 image->set_attr(
"norm_add",c0);
5316 if (!image->is_complex())
throw ImageFormatException(
"Fourier binary thresholding processor only works for complex images");
5321 float* d = image->get_data();
5322 for(
size_t i = 0; i < image->get_size()/2; ++i, d+=2) {
5323 if ( *d < threshold ) {
5329 image->set_ri(
true);
5341 float distance_sigma =
params[
"distance_sigma"];
5342 float value_sigma =
params[
"value_sigma"];
5343 int max_iter =
params[
"niter"];
5344 int half_width =
params[
"half_width"];
5346 if (half_width < distance_sigma) {
5347 LOGWARN(
"localwidth(=%d) should be larger than distance_sigma=(%f)\n",
5348 half_width, distance_sigma);
5351 distance_sigma *= distance_sigma;
5353 float image_sigma = image->get_attr(
"sigma");
5354 if (image_sigma > value_sigma) {
5355 LOGWARN(
"image sigma(=%f) should be smaller than value_sigma=(%f)\n",
5356 image_sigma, value_sigma);
5358 value_sigma *= value_sigma;
5360 int nx = image->get_xsize();
5361 int ny = image->get_ysize();
5362 int nz = image->get_zsize();
5365 int width=nx, height=ny;
5369 float tempfloat1,tempfloat2,tempfloat3;
5370 int index1,index2,index;
5372 int tempint1,tempint3;
5375 tempint3=width+2*half_width;
5377 float* mask=(
float*)calloc((2*half_width+1)*(2*half_width+1),
sizeof(
float));
5378 float* OrgImg=(
float*)calloc((2*half_width+width)*(2*half_width+height),
sizeof(