doxygen/processor_8cpp_source.html

/*

 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)

 * Copyright (c) 2000-2006 Baylor College of Medicine

 *

 * This software is issued under a joint BSD/GNU license. You may use the

 * source code in this file under either license. However, note that the

 * complete EMAN2 and SPARX software packages have some GPL dependencies,

 * so you are responsible for compliance with the licenses of these packages

 * if you opt to use BSD licensing. The warranty disclaimer below holds

 * in either instance.

 *

 * This complete copyright notice must be included in any revised version of the

 * source code. Additional authorship citations may be added, but existing

 * author citations must be preserved.

 *

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 2 of the License, or

 * (at your option) any later version.

 *

 * This program is distributed in the hope that it will be useful,

 * but WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

 * GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

 *

 * */


#include "processor.h"

#include "sparx/processor_sparx.h"

#include "plugins/processor_template.h"

#include "reconstructor_tools.h"

#include "cmp.h"

#include "ctf.h"

#include "xydata.h"

#include "emdata.h"

#include "emfft.h"

#include "emassert.h"

#include "randnum.h"

#include "symmetry.h"

#include "averager.h"

#include "util.h"


#include <gsl/gsl_randist.h>

#include <gsl/gsl_statistics.h>

#include <gsl/gsl_wavelet.h>

#include <gsl/gsl_wavelet2d.h>

#include <gsl/gsl_multimin.h>

#include <algorithm>

#include <gsl/gsl_fit.h>

#include <ctime>


#ifdef __APPLE__

        typedef unsigned int uint;

#endif  //__APPLE__


#ifdef _WIN32

        typedef unsigned int uint;

#endif  //_WIN32


#ifdef EMAN2_USING_CUDA

//#include "cuda/cuda_util.h"

#include "cuda/cuda_processor.h"

#endif // EMAN2_USING_CUDA


using namespace EMAN;

using std::reverse;


const string SNREvalProcessor::NAME = "eval.maskedsnr";

const string AmpweightFourierProcessor::NAME = "filter.ampweight";

const string Axis0FourierProcessor::NAME = "filter.xyaxes0";

const string ConvolutionProcessor::NAME = "math.convolution";

const string BispecSliceProcessor::NAME = "math.bispectrum.slice";

const string HarmonicProcessor::NAME = "math.harmonic";

const string XGradientProcessor::NAME = "math.edge.xgradient";

const string YGradientProcessor::NAME = "math.edge.ygradient";

const string ZGradientProcessor::NAME = "math.edge.zgradient";

const string ImageDivergenceProcessor::NAME = "math.divergence";

const string GradientMagnitudeProcessor::NAME = "math.gradient.magnitude";

const string GradientDirectionProcessor::NAME = "math.gradient.direction";

const string SDGDProcessor::NAME = "math.laplacian.sdgd";

const string LaplacianProcessor::NAME = "math.laplacian";

const string LaplacianDirectionProcessor::NAME = "math.laplacian.direction";

const string LaplacianMagnitudeProcessor::NAME = "math.laplacian.magnitude";

const string Wiener2DAutoAreaProcessor::NAME = "filter.wiener2dauto";

const string Wiener2DFourierProcessor::NAME = "filter.wiener2d";

const string CtfSimProcessor::NAME = "math.simulatectf";

const string LinearRampFourierProcessor::NAME = "filter.linearfourier";

const string LoGFourierProcessor::NAME = "filter.LoG";

const string DoGFourierProcessor::NAME = "filter.DoG";

const string AzSharpProcessor::NAME = "filter.azimuthal.contrast";

const string HighpassAutoPeakProcessor::NAME = "filter.highpass.autopeak";

const string LinearRampProcessor::NAME = "eman1.filter.ramp";

const string AbsoluteValueProcessor::NAME = "math.absvalue";

const string CCCSNRProcessor::NAME = "math.ccc_snr_wiener";

const string FloorValueProcessor::NAME = "math.floor";

const string BooleanProcessor::NAME = "threshold.notzero";

const string KmeansSegmentProcessor::NAME = "segment.kmeans";

const string GaussSegmentProcessor::NAME = "segment.gauss";

const string DistanceSegmentProcessor::NAME = "segment.distance";

const string WatershedProcessor::NAME = "segment.watershed";

const string SegmentSubunitProcessor::NAME = "segment.subunit";

const string RecipCarefullyProcessor::NAME = "math.reciprocal";

const string SubtractOptProcessor::NAME = "math.sub.optimal";

const string ValuePowProcessor::NAME = "math.pow";

const string ValueSquaredProcessor::NAME = "math.squared";

const string ValueSqrtProcessor::NAME = "math.sqrt";

const string DiscritizeProcessor::NAME = "threshold.discritize.sigma";

const string ToZeroProcessor::NAME = "threshold.belowtozero";

const string AboveToZeroProcessor::NAME="threshold.abovetozero";

const string OutlierProcessor::NAME="threshold.outlier.localmean";

const string Rotate180Processor::NAME = "math.rotate.180";

const string TransformProcessor::NAME = "xform";

const string IntTranslateProcessor::NAME = "xform.translate.int";

const string ScaleTransformProcessor::NAME = "xform.scale";

const string ApplySymProcessor::NAME = "xform.applysym";

const string ClampingProcessor::NAME = "threshold.clampminmax";

const string RangeZeroProcessor::NAME = "threshold.rangetozero";

const string NSigmaClampingProcessor::NAME = "threshold.clampminmax.nsigma";

const string ToMinvalProcessor::NAME = "threshold.belowtominval";

const string CutToZeroProcessor::NAME = "threshold.belowtozero_cut";

const string BinarizeProcessor::NAME = "threshold.binary";

//const string BinarizeAmpProcessor::NAME = "threshold.amp.binary";

const string BinarizeFourierProcessor::NAME = "threshold.binary.fourier";

const string CollapseProcessor::NAME = "threshold.compress";

const string LinearXformProcessor::NAME = "math.linear";

const string ExpProcessor::NAME = "math.exp";

const string FiniteProcessor::NAME = "math.finite";

const string RangeThresholdProcessor::NAME = "threshold.binaryrange";

const string SigmaProcessor::NAME = "math.sigma";

const string LogProcessor::NAME = "math.log";

const string MaskSharpProcessor::NAME = "mask.sharp";

const string MaskSoftProcessor::NAME = "mask.soft";

const string MaskEdgeMeanProcessor::NAME = "mask.ringmean";

const string MaskNoiseProcessor::NAME = "mask.noise";

const string MaskGaussProcessor::NAME = "mask.gaussian";

const string MaskAzProcessor::NAME = "mask.cylinder";

const string MaskGaussNonuniformProcessor::NAME = "mask.gaussian.nonuniform";

const string MaskGaussInvProcessor::NAME = "math.gausskernelfix";

const string GridKernelFixProcessor::NAME = "math.gridkernelfix";

const string LinearPyramidProcessor::NAME = "math.linearpyramid";

const string SetBitsProcessor::NAME = "math.setbits";

const string MakeRadiusSquaredProcessor::NAME = "math.toradiussqr";

const string MakeRadiusProcessor::NAME = "math.toradius";

const string ComplexNormPixel::NAME = "complex.normpixels";

const string ZeroConstantProcessor::NAME = "mask.contract"; // This is broken, it never worked. Somebody didn't think it through properly

const string BoxMedianProcessor::NAME = "eman1.filter.median";

const string BoxSigmaProcessor::NAME = "math.localsigma";

const string NonConvexProcessor::NAME = "math.nonconvex";

const string BoxMaxProcessor::NAME = "math.localmax";

const string LocalMinAbsProcessor::NAME = "math.localminabs";

const string MinusPeakProcessor::NAME = "math.submax";

const string PeakOnlyProcessor::NAME = "mask.onlypeaks";

const string DiffBlockProcessor::NAME = "eman1.filter.blockrange";

const string CutoffBlockProcessor::NAME = "eman1.filter.blockcutoff";

const string MaxShrinkProcessor::NAME = "math.maxshrink";

const string MinShrinkProcessor::NAME = "math.minshrink";

const string MeanShrinkProcessor::NAME = "math.meanshrink";

const string MedianShrinkProcessor::NAME = "math.medianshrink";

const string FFTResampleProcessor::NAME = "math.fft.resample";

const string GradientRemoverProcessor::NAME = "math.lineargradientfix";

const string GradientPlaneRemoverProcessor::NAME = "filter.gradientPlaneRemover";

const string FlattenBackgroundProcessor::NAME = "filter.flattenbackground";

const string RampProcessor::NAME = "filter.ramp";

const string VerticalStripeProcessor::NAME = "math.verticalstripefix";

const string RealToFFTProcessor::NAME = "math.realtofft";

const string SigmaZeroEdgeProcessor::NAME = "mask.zeroedgefill";

const string WedgeFillProcessor::NAME = "mask.wedgefill";

const string FFTPeakProcessor::NAME = "mask.fft.peak";

const string FFTConeProcessor::NAME = "mask.fft.cone";

const string FFTWedgeProcessor::NAME = "mask.fft.wedge";

const string BeamstopProcessor::NAME = "mask.beamstop";

const string MeanZeroEdgeProcessor::NAME = "mask.dampedzeroedgefill";

const string AverageXProcessor::NAME = "math.averageovery";

const string DecayEdgeProcessor::NAME = "mask.decayedge2d";

const string ZeroEdgeRowProcessor::NAME = "mask.zeroedge2d";

const string ZeroEdgePlaneProcessor::NAME = "mask.zeroedge3d";

const string BilateralProcessor::NAME = "filter.bilateral";

const string NormalizeUnitProcessor::NAME = "normalize.unitlen";

const string NormalizeUnitSumProcessor::NAME = "normalize.unitsum";

const string NormalizeStdProcessor::NAME = "normalize";

const string NormalizeMaskProcessor::NAME = "normalize.mask";

const string NormalizeRampNormVar::NAME = "normalize.ramp.normvar";

const string NormalizeByMassProcessor::NAME = "normalize.bymass";

const string NormalizeEdgeMeanProcessor::NAME = "normalize.edgemean";

const string NormalizeCircleMeanProcessor::NAME = "normalize.circlemean";

const string NormalizeLREdgeMeanProcessor::NAME = "normalize.lredge";

const string NormalizeHistPeakProcessor::NAME = "normalize.histpeak";

const string NormalizeMaxMinProcessor::NAME = "normalize.maxmin";

const string NormalizeRowProcessor::NAME = "normalize.rows";

const string NormalizeToLeastSquareProcessor::NAME = "normalize.toimage";

const string RotationalAverageProcessor::NAME = "math.rotationalaverage";

const string RotationalSubstractProcessor::NAME = "math.rotationalsubtract";

const string TransposeProcessor::NAME = "xform.transpose";

const string FlipProcessor::NAME = "xform.flip";

const string AddNoiseProcessor::NAME = "math.addnoise";

const string AddSigmaNoiseProcessor::NAME = "math.addsignoise";

const string AddRandomNoiseProcessor::NAME = "math.addspectralnoise";

const string FourierToCornerProcessor::NAME = "xform.fourierorigin.tocorner";

const string FourierToCenterProcessor::NAME = "xform.fourierorigin.tocenter";

const string PhaseToCenterProcessor::NAME = "xform.phaseorigin.tocenter";

const string PhaseToCornerProcessor::NAME = "xform.phaseorigin.tocorner";

const string AutoMask2DProcessor::NAME = "mask.auto2d";

const string AutoMaskAsymUnit::NAME = "mask.asymunit";

const string AutoMask3DProcessor::NAME = "mask.auto3d.thresh";

const string AutoMask3D2Processor::NAME = "mask.auto3d";

const string AutoMaskDustProcessor::NAME = "mask.dust3d";

const string AddMaskShellProcessor::NAME = "mask.addshells";

const string IterMultiMaskProcessor::NAME = "mask.addshells.multilevel";

const string IterBinMaskProcessor::NAME = "mask.addshells.gauss";

const string PhaseToMassCenterProcessor::NAME = "xform.phasecenterofmass";

const string ToMassCenterProcessor::NAME = "xform.centerofmass";

const string ToCenterProcessor::NAME = "xform.center";

const string ACFCenterProcessor::NAME = "xform.centeracf";

const string CTFCorrProcessor::NAME = "filter.ctfcorr.simple";

const string SNRProcessor::NAME = "eman1.filter.snr";

const string FileFourierProcessor::NAME = "eman1.filter.byfile";

const string FSCFourierProcessor::NAME = "filter.wiener.byfsc";

const string SymSearchProcessor::NAME = "misc.symsearch";

const string MaskPackProcessor::NAME = "misc.mask.pack";

const string LocalNormProcessor::NAME = "normalize.local";

const string StripeXYProcessor::NAME = "math.xystripefix";

const string BadLineXYProcessor::NAME = "math.xybadlines";

const string IndexMaskFileProcessor::NAME = "mask.fromfile";

// const string CoordinateMaskFileProcessor::NAME = "mask.fromfile.sizediff";

const string PaintProcessor::NAME = "mask.paint";

const string DirectionalSumProcessor::NAME = "misc.directional_sum";

template<> const string BinaryOperateProcessor<MaxPixelOperator>::NAME = "math.max";            // These 2 should not really be processors

template<> const string BinaryOperateProcessor<MinPixelOperator>::NAME = "math.min";

const string MaxPixelOperator::NAME = "math.max";

const string MinPixelOperator::NAME = "math.min";

const string MatchSFProcessor::NAME = "filter.matchto";

const string SetSFProcessor::NAME = "filter.setstrucfac";

const string SetIsoPowProcessor::NAME = "filter.setisotropicpow";

const string SmartMaskProcessor::NAME = "mask.smart";

const string TestImagePureGaussian::NAME = "testimage.puregaussian";

const string TestImageFourierGaussianBand::NAME = "testimage.fourier.gaussianband";

const string TestImageFourierNoiseGaussian::NAME = "testimage.noise.fourier.gaussian";

const string TestImageFourierNoiseProfile::NAME = "testimage.noise.fourier.profile";

const string CTFSNRWeightProcessor::NAME = "ctf.snr.weight";

const string TestImageLineWave::NAME = "testimage.linewave";

const string TestTomoImage::NAME = "testimage.tomo.objects";

const string TestImageGradient::NAME = "testimage.gradient";

const string TestImageAxes::NAME = "testimage.axes";

const string TestImageGaussian::NAME = "testimage.gaussian";

const string TestImageScurve::NAME = "testimage.scurve";

const string TestImageSphericalWave::NAME = "testimage.sphericalwave";

const string TestImageSinewave::NAME = "testimage.sinewave";

const string TestImageSinewaveCircular::NAME = "testimage.sinewave.circular";

const string TestImageSquarecube::NAME = "testimage.squarecube";

const string TestImageEllipse::NAME = "testimage.ellipsoid";

const string TestImageHollowEllipse::NAME = "testimage.ellipsoid.hollow";

const string TestImageCirclesphere::NAME = "testimage.circlesphere";

const string TestImageNoiseUniformRand::NAME = "testimage.noise.uniform.rand";

const string TestImageNoiseGauss::NAME = "testimage.noise.gauss";

const string TestImageCylinder::NAME = "testimage.cylinder";

const string TestImageDisc::NAME = "testimage.disc";

const string CCDNormProcessor::NAME = "filter.ccdnorm";

const string WaveletProcessor::NAME = "basis.wavelet";

const string EnhanceProcessor::NAME = "filter.enhance";

const string TomoTiltEdgeMaskProcessor::NAME = "tomo.tiltedgemask";

const string TomoTiltAngleWeightProcessor::NAME = "tomo.tiltangleweight";

const string FFTProcessor::NAME = "basis.fft";

const string RadialProcessor::NAME = "mask.radialprofile";

const string HistogramBin::NAME = "histogram.bin";

const string ModelEMCylinderProcessor::NAME = "math.model_em_cylinder";

const string ApplyPolynomialProfileToHelix::NAME = "math.poly_radial_profile";

const string BinarySkeletonizerProcessor::NAME="gorgon.binary_skel";

const string MirrorProcessor::NAME = "xform.mirror";

const string ReverseProcessor::NAME = "xform.reverse";

const string NewLowpassTopHatProcessor::NAME = "filter.lowpass.tophat";

const string NewHighpassTopHatProcessor::NAME = "filter.highpass.tophat";

const string NewBandpassTopHatProcessor::NAME = "filter.bandpass.tophat";

const string NewHomomorphicTopHatProcessor::NAME = "filter.homomorphic.tophat";

const string NewLowpassGaussProcessor::NAME = "filter.lowpass.gauss";

const string LowpassAutoBProcessor::NAME="filter.lowpass.autob";

const string NewHighpassGaussProcessor::NAME = "filter.highpass.gauss";

const string NewBandpassGaussProcessor::NAME = "filter.bandpass.gauss";

const string NewHomomorphicGaussProcessor::NAME = "filter.homomorphic.gauss";

const string NewInverseGaussProcessor::NAME = "filter.gaussinverse";

const string GaussZFourierProcessor::NAME = "filter.lowpass.gaussz";

const string SHIFTProcessor::NAME = "filter.shift";

const string InverseKaiserI0Processor::NAME = "filter.kaiser_io_inverse";

const string InverseKaiserSinhProcessor::NAME = "filter.kaisersinhinverse";

const string NewRadialTableProcessor::NAME = "filter.radialtable";

const string LowpassRandomPhaseProcessor::NAME = "filter.lowpass.randomphase";

const string NewLowpassButterworthProcessor::NAME = "filter.lowpass.butterworth";

const string NewHighpassButterworthProcessor::NAME = "filter.highpass.butterworth";

const string NewHomomorphicButterworthProcessor::NAME = "filter.homomorphic.butterworth";

const string NewLowpassTanhProcessor::NAME = "filter.lowpass.tanh";

const string NewHighpassTanhProcessor::NAME = "filter.highpass.tanh";

const string NewHomomorphicTanhProcessor::NAME = "filter.homomorphic.tanh";

const string NewBandpassTanhProcessor::NAME = "filter.bandpass.tanh";

const string CTF_Processor::NAME = "filter.CTF_";

const string ConvolutionKernelProcessor::NAME = "filter.convolution.kernel";

const string RotateInFSProcessor::NAME = "rotateinfs";

const string CircularAverageBinarizeProcessor::NAME = "threshold.binary.circularmean";

const string ObjDensityProcessor::NAME = "morph.object.density";

const string ObjLabelProcessor::NAME = "morph.object.label";

const string BwThinningProcessor::NAME = "morph.thin";

const string BwMajorityProcessor::NAME = "morph.majority";

const string PruneSkeletonProcessor::NAME = "morph.prune";

const string ManhattanDistanceProcessor::NAME = "math.distance.manhattan";

const string BinaryDilationProcessor::NAME = "morph.dilate.binary";

const string BinaryErosionProcessor::NAME = "morph.erode.binary";

const string BinaryOpeningProcessor::NAME = "morph.open.binary";

const string BinaryClosingProcessor::NAME = "morph.close.binary";

const string BinaryMorphGradientProcessor::NAME = "morph.gradient.binary";

const string BinaryExternalGradientProcessor::NAME = "morph.ext_grad.binary";

const string BinaryInternalGradientProcessor::NAME = "morph.int_grad.binary";

const string BinaryTopHatProcessor::NAME = "morph.tophat.binary";

const string BinaryBlackHatProcessor::NAME = "morph.blackhat.binary";

const string GrowSkeletonProcessor::NAME = "morph.grow";

const string FixSignProcessor::NAME = "math.fixmode";

const string ZThicknessProcessor::NAME = "misc.zthick";

const string ReplaceValuefromListProcessor::NAME = "misc.colorlabel";

const string PolyMaskProcessor::NAME = "mask.poly";

const string AmpMultProcessor::NAME = "math.multamplitude";


//#ifdef EMAN2_USING_CUDA

//const string CudaMultProcessor::NAME = "cuda.math.mult";

//const string CudaCorrelationProcessor::NAME = "cuda.correlate";

//#endif //EMAN2_USING_CUDA


#if 0

//const string XYZProcessor::NAME = "XYZ";

#endif  //0


template <> Factory < Processor >::Factory()

