EMAN2
processor.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#include "processor.h"
33#include "sparx/processor_sparx.h"
35#include "reconstructor_tools.h"
36#include "cmp.h"
37#include "ctf.h"
38#include "xydata.h"
39#include "emdata.h"
40#include "emfft.h"
41#include "emassert.h"
42#include "randnum.h"
43#include "symmetry.h"
44#include "averager.h"
45#include "util.h"
46
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>
52#include <algorithm>
53#include <gsl/gsl_fit.h>
54#include <ctime>
55
56#ifdef __APPLE__
57 typedef unsigned int uint;
58#endif //__APPLE__
59
60#ifdef _WIN32
61 typedef unsigned int uint;
62#endif //_WIN32
63
64#ifdef EMAN2_USING_CUDA
65//#include "cuda/cuda_util.h"
66#include "cuda/cuda_processor.h"
67#endif // EMAN2_USING_CUDA
68
69using namespace EMAN;
70using std::reverse;
71
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";
126//const string BinarizeAmpProcessor::NAME = "threshold.amp.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"; // This is broken, it never worked. Somebody didn't think it through properly
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";
228// const string CoordinateMaskFileProcessor::NAME = "mask.fromfile.sizediff";
229const string PaintProcessor::NAME = "mask.paint";
230const string DirectionalSumProcessor::NAME = "misc.directional_sum";
231template<> const string BinaryOperateProcessor<MaxPixelOperator>::NAME = "math.max"; // These 2 should not really be processors
232template<> const string BinaryOperateProcessor<MinPixelOperator>::NAME = "math.min";
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";
322
323//#ifdef EMAN2_USING_CUDA
324//const string CudaMultProcessor::NAME = "cuda.math.mult";
325//const string CudaCorrelationProcessor::NAME = "cuda.correlate";
326//#endif //EMAN2_USING_CUDA
327
328#if 0
329//const string XYZProcessor::NAME = "XYZ";
330#endif //0
331
332
334{
335 force_add<HighpassAutoPeakProcessor>();
336 force_add<LinearRampFourierProcessor>();
337 force_add<LoGFourierProcessor>();
338 force_add<DoGFourierProcessor>();
339// force_add<AzSharpProcessor>();
340 force_add<FixSignProcessor>();
341
342 force_add<AmpweightFourierProcessor>();
343 force_add<Axis0FourierProcessor>();
344 force_add<Wiener2DFourierProcessor>();
345 force_add<LowpassAutoBProcessor>();
346 force_add<CtfSimProcessor>();
347
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>();
368
369 force_add<ClampingProcessor>();
370 force_add<NSigmaClampingProcessor>();
371
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>();
379// force_add<BinarizeAmpProcessor>();
380 force_add<BinarizeFourierProcessor>();
381 force_add<CollapseProcessor>();
382 force_add<LinearXformProcessor>();
383 force_add<SetBitsProcessor>();
384
385 force_add<CCCSNRProcessor>();
386 force_add<ExpProcessor>();
387 force_add<RangeThresholdProcessor>();
388 force_add<SigmaProcessor>();
389 force_add<LogProcessor>();
390 force_add<FiniteProcessor>();
391
392 force_add< BinaryOperateProcessor<MaxPixelOperator> >();
393 force_add< BinaryOperateProcessor<MinPixelOperator> >();
394
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>();
406
407 force_add<MaxShrinkProcessor>();
408 force_add<MinShrinkProcessor>();
409 force_add<MeanShrinkProcessor>();
410 force_add<MedianShrinkProcessor>();
411 force_add<FFTResampleProcessor>();
412 force_add<NonConvexProcessor>();
413
414 force_add<MakeRadiusSquaredProcessor>();
415 force_add<MakeRadiusProcessor>();
416
417 force_add<ComplexNormPixel>();
418
419 force_add<LaplacianProcessor>();
420// force_add<ZeroConstantProcessor>(); // this is badly written, it does not work and never did. Who wrote this !?!?
421
422 force_add<BoxMedianProcessor>();
423 force_add<BoxSigmaProcessor>();
424 force_add<BoxMaxProcessor>();
425 force_add<LocalMinAbsProcessor>();
426
427 force_add<MinusPeakProcessor>();
428 force_add<PeakOnlyProcessor>();
429 force_add<DiffBlockProcessor>();
430
431 force_add<CutoffBlockProcessor>();
432// force_add<GradientRemoverProcessor>();
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>();
443
444 force_add<BeamstopProcessor>();
445 force_add<MeanZeroEdgeProcessor>();
446 force_add<AverageXProcessor>();
447 force_add<DecayEdgeProcessor>();
448 force_add<ZeroEdgeRowProcessor>();
449 force_add<ZeroEdgePlaneProcessor>();
450
451 force_add<BilateralProcessor>();
452
453 force_add<ConvolutionProcessor>();
454 force_add<BispecSliceProcessor>();
455 force_add<HarmonicProcessor>();
456
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>();
469
470 force_add<HistogramBin>();
471
472 force_add<NormalizeToLeastSquareProcessor>();
473
474 force_add<RotationalAverageProcessor>();
475 force_add<RotationalSubstractProcessor>();
476 force_add<FlipProcessor>();
477 force_add<TransposeProcessor>();
478 force_add<MirrorProcessor>();
479 force_add<ReverseProcessor>();
480
481 force_add<AddNoiseProcessor>();
482 force_add<AddSigmaNoiseProcessor>();
483 force_add<AddRandomNoiseProcessor>();
484
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>();
497
498 force_add<CTFSNRWeightProcessor>();
499
500 force_add<ToMassCenterProcessor>();
501 force_add<ToCenterProcessor>();
502 force_add<PhaseToMassCenterProcessor>();
503 force_add<ACFCenterProcessor>();
504// force_add<SNRProcessor>();
505 force_add<CTFCorrProcessor>();
506 force_add<FSCFourierProcessor>();
507
508 force_add<XGradientProcessor>();
509 force_add<YGradientProcessor>();
510 force_add<ZGradientProcessor>();
511
512 force_add<ImageDivergenceProcessor>();
513 force_add<GradientMagnitudeProcessor>();
514 force_add<GradientDirectionProcessor>();
515 force_add<LaplacianMagnitudeProcessor>();
516 force_add<LaplacianDirectionProcessor>();
517 force_add<SDGDProcessor>();
518
519// force_add<FileFourierProcessor>();
520
521 force_add<SymSearchProcessor>();
522 force_add<MaskPackProcessor>();
523 force_add<StripeXYProcessor>();
524 force_add<BadLineXYProcessor>();
525 force_add<LocalNormProcessor>();
526
527 force_add<IndexMaskFileProcessor>();
528// force_add<CoordinateMaskFileProcessor>();
529 force_add<SetIsoPowProcessor>();
530 force_add<SetSFProcessor>();
531 force_add<MatchSFProcessor>();
532
533 force_add<SmartMaskProcessor>();
534
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>();
556
557 force_add<TomoTiltEdgeMaskProcessor>();
558 force_add<TomoTiltAngleWeightProcessor>();
559
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>();
584
585// force_add<WaveletProcessor>();
586 force_add<EnhanceProcessor>();
587 force_add<FFTProcessor>();
588 force_add<RadialProcessor>();
589
590 force_add<DirectionalSumProcessor>();
591 force_add<ConvolutionKernelProcessor>();
592
593 //Gorgon-related processors
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>();
605
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>();
620
621//#ifdef EMAN2_USING_CUDA
622// force_add<CudaMultProcessor>();
623// force_add<CudaCorrelationProcessor>();
624//#endif // EMAN2_USING_CUDA
625
626// force_add<XYZProcessor>();
627}
628
629void FiniteProcessor::process_pixel(float *x) const
630{
631 if ( !Util::goodf(x) ) {
632 *x = to;
633 }
634}
635
636
637EMData* Processor::process(const EMData * const image)
638{
639 EMData * result = image->copy();
640// printf("Default copy called\n");
641 process_inplace(result);
642 return result;
643}
644
646{
647 if (!image) {
648 LOGWARN("NULL image");
649 return;
650 }
651
652 EMData *processor_image = create_processor_image();
653
654 if (image->is_complex()) {
655 (*image) *= *processor_image;
656 }
657 else {
658 EMData *fft = image->do_fft();
659 (*fft) *= (*processor_image);
660 EMData *ift = fft->do_ift();
661
662 float *data = image->get_data();
663 float *t = data;
664 float *ift_data = ift->get_data();
665
666 data = ift_data;
667 ift_data = t;
668
669 ift->update();
670
671 if( fft )
672 {
673 delete fft;
674 fft = 0;
675 }
676
677 if( ift )
678 {
679 delete ift;
680 ift = 0;
681 }
682 }
683
684 image->update();
685}
686
687#define FFTRADIALOVERSAMPLE 4
689{
690 if (!image) {
691 LOGWARN("NULL Image");
692 return;
693 }
694
695 preprocess(image);
696
697 int array_size = FFTRADIALOVERSAMPLE * image->get_ysize();
698 float step=0.5f/array_size;
699
700 bool return_radial=(bool)params.set_default("return_radial",0);
701 vector < float >yarray(array_size);
702
703 create_radial_func(yarray);
704
705 if (image->is_complex()) {
706 image->apply_radial_func(0, step, yarray);
707 }
708 else {
709 EMData *fft = image->do_fft();
710 fft->apply_radial_func(0, step, yarray);
711 EMData *ift = fft->do_ift();
712
713 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));
714
715 //ift->update(); Unecessary
716
717 if( fft )
718 {
719 delete fft;
720 fft = 0;
721 }
722
723 if( ift )
724 {
725 delete ift;
726 ift = 0;
727 }
728 }
729 if (return_radial) image->set_attr("filter_curve",yarray);
730
731 image->update();
732}
733
735{
736 if (!image) {
737 LOGWARN("NULL Image");
738 return;
739 }
740
741 preprocess(image);
742
743// int array_size = FFTRADIALOVERSAMPLE * image->get_ysize();
744// float step=0.5f/array_size;
745//
746// vector < float >yarray(array_size);
747
748 bool return_radial=(bool)params.set_default("return_radial",0);
749 bool interpolate=(bool)params.set_default("interpolate",1);
750
751 float cornerscale;
752 if (image->get_zsize()>1) cornerscale=sqrt(3.0);
753 else cornerscale=sqrt(2.0);
754
755 if (image->is_complex()) {
756 vector <float>yarray = image->calc_radial_dist(floor(image->get_ysize()*cornerscale/2),0,1.0,1);
757 create_radial_func(yarray,image);
758 image->apply_radial_func(0, 1.0f/image->get_ysize(), yarray,interpolate);
759 if (return_radial) image->set_attr("filter_curve",yarray);
760 }
761 else {
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);
764 create_radial_func(yarray,fft);
765 // 4/30/10 stevel turned off interpolation to fix problem with matched filter
766 // 9/12/14 stevel, not sure why I turned off interp. Seems to cause rather than fix problems. Adding option to enable. Afraid to turn it on
767 fft->apply_radial_func(0, 1.0f/image->get_ysize(), yarray,interpolate);
768 EMData *ift = fft->do_ift();
769
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);
772
773// for (int i=0; i<yarray.size(); i++) printf("%d\t%f\n",i,yarray[i]);
774
775 //ift->update(); Unecessary
776
777 delete fft;
778 delete ift;
779
780 }
781
782 image->update();
783}
784
786 EMData *mask=(EMData *)params["mask"];
787 bool unpack = (int)params.set_default("unpack",0);
788
789 if ((float)mask->get_attr("mean_nonzero")!=1.0f) throw InvalidParameterException("MaskPackProcessor requires a binary mask");
790
791 int n_nz=(int)mask->get_attr("square_sum"); // true only for a binary mask
792
793 EMData *ret = 0;
794 if (unpack) {
795 ret = new EMData(mask->get_xsize(),mask->get_ysize(),mask->get_zsize());
796 ret->to_zero();
797
798 size_t ix=0;
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++));
801 }
802
803 } else {
804 ret = new EMData(n_nz,1,1);
805 ret->to_zero();
806
807 size_t ix=0;
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));
810 }
811 }
812
813 ret->update();
814 return ret;
815}
816
817// Looks like this hasn't actually been written yet...
819{
820 if (!image) {
821 LOGWARN("NULL Image");
822 return;
823 }
824
825 string mode=(string)params["mode"];
826
827 int nx = image->get_xsize();
828 int ny = image->get_ysize();
829 int nz = image->get_zsize();
830
831 // make an empty FFT object.
832 EMData *fft=new EMData((nx&0xfffffffe)+2,ny,nz);
833 fft->set_complex(1);
834 fft->set_ri(1);
835 if (nx&1) fft->set_fftpad(1);
836 fft->to_zero();
837
838 // Copy the kernel. Note that the kernel is 3x oversampled, so we are just computing an approximate kernel locally
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++) {
843 fft->set_complex_at(x,y,z,FourierInserter3DMode7::kernel[x*3][abs(y)*3][abs(z)*3]);
844 }
845 }
846 }
847 }
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++) {
852 fft->set_complex_at(x,y,z,FourierInserter3DMode11::kernel[x*3][abs(y)*3][abs(z)*3]);
853 }
854 }
855 }
856 }
857 else throw InvalidParameterException("Gridding kernel correction of unknown mode, only gridding_5 or gridding_7 allowed");
858
859 EMData *real=fft->do_ift(); // this is the kernel ift
860 real->process_inplace("xform.phaseorigin.tocenter");
861 real->mult(2.0f/(float)real->get_attr("maximum"));
862 real->process_inplace("math.reciprocal"); // reciprocal to make a correction volume
863 real->process_inplace("threshold.clampminmax",Dict("minval",-4.0f,"maxval",4.0f)); // block overcorrection of noise near the edges
864// real->write_image("invkernel.hdf");
865
866 image->mult(*real); // apply the correction
867 delete real;
868 delete fft;
869}
870
871// Looks like this hasn't actually been written yet...
873{
874 if (!image) {
875 LOGWARN("NULL Image");
876 return;
877 }
878
879 float az_scale=(float)params.set_default("az_scale",1.0);
880 float cornerscale;
881 EMData *fft;
882
883 if (image->is_complex()) fft=image;
884 else EMData *fft = image->do_fft();
885
886 int nx=fft->get_xsize();
887 int ny=fft->get_ysize();
888 int nz=fft->get_zsize();
889 if (nz!=1) {
890 }
891 else {
892
893 }
894 if (image->is_complex()) {
895 EMData *ift = fft->do_ift();
896
897 memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));
898
899 delete fft;
900 delete ift;
901 }
902
903 image->update();
904}
905
907{
908 EMData *fft;
909
910 if (!image) throw InvalidParameterException("FFTPeakProcessor: no image provided");
911 if (!image->is_complex()) fft = image->do_fft();
912 else fft = image;
913
914 int nx=fft->get_xsize();
915 int ny=fft->get_ysize();
916 int nz=fft->get_zsize();
917 float thresh_sigma = (float)params.set_default("thresh_sigma", 1.0);
918 bool removepeaks = (bool)params.set_default("removepeaks",0);
919 bool to_mean = (bool)params.set_default("to_mean",0);
920
921 vector<float> amp=fft->calc_radial_dist(nx/2,0,1,0); // amplitude
922 vector<float> thr=fft->calc_radial_dist(nx/2,0,1,1); // start with inten (avg of amp^2)
923
924 for (int i=0; i<nx/2; i++) thr[i]=thresh_sigma*sqrt(thr[i]-amp[i]*amp[i])+amp[i]; // sqrt(inten-amp^2) is st
925
926 if (nz>1) {
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++) {
930 float r2=Util::hypot3(x,y,z);
931 int r=int(r2);
932 if (r>=nx/2) continue;
933 complex<float> v=fft->get_complex_at(x,y,z);
934 float va=std::abs(v);
935 if (va>=thr[r]) {
936 if (removepeaks) {
937 if (to_mean) {
938 v*=amp[r]/va;
939 fft->set_complex_at(x,y,z,v);
940 }
941 else fft->set_complex_at(x,y,z,0.0);
942 }
943 }
944 else {
945 if (!removepeaks) {
946 if (to_mean) {
947 v*=amp[r]/va;
948 fft->set_complex_at(x,y,z,v);
949 }
950 else fft->set_complex_at(x,y,z,0.0);
951 }
952 }
953 }
954 }
955 }
956 }
957 else {
958 for (int y=-ny/2; y<ny/2; y++) {
959 for (int x=0; x<nx/2; x++) {
960 int r=Util::hypot_fast_int(x,y);
961 if (r>=nx/2) continue;
962 complex<float> v=fft->get_complex_at(x,y);
963 float va=std::abs(v);
964 if (va>=thr[r]) {
965 if (removepeaks) {
966 if (to_mean) {
967 v*=amp[r]/va;
968 fft->set_complex_at(x,y,v);
969 }
970 else fft->set_complex_at(x,y,0.0);
971 }
972 }
973 else {
974 if (!removepeaks) {
975 if (to_mean) {
976 v*=amp[r]/va;
977 fft->set_complex_at(x,y,v);
978 }
979 else fft->set_complex_at(x,y,0.0);
980 }
981 }
982 }
983 }
984 }
985
986 if (fft!=image) {
987 EMData *ift=fft->do_ift();
988 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));
989 delete fft;
990 delete ift;
991 }
992 image->update();
993
994// image->update();
995}
996
998{
999 EMData *fft;
1000
1001 if (!image) throw InvalidParameterException("FFTConeProcessor: no image provided");
1002 if (!image->is_complex()) fft = image->do_fft();
1003 else fft = image;
1004
1005
1006 int nx=fft->get_xsize();
1007 int ny=fft->get_ysize();
1008 int nz=fft->get_zsize();
1009 if (nz==1) throw ImageDimensionException("FFTConeProcessor only works on 3-D images");
1010
1011 float angle = (float)params.set_default("angle",15.0f);
1012 float rmin = (float)params.set_default("rmin",1.0f);
1013
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++) {
1017 float ang=0;
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;
1019 if (ang>angle || Util::hypot3(x,y,z)<rmin) continue;
1020 fft->set_complex_at(x,y,z,0.0f);
1021 }
1022 }
1023 }
1024
1025 if (fft!=image) {
1026 EMData *ift=fft->do_ift();
1027 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));
1028 delete fft;
1029 delete ift;
1030 }
1031 image->update();
1032
1033// image->update();
1034}
1035
1037{
1038 EMData *fft;
1039
1040 if (!image) throw InvalidParameterException("FFTWedgeProcessor: no image provided");
1041 if (!image->is_complex()) fft = image->do_fft();
1042 else fft = image;
1043
1044
1045 int nx=fft->get_xsize();
1046 int ny=fft->get_ysize();
1047 int nz=fft->get_zsize();
1048 if (nz==1) throw ImageDimensionException("FFTWedgeProcessor only works on 3-D images");
1049
1050 float anglemin = (float)params.set_default("anglemin",-30.0f);
1051 float anglemax = (float)params.set_default("anglemax",30.0f);
1052 float rmin = (float)params.set_default("rmin",1.0f);
1053
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++) {
1057 float ang=90.0f;
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);
1061 }
1062 }
1063 }
1064
1065 if (fft!=image) {
1066 EMData *ift=fft->do_ift();
1067 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));
1068 delete fft;
1069 delete ift;
1070 }
1071 image->update();
1072
1073// image->update();
1074}
1075
1076
1078{
1079 if (!image) throw InvalidParameterException("WedgeFillProcessor: no image provided");
1080 if (!image->is_complex()) throw ImageFormatException("WedgeFillProcessor: target image must be complex");
1081
1082 EMData *source=(EMData *)params.set_default("fillsource",(EMData *)NULL);
1083// if (!source) throw InvalidParameterException("WedgeFillProcessor: fillsource required");
1084
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");
1087
1088 int nx=image->get_xsize();
1089 int ny=image->get_ysize();
1090 int nz=image->get_zsize();
1091 float thresh_sigma = (float)params.set_default("thresh_sigma", 0.5);
1092 bool dosigma = 1 ? thresh_sigma>0.0 : 0;
1093
1094 float maxtilt = (float)params.set_default("maxtilt", 90.0);
1095 bool dotilt = 1 ? maxtilt <90.0 : 0;
1096 maxtilt*=M_PI/180.0;
1097
1098 vector<float> sigmaimg;
1099 if (dosigma) {
1100 sigmaimg=image->calc_radial_dist(nx/2,0,1,4);
1101 for (int i=0; i<nx/2; i++) sigmaimg[i]*=sigmaimg[i]*thresh_sigma; // anything less than 1/10 sigma is considered to be missing
1102 }
1103
1104 vector<int> realpixel(nx/2);
1105 for (int i=0; i<nx/2; i++) realpixel[i]=0;
1106
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) {
1110 float r2=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z); // origin at 0,0; periodic
1111 int r=int(r2);
1112 if (r<3) continue; // too few points at r<3 to consider any "missing"
1113
1114 float tilt = 0.0;
1115 if (dotilt) tilt=atan2((float)(z<nz/2?z:nz-z),(float)(x/2));
1116
1117 float v1r=image->get_value_at(x,y,z);
1118 float v1i=image->get_value_at(x+1,y,z);
1119 float v1=Util::square_sum(v1r,v1i);
1120// if (r<10) printf("%d %d %d %d\t%1.3g %1.3g\n",x,y,z,r,v1,sigmaimg[r]);
1121 if ((!dosigma || v1>sigmaimg[r]) && r<nx/2 && (!dotilt || fabs(tilt)<maxtilt)){
1122 realpixel[r]++;
1123 continue;
1124 }
1125
1126 if (!source) {
1127 image->set_value_at_fast(x,y,z,0);
1128 image->set_value_at_fast(x+1,y,z,0);
1129 }
1130 else {
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));
1133 }
1134 }
1135 }
1136 }
1137
1138 image->set_attr("real_pixels", realpixel);
1139 image->update();
1140}
1141
1142
1143void LowpassAutoBProcessor::create_radial_func(vector < float >&radial_mask,EMData *image) const{
1144 float apix=(float)image->get_attr("apix_x");
1145 int verbose=(int)params["verbose"];
1146// int adaptnoise=params.set_default("adaptnoise",0);
1147 float noisecutoff=(float)params.set_default("noisecutoff",1.0/6.0);
1148 if (apix<=0 || apix>7.0f) throw ImageFormatException("apix_x > 7.0 or <0");
1149 float ds=1.0f/(apix*image->get_ysize()); // 0.5 is because radial mask is 2x oversampled
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;
1156 }
1157 if (end<start+2) {
1158 printf("WARNING: noise cutoff too close to 15 A ! Results will not be good...");
1159 start=end-5;
1160 }
1161
1162 FILE *out=NULL;
1163 if (verbose>2) {
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");
1166 }
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++ ) { // -2 is arbitrary because sometimes the last pixel or two have funny values
1172 x[i-start]=ds*ds*i*i;
1173 if (radial_mask[i]>0) y[i-start]=log(radial_mask[i]); // ok
1174 else if (i>start) y[i-start]=y[i-start-1]; // not good
1175 else y[i-start]=0.0; // bad
1176 if (i>start &&i<radial_mask.size()-3) dy[i-start]=y[i-start]-y[i-start-1]; // creates a 'derivative' of sorts, for use in adaptnoise
1177 if (out) fprintf(out,"%f\t%f\n",x[i-start],y[i-start]);
1178 }
1179 if (out) fclose(out);
1180
1181 float slope=0,intercept=0;
1182 Util::calc_least_square_fit(end-start, x,y,&slope,&intercept,1);
1183
1184 if (verbose) printf("slope=%f intercept=%f\n",slope,intercept);
1185
1186 float B=(float)params["bfactor"]+slope*4.0f/2.0f; // *4 is for Henderson definition, 2.0 is for intensity vs amplitude
1187 float B2=(float)params["bfactor"];
1188
1189 if (verbose) printf("User B = %1.2f Corrective B = %1.2f Total B=%1.3f\n",(float)params["bfactor"],slope*2.0,B);
1190
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);
1195 }
1196 if (verbose>1) Util::save_data(0,ds,radial_mask,"filter.txt");
1197
1198 free(x);
1199 free(y);
1200 free(dy);
1201 }
1202
1204 size_t nxyz = image->get_size();
1205 int mode=params.set_default("wiener",0);
1206 float scale=params.set_default("scalesnr",2.0f);
1207
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);
1214 }
1215}
1216
1217
1219{
1220 int ys=image->get_ysize();
1221
1222 EMData *mask1=new EMData(ys,ys,1);
1223 mask1->process_inplace("mask.gaussian",Dict("outer_radius", ys/2.0));
1224 EMData *mask2=mask1->copy();
1225 mask2->mult(-1.0f);
1226 mask2->add(1.0);
1227 mask2->process_inplace("mask.decayedge2d",Dict("width",4));
1228
1229/*
1230
1231
1232
1233mask1=EMData(ys2,ys2,1)
1234 mask1.to_one()
1235 mask1.process_inplace("mask.gaussian",{"outer_radius":radius})
1236 mask2=mask1.copy()*-1+1
1237# mask1.process_inplace("mask.decayedge2d",{"width":4})
1238 mask2.process_inplace("mask.decayedge2d",{"width":4})
1239 mask1.clip_inplace(Region(-(ys2*(oversamp-1)/2),-(ys2*(oversamp-1)/2),ys,ys))
1240 mask2.clip_inplace(Region(-(ys2*(oversamp-1)/2),-(ys2*(oversamp-1)/2),ys,ys))
1241 ratio1=mask1.get_attr("square_sum")/(ys*ys) #/1.035
1242 ratio2=mask2.get_attr("square_sum")/(ys*ys)
1243 masks[(ys,radius)]=(mask1,ratio1,mask2,ratio2)
1244
1245
1246
1247*/
1248}
1249
1251{
1252 EMData *fft;
1253 float *fftd;
1254 int f=0;
1255// static float sum1=0,sum1a=0;
1256// static double sum2=0,sum2a=0;
1257
1258 if (!image) {
1259 LOGWARN("NULL Image");
1260 return;
1261 }
1262
1263 if (!image->is_complex()) fft = image->do_fft();
1264 else fft = image;
1265
1266 int neighbor=params.set_default("neighbor",0);
1267 float neighbornorm=params.set_default("neighbornorm",sqrt(2.0f));
1268
1269 int nx=fft->get_xsize()/2;
1270 int ny=fft->get_ysize();
1271 if (params.set_default("x",1)) {
1272 if (neighbor) {
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); // sqrt is because the pixels are normally pretty noisy
1274 }
1275 else {
1276 for (int x=1; x<nx; x++) fft->set_complex_at(x,0,0.0f);
1277 }
1278 }
1279 if (params.set_default("y",1)) {
1280 if (neighbor) {
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); // sqrt is because the pixels are normally pretty noisy
1282 }
1283 else {
1284 for (int y=1; y<ny; y++) fft->set_complex_at(0,y,0.0f);
1285 }
1286 }
1287
1288 if (fft!=image) {
1289// fft->update();
1290 EMData *ift=fft->do_ift();
1291 memcpy(image->get_data(),ift->get_data(),(nx*2-2)*ny*sizeof(float));
1292 delete fft;
1293 delete ift;
1294 }
1295
1296// image->update();
1297
1298}
1299
1301{
1302 if(params.has_key("apix")) {
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"]);
1306 }
1307
1308 int xynoz = params.set_default("xynoz",0);
1309
1310 const Dict dict = image->get_attr_dict();
1311 if( params.has_key("cutoff_freq") ) {
1312 float val = (float)params["cutoff_freq"] * (float)dict["apix_x"];
1313 params["cutoff_abs"] = val;
1314 }
1315 else if( params.has_key("cutoff_pixels") ) {
1316 float val = ((float)params["cutoff_pixels"] / (float)dict["nx"]);
1317 params["cutoff_abs"] = val;
1318 }
1319
1320 float omega = params["cutoff_abs"];
1321 float zcenter=params.set_default("centerfreq",0.0f);
1322 float hppix = params.set_default("hppix",-1.0f);
1323
1324 omega = (omega<0?-1.0:1.0)*0.5f/omega/omega;
1325
1326 EMData *fft;
1327 int f=0;
1328
1329 if (!image) {
1330 LOGWARN("NULL Image");
1331 return;
1332 }
1333
1334 if (!image->is_complex()) {
1335 fft = image->do_fft();
1336 f=1;
1337 }
1338 else {
1339 fft=image;
1340 }
1341
1342 int nx=fft->get_xsize();
1343 int ny=fft->get_ysize();
1344 int nz=fft->get_ysize();
1345 omega /=(nz*nz)/4;
1346 zcenter=zcenter*(float)dict["apix_x"]*nz;
1347
1348 if (xynoz) {
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);
1353 float r=Util::hypot_fast(x,y);
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));
1355 }
1356 }
1357 }
1358 }
1359 else {
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));
1365 }
1366 }
1367 }
1368 }
1369
1370 if (f) {
1371 EMData *ift=fft->do_ift();
1372 memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));
1373 delete fft;
1374 delete ift;
1375 }
1376
1377}
1378
1379
1381{
1382 EMData *fft;
1383 float *fftd;
1384 int f=0;
1385// static float sum1=0,sum1a=0;
1386// static double sum2=0,sum2a=0;
1387
1388 if (!image) {
1389 LOGWARN("NULL Image");
1390 return;
1391 }
1392
1393 if (!image->is_complex()) {
1394 fft = image->do_fft();
1395 fftd = fft->get_data();
1396 f=1;
1397 }
1398 else {
1399 fft=image;
1400 fftd=image->get_data();
1401 }
1402 float *sumd = NULL;
1403 if (sum) sumd=sum->get_data();
1404//printf("%d %d %d %d\n",fft->get_xsize(),fft->get_ysize(),sum->get_xsize(),sum->get_ysize());
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) {
1407 float c;
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; // prevents divide by zero in normalization
1411 fftd[i]*=c;
1412 fftd[i+1]*=c;
1413 if (sumd) { sumd[i]+=c; sumd[i+1]+=0; }
1414
1415 // debugging
1416/* if (i==290*1+12) {
1417 sum1+=fftd[i];
1418 sum2+=fftd[i];
1419 printf("%f\t%f\t%f\t%f\t%f\t%f\n",sum1,sum2,fftd[i],fftd[i+1],sumd[i],c);
1420 }
1421 if (i==290*50+60) {
1422 sum1a+=fftd[i];
1423 sum2a+=fftd[i];
1424 printf("%f\t%f\t%f\t%f\t%f\t%f\n",sum1a,sum2a,fftd[i],fftd[i+1],sumd[i],c);
1425 }*/
1426 }
1427
1428 if (f) {
1429 fft->update();
1430 EMData *ift=fft->do_ift();
1431 memcpy(image->get_data(),ift->get_data(),n*sizeof(float));
1432 delete fft;
1433 delete ift;
1434 }
1435
1436 sum->update();
1437 image->update();
1438
1439}
1440
1442{
1443 printf("Process inplace not implemented. Please use process.\n");
1444 return;
1445}
1446
1447
1449{
1450
1451 float minratio = params.set_default("minratio",0.5f);
1452 int maxnseg = params.set_default("maxnseg",0);
1453 int skipseg = params.set_default("skipseg",0);
1454 float width = params.set_default("width",10.0f);
1455 int verbose = params.set_default("verbose",0);
1456 EMData *mask= (EMData *)params.set_default("mask",(EMData *)NULL);
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();
1461
1462 EMData * result = image->process("threshold.belowtozero",Dict("minval",0.0f));
1463 // The intent of these filters is to insure that the map effectively
1464// result->process_inplace("filter.lowpass.gauss",Dict("cutoff_freq",0.45/width)); // 0.45 = sqrt(2)/pi, resolvability threshold
1465// result->process_inplace("filter.lowpass.gauss",Dict("cutoff_freq",0.637/width)); // 0.637 = 2/pi, (may be Rayleigh criterion?)
1466 if (mask!=NULL) {
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));
1469 }
1470
1471 XYData gsf;
1472 gsf.make_gauss(nx*4,1.0f/apix,0.45/width);
1473// gsf.make_gauss(nx*4,1.0f/apix,0.637/width);
1474 result->process_inplace("filter.setstrucfac",Dict("apix",apix,"strucfac",&gsf));
1475
1476 // Generate a Gaussian we can subtract out of the volume quickly
1477 int gs=2*width/(float)image->get_attr("apix_x");
1478 EMData gsub(gs,gs,gs);
1479 gsub.to_one();
1480 gsub.process_inplace("mask.gaussian",Dict("outer_radius",int(width/(2.0*apix))));
1481
1482 if (verbose>2) {
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"));
1486 result->write_image("segfilt.hdf",0,EMUtil::IMAGE_UNKNOWN,false,0,EMUtil::EM_COMPRESSED);
1487 }
1488
1489 EMData *cache = NULL;
1490 if (skipseg==2) cache=result->copy();
1491
1492 // original version of the code. A bit slow but works very well
1493// vector<float> centers;
1494// vector<float> amps;
1495// if (maxnseg==0) maxnseg=2000000000;
1496// float maxval=0.0f;
1497//
1498// while (amps.size()<maxnseg) {
1499// IntPoint p = result->calc_max_location();
1500// FloatPoint fp(p[0],p[1],p[2]);
1501//
1502// float amp=result->get_value_at(p[0],p[1],p[2]);
1503// if (amp<maxval*minratio) break;
1504// amps.push_back(amp);
1505// centers.push_back(p[0]);
1506// centers.push_back(p[1]);
1507// centers.push_back(p[2]);
1508// if (amp>maxval) maxval=amp; // really this should only happen once...
1509//
1510// result->insert_scaled_sum(&gsub,fp,1.0,-amp);
1511// }
1512
1513 // This was an alternate version, to improve speed. While it is faster, the approximations its making
1514 // produce slightly less optimal results.
1515 vector<float> centers;
1516 vector<float> amps;
1517 if (maxnseg==0) maxnseg=2000000000;
1518
1519 float thr=(float)result->get_attr("maximum")*minratio;
1520 while (amps.size()<maxnseg) {
1521 if ((float)result->get_attr("maximum")<=thr) break;
1522 // We find all peak values greater than our threshold
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());
1526
1527 for (int i=0; i<pixels.size(); i++) {
1528 IntPoint p = pixels[i].get_point();
1529 FloatPoint fp(p[0],p[1],p[2]);
1530
1531 float amp=result->get_value_at(p[0],p[1],p[2]);
1532 if (amp!=pixels[i].get_value()) continue; // skip any values near a value we've just changed, should come back in the next round
1533 // printf("%d %d %d %f %f %f\n",p[0],p[1],p[2],amp,pixels[i].get_value(),pixels[i+1].get_value());
1534// if (i<pixels.size()-1 && amp<pixels[i+1].get_value()) continue; // if one of the identified peaks has been reduced by a nearby peak too much
1535 if (amp<thr) continue;
1536// if (amp<thr) { maxnseg=amps.size(); break; }
1537 amps.push_back(amp);
1538 centers.push_back(p[0]);
1539 centers.push_back(p[1]);
1540 centers.push_back(p[2]);
1541
1542 result->insert_scaled_sum(&gsub,fp,1.0,-amp);
1543 }
1544 }
1545
1546 if (verbose) printf("Found %d centers\n",amps.size());
1547
1548 if (verbose>2) {
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"));
1552 result->write_image("segfilt.hdf",1,EMUtil::IMAGE_UNKNOWN,false,0,EMUtil::EM_COMPRESSED);
1553 }
1554 if (!skipseg) {
1555 // after we have our list of centers classify pixels
1556 for (int z=0; z<nz; z++) {
1557 for (int y=0; y<ny; y++) {
1558 for (int x=0; x<nz; x++) {
1559// if (image->get_value_at(x,y,z)<thr) {
1560// result->set_value_at(x,y,z,-1.0); //below threshold -> -1 (unclassified)
1561// continue;
1562// }
1563 int bcls=-1; // best matching class
1564 float bdist=(float)(nx+ny+nz); // distance for best class
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; }
1568 }
1569 result->set_value_at(x,y,z,(float)bcls); // set the pixel to the class number
1570 }
1571 }
1572 }
1573 if (verbose) printf("segmented\n");
1574 }
1575
1576 //if skipseg is set to 2, we return the filtered, unsegmented map (with valid centers)
1577 if (skipseg==2) {
1578 delete result;
1579 result=cache;
1580 }
1581
1582 result->set_attr("segment_centers",centers);
1583 result->set_attr("segment_amps",amps);
1584
1585 return result;
1586}
1587
1589{
1590 printf("Process inplace not implemented. Please use process.\n");
1591 return;
1592}
1593
1594
1596{
1597 EMData * result = image->copy();
1598
1599 float thr = params.set_default("thr",0.9f);
1600 float minsegsep = params.set_default("minsegsep",5.0f);
1601 float maxsegsep = params.set_default("maxsegsep",5.1f);
1602 int verbose = params.set_default("verbose",0);
1603
1604 vector<Pixel> pixels=image->calc_highest_locations(thr);
1605
1606 vector<float> centers(3); // only 1 to start
1607 int nx=image->get_xsize();
1608 int ny=image->get_ysize();
1609 int nz=image->get_zsize();
1610// int nxy=nx*ny;
1611
1612 // seed the process with the highest valued point
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());
1617
1618 // outer loop. We add one center per iteration
1619 // This is NOT a very efficient algorithm, it assumes points are fairly sparse
1620 while (pixels.size()>0) {
1621 // iterate over pixels until we find a new center (then stop), delete any 'bad' pixels
1622 // no iterators because we remove elements
1623
1624 for (unsigned int i=0; i<pixels.size(); i++) {
1625
1626 Pixel p=pixels[i];
1627 // iterate over existing centers to see if this pixel should be removed ... technically we only should need to check the last center
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);
1630 if (d<minsegsep) { // conflicts with existing center, erase
1631 pixels.erase(pixels.begin()+i);
1632 i--;
1633 break;
1634 }
1635 }
1636 }
1637
1638 int found=0;
1639 for (unsigned int i=0; i<pixels.size() && found==0; i++) {
1640 Pixel p=pixels[i];
1641
1642 // iterate again to see if this may be a new valid center. Start at the end so we tend to build chains
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);
1645 if (d<maxsegsep) { // we passed minsegsep question already, so we know we're in the 'good' range
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); // in the centers list now, don't need it any more
1650 found=1;
1651 break;
1652 }
1653 }
1654 }
1655
1656 // If we went through the whole list and didn't find one, we need to reseed again
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());
1663 }
1664
1665 if (verbose) printf("%d points found\n",(int)(centers.size()/3));
1666 }
1667
1668 // after we have our list of centers classify pixels
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); //below threshold -> -1 (unclassified)
1674 continue;
1675 }
1676 int bcls=-1; // best matching class
1677 float bdist=(float)(nx+ny+nz); // distance for best class
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; }
1681 }
1682 result->set_value_at(x,y,z,(float)bcls); // set the pixel to the class number
1683 }
1684 }
1685 }
1686
1687 result->set_attr("segment_centers",centers);
1688
1689 return result;
1690}
1691
1693{
1694 string s=(string)params.set_default("sym","c1");
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();
1698
1699 Averager* imgavg = Factory<Averager>::get((string)params.set_default("averager","mean"));
1700
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");
1703 if (n<=0) throw InvalidValueException(n,"xform.applysym: Cn symmetry, n>0");
1704
1705 for (int i=0; i<n; i++) {
1706 Transform t(Dict("type","2d","alpha",(float)(i*360.0f/n)));
1707 EMData* transformed = image->process("xform",Dict("transform",&t));
1708 imgavg->add_image(transformed);
1709 delete transformed;
1710 }
1711 EMData *ret=imgavg->finish();
1712 delete imgavg;
1713 return ret;
1714 }
1715
1716 Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));
1717 vector<Transform> transforms = sym->get_syms();
1718
1719 for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {
1720 Transform t = *trans_it;
1721 EMData* transformed = image->process("xform",Dict("transform",&t));
1722 imgavg->add_image(transformed);
1723 delete transformed;
1724 }
1725 EMData *ret=imgavg->finish();
1726 delete imgavg;
1727 return ret;
1728}
1729
1731{
1732 EMData *tmp=process(image);
1733 memcpy(image->get_data(),tmp->get_data(),(size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
1734 delete tmp;
1735 image->update();
1736 return;
1737}
1738
1740{
1741 EMData * result = image->copy();
1742
1743 int nseg = params.set_default("nseg",12);
1744 float thr = params.set_default("thr",-1.0e30f);
1745 int ampweight = params.set_default("ampweight",1);
1746 float maxsegsize = params.set_default("maxsegsize",10000.0f);
1747 float minsegsep = params.set_default("minsegsep",0.0f);
1748 int maxiter = params.set_default("maxiter",100);
1749 int maxvoxmove = params.set_default("maxvoxmove",25);
1750 int verbose = params.set_default("verbose",0);
1751 bool pseudoatom = params.set_default("pseudoatom",0);
1752 float sep = params.set_default("sep",3.78f);
1753
1754 int nx=image->get_xsize();
1755 int ny=image->get_ysize();
1756 int nz=image->get_zsize();
1757// int nxy=nx*ny;
1758// image->process_inplace("threshold.belowtozero");
1759 if (thr==-1.0e30f){
1760 thr=float(image->get_attr("mean_nonzero"))+ 1.0 *float(image->get_attr("sigma_nonzero"));
1761 printf("Estimated map threshold: %4f\n", thr);
1762 }
1763 // seed
1764 vector<float> centers(nseg*3);
1765 vector<float> count(nseg);
1766 // Alternative seeding method for paudoatom generation. Seed on the gird.
1767 if (pseudoatom){
1768 float ax=image->get_attr("apix_x");
1769 sep/=ax;
1770 if (verbose) printf("Seeding .....\n");
1771 int sx=int(nx/sep)+1,sy=int(ny/sep)+1,sz=int(nz/sep)+1;
1772 EMData m(sx,sy,sz);
1773 EMData mcount(sx,sy,sz);
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);
1779 if (v>thr){
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));
1782 }
1783 }
1784 }
1785 }
1786 m.div((nx/sx)*(ny/sy)*(nz/sz));
1787 int nsum=0;
1788 float l=image->get_attr("minimum"),r=image->get_attr("maximum"),th=0;
1789 while (abs(nsum-nseg)>0){
1790 th=(l+r)/2;
1791 nsum=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;
1796 }
1797 }
1798 }
1799// if (verbose) printf("%3f\t %3f\t %3f,\t %4d\t %4d\n", l,th,r,nsum,nseg);
1800 if (nsum>nseg) l=th;
1801 if (nsum<nseg) r=th;
1802 if ((r-l)<.0001)break;
1803 }
1804// nseg=nsum;
1805 if (verbose) printf("%3d pseudoatoms seeded at threshold value %3f\n", nsum, th);
1806 int q=0;
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){
1811 if(q<nseg*3){
1812 centers[q]= float(i+.5)*sep;
1813 centers[q+1]=float(j+.5)*sep;
1814 centers[q+2]=float(k+.5)*sep;
1815 q+=3;
1816 }
1817 }
1818 }
1819 }
1820 }
1821
1822 }
1823 // Default: random seeding.
1824 else{
1825 for (int i=0; i<nseg*3; i+=3) {
1826 centers[i]= Util::get_frand(0.0f,(float)nx);
1827 centers[i+1]=Util::get_frand(0.0f,(float)ny);
1828 centers[i+2]=Util::get_frand(0.0f,(float)nz);
1829 }
1830 }
1831
1832 for (int iter=0; iter<maxiter; iter++) {
1833 // **** classify
1834 size_t pixmov=0; // count of moved pixels
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); //below threshold -> -1 (unclassified)
1840 continue;
1841 }
1842 int bcls=-1; // best matching class
1843 float bdist=(float)(nx+ny+nz); // distance for best class
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; }
1847 }
1848 if ((int)result->get_value_at(x,y,z)!=bcls) pixmov++;
1849 if (bdist>maxsegsize) result->set_value_at(x,y,z,-1); // pixel is too far from any center
1850 else result->set_value_at(x,y,z,(float)bcls); // set the pixel to the class number
1851 }
1852 }
1853 }
1854
1855 // **** adjust centers
1856 for (int i=0; i<nseg*3; i++) centers[i]=0;
1857 for (int i=0; i<nseg; i++) count[i]=0;
1858 // weighted sums
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;
1864 float w=1.0;
1865 if (ampweight) w=image->get_value_at(x,y,z);
1866
1867 centers[cls*3]+=x*w;
1868 centers[cls*3+1]+=y*w;
1869 centers[cls*3+2]+=z*w;
1870 count[cls]+=w;
1871 }
1872 }
1873 }
1874 // now each becomes center of mass, or gets randomly reseeded
1875 int nreseed=0;
1876 for (int c=0; c<nseg; c++) {
1877 // reseed
1878 if (count[c]==0) {
1879 nreseed++;
1880 do {
1881 centers[c*3]= Util::get_frand(0.0f,(float)nx);
1882 centers[c*3+1]=Util::get_frand(0.0f,(float)ny);
1883 centers[c*3+2]=Util::get_frand(0.0f,(float)nz);
1884 } while (image->get_value_at((int)centers[c*3],(int)centers[c*3+1],(int)centers[c*3+2])<thr); // This makes sure the new point is inside density
1885 }
1886 // center of mass
1887 else {
1888 centers[c*3]/=count[c];
1889 centers[c*3+1]/=count[c];
1890 centers[c*3+2]/=count[c];
1891 }
1892 }
1893
1894 // with minsegsep, check separation
1895 if (minsegsep>0) {
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) {
1899 nreseed++;
1900 do {
1901 centers[c1*3]= Util::get_frand(0.0f,(float)nx);
1902 centers[c1*3+1]=Util::get_frand(0.0f,(float)ny);
1903 centers[c1*3+2]=Util::get_frand(0.0f,(float)nz);
1904 } while (image->get_value_at((int)centers[c1*3],(int)centers[c1*3+1],(int)centers[c1*3+2])<thr);
1905 }
1906 }
1907 }
1908 }
1909
1910
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; // termination conditions met
1913 }
1914
1915 result->set_attr("segment_centers",centers);
1916
1917 return result;
1918}
1919
1921{
1922 printf("Process inplace not implemented. Please use process.\n");
1923 return;
1924}
1925
1926
1928
1929// if (image->get_zsize()!=1) { throw ImageDimensionException("Only 2-D images supported"); }
1930
1931 int nx=image->get_xsize();
1932 int ny=image->get_ysize();
1933 int nz=image->get_zsize();
1934 float x0 = params.set_default("x0",nx/2);
1935 float y0 = params.set_default("y0",ny/2);
1936 float z0 = params.set_default("z0",nz/2);
1937 float xwidth = params.set_default("xwidth",nx);
1938 float ywidth = params.set_default("ywidth",ny);
1939 float zwidth = params.set_default("zwidth",nz);
1940
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;
1947 float val=0.0;
1948 if (xf>0&&yf>0&&zf>0) val=xf*yf*zf;
1949 image->mult_value_at_fast(x,y,z,val);
1950 }
1951 }
1952 }
1953}
1954
1956{
1957// TODO NOT IMPLEMENTED YET !!!
1958 EMData *ret = 0;
1959 const EMData *fft;
1960// float *fftd;
1961// int f=0;
1962
1963 if (!image) {
1964 LOGWARN("NULL Image");
1965 return ret;
1966 }
1967 throw NullPointerException("Processor not yet implemented");
1968
1969// if (!image->is_complex()) {
1970// fft = image->do_fft();
1971// fftd = fft->get_data();
1972// f=1;
1973// }
1974// else {
1975// fft=image;
1976// fftd=image->get_data();
1977// }
1978
1979 return ret;
1980}
1981
1983 EMData *tmp=process(image);
1984 memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
1985 delete tmp;
1986 image->update();
1987 return;
1988}
1989
1990
1992{
1993 const EMData *in2 = 0;
1994 if (in->is_complex()) in2=in;
1995 else in=in->do_fft();
1996
1997 EMData *filt = in->copy_head();
1998 Ctf *ictf = ctf;
1999
2000 if (!ictf) ctf=(Ctf *)in->get_attr("ctf");
2001
2003 filt->mult(*in2);
2004 EMData *ret=filt->do_ift();
2005
2006 delete filt;
2007 if (!in->is_complex()) delete in2;
2008
2009 if(!ictf && ctf) {delete ctf; ctf=0;}
2010 return(ret);
2011/* const EMData *fft;
2012 float *fftd;
2013 int f=0;
2014
2015 if (!image) {
2016 LOGWARN("NULL Image");
2017 return ret;
2018 }
2019
2020 if (!image->is_complex()) {
2021 fft = image->do_fft();
2022 fftd = fft->get_data();
2023 f=1;
2024 }
2025 else {
2026 fft=image;
2027 fftd=image->get_data();
2028 }
2029 powd=image->get_data();
2030
2031 int bad=0;
2032 for (int i=0; i<image->get_xsize()*image->get_ysize(); i+=2) {
2033 snr=(fftd[i]*fftd[i]+fftd[i+1]*fftd[i+1]-powd[i])/powd[i];
2034 if (snr<0) { bad++; snr=0; }
2035
2036 }
2037
2038 print("%d bad pixels\n",snr);
2039*/ return ret;
2040
2041}
2042
2044 EMData *tmp=process(image);
2045 memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
2046 delete tmp;
2047 image->update();
2048 return;
2049}
2050
2051void LinearRampFourierProcessor::create_radial_func(vector < float >&radial_mask) const
2052{
2053 Assert(radial_mask.size() > 0);
2054 for (size_t i = 0; i < radial_mask.size(); i++) {
2055 radial_mask[i] = (float)i;
2056 }
2057}
2058
2059void LowpassRandomPhaseProcessor::create_radial_func(vector < float >&radial_mask) const { };
2060
2062{
2063 float cutoff=0;
2064 preprocess(image);
2065 if( params.has_key("cutoff_abs") ) {
2066 cutoff=(float)params["cutoff_abs"];
2067 }
2068 else {
2069 printf("A cutoff_* parameter is required by filter.lowpass.randomphase\n");
2070 return;
2071 }
2072
2073
2074 if (image->get_zsize()==1) {
2075 int flag=0;
2076 if (!image->is_complex()) { image->do_fft_inplace(); flag=1; }
2077 image->ri2ap();
2078 int nx=image->get_xsize();
2079 int ny=image->get_ysize();
2080
2081 int z=0;
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); // location of the amplitude
2087 data[idx+1]=Util::get_frand(0.0f,(float)(M_PI*2.0));
2088 }
2089 }
2090 }
2091
2092 image->ap2ri();
2093
2094 if (flag) {
2095 image->do_ift_inplace();
2096 image->depad();
2097 }
2098 }
2099 else { // 3D
2100 int flag=0;
2101 if (!image->is_complex()) { image->do_fft_inplace(); flag=1; }
2102 image->ri2ap();
2103 int nx=image->get_xsize();
2104 int ny=image->get_ysize();
2105 int nz=image->get_zsize();
2106
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++) {
2111 if (Util::hypot3(x/float(nx),y/float(ny),z/float(nz))>=cutoff) {
2112 size_t idx=image->get_complex_index_fast(x,y,z); // location of the amplitude
2113 data[idx+1]=Util::get_frand(0.0f,(float)(M_PI*2.0));
2114 }
2115 }
2116 }
2117 }
2118 image->ap2ri();
2119
2120 if (flag) {
2121 image->do_ift_inplace();
2122 image->depad();
2123 }
2124 }
2125}
2126
2128{
2129 if(params.has_key("apix")) {
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"]);
2133 }
2134
2135 const Dict dict = image->get_attr_dict();
2136
2137 if( params.has_key("cutoff_abs") ) {
2138 highpass = params["cutoff_abs"];
2139 }
2140 else if( params.has_key("cutoff_freq") ) {
2141 highpass = (float)params["cutoff_freq"] * (float)dict["apix_x"] * (float)dict["ny"] / 2.0f;
2142 }
2143 else if( params.has_key("cutoff_pixels") ) {
2144 highpass = (float)params["cutoff_pixels"] / (float)dict["nx"];
2145 }
2146}
2147
2148void HighpassAutoPeakProcessor::create_radial_func(vector < float >&radial_mask, EMData *) const
2149{
2150 unsigned int c;
2151
2152// for (unsigned int i=0; i<radial_mask.size(); i++) printf("%d\t%f\n",i,radial_mask[i]);
2153 for (c=2; c<radial_mask.size(); c++) if (radial_mask[c-1]<=radial_mask[c]) break;
2154 if (c>highpass) c=(unsigned int)highpass; // the *2 is for the 2x oversampling
2155
2156 radial_mask[0]=0.0;
2157// for (int i=1; i<radial_mask.size(); i++) radial_mask[i]=(i<=c?radial_mask[c+1]/radial_mask[i]:1.0);
2158 for (unsigned int i=1; i<radial_mask.size(); i++) radial_mask[i]=(i<=c?0.0f:1.0f);
2159
2160 printf("%f %d\n",highpass,c);
2161// for (unsigned int i=0; i<radial_mask.size(); i++) printf("%d\t%f\n",i,radial_mask[i]);
2162
2163}
2164
2165void LinearRampProcessor::create_radial_func(vector < float >&radial_mask) const
2166{
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++) {
2171 radial_mask[i] = intercept + ((slope - intercept) * i) / (size - 1.0f);
2172 x += step;
2173 }
2174}
2175
2176void LoGFourierProcessor::create_radial_func(vector < float >&radial_mask) const
2177{
2178
2179 Assert(radial_mask.size() > 0);
2180 float x = 0.0f , nqstep = 0.5f/radial_mask.size();
2181 size_t size=radial_mask.size();
2182 float var = sigma*sigma;
2183 for (size_t i = 0; i < size; i++) {
2184 radial_mask[i] = ((x*x - var)/var*var)*exp(-x*x/2*var);
2185 x += nqstep;
2186 }
2187}
2188
2189void DoGFourierProcessor::create_radial_func(vector < float >&radial_mask) const
2190{
2191
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++) {
2197 radial_mask[i] = norm*((1.0f/sigma1*exp(-x*x/(2.0f*sigma1*sigma1))) - (1.0f/sigma2*exp(-x*x/(2.0f*sigma2*sigma2))));
2198 x += nqstep;
2199 }
2200}
2201
2203{
2204 if (!image) {
2205 LOGWARN("NULL Image");
2206 return;
2207 }
2208
2209 float gauss_width = params.set_default("gauss_width",0.0f);
2210
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();
2216
2217 for (size_t i = 0; i < size; ++i) {
2218 if (data[i]>=minval && data[i]<=maxval) data[i]=0.0f;
2219 }
2220 image->update();
2221 }
2222 else {
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);
2227
2228 for (int z=0; z<nz; z++) {
2229 for (int y=0; y<ny; y++) {
2230 for (int x=0; x<nx; x++) {
2231 float cor=exp(-Util::hypot3sq(x-nx/2,y-ny/2,z-nz/2)*wid);
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);
2233 }
2234 }
2235 }
2236 image->update();
2237 }
2238
2239}
2240
2241
2242
2244{
2245 if (!image) {
2246 LOGWARN("NULL Image");
2247 return;
2248 }
2249
2250 maxval = image->get_attr("maximum");
2251 mean = image->get_attr("mean");
2252 sigma = image->get_attr("sigma");
2253
2254 calc_locals(image);
2255
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();
2260
2261 for (size_t i = 0; i < size; ++i) {
2262 process_pixel(&data[i]);
2263 }
2264 image->update();
2265}
2266
2267
2269{
2270 if (!image) {
2271 LOGWARN("NULL Image");
2272 return;
2273 }
2274
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();
2281 is_complex = image->is_complex();
2282
2283 calc_locals(image);
2284
2285
2286 if (!is_valid()) {
2287 return;
2288 }
2289
2290 float *data = image->get_data();
2291 size_t i = 0;
2292
2293 for (int z = 0; z < nz; z++) {
2294 for (int y = 0; y < ny; y++) {
2295 for (int x = 0; x < nx; x++) {
2296 process_pixel(&data[i], x, y, z);
2297 ++i;
2298 }
2299 }
2300 }
2301 image->update();
2302}
2303
2305 int nx=image->get_xsize();
2306 int ny=image->get_ysize();
2307 int nz=image->get_zsize();
2308
2309 if (nz==1) {
2310 float r;
2311 for (int j=(y<r2?0:y-r2); j<(y+r2>ny?ny:y+r2); j++) {
2312 for (int i=(x<r2?0:x-r2); i<(x+r2>nx?nx:x+r2); i++) {
2313 r=sqrt(float(Util::square(i-x)+Util::square(j-y)));
2314 if (r>r2 && r>r1) continue;
2315 if (r>r1) image->set_value_at(i,j,0,v2*(r-r1)/(r2-r1)+v1*(r2-r)/(r2-r1));
2316 else image->set_value_at(i,j,0,v1);
2317 }
2318 }
2319 }
2320 else {
2321 float r;
2322 for (int k=(z<r2?0:z-r2); k<(z+r2>nz?nz:z+r2); k++) {
2323 for (int j=(y<r2?0:y-r2); j<(y+r2>ny?ny:y+r2); j++) {
2324 for (int i=(x<r2?0:x-r2); i<(x+r2>nx?nx:x+r2); i++) {
2325 r=sqrt(float(Util::square(i-x)+Util::square(j-y)+Util::square(k-z)));
2326 if (r>r2 && r>r1) continue;
2327 if (r>r1) image->set_value_at(i,j,k,v2*(r-r1)/(r2-r1)+v1*(r2-r)/(r2-r1));
2328 else image->set_value_at(i,j,k,v1);
2329 }
2330 }
2331 }
2332 }
2333 image->update();
2334}
2335
2337 if (image->is_complex()) throw ImageFormatException("MaskAzProcessor: target image must be real");
2338
2339 int nx=image->get_xsize();
2340 int ny=image->get_ysize();
2341 int nz=image->get_zsize();
2342
2343 float phicen = params.set_default("phicen",0.0f)*M_PI/180.0;
2344 phicen=Util::angle_norm_pi(phicen);
2345 float phirange = params.set_default("phirange",180.0f)*M_PI/180.0;
2346 float phitrirange = params.set_default("phitrirange",0.0f)*M_PI/180.0;
2347 int phitriangle = params.set_default("phitriangle",0);
2348 float cx = params.set_default("cx",nx/2);
2349 float cy = params.set_default("cy",ny/2);
2350 float zmin = params.set_default("zmin",0);
2351 float zmax = params.set_default("zmax",nz);
2352 float ztri = params.set_default("ztriangle",0.0f);
2353 float inner_radius = params.set_default("inner_radius",0.0f);
2354 float outer_radius = params.set_default("outer_radius",nx+ny);
2355
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);
2362 float val=0.0f;
2363 if (r>inner_radius&&r<=outer_radius) {
2364 float as=Util::angle_sub_2pi(az,phicen);
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);
2368 }
2369 if (r==0 && inner_radius<=0) val=1.0;
2370
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));
2376 }
2377 }
2378 }
2379
2380}
2381
2383{
2384 xc = Util::fast_floor(nx/2.0f) + dx;
2385 yc = Util::fast_floor(ny/2.0f) + dy;
2386 zc = Util::fast_floor(nz/2.0f) + dz;
2387
2388 if (outer_radius < 0) {
2389 outer_radius = nx / 2 + outer_radius +1;
2391 }
2392
2393 if (inner_radius <= 0) {
2395 }
2396}
2397
2398
2400{
2401 if (!image) {
2402 throw NullPointerException("NULL image");
2403 }
2404 int nitems = 0;
2405 float sum = 0;
2406 float *data = image->get_data();
2407 size_t i = 0;
2408
2409 xc = Util::fast_floor(nx/2.0f) + dx;
2410 yc = Util::fast_floor(ny/2.0f) + dy;
2411 zc = Util::fast_floor(nz/2.0f) + dz;
2412
2413 for (int z = 0; z < nz; ++z) {
2414 for (int y = 0; y < ny; ++y) {
2415 for (int x = 0; x < nx; ++x) {
2416 float x1 = sqrt((x - xc) * (x - xc) + (y - yc) * (y - yc) + (z - zc) * (z - zc));
2417 if (x1 <= outer_radius + ring_width && x1 >= outer_radius - ring_width) {
2418 sum += data[i];
2419 ++nitems;
2420 }
2421 ++i;
2422 }
2423 }
2424 }
2425
2426 ring_avg = sum / nitems;
2427}
2428
2430{
2431 if (!image) {
2432 LOGWARN("NULL Image");
2433 return;
2434 }
2435
2436 float minval = params.set_default("minval",0.0f);
2437 float newval = params.set_default("newval",minval);
2438
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();
2443
2444
2445
2446 for (size_t i = 0; i < size; ++i) {
2447 if (data[i]<minval) data[i]=newval;
2448 }
2449 image->update();
2450}
2451
2452
2454{
2455 if (!image) {
2456 LOGWARN("NULL image");
2457 return;
2458 }
2459 if (!image->is_complex()) {
2460 LOGWARN("cannot apply complex processor on a real image. Nothing is done.");
2461 return;
2462 }
2463
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();
2468
2469 image->ri2ap();
2470
2471 for (size_t i = 0; i < size; i += 2) {
2472 process_pixel(data);
2473 data += 2;
2474 }
2475
2476 image->update();
2477 image->ap2ri();
2478}
2479
2480
2481
2483{
2484 if (!image) {
2485 LOGWARN("NULL Image");
2486 return;
2487 }
2488
2489 float *data = image->get_data();
2490
2491 nx = image->get_xsize();
2492 ny = image->get_ysize();
2493 nz = image->get_zsize();
2494
2495 int n = (areasize - 1) / 2;
2497
2498 if (nz > 1) {
2500 }
2501
2502 float *matrix = new float[matrix_size];
2503 kernel = new float[matrix_size];
2504
2505 size_t cpysize = areasize * sizeof(float);
2506 size_t start = (nx * ny + nx + 1) * n;
2507
2508 int xend = nx - n;
2509 int yend = ny - n;
2510
2511 int zstart = n;
2512 int zend = nz - n;
2513
2514 int zbox_start = 0;
2515 int zbox_end = areasize;
2516
2517 if (nz == 1) {
2518 zstart = 0;
2519 zend = 1;
2520 zbox_end = 1;
2521 }
2522
2523 size_t nsec = (size_t)nx * (size_t)ny;
2524 int box_nsec = areasize * areasize;
2525
2526 create_kernel();
2527
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));
2531
2532 size_t k;
2533 for (int z = zstart; z < zend; z++) {
2534 for (int y = n; y < yend; y++) {
2535 for (int x = n; x < xend; x++) {
2536
2537 k = (size_t)z * nsec + y * nx + x;
2538
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);
2543 }
2544 }
2545
2546 process_pixel(&data[k], (float) x, (float) y, (float) z, matrix);
2547 }
2548 }
2549 }
2550
2551 if( matrix )
2552 {
2553 delete[]matrix;
2554 matrix = 0;
2555 }
2556
2557 if( kernel )
2558 {
2559 delete[]kernel;
2560 kernel = 0;
2561 }
2562 image->update();
2563}
2564
2566{
2567 EMData* d = new EMData();
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");
2572 image->add(*d);
2573 image->process_inplace("normalize");
2574 delete d;
2575}
2576
2578{
2579 if (nz == 1) {
2580 memset(kernel, 0, areasize * areasize);
2581 kernel[1] = -0.25f;
2582 kernel[3] = -0.25f;
2583 kernel[5] = -0.25f;
2584 kernel[7] = -0.25f;
2585 kernel[4] = 1;
2586 }
2587 else {
2588 memset(kernel, 0, (size_t)areasize * areasize * areasize);
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;
2595 kernel[13] = 1;
2596 }
2597}
2598
2600{
2601 if (!image) return;
2602
2603 float nsigma = params.set_default("nsigma",3.0f);
2604 int bits = params.set_default("bits",5);
2605 int floatbits = params.set_default("floatbits",-1);
2606 float sigma = image->get_attr("sigma");
2607 float mean = image->get_attr("mean");
2608 float bitmax=pow(2.,bits);
2609
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;
2614
2615 if (floatbits<=0) {
2616 float min=image->get_attr("minimum");
2617 float max=image->get_attr("maximum");
2618
2619 // our max/min are either the actual max/min, or mean+-nsigma*sigma
2620 // whichever produces the smaller range. That way by specifying a large
2621 // nsigma you can get the actual min/max of the image easily
2622 if (mean-nsigma*sigma>min) min=mean-nsigma*sigma;
2623 if (mean+nsigma*sigma<max) max=mean+nsigma*sigma;
2624
2625
2626 for (size_t i=0; i<total_size; i++) {
2627 float newval=floor((image->get_value_at_index(i)-min)*bitmax/(max-min)
2628
2629 );
2630 if (newval<0) newval=0;
2631 if (newval>=bitmax) newval=bitmax-1.0f;
2632 image->set_value_at_index(i,newval);
2633 }
2634 }
2635 else {
2636 // In this version, we keep a specified number of significant bits in the floating point significand, leaving the exponent and sign alone
2637 if (floatbits>22) throw InvalidValueException(floatbits,"floatbits must be <23");
2638 unsigned int bitmask= 0xffffffff << (23-floatbits); // this will leave 1 in the sign and exponent bits, along with the specified number of significand bits
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;
2642 *newvali&=bitmask;
2643 image->set_value_at_index(i,newval);
2644 }
2645 }
2646}
2647
2648
2650 EMData *cpy = process(image);
2651 memcpy(image->get_data(),cpy->get_data(),image->get_size()*sizeof(float));
2652 delete cpy;
2653}
2654
2655// mirrors coordinates at box edges
2656inline int MIRE(int x,int nx) { return x<0?-x:(x>=nx?nx-(x-nx+1):x); }
2657
2659{
2660 if (!image) {
2661 LOGWARN("NULL Image");
2662 return NULL;
2663 }
2664
2665
2666 int nx = image->get_xsize();
2667 int ny = image->get_ysize();
2668 int nz = image->get_zsize();
2669
2670 int dx=1,dy=1,dz=1;
2671 if (params.has_key("radius")) dx=dy=dz=params["radius"];
2672 if (params.has_key("xsize")) dx=params["xsize"];
2673 if (params.has_key("ysize")) dy=params["ysize"];
2674 if (params.has_key("zsize")) dz=params["zsize"];
2675 if (nz==1) dz=0;
2676
2677 int matrix_size = (2*dx+1)*(2*dy+1)*(2*dz+1);
2678
2679 vector<float> array(matrix_size);
2680// image->process_inplace("normalize");
2681
2682 EMData *ret = image->copy_head();
2683
2684 // The old version of this code had a lot of hand optimizations, which likely weren't accomplishing much
2685 // This is much simpler, but relies on the compiler to optimize. May be some cost associated with the new
2686 // edge-mirroring policy which could be hand optomized if necessary
2687 for (int k=0; k<nz; k++) {
2688 for (int j=0; j<ny; j++) {
2689 for (int i=0; i<nx; i++) {
2690
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));
2695 }
2696 }
2697 }
2698 float newv=image->get_value_at(i,j,k);
2699 process_pixel(&newv,array.data(),matrix_size);
2700 ret->set_value_at(i,j,k,newv);
2701 }
2702 }
2703 }
2704
2705 ret->update();
2706
2707 return ret;
2708}
2709
2711{
2712 if (!image) {
2713 LOGWARN("NULL Image");
2714 return;
2715 }
2716
2717 int nz = image->get_zsize();
2718
2719 if (nz > 1) {
2720 LOGERR("%s Processor doesn't support 3D", get_name().c_str());
2721 throw ImageDimensionException("3D model not supported");
2722 }
2723
2724 int nx = image->get_xsize();
2725 int ny = image->get_ysize();
2726
2727 int v1 = params["cal_half_width"];
2728 int v2 = params["fill_half_width"];
2729
2730 int v0 = v1 > v2 ? v1 : v2;
2731
2732 if (v2 <= 0) {
2733 v2 = v1;
2734 }
2735
2736 float *data = image->get_data();
2737
2738 for (int y = v0; y <= ny - v0 - 1; y += v2) {
2739 for (int x = v0; x <= nx - v0 - 1; x += v2) {
2740
2741 float sum = 0;
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];
2745 }
2746 }
2747 float mean = sum / ((v1 * 2 + 1) * (v1 * 2 + 1));
2748
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;
2752 }
2753 }
2754 }
2755 }
2756
2757 image->update();
2758}
2759
2760
2762{
2763 if (!image) {
2764 LOGWARN("NULL Image");
2765 return;
2766 }
2767 int nz = image->get_zsize();
2768
2769 if (nz > 1) {
2770 LOGERR("%s Processor doesn't support 3D", get_name().c_str());
2771 throw ImageDimensionException("3D model not supported");
2772 }
2773
2774 int nx = image->get_xsize();
2775 int ny = image->get_ysize();
2776
2777 float value1 = params["value1"];
2778 float value2 = params["value2"];
2779
2780 int v1 = (int) value1;
2781 int v2 = (int) value2;
2782 if (v2 > v1 / 2) {
2783 LOGERR("invalid value2 '%f' in CutoffBlockProcessor", value2);
2784 return;
2785 }
2786
2787 if (v2 <= 0) {
2788 v2 = v1;
2789 }
2790
2791 float *data = image->get_data();
2792 int y = 0, x = 0;
2793 for (y = 0; y <= ny - v1; y += v1) {
2794 for (x = 0; x <= nx - v1; x += v1) {
2795
2796 EMData *clip = image->get_clip(Region(x, y, v1, v1));
2797 EMData *fft = clip->do_fft();
2798
2799 float *fft_data = fft->get_data();
2800 float sum = 0;
2801 int nitems = 0;
2802
2803 for (int i = -v2; i < v2; i++) {
2804 for (int j = 0; j < v2; j++) {
2805 if (j == 0 && i == 0) {
2806 continue;
2807 }
2808
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]);
2812 nitems++;
2813 }
2814 }
2815 }
2816
2817 if( clip )
2818 {
2819 delete clip;
2820 clip = 0;
2821 }
2822
2823 float mean = sum / nitems;
2824
2825 for (int i = y; i < y + v1; i++) {
2826 for (int j = x; j < x + v1; j++) {
2827 data[i * nx + j] = mean;
2828 }
2829 }
2830 }
2831 }
2832
2833 memset(&data[y * nx], 0, (ny - y) * nx * sizeof(float));
2834
2835 for (int i = 0; i < ny; i++) {
2836 memset(&data[i * nx + x], 0, (nx - x) * sizeof(float));
2837 }
2838
2839 image->update();
2840}
2841
2843{
2844 if (image->is_complex()) throw ImageFormatException("Error, the median shrink processor does not work on complex images");
2845
2846 int shrink_factor = params.set_default("n",0);
2847 if (shrink_factor <= 1) {
2848 throw InvalidValueException(shrink_factor,
2849 "median shrink: shrink factor must > 1");
2850 }
2851
2852 int nx = image->get_xsize();
2853 int ny = image->get_ysize();
2854 int nz = image->get_zsize();
2855
2856// if ((nx % shrink_factor != 0) || (ny % shrink_factor != 0) || (nz > 1 && (nz % shrink_factor != 0))) {
2857// throw InvalidValueException(shrink_factor, "Image size not divisible by shrink factor");
2858// }
2859
2860
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;
2865
2866 EMData* copy = image->copy();
2867 image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
2868 accrue_median(image,copy,shrink_factor);
2869 image->update();
2870 if( copy )
2871 {
2872 delete copy;
2873 copy = 0;
2874 }
2875}
2876
2877//
2879{
2880 if (image->is_complex()) throw ImageFormatException("Error, the median shrink processor does not work on complex images");
2881
2882 int shrink_factor = params.set_default("n",0);
2883 if (shrink_factor <= 1) {
2884 throw InvalidValueException(shrink_factor,
2885 "median shrink: shrink factor must > 1");
2886 }
2887 int nx = image->get_xsize();
2888 int ny = image->get_ysize();
2889 int nz = image->get_zsize();
2890
2891
2892// if ((nx % shrink_factor != 0) || (ny % shrink_factor != 0) || (nz > 1 && (nz % shrink_factor != 0))) {
2893// throw InvalidValueException(shrink_factor, "Image size not divisible by shrink factor");
2894// }
2895
2896
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;
2901
2902// EMData* ret = new EMData(shrunken_nx, shrunken_ny, shrunken_nz);
2903 EMData *ret = image->copy_head();
2904 ret->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
2905
2906 accrue_median(ret,image,shrink_factor);
2907 ret->update();
2908 return ret;
2909}
2910
2911void MedianShrinkProcessor::accrue_median(EMData* to, const EMData* const from,const int shrink_factor)
2912{
2913
2914 int nx_old = from->get_xsize();
2915 int ny_old = from->get_ysize();
2916
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;
2922 }
2923
2924 float *mbuf = new float[threed_shrink_factor];
2925
2926
2927 int nxy_old = nx_old * ny_old;
2928
2929 int nx = to->get_xsize();
2930 int ny = to->get_ysize();
2931 int nz = to->get_zsize();
2932 int nxy_new = nx * ny;
2933
2934 float * rdata = to->get_data();
2935 const float *const data_copy = from->get_const_data();
2936
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;
2941
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;
2946
2947 for (int i = 0; i < nx; i++) {
2948 int i_min = i * shrink_factor;
2949 int i_max = (i + 1) * shrink_factor;
2950
2951 size_t k = 0;
2952 for (int l2 = l_min; l2 < l_max; l2++) {
2953 size_t cur_l2 = l2 * nxy_old;
2954
2955 for (int j2 = j_min; j2 < j_max; j2++) {
2956 size_t cur_j2 = j2 * nx_old + cur_l2;
2957
2958 for (int i2 = i_min; i2 < i_max; i2++) {
2959 mbuf[k] = data_copy[i2 + cur_j2];
2960 ++k;
2961 }
2962 }
2963 }
2964
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]) {
2968 float f = mbuf[i2];
2969 mbuf[i2] = mbuf[k];
2970 mbuf[k] = f;
2971 }
2972 }
2973 }
2974
2975 rdata[i + cur_j] = mbuf[threed_shrink_factor / 2];
2976 }
2977 }
2978 }
2979
2980 if( mbuf )
2981 {
2982 delete[]mbuf;
2983 mbuf = 0;
2984 }
2985
2986 to->scale_pixel((float)shrink_factor);
2987}
2988
2990{
2991 if (image->is_complex()) throw ImageFormatException("FFTResampleProcessor does not support complex images.");
2992
2993 float sample_rate = params.set_default("n",0.0f);
2994 if (sample_rate <= 0.0F ) {
2995 throw InvalidValueException(sample_rate, "downsampling must be >0 ");
2996 }
2997
2998 int nx=image->get_xsize();
2999 int ny=image->get_ysize();
3000 int nz=image->get_zsize();
3001
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));
3005
3006 if (nnx==nx && nny==ny && nnz==nz) return image->copy();
3007
3008 EMData *result;
3009 //downsample
3010 if (nnx<nx||nny<ny||nnz<nz) {
3011 // the type casting here is because FourTruncate was not defined to be const (but it is)
3012 result=((EMData *)image)->FourTruncate(nnx, nny, nnz, 1, 0); // nnx,nny,nnz,returnreal,normalize
3013 } //upscale
3014 else {
3015 // the type casting here is because FourTruncate was not defined to be const (but it is)
3016 result=((EMData *)image)->FourInterpol(nnx, nny, nnz, 1, 0); // nnx,nny,nnz,returnreal,normalize
3017 }
3018 result->scale_pixel((float)nx/(float)nnx);
3019 result->update();
3020 return result;
3021
3022
3023// EMData* result;
3024// if (image->is_complex()) result = image->copy();
3025// else result = image->do_fft();
3026// fft_resample(result,image,sample_rate);
3027 // The image may have been padded - we should shift it so that the phase origin is where FFTW expects it
3028// result->update();
3029// result->scale_pixel(sample_rate);
3030// return result;
3031}
3032
3034{
3035 float sample_rate = params.set_default("n",0.0f);
3036 if (sample_rate <= 0.0F ) {
3037 throw InvalidValueException(sample_rate, "sample rate (n) must be >0 ");
3038 }
3039 //if (image->is_complex()) throw ImageFormatException("Error, the fft resampling processor does not work on complex images");
3040
3041 int nx=image->get_xsize();
3042 int ny=image->get_ysize();
3043 int nz=image->get_zsize();
3044
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));
3048
3049 if (nnx==nx && nny==ny && nnz==nz) return;
3050
3051 EMData *result;
3052 //downsample
3053 if (nnx<nx||nny<ny||nnz<nz) {
3054 // the type casting here is because FourTruncate was not defined to be const (but it is)
3055 result=((EMData *)image)->FourTruncate(nnx, nny, nnz, 1, 0); // nnx,nny,nnz,returnreal,normalize
3056 } //upscale
3057 else {
3058 // the type casting here is because FourTruncate was not defined to be const (but it is)
3059 result=((EMData *)image)->FourInterpol(nnx, nny, nnz, 1, 0); // nnx,nny,nnz,returnreal,normalize
3060 }
3061
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);
3065 image->update();
3066 delete result;
3067
3068 //fft_resample(image,image,sample_rate);
3069
3070 //image->scale_pixel(sample_rate);
3071
3072}
3073
3074void FFTResampleProcessor::fft_resample(EMData* to, const EMData *const from, const float& sample_rate) {
3075 int nx = from->get_xsize();
3076 int ny = from->get_ysize();
3077 int nz = from->get_zsize();
3078
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);
3082
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;
3086
3087 int ndim = from->get_ndim();
3088 if ( ndim < 3 ) {
3089 new_nz = 1;
3090 }
3091 if ( ndim < 2 ) {
3092 new_ny = 1;
3093 }
3094
3095 int fft_x_correction = 1;
3096 if (new_nx % 2 == 0) fft_x_correction = 2;
3097
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;
3101
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;
3105
3106 if ( ! to->is_complex()) to->do_fft_inplace();
3107
3108 if (ndim != 1) to->process_inplace("xform.fourierorigin.tocenter");
3109
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);
3111 to->clip_inplace(clip);
3112
3113 if (fft_x_correction == 1) to->set_fftodd(true);
3114 else to->set_fftodd(false);
3115
3116 if (ndim != 1) to->process_inplace("xform.fourierorigin.tocorner");
3117
3118 to->do_ift_inplace();
3119 to->depad_corner();
3120
3121}
3122
3123
3125{
3126 if (image->is_complex()) throw ImageFormatException("Error, the mean shrink processor does not work on complex images");
3127
3128 if (image->get_ndim() == 1) { throw ImageDimensionException("Error, mean shrink works only for 2D & 3D images"); }
3129
3130 float shrink_factor0 = params.set_default("n",0.0f);
3131 int shrink_factor = int(shrink_factor0);
3132 if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {
3133 throw InvalidValueException(shrink_factor0,
3134 "mean shrink: shrink factor must be >1 integer or 1.5");
3135 }
3136
3137 int nx = image->get_xsize();
3138 int ny = image->get_ysize();
3139 int nz = image->get_zsize();
3140
3141
3142 // here handle the special averaging by 1.5 for 2D case
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");
3145
3146 int shrunken_nx = (int(nx / 1.5)+1)/2*2; // make sure the output size is even
3147 int shrunken_ny = (int(ny / 1.5)+1)/2*2;
3148 EMData* result = new EMData(shrunken_nx,shrunken_ny,1);
3149
3150 accrue_mean_one_p_five(result,image);
3151 result->update();
3152
3153 return result;
3154 }
3155
3156 int shrunken_nx = nx / shrink_factor;
3157 int shrunken_ny = ny / shrink_factor;
3158 int shrunken_nz = 1;
3159
3160 if (nz > 1) {
3161 shrunken_nz = nz / shrink_factor;
3162 }
3163
3164// EMData* result = new EMData(shrunken_nx,shrunken_ny,shrunken_nz);
3165 EMData* result = image->copy_head();
3166 result->set_size(shrunken_nx,shrunken_ny,shrunken_nz);
3167 accrue_mean(result,image,shrink_factor);
3168
3169 result->update();
3170
3171 return result;
3172}
3173
3175{
3176 if (image->is_complex()) throw ImageFormatException("Error, the mean shrink processor does not work on complex images");
3177
3178 if (image->get_ndim() == 1) { throw ImageDimensionException("Error, mean shrink works only for 2D & 3D images"); }
3179
3180 float shrink_factor0 = params.set_default("n",0.0f);
3181 int shrink_factor = int(shrink_factor0);
3182 if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {
3183 throw InvalidValueException(shrink_factor0,
3184 "mean shrink: shrink factor must be >1 integer or 1.5");
3185 }
3186
3187/* if ((nx % shrink_factor != 0) || (ny % shrink_factor != 0) ||
3188 (nz > 1 && (nz % shrink_factor != 0))) {
3189 throw InvalidValueException(shrink_factor,
3190 "Image size not divisible by shrink factor");
3191}*/
3192
3193 int nx = image->get_xsize();
3194 int ny = image->get_ysize();
3195 int nz = image->get_zsize();
3196 // here handle the special averaging by 1.5 for 2D case
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");
3199
3200 int shrunken_nx = (int(nx / 1.5)+1)/2*2; // make sure the output size is even
3201 int shrunken_ny = (int(ny / 1.5)+1)/2*2;
3202
3203 EMData* orig = image->copy();
3204 image->set_size(shrunken_nx, shrunken_ny, 1); // now nx = shrunken_nx, ny = shrunken_ny
3205 image->to_zero();
3206
3207 accrue_mean_one_p_five(image,orig);
3208
3209 if( orig ) {
3210 delete orig;
3211 orig = 0;
3212 }
3213 image->update();
3214
3215 return;
3216 }
3217
3218 accrue_mean(image,image,shrink_factor);
3219
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;
3224
3225 image->update();
3226 image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);
3227}
3228
3229void MeanShrinkProcessor::accrue_mean(EMData* to, const EMData* const from,const int shrink_factor)
3230{
3231 const float * const data = from->get_const_data();
3232 float* rdata = to->get_data();
3233
3234 size_t nx = from->get_xsize();
3235 size_t ny = from->get_ysize();
3236 size_t nz = from->get_zsize();
3237 size_t nxy = nx*ny;
3238
3239
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;
3244
3245 int normalize_shrink_factor = shrink_factor * shrink_factor;
3246 int z_shrink_factor = 1;
3247
3248 if (nz > 1) {
3249 shrunken_nz = nz / shrink_factor;
3250 normalize_shrink_factor *= shrink_factor;
3251 z_shrink_factor = shrink_factor;
3252 }
3253
3254 float invnormfactor = 1.0f/(float)normalize_shrink_factor;
3255
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;
3260
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;
3265
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;
3269
3270 float sum = 0;
3271 for (size_t kk = k_min; kk < k_max; kk++) {
3272 size_t cur_kk = kk * nxy;
3273
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];
3278 }
3279 }
3280 }
3281 rdata[i + cur_j] = sum * invnormfactor;
3282 }
3283 }
3284 }
3285 to->scale_pixel((float)shrink_factor);
3286}
3287
3288
3290{
3291 int nx0 = from->get_xsize(), ny0 = from->get_ysize(); // the original size
3292
3293 int nx = to->get_xsize(), ny = to->get_ysize();
3294
3295 float *data = to->get_data();
3296 const float * const data0 = from->get_const_data();
3297
3298 for (int j = 0; j < ny; j++) {
3299 int jj = int(j * 1.5);
3300 float jw0 = 1.0F, jw1 = 0.5F; // 3x3 -> 2x2, so each new pixel should have 2.25 of the old pixels
3301 if ( j%2 ) {
3302 jw0 = 0.5F;
3303 jw1 = 1.0F;
3304 }
3305 for (int i = 0; i < nx; i++) {
3306 int ii = int(i * 1.5);
3307 float iw0 = 1.0F, iw1 = 0.5F;
3308 float w = 0.0F;
3309
3310 if ( i%2 ) {
3311 iw0 = 0.5F;
3312 iw1 = 1.0F;
3313 }
3314 if ( jj < ny0 ) {
3315 if ( ii < nx0 ) {
3316 data[j * nx + i] = data0[ jj * nx0 + ii ] * jw0 * iw0 ;
3317 w += jw0 * iw0 ;
3318 if ( ii+1 < nx0 ) {
3319 data[j * nx + i] += data0[ jj * nx0 + ii + 1] * jw0 * iw1;
3320 w += jw0 * iw1;
3321 }
3322 }
3323 if ( jj +1 < ny0 ) {
3324 if ( ii < nx0 ) {
3325 data[j * nx + i] += data0[ (jj+1) * nx0 + ii ] * jw1 * iw0;
3326 w += jw1 * iw0;
3327 if ( ii+1 < nx0 ) {
3328 data[j * nx + i] += data0[ (jj+1) * nx0 + ii + 1] * jw1 * iw1;
3329 w += jw1 * iw1;
3330 }
3331 }
3332 }
3333 }
3334 if ( w>0 ) data[j * nx + i] /= w;
3335 }
3336 }
3337
3338 to->update();
3339 to->scale_pixel((float)1.5);
3340}
3341
3342// This would have to be moved into the header if it were required in other source files
3343template<class LogicOp>
3345{
3346 // The basic idea of this code is to iterate through each pixel in the output image
3347 // determining its value by investigation a region of the input image
3348
3349 if (!image) throw NullPointerException("Attempt to max shrink a null image");
3350
3351 if (image->is_complex() ) throw ImageFormatException("Can not max shrink a complex image");
3352
3353
3354 int shrink = params.set_default("n",2);
3355 int search = params.set_default("search",2);
3356
3357 if ( shrink < 0 ) throw InvalidValueException(shrink, "Can not shrink by a value less than 0");
3358
3359
3360 int nz = image->get_zsize();
3361 int ny = image->get_ysize();
3362 int nx = image->get_xsize();
3363
3364 if (nx == 1 && ny == 1 && nz == 1 ) return image->copy();
3365
3366 LogicOp op;
3367 EMData* return_image = new EMData();
3368
3369 int shrinkx = shrink;
3370 int shrinky = shrink;
3371 int shrinkz = shrink;
3372
3373 int searchx = search;
3374 int searchy = search;
3375 int searchz = search;
3376
3377 // Clamping the shrink values to the dimension lengths
3378 // ensures that the return image has non zero dimensions
3379 if ( shrinkx > nx ) shrinkx = nx;
3380 if ( shrinky > ny ) shrinky = ny;
3381 if ( shrinkz > nz ) shrinkz = nz;
3382
3383 if ( nz == 1 && ny == 1 )
3384 {
3385 return_image->set_size(nx/shrinkx);
3386 for(int i = 0; i < nx/shrinkx; ++i)
3387 {
3388 float tmp = op.get_start_val();
3389 for(int s=0; s < searchx; ++s)
3390 {
3391 int idx = shrinkx*i+s;
3392 // Don't ask for memory beyond limits
3393 if ( idx > nx ) break;
3394 else
3395 {
3396 float val = image->get_value_at(idx);
3397 if ( op( val,tmp) ) tmp = val;
3398 }
3399 }
3400 return_image->set_value_at(i,tmp);
3401 }
3402 }
3403 else if ( nz == 1 )
3404 {
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;
3417
3418 float val = image->get_value_at(xidx,yidx);
3419 if ( op( val,tmp) ) tmp = val;
3420 }
3421 }
3422 return_image->set_value_at(x,y,tmp);
3423 }
3424 }
3425 }
3426 else
3427 {
3428 int tz = nz/shrinkz;
3429 int ty = ny/shrinky;
3430 int tx = nx/shrinkx;
3431
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();
3437
3438 for(int sz=0; sz < searchz; ++sz) {
3439 int zidx = shrinkz*z+sz;
3440 if ( zidx >= nz) break;
3441
3442 for(int sy=0; sy < searchy; ++sy) {
3443 int yidx = shrinky*y+sy;
3444 if ( yidx >= ny) break;
3445
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;
3451 }
3452 }
3453 }
3454 return_image->set_value_at(x,y,z,tmp);
3455 }
3456 }
3457 }
3458 }
3459 return_image->update();
3460
3461 return return_image;
3462}
3463
3464template<class LogicOp>
3466{
3467 // The basic idea of this code is to iterate through each pixel in the output image
3468 // determining its value by investigation a region of the input image
3469 if (!image) throw NullPointerException("Attempt to max shrink a null image");
3470
3471 if (image->is_complex() ) throw ImageFormatException("Can not max shrink a complex image");
3472
3473
3474 int shrink = params.set_default("shrink",2);
3475 int search = params.set_default("search",2);
3476
3477 if ( shrink < 0 ) throw InvalidValueException(shrink, "Can not shrink by a value less than 0");
3478
3479
3480 int nz = image->get_zsize();
3481 int ny = image->get_ysize();
3482 int nx = image->get_xsize();
3483
3484 LogicOp op;
3485
3486 int shrinkx = shrink;
3487 int shrinky = shrink;
3488 int shrinkz = shrink;
3489
3490 int searchx = search;
3491 int searchy = search;
3492 int searchz = search;
3493
3494 // Clamping the shrink values to the dimension lengths
3495 // ensures that the return image has non zero dimensions
3496 if ( shrinkx > nx ) shrinkx = nx;
3497 if ( shrinky > ny ) shrinky = ny;
3498 if ( shrinkz > nz ) shrinkz = nz;
3499
3500 if (nx == 1 && ny == 1 && nz == 1 ) return;
3501
3502 if ( nz == 1 && ny == 1 )
3503 {
3504 for(int i = 0; i < nx/shrink; ++i)
3505 {
3506 float tmp = op.get_start_val();
3507 for(int s=0; s < searchx; ++s)
3508 {
3509 int idx = shrinkx*i+s;
3510 if ( idx > nx ) break;
3511 else
3512 {
3513 float val = image->get_value_at(idx);
3514 if ( op( val,tmp) ) tmp = val;
3515 }
3516 }
3517 image->set_value_at(i,tmp);
3518 }
3519
3520 image->set_size(nx/shrinkx);
3521 }
3522 else if ( nz == 1 )
3523 {
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;
3535
3536 float val = image->get_value_at(xidx,yidx);
3537 if ( op( val,tmp) ) tmp = val;
3538 }
3539 }
3540 (*image)(x+tx*y) = tmp;
3541 }
3542 }
3543 image->set_size(tx,ty);
3544 }
3545 else
3546 {
3547 int tnxy = nx/shrinkx*ny/shrinky;
3548 int tz = nz/shrinkz;
3549 int ty = ny/shrinky;
3550 int tx = nx/shrinkx;
3551
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;
3565
3566 float val = image->get_value_at(xidx,yidx,zidx);
3567 if ( op( val,tmp) ) tmp = val;
3568 }
3569 }
3570 }
3571 (*image)(x+tx*y+tnxy*z) = tmp;
3572 }
3573 }
3574 }
3575 image->set_size(tx,ty,tz);
3576 }
3577 image->update();
3578}
3579
3580
3581
3583{
3584 if (!image) {
3585 LOGWARN("NULL Image");
3586 return;
3587 }
3588
3589 int nz = image->get_zsize();
3590 if (nz > 1) {
3591 LOGERR("%s Processor doesn't support 3D model", get_name().c_str());
3592 throw ImageDimensionException("3D model not supported");
3593 }
3594
3595 int nx = image->get_xsize();
3596 int ny = image->get_ysize();
3597 float *dy = new float[ny];
3598 float m = 0;
3599 float b = 0;
3600 float sum_y = 0;
3601 float *data = image->get_data();
3602
3603 for (int i = 0; i < ny; i++) {
3604 Util::calc_least_square_fit(nx, 0, data + i * nx, &m, &b, false);
3605 dy[i] = b;
3606 sum_y += m;
3607 }
3608
3609 float mean_y = sum_y / ny;
3610 float sum_x = 0;
3611 Util::calc_least_square_fit(ny, 0, dy, &sum_x, &b, false);
3612
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;
3616 }
3617 }
3618
3619 image->update();
3620}
3621
3623{
3624
3625 EMData* mask = params.set_default("mask",(EMData*)0);
3626 int radius = params.set_default("radius",0);
3627
3628 if (radius != 0 && mask != 0) throw InvalidParameterException("Error - the mask and radius parameters are mutually exclusive.");
3629
3630 if (mask == 0 && radius == 0) throw InvalidParameterException("Error - you must specify either the mask or the radius parameter.");
3631
3632 // If the radius isn't 0, then turn the mask into the thing we want...
3633 bool deletemask = false;
3634 if (radius != 0) {
3635 mask = new EMData;
3636 int n = image->get_ndim();
3637 if (n==1){
3638 mask->set_size(2*radius+1);
3639 } else if (n==2) {
3640 mask->set_size(2*radius+1,2*radius+1);
3641 }
3642 else /*n==3*/ {
3643 mask->set_size(2*radius+1,2*radius+1,2*radius+1);
3644 }
3645 // assuming default behavior is to make a circle/sphere with using the radius of the mask
3646 mask->process_inplace("testimage.circlesphere");
3647 }
3648
3649 // Double check that that mask isn't too big
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; // Sanity check
3654 if (ny == 1) nyc = 1; // Sanity check
3655
3656 if ( mnx > nx || mny > ny || mnz > nz)
3657 throw ImageDimensionException("Can not flatten using a mask that is larger than the image.");
3658
3659 // Get the normalization factor
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);
3663 }
3664 // If the sum is zero the user probably doesn't understand that determining a measure of the mean requires
3665 // strictly positive numbers. The user has specified a mask that consists entirely of zeros, or the mask
3666 // has a mean of zero.
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;
3669
3670 // The mask can now be automatically resized to the dimensions of the image
3671// bool undoclip = false;
3672
3673 Region r;
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);
3677 mask->clip_inplace(r,0);
3678// undoclip = true;
3679// if ( mnx < nx || mny < ny || mnz < nz) {
3680// Region r((mnx-nx)/2, (mny-ny)/2,(mnz-nz)/2,nx,ny,nz);
3681// mask->clip_inplace(r);
3682// undoclip = true;
3683// }
3684
3685 Region r2;
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);
3689 image->clip_inplace(r2,image->get_edge_mean());
3690 // Finally do the convolution
3691 EMData* m = image->convolute(mask);
3692 // Normalize so that m is truly the local mean
3693 m->mult(normfac);
3694 // Before we can subtract, the mean must be phase shifted
3695 m->process_inplace("xform.phaseorigin.tocenter");
3696 // Subtract the local mean
3697// image->write_image("a.mrc");
3698// m->write_image("b.mrc");
3699 image->sub(*m); // WE'RE DONE!
3700 delete m;
3701
3702 if (deletemask) {
3703 delete mask;
3704 } else { // I clipped it inplace, so undo this clipping so the user gets back what the put in
3705 Region r;
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);
3709 mask->clip_inplace(r);
3710 }
3711
3712 Region r3;
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);
3716 image->clip_inplace(r3);
3717// if ( undoclip ) {
3718// Region r((nx-mnx)/2, (ny-mny)/2, (nz-mnz)/2,mnx,mny,mnz);
3719// mask->clip_inplace(r);
3720// }
3721
3722}
3723
3725 if (!image) { LOGWARN("NULL IMAGE"); return; }
3726 //int isinten=image->get_attr_default("is_intensity",0);
3727
3728 // 1-D
3729 if (image->get_ysize()==1) {
3730
3731 }
3732 // 2-D
3733 else if (image->get_zsize()==1) {
3734// if (!isinten) throw ImageDimensionException("Only complex intensity images currently supported by NonConvexProcessor");
3735 int nx2=image->get_xsize()/2;
3736 int ny2=image->get_ysize()/2;
3737 vector<float> rdist = image->calc_radial_dist(nx2*1.5,0,1,false); // radial distribution to make sure nonconvex values decrease radially
3738 // Make sure rdist is decreasing (or flat)
3739 for (int i=1; i<nx2; i++) {
3740 if (rdist[i]>rdist[i-1]) rdist[i]=rdist[i-1];
3741 }
3742
3743 image->process_inplace("xform.fourierorigin.tocenter");
3744 EMData* binary=image->copy();
3745
3746 // First we eliminate convex points from the input image (set to zero)
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])) { // point is considered nonconvex if lower than surrounding values and lower than mean
3752 image->set_value_at_fast(x/2+nx2,y,0.0); // we are turning image into a full real-space intensity image for now
3753 image->set_value_at_fast(nx2-x/2,ny2*2-y-1,0.0);
3754 }
3755 else {
3756 image->set_value_at_fast(x/2+nx2,y,cen); // we are turning image into a full real-space intensity image for now
3757 image->set_value_at_fast(nx2-x/2,ny2*2-y-1,cen); // It will contain non-zero values only for nonconvex points
3758 }
3759 }
3760 }
3761 image->set_value_at_fast(nx2+1,ny2,(*binary)(2,ny2)); // We keep the points near the Fourier origin as a central anchor even though it's convex
3762 image->set_value_at_fast(nx2-1,ny2,(*binary)(2,ny2)); // We keep the points near the Fourier origin as a central anchor even though it's convex
3763 image->set_value_at_fast(nx2,ny2+1,(*binary)(0,ny2+1)); // We keep the points near the Fourier origin as a central anchor even though it's convex
3764 image->set_value_at_fast(nx2,ny2-1,(*binary)(0,ny2-1)); // We keep the points near the Fourier origin as a central anchor even though it's convex
3765 for (int y=0; y<ny2*2; y++) image->set_value_at_fast(0,y,0.0f);
3766
3767 // Now make a binary version of the convex points
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; // binary version of the convex points in image
3773 }
3774 binary->update();
3775
3776 // We now use a Gaussian filter on both images, to use Gaussian interpolation to fill in zero values
3777 image->set_complex(false); // so we can use a Gaussian filter on it
3778 binary->set_complex(false);
3779
3780/* image->write_image("con.hdf",0);*/
3781 image->set_fftpad(false);
3782 binary->set_fftpad(false);
3783
3784 // Gaussian blur of both images
3785 image->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.04f));
3786 binary->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.04f));
3787
3788/* image->write_image("con.hdf",1);
3789 binary->write_image("con.hdf",2);*/
3790
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);
3796 }
3797 }
3798 image->set_complex(true);
3799 image->set_fftpad(true);
3800 image->process_inplace("xform.fourierorigin.tocorner");
3801 delete binary;
3802 }
3803 else throw ImageDimensionException("3D maps not yet supported by NonConvexProcessor");
3804
3805}
3806
3807
3808#include <gsl/gsl_linalg.h>
3810{
3811 if (!image) {
3812 LOGWARN("NULL Image");
3813 return;
3814 }
3815
3816 int nz = image->get_zsize();
3817 if (nz > 1) {
3818 LOGERR("%s Processor doesn't support 3D model", get_name().c_str());
3819 throw ImageDimensionException("3D map not supported");
3820 }
3821
3822 int nx = image->get_xsize();
3823 int ny = image->get_ysize();
3824 float *d = image->get_data();
3825 EMData *mask = 0;
3826 float *dm = 0;
3827 if (params.has_key("mask")) {
3828 mask = params["mask"];
3829 if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) {
3830 LOGERR("%s Processor requires same size mask image", get_name().c_str());
3831 throw ImageDimensionException("wrong size mask image");
3832 }
3833 dm = mask->get_data();
3834 }
3835 int count = 0;
3836 if (dm) {
3837 for(int i=0; i<nx*ny; i++) {
3838 if(dm[i]) count++;
3839 }
3840 }
3841 else {
3842 count = nx * ny;
3843 }
3844 if(count<3) {
3845 LOGERR("%s Processor requires at least 3 pixels to fit a plane", get_name().c_str());
3846 throw ImageDimensionException("too few usable pixels to fit a plane");
3847 }
3848 // Allocate the working space
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);
3852
3853 double m[3] = {0, 0, 0};
3854 int index=0;
3855 if (dm) {
3856 for(int j=0; j<ny; j++){
3857 for(int i=0; i<nx; i++){
3858 int ij=j*nx+i;
3859 if(dm[ij]) {
3860 m[0]+=i; // x
3861 m[1]+=j; // y
3862 m[2]+=d[ij]; // z
3863 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \
3864 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/
3865 index++;
3866 }
3867 }
3868 }
3869 }
3870 else {
3871 for(int j=0; j<ny; j++){
3872 for(int i=0; i<nx; i++){
3873 int ij=j*nx+i;
3874 m[0]+=i; // x
3875 m[1]+=j; // y
3876 m[2]+=d[ij]; // z
3877 /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \
3878 index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/
3879 index++;
3880 }
3881 }
3882 }
3883
3884 for(int i=0; i<3; i++) m[i]/=count; // compute center of the plane
3885
3886 index=0;
3887 if (dm) {
3888 for(int j=0; j<ny; j++){
3889 for(int i=0; i<nx; i++){
3890 int ij=j*nx+i;
3891 if(dm[ij]) {
3892 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]);
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]);
3896 index++;
3897 }
3898 }
3899 }
3900 mask->update();
3901 }
3902 else {
3903 for(int j=0; j<ny; j++){
3904 for(int i=0; i<nx; i++){
3905 int ij=j*nx+i;
3906 //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]);
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]);
3910 index++;
3911 }
3912 }
3913 }
3914
3915 // SVD decomposition and use the V vector associated with smallest singular value as the plan normal
3916 gsl_linalg_SV_decomp_jacobi(A, V, S);
3917
3918 double n[3];
3919 for(int i=0; i<3; i++) n[i] = gsl_matrix_get(V, i, 2);
3920
3921 #ifdef DEBUG
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]);
3927 #endif
3928
3929 int changeZero = 0;
3930 if (params.has_key("changeZero")) changeZero = params["changeZero"];
3931 if (changeZero) {
3932 for(int j=0; j<nx; j++){
3933 for(int i=0; i<ny; i++){
3934 int ij = j*nx+i;
3935 d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);
3936 }
3937 }
3938 }
3939 else {
3940 for(int j=0; j<nx; j++){
3941 for(int i=0; i<ny; i++){
3942 int ij = j*nx+i;
3943 if(d[ij]) d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);
3944 }
3945 }
3946 }
3947 image->update();
3948 // set return plane parameters
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]);
3953 params["planeParam"]=EMObject(planeParam);
3954}
3955
3957{
3958 if (!image) {
3959 LOGWARN("NULL Image");
3960 return;
3961 }
3962
3963 int nx = image->get_xsize();
3964 int ny = image->get_ysize();
3965 int nz = image->get_zsize();
3966
3967 float *data = image->get_data();
3968 float sigma = image->get_attr("sigma");
3969
3970 for (int k = 0; k < nz; k++) {
3971 for (int i = 0; i < nx; i++) {
3972 double sum = 0;
3973 for (int j = ny / 4; j < 3 * ny / 4; j++) {
3974 sum += data[i + j * nx];
3975 }
3976
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;
3980 }
3981 }
3982 }
3983
3984 image->update();
3985}
3986
3988{
3989 if (!image) {
3990 LOGWARN("NULL Image");
3991 return;
3992 }
3993
3994 //Note : real image only!
3995 if(image->is_complex()) {
3996 LOGERR("%s Processor only operates on real images", get_name().c_str());
3997 throw ImageFormatException("apply to real image only");
3998 }
3999
4000 // Note : 2D only!
4001 int nz = image->get_zsize();
4002 if (nz > 1) {
4003 LOGERR("%s Processor doesn't support 3D models", get_name().c_str());
4004 throw ImageDimensionException("3D model not supported");
4005 }
4006
4007 EMData *ff=image->do_fft();
4008 ff->ri2ap();
4009
4010 int nx=image->get_xsize();
4011 int ny=image->get_ysize();
4012
4013 int x,y;
4014 float norm=static_cast<float>(nx*ny);
4015
4016 for (y=0; y<ny; y++) image->set_value_at(0,y,0);
4017
4018 for (x=1; x<nx/2; x++) {
4019 for (y=0; y<ny; y++) {
4020 int y2;
4021 if (y<ny/2) y2=y+ny/2;
4022 else if (y==ny/2) y2=ny;
4023 else y2=y-ny/2;
4024 image->set_value_at(x,y,ff->get_value_at(nx-x*2,ny-y2)/norm);
4025 }
4026 }
4027
4028 for (x=nx/2; x<nx; x++) {
4029 for (y=0; y<ny; y++) {
4030 int y2;
4031 if (y<ny/2) y2=y+ny/2;
4032 else y2=y-ny/2;
4033 image->set_value_at(x,y,ff->get_value_at(x*2-nx,y2)/norm);
4034 }
4035 }
4036
4037 image->update();
4038 if( ff )
4039 {
4040 delete ff;
4041 ff = 0;
4042 }
4043}
4044
4046{
4047 if (!image) {
4048 LOGWARN("NULL Image");
4049 return;
4050 }
4051
4052 if (image->get_zsize() > 1) {
4053 LOGERR("%s Processor doesn't support 3D model", get_name().c_str());
4054 throw ImageDimensionException("3D model not supported");
4055 }
4056
4057 float nonzero = params.set_default("nonzero",false);
4058
4059 float *d = image->get_data();
4060 int i = 0;
4061 int j = 0;
4062
4063 int nx = image->get_xsize();
4064 int ny = image->get_ysize();
4065
4066 float zval=9.99e23f; // we're assuming we won't run into this exact value for an edge, not great programming, but good enough
4067 if (nonzero) {
4068 int x,y;
4069 size_t corn=nx*ny-1;
4070
4071 // this set of 4 tests looks for any edges with exactly the same value
4072 for (x=1; x<nx; x++) { if (d[x]!=d[0]) break;}
4073 if (x==nx) zval=d[0];
4074
4075 for (y=1; y<ny; y++) { if (d[y*nx]!=d[0]) break; }
4076 if (y==ny) zval=d[0];
4077
4078 for (x=1; x<nx; x++) { if (d[corn-x]!=d[corn]) break;}
4079 if (x==nx) zval=d[corn];
4080
4081 for (y=1; y<ny; y++) { if (d[corn-y*nx]!=d[corn]) break; }
4082 if (y==ny) zval=d[corn];
4083
4084 if (zval!=9.99e23f) {
4085 image->set_attr("hadzeroedge",1);
4086// printf("zeroedge %f\n",zval);
4087
4088 }
4089 else image->set_attr("hadzeroedge",0);
4090
4091 // This tries to detect images where the edges have been filled with the nearest non-zero value. The filter does nothing, but we set the tag.
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;
4094 }
4095 if (x==nx/2+5) image->set_attr("hadzeroedge",2);
4096
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;
4099 }
4100 if (x==nx/2+5) image->set_attr("hadzeroedge",2);
4101
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;
4104 }
4105 if (y==ny/2+5) image->set_attr("hadzeroedge",2);
4106
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;
4109 }
4110 if (y==ny/2+5) image->set_attr("hadzeroedge",2);
4111
4112 }
4113 if (zval==9.99e23f) zval=0;
4114
4115 for (j = 0; j < ny; j++) {
4116 for (i = 0; i < nx - 1; i++) {
4117 if (d[i + j * nx] != zval) {
4118 break;
4119 }
4120 }
4121
4122 float v = d[i + j * nx];
4123 while (i >= 0) {
4124 d[i + j * nx] = v;
4125 i--;
4126 }
4127
4128 for (i = nx - 1; i > 0; i--) {
4129 if (d[i + j * nx] != zval)
4130 break;
4131 }
4132 v = d[i + j * nx];
4133 while (i < nx) {
4134 d[i + j * nx] = v;
4135 i++;
4136 }
4137 }
4138
4139 for (i = 0; i < nx; i++) {
4140 for (j = 0; j < ny; j++) {
4141 if (d[i + j * nx] != zval)
4142 break;
4143 }
4144
4145 float v = d[i + j * nx];
4146 while (j >= 0) {
4147 d[i + j * nx] = v;
4148 j--;
4149 }
4150
4151 for (j = ny - 1; j > 0; j--) {
4152 if (d[i + j * nx] != zval)
4153 break;
4154 }
4155 v = d[i + j * nx];
4156 while (j < ny) {
4157 d[i + j * nx] = v;
4158 j++;
4159 }
4160 }
4161
4162
4163 image->update();
4164}
4165
4167{
4168 if (!image) {
4169 LOGWARN("NULL Image");
4170 return;
4171 }
4172
4173 bool fix_zero=(bool)params.set_default("fix_zero",0);
4174 float sigmamult=(float)params.set_default("sigma",3.0);
4175
4176 if (sigmamult<=0.0) throw InvalidValueException(sigmamult,"threshold.outlier.localmean: sigma must be >0");
4177
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;
4180
4181 int nx=image->get_xsize();
4182 int ny=image->get_ysize();
4183 int nz=image->get_zsize();
4184
4185 // This isn't optimally efficient
4186 EMData *im[2];
4187 im[0]=image;
4188 im[1]=image->copy_head();
4189
4190 if (nz==1) {
4191 int repeat=1;
4192 while (repeat) {
4193 memcpy(im[1]->get_data(),im[0]->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
4194 repeat=0;
4195 for (int y=0; y<ny; y++) {
4196 for (int x=0; x<nx; x++) {
4197 // if the pixel is an outlier
4198 float pix=im[1]->get_value_at(x,y);
4199 if (pix>hithr || pix<lothr || (pix==0 && fix_zero)) {
4200 int y0=0>y-1?0:y-1;
4201 int y1=y>=ny-1?ny-1:y+1;
4202 int x0=0>x-1?0:x-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;
4209 c+=lpix;
4210 nc++;
4211 }
4212 }
4213 if (nc!=0) im[0]->set_value_at(x,y,c/nc);
4214 else repeat=1;
4215 }
4216 }
4217 }
4218 }
4219 }
4220 else {
4221 throw ImageDimensionException("threshold.outlier.localmean: 3D not yet implemented");
4222 }
4223 delete im[1];
4224}
4225
4226
4228{
4229 if (!image) {
4230 LOGWARN("NULL Image");
4231 return;
4232 }
4233 if (image->get_zsize() > 1) {
4234 LOGERR("BeamstopProcessor doesn't support 3D model");
4235 throw ImageDimensionException("3D model not supported");
4236 }
4237
4238 float value1 = params["value1"];
4239 float value2 = params["value2"];
4240 float value3 = params["value3"];
4241
4242 float thr = fabs(value1);
4243 float *data = image->get_data();
4244 int cenx = (int) value2;
4245 int ceny = (int) value3;
4246
4247 int nx = image->get_xsize();
4248 int ny = image->get_ysize();
4249
4250 if (cenx <= 0) {
4251 cenx = nx / 2;
4252 }
4253
4254 if (ceny <= 0) {
4255 ceny = ny / 2;
4256 }
4257
4258 int mxr = (int) floor(sqrt(2.0f) * nx / 2);
4259
4260 float *mean_values = new float[mxr];
4261 float *sigma_values = new float[mxr];
4262 double sum = 0;
4263 int count = 0;
4264 double square_sum = 0;
4265
4266 for (int i = 0; i < mxr; i++) {
4267 sum = 0;
4268 count = 0;
4269 square_sum = 0;
4270 int nitems = 6 * i + 2;
4271
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);
4276
4277 if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny) {
4278 continue;
4279 }
4280
4281 float f = data[x0 + y0 * nx];
4282 sum += f;
4283 square_sum += f * f;
4284 count++;
4285 }
4286
4287 mean_values[i] = (float)sum / count;
4288 sigma_values[i] = (float) sqrt(square_sum / count - mean_values[i] * mean_values[i]);
4289 }
4290
4291
4292 for (int k = 0; k < 5; k++) {
4293 for (int i = 0; i < mxr; i++) {
4294 sum = 0;
4295 count = 0;
4296 square_sum = 0;
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];
4300
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);
4305
4306 if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny ||
4307 data[x0 + y0 * nx] < thr1 || data[x0 + y0 * nx] > thr2) {
4308 continue;
4309 }
4310
4311 sum += data[x0 + y0 * nx];
4312 square_sum += data[x0 + y0 * nx] * data[x0 + y0 * nx];
4313 count++;
4314 }
4315
4316 mean_values[i] = (float) sum / count;
4317 sigma_values[i] = (float) sqrt(square_sum / count - mean_values[i] * mean_values[i]);
4318 }
4319 }
4320
4321 for (int i = 0; i < nx; i++) {
4322 for (int j = 0; j < ny; j++) {
4323
4324 int r = Util::round(hypot((float) i - cenx, (float) j - ceny));
4325
4326 if (value1 < 0) {
4327 if (data[i + j * nx] < (mean_values[r] - sigma_values[r] * thr)) {
4328 data[i + j * nx] = 0;
4329 }
4330 else {
4331 data[i + j * nx] -= mean_values[r];
4332 }
4333 continue;
4334 }
4335 if (data[i + j * nx] > (mean_values[r] - sigma_values[r] * thr)) {
4336 continue;
4337 }
4338 data[i + j * nx] = mean_values[r];
4339 }
4340 }
4341
4342 if( mean_values )
4343 {
4344 delete[]mean_values;
4345 mean_values = 0;
4346 }
4347
4348 if( sigma_values )
4349 {
4350 delete[]sigma_values;
4351 sigma_values = 0;
4352 }
4353
4354 image->update();
4355}
4356
4357
4358
4360{
4361 if (!image) {
4362 LOGWARN("NULL Image");
4363 return;
4364 }
4365 if (image->get_zsize() > 1) {
4366 LOGERR("MeanZeroEdgeProcessor doesn't support 3D model");
4367 throw ImageDimensionException("3D model not supported");
4368 }
4369
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");
4374
4375 float *d = image->get_data();
4376 int i = 0;
4377 int j = 0;
4378
4379 for (j = 0; j < ny; j++) {
4380 for (i = 0; i < nx - 1; i++) {
4381 if (d[i + j * nx] != 0) {
4382 break;
4383 }
4384 }
4385
4386 if (i == nx - 1) {
4387 i = -1;
4388 }
4389
4390 float v = d[i + j * nx] - mean_nonzero;
4391
4392 while (i >= 0) {
4393 v *= 0.9f;
4394 d[i + j * nx] = v + mean_nonzero;
4395 i--;
4396 }
4397
4398
4399 for (i = nx - 1; i > 0; i--) {
4400 if (d[i + j * nx] != 0) {
4401 break;
4402 }
4403 }
4404
4405 if (i == 0) {
4406 i = nx;
4407 }
4408
4409 v = d[i + j * nx] - mean_nonzero;
4410
4411 while (i < nx) {
4412 v *= .9f;
4413 d[i + j * nx] = v + mean_nonzero;
4414 i++;
4415 }
4416 }
4417
4418
4419 for (i = 0; i < nx; i++) {
4420 for (j = 0; j < ny; j++) {
4421 if (d[i + j * nx] != 0)
4422 break;
4423 }
4424
4425 float v = d[i + j * nx] - mean_nonzero;
4426
4427 while (j >= 0) {
4428 v *= .9f;
4429 d[i + j * nx] = v + mean_nonzero;
4430 j--;
4431 }
4432
4433 for (j = ny - 1; j > 0; j--) {
4434 if (d[i + j * nx] != 0)
4435 break;
4436 }
4437
4438 v = d[i + j * nx] - mean_nonzero;
4439
4440 while (j < ny) {
4441 v *= .9f;
4442 d[i + j * nx] = v + mean_nonzero;
4443 j++;
4444 }
4445 }
4446
4447 image->update();
4448}
4449
4450
4451
4453{
4454 if (!image) {
4455 LOGWARN("NULL Image");
4456 return;
4457 }
4458
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;
4464
4465 size_t idx;
4466 for (int z = 0; z < nz; z++) {
4467 for (int x = 0; x < nx; x++) {
4468 double sum = 0;
4469 for (int y = 0; y < ny; y++) {
4470 idx = x + y * nx + z * nxy;
4471 sum += data[idx];
4472 }
4473 float mean = (float) sum / ny;
4474
4475 for (int y = 0; y < ny; y++) {
4476 idx = x + y * nx + z * nxy;
4477 data[idx] = mean;
4478 }
4479 }
4480 }
4481
4482 image->update();
4483}
4484
4486{
4487 if (!image) {
4488 LOGWARN("NULL Image");
4489 return;
4490 }
4491
4492 int nx = image->get_xsize();
4493 int ny = image->get_ysize();
4494
4495 float *d = image->get_data();
4496 int width = params["width"];
4497
4498 if (width > min(nx,ny)/2.){
4499 LOGERR("width parameter cannot be greater than min(nx,ny)/2");
4500 throw InvalidParameterException("width cannot be greater than min(nx,ny)/2");
4501 }
4502
4503 if (image->get_zsize() > 1){
4504 for (int k=0; k<image->get_zsize(); k++){
4505 int zidx = k*nx*ny;
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;
4511 }
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;
4515 }
4516 }
4517 }
4518 }
4519 else {
4520
4521 for (int i=0; i<width; i++) {
4522 float frac=i/(float)width;
4523 for (int j=0; j<nx; j++) {
4524 d[j+i*nx]*=frac;
4525 d[nx*ny-j-i*nx-1]*=frac;
4526 }
4527 for (int j=0; j<ny; j++) {
4528 d[j*nx+i]*=frac;
4529 d[nx*ny-j*nx-i-1]*=frac;
4530 }
4531 }
4532 }
4533
4534 image->update();
4535}
4536
4538{
4539 if (!image) {
4540 LOGWARN("NULL Image");
4541 return;
4542 }
4543
4544 if (image->get_zsize() > 1) {
4545 LOGERR("ZeroEdgeRowProcessor is not supported in 3D models");
4546 throw ImageDimensionException("3D model not supported");
4547 }
4548
4549 int nx = image->get_xsize();
4550 int ny = image->get_ysize();
4551
4552 float *d = image->get_data();
4553 int apodize = params.set_default("apodize",0);
4554 if (apodize>0) throw InvalidValueException(apodize," apodize not yet supported.");
4555
4556 int top_nrows = std::max(0,params.set_default("y0",0)-apodize);
4557 int bottom_nrows = std::max(0,params.set_default("y1",top_nrows)-apodize);
4558
4559 int left_ncols = std::max(0,params.set_default("x0",0)-apodize);
4560 int right_ncols = std::max(0,params.set_default("x1",left_ncols)-apodize);
4561
4562
4563
4564 size_t row_size = nx * sizeof(float);
4565
4566 memset(d, 0, top_nrows * row_size);
4567 memset(d + (ny - bottom_nrows) * nx, 0, bottom_nrows * row_size);
4568
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));
4572 }
4573 image->update();
4574}
4575
4577{
4578 if (!image) {
4579 LOGWARN("NULL Image");
4580 return;
4581 }
4582
4583 if (image->get_zsize() <= 1) {
4584 LOGERR("ZeroEdgePlaneProcessor only support 3D models");
4585 throw ImageDimensionException("3D model only");
4586 }
4587
4588 int nx = image->get_xsize();
4589 int ny = image->get_ysize();
4590 int nz = image->get_zsize();
4591
4592 float *d = image->get_data();
4593
4594 int x0=params["x0"];
4595 int x1=params.set_default("x1",x0);
4596 int y0=params["y0"];
4597 int y1=params.set_default("y1",y0);
4598 int z0=params["z0"];
4599 int z1=params.set_default("z1",z0);
4600
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;
4606 int max_y = ny-y1;
4607 size_t x0size = x0*sizeof(float);
4608 size_t x1size = x1*sizeof(float);
4609
4610 memset(d,0,z0*sec_size); // zero -z
4611 memset(d+(nxy*(nz-z1)),0,sec_size*z1); // zero +z
4612
4613 for (int z=z0; z<nz-z1; z++) {
4614 memset(d+z*nxy,0,y0row); // zero -y
4615 memset(d+z*nxy+(ny-y1)*nx,0,y1row); // zero +y
4616
4617 int znxy = z * nxy;
4618 int znxy2 = znxy + nx - x1;
4619
4620 for (int y=y0; y<max_y; y++) {
4621 memset(d+znxy+y*nx,0,x0size); // zero -x
4622 memset(d+znxy2+y*nx,0,x1size); // zero +x
4623 }
4624 }
4625
4626 image->update();
4627}
4628
4629
4631{
4632 return image->get_attr("sigma");
4633}
4634
4636{
4637 if (!image) {
4638 LOGWARN("cannot do normalization on NULL image");
4639 return;
4640 }
4641
4642 if (image->is_complex()) {
4643 LOGWARN("cannot do normalization on complex image");
4644 return;
4645 }
4646
4647 float sigma = calc_sigma(image);
4648 if (sigma == 0 || !Util::goodf(&sigma)) {
4649 LOGWARN("cannot do normalization on image with sigma = 0");
4650 return;
4651 }
4652
4653 float mean = calc_mean(image);
4654
4655 size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
4656 float *data = image->get_data();
4657
4658 for (size_t i = 0; i < size; ++i) {
4659 data[i] = (data[i] - mean) / sigma;
4660 }
4661
4662 image->update();
4663}
4664
4666{
4667 if (!image) {
4668 LOGWARN("NULL Image");
4669 return 0;
4670 }
4671 float ret=sqrt((float)image->get_attr("square_sum"));
4672 return ret==0.0f?1.0f:ret;
4673}
4674
4676{
4677 if (!image) {
4678 LOGWARN("NULL Image");
4679 return 0;
4680 }
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;
4683}
4684
4686{
4687 if (!image) {
4688 LOGWARN("NULL Image");
4689 return;
4690 }
4691 EMData *mask = params["mask"];
4692 int no_sigma = params.set_default("no_sigma",0);
4693 int apply_mask = params.set_default("apply_mask",0);
4694
4695 int nx=image->get_xsize();
4696 int ny=image->get_ysize();
4697 int nz=image->get_zsize();
4698
4699 double sum=0;
4700 double sumsq=0;
4701 double nnz=0;
4702
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);
4707 if (m==0) continue;
4708 if (!apply_mask) m=1.0; // If not applying the mask (which may have values between 0 and 1) we don't want to alter the values
4709 double v=(double)image->get_value_at(x,y,z)*m;
4710 sum+=v;
4711 sumsq+=v*v;
4712 nnz+=1.0;
4713 }
4714 }
4715 }
4716
4717 float mean=(float)(sum/nnz);
4718 float sigma=sqrt((float)(sumsq-sum*sum/nnz)/nnz);
4719 if (no_sigma) sigma=1.0f;
4720
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; // If not applying the mask (which may have values between 0 and 1) we don't want to alter the values
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);
4728 }
4729 }
4730 }
4731}
4732
4734{
4735 if (!image) {
4736 LOGWARN("cannot do normalization on NULL image");
4737 return;
4738 }
4739
4740 if (image->is_complex()) {
4741 LOGWARN("cannot do normalization on complex image");
4742 return;
4743 }
4744
4745 image->process_inplace( "filter.ramp" );
4746 int nx = image->get_xsize();
4747 EMData mask(nx,nx);
4748 mask.process_inplace("testimage.circlesphere", Dict("radius",nx/2-2,"fill",1));
4749
4750 vector<float> rstls = Util::infomask( image, &mask, false);
4751 image->add((float)-rstls[0]);
4752 image->mult((float)1.0/rstls[1]);
4753 image->update();
4754}
4755
4757{
4758 float mass = params.set_default("mass",-1.0f);
4759 int verbose = params.set_default("verbose",0);
4760
4761 if (mass <= 0) throw InvalidParameterException("You must specify a positive non zero mass");
4762
4763 float tthr = params.set_default("thr",(float)image->get_attr("mean")+(float)image->get_attr("sigma"));
4764
4765 float apix = image->get_attr_default("apix_x",1.0f);
4766 apix = params.set_default("apix",apix);
4767
4768 if (apix <= 0) throw InvalidParameterException("You must specify a positive non zero apix");
4769
4770 float step = ((float)image->get_attr("sigma"))/5.0f;
4771
4772 if (step==0) throw InvalidParameterException("This image has sigma=0, cannot give it mass");
4773
4774
4775 size_t n = image->get_size();
4776 float* d = image->get_data();
4777
4778
4779 float thr=(float)image->get_attr("mean")+(float)image->get_attr("sigma")/2.0;
4780 int count=0;
4781 for (size_t i=0; i<n; ++i) {
4782 if (d[i]>=thr) ++count;
4783 }
4784 if (verbose) printf("apix=%1.3f\tmass=%1.1f\tthr=%1.2f\tstep=%1.3g\n",apix,mass,thr,step);
4785
4786 float max = image->get_attr("maximum");
4787 float min = image->get_attr("minimum");
4788 for (int j=0; j<4; j++) {
4789 int err=0;
4790 while (thr<max && count*apix*apix*apix*.81/1000.0>mass) {
4791 thr+=step;
4792 count=0;
4793 for (size_t i=0; i<n; ++i) {
4794 if (d[i]>=thr) ++count;
4795 }
4796 err+=1;
4797 if (err>1000) throw InvalidParameterException("Specified mass could not be achieved");
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);
4799 }
4800
4801 step/=4.0;
4802
4803 while (thr>min && count*apix*apix*apix*.81/1000.0<mass) {
4804 thr-=step;
4805 count=0;
4806 for (size_t i=0; i<n; ++i) {
4807 if (d[i]>=thr) ++count;
4808 }
4809 err+=1;
4810 if (err>1000) throw InvalidParameterException("Specified mass could not be achieved");
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);
4812
4813 }
4814 step/=4.0;
4815 }
4816
4817 // We don't adjust the map if we get a negative value
4818 if ((float)tthr/thr>0) image->mult((float)tthr/thr);
4819 else printf("WARNING: could not normalize map to specified mass.");
4820 image->update();
4821}
4822
4824{
4825 if (!image) {
4826 LOGWARN("NULL Image");
4827 return 0;
4828 }
4829 return image->get_edge_mean();
4830}
4831
4833{
4834 if (!image) {
4835 LOGWARN("NULL Image");
4836 return 0;
4837 }
4838// return image->get_circle_mean();
4839 int nx=image->get_xsize();
4840 int ny=image->get_ysize();
4841 int nz=image->get_zsize();
4842
4843 float radius = params.set_default("radius",((float)ny/2-2));
4844 float width = params.set_default("width",2.0f);
4845 if (radius<0) radius=ny/2+radius;
4846
4847 static bool busy = false; // reduce problems with threads and different image sizes
4848 static EMData *mask = 0;
4849 static int oldradius=radius;
4850 static int oldwidth=width;
4851
4852 if (!mask || !EMUtil::is_same_size(image, mask)||radius!=oldradius||width!=oldwidth) {
4853 while (busy) ;
4854 busy=true;
4855 if (!mask) {
4856 mask = new EMData();
4857 }
4858 mask->set_size(nx, ny, nz);
4859 mask->to_one();
4860
4861 mask->process_inplace("mask.sharp", Dict("inner_radius", radius,
4862 "outer_radius", radius + width));
4863
4864 }
4865 busy=true;
4866 double n = 0,s=0;
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]; }
4872 }
4873
4874 float result = (float)(s/n);
4875// printf("cmean=%f\n",result);
4876 busy=false;
4877
4878 return result;
4879}
4880
4881
4883{
4884 if (!image) {
4885 LOGWARN("NULL Image");
4886 return 0;
4887 }
4888 float maxval = image->get_attr("maximum");
4889 float minval = image->get_attr("minimum");
4890 return (maxval + minval) / 2;
4891}
4892
4894{
4895 if (!image) {
4896 LOGWARN("NULL Image");
4897 return 0;
4898 }
4899 float maxval = image->get_attr("maximum");
4900 float minval = image->get_attr("minimum");
4901 return (maxval - minval) / 2;
4902}
4903
4905{
4906 if (!image) {
4907 LOGWARN("NULL Image");
4908 return 0;
4909 }
4910 double sum = 0;
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;
4916
4917 for (size_t i = 0; i < nyz; i++) {
4918 size_t l = i * nx;
4919 size_t r = l + nx - 2;
4920 sum += d[l] + d[l + 1] + d[r] + d[r + 1];
4921 }
4922 float mean = (float) sum / (4 * nyz);
4923 return mean;
4924}
4925
4927{
4928 if (!image) {
4929 LOGWARN("NULL Image");
4930 return;
4931 }
4932
4933
4934
4935
4936 float *rdata = image->get_data();
4937 int nx = image->get_xsize();
4938 int ny = image->get_ysize();
4939 int nz = image->get_zsize();
4940
4941 if (nz > 1) {
4942 // Muyuan 08/2018: normalize by slice when nz>1
4943 // since this processor is not used anyways..
4944 for (int z = 0; z < nz; z++) {
4945 double row_sum = 0;
4946 double row_sqr = 0;
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];
4950 row_sum += d;
4951 d*=d;
4952 row_sqr += d;
4953 }
4954 }
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;
4962 }
4963 }
4964
4965 }
4966
4967
4968// LOGERR("row normalize only works for 2D image");
4969// return;
4970 }
4971
4972 else {
4973 if (params.has_key("unitlen")) {
4974 for (int y = 0; y < ny; y++) {
4975 double row_len = 0;
4976 for (int x = 0; x < nx; x++) {
4977 row_len += pow((double)rdata[x + y * nx],2.0);
4978 }
4979 row_len=sqrt(row_len);
4980 if (row_len==0) row_len=1.0;
4981
4982 for (int x = 0; x < nx; x++) {
4983 rdata[x + y * nx] /= (float)row_len;
4984 }
4985 }
4986 image->update();
4987 return;
4988 }
4989
4990 for (int y = 0; y < ny; y++) {
4991 double row_sum = 0;
4992 for (int x = 0; x < nx; x++) {
4993 row_sum += rdata[x + y * nx];
4994 }
4995
4996 double row_mean = row_sum / nx;
4997 if (row_mean <= 0) {
4998 row_mean = 1;
4999 }
5000
5001 for (int x = 0; x < nx; x++) {
5002 rdata[x + y * nx] /= (float)row_mean;
5003 }
5004 }
5005 }
5006 image->update();
5007}
5008
5010{
5011 if (!image) {
5012 LOGWARN("NULL Image");
5013 return 0;
5014 }
5015 return image->get_attr("mean");
5016}
5017
5019{
5020 if (!image) {
5021 LOGWARN("NULL Image");
5022 return 0;
5023 }
5024 float mean = image->get_attr("mean");
5025 float sig = image->get_attr("sigma");
5026 size_t n=image->get_size();
5027
5028 // We want at least 100 values per bin on average
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);
5033
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);
5040 hist[bin]++;
5041 }
5042
5043 // yes, probably could use fancy vector methods...
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; }
5047 }
5048
5049 // weighted average position of the peak and nearest neighbors
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);
5053
5054 return ret;
5055}
5056
5057
5058
5060{
5061 if (!image) {
5062 LOGWARN("NULL Image");
5063 return NULL;
5064 }
5065
5066 EMData *refr = params["ref"];
5067 EMData *actual = params.set_default("actual",(EMData*)NULL);
5068 EMData *ref;
5069 bool return_radial = params.set_default("return_radial",false);
5070 bool return_fft = params.set_default("return_fft",false);
5071 bool ctfweight = params.set_default("ctfweight",false);
5072 bool return_presigma = params.set_default("return_presigma",false);
5073 bool return_subim = params.set_default("return_subim",false);
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()); // include the corners unless explicitly excluded
5076
5077 // We will be modifying imf, so it needs to be a copy
5078 EMData *imf;
5079 if (image->is_complex()) imf=image->copy();
5080 else imf=image->do_fft();
5081
5082 if (ctfweight) {
5083 EMData *ctfi=imf->copy_head();
5084 Ctf *ctf;
5085// if (image->has_attr("ctf"))
5086 ctf=(Ctf *)(image->get_attr("ctf"));
5087// else ctf=(Ctf *)(ref->get_attr("ctf"));
5089 imf->mult(*ctfi);
5090 delete ctfi;
5091 }
5092
5093 // Make sure ref is complex
5094 if (refr->is_complex()) ref=refr;
5095 else ref=refr->do_fft();
5096
5097 EMData *actf;
5098 if (actual==NULL) actf=ref;
5099 else {
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();
5103 }
5104
5105 int ny2=(int)(image->get_ysize()*sqrt(2.0)/2);
5106 vector <double>rad(ny2+1);
5107 vector <double>norm(ny2+1);
5108
5109 // We are essentially computing an FSC here, but while the reference (the image
5110 // we plan to subtract) is normalized, the other image is not. This gives us a filter
5111 // to apply to the reference to optimally eliminate its contents from 'image'.
5112 for (int y=-ny2; y<ny2; y++) {
5113 for (int x=0; x<ny2; x++) {
5114 int r=int(Util::hypot_fast(x,y));
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());
5119// norm[r]+=v2.real()*v2.real()+v2.imag()*v2.imag()+v1.real()*v1.real()+v1.imag()*v1.imag();
5120 norm[r]+=(double)(v2.real()*v2.real()+v2.imag()*v2.imag());
5121 }
5122 }
5123 for (int i=1; i<ny2; i++) rad[i]/=norm[i];
5124 rad[0]=0;
5125// FILE *out=fopen("dbug.txt","w");
5126// for (int i=0; i<ny2; i++) fprintf(out,"%lf\t%lf\t%lf\n",(float)i,rad[i],norm[i]);
5127// fclose(out);
5128
5129 float oldsig=-1.0;
5130 // This option computes the real-space sigma on the input-image after the specified filter
5131 // This is an expensive option, but more efficient than computing the same using other means
5132 if (return_presigma) {
5133 for (int y=-ny2; y<ny2; y++) {
5134 for (int x=0; x<imf->get_xsize()/2; x++) {
5135 int r=int(Util::hypot_fast(x,y));
5136 if (r>=ny2 || r>=si1 || r<si0) {
5137 imf->set_complex_at(x,y,0);
5138 continue;
5139 }
5140 std::complex<float> v1=imf->get_complex_at(x,y);
5141 imf->set_complex_at(x,y,v1);
5142 }
5143 }
5144 EMData *tmp=imf->do_ift();
5145 oldsig=(float)tmp->get_attr("sigma");
5146 delete tmp;
5147 }
5148
5149 for (int y=-ny2; y<ny2; y++) {
5150 for (int x=0; x<imf->get_xsize()/2; x++) {
5151 int r=int(Util::hypot_fast(x,y));
5152 if (r>=ny2 || r>=si1 || r<si0) {
5153 imf->set_complex_at(x,y,0);
5154 continue;
5155 }
5156 std::complex<float> v1=imf->get_complex_at(x,y);
5157 std::complex<float> v2=actf->get_complex_at(x,y);
5158 v2*=(float)rad[r];
5159 if (return_subim) imf->set_complex_at(x,y,v2);
5160 else imf->set_complex_at(x,y,v1-v2);
5161 }
5162 }
5163
5164 if (!refr->is_complex()) delete ref;
5165 if (actual!=NULL && !actual->is_complex()) delete actf;
5166
5167 vector <float>radf;
5168 if (return_radial) {
5169 radf.resize(ny2);
5170 for (int i=0; i<ny2; i++) radf[i]=(float)rad[i];
5171 }
5172
5173 if (!return_fft) {
5174 EMData *ret=imf->do_ift();
5175 delete imf;
5176 if (return_radial) ret->set_attr("filter_curve",radf);
5177 if (return_presigma) {
5178 ret->set_attr("sigma_presub",oldsig);
5179// printf("Set %f\n",(float)image->get_attr("sigma_presub"));
5180 }
5181 return ret;
5182 }
5183 if (return_radial) imf->set_attr("filter_curve",radf);
5184 if (return_presigma) imf->set_attr("sigma_presub",oldsig);
5185 return imf;
5186}
5187
5189{
5190 if (!image) {
5191 LOGWARN("NULL image");
5192 return;
5193 }
5194
5195 EMData *tmp=process(image);
5196 memcpy(image->get_data(),tmp->get_data(),(size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));
5197 delete tmp;
5198 image->update();
5199 return;
5200}
5201
5202
5204{
5205 if (!image) {
5206 LOGWARN("NULL Image");
5207 return;
5208 }
5209
5210 EMData *to = params["to"];
5211
5212 bool ignore_zero = params.set_default("ignore_zero",true);
5213 bool fourieramp = params.set_default("fourieramp",false);
5214 bool debug = params.set_default("debug",false);
5215 float ignore_lowsig = params.set_default("ignore_lowsig",-1.0);
5216 float low_threshold = params.set_default("low_threshold",-FLT_MAX);
5217 float high_threshold = params.set_default("high_threshold",FLT_MAX);
5218
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;
5223 size_t size2 = 0;
5224
5225 EMData *fim = NULL;
5226 EMData *fto = NULL;
5227 float *dimage; // unthresholded image data
5228 float *dto;
5229 int step=1; // used in Fourier mode
5230
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;
5235
5236 if (fourieramp) {
5237 fim=image->do_fft();
5238 fto=to->do_fft();
5239 fim->ri2ap();
5240 fto->ri2ap();
5241 dimage=fim->get_data();
5242 dto=fto->get_data();
5243 step=2;
5244 size2=(size_t)(nx+2) * ny * nz;
5245
5246 // sigma for thresholding
5247 meani=meant=0; // Fourier amplitude >=0, so we reference everything to this point
5248 sigi=sigt=0;
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));
5252 }
5253 else {
5254 dimage = image->get_data();
5255 dto = to->get_data();
5256 size2=size;
5257 }
5258
5259
5260 // rewrote this to just use GSL and get rid of David's old code.
5261 // The two passes are just to make sure we don't eat too much RAM if we do 3D
5262 if (ignore_lowsig<0) ignore_lowsig=0;
5263
5264 size_t count=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))) {
5270 count++;
5271 }
5272 }
5273
5274
5275 double *x=(double *)malloc(count*sizeof(double));
5276 double *y=(double *)malloc(count*sizeof(double));
5277 count=0;
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))) {
5283 x[count]=dimage[i];
5284 y[count]=dto[i];
5285 count++;
5286 }
5287 }
5288 double c0,c1;
5289 double cov00,cov01,cov11,sumsq;
5290 gsl_fit_linear (x, 1, y, 1, count, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
5291
5292 if (debug) {
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]);
5296 }
5297 fclose(out);
5298 printf("add %lf\tmul %lf\t%lf\t%lf\t%lf\t%lf\n",c0,c1,cov00,cov01,cov11,sumsq);
5299 }
5300
5301 free(x);
5302 free(y);
5303 if (fim!=NULL) delete fim;
5304 if (fto!=NULL) delete fto;
5305
5306 if (fourieramp) c0=0; // c0 makes no sense in this context. Really just a measure of noise
5307 dimage = image->get_data();
5308 for (size_t i = 0; i < size; ++i) dimage[i]=dimage[i]*c1+c0;
5309 image->update();
5310 image->set_attr("norm_mult",c1);
5311 image->set_attr("norm_add",c0);
5312}
5313
5315 ENTERFUNC;
5316 if (!image->is_complex()) throw ImageFormatException("Fourier binary thresholding processor only works for complex images");
5317
5318 float threshold = params.set_default("value",0.0f);
5319 image->ri2ap(); // works for cuda
5320
5321 float* d = image->get_data();
5322 for( size_t i = 0; i < image->get_size()/2; ++i, d+=2) {
5323 if ( *d < threshold ) {
5324 *d = 0;
5325 *(d+1) = 0;
5326 }
5327 }
5328 image->ap2ri();
5329 image->set_ri(true); // So it can be used for fourier multiplaction, for example
5330 image->update();
5331 EXITFUNC;
5332}
5333
5335{
5336 if (!image) {
5337 LOGWARN("NULL Image");
5338 return;
5339 }
5340
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"];
5345
5346 if (half_width < distance_sigma) {
5347 LOGWARN("localwidth(=%d) should be larger than distance_sigma=(%f)\n",
5348 half_width, distance_sigma);
5349 }
5350
5351 distance_sigma *= distance_sigma;
5352
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);
5357 }
5358 value_sigma *= value_sigma;
5359
5360 int nx = image->get_xsize();
5361 int ny = image->get_ysize();
5362 int nz = image->get_zsize();
5363
5364 if(nz==1) { //for 2D image
5365 int width=nx, height=ny;
5366
5367 int i,j,m,n;
5368
5369 float tempfloat1,tempfloat2,tempfloat3;
5370 int index1,index2,index;
5371 int Iter;
5372 int tempint1,tempint3;
5373
5374 tempint1=width;
5375 tempint3=width+2*half_width;
5376
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(