EMAN2
proc2d.cpp
Go to the documentation of this file.
00001 /*
00002  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00003  * Copyright (c) 2000-2006 Baylor College of Medicine
00004  * 
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  * 
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  * 
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  * 
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00024  * GNU General Public License for more details.
00025  * 
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00029  * 
00030  * */
00031 
00032 #include "emdata.h"
00033 #include "ctf.h"
00034 
00035 #include <cctype>
00036 
00037 using namespace EMAN;
00038 using std::map;
00039 
00040 void usage(const char *progname)
00041 {
00042     printf("\n%s inputfile outputfile [-vN] [first=<n>] [last=<n>] [inplace] [edgenorm] [norm[=<mean>,<sigma>]] [invert] [flip] [rot=<angle>] [trans=<dx>,<dy>] [center] [acfcenter] [scale=<sca>] [clip=<x,y>] [shrink=<n>] [meanshrink=<n>] [apix=<A/pix>] [lp=<filt r>] [hp=<filt r>]  [blfilt=<sigma1,sigma2,iter,localwidth>] [mask=<radius>] [imask=<radius>] [noisemask]  [mrc]  [automask=x] [pif] [hdf] [png] [em] [spider] [spider-single] [pgm[=<low>,<hi>]] [vtk] [sharphp=<pixels>] [norefs] [rotav] [average]  [split=n] [ctfsplit] [interlv=<file 2 interleave>] [sym=<Cn>] [plt[=<plt file>]] [setsfpairs] [addnoise=<level>] [randomize=<dr>,<da>,<flip>]  [selfcl[=<steps>][,<mode>]] [radon] [rfp] [mraprep] [phot] [rsub] rfilt=<filtername><:param1=value1><...>\n", progname);
00043 }
00044 
00045 int main(int argc, char *argv[])
00046 {
00047 
00048     if (argc < 3) {
00049                 usage(argv[0]);
00050                 exit(0);
00051     }
00052 
00053     Util::set_log_level(argc, argv);
00054 
00055     const int MAXMICROCTF = 1000;
00056     float *defocus_val = new float[MAXMICROCTF];
00057     float *bfactor_val = new float[MAXMICROCTF];
00058 
00059     string center = "center";
00060     string acfcenter = "acfcenter";
00061     string invert = "invert";
00062     string phot = "phot";
00063 
00064     string avgonly = "average";
00065     EMData *average = 0;
00066 
00067     string rfp = "rfp";
00068     string noisemask = "noisemask";
00069     string mraprep = "mraprep";
00070     string ctfsplit = "ctfsplit";
00071     string setsfpairs = "setsfpairs";
00072     string norefs = "norefs";
00073     string mrc8bit = "mrc8bit";
00074 
00075     string mrc = "mrc";
00076     string pif = "pif";
00077     string png = "png";
00078     string hdf = "hdf";
00079         string vtk = "vtk";
00080         
00081     string em = "em";
00082     string spidersingle = "spider-single";
00083 
00084     string spider = "spider";
00085     string rotav = "rotav";
00086     string inplace = "inplace";
00087     string rsub = "rsub";
00088 
00089     string flip = "flip";
00090     string edgenorm = "edgenorm";
00091     string radon = "radon";
00092 
00093     float blsigma1 = 0;
00094     float blsigma2 = 0;
00095     int bliter = 0;
00096     int blwidth = 0;
00097 
00098     map < string, bool > argdict;
00099     map < string, string > argfilters;
00100 
00101     float apix = 0;
00102     float lp = 0;
00103     float tlp = 0;
00104 
00105     float hp = 0;
00106     float sharphp = 0;
00107 
00108     int norm = 0;
00109     float nmean = 0;
00110     float nsig = 0;
00111 
00112     float ramp = 0;
00113     int mask = 0;
00114     int imask = 0;
00115     float automask = 0;
00116     
00117     float anoise = 0;
00118     float rot = 0;
00119     float dx = 0;
00120     float dy = 0;
00121 
00122     int rize = 0;
00123     int rizef = 0;
00124     float rizedx = 0;
00125     float rizeda = 0;
00126 
00127     float scale = 1;
00128     int clipx = -1;
00129     int clipy = -1;
00130 
00131     int shrink = 1;
00132     int csym = 0;
00133 
00134     int scl = 0;
00135     int sclmd = 0;
00136 
00137     int pgm = 0;
00138     int pgmlo = 0;
00139     int pgmhi = 0;
00140 
00141     bool rfilt = false;
00142     
00143     const char *interleaved_file = 0;
00144 
00145     int n0 = 0;
00146     int n1 = -1;
00147 
00148     argfilters[noisemask] = "MakeRadiusSquared";
00149     string filtername;
00150     Dict params_dict;
00151     
00152     for (int i = 3; i < argc; i++) {
00153                 const char *arg = argv[i];
00154 
00155                 if (Util::get_str_float(arg, "apix=", &apix)) {
00156                 }
00157                 else if (Util::get_str_float(arg, "lp=", &lp)) {
00158                 }
00159                 else if (Util::get_str_float(arg, "tlp=", &tlp)) {
00160                 }
00161                 else if (Util::get_str_float(arg, "hp=", &hp)) {
00162                 }
00163                 else if (Util::get_str_float(arg, "sharphp=", &sharphp)) {
00164                 }
00165                 else if (Util::get_str_int(arg, "first=", &n0)) {
00166                 }
00167                 else if (Util::get_str_int(arg, "last=", &n1)) {
00168                 }
00169                 else if (Util::get_str_float(arg, "ramp=", &ramp)) {
00170                 }
00171                 else if (Util::get_str_int(arg, "mask=", &mask)) {
00172                 }
00173                 else if (Util::get_str_int(arg, "imask=", &imask)) {
00174                 }
00175                 else if (Util::get_str_float(arg, "addnoise=", &anoise)) {
00176                 }
00177                 else if (Util::get_str_float(arg, "rot=", &rot)) {
00178                 }
00179                 else if (Util::get_str_float(arg, "automask=", &automask)) {
00180                 }
00181                 else if (Util::get_str_float(arg, "trans=", &dx, &dy)) {
00182                 }
00183                 else if (Util::get_str_float(arg, "scale=", &scale)) {
00184                 }
00185                 else if (Util::get_str_int(arg, "clip=", &clipx, &clipy)) {
00186                 }
00187                 else if (Util::get_str_int(arg, "shrink=", &shrink)) {
00188                 }
00189                 else if (strncmp(argv[i], "blfilt=", 7) == 0) {
00190                         sscanf(argv[i] + 7, "%f,%f,%d,%d", &blsigma1, &blsigma2, &bliter, &blwidth);
00191                 }
00192                 else if (strncmp(arg, "sym=", 4) == 0) {
00193                         if (tolower(argv[i][4]) != 'c') {
00194                                 LOGERR("Only C symmetry currently supported!\n");
00195                                 exit(1);
00196                         }
00197                         csym = atoi(&argv[i][5]);
00198                 }
00199                 else if (strncmp(arg, "randomize=", 10) == 0) {
00200                         rize = 1;
00201                         sscanf(&argv[i][10], "%f,%f,%d", &rizedx, &rizeda, &rizef);
00202                         rizeda *= M_PI / 180.0f;
00203                 }
00204                 else if (Util::get_str_float(arg, "norm", &norm, &nmean, &nsig)) {
00205                         if (norm != 2) {
00206                                 norm = 1;
00207                         }
00208                 }
00209                 else if (strncmp(argv[i], "selfcl", 6) == 0) {
00210                         if (argv[i][6] == '=') {
00211                                 sscanf(argv[i] + 7, "%d,%d", &scl, &sclmd);
00212                                 scl /= 2;
00213                         }
00214                         else {
00215                                 scl = 90;
00216                         }
00217                 }
00218                 else if (strncmp(arg, "interlv=", 8) == 0) {
00219                         interleaved_file = arg + 8;
00220                 }
00221                 else if (Util::get_str_int(arg, "pgm", &pgm, &pgmlo, &pgmhi)) {
00222                         if (pgm != 2) {
00223                                 pgm = 1;
00224                         }
00225                 }
00226                 else if (Util::sstrncmp(arg, noisemask.c_str())) {
00227                         argfilters[noisemask] = "NoiseMask";
00228                 }
00229                 else if (strncmp(arg, "rfilt=", 6) == 0) {
00230                         rfilt = true;
00231                         string s(&arg[6]);
00232                         size_t cpos = s.find(":");
00233             
00234                         if (cpos == string::npos) {
00235                                 filtername = s;
00236                         }
00237                         else {
00238                                 filtername = s.substr(0, cpos);
00239                                 size_t slen = s.length();
00240                 
00241                                 while (cpos < slen) {
00242                                         size_t eqpos = s.find("=", cpos);
00243                                         string param1 = s.substr(cpos+1, eqpos-cpos-1);
00244                     
00245                                         cpos = s.find(":", cpos+1);
00246                                         string val1;
00247                                         if (cpos == string::npos) {
00248                                                 val1 = s.substr(eqpos+1);
00249                                                 cpos = slen;
00250                                         }
00251                                         else {
00252                                                 val1 = s.substr(eqpos+1, cpos-eqpos-1);
00253                                         }
00254                     
00255                                         params_dict[param1] = atof(val1.c_str());
00256                                 }
00257                         }
00258             
00259                 }
00260                 else {
00261                         argdict[arg] = true;
00262                 }
00263     }
00264 
00265 
00266 
00267     EMData *d = new EMData();
00268 
00269     int nimg = EMUtil::get_image_count(argv[1]);
00270     if (nimg <= n1 || n1 < 0) {
00271                 n1 = nimg - 1;
00272     }
00273     
00274     char outfile[1024];
00275     sprintf(outfile, "%s", argv[2]);
00276     char spiderformat[256];
00277 
00278     EMData *ld = new EMData();
00279     vector < float >sfcurve1;
00280     try {
00281                 for (int i = n0; i <= n1; i++) {
00282                         d->read_image(argv[1], i);
00283                         int nx = d->get_xsize();
00284         
00285                         if ((argdict[ctfsplit] != 0) && (i == n0 || !EMUtil::is_same_ctf(d, ld))) {
00286                                 Ctf *ctf = d->get_ctf();
00287 
00288                                 int j = 0;
00289                                 for (j = 1; j < argdict[ctfsplit]; j++) {
00290                                         if (defocus_val[j] == ctf->get_defocus() && bfactor_val[j] == ctf->get_bfactor()) {
00291                                                 break;
00292                                         }
00293                                 }
00294                                 if (argdict[ctfsplit] <= j) {
00295                                         argdict[ctfsplit] = j + 1;
00296                                         printf("New CTF at %d\n", i);
00297                                 }
00298                                 defocus_val[j] = ctf->get_defocus();
00299                                 bfactor_val[j] = ctf->get_bfactor();
00300                                 outfile[strlen(outfile) - 4] = 0;
00301                                 sprintf(outfile + strlen(outfile), ".%02d.img", j);
00302 
00303                                 if( ld )
00304                                 {
00305                                         delete ld;
00306                                         ld = 0;
00307                                 }
00308                                 ld = d->copy();
00309                         }
00310 
00311                         float sigma = d->get_attr("sigma");
00312                         if (!Util::goodf(&sigma)) {
00313                                 LOGWARN("Warning! bad Sigma for image %d in file '%s'", i, argv[1]);
00314                                 continue;
00315                         }
00316 #if 0
00317                         if (argdict[norefs] && d->get_nimg() < 0) {
00318                                 continue;
00319                         }
00320 #endif
00321                         if (argdict[edgenorm]) {
00322                                 d->process_inplace("normalize.circlemean");
00323                         }
00324 
00325                         if (norm == 1) {
00326                                 d->process_inplace("normalize");
00327                         }
00328                         else if (norm == 2) {
00329                                 (*d) *= nsig / sigma;
00330                                 (*d) += nmean - (float)d->get_attr("mean");
00331                         }
00332 
00333                         if (argdict[flip]) {
00334                                 d->process_inplace("xform.flip", Dict("axis", "y"));
00335                         }
00336 
00337                         if (argdict[invert]) {
00338                                 (*d) *= -1;
00339                         }
00340 
00341                         if (ramp) {
00342                                 d->process_inplace("eman1.filter.ramp", Dict("intercept", 1, "slope", ramp));       
00343                         }
00344 
00345                         int y = d->get_ysize();
00346 
00347                         if (argdict[setsfpairs]) {
00348                                 if (mask) {
00349                                         d->process_inplace(argfilters[noisemask], Dict("outer_radius", mask));
00350                                 }
00351                                 EMData *dataf = d->do_fft();
00352 
00353                                 float x0 = 0;
00354                                 float step = 0.5;
00355             
00356                                 if (i % 2 == 0) {
00357                                         sfcurve1 = dataf->calc_radial_dist(nx, x0, step);
00358                                 }
00359                                 else {
00360                                         vector < float >sfcurve2 = dataf->calc_radial_dist(nx, x0, step);
00361 
00362                                         for (int j = 0; j < nx; j++) {
00363                                                 if (sfcurve1[j] > 0 && sfcurve2[j] > 0) {
00364                                                         sfcurve2[j] = sqrt(sfcurve1[j] / sfcurve2[j]);
00365                                                 }
00366                                                 else {
00367                                                         sfcurve2[j] = 0;
00368                                                 }
00369                                         }
00370 
00371                                         dataf->apply_radial_func(x0, step, sfcurve2);
00372                                         if( d )
00373                                         {
00374                                                 delete d;
00375                                                 d = 0;
00376                                         }
00377                                         d = dataf->do_ift();
00378                                 }
00379                                 
00380                                 if( dataf )
00381                                 {
00382                                         delete dataf;
00383                                         dataf = 0;
00384                                 }
00385                         }
00386 
00387                         float Xlp = lp;
00388                         float Xtlp = tlp;
00389                         float Xhp = hp;
00390                         float Xsharphp = sharphp;
00391 
00392                         if (apix > 0 && lp) {
00393                                 Xlp = y * apix / lp;
00394                         }
00395 
00396                         if (apix > 0 && tlp) {
00397                                 Xtlp = y * apix / tlp;
00398                         }
00399 
00400                         if (apix > 0 && hp) {
00401                                 Xhp = y * apix / hp;
00402                         }
00403 
00404                         if (apix > 0 && sharphp) {
00405                                 Xsharphp = y * apix / sharphp;
00406                         }
00407 
00408                         if (Xlp || Xhp) {
00409                                 d->process_inplace("filter.lowpass.gauss", Dict("cutoff_abs", Xhp == 0 ? -10.0 : Xhp));
00410                                 d->process_inplace("filter.highpass.tanh", Dict("cutoff_abs", Xlp == 0 ? 100000.0 : Xlp));
00411                         }
00412 
00413                         if (Xtlp) {
00414                                 d->process_inplace("filter.lowpass.tanh", Dict("cutoff_abs", -10.0));
00415                                 d->process_inplace("filter.highpass.tanh", Dict("cutoff_abs", Xtlp));
00416                         }
00417 
00418                         if (Xsharphp) {
00419                                 d->process_inplace("filter.lowpass.tophat", Dict("cutoff_abs", Xsharphp));
00420                                 d->process_inplace("filter.highpass.tophat", Dict("cutoff_abs", 100000.0));
00421                         }
00422 
00423                         if (mask) {
00424                                 d->process_inplace(argfilters[noisemask], Dict("outer_radius", mask));
00425                         }
00426 
00427                         if (imask > 0) {
00428                                 d->process_inplace("mask.sharp", Dict("inner_radius", imask, "value", 0));
00429                         }
00430 
00431                         if (automask) {
00432                                 d->process_inplace("mask.auto2d", Dict("threshold", automask));
00433                         }
00434 
00435                         // uses correlation with 180 deg rot for centering
00436                         if (argdict[acfcenter]) {
00437                                 Dict params;
00438                                 params["intonly"] = 1;
00439                                 params["maxshift"] = d->get_xsize() / 4;
00440                                 d->align("translate", 0, params);
00441                                 //d->rotate_translate();
00442                         }
00443 
00444                         if (argdict[center]) {
00445                                 d->process_inplace("xform.centerofmass", Dict("int_shift_only", 1));
00446                         }
00447 
00448                         if (argdict[phot]) {
00449                                 d->process_inplace("xform.phaseorigin");
00450                         }
00451 
00452                         if (anoise) {
00453                                 d->process_inplace("math.addnoise");
00454                         }
00455 
00456                         if (argdict[rfp]) {
00457                                 EMData *e = d->make_rotational_footprint();
00458                                 e->append_image("rfp.hed");
00459                         }
00460 
00461 
00462                         if (rot || dx || dy || rize) {
00463                                 if (!rize) {
00464                                         d->rotate_translate(rot * M_PI / 180.0f, 0, 0, dx, dy, 0);
00465                                 }
00466                                 else {
00467                                         if (rizef && rand() % 2) {
00468                                                 d->process_inplace("xform.flip", Dict("axis", "y"));
00469                                         }
00470 
00471                                         if (rizeda > 0) {
00472                                                 d->rotate(Util::get_frand(-rizeda / 2.0f, rizeda / 2.0f), 0, 0);
00473                                         }
00474                 
00475                                         if (rizedx > 0) {
00476                                                 d->translate(Util::get_gauss_rand(0, rizedx),
00477                                                                          Util::get_gauss_rand(0, rizedx), 0.0f);
00478                                         }
00479                                 }
00480                         }
00481 
00482                         Dict rr = d->get_transform().get_rotation("eman");
00483                         //int nimg = d->get_nimg();
00484 
00485                         if (scale < 1.0) {
00486                                 Transform t;
00487                                 t.set_scale(scale);
00488                                 d->transform(t);
00489                         }
00490 
00491                         if (clipx > 0) {
00492                                 EMData *e = d->get_clip(Region((d->get_xsize() - clipx) / 2,
00493                                                                                            (d->get_ysize() - clipy) / 2,
00494                                                                                            clipx, clipy));
00495                                 if( d )
00496                                 {
00497                                         delete d;
00498                                         d = 0;
00499                                 }
00500                                 d = e;
00501                         }
00502 
00503                         if (scale > 1.0) {
00504                                 Transform t;
00505                                 t.set_scale(scale);
00506                                 d->transform(t);
00507                         }
00508 
00509                         d->set_rotation((float)rr["alt"], (float)rr["az"], (float)rr["phi"]);
00510                         //d->setNImg(nimg);
00511 
00512                         if (fabs((float)shrink) > 1) {
00513                                 d->median_shrink((int) fabs((float)shrink));
00514                         }
00515 
00516                         if (argdict[rotav]) {
00517                                 d->process_inplace("math.rotationalaverage");
00518                         }
00519 
00520                         if (csym > 1) {
00521                                 EMData *f = d->copy();
00522 
00523                                 EMData *e = f->copy();
00524                                 for (int j = 1; j < csym; j++) {
00525                                         e->rotate(j * M_PI * 2.0f / csym, 0, 0);
00526                                         (*d) += (*e);
00527                                 }
00528 
00529                                 *d *= (1.0f / csym);
00530 
00531                                 if( e )
00532                                 {
00533                                         delete e;
00534                                         e = 0;
00535                                 }
00536                                 if( f )
00537                                 {
00538                                         delete f;
00539                                         f = 0;
00540                                 }
00541                         }
00542 
00543                         if (argdict[rsub]) {
00544                                 d->process_inplace("math.rotationalsubtract");
00545                         }
00546 
00547                         if (scl) {
00548                                 EMData *sc = new EMData();
00549                                 if (sclmd == 0) {
00550                                         sc->common_lines_real(d, d, scl, true);
00551                                 }
00552                                 else {
00553                                         EMData *e = d->copy();
00554                                         e->process_inplace("xform.phaseorigin");
00555 
00556                                         if (sclmd == 1) {
00557                                                 sc->common_lines(e, e, sclmd, scl, true);
00558                                                 sc->process_inplace("math.linear", Dict("shift", -90.0, "scale", -1.0));
00559                                         }
00560                                         else if (sclmd == 2) {
00561                                                 sc->common_lines(e, e, sclmd, scl, true);
00562                                         }
00563                                         else {
00564                                                 LOGERR("Invalid common-lines mode");
00565                                                 exit(1);
00566                                         }
00567                                         if( e )
00568                                         {
00569                                                 delete e;
00570                                                 e = 0;
00571                                         }
00572                                 }
00573                                 if( d )
00574                                 {
00575                                         delete d;
00576                                         d = 0;
00577                                 }
00578                                         d = sc;
00579                                 
00580                         }
00581 
00582                         if (argdict[radon]) {
00583                                 EMData *r = d->do_radon();
00584                                 if( d )
00585                                 {
00586                                         delete d;
00587                                         d = 0;
00588                                 }
00589                                 d = r;
00590                         }
00591 
00592                         if (rfilt) {
00593                                 d->process_inplace(filtername, params_dict);
00594                         }
00595 
00596                         if (bliter > 0 && blwidth > 0) {
00597                                 Dict p;
00598                                 p["distance_sigma"] = blsigma1;
00599                                 p["value_sigma"] = blsigma2;
00600                                 p["niter"] = bliter;
00601                                 p["half_width"] = blwidth;
00602 
00603                                 d->process_inplace("bilateral", p);
00604                         }
00605 
00606 #if 0
00607                         if (filefilt) {
00608                                 EMData *d2 = d->do_fft();
00609                                 if( d )
00610                                 {
00611                                         delete d;
00612                                         d = 0;
00613                                 }
00614 
00615                                 d2->apply_radial_func(nxyd, xd[0], xd[1] - xd[0], yd, true);
00616                                 d = d2->do_ift();
00617                         }
00618 #endif
00619 
00620                         if (argdict[avgonly]) {
00621                                 if (!average) {
00622                                         average = d->copy();
00623                                 }
00624                                 else {
00625                                         (*average) += (*d);
00626                                 }
00627                                 continue;
00628                         }
00629 #if 0
00630                         if (fftavgfile) {
00631                                 if (!fftavg) {
00632                                         fftavg = new EMData;
00633                                         fftavg->setSize(d->get_xsize() + 2, d->get_ysize());
00634                                         fftavg->setComplex(1);
00635                                         fftavg->zero();
00636                                 }
00637                                 d->applyMask(-1, 6);
00638                                 d->normalize();
00639                                 EMData *df = d->do_fft();
00640                                 df->multConst(df->get_ysize()); // just for scaling of the intensity level
00641                                 fftavg->addIncoherent(df);
00642 
00643                                 if( df )
00644                                 {
00645                                         delete df;
00646                                         df = 0;
00647                                 }
00648                                 continue;               // no writing yet
00649                         }
00650 #endif
00651 
00652 #if 0
00653                         if (strlen(sfout)) {
00654             
00655                                 int y = d->get_ysize();
00656                                 float *curve = (float *) malloc(y * 4 * sizeof(float));
00657                                 EMData *dataf = d->do_fft();
00658 
00659                                 dataf->calc_radial_dist(y, 0, .5, curve);
00660                                 // single images go to outfile
00661                                 if (n1 == 0) {
00662                                         save_data(0, 1.0 / (apix * 2.0 * y), curve, y, sfout);
00663                                 }
00664                                 // multi-images get numbered
00665                                 else {
00666                                         sprintf(outfile2, "%s.%03d", sfout, i + 100);
00667                                         save_data(0, 1.0 / (apix * 2.0 * y), curve, y, outfile2);
00668                                 }
00669                                 if (sfout_n) {
00670                                         for (int j = 0; j < sfout_n; j++) {
00671                                                 dataf->calc_radial_dist(y, 0, .5, curve, j * 2 * M_PI / sfout_n,
00672                                                                                                 2 * M_PI / sfout_n);
00673                                                 // single images go to outfile
00674                                                 if (n1 == 0) {
00675                                                         sprintf(outfile2, "%s-%d-%d.pwr", file_basename(sfout), j, sfout_n);
00676                                                         save_data(0, 1.0 / (apix * 2.0 * y), curve, y, outfile2);
00677                                                 }
00678                                                 // multi-images get numbered
00679                                                 else {
00680                                                         sprintf(outfile2, "%s-%d-%d.pwr.%03d", file_basename(sfout), j, sfout_n,
00681                                                                         i + 100);
00682                                                         save_data(0, 1.0 / (apix * 2.0 * y), curve, y, outfile2);
00683                                                 }
00684                                         }
00685                                 }
00686                         }
00687 #endif
00688 
00689                         if (argdict[mrc]) {
00690                                 if (n1 != 0) {
00691                                         sprintf(outfile, "%03d.%s", i + 100, outfile);
00692                                 }
00693                                 d->write_image(outfile, 0, EMUtil::IMAGE_MRC);
00694                         }
00695                         else if (argdict[spidersingle]) {
00696                                 if (n1 != 0) {
00697                                         sprintf(spiderformat, "%s%%0%dd.spi", outfile, (int) (log10((float) i)) + 1);
00698                                         sprintf(outfile, spiderformat, i);
00699                                 }
00700 
00701                                 d->write_image(outfile, 0, EMUtil::IMAGE_SINGLE_SPIDER);
00702                         }
00703                         else if (argdict[hdf]) {
00704                                 d->append_image(outfile,  EMUtil::IMAGE_HDF);
00705                         }
00706                         else if (argdict[em]) {
00707                                 d->write_image(outfile, 0, EMUtil::IMAGE_EM);
00708                         }
00709                         else if (argdict[pif]) {
00710                                 if (n1 != 0) {
00711                                         sprintf(outfile, "%03d.%s", i + 100, outfile);
00712                                 }
00713 
00714                                 d->write_image(outfile, 0, EMUtil::IMAGE_PIF);
00715                         }
00716                         else if (argdict[png]) {
00717                                 if (n1 != 0) {
00718                                         sprintf(outfile, "%03d.%s", i + 100, outfile);
00719                                 }
00720 
00721                                 d->write_image(outfile, 0, EMUtil::IMAGE_PNG);
00722                         }
00723                         else if (pgm) {
00724                                 if (pgmlo >= pgmhi) {
00725                                         pgmlo = (int) d->get_attr("minimum");
00726                                         pgmhi = (int) d->get_attr("maximum");
00727                                 }
00728 
00729                                 d->set_attr("min_gray", pgmlo);
00730                                 d->set_attr("max_gray", pgmhi);
00731 
00732                                 if (n1 != 0) {
00733                                         sprintf(outfile, "%03d.%s", i + 100, outfile);
00734                                 }
00735                                 d->write_image(outfile, 0, EMUtil::IMAGE_PGM);
00736                         }
00737                         else if (argdict[spider]) {
00738                                 d->write_image(outfile, 0, EMUtil::IMAGE_SPIDER);
00739                         }
00740                         else if (argdict[vtk]) {
00741                                 d->write_image(outfile, 0, EMUtil::IMAGE_VTK);
00742                         }
00743                         else {
00744                                 if (argdict[inplace]) {
00745                                         d->write_image(outfile, i);
00746                                 }
00747                                 else if (argdict[mraprep]) {
00748                                         char nf[1024];
00749                                         sprintf(nf, "%s%04d.lst", outfile, i);
00750                                         d->write_image(nf, 0, EMUtil::IMAGE_LST);
00751                                 }
00752                                 else {
00753                                         d->append_image(outfile);
00754                                 }
00755                         }
00756 
00757                         if (interleaved_file) {
00758                                 d->read_image(interleaved_file, i);
00759                                 d->write_image(outfile, -1, EMUtil::IMAGE_IMAGIC);
00760                         }
00761 
00762                 }
00763         }
00764         catch(E2Exception &e ) {
00765                 printf("%s\n", e.what());
00766         }
00767         
00768     if (average) {
00769                 //average->setNImg(n1-n0+1);
00770                 average->process_inplace("normalize");
00771                 average->write_image(outfile, -1);
00772     }
00773 #if 0
00774     if (fftavg) {
00775                 fftavg->multConst(1. / sqrt((float) n1 - n0 + 1));
00776                 fftavg->writeMRC(fftavgfile);
00777     }
00778 #endif
00779     int n_outimg = EMUtil::get_image_count(argv[2]);
00780     printf("%d images\n", n_outimg);
00781 
00782     if (d) {
00783                 delete d;
00784                 d = 0;
00785     }
00786 
00787     if (ld) {
00788                 delete ld;
00789                 ld = 0;
00790     }
00791 
00792         if( defocus_val )
00793         {
00794         delete [] defocus_val;
00795         defocus_val = 0;
00796         }
00797 
00798         if( bfactor_val )
00799         {
00800         delete [] bfactor_val;
00801         bfactor_val = 0;
00802         }
00803     
00804     return 0;
00805 }