{

        force_add<HighpassAutoPeakProcessor>();

        force_add<LinearRampFourierProcessor>();

        force_add<LoGFourierProcessor>();

        force_add<DoGFourierProcessor>();

//      force_add<AzSharpProcessor>();

        force_add<FixSignProcessor>();


        force_add<AmpweightFourierProcessor>();

        force_add<Axis0FourierProcessor>();

        force_add<Wiener2DFourierProcessor>();

        force_add<LowpassAutoBProcessor>();

        force_add<CtfSimProcessor>();


        force_add<LinearPyramidProcessor>();

        force_add<LinearRampProcessor>();

        force_add<AbsoluteValueProcessor>();

        force_add<FloorValueProcessor>();

        force_add<BooleanProcessor>();

        force_add<KmeansSegmentProcessor>();

        force_add<GaussSegmentProcessor>();

        force_add<SegmentSubunitProcessor>();

        force_add<DistanceSegmentProcessor>();

        force_add<ValuePowProcessor>();

        force_add<ValueSquaredProcessor>();

        force_add<ValueSqrtProcessor>();

        force_add<DiscritizeProcessor>();

        force_add<Rotate180Processor>();

        force_add<TransformProcessor>();

        force_add<ScaleTransformProcessor>();

        force_add<ApplySymProcessor>();

        force_add<IntTranslateProcessor>();

        force_add<RecipCarefullyProcessor>();

        force_add<SubtractOptProcessor>();


        force_add<ClampingProcessor>();

        force_add<NSigmaClampingProcessor>();


        force_add<ToZeroProcessor>();

        force_add<RangeZeroProcessor>();

        force_add<AboveToZeroProcessor>();

        force_add<OutlierProcessor>();

        force_add<ToMinvalProcessor>();

        force_add<CutToZeroProcessor>();

        force_add<BinarizeProcessor>();

//      force_add<BinarizeAmpProcessor>();

        force_add<BinarizeFourierProcessor>();

        force_add<CollapseProcessor>();

        force_add<LinearXformProcessor>();

        force_add<SetBitsProcessor>();


        force_add<CCCSNRProcessor>();

        force_add<ExpProcessor>();

        force_add<RangeThresholdProcessor>();

        force_add<SigmaProcessor>();

        force_add<LogProcessor>();

        force_add<FiniteProcessor>();


        force_add< BinaryOperateProcessor<MaxPixelOperator> >();

        force_add< BinaryOperateProcessor<MinPixelOperator> >();


        force_add<PaintProcessor>();

        force_add<WatershedProcessor>();

        force_add<MaskSharpProcessor>();

        force_add<MaskSoftProcessor>();

        force_add<MaskEdgeMeanProcessor>();

        force_add<MaskNoiseProcessor>();

        force_add<MaskGaussProcessor>();

        force_add<MaskGaussNonuniformProcessor>();

        force_add<MaskGaussInvProcessor>();

        force_add<GridKernelFixProcessor>();

        force_add<MaskAzProcessor>();


        force_add<MaxShrinkProcessor>();

        force_add<MinShrinkProcessor>();

        force_add<MeanShrinkProcessor>();

        force_add<MedianShrinkProcessor>();

        force_add<FFTResampleProcessor>();

        force_add<NonConvexProcessor>();


        force_add<MakeRadiusSquaredProcessor>();

        force_add<MakeRadiusProcessor>();


        force_add<ComplexNormPixel>();


        force_add<LaplacianProcessor>();

//      force_add<ZeroConstantProcessor>();       // this is badly written, it does not work and never did. Who wrote this !?!?


        force_add<BoxMedianProcessor>();

        force_add<BoxSigmaProcessor>();

        force_add<BoxMaxProcessor>();

        force_add<LocalMinAbsProcessor>();


        force_add<MinusPeakProcessor>();

        force_add<PeakOnlyProcessor>();

        force_add<DiffBlockProcessor>();


        force_add<CutoffBlockProcessor>();

//      force_add<GradientRemoverProcessor>();

        force_add<GradientPlaneRemoverProcessor>();

        force_add<FlattenBackgroundProcessor>();

        force_add<VerticalStripeProcessor>();

        force_add<RealToFFTProcessor>();

        force_add<SigmaZeroEdgeProcessor>();

        force_add<WedgeFillProcessor>();

        force_add<FFTPeakProcessor>();

        force_add<FFTConeProcessor>();

        force_add<FFTWedgeProcessor>();

        force_add<RampProcessor>();


        force_add<BeamstopProcessor>();

        force_add<MeanZeroEdgeProcessor>();

        force_add<AverageXProcessor>();

        force_add<DecayEdgeProcessor>();

        force_add<ZeroEdgeRowProcessor>();

        force_add<ZeroEdgePlaneProcessor>();


        force_add<BilateralProcessor>();


        force_add<ConvolutionProcessor>();

        force_add<BispecSliceProcessor>();

        force_add<HarmonicProcessor>();


        force_add<NormalizeStdProcessor>();

        force_add<NormalizeHistPeakProcessor>();

        force_add<NormalizeUnitProcessor>();

        force_add<NormalizeUnitSumProcessor>();

        force_add<NormalizeMaskProcessor>();

        force_add<NormalizeEdgeMeanProcessor>();

        force_add<NormalizeCircleMeanProcessor>();

        force_add<NormalizeLREdgeMeanProcessor>();

        force_add<NormalizeMaxMinProcessor>();

        force_add<NormalizeByMassProcessor>();

        force_add<NormalizeRowProcessor>();

        force_add<NormalizeRampNormVar>();


        force_add<HistogramBin>();


        force_add<NormalizeToLeastSquareProcessor>();


        force_add<RotationalAverageProcessor>();

        force_add<RotationalSubstractProcessor>();

        force_add<FlipProcessor>();

        force_add<TransposeProcessor>();

        force_add<MirrorProcessor>();

        force_add<ReverseProcessor>();


        force_add<AddNoiseProcessor>();

        force_add<AddSigmaNoiseProcessor>();

        force_add<AddRandomNoiseProcessor>();


        force_add<PhaseToCenterProcessor>();

        force_add<PhaseToCornerProcessor>();

        force_add<FourierToCenterProcessor>();

        force_add<FourierToCornerProcessor>();

        force_add<AutoMask2DProcessor>();

        force_add<AutoMask3DProcessor>();

        force_add<AutoMask3D2Processor>();

        force_add<AutoMaskDustProcessor>();

        force_add<AddMaskShellProcessor>();

        force_add<IterMultiMaskProcessor>();

        force_add<IterBinMaskProcessor>();

        force_add<AutoMaskAsymUnit>();


        force_add<CTFSNRWeightProcessor>();


        force_add<ToMassCenterProcessor>();

        force_add<ToCenterProcessor>();

        force_add<PhaseToMassCenterProcessor>();

        force_add<ACFCenterProcessor>();

//      force_add<SNRProcessor>();

        force_add<CTFCorrProcessor>();

        force_add<FSCFourierProcessor>();


        force_add<XGradientProcessor>();

        force_add<YGradientProcessor>();

        force_add<ZGradientProcessor>();


        force_add<ImageDivergenceProcessor>();

        force_add<GradientMagnitudeProcessor>();

        force_add<GradientDirectionProcessor>();

        force_add<LaplacianMagnitudeProcessor>();

        force_add<LaplacianDirectionProcessor>();

        force_add<SDGDProcessor>();


//      force_add<FileFourierProcessor>();


        force_add<SymSearchProcessor>();

        force_add<MaskPackProcessor>();

        force_add<StripeXYProcessor>();

        force_add<BadLineXYProcessor>();

        force_add<LocalNormProcessor>();


        force_add<IndexMaskFileProcessor>();

//      force_add<CoordinateMaskFileProcessor>();

        force_add<SetIsoPowProcessor>();

        force_add<SetSFProcessor>();

        force_add<MatchSFProcessor>();


        force_add<SmartMaskProcessor>();


        force_add<TestImageGaussian>();

        force_add<TestImagePureGaussian>();

        force_add<TestImageSinewave>();

        force_add<TestImageSphericalWave>();

        force_add<TestImageSinewaveCircular>();

        force_add<TestImageSquarecube>();

        force_add<TestImageCirclesphere>();

        force_add<TestImageAxes>();

        force_add<TestImageNoiseUniformRand>();

        force_add<TestImageNoiseGauss>();

        force_add<TestImageScurve>();

        force_add<TestImageCylinder>();

        force_add<TestImageDisc>();

        force_add<TestImageGradient>();

        force_add<TestTomoImage>();

        force_add<TestImageLineWave>();

        force_add<TestImageEllipse>();

        force_add<TestImageHollowEllipse>();

        force_add<TestImageFourierGaussianBand>();

        force_add<TestImageFourierNoiseGaussian>();

        force_add<TestImageFourierNoiseProfile>();


        force_add<TomoTiltEdgeMaskProcessor>();

        force_add<TomoTiltAngleWeightProcessor>();


        force_add<NewLowpassTopHatProcessor>();

        force_add<NewHighpassTopHatProcessor>();

        force_add<NewBandpassTopHatProcessor>();

        force_add<NewHomomorphicTopHatProcessor>();

        force_add<NewLowpassGaussProcessor>();

        force_add<NewHighpassGaussProcessor>();

        force_add<NewBandpassGaussProcessor>();

        force_add<NewHomomorphicGaussProcessor>();

        force_add<NewInverseGaussProcessor>();

        force_add<GaussZFourierProcessor>();

        force_add<LowpassRandomPhaseProcessor>();

        force_add<NewLowpassButterworthProcessor>();

        force_add<NewHighpassButterworthProcessor>();

        force_add<NewHomomorphicButterworthProcessor>();

        force_add<NewLowpassTanhProcessor>();

        force_add<NewHighpassTanhProcessor>();

        force_add<NewBandpassTanhProcessor>();

        force_add<NewHomomorphicTanhProcessor>();

        force_add<NewRadialTableProcessor>();

        force_add<InverseKaiserI0Processor>();

        force_add<InverseKaiserSinhProcessor>();

        force_add<CCDNormProcessor>();

        force_add<CTF_Processor>();

        force_add<SHIFTProcessor>();


//      force_add<WaveletProcessor>();

        force_add<EnhanceProcessor>();

        force_add<FFTProcessor>();

        force_add<RadialProcessor>();


        force_add<DirectionalSumProcessor>();

        force_add<ConvolutionKernelProcessor>();


        //Gorgon-related processors

        force_add<ModelEMCylinderProcessor>();

        force_add<ApplyPolynomialProfileToHelix>();

        force_add<BinarySkeletonizerProcessor>();

        force_add<RotateInFSProcessor>();

        force_add<CircularAverageBinarizeProcessor>();

        force_add<ObjDensityProcessor>();

        force_add<ObjLabelProcessor>();

        force_add<BwThinningProcessor>();

        force_add<BwMajorityProcessor>();

        force_add<PruneSkeletonProcessor>();

        force_add<GrowSkeletonProcessor>();


        force_add<ManhattanDistanceProcessor>();

        force_add<BinaryDilationProcessor>();

        force_add<BinaryErosionProcessor>();

        force_add<BinaryClosingProcessor>();

        force_add<BinaryOpeningProcessor>();

        force_add<BinaryMorphGradientProcessor>();

        force_add<BinaryExternalGradientProcessor>();

        force_add<BinaryInternalGradientProcessor>();

        force_add<BinaryTopHatProcessor>();

        force_add<BinaryBlackHatProcessor>();

        force_add<ZThicknessProcessor>();

        force_add<ReplaceValuefromListProcessor>();

        force_add<PolyMaskProcessor>();

        force_add<AmpMultProcessor>();


//#ifdef EMAN2_USING_CUDA

//      force_add<CudaMultProcessor>();

//      force_add<CudaCorrelationProcessor>();

//#endif // EMAN2_USING_CUDA


//      force_add<XYZProcessor>();

}


void FiniteProcessor::process_pixel(float *x) const

{

        if ( !Util::goodf(x) ) {

                *x = to;

        }

}


EMData* Processor::process(const EMData * const image)

{

        EMData * result = image->copy();

//      printf("Default copy called\n");

        process_inplace(result);

        return result;

}


void ImageProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL image");

                return;

        }


        EMData *processor_image = create_processor_image();


        if (image->is_complex()) {

                (*image) *= *processor_image;

        }

        else {

                EMData *fft = image->do_fft();

                (*fft) *= (*processor_image);

                EMData *ift = fft->do_ift();


                float *data = image->get_data();

                float *t = data;

                float *ift_data = ift->get_data();


                data = ift_data;

                ift_data = t;


                ift->update();


                if( fft )

                {

                        delete fft;

                        fft = 0;

                }


                if( ift )

                {

                        delete ift;

                        ift = 0;

                }

        }


        image->update();

}


#define FFTRADIALOVERSAMPLE 4

void FourierProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        preprocess(image);


        int array_size = FFTRADIALOVERSAMPLE * image->get_ysize();

        float step=0.5f/array_size;


        bool return_radial=(bool)params.set_default("return_radial",0);

        vector < float >yarray(array_size);


        create_radial_func(yarray);


        if (image->is_complex()) {

                image->apply_radial_func(0, step, yarray);

        }

        else {

                EMData *fft = image->do_fft();

                fft->apply_radial_func(0, step, yarray);

                EMData *ift = fft->do_ift();


                memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));


                //ift->update(); Unecessary


                if( fft )

                {

                        delete fft;

                        fft = 0;

                }


                if( ift )

                {

                        delete ift;

                        ift = 0;

                }

        }

        if (return_radial) image->set_attr("filter_curve",yarray);


        image->update();

}


void FourierAnlProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        preprocess(image);


//      int array_size = FFTRADIALOVERSAMPLE * image->get_ysize();

//      float step=0.5f/array_size;

//

//      vector < float >yarray(array_size);


        bool return_radial=(bool)params.set_default("return_radial",0);

        bool interpolate=(bool)params.set_default("interpolate",1);


        float cornerscale;

        if (image->get_zsize()>1) cornerscale=sqrt(3.0);

        else cornerscale=sqrt(2.0);


        if (image->is_complex()) {

                vector <float>yarray = image->calc_radial_dist(floor(image->get_ysize()*cornerscale/2),0,1.0,1);

                create_radial_func(yarray,image);

                image->apply_radial_func(0, 1.0f/image->get_ysize(), yarray,interpolate);

                if (return_radial) image->set_attr("filter_curve",yarray);

        }

        else {

                EMData *fft = image->do_fft();

                vector <float>yarray = fft->calc_radial_dist((int)floor(fft->get_ysize()*cornerscale/2.0),0,1.0,1);

                create_radial_func(yarray,fft);

                // 4/30/10 stevel turned off interpolation to fix problem with matched filter

                // 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

                fft->apply_radial_func(0, 1.0f/image->get_ysize(), yarray,interpolate);

                EMData *ift = fft->do_ift();


                memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));

                if (return_radial) image->set_attr("filter_curve",yarray);


//              for (int i=0; i<yarray.size(); i++) printf("%d\t%f\n",i,yarray[i]);


                //ift->update(); Unecessary


                delete fft;

                delete ift;


        }


        image->update();

}


EMData* MaskPackProcessor::process(const EMData *image) {

        EMData *mask=(EMData *)params["mask"];

        bool unpack = (int)params.set_default("unpack",0);


        if ((float)mask->get_attr("mean_nonzero")!=1.0f) throw InvalidParameterException("MaskPackProcessor requires a binary mask");


        int n_nz=(int)mask->get_attr("square_sum");             // true only for a binary mask


        EMData *ret = 0;

        if (unpack) {

                ret = new EMData(mask->get_xsize(),mask->get_ysize(),mask->get_zsize());

                ret->to_zero();


                size_t ix=0;

                for (size_t i=0; i<ret->get_size(); i++) {

                        if (mask->get_value_at_index(i)) ret->set_value_at_index(i,image->get_value_at_index(ix++));

                }


        } else {

                ret = new EMData(n_nz,1,1);

                ret->to_zero();


                size_t ix=0;

                for (size_t i=0; i<image->get_size(); i++) {

                        if (mask->get_value_at_index(i)) ret->set_value_at_index(ix++,image->get_value_at_index(i));

                }

        }


        ret->update();

        return ret;

}


// Looks like this hasn't actually been written yet...

void GridKernelFixProcessor::process_inplace(EMData *image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        string mode=(string)params["mode"];


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        // make an empty FFT object.

        EMData *fft=new EMData((nx&0xfffffffe)+2,ny,nz);

        fft->set_complex(1);

        fft->set_ri(1);

        if (nx&1) fft->set_fftpad(1);

        fft->to_zero();


        // Copy the kernel. Note that the kernel is 3x oversampled, so we are just computing an approximate kernel locally

        if (mode=="gridding_5") {

                for (int z=-2; z<3; z++) {

                        for (int y=-2; y<3; y++) {

                                for (int x=0; x<3; x++) {

                                        fft->set_complex_at(x,y,z,FourierInserter3DMode7::kernel[x*3][abs(y)*3][abs(z)*3]);

                                }

                        }

                }

        }

        else if (mode=="gridding_7") {

                for (int z=-3; z<4; z++) {

                        for (int y=-3; y<4; y++) {

                                for (int x=0; x<4; x++) {

                                        fft->set_complex_at(x,y,z,FourierInserter3DMode11::kernel[x*3][abs(y)*3][abs(z)*3]);

                                }

                        }

                }

        }

        else throw InvalidParameterException("Gridding kernel correction of unknown mode, only gridding_5 or gridding_7 allowed");


        EMData *real=fft->do_ift();                                     // this is the kernel ift

        real->process_inplace("xform.phaseorigin.tocenter");

        real->mult(2.0f/(float)real->get_attr("maximum"));

        real->process_inplace("math.reciprocal");       // reciprocal to make a correction volume

        real->process_inplace("threshold.clampminmax",Dict("minval",-4.0f,"maxval",4.0f));      // block overcorrection of noise near the edges

//      real->write_image("invkernel.hdf");


        image->mult(*real);             // apply the correction

        delete real;

        delete fft;

}


// Looks like this hasn't actually been written yet...

void AzSharpProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float az_scale=(float)params.set_default("az_scale",1.0);

        float cornerscale;

        EMData *fft;


        if (image->is_complex()) fft=image;

        else EMData *fft = image->do_fft();


        int nx=fft->get_xsize();

        int ny=fft->get_ysize();

        int nz=fft->get_zsize();

        if (nz!=1) {

        }

        else {


        }

        if (image->is_complex()) {

                EMData *ift = fft->do_ift();


                memcpy(image->get_data(),ift->get_data(),ift->get_xsize()*ift->get_ysize()*ift->get_zsize()*sizeof(float));


                delete fft;

                delete ift;

        }


        image->update();

}


void FFTPeakProcessor::process_inplace(EMData * image)

{

        EMData *fft;


        if (!image) throw InvalidParameterException("FFTPeakProcessor: no image provided");

        if (!image->is_complex()) fft = image->do_fft();

        else fft = image;


        int nx=fft->get_xsize();

        int ny=fft->get_ysize();

        int nz=fft->get_zsize();

        float thresh_sigma = (float)params.set_default("thresh_sigma", 1.0);

        bool removepeaks = (bool)params.set_default("removepeaks",0);

        bool to_mean = (bool)params.set_default("to_mean",0);


        vector<float> amp=fft->calc_radial_dist(nx/2,0,1,0);    // amplitude

        vector<float> thr=fft->calc_radial_dist(nx/2,0,1,1);    // start with inten (avg of amp^2)


        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


        if (nz>1) {

                for (int z=-nz/2; z<nz/2; z++) {

                        for (int y=-ny/2; y<ny/2; y++) {

                                for (int x=0; x<nx/2; x++) {

                                        float r2=Util::hypot3(x,y,z);

                                        int r=int(r2);

                                        if (r>=nx/2) continue;

                                        complex<float> v=fft->get_complex_at(x,y,z);

                                        float va=std::abs(v);

                                        if (va>=thr[r]) {

                                                if (removepeaks) {

                                                        if (to_mean) {

                                                                v*=amp[r]/va;

                                                                fft->set_complex_at(x,y,z,v);

                                                        }

                                                        else fft->set_complex_at(x,y,z,0.0);

                                                }

                                        }

                                        else {

                                                if (!removepeaks) {

                                                        if (to_mean) {

                                                                v*=amp[r]/va;

                                                                fft->set_complex_at(x,y,z,v);

                                                        }

                                                        else fft->set_complex_at(x,y,z,0.0);

                                                }

                                        }

                                }

                        }

                }

        }

        else {

                for (int y=-ny/2; y<ny/2; y++) {

                        for (int x=0; x<nx/2; x++) {

                                int r=Util::hypot_fast_int(x,y);

                                if (r>=nx/2) continue;

                                complex<float> v=fft->get_complex_at(x,y);

                                float va=std::abs(v);

                                if (va>=thr[r]) {

                                        if (removepeaks) {

                                                if (to_mean) {

                                                        v*=amp[r]/va;

                                                        fft->set_complex_at(x,y,v);

                                                }

                                                else fft->set_complex_at(x,y,0.0);

                                        }

                                }

                                else {

                                        if (!removepeaks) {

                                                if (to_mean) {

                                                        v*=amp[r]/va;

                                                        fft->set_complex_at(x,y,v);

                                                }

                                                else fft->set_complex_at(x,y,0.0);

                                        }

                                }

                        }

                }

        }


        if (fft!=image) {

                EMData *ift=fft->do_ift();

                memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));

                delete fft;

                delete ift;

        }

        image->update();


//      image->update();

}


void FFTConeProcessor::process_inplace(EMData * image)

{

        EMData *fft;


        if (!image) throw InvalidParameterException("FFTConeProcessor: no image provided");

        if (!image->is_complex()) fft = image->do_fft();

        else fft = image;


        int nx=fft->get_xsize();

        int ny=fft->get_ysize();

        int nz=fft->get_zsize();

        if (nz==1) throw ImageDimensionException("FFTConeProcessor only works on 3-D images");


        float angle = (float)params.set_default("angle",15.0f);

        float rmin = (float)params.set_default("rmin",1.0f);


        for (int z=-nz/2; z<nz/2; z++) {

                for (int y=-ny/2; y<ny/2; y++) {

                        for (int x=0; x<nx/2; x++) {

                                float ang=0;

                                if (x!=0 ||y!=0) ang=90.0-atan(fabs(float(z)/nz)/hypot(float(x)/nx,float(y)/ny))*180.0/M_PI;

                                if (ang>angle || Util::hypot3(x,y,z)<rmin) continue;

                                fft->set_complex_at(x,y,z,0.0f);

                        }

                }

        }


        if (fft!=image) {

                EMData *ift=fft->do_ift();

                memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));

                delete fft;

                delete ift;

        }

        image->update();


//      image->update();

}


void FFTWedgeProcessor::process_inplace(EMData * image)

{

        EMData *fft;


        if (!image) throw InvalidParameterException("FFTWedgeProcessor: no image provided");

        if (!image->is_complex()) fft = image->do_fft();

        else fft = image;


        int nx=fft->get_xsize();

        int ny=fft->get_ysize();

        int nz=fft->get_zsize();

        if (nz==1) throw ImageDimensionException("FFTWedgeProcessor only works on 3-D images");


        float anglemin = (float)params.set_default("anglemin",-30.0f);

        float anglemax = (float)params.set_default("anglemax",30.0f);

        float rmin = (float)params.set_default("rmin",1.0f);


        for (int z=-nz/2; z<nz/2; z++) {

                for (int y=-ny/2; y<ny/2; y++) {

                        for (int x=0; x<nx/2; x++) {

                                float ang=90.0f;

                                if (z!=0) ang=atan((float(y)/ny)/fabs(float(z)/nz))*180.0/M_PI;

                                if (ang<anglemin||ang>anglemax || Util::hypot3(x,y,z)<rmin) continue;

                                fft->set_complex_at(x,y,z,0.0f);

                        }

                }

        }


        if (fft!=image) {

                EMData *ift=fft->do_ift();

                memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));

                delete fft;

                delete ift;

        }

        image->update();


//      image->update();

}


void WedgeFillProcessor::process_inplace(EMData * image)

{

        if (!image) throw InvalidParameterException("WedgeFillProcessor: no image provided");

        if (!image->is_complex()) throw ImageFormatException("WedgeFillProcessor: target image must be complex");


        EMData *source=(EMData *)params.set_default("fillsource",(EMData *)NULL);

//      if (!source) throw InvalidParameterException("WedgeFillProcessor: fillsource required");


        if (source && !source->is_complex()) throw ImageFormatException("WedgeFillProcessor: fillsource must be complex");

        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");


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();

        float thresh_sigma = (float)params.set_default("thresh_sigma", 0.5);

        bool dosigma = 1 ? thresh_sigma>0.0 : 0;


        float maxtilt = (float)params.set_default("maxtilt", 90.0);

        bool dotilt = 1 ? maxtilt <90.0 : 0;

        maxtilt*=M_PI/180.0;


        vector<float> sigmaimg;

        if (dosigma) {

                sigmaimg=image->calc_radial_dist(nx/2,0,1,4);

                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

        }


        vector<int> realpixel(nx/2);

        for (int i=0; i<nx/2; i++) realpixel[i]=0;


        for (int z=0; z<nz; z++) {

                for (int y=0; y<ny; y++) {

                        for (int x=0; x<nx; x+=2) {

                                float r2=Util::hypot3(x/2,y<ny/2?y:ny-y,z<nz/2?z:nz-z); // origin at 0,0; periodic

                                int r=int(r2);

                                if (r<3) continue;              // too few points at r<3 to consider any "missing"


                                float tilt = 0.0;

                                if (dotilt) tilt=atan2((float)(z<nz/2?z:nz-z),(float)(x/2));


                                float v1r=image->get_value_at(x,y,z);

                                float v1i=image->get_value_at(x+1,y,z);

                                float v1=Util::square_sum(v1r,v1i);

//                              if (r<10) printf("%d %d %d %d\t%1.3g %1.3g\n",x,y,z,r,v1,sigmaimg[r]);

                                if ((!dosigma || v1>sigmaimg[r]) && r<nx/2 && (!dotilt || fabs(tilt)<maxtilt)){

                                        realpixel[r]++;

                                        continue;

                                }


                                if (!source) {

                                        image->set_value_at_fast(x,y,z,0);

                                        image->set_value_at_fast(x+1,y,z,0);

                                }

                                else {

                                        image->set_value_at_fast(x,y,z,source->get_value_at(x,y,z));

                                        image->set_value_at_fast(x+1,y,z,source->get_value_at(x+1,y,z));

                                }

                        }

                }

        }


        image->set_attr("real_pixels", realpixel);

        image->update();

}


void LowpassAutoBProcessor::create_radial_func(vector < float >&radial_mask,EMData *image) const{

        float apix=(float)image->get_attr("apix_x");

        int verbose=(int)params["verbose"];

//      int adaptnoise=params.set_default("adaptnoise",0);

        float noisecutoff=(float)params.set_default("noisecutoff",1.0/6.0);

        if (apix<=0 || apix>7.0f) throw ImageFormatException("apix_x > 7.0 or <0");

        float ds=1.0f/(apix*image->get_ysize());        // 0.5 is because radial mask is 2x oversampled

        unsigned int start=(int)floor(1.0/(15.0*ds));

        unsigned int end=radial_mask.size()-2;

        if (noisecutoff>0) end=(int)floor(noisecutoff/ds);

        if (end>radial_mask.size()-2) {

                if (verbose) printf("WARNING: specified noisecutoff too close to Nyquist, reset !");

                end=radial_mask.size()-2;

        }

        if (end<start+2) {

                printf("WARNING: noise cutoff too close to 15 A ! Results will not be good...");

                start=end-5;

        }


        FILE *out=NULL;

        if (verbose>2)  {

                printf("Autob -> %d - %d  ds=%g apix=%g rdlmsk=%d\n",start,end,ds,apix,int(radial_mask.size()));

                out=fopen("fitplot.txt","w");

        }

        int N=(radial_mask.size()-start-2);

        float *x=(float *)malloc(N*sizeof(float));

        float *y=(float *)malloc(N*sizeof(float));

        float *dy=(float *)malloc(N*sizeof(float));

        for (unsigned int i=start; i<radial_mask.size()-2; i++ ) {              // -2 is arbitrary because sometimes the last pixel or two have funny values

                x[i-start]=ds*ds*i*i;

                if (radial_mask[i]>0) y[i-start]=log(radial_mask[i]); // ok

                else if (i>start) y[i-start]=y[i-start-1];              // not good

                else y[i-start]=0.0;                                                    // bad

                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

                if (out) fprintf(out,"%f\t%f\n",x[i-start],y[i-start]);

        }

        if (out) fclose(out);


        float slope=0,intercept=0;

        Util::calc_least_square_fit(end-start, x,y,&slope,&intercept,1);


        if (verbose) printf("slope=%f  intercept=%f\n",slope,intercept);


        float B=(float)params["bfactor"]+slope*4.0f/2.0f;               // *4 is for Henderson definition, 2.0 is for intensity vs amplitude

        float B2=(float)params["bfactor"];


        if (verbose) printf("User B = %1.2f      Corrective B = %1.2f  Total B=%1.3f\n",(float)params["bfactor"],slope*2.0,B);


        float cutval=exp(-B*pow(end*ds,2.0f)/4.0f)/exp(-B2*pow(end*ds,2.0f)/4.0f);

        for (unsigned int i=0; i<radial_mask.size(); i++) {

                if (i<=end) radial_mask[i]=exp(-B*pow(i*ds,2.0f)/4.0f);

                else radial_mask[i]=cutval*exp(-B2*pow(i*ds,2.0f)/4.0f);

        }

        if (verbose>1) Util::save_data(0,ds,radial_mask,"filter.txt");


        free(x);

        free(y);

        free(dy);

 }


void CCCSNRProcessor::process_inplace(EMData *image) {

        size_t nxyz = image->get_size();

        int mode=params.set_default("wiener",0);

        float scale=params.set_default("scalesnr",2.0f);


        for (size_t i=0; i<nxyz; i++) {

                float v=image->get_value_at_index(i);

                float snr=(v>=.9999)?10000.0f:scale*v/(1.0f-v);

                if (snr<0) snr=0.0f;

                if (mode) image->set_value_at_index(i,snr/(1.0f+snr));

                else image->set_value_at_index(i,snr);

        }

}


void SNREvalProcessor::process_inplace(EMData * image)

{

        int ys=image->get_ysize();


        EMData *mask1=new EMData(ys,ys,1);

        mask1->process_inplace("mask.gaussian",Dict("outer_radius", ys/2.0));

        EMData *mask2=mask1->copy();

        mask2->mult(-1.0f);

        mask2->add(1.0);

        mask2->process_inplace("mask.decayedge2d",Dict("width",4));


/*


mask1=EMData(ys2,ys2,1)

                mask1.to_one()

                mask1.process_inplace("mask.gaussian",{"outer_radius":radius})

                mask2=mask1.copy()*-1+1

#               mask1.process_inplace("mask.decayedge2d",{"width":4})

                mask2.process_inplace("mask.decayedge2d",{"width":4})

                mask1.clip_inplace(Region(-(ys2*(oversamp-1)/2),-(ys2*(oversamp-1)/2),ys,ys))

                mask2.clip_inplace(Region(-(ys2*(oversamp-1)/2),-(ys2*(oversamp-1)/2),ys,ys))

                ratio1=mask1.get_attr("square_sum")/(ys*ys) #/1.035

                ratio2=mask2.get_attr("square_sum")/(ys*ys)

                masks[(ys,radius)]=(mask1,ratio1,mask2,ratio2)


*/

}


void Axis0FourierProcessor::process_inplace(EMData * image)

{

        EMData *fft;

        float *fftd;

        int f=0;

//      static float sum1=0,sum1a=0;

//      static double sum2=0,sum2a=0;


        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        if (!image->is_complex()) fft = image->do_fft();

        else fft = image;


        int neighbor=params.set_default("neighbor",0);

        float neighbornorm=params.set_default("neighbornorm",sqrt(2.0f));


        int nx=fft->get_xsize()/2;

        int ny=fft->get_ysize();

        if (params.set_default("x",1)) {

                if (neighbor) {

                                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

                }

                else {

                                for (int x=1; x<nx; x++) fft->set_complex_at(x,0,0.0f);

                }

        }

        if (params.set_default("y",1)) {

                if (neighbor) {

                                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

                }

                else {

                                for (int y=1; y<ny; y++) fft->set_complex_at(0,y,0.0f);

                }

        }


        if (fft!=image) {

//              fft->update();

                EMData *ift=fft->do_ift();

                memcpy(image->get_data(),ift->get_data(),(nx*2-2)*ny*sizeof(float));

                delete fft;

                delete ift;

        }


//      image->update();


}


void GaussZFourierProcessor::process_inplace(EMData * image)

{

        if(params.has_key("apix")) {

                image->set_attr("apix_x", (float)params["apix"]);

                image->set_attr("apix_y", (float)params["apix"]);

                image->set_attr("apix_z", (float)params["apix"]);

        }


        int xynoz = params.set_default("xynoz",0);


        const Dict dict = image->get_attr_dict();

        if( params.has_key("cutoff_freq") ) {

                float val =  (float)params["cutoff_freq"] * (float)dict["apix_x"];

                params["cutoff_abs"] = val;

        }

        else if( params.has_key("cutoff_pixels") ) {

                float val = ((float)params["cutoff_pixels"] / (float)dict["nx"]);

                params["cutoff_abs"] = val;

        }


        float omega = params["cutoff_abs"];

        float zcenter=params.set_default("centerfreq",0.0f);

        float hppix = params.set_default("hppix",-1.0f);


        omega = (omega<0?-1.0:1.0)*0.5f/omega/omega;


        EMData *fft;

        int f=0;


        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        if (!image->is_complex()) {

                fft = image->do_fft();

                f=1;

        }

        else {

                fft=image;

        }


        int nx=fft->get_xsize();

        int ny=fft->get_ysize();

        int nz=fft->get_ysize();

        omega /=(nz*nz)/4;

        zcenter=zcenter*(float)dict["apix_x"]*nz;


        if (xynoz) {

                for (int x=0; x<nx/2; x++) {

                        for (int z=(x==0?0:-nz/2); z<nz/2; z++) {

                                for (int y=(x==0&&z==0?0:-ny/2); y<ny/2; y++) {

                                        std::complex <float> v=fft->get_complex_at(x,y,z);

                                        float r=Util::hypot_fast(x,y);

                                        fft->set_complex_at(x,y,z,v*exp(-omega*(r-zcenter)*(r-zcenter))*(r>hppix?1.0f:Util::hypot3(x,y,z)/hppix));

                                }

                        }

                }

        }

        else {

                for (int z=-nz/2; z<nz/2; z++) {

                        for (int y=-ny/2; y<ny/2; y++) {

                                for (int x=0; x<nx/2; x++) {

                                        std::complex <float> v=fft->get_complex_at(x,y,z);

                                        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));

                                }

                        }

                }

        }


        if (f) {

                EMData *ift=fft->do_ift();

                memcpy(image->get_data(),ift->get_data(),(nx-2)*ny*nz*sizeof(float));

                delete fft;

                delete ift;

        }


}


void AmpweightFourierProcessor::process_inplace(EMData * image)

{

        EMData *fft;

        float *fftd;

        int f=0;

//      static float sum1=0,sum1a=0;

//      static double sum2=0,sum2a=0;


        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        if (!image->is_complex()) {

                fft = image->do_fft();

                fftd = fft->get_data();

                f=1;

        }

        else {

                fft=image;

                fftd=image->get_data();

        }

        float *sumd = NULL;

        if (sum) sumd=sum->get_data();

//printf("%d %d    %d %d\n",fft->get_xsize(),fft->get_ysize(),sum->get_xsize(),sum->get_ysize());

        size_t n = (size_t)fft->get_xsize()*fft->get_ysize()*fft->get_zsize();

        for (size_t i=0; i<n; i+=2) {

                float c;

                if (dosqrt) c=pow(fftd[i]*fftd[i]+fftd[i+1]*fftd[i+1],0.25f);

                else c = static_cast<float>(hypot(fftd[i],fftd[i+1]));

                if (c==0) c=1.0e-30f;   // prevents divide by zero in normalization

                fftd[i]*=c;

                fftd[i+1]*=c;

                if (sumd) { sumd[i]+=c; sumd[i+1]+=0; }


                // debugging

/*              if (i==290*1+12) {

                        sum1+=fftd[i];

                        sum2+=fftd[i];

                        printf("%f\t%f\t%f\t%f\t%f\t%f\n",sum1,sum2,fftd[i],fftd[i+1],sumd[i],c);

                }

                if (i==290*50+60) {

                        sum1a+=fftd[i];

                        sum2a+=fftd[i];

                        printf("%f\t%f\t%f\t%f\t%f\t%f\n",sum1a,sum2a,fftd[i],fftd[i+1],sumd[i],c);

        }*/

        }


        if (f) {

                fft->update();

                EMData *ift=fft->do_ift();

                memcpy(image->get_data(),ift->get_data(),n*sizeof(float));

                delete fft;

                delete ift;

        }


        sum->update();

        image->update();


}


void GaussSegmentProcessor::process_inplace(EMData *)

{

        printf("Process inplace not implemented. Please use process.\n");

        return;

}


EMData *GaussSegmentProcessor::process(const EMData * const image)

{


        float minratio = params.set_default("minratio",0.5f);

        int maxnseg = params.set_default("maxnseg",0);

        int skipseg = params.set_default("skipseg",0);

        float width = params.set_default("width",10.0f);

        int verbose = params.set_default("verbose",0);

        EMData *mask= (EMData *)params.set_default("mask",(EMData *)NULL);

        float apix=image->get_attr("apix_x");

        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        EMData * result = image->process("threshold.belowtozero",Dict("minval",0.0f));

        // The intent of these filters is to insure that the map effectively

//      result->process_inplace("filter.lowpass.gauss",Dict("cutoff_freq",0.45/width)); // 0.45 = sqrt(2)/pi, resolvability threshold

//      result->process_inplace("filter.lowpass.gauss",Dict("cutoff_freq",0.637/width));        // 0.637 = 2/pi, (may be Rayleigh criterion?)

        if (mask!=NULL) {

                result->mult(*mask);

                result->process_inplace("normalize.local",Dict("radius",nx*apix/(3.0f*width),"threshold",(float)result->get_attr("sigma_nonzero")*1.2));

        }


        XYData gsf;

        gsf.make_gauss(nx*4,1.0f/apix,0.45/width);

//      gsf.make_gauss(nx*4,1.0f/apix,0.637/width);

        result->process_inplace("filter.setstrucfac",Dict("apix",apix,"strucfac",&gsf));


        // Generate a Gaussian we can subtract out of the volume quickly

        int gs=2*width/(float)image->get_attr("apix_x");

        EMData gsub(gs,gs,gs);

        gsub.to_one();

        gsub.process_inplace("mask.gaussian",Dict("outer_radius",int(width/(2.0*apix))));


        if (verbose>2) {

                result->set_attr("render_bits",12);

                result->set_attr("render_min",(float)result->get_attr("minimum"));

                result->set_attr("render_max",(float)result->get_attr("maximum"));

                result->write_image("segfilt.hdf",0,EMUtil::IMAGE_UNKNOWN,false,0,EMUtil::EM_COMPRESSED);

        }


        EMData *cache = NULL;

        if (skipseg==2) cache=result->copy();


        // original version of the code. A bit slow but works very well

//      vector<float> centers;

//      vector<float> amps;

//      if (maxnseg==0) maxnseg=2000000000;

//      float maxval=0.0f;

//

//      while (amps.size()<maxnseg) {

//              IntPoint p = result->calc_max_location();

//              FloatPoint fp(p[0],p[1],p[2]);

//

//              float amp=result->get_value_at(p[0],p[1],p[2]);

//              if (amp<maxval*minratio) break;

//              amps.push_back(amp);

//              centers.push_back(p[0]);

//              centers.push_back(p[1]);

//              centers.push_back(p[2]);

//              if (amp>maxval) maxval=amp;             // really this should only happen once...

//

//              result->insert_scaled_sum(&gsub,fp,1.0,-amp);

//      }


        // This was an alternate version, to improve speed. While it is faster, the approximations its making

        // produce slightly less optimal results.

        vector<float> centers;

        vector<float> amps;

        if (maxnseg==0) maxnseg=2000000000;


        float thr=(float)result->get_attr("maximum")*minratio;

        while (amps.size()<maxnseg) {

                if ((float)result->get_attr("maximum")<=thr) break;

                // We find all peak values greater than our threshold

                vector<Pixel> pixels=result->calc_highest_locations(thr);

                if (pixels.size()==0) break;

                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());


                for (int i=0; i<pixels.size(); i++) {

                        IntPoint p = pixels[i].get_point();

                        FloatPoint fp(p[0],p[1],p[2]);


                        float amp=result->get_value_at(p[0],p[1],p[2]);

                        if (amp!=pixels[i].get_value()) continue;       // skip any values near a value we've just changed, should come back in the next round

                        //                      printf("%d %d %d     %f  %f   %f\n",p[0],p[1],p[2],amp,pixels[i].get_value(),pixels[i+1].get_value());

//                      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

                        if (amp<thr) continue;

//                      if (amp<thr) { maxnseg=amps.size(); break; }

                        amps.push_back(amp);

                        centers.push_back(p[0]);

                        centers.push_back(p[1]);

                        centers.push_back(p[2]);


                        result->insert_scaled_sum(&gsub,fp,1.0,-amp);

                }

        }


        if (verbose) printf("Found %d centers\n",amps.size());


        if (verbose>2) {

                result->set_attr("render_bits",12);

                result->set_attr("render_min",(float)result->get_attr("minimum"));

                result->set_attr("render_max",(float)result->get_attr("maximum"));

                result->write_image("segfilt.hdf",1,EMUtil::IMAGE_UNKNOWN,false,0,EMUtil::EM_COMPRESSED);

        }

        if (!skipseg) {

        // after we have our list of centers classify pixels

        for (int z=0; z<nz; z++) {

                for (int y=0; y<ny; y++) {

                        for (int x=0; x<nz; x++) {

//                              if (image->get_value_at(x,y,z)<thr) {

//                                      result->set_value_at(x,y,z,-1.0);               //below threshold -> -1 (unclassified)

//                                      continue;

//                              }

                                int bcls=-1;                    // best matching class

                                float bdist=(float)(nx+ny+nz);  // distance for best class

                                for (unsigned int c=0; c<centers.size()/3; c++) {

                                        float d=Util::hypot3(x-centers[c*3],y-centers[c*3+1],z-centers[c*3+2]);

                                        if (d<bdist) { bdist=d; bcls=c; }

                                }

                                result->set_value_at(x,y,z,(float)bcls);                // set the pixel to the class number

                        }

                }

        }

        if (verbose) printf("segmented\n");

        }


        //if skipseg is set to 2, we return the filtered, unsegmented map (with valid centers)

        if (skipseg==2) {

                delete result;

                result=cache;

        }


        result->set_attr("segment_centers",centers);

        result->set_attr("segment_amps",amps);


        return result;

}


void DistanceSegmentProcessor::process_inplace(EMData *)

{

        printf("Process inplace not implemented. Please use process.\n");

        return;

}


EMData *DistanceSegmentProcessor::process(const EMData * const image)

{

        EMData * result = image->copy();


        float thr = params.set_default("thr",0.9f);

        float minsegsep = params.set_default("minsegsep",5.0f);

        float maxsegsep = params.set_default("maxsegsep",5.1f);

        int verbose = params.set_default("verbose",0);


        vector<Pixel> pixels=image->calc_highest_locations(thr);


        vector<float> centers(3);       // only 1 to start

        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();

//      int nxy=nx*ny;


        // seed the process with the highest valued point

        centers[0]=(float)pixels[0].x;

        centers[1]=(float)pixels[0].y;

        centers[2]=(float)pixels[0].z;

        pixels.erase(pixels.begin());


        // outer loop. We add one center per iteration

        // This is NOT a very efficient algorithm, it assumes points are fairly sparse

        while (pixels.size()>0) {

                // iterate over pixels until we find a new center (then stop), delete any 'bad' pixels

                // no iterators because we remove elements


                for (unsigned int i=0; i<pixels.size(); i++) {


                        Pixel p=pixels[i];

                        // iterate over existing centers to see if this pixel should be removed ... technically we only should need to check the last center

                        for (unsigned int j=0; j<centers.size(); j+=3) {

                                float d=Util::hypot3(centers[j]-p.x,centers[j+1]-p.y,centers[j+2]-p.z);

                                if (d<minsegsep) {              // conflicts with existing center, erase

                                        pixels.erase(pixels.begin()+i);

                                        i--;

                                        break;

                                }

                        }

                }


                int found=0;

                for (unsigned int i=0; i<pixels.size() && found==0; i++) {

                        Pixel p=pixels[i];


                        // iterate again to see if this may be a new valid center. Start at the end so we tend to build chains

                        for (unsigned int j=centers.size()-3; j>0; j-=3) {

                                float d=Util::hypot3(centers[j]-p.x,centers[j+1]-p.y,centers[j+2]-p.z);

                                if (d<maxsegsep) {              // we passed minsegsep question already, so we know we're in the 'good' range

                                        centers.push_back((float)p.x);

                                        centers.push_back((float)p.y);

                                        centers.push_back((float)p.z);

                                        pixels.erase(pixels.begin()+i); // in the centers list now, don't need it any more

                                        found=1;

                                        break;

                                }

                        }

                }


                // If we went through the whole list and didn't find one, we need to reseed again

                if (!found && pixels.size()) {

                        if (verbose) printf("New chain\n");

                        centers.push_back((float)pixels[0].x);

                        centers.push_back((float)pixels[0].y);

                        centers.push_back((float)pixels[0].z);

                        pixels.erase(pixels.begin());

                }


                if (verbose) printf("%d points found\n",(int)(centers.size()/3));

        }


        // after we have our list of centers classify pixels

        for (int z=0; z<nz; z++) {

                for (int y=0; y<ny; y++) {

                        for (int x=0; x<nz; x++) {

                                if (image->get_value_at(x,y,z)<thr) {

                                        result->set_value_at(x,y,z,-1.0);               //below threshold -> -1 (unclassified)

                                        continue;

                                }

                                int bcls=-1;                    // best matching class

                                float bdist=(float)(nx+ny+nz);  // distance for best class

                                for (unsigned int c=0; c<centers.size()/3; c++) {

                                        float d=Util::hypot3(x-centers[c*3],y-centers[c*3+1],z-centers[c*3+2]);

                                        if (d<bdist) { bdist=d; bcls=c; }

                                }

                                result->set_value_at(x,y,z,(float)bcls);                // set the pixel to the class number

                        }

                }

        }


        result->set_attr("segment_centers",centers);


        return result;

}


EMData* ApplySymProcessor::process(const EMData * const image)

{

        string s=(string)params.set_default("sym","c1");

        if (s.length()<2) return image->copy();

        int n=atoi(s.c_str()+1);

        if ((s[0]=='c' || s[0]=='C') && n==1) return image->copy();


        Averager* imgavg = Factory<Averager>::get((string)params.set_default("averager","mean"));


        if (image->get_zsize()==1) {

                if (s[0]!='c' && s[0]!='C') throw ImageDimensionException("xform.applysym: Cn symmetry required for 2-D symmetrization");

                if (n<=0) throw InvalidValueException(n,"xform.applysym: Cn symmetry, n>0");


                for (int i=0; i<n; i++) {

                        Transform t(Dict("type","2d","alpha",(float)(i*360.0f/n)));

                        EMData* transformed = image->process("xform",Dict("transform",&t));

                        imgavg->add_image(transformed);

                        delete transformed;

                }

                EMData *ret=imgavg->finish();

                delete imgavg;

                return ret;

        }


        Symmetry3D* sym = Factory<Symmetry3D>::get((string)params.set_default("sym","c1"));

        vector<Transform> transforms = sym->get_syms();


        for(vector<Transform>::const_iterator trans_it = transforms.begin(); trans_it != transforms.end(); trans_it++) {

                Transform t = *trans_it;

                EMData* transformed = image->process("xform",Dict("transform",&t));

                imgavg->add_image(transformed);

                delete transformed;

        }

        EMData *ret=imgavg->finish();

        delete imgavg;

        return ret;

}


void ApplySymProcessor::process_inplace(EMData* image)

{

        EMData *tmp=process(image);

        memcpy(image->get_data(),tmp->get_data(),(size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));

        delete tmp;

        image->update();

        return;

}


EMData* KmeansSegmentProcessor::process(const EMData * const image)

{

        EMData * result = image->copy();


        int nseg = params.set_default("nseg",12);

        float thr = params.set_default("thr",-1.0e30f);

        int ampweight = params.set_default("ampweight",1);

        float maxsegsize = params.set_default("maxsegsize",10000.0f);

        float minsegsep = params.set_default("minsegsep",0.0f);

        int maxiter = params.set_default("maxiter",100);

        int maxvoxmove = params.set_default("maxvoxmove",25);

        int verbose = params.set_default("verbose",0);

        bool pseudoatom = params.set_default("pseudoatom",0);

        float sep = params.set_default("sep",3.78f);


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();

//      int nxy=nx*ny;

//      image->process_inplace("threshold.belowtozero");

        if (thr==-1.0e30f){

                thr=float(image->get_attr("mean_nonzero"))+ 1.0 *float(image->get_attr("sigma_nonzero"));

                printf("Estimated map threshold: %4f\n", thr);

        }

        // seed

        vector<float> centers(nseg*3);

        vector<float> count(nseg);

        // Alternative seeding method for paudoatom generation. Seed on the gird.

        if (pseudoatom){

                float ax=image->get_attr("apix_x");

                sep/=ax;

                if (verbose) printf("Seeding .....\n");

                int sx=int(nx/sep)+1,sy=int(ny/sep)+1,sz=int(nz/sep)+1;

                EMData m(sx,sy,sz);

                EMData mcount(sx,sy,sz);

                for (int i=0; i<nx; i++){

                        for (int j=0; j<ny; j++){

                                for (int k=0; k<nz; k++){

                                        int ni=(i/sep),nj=(j/sep),nk=(k/sep);

                                        float v=image->get_value_at(i,j,k);

                                        if (v>thr){

                                                m.set_value_at(ni,nj,nk,(m.get_value_at(ni,nj,nk)+v));

                                                mcount.set_value_at(ni,nj,nk,(mcount.get_value_at(ni,nj,nk)+1));

                                        }

                                }

                        }

                }

                m.div((nx/sx)*(ny/sy)*(nz/sz));

                int nsum=0;

                float l=image->get_attr("minimum"),r=image->get_attr("maximum"),th=0;

                while (abs(nsum-nseg)>0){

                        th=(l+r)/2;

                        nsum=0;

                        for (int i=0; i<sx; i++){

                                for (int j=0; j<sy; j++){

                                        for (int k=0; k<sz; k++){

                                                if (m.get_value_at(i,j,k)>th)  nsum+=1;

                                        }

                                }

                        }

//                      if (verbose) printf("%3f\t %3f\t %3f,\t %4d\t %4d\n", l,th,r,nsum,nseg);

                        if (nsum>nseg) l=th;

                        if (nsum<nseg) r=th;

                        if ((r-l)<.0001)break;

                }

//              nseg=nsum;

                if (verbose) printf("%3d pseudoatoms seeded at threshold value %3f\n", nsum, th);

                int q=0;

                for (int i=0; i<sx; i++){

                        for (int j=0; j<sy; j++){

                                for (int k=0; k<sz; k++){

                                        if (m.get_value_at(i,j,k)>th){

                                                if(q<nseg*3){

                                                        centers[q]=      float(i+.5)*sep;

                                                        centers[q+1]=float(j+.5)*sep;

                                                        centers[q+2]=float(k+.5)*sep;

                                                        q+=3;

                                                }

                                        }

                                }

                        }

                }


        }

        // Default: random seeding.

        else{

                for (int i=0; i<nseg*3; i+=3) {

                        centers[i]=      Util::get_frand(0.0f,(float)nx);

                        centers[i+1]=Util::get_frand(0.0f,(float)ny);

                        centers[i+2]=Util::get_frand(0.0f,(float)nz);

                }

        }


        for (int iter=0; iter<maxiter; iter++) {

                // **** classify

                size_t pixmov=0;                // count of moved pixels

                for (int z=0; z<nz; z++) {

                        for (int y=0; y<ny; y++) {

                                for (int x=0; x<nx; x++) {

                                        if (image->get_value_at(x,y,z)<thr) {

                                                result->set_value_at(x,y,z,-1.0);               //below threshold -> -1 (unclassified)

                                                continue;

                                        }

                                        int bcls=-1;                    // best matching class

                                        float bdist=(float)(nx+ny+nz);  // distance for best class

                                        for (int c=0; c<nseg; c++) {

                                                float d=Util::hypot3(x-centers[c*3],y-centers[c*3+1],z-centers[c*3+2]);

                                                if (d<bdist) { bdist=d; bcls=c; }

                                        }

                                        if ((int)result->get_value_at(x,y,z)!=bcls) pixmov++;

                                        if (bdist>maxsegsize) result->set_value_at(x,y,z,-1);           // pixel is too far from any center

                                        else result->set_value_at(x,y,z,(float)bcls);           // set the pixel to the class number

                                }

                        }

                }


                // **** adjust centers

                for (int i=0; i<nseg*3; i++) centers[i]=0;

                for (int i=0; i<nseg; i++) count[i]=0;

                // weighted sums

                for (int z=0; z<nz; z++) {

                        for (int y=0; y<ny; y++) {

                                for (int x=0; x<nx; x++) {

                                        int cls = (int)result->get_value_at(x,y,z);

                                        if (cls==-1) continue;

                                        float w=1.0;

                                        if (ampweight) w=image->get_value_at(x,y,z);


                                        centers[cls*3]+=x*w;

                                        centers[cls*3+1]+=y*w;

                                        centers[cls*3+2]+=z*w;

                                        count[cls]+=w;

                                }

                        }

                }

                // now each becomes center of mass, or gets randomly reseeded

                int nreseed=0;

                for (int c=0; c<nseg; c++) {

                        // reseed

                        if (count[c]==0) {

                                nreseed++;

                                do {

                                        centers[c*3]=  Util::get_frand(0.0f,(float)nx);

                                        centers[c*3+1]=Util::get_frand(0.0f,(float)ny);

                                        centers[c*3+2]=Util::get_frand(0.0f,(float)nz);

                                } 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

                        }

                        // center of mass

                        else {

                                centers[c*3]/=count[c];

                                centers[c*3+1]/=count[c];

                                centers[c*3+2]/=count[c];

                        }

                }


                // with minsegsep, check separation

                if (minsegsep>0) {

                        for (int c1=0; c1<nseg-1; c1++) {

                                for (int c2=c1+1; c2<nseg; c2++) {

                                        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) {

                                                nreseed++;

                                                do {

                                                        centers[c1*3]=  Util::get_frand(0.0f,(float)nx);

                                                        centers[c1*3+1]=Util::get_frand(0.0f,(float)ny);

                                                        centers[c1*3+2]=Util::get_frand(0.0f,(float)nz);

                                                } while (image->get_value_at((int)centers[c1*3],(int)centers[c1*3+1],(int)centers[c1*3+2])<thr);

                                        }

                                }

                        }

                }


                if (verbose) printf("Iteration %3d: %6ld voxels moved, %3d classes reseeded\n",iter,pixmov,nreseed);

                if (nreseed==0 && pixmov<(size_t)maxvoxmove) break;             // termination conditions met

        }


        result->set_attr("segment_centers",centers);


        return result;

}


void KmeansSegmentProcessor::process_inplace(EMData *)

{

        printf("Process inplace not implemented. Please use process.\n");

        return;

}


void LinearPyramidProcessor::process_inplace(EMData *image) {


//      if (image->get_zsize()!=1) { throw ImageDimensionException("Only 2-D images supported"); }


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();

        float x0 = params.set_default("x0",nx/2);

        float y0 = params.set_default("y0",ny/2);

        float z0 = params.set_default("z0",nz/2);

        float xwidth = params.set_default("xwidth",nx);

        float ywidth = params.set_default("ywidth",ny);

        float zwidth = params.set_default("zwidth",nz);


        for (int z=0; z<nz; z++) {

                for (int y=0; y<ny; y++) {

                        for (int x=0; x<nx; x++) {

                                float xf=1.0-fabs(x-x0)*2.0/xwidth;

                                float yf=1.0-fabs(y-y0)*2.0/ywidth;

                                float zf=1.0-fabs(z-z0)*2.0/zwidth;

                                float val=0.0;

                                if (xf>0&&yf>0&&zf>0) val=xf*yf*zf;

                                image->mult_value_at_fast(x,y,z,val);

                        }

                }

        }

}


EMData * Wiener2DAutoAreaProcessor::process(const EMData * image)

{

// TODO NOT IMPLEMENTED YET !!!

        EMData *ret = 0;

        const EMData *fft;

//      float *fftd;

//      int f=0;


        if (!image) {

                LOGWARN("NULL Image");

                return ret;

        }

        throw NullPointerException("Processor not yet implemented");


//      if (!image->is_complex()) {

//              fft = image->do_fft();

//              fftd = fft->get_data();

//              f=1;

//      }

//      else {

//              fft=image;

//              fftd=image->get_data();

//      }


        return ret;

}


void Wiener2DAutoAreaProcessor::process_inplace(EMData *image) {

        EMData *tmp=process(image);

        memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));

        delete tmp;

        image->update();

        return;

}


EMData * Wiener2DFourierProcessor::process(const EMData *in)

{

        const EMData *in2 = 0;

        if (in->is_complex()) in2=in;

        else in=in->do_fft();


        EMData *filt = in->copy_head();

        Ctf *ictf = ctf;


        if (!ictf) ctf=(Ctf *)in->get_attr("ctf");


        ictf->compute_2d_complex(filt,Ctf::CTF_WIENER_FILTER);

        filt->mult(*in2);

        EMData *ret=filt->do_ift();


        delete filt;

        if (!in->is_complex()) delete in2;


        if(!ictf && ctf) {delete ctf; ctf=0;}

        return(ret);

/*      const EMData *fft;

        float *fftd;

        int f=0;


        if (!image) {

                LOGWARN("NULL Image");

                return ret;

        }


        if (!image->is_complex()) {

                fft = image->do_fft();

                fftd = fft->get_data();

                f=1;

        }

        else {

                fft=image;

                fftd=image->get_data();

        }

        powd=image->get_data();


        int bad=0;

        for (int i=0; i<image->get_xsize()*image->get_ysize(); i+=2) {

                snr=(fftd[i]*fftd[i]+fftd[i+1]*fftd[i+1]-powd[i])/powd[i];

                if (snr<0) { bad++; snr=0; }


        }


        print("%d bad pixels\n",snr);

*/      return ret;


}


void Wiener2DFourierProcessor::process_inplace(EMData *image) {

        EMData *tmp=process(image);

        memcpy(image->get_data(),tmp->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));

        delete tmp;

        image->update();

        return;

}


void LinearRampFourierProcessor::create_radial_func(vector < float >&radial_mask) const

{

        Assert(radial_mask.size() > 0);

        for (size_t i = 0; i < radial_mask.size(); i++) {

                radial_mask[i] = (float)i;

        }

}


void LowpassRandomPhaseProcessor::create_radial_func(vector < float >&radial_mask) const { };


void LowpassRandomPhaseProcessor::process_inplace(EMData *image)

{

        float cutoff=0;

        preprocess(image);

        if( params.has_key("cutoff_abs") ) {

                cutoff=(float)params["cutoff_abs"];

        }

        else {

                printf("A cutoff_* parameter is required by filter.lowpass.randomphase\n");

                return;

        }


        if (image->get_zsize()==1) {

                int flag=0;

                if (!image->is_complex()) { image->do_fft_inplace(); flag=1; }

                image->ri2ap();

                int nx=image->get_xsize();

                int ny=image->get_ysize();


                int z=0;

                float *data=image->get_data();

                for (int y=-ny/2; y<ny/2; y++) {

                        for (int x=0; x<nx/2+1; x++) {

                                if (hypot(x/float(nx),y/float(ny))>=cutoff) {

                                        size_t idx=image->get_complex_index_fast(x,y,z);                // location of the amplitude

                                        data[idx+1]=Util::get_frand(0.0f,(float)(M_PI*2.0));

                                }

                        }

                }


                image->ap2ri();


                if (flag) {

                        image->do_ift_inplace();

                        image->depad();

                }

        }

        else {          // 3D

                int flag=0;

                if (!image->is_complex()) { image->do_fft_inplace(); flag=1; }

                image->ri2ap();

                int nx=image->get_xsize();

                int ny=image->get_ysize();

                int nz=image->get_zsize();


                float *data=image->get_data();

                for (int z=-nz/2; z<nz/2; z++) {

                        for (int y=-ny/2; y<ny/2; y++) {

                                for (int x=0; x<nx/2; x++) {

                                        if (Util::hypot3(x/float(nx),y/float(ny),z/float(nz))>=cutoff) {

                                                size_t idx=image->get_complex_index_fast(x,y,z);                // location of the amplitude

                                                data[idx+1]=Util::get_frand(0.0f,(float)(M_PI*2.0));

                                        }

                                }

                        }

                }

                image->ap2ri();


                if (flag) {

                        image->do_ift_inplace();

                        image->depad();

                }

        }

}


void HighpassAutoPeakProcessor::preprocess(EMData * image)

{

        if(params.has_key("apix")) {

                image->set_attr("apix_x", (float)params["apix"]);

                image->set_attr("apix_y", (float)params["apix"]);

                image->set_attr("apix_z", (float)params["apix"]);

        }


        const Dict dict = image->get_attr_dict();


        if( params.has_key("cutoff_abs") ) {

                highpass = params["cutoff_abs"];

        }

        else if( params.has_key("cutoff_freq") ) {

                highpass = (float)params["cutoff_freq"] * (float)dict["apix_x"] * (float)dict["ny"] / 2.0f;

        }

        else if( params.has_key("cutoff_pixels") ) {

                highpass = (float)params["cutoff_pixels"] / (float)dict["nx"];

        }

}


void HighpassAutoPeakProcessor::create_radial_func(vector < float >&radial_mask, EMData *) const

{

        unsigned int c;


//      for (unsigned int i=0; i<radial_mask.size(); i++) printf("%d\t%f\n",i,radial_mask[i]);

        for (c=2; c<radial_mask.size(); c++) if (radial_mask[c-1]<=radial_mask[c]) break;

        if (c>highpass) c=(unsigned int)highpass;               // the *2 is for the 2x oversampling


        radial_mask[0]=0.0;

//      for (int i=1; i<radial_mask.size(); i++) radial_mask[i]=(i<=c?radial_mask[c+1]/radial_mask[i]:1.0);

        for (unsigned int i=1; i<radial_mask.size(); i++) radial_mask[i]=(i<=c?0.0f:1.0f);


        printf("%f %d\n",highpass,c);

//      for (unsigned int i=0; i<radial_mask.size(); i++) printf("%d\t%f\n",i,radial_mask[i]);


}


void LinearRampProcessor::create_radial_func(vector < float >&radial_mask) const

{

        Assert(radial_mask.size() > 0);

        float x = 0.0f , step = 0.5f/radial_mask.size();

        size_t size=radial_mask.size();

        for (size_t i = 0; i < size; i++) {

                radial_mask[i] = intercept + ((slope - intercept) * i) / (size - 1.0f);

                x += step;

        }

}


void LoGFourierProcessor::create_radial_func(vector < float >&radial_mask) const

{


        Assert(radial_mask.size() > 0);

        float x = 0.0f , nqstep = 0.5f/radial_mask.size();

        size_t size=radial_mask.size();

        float var = sigma*sigma;

        for (size_t i = 0; i < size; i++) {

                radial_mask[i] = ((x*x - var)/var*var)*exp(-x*x/2*var);

                x += nqstep;

        }

}


void DoGFourierProcessor::create_radial_func(vector < float >&radial_mask) const

{


        Assert(radial_mask.size() > 0);

        float x = 0.0f , nqstep = 0.5f/radial_mask.size();

        size_t size=radial_mask.size();

        float norm = 1.0f/sqrt(2*M_PI);

        for (size_t i = 0; i < size; i++) {

                radial_mask[i] = norm*((1.0f/sigma1*exp(-x*x/(2.0f*sigma1*sigma1))) - (1.0f/sigma2*exp(-x*x/(2.0f*sigma2*sigma2))));

                x += nqstep;

        }

}


void RangeZeroProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float gauss_width = params.set_default("gauss_width",0.0f);


        if (gauss_width<=0) {

                size_t size = (size_t)image->get_xsize() *

                                        (size_t)image->get_ysize() *

                                        (size_t)image->get_zsize();

                float *data = image->get_data();


                for (size_t i = 0; i < size; ++i) {

                        if (data[i]>=minval && data[i]<=maxval) data[i]=0.0f;

                }

                image->update();

        }

        else {

                int nx = image->get_xsize();

                int ny = image->get_ysize();

                int nz = image->get_zsize();

                float wid=gauss_width/(ny*ny);


                for (int z=0; z<nz; z++) {

                        for (int y=0; y<ny; y++) {

                                for (int x=0; x<nx; x++) {

                                        float cor=exp(-Util::hypot3sq(x-nx/2,y-ny/2,z-nz/2)*wid);

                                        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);

                                }

                        }

                }

                image->update();

        }


}


void RealPixelProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        maxval = image->get_attr("maximum");

        mean = image->get_attr("mean");

        sigma = image->get_attr("sigma");


        calc_locals(image);


        size_t size = (size_t)image->get_xsize() *

                                  (size_t)image->get_ysize() *

                                  (size_t)image->get_zsize();

        float *data = image->get_data();


        for (size_t i = 0; i < size; ++i) {

                process_pixel(&data[i]);

        }

        image->update();

}


void CoordinateProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        maxval = image->get_attr("maximum");

        mean = image->get_attr("mean");

        sigma = image->get_attr("sigma");

        nx = image->get_xsize();

        ny = image->get_ysize();

        nz = image->get_zsize();

        is_complex = image->is_complex();


        calc_locals(image);


        if (!is_valid()) {

                return;

        }


        float *data = image->get_data();

        size_t i = 0;


        for (int z = 0; z < nz; z++) {

                for (int y = 0; y < ny; y++) {

                        for (int x = 0; x < nx; x++) {

                                process_pixel(&data[i], x, y, z);

                                ++i;

                        }

                }

        }

        image->update();

}


void PaintProcessor::process_inplace(EMData *image) {

        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        if (nz==1) {

                float r;

                for (int j=(y<r2?0:y-r2); j<(y+r2>ny?ny:y+r2); j++) {

                        for (int i=(x<r2?0:x-r2); i<(x+r2>nx?nx:x+r2); i++) {

                                r=sqrt(float(Util::square(i-x)+Util::square(j-y)));

                                if (r>r2 && r>r1) continue;

                                if (r>r1) image->set_value_at(i,j,0,v2*(r-r1)/(r2-r1)+v1*(r2-r)/(r2-r1));

                                else image->set_value_at(i,j,0,v1);

                        }

                }

        }

        else {

                float r;

                for (int k=(z<r2?0:z-r2); k<(z+r2>nz?nz:z+r2); k++) {

                        for (int j=(y<r2?0:y-r2); j<(y+r2>ny?ny:y+r2); j++) {

                                for (int i=(x<r2?0:x-r2); i<(x+r2>nx?nx:x+r2); i++) {

                                r=sqrt(float(Util::square(i-x)+Util::square(j-y)+Util::square(k-z)));

                                if (r>r2 && r>r1) continue;

                                if (r>r1) image->set_value_at(i,j,k,v2*(r-r1)/(r2-r1)+v1*(r2-r)/(r2-r1));

                                else image->set_value_at(i,j,k,v1);

                                }

                        }

                }

        }

        image->update();

}


void MaskAzProcessor::process_inplace(EMData *image) {

        if (image->is_complex()) throw ImageFormatException("MaskAzProcessor: target image must be real");


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        float phicen = params.set_default("phicen",0.0f)*M_PI/180.0;

        phicen=Util::angle_norm_pi(phicen);

        float phirange = params.set_default("phirange",180.0f)*M_PI/180.0;

        float phitrirange = params.set_default("phitrirange",0.0f)*M_PI/180.0;

        int phitriangle = params.set_default("phitriangle",0);

        float cx = params.set_default("cx",nx/2);

        float cy = params.set_default("cy",ny/2);

        float zmin = params.set_default("zmin",0);

        float zmax = params.set_default("zmax",nz);

        float ztri = params.set_default("ztriangle",0.0f);

        float inner_radius = params.set_default("inner_radius",0.0f);

        float outer_radius = params.set_default("outer_radius",nx+ny);


        for (int x=0; x<nx; x++) {

                for (int y=0; y<ny; y++) {

                        float az=atan2(y-cy,x-cx);

                        float az2=az-M_PI*2.0f;

                        float az3=az+M_PI*2.0f;

                        float r=hypot(y-cy,x-cx);

                        float val=0.0f;

                        if (r>inner_radius&&r<=outer_radius) {

                                float as=Util::angle_sub_2pi(az,phicen);

                                if (as<phirange) val=1.0f;

                                else if (!phitriangle || as>phirange+phitrirange) val=0.0f;

                                else val=1.0-(as-phirange)/float(phitrirange);

                        }

                        if (r==0 && inner_radius<=0)  val=1.0;


                        for (int z=0; z<nz; z++) {

                                if (z<zmin-ztri || z>zmax+ztri) image->mult_value_at_fast(x,y,z,0);

                                else if (z>=zmin+ztri && z<=zmax-ztri) image->mult_value_at_fast(x,y,z,val);

                                else if (z>=zmin-ztri && z<=zmin+ztri) image->mult_value_at_fast(x,y,z,val*((z-zmin)/(2.0f*ztri)+0.5));

                                else image->mult_value_at_fast(x,y,z,val*((zmax-z)/(2.0f*ztri)+0.5));

                        }

                }

        }


}


void CircularMaskProcessor::calc_locals(EMData *)

{

        xc = Util::fast_floor(nx/2.0f) + dx;

        yc = Util::fast_floor(ny/2.0f) + dy;

        zc = Util::fast_floor(nz/2.0f) + dz;


        if (outer_radius < 0) {

                outer_radius = nx / 2 + outer_radius +1;

                outer_radius_square = outer_radius * outer_radius;

        }


        if (inner_radius <= 0) {

                inner_radius_square = 0;

        }

}


void MaskEdgeMeanProcessor::calc_locals(EMData * image)

{

        if (!image) {

                throw NullPointerException("NULL image");

        }

        int nitems = 0;

        float sum = 0;

        float *data = image->get_data();

        size_t i = 0;


        xc = Util::fast_floor(nx/2.0f) + dx;

        yc = Util::fast_floor(ny/2.0f) + dy;

        zc = Util::fast_floor(nz/2.0f) + dz;


        for (int z = 0; z < nz; ++z) {

                for (int y = 0; y < ny; ++y) {

                        for (int x = 0; x < nx; ++x) {

                                float x1 = sqrt((x - xc) * (x - xc) + (y - yc) * (y - yc) + (z - zc) * (z - zc));

                                if (x1 <= outer_radius + ring_width && x1 >= outer_radius - ring_width) {

                                        sum += data[i];

                                        ++nitems;

                                }

                                ++i;

                        }

                }

        }


        ring_avg = sum / nitems;

}


void ToMinvalProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float minval = params.set_default("minval",0.0f);

        float newval = params.set_default("newval",minval);


        size_t size = (size_t)image->get_xsize() *

                                  (size_t)image->get_ysize() *

                                  (size_t)image->get_zsize();

        float *data = image->get_data();


        for (size_t i = 0; i < size; ++i) {

                if (data[i]<minval) data[i]=newval;

        }

        image->update();

}


void ComplexPixelProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL image");

                return;

        }

        if (!image->is_complex()) {

                LOGWARN("cannot apply complex processor on a real image. Nothing is done.");

                return;

        }


        size_t size = (size_t)image->get_xsize() *

                                  (size_t)image->get_ysize() *

                                  (size_t)image->get_zsize();

        float *data = image->get_data();


        image->ri2ap();


        for (size_t i = 0; i < size; i += 2) {

                process_pixel(data);

                data += 2;

        }


        image->update();

        image->ap2ri();

}


void AreaProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float *data = image->get_data();


        nx = image->get_xsize();

        ny = image->get_ysize();

        nz = image->get_zsize();


        int n = (areasize - 1) / 2;

        matrix_size = areasize * areasize;


        if (nz > 1) {

                matrix_size *= areasize;

        }


        float *matrix = new float[matrix_size];

        kernel = new float[matrix_size];


        size_t cpysize = areasize * sizeof(float);

        size_t start = (nx * ny + nx + 1) * n;


        int xend = nx - n;

        int yend = ny - n;


        int zstart = n;

        int zend = nz - n;


        int zbox_start = 0;

        int zbox_end = areasize;


        if (nz == 1) {

                zstart = 0;

                zend = 1;

                zbox_end = 1;

        }


        size_t nsec = (size_t)nx * (size_t)ny;

        int box_nsec = areasize * areasize;


        create_kernel();


        size_t total_size = (size_t)nx * (size_t)ny * (size_t)nz;

        float *data2 = new float[total_size];

        memcpy(data2, data, total_size * sizeof(float));


        size_t k;

        for (int z = zstart; z < zend; z++) {

                for (int y = n; y < yend; y++) {

                        for (int x = n; x < xend; x++) {


                                k = (size_t)z * nsec + y * nx + x;


                                for (int bz = zbox_start; bz < zbox_end; bz++) {

                                        for (int by = 0; by < areasize; by++) {

                                                memcpy(&matrix[(size_t)bz * box_nsec + by * areasize],

                                                           &data2[k - start + bz * nsec + by * nx], cpysize);

                                        }

                                }


                                process_pixel(&data[k], (float) x, (float) y, (float) z, matrix);

                        }

                }

        }


        if( matrix )

        {

                delete[]matrix;

                matrix = 0;

        }


        if( kernel )

        {

                delete[]kernel;

                kernel = 0;

        }

        image->update();

}


void LaplacianProcessor::process_inplace(EMData* image)

{

        EMData* d = new EMData();

        d = image->process("math.edge.xgradient");

        d->process_inplace("math.edge.xgradient");

        image->process_inplace("math.edge.ygradient");

        image->process_inplace("math.edge.xgradient");

        image->add(*d);

        image->process_inplace("normalize");

        delete d;

}


void LaplacianProcessor::create_kernel() const

{

        if (nz == 1) {

                memset(kernel, 0, areasize * areasize);

                kernel[1] = -0.25f;

                kernel[3] = -0.25f;

                kernel[5] = -0.25f;

                kernel[7] = -0.25f;

                kernel[4] = 1;

        }

        else {

                memset(kernel, 0, (size_t)areasize * areasize * areasize);

                kernel[4] = -1.0f / 6.0f;

                kernel[10] = -1.0f / 6.0f;

                kernel[12] = -1.0f / 6.0f;

                kernel[14] = -1.0f / 6.0f;

                kernel[16] = -1.0f / 6.0f;

                kernel[22] = -1.0f / 6.0f;

                kernel[13] = 1;

        }

}


void SetBitsProcessor::process_inplace(EMData *image)

{

        if (!image) return;


        float nsigma = params.set_default("nsigma",3.0f);

        int bits = params.set_default("bits",5);

        int floatbits = params.set_default("floatbits",-1);

        float sigma = image->get_attr("sigma");

        float mean = image->get_attr("mean");

        float bitmax=pow(2.,bits);


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();

        size_t total_size = (size_t)nx * (size_t)ny * (size_t)nz;


        if (floatbits<=0) {

                float min=image->get_attr("minimum");

                float max=image->get_attr("maximum");


                // our max/min are either the actual max/min, or mean+-nsigma*sigma

                // whichever produces the smaller range. That way by specifying a large

                // nsigma you can get the actual min/max of the image easily

                if (mean-nsigma*sigma>min) min=mean-nsigma*sigma;

                if (mean+nsigma*sigma<max) max=mean+nsigma*sigma;


                for (size_t i=0; i<total_size; i++) {

                        float newval=floor((image->get_value_at_index(i)-min)*bitmax/(max-min)


                        );

                        if (newval<0) newval=0;

                        if (newval>=bitmax) newval=bitmax-1.0f;

                        image->set_value_at_index(i,newval);

                }

        }

        else {

                // In this version, we keep a specified number of significant bits in the floating point significand, leaving the exponent and sign alone

                if (floatbits>22) throw InvalidValueException(floatbits,"floatbits must be <23");

                unsigned int bitmask= 0xffffffff << (23-floatbits); // this will leave 1 in the sign and exponent bits, along with the specified number of significand bits

                for (size_t i=0; i<total_size; i++) {

                        float newval=image->get_value_at_index(i);

                        unsigned int *newvali = (unsigned int*)&newval;

                        *newvali&=bitmask;

                        image->set_value_at_index(i,newval);

                }

        }

}


void BoxStatProcessor::process_inplace(EMData *image) {

        EMData *cpy = process(image);

        memcpy(image->get_data(),cpy->get_data(),image->get_size()*sizeof(float));

        delete cpy;

}


// mirrors coordinates at box edges

inline int MIRE(int x,int nx) { return x<0?-x:(x>=nx?nx-(x-nx+1):x); }


EMData *BoxStatProcessor::process(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return NULL;

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        int dx=1,dy=1,dz=1;

        if (params.has_key("radius")) dx=dy=dz=params["radius"];

        if (params.has_key("xsize")) dx=params["xsize"];

        if (params.has_key("ysize")) dy=params["ysize"];

        if (params.has_key("zsize")) dz=params["zsize"];

        if (nz==1) dz=0;


        int matrix_size = (2*dx+1)*(2*dy+1)*(2*dz+1);


        vector<float> array(matrix_size);

//      image->process_inplace("normalize");


        EMData *ret = image->copy_head();


        // The old version of this code had a lot of hand optimizations, which likely weren't accomplishing much

        // This is much simpler, but relies on the compiler to optimize. May be some cost associated with the new

        // edge-mirroring policy which could be hand optomized if necessary

        for (int k=0; k<nz; k++) {

                for (int j=0; j<ny; j++) {

                        for (int i=0; i<nx; i++) {


                                for (int kk=k-dz,s=0; kk<=k+dz; kk++) {

                                        for (int jj=j-dy; jj<=j+dy; jj++) {

                                                for (int ii=i-dx; ii<=i+dx; ii++,s++) {

                                                        array[s]=image->get_value_at(MIRE(ii,nx),MIRE(jj,ny),MIRE(kk,nz));

                                                }

                                        }

                                }

                                float newv=image->get_value_at(i,j,k);

                                process_pixel(&newv,array.data(),matrix_size);

                                ret->set_value_at(i,j,k,newv);

                        }

                }

        }


        ret->update();


        return ret;

}


void DiffBlockProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        int nz = image->get_zsize();


        if (nz > 1) {

                LOGERR("%s Processor doesn't support 3D", get_name().c_str());

                throw ImageDimensionException("3D model not supported");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();


        int v1 = params["cal_half_width"];

        int v2 = params["fill_half_width"];


        int v0 = v1 > v2 ? v1 : v2;


        if (v2 <= 0) {

                v2 = v1;

        }


        float *data = image->get_data();


        for (int y = v0; y <= ny - v0 - 1; y += v2) {

                for (int x = v0; x <= nx - v0 - 1; x += v2) {


                        float sum = 0;

                        for (int y1 = y - v1; y1 <= y + v1; y1++) {

                                for (int x1 = x - v1; x1 <= x + v1; x1++) {

                                        sum += data[x1 + y1 * nx];

                                }

                        }

                        float mean = sum / ((v1 * 2 + 1) * (v1 * 2 + 1));


                        for (int j = y - v2; j <= y + v2; j++) {

                                for (int i = x - v2; i <= x + v2; i++) {

                                        data[i + j * nx] = mean;

                                }

                        }

                }

        }


        image->update();

}


void CutoffBlockProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }

        int nz = image->get_zsize();


        if (nz > 1) {

                LOGERR("%s Processor doesn't support 3D", get_name().c_str());

                throw ImageDimensionException("3D model not supported");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();


        float value1 = params["value1"];

        float value2 = params["value2"];


        int v1 = (int) value1;

        int v2 = (int) value2;

        if (v2 > v1 / 2) {

                LOGERR("invalid value2 '%f' in CutoffBlockProcessor", value2);

                return;

        }


        if (v2 <= 0) {

                v2 = v1;

        }


        float *data = image->get_data();

        int y = 0, x = 0;

        for (y = 0; y <= ny - v1; y += v1) {

                for (x = 0; x <= nx - v1; x += v1) {


                        EMData *clip = image->get_clip(Region(x, y, v1, v1));

                        EMData *fft = clip->do_fft();


                        float *fft_data = fft->get_data();

                        float sum = 0;

                        int nitems = 0;


                        for (int i = -v2; i < v2; i++) {

                                for (int j = 0; j < v2; j++) {

                                        if (j == 0 && i == 0) {

                                                continue;

                                        }


                                        if (hypot(j, i) < value2) {

                                                int t = j * 2 + (i + v1 / 2) * (v1 + 2);

                                                sum += (fft_data[t] * fft_data[t] + fft_data[t + 1] * fft_data[t + 1]);

                                                nitems++;

                                        }

                                }

                        }


                        if( clip )

                        {

                                delete clip;

                                clip = 0;

                        }


                        float mean = sum / nitems;


                        for (int i = y; i < y + v1; i++) {

                                for (int j = x; j < x + v1; j++) {

                                        data[i * nx + j] = mean;

                                }

                        }

                }

        }


        memset(&data[y * nx], 0, (ny - y) * nx * sizeof(float));


        for (int i = 0; i < ny; i++) {

                memset(&data[i * nx + x], 0, (nx - x) * sizeof(float));

        }


        image->update();

}


void MedianShrinkProcessor::process_inplace(EMData * image)

{

        if (image->is_complex()) throw ImageFormatException("Error, the median shrink processor does not work on complex images");


        int shrink_factor =      params.set_default("n",0);

        if (shrink_factor <= 1) {

                throw InvalidValueException(shrink_factor,

                                                                        "median shrink: shrink factor must > 1");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


//      if ((nx % shrink_factor != 0) || (ny % shrink_factor != 0) || (nz > 1 && (nz % shrink_factor != 0))) {

//              throw InvalidValueException(shrink_factor, "Image size not divisible by shrink factor");

//      }


        int shrunken_nx = nx / shrink_factor;

        int shrunken_ny = ny / shrink_factor;

        int shrunken_nz = 1;

        if (nz > 1) shrunken_nz = nz / shrink_factor;


        EMData* copy = image->copy();

        image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);

        accrue_median(image,copy,shrink_factor);

        image->update();

        if( copy )

        {

                delete copy;

                copy = 0;

        }

}


//

EMData* MedianShrinkProcessor::process(const EMData *const image)

{

        if (image->is_complex()) throw ImageFormatException("Error, the median shrink processor does not work on complex images");


        int shrink_factor =      params.set_default("n",0);

        if (shrink_factor <= 1) {

                throw InvalidValueException(shrink_factor,

                                                                        "median shrink: shrink factor must > 1");

        }

        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


//      if ((nx % shrink_factor != 0) || (ny % shrink_factor != 0) || (nz > 1 && (nz % shrink_factor != 0))) {

//              throw InvalidValueException(shrink_factor, "Image size not divisible by shrink factor");

//      }


        int shrunken_nx = nx / shrink_factor;

        int shrunken_ny = ny / shrink_factor;

        int shrunken_nz = 1;

        if (nz > 1) shrunken_nz = nz / shrink_factor;


//      EMData* ret = new EMData(shrunken_nx, shrunken_ny, shrunken_nz);

        EMData *ret = image->copy_head();

        ret->set_size(shrunken_nx, shrunken_ny, shrunken_nz);


        accrue_median(ret,image,shrink_factor);

        ret->update();

        return ret;

}


void MedianShrinkProcessor::accrue_median(EMData* to, const EMData* const from,const int shrink_factor)

{


        int nx_old = from->get_xsize();

        int ny_old = from->get_ysize();


        int threed_shrink_factor = shrink_factor * shrink_factor;

        int z_shrink_factor = 1;

        if (from->get_zsize() > 1) {

                threed_shrink_factor *= shrink_factor;

                z_shrink_factor = shrink_factor;

        }


        float *mbuf = new float[threed_shrink_factor];


        int nxy_old = nx_old * ny_old;


        int nx = to->get_xsize();

        int ny = to->get_ysize();

        int nz = to->get_zsize();

        int nxy_new = nx * ny;


        float * rdata = to->get_data();

        const float *const data_copy = from->get_const_data();


        for (int l = 0; l < nz; l++) {

                int l_min = l * shrink_factor;

                int l_max = l * shrink_factor + z_shrink_factor;

                size_t cur_l = (size_t)l * nxy_new;


                for (int j = 0; j < ny; j++) {

                        int j_min = j * shrink_factor;

                        int j_max = (j + 1) * shrink_factor;

                        size_t cur_j = j * nx + cur_l;


                        for (int i = 0; i < nx; i++) {

                                int i_min = i * shrink_factor;

                                int i_max = (i + 1) * shrink_factor;


                                size_t k = 0;

                                for (int l2 = l_min; l2 < l_max; l2++) {

                                        size_t cur_l2 = l2 * nxy_old;


                                        for (int j2 = j_min; j2 < j_max; j2++) {

                                                size_t cur_j2 = j2 * nx_old + cur_l2;


                                                for (int i2 = i_min; i2 < i_max; i2++) {

                                                        mbuf[k] = data_copy[i2 + cur_j2];

                                                        ++k;

                                                }

                                        }

                                }


                                for (k = 0; k < size_t(threed_shrink_factor / 2 + 1); k++) {

                                        for (int i2 = k + 1; i2 < threed_shrink_factor; i2++) {

                                                if (mbuf[i2] < mbuf[k]) {

                                                        float f = mbuf[i2];

                                                        mbuf[i2] = mbuf[k];

                                                        mbuf[k] = f;

                                                }

                                        }

                                }


                                rdata[i + cur_j] = mbuf[threed_shrink_factor / 2];

                        }

                }

        }


        if( mbuf )

        {

                delete[]mbuf;

                mbuf = 0;

        }


        to->scale_pixel((float)shrink_factor);

}


EMData* FFTResampleProcessor::process(const EMData *const image)

{

        if (image->is_complex()) throw ImageFormatException("FFTResampleProcessor does not support complex images.");


        float sample_rate = params.set_default("n",0.0f);

        if (sample_rate <= 0.0F  )      {

                throw InvalidValueException(sample_rate,        "downsampling must be >0 ");

        }


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        int nnx=(int)floor(nx/sample_rate+0.5f);

        int nny=(ny==1?1:(int)floor(ny/sample_rate+0.5f));

        int nnz=(nz==1?1:(int)floor(nz/sample_rate+0.5f));


        if (nnx==nx && nny==ny && nnz==nz) return image->copy();


        EMData *result;

        //downsample

        if (nnx<nx||nny<ny||nnz<nz) {

                // the type casting here is because FourTruncate was not defined to be const (but it is)

                result=((EMData *)image)->FourTruncate(nnx, nny, nnz, 1, 0);    // nnx,nny,nnz,returnreal,normalize

        } //upscale

        else {

                // the type casting here is because FourTruncate was not defined to be const (but it is)

                result=((EMData *)image)->FourInterpol(nnx, nny, nnz, 1, 0);    // nnx,nny,nnz,returnreal,normalize

        }

        result->scale_pixel((float)nx/(float)nnx);

        result->update();

        return result;


//      EMData* result;

//      if (image->is_complex()) result = image->copy();

//      else result = image->do_fft();

//      fft_resample(result,image,sample_rate);

        // The image may have been padded - we should shift it so that the phase origin is where FFTW expects it

//      result->update();

//      result->scale_pixel(sample_rate);

//      return result;

}


void FFTResampleProcessor::process_inplace(EMData * image)

{

        float sample_rate = params.set_default("n",0.0f);

        if (sample_rate <= 0.0F  )      {

                throw InvalidValueException(sample_rate,        "sample rate (n) must be >0 ");

        }

        //if (image->is_complex()) throw ImageFormatException("Error, the fft resampling processor does not work on complex images");


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        int nnx=(int)floor(nx/sample_rate+0.5f);

        int nny=(ny==1?1:(int)floor(ny/sample_rate+0.5f));

        int nnz=(nz==1?1:(int)floor(nz/sample_rate+0.5f));


        if (nnx==nx && nny==ny && nnz==nz) return;


        EMData *result;

        //downsample

        if (nnx<nx||nny<ny||nnz<nz) {

                // the type casting here is because FourTruncate was not defined to be const (but it is)

                result=((EMData *)image)->FourTruncate(nnx, nny, nnz, 1, 0);    // nnx,nny,nnz,returnreal,normalize

        } //upscale

        else {

                // the type casting here is because FourTruncate was not defined to be const (but it is)

                result=((EMData *)image)->FourInterpol(nnx, nny, nnz, 1, 0);    // nnx,nny,nnz,returnreal,normalize

        }


        image->set_size(nnx,nny,nnz);

        memcpy(image->get_data(),result->get_data(),nnx*nny*nnz*sizeof(float));

        image->scale_pixel((float)nx/(float)nnx);

        image->update();

        delete result;


        //fft_resample(image,image,sample_rate);


        //image->scale_pixel(sample_rate);


}


void FFTResampleProcessor::fft_resample(EMData* to, const EMData *const from, const float& sample_rate) {

        int nx = from->get_xsize();

        int ny = from->get_ysize();

        int nz = from->get_zsize();


        int new_nx = static_cast<int>( static_cast<float> (nx) / sample_rate);

        int new_ny = static_cast<int>( static_cast<float> (ny) / sample_rate);

        int new_nz = static_cast<int>( static_cast<float> (nz) / sample_rate);


        if (new_nx == 0) throw UnexpectedBehaviorException("The resample rate causes the pixel dimensions in the x direction to go to zero");

        if (new_ny == 0) new_ny = 1;

        if (new_nz == 0) new_nz = 1;


        int ndim = from->get_ndim();

        if ( ndim < 3 ) {

                new_nz = 1;

        }

        if ( ndim < 2 ) {

                new_ny = 1;

        }


        int fft_x_correction = 1;

        if (new_nx % 2 == 0) fft_x_correction = 2;


        int fft_y_correction = 0;

        if (ny != 1 && new_ny % 2 == 0 && ny % 2 == 1) fft_y_correction = 1;

        else if (ny != 1 && new_ny % 2 == 1 && ny % 2 == 0) fft_y_correction = -1;


        int fft_z_correction = 0;

        if (nz != 1 && new_nz % 2 == 0 && nz % 2 == 1) fft_z_correction = 1;

        else if (nz != 1 && new_nz % 2 == 1 && nz % 2 == 0) fft_z_correction = -1;


        if ( ! to->is_complex()) to->do_fft_inplace();


        if (ndim != 1) to->process_inplace("xform.fourierorigin.tocenter");


        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);

        to->clip_inplace(clip);


        if (fft_x_correction == 1) to->set_fftodd(true);

        else to->set_fftodd(false);


        if (ndim != 1) to->process_inplace("xform.fourierorigin.tocorner");


        to->do_ift_inplace();

        to->depad_corner();


}


EMData* MeanShrinkProcessor::process(const EMData *const image)

{

        if (image->is_complex()) throw ImageFormatException("Error, the mean shrink processor does not work on complex images");


        if (image->get_ndim() == 1) { throw ImageDimensionException("Error, mean shrink works only for 2D & 3D images"); }


        float shrink_factor0 = params.set_default("n",0.0f);

        int shrink_factor = int(shrink_factor0);

        if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {

                throw InvalidValueException(shrink_factor0,

                                                                        "mean shrink: shrink factor must be >1 integer or 1.5");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        // here handle the special averaging by 1.5 for 2D case

        if (shrink_factor0==1.5 ) {

                if (nz > 1 ) throw InvalidValueException(shrink_factor0, "mean shrink: only support 2D images for shrink factor = 1.5");


                int shrunken_nx = (int(nx / 1.5)+1)/2*2;        // make sure the output size is even

                int shrunken_ny = (int(ny / 1.5)+1)/2*2;

                EMData* result = new EMData(shrunken_nx,shrunken_ny,1);


                accrue_mean_one_p_five(result,image);

                result->update();


                return result;

        }


        int shrunken_nx = nx / shrink_factor;

        int shrunken_ny = ny / shrink_factor;

        int shrunken_nz = 1;


        if (nz > 1) {

                shrunken_nz = nz / shrink_factor;

        }


//      EMData* result = new EMData(shrunken_nx,shrunken_ny,shrunken_nz);

        EMData* result = image->copy_head();

        result->set_size(shrunken_nx,shrunken_ny,shrunken_nz);

        accrue_mean(result,image,shrink_factor);


        result->update();


        return result;

}


void MeanShrinkProcessor::process_inplace(EMData * image)

{

        if (image->is_complex()) throw ImageFormatException("Error, the mean shrink processor does not work on complex images");


        if (image->get_ndim() == 1) { throw ImageDimensionException("Error, mean shrink works only for 2D & 3D images"); }


        float shrink_factor0 = params.set_default("n",0.0f);

        int shrink_factor = int(shrink_factor0);

        if (shrink_factor0 <= 1.0F || ((shrink_factor0 != shrink_factor) && (shrink_factor0 != 1.5F) ) ) {

                throw InvalidValueException(shrink_factor0,

                                                                        "mean shrink: shrink factor must be >1 integer or 1.5");

        }


/*      if ((nx % shrink_factor != 0) || (ny % shrink_factor != 0) ||

        (nz > 1 && (nz % shrink_factor != 0))) {

        throw InvalidValueException(shrink_factor,

        "Image size not divisible by shrink factor");

}*/


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();

        // here handle the special averaging by 1.5 for 2D case

        if (shrink_factor0==1.5 ) {

                if (nz > 1 ) throw InvalidValueException(shrink_factor0, "mean shrink: only support 2D images for shrink factor = 1.5");


                int shrunken_nx = (int(nx / 1.5)+1)/2*2;        // make sure the output size is even

                int shrunken_ny = (int(ny / 1.5)+1)/2*2;


                EMData* orig = image->copy();

                image->set_size(shrunken_nx, shrunken_ny, 1);   // now nx = shrunken_nx, ny = shrunken_ny

                image->to_zero();


                accrue_mean_one_p_five(image,orig);


                if( orig ) {

                        delete orig;

                        orig = 0;

                }

                image->update();


                return;

        }


        accrue_mean(image,image,shrink_factor);


        int shrunken_nx = nx / shrink_factor;

        int shrunken_ny = ny / shrink_factor;

        int shrunken_nz = 1;

        if (nz > 1) shrunken_nz = nz / shrink_factor;


        image->update();

        image->set_size(shrunken_nx, shrunken_ny, shrunken_nz);

}


void MeanShrinkProcessor::accrue_mean(EMData* to, const EMData* const from,const int shrink_factor)

{

        const float * const data = from->get_const_data();

        float* rdata = to->get_data();


        size_t nx = from->get_xsize();

        size_t ny = from->get_ysize();

        size_t nz = from->get_zsize();

        size_t nxy = nx*ny;


        size_t shrunken_nx = nx / shrink_factor;

        size_t shrunken_ny = ny / shrink_factor;

        size_t shrunken_nz = 1;

        size_t shrunken_nxy = shrunken_nx * shrunken_ny;


        int normalize_shrink_factor = shrink_factor * shrink_factor;

        int z_shrink_factor = 1;


        if (nz > 1) {

                shrunken_nz = nz / shrink_factor;

                normalize_shrink_factor *= shrink_factor;

                z_shrink_factor = shrink_factor;

        }


        float invnormfactor = 1.0f/(float)normalize_shrink_factor;


        for (size_t k = 0; k < shrunken_nz; k++) {

                size_t k_min = k * shrink_factor;

                size_t k_max = k * shrink_factor + z_shrink_factor;

                size_t cur_k = k * shrunken_nxy;


                for (size_t j = 0; j < shrunken_ny; j++) {

                        size_t j_min = j * shrink_factor;

                        size_t j_max = j * shrink_factor + shrink_factor;

                        size_t cur_j = j * shrunken_nx + cur_k;


                        for (size_t i = 0; i < shrunken_nx; i++) {

                                size_t i_min = i * shrink_factor;

                                size_t i_max = i * shrink_factor + shrink_factor;


                                float sum = 0;

                                for (size_t kk = k_min; kk < k_max; kk++) {

                                        size_t cur_kk = kk * nxy;


                                        for (size_t jj = j_min; jj < j_max; jj++) {

                                                size_t cur_jj = jj * nx + cur_kk;

                                                for (size_t ii = i_min; ii < i_max; ii++) {

                                                        sum += data[ii + cur_jj];

                                                }

                                        }

                                }

                                rdata[i + cur_j] = sum * invnormfactor;

                        }

                }

        }

        to->scale_pixel((float)shrink_factor);

}


void MeanShrinkProcessor::accrue_mean_one_p_five(EMData* to, const EMData * const from)

{

        int nx0 = from->get_xsize(), ny0 = from->get_ysize();   // the original size


        int nx = to->get_xsize(), ny = to->get_ysize();


        float *data = to->get_data();

        const float * const data0 = from->get_const_data();


        for (int j = 0; j < ny; j++) {

                int jj = int(j * 1.5);

                float jw0 = 1.0F, jw1 = 0.5F;   // 3x3 -> 2x2, so each new pixel should have 2.25 of the old pixels

                if ( j%2 ) {

                        jw0 = 0.5F;

                        jw1 = 1.0F;

                }

                for (int i = 0; i < nx; i++) {

                        int ii = int(i * 1.5);

                        float iw0 = 1.0F, iw1 = 0.5F;

                        float w = 0.0F;


                        if ( i%2 ) {

                                iw0 = 0.5F;

                                iw1 = 1.0F;

                        }

                        if ( jj < ny0 ) {

                                if ( ii < nx0 ) {

                                        data[j * nx + i] = data0[ jj * nx0 + ii ] * jw0 * iw0 ;

                                        w += jw0 * iw0 ;

                                        if ( ii+1 < nx0 ) {

                                                data[j * nx + i] += data0[ jj * nx0 + ii + 1] * jw0 * iw1;

                                                w += jw0 * iw1;

                                        }

                                }

                                if ( jj +1 < ny0 ) {

                                        if ( ii < nx0 ) {

                                                data[j * nx + i] += data0[ (jj+1) * nx0 + ii ] * jw1 * iw0;

                                                w += jw1 * iw0;

                                                if ( ii+1 < nx0 ) {

                                                        data[j * nx + i] += data0[ (jj+1) * nx0 + ii + 1] * jw1 * iw1;

                                                        w += jw1 * iw1;

                                                }

                                        }

                                }

                        }

                        if ( w>0 ) data[j * nx + i] /= w;

                }

        }


        to->update();

        to->scale_pixel((float)1.5);

}


// This would have to be moved into the header if it were required in other source files

template<class LogicOp>

EMData* BooleanShrinkProcessor::process(const EMData *const image, Dict& params)

{

        // The basic idea of this code is to iterate through each pixel in the output image

        // determining its value by investigation a region of the input image


        if (!image) throw NullPointerException("Attempt to max shrink a null image");


        if (image->is_complex() ) throw ImageFormatException("Can not max shrink a complex image");


        int shrink = params.set_default("n",2);

        int search = params.set_default("search",2);


        if ( shrink < 0 ) throw InvalidValueException(shrink, "Can not shrink by a value less than 0");


        int nz = image->get_zsize();

        int ny = image->get_ysize();

        int nx = image->get_xsize();


        if (nx == 1 && ny == 1 && nz == 1 ) return image->copy();


        LogicOp op;

        EMData* return_image = new EMData();


        int shrinkx = shrink;

        int shrinky = shrink;

        int shrinkz = shrink;


        int searchx = search;

        int searchy = search;

        int searchz = search;


        // Clamping the shrink values to the dimension lengths

        // ensures that the return image has non zero dimensions

        if ( shrinkx > nx ) shrinkx = nx;

        if ( shrinky > ny ) shrinky = ny;

        if ( shrinkz > nz ) shrinkz = nz;


        if ( nz == 1 && ny == 1 )

        {

                return_image->set_size(nx/shrinkx);

                for(int i = 0; i < nx/shrinkx; ++i)

                {

                        float tmp = op.get_start_val();

                        for(int s=0; s < searchx; ++s)

                        {

                                int idx = shrinkx*i+s;

                                // Don't ask for memory beyond limits

                                if ( idx > nx ) break;

                                else

                                {

                                        float val = image->get_value_at(idx);

                                        if ( op( val,tmp) ) tmp = val;

                                }

                        }

                        return_image->set_value_at(i,tmp);

                }

        }

        else if ( nz == 1 )

        {

                int ty = ny/shrinky;

                int tx = nx/shrinkx;

                return_image->set_size(tx,ty);

                for(int y = 0; y < ty; ++y) {

                        for(int x = 0; x < tx; ++x) {

                                float tmp = op.get_start_val();

                                for(int sy=0; sy < searchy; ++sy) {

                                        int yidx = shrinky*y+sy;

                                        if ( yidx >= ny) break;

                                        for(int sx=0; sx < searchx; ++sx) {

                                                int xidx = shrinkx*x+sx;

                                                if ( xidx >= nx) break;


                                                float val = image->get_value_at(xidx,yidx);

                                                if ( op( val,tmp) ) tmp = val;

                                        }

                                }

                                return_image->set_value_at(x,y,tmp);

                        }

                }

        }

        else

        {

                int tz = nz/shrinkz;

                int ty = ny/shrinky;

                int tx = nx/shrinkx;


                return_image->set_size(tx,ty,tz);

                for(int z = 0; z < tz; ++z) {

                        for(int y = 0; y < ty; ++y) {

                                for(int x = 0; x < tx; ++x) {

                                        float tmp = op.get_start_val();


                                        for(int sz=0; sz < searchz; ++sz) {

                                                int zidx = shrinkz*z+sz;

                                                if ( zidx >= nz) break;


                                                for(int sy=0; sy < searchy; ++sy) {

                                                        int yidx = shrinky*y+sy;

                                                        if ( yidx >= ny) break;


                                                        for(int sx=0; sx < searchx; ++sx) {

                                                                int xidx = shrinkx*x+sx;

                                                                if ( xidx >= nx) break;

                                                                float val = image->get_value_at(xidx,yidx,zidx);

                                                                if ( op( val,tmp) ) tmp = val;

                                                        }

                                                }

                                        }

                                        return_image->set_value_at(x,y,z,tmp);

                                }

                        }

                }

        }

        return_image->update();


        return return_image;

}


template<class LogicOp>

void BooleanShrinkProcessor::process_inplace(EMData * image, Dict& params)

{

        // The basic idea of this code is to iterate through each pixel in the output image

        // determining its value by investigation a region of the input image

        if (!image) throw NullPointerException("Attempt to max shrink a null image");


        if (image->is_complex() ) throw ImageFormatException("Can not max shrink a complex image");


        int shrink = params.set_default("shrink",2);

        int search = params.set_default("search",2);


        if ( shrink < 0 ) throw InvalidValueException(shrink, "Can not shrink by a value less than 0");


        int nz = image->get_zsize();

        int ny = image->get_ysize();

        int nx = image->get_xsize();


        LogicOp op;


        int shrinkx = shrink;

        int shrinky = shrink;

        int shrinkz = shrink;


        int searchx = search;

        int searchy = search;

        int searchz = search;


        // Clamping the shrink values to the dimension lengths

        // ensures that the return image has non zero dimensions

        if ( shrinkx > nx ) shrinkx = nx;

        if ( shrinky > ny ) shrinky = ny;

        if ( shrinkz > nz ) shrinkz = nz;


        if (nx == 1 && ny == 1 && nz == 1 ) return;


        if ( nz == 1 && ny == 1 )

        {

                for(int i = 0; i < nx/shrink; ++i)

                {

                        float tmp = op.get_start_val();

                        for(int s=0; s < searchx; ++s)

                        {

                                int idx = shrinkx*i+s;

                                if ( idx > nx ) break;

                                else

                                {

                                        float val = image->get_value_at(idx);

                                        if ( op( val,tmp) ) tmp = val;

                                }

                        }

                        image->set_value_at(i,tmp);

                }


                image->set_size(nx/shrinkx);

        }

        else if ( nz == 1 )

        {

                int ty = ny/shrinky;

                int tx = nx/shrinkx;

                for(int y = 0; y < ty; ++y) {

                        for(int x = 0; x < tx; ++x) {

                                float tmp = op.get_start_val();

                                for(int sy=0; sy < searchy; ++sy) {

                                        int yidx = shrinky*y+sy;

                                        if ( yidx >= ny) break;

                                        for(int sx=0; sx < searchx; ++sx) {

                                                int xidx = shrinkx*x+sx;

                                                if ( xidx >= nx) break;


                                                float val = image->get_value_at(xidx,yidx);

                                                if ( op( val,tmp) ) tmp = val;

                                        }

                                }

                                (*image)(x+tx*y) = tmp;

                        }

                }

                image->set_size(tx,ty);

        }

        else

        {

                int tnxy = nx/shrinkx*ny/shrinky;

                int tz = nz/shrinkz;

                int ty = ny/shrinky;

                int tx = nx/shrinkx;


                for(int z = 0; z < tz; ++z) {

                        for(int y = 0; y < ty; ++y) {

                                for(int x = 0; x < tx; ++x) {

                                        float tmp = op.get_start_val();

                                        for(int sz=0; sz < searchz; ++sz) {

                                                int zidx = shrinkz*z+sz;

                                                if ( zidx >= nz) break;

                                                for(int sy=0; sy < searchy; ++sy) {

                                                        int yidx = shrinky*y+sy;

                                                        if ( yidx >= ny) break;

                                                        for(int sx=0; sx < shrinkx; ++sx) {

                                                                int xidx = shrinkx*x+sx;

                                                                if ( xidx >= nx) break;


                                                                float val = image->get_value_at(xidx,yidx,zidx);

                                                                if ( op( val,tmp) ) tmp = val;

                                                        }

                                                }

                                        }

                                        (*image)(x+tx*y+tnxy*z) = tmp;

                                }

                        }

                }

                image->set_size(tx,ty,tz);

        }

        image->update();

}


void GradientRemoverProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        int nz = image->get_zsize();

        if (nz > 1) {

                LOGERR("%s Processor doesn't support 3D model", get_name().c_str());

                throw ImageDimensionException("3D model not supported");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        float *dy = new float[ny];

        float m = 0;

        float b = 0;

        float sum_y = 0;

        float *data = image->get_data();


        for (int i = 0; i < ny; i++) {

                Util::calc_least_square_fit(nx, 0, data + i * nx, &m, &b, false);

                dy[i] = b;

                sum_y += m;

        }


        float mean_y = sum_y / ny;

        float sum_x = 0;

        Util::calc_least_square_fit(ny, 0, dy, &sum_x, &b, false);


        for (int j = 0; j < ny; j++) {

                for (int i = 0; i < nx; i++) {

                        data[i + j * nx] -= i * sum_x + j * mean_y + b;

                }

        }


        image->update();

}


void FlattenBackgroundProcessor::process_inplace(EMData * image)

{


        EMData* mask = params.set_default("mask",(EMData*)0);

        int radius = params.set_default("radius",0);


        if (radius != 0 && mask != 0) throw InvalidParameterException("Error - the mask and radius parameters are mutually exclusive.");


        if (mask == 0 && radius == 0) throw InvalidParameterException("Error - you must specify either the mask or the radius parameter.");


        // If the radius isn't 0, then turn the mask into the thing we want...

        bool deletemask = false;

        if (radius != 0) {

                mask = new EMData;

                int n = image->get_ndim();

                if (n==1){

                        mask->set_size(2*radius+1);

                } else if (n==2) {

                        mask->set_size(2*radius+1,2*radius+1);

                }

                else /*n==3*/ {

                        mask->set_size(2*radius+1,2*radius+1,2*radius+1);

                }

                // assuming default behavior is to make a circle/sphere with using the radius of the mask

                mask->process_inplace("testimage.circlesphere");

        }


        // Double check that that mask isn't too big

        int mnx = mask->get_xsize(); int mny = mask->get_ysize(); int mnz = mask->get_zsize();

        int nx = image->get_xsize(); int ny = image->get_ysize(); int nz = image->get_zsize();

        int nxc = nx+mnx; int nyc = ny+mny; int nzc = nz+mnz;

        if (nz == 1) nzc = 1; // Sanity check

        if (ny == 1) nyc = 1; // Sanity check


        if ( mnx > nx || mny > ny || mnz > nz)

                throw ImageDimensionException("Can not flatten using a mask that is larger than the image.");


        // Get the normalization factor

        float normfac = 0.0;

        for (int i=0; i<mask->get_xsize()*mask->get_ysize()*mask->get_zsize(); ++i){

                normfac += mask->get_value_at(i);

        }

        // If the sum is zero the user probably doesn't understand that determining a measure of the mean requires

        // strictly positive numbers. The user has specified a mask that consists entirely of zeros, or the mask

        // has a mean of zero.

        if (normfac == 0) throw InvalidParameterException("Error - the pixels in the mask sum to zero. This breaks the flattening procedure");

        normfac = 1.0f/normfac;


        // The mask can now be automatically resized to the dimensions of the image

//      bool undoclip = false;


        Region r;

        if (ny == 1) r = Region((mnx-nxc)/2,nxc);

        else if (nz == 1) r = Region((mnx-nxc)/2, (mny-nyc)/2,nxc,nyc);

        else r = Region((mnx-nxc)/2, (mny-nyc)/2,(mnz-nzc)/2,nxc,nyc,nzc);

        mask->clip_inplace(r,0);

//      undoclip = true;

//      if ( mnx < nx || mny < ny || mnz < nz) {

//              Region r((mnx-nx)/2, (mny-ny)/2,(mnz-nz)/2,nx,ny,nz);

//              mask->clip_inplace(r);

//              undoclip = true;

//      }


        Region r2;

        if (ny == 1) r2 = Region((nx-nxc)/2,nxc);

        else if (nz == 1) r2 = Region((nx-nxc)/2, (ny-nyc)/2,nxc,nyc);

        else r2 = Region((nx-nxc)/2, (ny-nyc)/2,(nz-nzc)/2,nxc,nyc,nzc);

        image->clip_inplace(r2,image->get_edge_mean());

        // Finally do the convolution

        EMData* m = image->convolute(mask);

        // Normalize so that m is truly the local mean

        m->mult(normfac);

        // Before we can subtract, the mean must be phase shifted

        m->process_inplace("xform.phaseorigin.tocenter");

        // Subtract the local mean

//      image->write_image("a.mrc");

//      m->write_image("b.mrc");

        image->sub(*m); // WE'RE DONE!

        delete m;


        if (deletemask) {

                delete mask;

        } else { // I clipped it inplace, so undo this clipping so the user gets back what the put in

                Region r;

                if (ny == 1) r = Region((nxc-mnx)/2,mnx);

                else if (nz == 1) r = Region((nxc-mnx)/2, (nyc-mny)/2,mnx,mny);

                else r = Region((nxc-mnx)/2, (nyc-mny)/2,(nzc-mnz)/2,mnx,mny,mnz);

                mask->clip_inplace(r);

        }


        Region r3;

        if (ny == 1) r3 = Region((nxc-nx)/2,nx);

        else if (nz == 1) r3 = Region((nxc-nx)/2, (nyc-ny)/2,nx,ny);

        else r3 = Region((nxc-nx)/2, (nyc-ny)/2,(nzc-nz)/2,nx,ny,nz);

        image->clip_inplace(r3);

//      if ( undoclip ) {

//              Region r((nx-mnx)/2, (ny-mny)/2, (nz-mnz)/2,mnx,mny,mnz);

//              mask->clip_inplace(r);

//      }


}


void NonConvexProcessor::process_inplace(EMData * image) {

        if (!image) { LOGWARN("NULL IMAGE"); return; }

        //int isinten=image->get_attr_default("is_intensity",0);


        // 1-D

        if (image->get_ysize()==1) {


        }

        // 2-D

        else if (image->get_zsize()==1) {

//              if (!isinten) throw ImageDimensionException("Only complex intensity images currently supported by NonConvexProcessor");

                int nx2=image->get_xsize()/2;

                int ny2=image->get_ysize()/2;

                vector<float> rdist = image->calc_radial_dist(nx2*1.5,0,1,false);               // radial distribution to make sure nonconvex values decrease radially

                // Make sure rdist is decreasing (or flat)

                for (int i=1; i<nx2; i++) {

                        if (rdist[i]>rdist[i-1]) rdist[i]=rdist[i-1];

                }


                image->process_inplace("xform.fourierorigin.tocenter");

                EMData* binary=image->copy();


                // First we eliminate convex points from the input image (set to zero)

                for (int x=0; x<image->get_xsize(); x+=2) {

                        for (int y=1; y<image->get_ysize()-1; y++) {

                                int r=(int)hypot((float)(x/2),(float)(y-ny2));

                                float cen=(*binary)(x,y);

                                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

                                        image->set_value_at_fast(x/2+nx2,y,0.0);        // we are turning image into a full real-space intensity image for now

                                        image->set_value_at_fast(nx2-x/2,ny2*2-y-1,0.0);

                                }

                                else {

                                        image->set_value_at_fast(x/2+nx2,y,cen);        // we are turning image into a full real-space intensity image for now

                                        image->set_value_at_fast(nx2-x/2,ny2*2-y-1,cen);        // It will contain non-zero values only for nonconvex points

                                }

                        }

                }

                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

                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

                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

                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

                for (int y=0; y<ny2*2; y++) image->set_value_at_fast(0,y,0.0f);


                // Now make a binary version of the convex points

                float *idat=image->get_data();

                float *bdat=binary->get_data();

                int nxy=(nx2*ny2*4);

                for (int i=0; i<nxy; i++) {

                        bdat[i]=idat[i]==0?0:1.0f;              // binary version of the convex points in image

                }

                binary->update();


                // We now use a Gaussian filter on both images, to use Gaussian interpolation to fill in zero values

                image->set_complex(false);              // so we can use a Gaussian filter on it

                binary->set_complex(false);


/*              image->write_image("con.hdf",0);*/

                image->set_fftpad(false);

                binary->set_fftpad(false);


                // Gaussian blur of both images

                image->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.04f));

                binary->process_inplace("filter.lowpass.gauss",Dict("cutoff_abs",0.04f));


/*              image->write_image("con.hdf",1);

                binary->write_image("con.hdf",2);*/


                for (int x=0; x<image->get_xsize(); x+=2) {

                        for (int y=0; y<image->get_ysize(); y++) {

                                float bv=binary->get_value_at(x/2+nx2,y);

                                image->set_value_at_fast(x,y,image->get_value_at(x/2+nx2,y)/(bv<=0?1.0f:bv));

                                image->set_value_at_fast(x+1,y,0.0);

                        }

                }

                image->set_complex(true);

                image->set_fftpad(true);

                image->process_inplace("xform.fourierorigin.tocorner");

                delete binary;

        }

        else throw ImageDimensionException("3D maps not yet supported by NonConvexProcessor");


}


#include <gsl/gsl_linalg.h>

void GradientPlaneRemoverProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        int nz = image->get_zsize();

        if (nz > 1) {

                LOGERR("%s Processor doesn't support 3D model", get_name().c_str());

                throw ImageDimensionException("3D map not supported");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        float *d = image->get_data();

        EMData *mask = 0;

        float *dm = 0;

        if (params.has_key("mask")) {

                mask = params["mask"];

                if (nx!=mask->get_xsize() || ny!=mask->get_ysize()) {

                        LOGERR("%s Processor requires same size mask image", get_name().c_str());

                        throw ImageDimensionException("wrong size mask image");

                }

                dm = mask->get_data();

        }

        int count = 0;

        if (dm) {

                for(int i=0; i<nx*ny; i++) {

                        if(dm[i]) count++;

                }

        }

        else {

                count = nx * ny;

        }

        if(count<3) {

                LOGERR("%s Processor requires at least 3 pixels to fit a plane", get_name().c_str());

                throw ImageDimensionException("too few usable pixels to fit a plane");

        }

        // Allocate the working space

        gsl_vector *S=gsl_vector_calloc(3);

        gsl_matrix *A=gsl_matrix_calloc(count,3);

        gsl_matrix *V=gsl_matrix_calloc(3,3);


        double m[3] = {0, 0, 0};

        int index=0;

        if (dm) {

                for(int j=0; j<ny; j++){

                        for(int i=0; i<nx; i++){

                                int ij=j*nx+i;

                                if(dm[ij]) {

                                        m[0]+=i;        // x

                                        m[1]+=j;        // y

                                        m[2]+=d[ij];    // z

                                        /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \

                                                        index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/

                                        index++;

                                }

                        }

                }

        }

        else {

                for(int j=0; j<ny; j++){

                        for(int i=0; i<nx; i++){

                                int ij=j*nx+i;

                                        m[0]+=i;        // x

                                        m[1]+=j;        // y

                                        m[2]+=d[ij];    // z

                                        /*printf("index=%d/%d\ti,j=%d,%d\tval=%g\txm,ym,zm=%g,%g,%g\n", \

                                                        index,count,i,j,d[ij],m[0]/(index+1),m[1]/(index+1),m[2]/(index+1));*/

                                        index++;

                        }

                }

        }


        for(int i=0; i<3; i++) m[i]/=count; // compute center of the plane


        index=0;

        if (dm) {

                for(int j=0; j<ny; j++){

                        for(int i=0; i<nx; i++){

                                int ij=j*nx+i;

                                if(dm[ij]) {

                                        //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]);

                                        gsl_matrix_set(A,index,0,i-m[0]);

                                        gsl_matrix_set(A,index,1,j-m[1]);

                                        gsl_matrix_set(A,index,2,d[ij]-m[2]);

                                        index++;

                                }

                        }

                }

                mask->update();

        }

        else {

                for(int j=0; j<ny; j++){

                        for(int i=0; i<nx; i++){

                                int ij=j*nx+i;

                                        //printf("index=%d/%d\ti,j=%d,%d\tval=%g\n",index,count,i,j,d[index]);

                                        gsl_matrix_set(A,index,0,i-m[0]);

                                        gsl_matrix_set(A,index,1,j-m[1]);

                                        gsl_matrix_set(A,index,2,d[ij]-m[2]);

                                        index++;

                        }

                }

        }


        // SVD decomposition and use the V vector associated with smallest singular value as the plan normal

        gsl_linalg_SV_decomp_jacobi(A, V, S);


        double n[3];

        for(int i=0; i<3; i++) n[i] = gsl_matrix_get(V, i, 2);


        #ifdef DEBUG

        printf("S=%g,%g,%g\n",gsl_vector_get(S,0), gsl_vector_get(S,1), gsl_vector_get(S,2));

        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));

        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));

        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));

        printf("Fitted plane: p0=%g,%g,%g\tn=%g,%g,%g\n",m[0],m[1],m[2],n[0],n[1],n[2]);

        #endif


        int changeZero = 0;

        if (params.has_key("changeZero")) changeZero = params["changeZero"];

        if (changeZero) {

                for(int j=0; j<nx; j++){

                        for(int i=0; i<ny; i++){

                                int ij = j*nx+i;

                                d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);

                        }

                }

        }

        else {

                for(int j=0; j<nx; j++){

                        for(int i=0; i<ny; i++){

                                int ij = j*nx+i;

                                if(d[ij]) d[ij]-=static_cast<float>(-((i-m[0])*n[0]+(j-m[1])*n[1])/n[2]+m[2]);

                        }

                }

        }

        image->update();

        // set return plane parameters

        vector< float > planeParam;

        planeParam.resize(6);

        for(int i=0; i<3; i++) planeParam[i] = static_cast<float>(n[i]);

        for(int i=0; i<3; i++) planeParam[i+3] = static_cast<float>(m[i]);

        params["planeParam"]=EMObject(planeParam);

}


void VerticalStripeProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        float *data = image->get_data();

        float sigma = image->get_attr("sigma");


        for (int k = 0; k < nz; k++) {

                for (int i = 0; i < nx; i++) {

                        double sum = 0;

                        for (int j = ny / 4; j < 3 * ny / 4; j++) {

                                sum += data[i + j * nx];

                        }


                        float mean = (float)sum / (ny / 2);

                        for (int j = 0; j < ny; j++) {

                                data[i + j * nx] = (data[i + j * nx] - mean) / sigma;

                        }

                }

        }


        image->update();

}


void RealToFFTProcessor::process_inplace(EMData *image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        //Note : real image only!

        if(image->is_complex()) {

                LOGERR("%s Processor only operates on real images", get_name().c_str());

                throw ImageFormatException("apply to real image only");

        }


        // Note : 2D only!

        int nz = image->get_zsize();

        if (nz > 1) {

                LOGERR("%s Processor doesn't support 3D models", get_name().c_str());

                throw ImageDimensionException("3D model not supported");

        }


        EMData *ff=image->do_fft();

        ff->ri2ap();


        int nx=image->get_xsize();

        int ny=image->get_ysize();


        int x,y;

        float norm=static_cast<float>(nx*ny);


        for (y=0; y<ny; y++) image->set_value_at(0,y,0);


        for (x=1; x<nx/2; x++) {

                for (y=0; y<ny; y++) {

                        int y2;

                        if (y<ny/2) y2=y+ny/2;

                        else if (y==ny/2) y2=ny;

                        else y2=y-ny/2;

                        image->set_value_at(x,y,ff->get_value_at(nx-x*2,ny-y2)/norm);

                }

        }


        for (x=nx/2; x<nx; x++) {

                for (y=0; y<ny; y++) {

                        int y2;

                        if (y<ny/2) y2=y+ny/2;

                        else y2=y-ny/2;

                        image->set_value_at(x,y,ff->get_value_at(x*2-nx,y2)/norm);

                }

        }


        image->update();

        if( ff )

        {

                delete ff;

                ff = 0;

        }

}


void SigmaZeroEdgeProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        if (image->get_zsize() > 1) {

                LOGERR("%s Processor doesn't support 3D model", get_name().c_str());

                throw ImageDimensionException("3D model not supported");

        }


        float nonzero = params.set_default("nonzero",false);


        float *d = image->get_data();

        int i = 0;

        int j = 0;


        int nx = image->get_xsize();

        int ny = image->get_ysize();


        float zval=9.99e23f;            // we're assuming we won't run into this exact value for an edge, not great programming, but good enough

        if (nonzero) {

                int x,y;

                size_t corn=nx*ny-1;


                // this set of 4 tests looks for any edges with exactly the same value

                for (x=1; x<nx; x++) { if (d[x]!=d[0]) break;}

                if (x==nx) zval=d[0];


                for (y=1; y<ny; y++) { if (d[y*nx]!=d[0]) break; }

                if (y==ny) zval=d[0];


                for (x=1; x<nx; x++) { if (d[corn-x]!=d[corn]) break;}

                if (x==nx) zval=d[corn];


                for (y=1; y<ny; y++) { if (d[corn-y*nx]!=d[corn]) break; }

                if (y==ny) zval=d[corn];


                if (zval!=9.99e23f) {

                        image->set_attr("hadzeroedge",1);

//                      printf("zeroedge %f\n",zval);


                }

                else image->set_attr("hadzeroedge",0);


                // 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.

                for (x=nx/2-5; x<nx/2+5; x++) {

                        if (d[x]!=d[x+nx] || d[x]!=d[x+nx*2] ) break;

                }

                if (x==nx/2+5) image->set_attr("hadzeroedge",2);


                for (x=nx/2-5; x<nx/2+5; x++) {

                        if (d[corn-x]!=d[corn-x-nx] || d[corn-x]!=d[corn-x-nx*2]) break;

                }

                if (x==nx/2+5) image->set_attr("hadzeroedge",2);


                for (y=ny/2-5; y<ny/2+5; y++) {

                        if (d[y*nx]!=d[y*nx+1] || d[y*nx]!=d[y*nx+2] ) break;

                }

                if (y==ny/2+5) image->set_attr("hadzeroedge",2);


                for (y=ny/2-5; y<ny/2+5; y++) {

                        if (d[corn-y*nx]!=d[corn-y*nx-1] || d[corn-y*nx]!=d[corn-y*nx-2]) break;

                }

                if (y==ny/2+5) image->set_attr("hadzeroedge",2);


        }

        if (zval==9.99e23f) zval=0;


        for (j = 0; j < ny; j++) {

                for (i = 0; i < nx - 1; i++) {

                        if (d[i + j * nx] != zval) {

                                break;

                        }

                }


                float v = d[i + j * nx];

                while (i >= 0) {

                        d[i + j * nx] = v;

                        i--;

                }


                for (i = nx - 1; i > 0; i--) {

                        if (d[i + j * nx] != zval)

                                break;

                }

                v = d[i + j * nx];

                while (i < nx) {

                        d[i + j * nx] = v;

                        i++;

                }

        }


        for (i = 0; i < nx; i++) {

                for (j = 0; j < ny; j++) {

                        if (d[i + j * nx] != zval)

                                break;

                }


                float v = d[i + j * nx];

                while (j >= 0) {

                        d[i + j * nx] = v;

                        j--;

                }


                for (j = ny - 1; j > 0; j--) {

                        if (d[i + j * nx] != zval)

                                break;

                }

                v = d[i + j * nx];

                while (j < ny) {

                        d[i + j * nx] = v;

                        j++;

                }

        }


        image->update();

}


void OutlierProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        bool fix_zero=(bool)params.set_default("fix_zero",0);

        float sigmamult=(float)params.set_default("sigma",3.0);


        if (sigmamult<=0.0) throw InvalidValueException(sigmamult,"threshold.outlier.localmean: sigma must be >0");


        float hithr=(float)image->get_attr("mean")+(float)(image->get_attr("sigma"))*sigmamult;

        float lothr=(float)image->get_attr("mean")-(float)(image->get_attr("sigma"))*sigmamult;


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        // This isn't optimally efficient

        EMData *im[2];

        im[0]=image;

        im[1]=image->copy_head();


        if (nz==1) {

                int repeat=1;

                while (repeat) {

                        memcpy(im[1]->get_data(),im[0]->get_data(),image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));

                        repeat=0;

                        for (int y=0; y<ny; y++) {

                                for (int x=0; x<nx; x++) {

                                        // if the pixel is an outlier

                                        float pix=im[1]->get_value_at(x,y);

                                        if (pix>hithr || pix<lothr || (pix==0 && fix_zero)) {

                                                int y0=0>y-1?0:y-1;

                                                int y1=y>=ny-1?ny-1:y+1;

                                                int x0=0>x-1?0:x-1;

                                                int x1=x>=nx-1?nx-1:x+1;

                                                float c=0.0f,nc=0.0f;

                                                for (int yy=y0; yy<=y1; yy++) {

                                                        for (int xx=x0; xx<=x1; xx++) {

                                                                float lpix=im[1]->get_value_at(xx,yy);

                                                                if (lpix>hithr || lpix<lothr || (lpix==0 && fix_zero)) continue;

                                                                c+=lpix;

                                                                nc++;

                                                        }

                                                }

                                                if (nc!=0) im[0]->set_value_at(x,y,c/nc);

                                                else repeat=1;

                                        }

                                }

                        }

                }

        }

        else {

                throw ImageDimensionException("threshold.outlier.localmean: 3D not yet implemented");

        }

        delete im[1];

}


void BeamstopProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }

        if (image->get_zsize() > 1) {

                LOGERR("BeamstopProcessor doesn't support 3D model");

                throw ImageDimensionException("3D model not supported");

        }


        float value1 = params["value1"];

        float value2 = params["value2"];

        float value3 = params["value3"];


        float thr = fabs(value1);

        float *data = image->get_data();

        int cenx = (int) value2;

        int ceny = (int) value3;


        int nx = image->get_xsize();

        int ny = image->get_ysize();


        if (cenx <= 0) {

                cenx = nx / 2;

        }


        if (ceny <= 0) {

                ceny = ny / 2;

        }


        int mxr = (int) floor(sqrt(2.0f) * nx / 2);


        float *mean_values = new float[mxr];

        float *sigma_values = new float[mxr];

        double sum = 0;

        int count = 0;

        double square_sum = 0;


        for (int i = 0; i < mxr; i++) {

                sum = 0;

                count = 0;

                square_sum = 0;

                int nitems = 6 * i + 2;


                for (int j = 0; j < nitems; j++) {

                        float ang = j * 2 * M_PI / nitems;

                        int x0 = (int) floor(cos(ang) * i + cenx);

                        int y0 = (int) floor(sin(ang) * i + ceny);


                        if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny) {

                                continue;

                        }


                        float f = data[x0 + y0 * nx];

                        sum += f;

                        square_sum += f * f;

                        count++;

                }


                mean_values[i] = (float)sum / count;

                sigma_values[i] = (float) sqrt(square_sum / count - mean_values[i] * mean_values[i]);

        }


        for (int k = 0; k < 5; k++) {

                for (int i = 0; i < mxr; i++) {

                        sum = 0;

                        count = 0;

                        square_sum = 0;

                        int nitems = 6 * i + 2;

                        double thr1 = mean_values[i] - sigma_values[i] * thr;

                        double thr2 = mean_values[i] + sigma_values[i];


                        for (int j = 0; j < nitems; j++) {

                                float ang = j * 2 * M_PI / nitems;

                                int x0 = (int) floor(cos(ang) * i + cenx);

                                int y0 = (int) floor(sin(ang) * i + ceny);


                                if (x0 < 0 || y0 < 0 || x0 >= nx || y0 >= ny ||

                                        data[x0 + y0 * nx] < thr1 || data[x0 + y0 * nx] > thr2) {

                                        continue;

                                }


                                sum += data[x0 + y0 * nx];

                                square_sum += data[x0 + y0 * nx] * data[x0 + y0 * nx];

                                count++;

                        }


                        mean_values[i] = (float) sum / count;

                        sigma_values[i] = (float) sqrt(square_sum / count - mean_values[i] * mean_values[i]);

                }

        }


        for (int i = 0; i < nx; i++) {

                for (int j = 0; j < ny; j++) {


                        int r = Util::round(hypot((float) i - cenx, (float) j - ceny));


                        if (value1 < 0) {

                                if (data[i + j * nx] < (mean_values[r] - sigma_values[r] * thr)) {

                                        data[i + j * nx] = 0;

                                }

                                else {

                                        data[i + j * nx] -= mean_values[r];

                                }

                                continue;

                        }

                        if (data[i + j * nx] > (mean_values[r] - sigma_values[r] * thr)) {

                                continue;

                        }

                        data[i + j * nx] = mean_values[r];

                }

        }


        if( mean_values )

        {

                delete[]mean_values;

                mean_values = 0;

        }


        if( sigma_values )

        {

                delete[]sigma_values;

                sigma_values = 0;

        }


        image->update();

}


void MeanZeroEdgeProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }

        if (image->get_zsize() > 1) {

                LOGERR("MeanZeroEdgeProcessor doesn't support 3D model");

                throw ImageDimensionException("3D model not supported");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        Dict dict = image->get_attr_dict();

        float mean_nonzero = dict.get("mean_nonzero");


        float *d = image->get_data();

        int i = 0;

        int j = 0;


        for (j = 0; j < ny; j++) {

                for (i = 0; i < nx - 1; i++) {

                        if (d[i + j * nx] != 0) {

                                break;

                        }

                }


                if (i == nx - 1) {

                        i = -1;

                }


                float v = d[i + j * nx] - mean_nonzero;


                while (i >= 0) {

                        v *= 0.9f;

                        d[i + j * nx] = v + mean_nonzero;

                        i--;

                }


                for (i = nx - 1; i > 0; i--) {

                        if (d[i + j * nx] != 0) {

                                break;

                        }

                }


                if (i == 0) {

                        i = nx;

                }


                v = d[i + j * nx] - mean_nonzero;


                while (i < nx) {

                        v *= .9f;

                        d[i + j * nx] = v + mean_nonzero;

                        i++;

                }

        }


        for (i = 0; i < nx; i++) {

                for (j = 0; j < ny; j++) {

                        if (d[i + j * nx] != 0)

                                break;

                }


                float v = d[i + j * nx] - mean_nonzero;


                while (j >= 0) {

                        v *= .9f;

                        d[i + j * nx] = v + mean_nonzero;

                        j--;

                }


                for (j = ny - 1; j > 0; j--) {

                        if (d[i + j * nx] != 0)

                                break;

                }


                v = d[i + j * nx] - mean_nonzero;


                while (j < ny) {

                        v *= .9f;

                        d[i + j * nx] = v + mean_nonzero;

                        j++;

                }

        }


        image->update();

}


void AverageXProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float *data = image->get_data();

        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();

        size_t nxy = (size_t)nx * ny;


        size_t idx;

        for (int z = 0; z < nz; z++) {

                for (int x = 0; x < nx; x++) {

                        double sum = 0;

                        for (int y = 0; y < ny; y++) {

                                idx = x + y * nx + z * nxy;

                                sum += data[idx];

                        }

                        float mean = (float) sum / ny;


                        for (int y = 0; y < ny; y++) {

                                idx = x + y * nx + z * nxy;

                                data[idx] = mean;

                        }

                }

        }


        image->update();

}


void DecayEdgeProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();


        float *d = image->get_data();

        int width = params["width"];


        if (width > min(nx,ny)/2.){

                LOGERR("width parameter cannot be greater than min(nx,ny)/2");

                throw InvalidParameterException("width cannot be greater than min(nx,ny)/2");

        }


        if (image->get_zsize() > 1){

                for (int k=0; k<image->get_zsize(); k++){

                        int zidx = k*nx*ny;

                        for (int i=0; i<width; i++) {

                                float frac=i/(float)width;

                                for (int j=0; j<nx; j++) {

                                        d[zidx+j+i*nx]*=frac;

                                        d[zidx+nx*ny-j-i*nx-1]*=frac;

                                }

                                for (int j=0; j<ny; j++) {

                                        d[zidx+j*nx+i]*=frac;

                                        d[zidx+nx*ny-j*nx-i-1]*=frac;

                                }

                        }

                }

        }

        else {


                for (int i=0; i<width; i++) {

                        float frac=i/(float)width;

                        for (int j=0; j<nx; j++) {

                                d[j+i*nx]*=frac;

                                d[nx*ny-j-i*nx-1]*=frac;

                        }

                        for (int j=0; j<ny; j++) {

                                d[j*nx+i]*=frac;

                                d[nx*ny-j*nx-i-1]*=frac;

                        }

                }

        }


        image->update();

}


void ZeroEdgeRowProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        if (image->get_zsize() > 1) {

                LOGERR("ZeroEdgeRowProcessor is not supported in 3D models");

                throw ImageDimensionException("3D model not supported");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();


        float *d = image->get_data();

        int apodize = params.set_default("apodize",0);

        if (apodize>0) throw InvalidValueException(apodize," apodize not yet supported.");


        int top_nrows = std::max(0,params.set_default("y0",0)-apodize);

        int bottom_nrows = std::max(0,params.set_default("y1",top_nrows)-apodize);


        int left_ncols = std::max(0,params.set_default("x0",0)-apodize);

        int right_ncols = std::max(0,params.set_default("x1",left_ncols)-apodize);


        size_t row_size = nx * sizeof(float);


        memset(d, 0, top_nrows * row_size);

        memset(d + (ny - bottom_nrows) * nx, 0, bottom_nrows * row_size);


        for (int i = top_nrows; i < ny - bottom_nrows; i++) {

                memset(d + i * nx, 0, left_ncols * sizeof(float));

                memset(d + i * nx + nx - right_ncols, 0, right_ncols * sizeof(float));

        }

        image->update();

}


void ZeroEdgePlaneProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        if (image->get_zsize() <= 1) {

                LOGERR("ZeroEdgePlaneProcessor only support 3D models");

                throw ImageDimensionException("3D model only");

        }


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        float *d = image->get_data();


        int x0=params["x0"];

        int x1=params.set_default("x1",x0);

        int y0=params["y0"];

        int y1=params.set_default("y1",y0);

        int z0=params["z0"];

        int z1=params.set_default("z1",z0);


        size_t row_size = nx * sizeof(float);

        size_t nxy = nx * ny;

        size_t sec_size = nxy * sizeof(float);

        size_t y0row = y0 * row_size;

        size_t y1row = y1 * row_size;

        int max_y = ny-y1;

        size_t x0size = x0*sizeof(float);

        size_t x1size = x1*sizeof(float);


        memset(d,0,z0*sec_size);                                        // zero -z

        memset(d+(nxy*(nz-z1)),0,sec_size*z1);          // zero +z


        for (int z=z0; z<nz-z1; z++) {

                memset(d+z*nxy,0,y0row);                        // zero -y

                memset(d+z*nxy+(ny-y1)*nx,0,y1row); // zero +y


                int znxy = z * nxy;

                int znxy2 = znxy + nx - x1;


                for (int y=y0; y<max_y; y++) {

                        memset(d+znxy+y*nx,0,x0size);   // zero -x

                        memset(d+znxy2+y*nx,0,x1size);  // zero +x

                }

        }


        image->update();

}


float NormalizeProcessor::calc_sigma(EMData * image) const

{

        return image->get_attr("sigma");

}


void NormalizeProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("cannot do normalization on NULL image");

                return;

        }


        if (image->is_complex()) {

                LOGWARN("cannot do normalization on complex image");

                return;

        }


        float sigma = calc_sigma(image);

        if (sigma == 0 || !Util::goodf(&sigma)) {

                LOGWARN("cannot do normalization on image with sigma = 0");

                return;

        }


        float mean = calc_mean(image);


        size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();

        float *data = image->get_data();


        for (size_t i = 0; i < size; ++i) {

                data[i] = (data[i] - mean) / sigma;

        }


        image->update();

}


float NormalizeUnitProcessor::calc_sigma(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        float ret=sqrt((float)image->get_attr("square_sum"));

        return ret==0.0f?1.0f:ret;

}


float NormalizeUnitSumProcessor::calc_sigma(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        float ret=(float)image->get_attr("mean")*image->get_xsize()*image->get_ysize()*image->get_zsize();

        return ret==0.0f?1.0f:ret;

}


void NormalizeMaskProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }

        EMData *mask = params["mask"];

        int no_sigma = params.set_default("no_sigma",0);

        int apply_mask = params.set_default("apply_mask",0);


        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        double sum=0;

        double sumsq=0;

        double nnz=0;


        for (int z=0; z<nz; z++) {

                for (int y=0; y<ny; y++) {

                        for (int x=0; x<nx; x++) {

                                double m=mask->get_value_at(x,y,z);

                                if (m==0) continue;

                                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

                                double v=(double)image->get_value_at(x,y,z)*m;

                                sum+=v;

                                sumsq+=v*v;

                                nnz+=1.0;

                        }

                }

        }


        float mean=(float)(sum/nnz);

        float sigma=sqrt((float)(sumsq-sum*sum/nnz)/nnz);

        if (no_sigma) sigma=1.0f;


        for (int z=0; z<nz; z++) {

                for (int y=0; y<ny; y++) {

                        for (int x=0; x<nx; x++) {

                                float m=mask->get_value_at(x,y,z);

                                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

                                if (m==0) image->set_value_at(x,y,z,0);

                                else image->set_value_at(x,y,z,(image->get_value_at(x,y,z)*m-mean)/sigma);

                        }

                }

        }

}


void NormalizeRampNormVar::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("cannot do normalization on NULL image");

                return;

        }


        if (image->is_complex()) {

                LOGWARN("cannot do normalization on complex image");

                return;

        }


        image->process_inplace( "filter.ramp" );

        int nx = image->get_xsize();

        EMData mask(nx,nx);

        mask.process_inplace("testimage.circlesphere", Dict("radius",nx/2-2,"fill",1));


        vector<float> rstls = Util::infomask( image, &mask, false);

        image->add((float)-rstls[0]);

        image->mult((float)1.0/rstls[1]);

        image->update();

}


void NormalizeByMassProcessor::process_inplace(EMData * image)

{

        float mass = params.set_default("mass",-1.0f);

        int verbose = params.set_default("verbose",0);


        if (mass <= 0) throw InvalidParameterException("You must specify a positive non zero mass");


        float tthr = params.set_default("thr",(float)image->get_attr("mean")+(float)image->get_attr("sigma"));


        float apix = image->get_attr_default("apix_x",1.0f);

        apix = params.set_default("apix",apix);


        if (apix <= 0) throw InvalidParameterException("You must specify a positive non zero apix");


        float step = ((float)image->get_attr("sigma"))/5.0f;


        if (step==0) throw InvalidParameterException("This image has sigma=0, cannot give it mass");


        size_t n = image->get_size();

        float* d = image->get_data();


        float thr=(float)image->get_attr("mean")+(float)image->get_attr("sigma")/2.0;

        int count=0;

        for (size_t i=0; i<n; ++i) {

                if (d[i]>=thr) ++count;

        }

        if (verbose) printf("apix=%1.3f\tmass=%1.1f\tthr=%1.2f\tstep=%1.3g\n",apix,mass,thr,step);


        float max = image->get_attr("maximum");

        float min = image->get_attr("minimum");

        for (int j=0; j<4; j++) {

                int err=0;

                while (thr<max && count*apix*apix*apix*.81/1000.0>mass) {

                        thr+=step;

                        count=0;

                        for (size_t i=0; i<n; ++i) {

                                if (d[i]>=thr) ++count;

                        }

                        err+=1;

                        if (err>1000) throw InvalidParameterException("Specified mass could not be achieved");

                        if (verbose>1) printf("%d\t%d\t%1.3f\t%1.2f\n",err,count,thr,count*apix*apix*apix*.81/1000.0);

                }


                step/=4.0;


                while (thr>min && count*apix*apix*apix*.81/1000.0<mass) {

                        thr-=step;

                        count=0;

                        for (size_t i=0; i<n; ++i) {

                                if (d[i]>=thr) ++count;

                        }

                        err+=1;

                        if (err>1000) throw InvalidParameterException("Specified mass could not be achieved");

                        if (verbose>1) printf("%d\t%d\t%1.3f\t%1.2f\n",err,count,thr,count*apix*apix*apix*.81/1000.0);


                }

                step/=4.0;

        }


        // We don't adjust the map if we get a negative value

        if ((float)tthr/thr>0) image->mult((float)tthr/thr);

        else printf("WARNING: could not normalize map to specified mass.");

        image->update();

}


float NormalizeEdgeMeanProcessor::calc_mean(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        return image->get_edge_mean();

}


float NormalizeCircleMeanProcessor::calc_mean(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

//      return image->get_circle_mean();

        int nx=image->get_xsize();

        int ny=image->get_ysize();

        int nz=image->get_zsize();


        float radius = params.set_default("radius",((float)ny/2-2));

        float width = params.set_default("width",2.0f);

        if (radius<0) radius=ny/2+radius;


        static bool busy = false;               // reduce problems with threads and different image sizes

        static EMData *mask = 0;

        static int oldradius=radius;

        static int oldwidth=width;


        if (!mask || !EMUtil::is_same_size(image, mask)||radius!=oldradius||width!=oldwidth) {

                while (busy) ;

                busy=true;

                if (!mask) {

                        mask = new EMData();

                }

                mask->set_size(nx, ny, nz);

                mask->to_one();


                mask->process_inplace("mask.sharp", Dict("inner_radius", radius,

                                                         "outer_radius", radius + width));


        }

        busy=true;

        double n = 0,s=0;

        float *d = mask->get_data();

        float * data = image->get_data();

        size_t size = (size_t)nx*ny*nz;

        for (size_t i = 0; i < size; ++i) {

                if (d[i]!=0.0f) { n+=1.0; s+=data[i]; }

        }


        float result = (float)(s/n);

//      printf("cmean=%f\n",result);

        busy=false;


        return result;

}


float NormalizeMaxMinProcessor::calc_sigma(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        float maxval = image->get_attr("maximum");

        float minval = image->get_attr("minimum");

        return (maxval + minval) / 2;

}


float NormalizeMaxMinProcessor::calc_mean(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        float maxval = image->get_attr("maximum");

        float minval = image->get_attr("minimum");

        return (maxval - minval) / 2;

}


float NormalizeLREdgeMeanProcessor::calc_mean(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        double sum = 0;

        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();

        float *d = image->get_data();

        size_t nyz = ny * nz;


        for (size_t i = 0; i < nyz; i++) {

                size_t l = i * nx;

                size_t r = l + nx - 2;

                sum += d[l] + d[l + 1] + d[r] + d[r + 1];

        }

        float mean = (float) sum / (4 * nyz);

        return mean;

}


void NormalizeRowProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float *rdata = image->get_data();

        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        if (nz > 1) {

                // Muyuan 08/2018: normalize by slice when nz>1

                // since this processor is not used anyways..

                for (int z = 0; z < nz; z++) {

                        double row_sum = 0;

                        double row_sqr = 0;

                        for (int y = 0; y < ny; y++) {

                                for (int x = 0; x < nx; x++) {

                                        float d = rdata[x + y * nx + z*nx*ny];

                                        row_sum += d;

                                        d*=d;

                                        row_sqr += d;

                                }

                        }

                        double row_mean = row_sum / nx / ny;

                        double row_std = row_sqr / nx / ny;

                        row_std = sqrt(row_std - row_mean*row_mean);

                        for (int y = 0; y < ny; y++) {

                                for (int x = 0; x < nx; x++) {

                                        rdata[x + y * nx+ z*nx*ny] -= (float)row_mean;

                                        rdata[x + y * nx+ z*nx*ny] /= (float)row_std;

                                }

                        }


                }


//              LOGERR("row normalize only works for 2D image");

//              return;

        }


        else {

                if (params.has_key("unitlen")) {

                        for (int y = 0; y < ny; y++) {

                                double row_len = 0;

                                for (int x = 0; x < nx; x++) {

                                        row_len += pow((double)rdata[x + y * nx],2.0);

                                }

                                row_len=sqrt(row_len);

                                if (row_len==0) row_len=1.0;


                                for (int x = 0; x < nx; x++) {

                                        rdata[x + y * nx] /= (float)row_len;

                                }

                        }

                        image->update();

                        return;

                }


                for (int y = 0; y < ny; y++) {

                        double row_sum = 0;

                        for (int x = 0; x < nx; x++) {

                                row_sum += rdata[x + y * nx];

                        }


                        double row_mean = row_sum / nx;

                        if (row_mean <= 0) {

                                row_mean = 1;

                        }


                        for (int x = 0; x < nx; x++) {

                                rdata[x + y * nx] /= (float)row_mean;

                        }

                }

        }

        image->update();

}


float NormalizeStdProcessor::calc_mean(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        return image->get_attr("mean");

}


float NormalizeHistPeakProcessor::calc_mean(EMData * image) const

{

        if (!image) {

                LOGWARN("NULL Image");

                return 0;

        }

        float mean = image->get_attr("mean");

        float sig = image->get_attr("sigma");

        size_t n=image->get_size();


        // We want at least 100 values per bin on average

        int hist_size = image->get_size()/100;

        if (hist_size>1000) hist_size=1000;

        if (hist_size<5) hist_size=5;

        vector<float> hist(hist_size,0.0f);


        float histmin=mean-2.0f*sig;

        float histmax=mean+2.0f*sig;

        float step=(histmax-histmin)/(hist_size-1);

        for (size_t i=0; i<n; i++) {

                int bin = (image->get_value_at_index(i)-histmin)/step;

                bin=bin<0?0:(bin>=hist_size?hist_size-1:bin);

                hist[bin]++;

        }


        // yes, probably could use fancy vector methods...

        int maxv=hist[1],maxn=1;

        for (int i=2; i<hist_size-1; i++) {

                if (hist[i]>maxv) { maxv=hist[i]; maxn=i; }

        }


        // weighted average position of the peak and nearest neighbors

        float v1=hist[maxn-1],v2=hist[maxn],v3=hist[maxn+1];

        float ret= (v1*(maxn-1)+v2*maxn+v3*(maxn+1))/(v1+v2+v3)*step+histmin;

        printf("%f %f %f %f %f %f\n",v1,v2,v3,mean,sig,ret);


        return ret;

}


EMData *SubtractOptProcessor::process(const EMData * const image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return NULL;

        }


        EMData *refr = params["ref"];

        EMData *actual = params.set_default("actual",(EMData*)NULL);

        EMData *ref;

        bool return_radial = params.set_default("return_radial",false);

        bool return_fft = params.set_default("return_fft",false);

        bool ctfweight = params.set_default("ctfweight",false);

        bool return_presigma = params.set_default("return_presigma",false);

        bool return_subim = params.set_default("return_subim",false);

        int si0=(int)floor(params.set_default("low_cutoff_frequency",0.0f)*image->get_ysize());

        int si1=(int)ceil(params.set_default("high_cutoff_frequency",0.7071f)*image->get_ysize());              // include the corners unless explicitly excluded


        // We will be modifying imf, so it needs to be a copy

        EMData *imf;

        if (image->is_complex()) imf=image->copy();

        else imf=image->do_fft();


        if (ctfweight) {

                EMData *ctfi=imf->copy_head();

                Ctf *ctf;

//              if (image->has_attr("ctf"))

                        ctf=(Ctf *)(image->get_attr("ctf"));

//              else ctf=(Ctf *)(ref->get_attr("ctf"));

                ctf->compute_2d_complex(ctfi,Ctf::CTF_INTEN);

                imf->mult(*ctfi);

                delete ctfi;

        }


        // Make sure ref is complex

        if (refr->is_complex()) ref=refr;

        else ref=refr->do_fft();


        EMData *actf;

        if (actual==NULL) actf=ref;

        else {

                if (ctfweight) throw InvalidCallException("math.sub.optimal: Sorry, cannot use ctfweight in combination with actual");

                if (actual->is_complex()) actf=actual;

                else actf=actual->do_fft();

        }


        int ny2=(int)(image->get_ysize()*sqrt(2.0)/2);

        vector <double>rad(ny2+1);

        vector <double>norm(ny2+1);


        // We are essentially computing an FSC here, but while the reference (the image

        // we plan to subtract) is normalized, the other image is not. This gives us a filter

        // to apply to the reference to optimally eliminate its contents from 'image'.

        for (int y=-ny2; y<ny2; y++) {

                for (int x=0; x<ny2; x++) {

                        int r=int(Util::hypot_fast(x,y));

                        if (r>ny2) continue;

                        std::complex<float> v1=imf->get_complex_at(x,y);

                        std::complex<float> v2=ref->get_complex_at(x,y);

                        rad[r]+=(double)(v1.real()*v2.real()+v1.imag()*v2.imag());

//                      norm[r]+=v2.real()*v2.real()+v2.imag()*v2.imag()+v1.real()*v1.real()+v1.imag()*v1.imag();

                        norm[r]+=(double)(v2.real()*v2.real()+v2.imag()*v2.imag());

                }

        }

        for (int i=1; i<ny2; i++) rad[i]/=norm[i];

        rad[0]=0;

//      FILE *out=fopen("dbug.txt","w");

//      for (int i=0; i<ny2; i++) fprintf(out,"%lf\t%lf\t%lf\n",(float)i,rad[i],norm[i]);

//      fclose(out);


        float oldsig=-1.0;

        // This option computes the real-space sigma on the input-image after the specified filter

        // This is an expensive option, but more efficient than computing the same using other means

        if (return_presigma) {

                for (int y=-ny2; y<ny2; y++) {

                        for (int x=0; x<imf->get_xsize()/2; x++) {

                                int r=int(Util::hypot_fast(x,y));

                                if (r>=ny2 || r>=si1 || r<si0) {

                                        imf->set_complex_at(x,y,0);

                                        continue;

                                }

                                std::complex<float> v1=imf->get_complex_at(x,y);

                                imf->set_complex_at(x,y,v1);

                        }

                }

                EMData *tmp=imf->do_ift();

                oldsig=(float)tmp->get_attr("sigma");

                delete tmp;

        }


        for (int y=-ny2; y<ny2; y++) {

                for (int x=0; x<imf->get_xsize()/2; x++) {

                        int r=int(Util::hypot_fast(x,y));

                        if (r>=ny2 || r>=si1 || r<si0) {

                                imf->set_complex_at(x,y,0);

                                continue;

                        }

                        std::complex<float> v1=imf->get_complex_at(x,y);

                        std::complex<float> v2=actf->get_complex_at(x,y);

                        v2*=(float)rad[r];

                        if (return_subim) imf->set_complex_at(x,y,v2);

                        else imf->set_complex_at(x,y,v1-v2);

                }

        }


        if (!refr->is_complex()) delete ref;

        if (actual!=NULL && !actual->is_complex()) delete actf;


        vector <float>radf;

        if (return_radial) {

                radf.resize(ny2);

                for (int i=0; i<ny2; i++) radf[i]=(float)rad[i];

        }


        if (!return_fft) {

                EMData *ret=imf->do_ift();

                delete imf;

                if (return_radial) ret->set_attr("filter_curve",radf);

                if (return_presigma) {

                        ret->set_attr("sigma_presub",oldsig);

//                      printf("Set %f\n",(float)image->get_attr("sigma_presub"));

                }

                return ret;

        }

        if (return_radial) imf->set_attr("filter_curve",radf);

        if (return_presigma) imf->set_attr("sigma_presub",oldsig);

        return imf;

}


void SubtractOptProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL image");

                return;

        }


        EMData *tmp=process(image);

        memcpy(image->get_data(),tmp->get_data(),(size_t)image->get_xsize()*image->get_ysize()*image->get_zsize()*sizeof(float));

        delete tmp;

        image->update();

        return;

}


void NormalizeToLeastSquareProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        EMData *to = params["to"];


        bool ignore_zero = params.set_default("ignore_zero",true);

        bool fourieramp = params.set_default("fourieramp",false);

        bool debug = params.set_default("debug",false);

        float ignore_lowsig = params.set_default("ignore_lowsig",-1.0);

        float low_threshold = params.set_default("low_threshold",-FLT_MAX);

        float high_threshold = params.set_default("high_threshold",FLT_MAX);


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();

        size_t size = (size_t)nx * ny * nz;

        size_t size2 = 0;


        EMData *fim = NULL;

        EMData *fto = NULL;

        float *dimage;          // unthresholded image data

        float *dto;

        int step=1;                     // used in Fourier mode


        float meani=(float)image->get_attr("mean");

        float meant=(float)to->get_attr("mean");

        float sigi=(float)image->get_attr("sigma")*ignore_lowsig;

        float sigt=(float)to->get_attr("sigma")*ignore_lowsig;


        if (fourieramp) {

                fim=image->do_fft();

                fto=to->do_fft();

                fim->ri2ap();

                fto->ri2ap();

                dimage=fim->get_data();

                dto=fto->get_data();

                step=2;

                size2=(size_t)(nx+2) * ny * nz;


                // sigma for thresholding

                meani=meant=0;                                          // Fourier amplitude >=0, so we reference everything to this point

                sigi=sigt=0;

                for (size_t i=0; i<size2; i+=2) { sigi+=pow(dimage[i],2.0f); sigt+=pow(dto[i],2.0f); }

                sigi=ignore_lowsig*sqrt(sigi/(size2/2));

                sigt=ignore_lowsig*sqrt(sigt/(size2/2));

        }

        else {

                dimage = image->get_data();

                dto = to->get_data();

                size2=size;

        }


        // rewrote this to just use GSL and get rid of David's old code.

        // The two passes are just to make sure we don't eat too much RAM if we do 3D

        if (ignore_lowsig<0) ignore_lowsig=0;


        size_t count=0;

        for (size_t i = 0; i < size2; i+=step) {

                if (dto[i] >= low_threshold && dto[i] <= high_threshold

                        && (dto[i]>=meant+sigt || dto[i]<=meant-sigt)

                        && (dimage[i]>=meani+sigi || dimage[i]<=meani-sigi)

                        && (!ignore_zero ||(dto[i] != 0.0f && dimage[i] != 0.0f))) {

                        count++;

                }

        }


        double *x=(double *)malloc(count*sizeof(double));

        double *y=(double *)malloc(count*sizeof(double));

        count=0;

        for (size_t i = 0; i < size2; i+=step) {

                if (dto[i] >= low_threshold && dto[i] <= high_threshold

                        && (dto[i]>=meant+sigt || dto[i]<=meant-sigt)

                        && (dimage[i]>=meani+sigi || dimage[i]<=meani-sigi)

                        && (!ignore_zero ||(dto[i] != 0.0f && dimage[i] != 0.0f))) {

                        x[count]=dimage[i];

                        y[count]=dto[i];

                        count++;

                }

        }

        double c0,c1;

        double cov00,cov01,cov11,sumsq;

        gsl_fit_linear (x, 1, y, 1, count, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);


        if (debug) {

                FILE*out=fopen("debug.txt","w");

                for (size_t i = 0; i < count; i++) {

                        fprintf(out,"%lf\t%lf\n",x[i],y[i]);

                }

                fclose(out);

                printf("add %lf\tmul %lf\t%lf\t%lf\t%lf\t%lf\n",c0,c1,cov00,cov01,cov11,sumsq);

        }


        free(x);

        free(y);

        if (fim!=NULL) delete fim;

        if (fto!=NULL) delete fto;


        if (fourieramp) c0=0;   // c0 makes no sense in this context. Really just a measure of noise

        dimage = image->get_data();

        for (size_t i = 0; i < size; ++i) dimage[i]=dimage[i]*c1+c0;

        image->update();

        image->set_attr("norm_mult",c1);

        image->set_attr("norm_add",c0);

}


void BinarizeFourierProcessor::process_inplace(EMData* image) {

        ENTERFUNC;

        if (!image->is_complex()) throw ImageFormatException("Fourier binary thresholding processor only works for complex images");


        float threshold = params.set_default("value",0.0f);

        image->ri2ap(); //      works for cuda


        float* d = image->get_data();

        for( size_t i = 0; i < image->get_size()/2; ++i, d+=2) {

                if ( *d < threshold ) {

                        *d = 0;

                        *(d+1) = 0;

                }

        }

                image->ap2ri();

        image->set_ri(true); // So it can be used for fourier multiplaction, for example

        image->update();

        EXITFUNC;

}


void BilateralProcessor::process_inplace(EMData * image)

{

        if (!image) {

                LOGWARN("NULL Image");

                return;

        }


        float distance_sigma = params["distance_sigma"];

        float value_sigma = params["value_sigma"];

        int max_iter = params["niter"];

        int half_width = params["half_width"];


        if (half_width < distance_sigma) {

                LOGWARN("localwidth(=%d) should be larger than distance_sigma=(%f)\n",

                                                        half_width, distance_sigma);

        }


        distance_sigma *= distance_sigma;


        float image_sigma = image->get_attr("sigma");

        if (image_sigma > value_sigma) {

                LOGWARN("image sigma(=%f) should be smaller than value_sigma=(%f)\n",

                                                        image_sigma, value_sigma);

        }

        value_sigma *= value_sigma;


        int nx = image->get_xsize();

        int ny = image->get_ysize();

        int nz = image->get_zsize();


        if(nz==1) { //for 2D image

                int width=nx, height=ny;


                int i,j,m,n;


                float tempfloat1,tempfloat2,tempfloat3;

                int       index1,index2,index;

                int       Iter;

                int       tempint1,tempint3;


                tempint1=width;

                tempint3=width+2*half_width;


                float* mask=(float*)calloc((2*half_width+1)*(2*half_width+1),sizeof(float));

                float* OrgImg=(float*)calloc((2*half_width+width)*(2*half_width+height),sizeof(