Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

projector.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "projector.h"
00037 #include "emdata.h"
00038 #include "interp.h"
00039 #include "emutil.h"
00040 #include "plugins/projector_template.h"
00041 
00042 #ifdef WIN32
00043         #define M_PI 3.14159265358979323846f
00044 #endif  //WIN32
00045 
00046 #ifdef EMAN2_USING_CUDA
00047 #include "cuda/cuda_util.h"
00048 #include "cuda/cuda_projector.h"
00049 #endif
00050 
00051 using namespace std;
00052 using namespace EMAN;
00053 
00054 const string GaussFFTProjector::NAME = "gauss_fft";
00055 const string FourierGriddingProjector::NAME = "fourier_gridding";
00056 const string PawelProjector::NAME = "pawel";
00057 const string StandardProjector::NAME = "standard";
00058 const string MaxValProjector::NAME = "maxval";
00059 const string ChaoProjector::NAME = "chao";
00060 
00061 template <> Factory < Projector >::Factory()
00062 {
00063         force_add<GaussFFTProjector>();
00064         force_add<PawelProjector>();
00065         force_add<StandardProjector>();
00066         force_add<MaxValProjector>();
00067         force_add<FourierGriddingProjector>();
00068         force_add<ChaoProjector>();
00069         
00070 //      force_add<XYZProjector>();
00071 }
00072 
00073 EMData *GaussFFTProjector::project3d(EMData * image) const
00074 {
00075         Transform* t3d = params["transform"];
00076         if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00077         if ( image->get_ndim() != 3 ) throw ImageDimensionException("Error, the projection volume must be 3D");
00078 
00079         EMData *f = image;
00080         if (!image->is_complex()) {
00081                 image->process_inplace("xform.phaseorigin.tocorner");
00082                 f = image->do_fft();
00083                 f->process_inplace("xform.fourierorigin.tocenter");
00084                 image->process_inplace("xform.phaseorigin.tocenter");
00085         }
00086 
00087         int f_nx = f->get_xsize();
00088         int f_ny = f->get_ysize();
00089         int f_nz = f->get_zsize();
00090 
00091         if (!f->is_complex() || f_nz != f_ny || f_nx != f_ny + 2) {
00092                 LOGERR("Cannot project this image");
00093                 return 0;
00094         }
00095 
00096         f->ap2ri();
00097 
00098         EMData *tmp = new EMData();
00099         tmp->set_size(f_nx, f_ny, 1);
00100         tmp->set_complex(true);
00101         tmp->set_ri(true);
00102 
00103         float *data = tmp->get_data();
00104 
00105         Transform r = t3d->get_rotation_transform();
00106         r.invert();
00107         float scale = t3d->get_scale();
00108 
00109         int mode = params["mode"];
00110         float gauss_width = 1;
00111         if ( mode == 0 ) mode = 2;
00112         if (mode == 2 ) {
00113                 gauss_width = EMConsts::I2G;
00114         }
00115         else if (mode == 3) {
00116                 gauss_width = EMConsts::I3G;
00117         }
00118         else if (mode == 4) {
00119                 gauss_width = EMConsts::I4G;
00120         }
00121         else if (mode == 6 || mode == 7) {
00122                 gauss_width = EMConsts::I5G;
00123         }
00124 
00125         for (int y = 0; y < f_ny; y++) {
00126                 for (int x = 0; x < f_nx / 2; x++) {
00127                         int ii = x * 2 + y * f_nx;
00128 #ifdef  _WIN32
00129                         if (_hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
00130 #else
00131                         if (hypot(x, y - f_ny / 2) >= f_ny / 2 - 1) {
00132 #endif  //_WIN32
00133                                 data[ii] = 0;
00134                                 data[ii + 1] = 0;
00135                         }
00136                         else {
00137                                 Vec3f coord(x,(y - f_ny / 2),0);
00138                                 Vec3f soln = r*coord;
00139                                 float xx = soln[0];
00140                                 float yy = soln[1];
00141                                 float zz = soln[2];
00142 
00143                                 int cc = 1;
00144                                 if (xx < 0) {
00145                                         xx = -xx;
00146                                         yy = -yy;
00147                                         zz = -zz;
00148                                         cc = -1;
00149                                 }
00150 
00151                                 if (scale != 1.0) {
00152                                         xx *= scale;
00153                                         yy *= scale;
00154                                         zz *= scale;
00155                                 }
00156 
00157                                 yy += f_ny / 2;
00158                                 zz += f_nz / 2;
00159 
00160                                 if (yy < 0 || xx < 0 || zz < 0) {
00161                                         data[ii] = 0;
00162                                         data[ii+1] = 0;
00163                                         continue;
00164                                 }
00165 
00166                                 if (interp_ft_3d(mode, f, xx, yy, zz, data + ii, gauss_width)) {
00167                                         data[ii + 1] *= cc;
00168                                 } else {
00169                                         data[ii] = 0;
00170                                         data[ii+1] = 0;
00171                                 }
00172                         }
00173                 }
00174         }
00175 
00176         f->update();
00177         tmp->update();
00178 
00179         tmp->process_inplace("xform.fourierorigin.tocorner");
00180         EMData *ret = tmp->do_ift();
00181         ret->process_inplace("xform.phaseorigin.tocenter");
00182 
00183         ret->translate(t3d->get_trans());
00184 
00185         if (t3d->get_mirror() ) ret->process_inplace("xform.flip",Dict("axis","x"));
00186 
00187         Dict filter_d;
00188         filter_d["gauss_width"] = gauss_width;
00189         filter_d["ring_width"] = ret->get_xsize() / 2;
00190         ret->process_inplace("math.gausskernelfix", filter_d);
00191 
00192         if( tmp )
00193         {
00194                 delete tmp;
00195                 tmp = 0;
00196         }
00197 
00198         ret->set_attr("xform.projection",t3d);
00199         ret->update();
00200 
00201         if(t3d) {delete t3d; t3d=0;}
00202 
00203         return ret;
00204 }
00205 
00206 
00207 
00208 bool GaussFFTProjector::interp_ft_3d(int mode, EMData * image, float x, float y,
00209                                                                          float z, float *data, float gw) const
00210 {
00211         float *rdata = image->get_data();
00212         int nx = image->get_xsize();
00213         int ny = image->get_ysize();
00214         int nz = image->get_zsize();
00215 
00216         if ( mode == 0 ) mode = 2;
00217 
00218         if (mode == 1) {
00219                 int x0 = 2 * (int) floor(x + 0.5f);
00220                 int y0 = (int) floor(y + 0.5f);
00221                 int z0 = (int) floor(z + 0.5f);
00222 
00223                 data[0] = rdata[x0 + y0 * nx + z0 * nx * ny];
00224                 data[1] = rdata[x0 + y0 * nx + z0 * nx * ny + 1];
00225                 return true;
00226         }
00227         else if (mode == 2) {
00228                 int x0 = (int) floor(x);
00229                 int y0 = (int) floor(y);
00230                 int z0 = (int) floor(z);
00231 
00232                 float dx = x - x0;
00233                 float dy = y - y0;
00234                 float dz = z - z0;
00235 
00236                 if (x0 <= nx/2- 2 && y0 <= ny - 1 && z0 <= nz - 1) {
00237                         int i = (int) (x0 * 2 + y0 * nx + z0 * nx * ny);
00238 
00239                         float n = Util::agauss(1, dx, dy, dz, gw) +
00240                                 Util::agauss(1, 1 - dx, dy, dz, gw) +
00241                                 Util::agauss(1, dx, 1 - dy, dz, gw) +
00242                                 Util::agauss(1, 1 - dx, 1 - dy, dz, gw) +
00243                                 Util::agauss(1, dx, dy, 1 - dz, gw) +
00244                                 Util::agauss(1, 1 - dx, dy, 1 - dz, gw) +
00245                                 Util::agauss(1, dx, 1 - dy, 1 - dz, gw) + Util::agauss(1, 1 - dx, 1 - dy, 1 - dz, gw);
00246 
00247                         data[0] = Util::agauss(rdata[i], dx, dy, dz, gw) +
00248                                 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
00249                                 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
00250                                 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
00251                                 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
00252                                 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
00253                                 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
00254                                 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
00255 
00256                         i++;
00257 
00258                         data[1] = Util::agauss(rdata[i], dx, dy, dz, gw) +
00259                                 Util::agauss(rdata[i + 2], 1 - dx, dy, dz, gw) +
00260                                 Util::agauss(rdata[i + nx], dx, 1 - dy, dz, gw) +
00261                                 Util::agauss(rdata[i + nx + 2], 1 - dx, 1 - dy, dz, gw) +
00262                                 Util::agauss(rdata[i + nx * ny], dx, dy, 1 - dz, gw) +
00263                                 Util::agauss(rdata[i + 2 + nx * ny], 1 - dx, dy, 1 - dz, gw) +
00264                                 Util::agauss(rdata[i + nx + nx * ny], dx, 1 - dy, 1 - dz, gw) +
00265                                 Util::agauss(rdata[i + 2 + nx + nx * ny], 1 - dx, 1 - dy, 1 - dz, gw) / n;
00266                         return true;
00267                 }
00268                 return false;
00269         }
00270         else if (mode == 3) {
00271                 int x0 = 2 * (int) floor(x + .5);
00272                 int y0 = (int) floor(y + .5);
00273                 int z0 = (int) floor(z + .5);
00274 
00275                 float *supp = image->setup4slice();
00276 
00277                 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00278                         float n = 0;
00279 
00280                         if (x0 == 0) {
00281                                 x0 += 4;
00282                                 size_t idx;
00283                                 float r, g;
00284                                 for (int k = z0 - 1; k <= z0 + 1; k++) {
00285                                         for (int j = y0 - 1; j <= y0 + 1; j++) {
00286                                                 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
00287                                                         r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00288                                                         g = exp(-r / gw);
00289                                                         n += g;
00290                                                         idx = i + j * 12 + (size_t)k * 12 * ny;
00291                                                         data[0] += supp[idx] * g;
00292                                                         data[1] += supp[idx + 1] * g;
00293                                                 }
00294                                         }
00295                                 }
00296                         }
00297                         else {
00298                                 size_t idx;
00299                                 float r, g;
00300                                 for (int k = z0 - 1; k <= z0 + 1; k++) {
00301                                         for (int j = y0 - 1; j <= y0 + 1; j++) {
00302                                                 for (int i = x0 - 2; i <= x0 + 2; i += 2) {
00303                                                         r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00304                                                         g = exp(-r / gw);
00305                                                         n += g;
00306                                                         idx = i + j * nx + (size_t)k * nx * ny;
00307                                                         data[0] += rdata[idx] * g;
00308                                                         data[1] += rdata[idx + 1] * g;
00309                                                 }
00310                                         }
00311                                 }
00312                         }
00313 
00314                         data[0] /= n;
00315                         data[1] /= n;
00316                         return true;
00317                 }
00318                 return false;
00319         }
00320         else if (mode == 4) {
00321                 int x0 = 2 * (int) floor(x);
00322                 int y0 = (int) floor(y);
00323                 int z0 = (int) floor(z);
00324 
00325                 float *supp = image->setup4slice();
00326 
00327                 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00328                         float n = 0;
00329 
00330                         if (x0 == 0) {
00331                                 x0 += 4;
00332                                 size_t idx;
00333                                 float r, g;
00334                                 for (int k = z0 - 1; k <= z0 + 2; k++) {
00335                                         for (int j = y0 - 1; j <= y0 + 2; j++) {
00336                                                 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
00337                                                         r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00338                                                         g = exp(-r / gw);
00339                                                         n += g;
00340                                                         idx = i + j * 12 + (size_t)k * 12 * ny;
00341                                                         data[0] += supp[idx] * g;
00342                                                         data[1] += supp[idx + 1] * g;
00343                                                 }
00344                                         }
00345                                 }
00346                         }
00347                         else {
00348                                 float r, g;
00349                                 size_t idx;
00350                                 for (int k = z0 - 1; k <= z0 + 2; k++) {
00351                                         for (int j = y0 - 1; j <= y0 + 2; j++) {
00352                                                 for (int i = x0 - 2; i <= x0 + 4; i += 2) {
00353                                                         r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00354                                                         g = exp(-r / gw);
00355                                                         n += g;
00356                                                         idx = i + j * nx + (size_t)k * nx * ny;
00357                                                         data[0] += rdata[idx] * g;
00358                                                         data[1] += rdata[idx + 1] * g;
00359                                                 }
00360                                         }
00361                                 }
00362                         }
00363                         data[0] /= n;
00364                         data[1] /= n;
00365                         return true;
00366                 }
00367                 return false;
00368         }
00369         else if (mode == 5) {
00370                 int x0 = 2 * (int) floor(x + .5);
00371                 int y0 = (int) floor(y + .5);
00372                 int z0 = (int) floor(z + .5);
00373 
00374                 float *supp = image->setup4slice();
00375                 float *gimx = Interp::get_gimx();
00376 
00377                 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00378                         int mx0 = -(int) floor((x - x0 / 2) * 39.0 + .5) - 78;
00379                         int my0 = -(int) floor((y - y0) * 39.0 + .5) - 78;
00380                         int mz0 = -(int) floor((z - z0) * 39.0 + .5) - 78;
00381 
00382                         float n = 0;
00383                         int mmz = mz0;
00384                         int mmy = my0;
00385                         int mmx = mx0;
00386 
00387                         if (x0 < 4) {
00388                                 x0 += 4;
00389 
00390                                 size_t idx;
00391                                 float g;
00392                                 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
00393                                         for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
00394                                                 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
00395                                                         g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
00396                                                         n += g;
00397                                                         idx = i + j * 12 + (size_t)k * 12 * ny;
00398                                                         data[0] += supp[idx] * g;
00399                                                         data[1] += supp[idx + 1] * g;
00400                                                 }
00401                                         }
00402                                 }
00403                         }
00404                         else {
00405                                 size_t ii;
00406                                 float g;
00407                                 for (int k = z0 - 2; k <= z0 + 2; k++, mmz += 39) {
00408                                         for (int j = y0 - 2; j <= y0 + 2; j++, mmy += 39) {
00409                                                 for (int i = x0 - 4; i <= x0 + 4; i += 2, mmx += 39) {
00410                                                         ii = i + j * nx + (size_t)k * nx * ny;
00411                                                         g = gimx[abs(mmx) + abs(mmy) * 100 + abs(mmz) * 10000];
00412                                                         n += g;
00413                                                         data[0] += rdata[ii] * g;
00414                                                         data[1] += rdata[ii + 1] * g;
00415                                                 }
00416                                         }
00417                                 }
00418                         }
00419 
00420                         data[0] /= n;
00421                         data[1] /= n;
00422                         return true;
00423                 }
00424                 return false;
00425         }
00426         else if (mode == 6) {
00427 
00428                 int x0 = 2 * (int) floor(x + .5);
00429                 int y0 = (int) floor(y + .5);
00430                 int z0 = (int) floor(z + .5);
00431 
00432                 float *supp = image->setup4slice();
00433 
00434                 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00435                         float n = 0;
00436 
00437                         if (x0 < 4) {
00438                                 x0 += 4;
00439                                 float r, g;
00440                                 size_t idx;
00441                                 for (int k = z0 - 2; k <= z0 + 2; k++) {
00442                                         for (int j = y0 - 2; j <= y0 + 2; j++) {
00443                                                 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00444                                                         r = Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z);
00445                                                         g = exp(-r / gw);
00446                                                         n += g;
00447                                                         idx = i + j * 12 + (size_t)k * 12 * ny;
00448                                                         data[0] += supp[idx] * g;
00449                                                         data[1] += supp[idx + 1] * g;
00450                                                 }
00451                                         }
00452                                 }
00453                         }
00454                         else {
00455                                 float r, g;
00456                                 size_t idx;
00457                                 for (int k = z0 - 2; k <= z0 + 2; k++) {
00458                                         for (int j = y0 - 2; j <= y0 + 2; j++) {
00459                                                 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00460                                                         r = Util::hypot3sq(i / 2.0f - x, j - y, k - z);
00461                                                         g = exp(-r / gw);
00462                                                         n += g;
00463                                                         idx = i + j * nx + (size_t)k * nx * ny;
00464                                                         data[0] += rdata[idx] * g;
00465                                                         data[1] += rdata[idx + 1] * g;
00466                                                 }
00467 
00468                                         }
00469                                 }
00470                         }
00471 
00472                         data[0] /= n;
00473                         data[1] /= n;
00474 
00475                         return true;
00476                 }
00477                 return false;
00478         }
00479         else if (mode == 7) {
00480                 int x0 = 2 * (int) floor(x + .5);
00481                 int y0 = (int) floor(y + .5);
00482                 int z0 = (int) floor(z + .5);
00483 
00484                 float *supp = image->setup4slice();
00485 
00486                 if (x0 < nx - 4 && y0 <= ny - 3 && z0 <= nz - 3 && y0 >= 2 && z0 >= 2) {
00487                         float n = 0;
00488                         if (x0 < 4) {
00489                                 x0 += 4;
00490                                 float r, g;
00491                                 size_t idx;
00492                                 for (int k = z0 - 2; k <= z0 + 2; k++) {
00493                                         for (int j = y0 - 2; j <= y0 + 2; j++) {
00494                                                 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00495                                                         r = sqrt(Util::hypot3sq(i / 2.0f - x - 2.0f, j - y, k - z));
00496                                                         g = Interp::hyperg(r);
00497                                                         n += g;
00498                                                         idx = i + j * 12 + (size_t)k * 12 * ny;
00499                                                         data[0] += supp[idx] * g;
00500                                                         data[1] += supp[idx + 1] * g;
00501                                                 }
00502                                         }
00503                                 }
00504                         }
00505                         else {
00506                                 float r, g;
00507                                 size_t idx;
00508                                 for (int k = z0 - 2; k <= z0 + 2; k++) {
00509                                         for (int j = y0 - 2; j <= y0 + 2; j++) {
00510                                                 for (int i = x0 - 4; i <= x0 + 4; i += 2) {
00511                                                         r = sqrt(Util::hypot3sq(i / 2.0f - x, j - y, k - z));
00512                                                         g = Interp::hyperg(r);
00513                                                         n += g;
00514                                                         idx = i + j * nx + (size_t)k * nx * ny;
00515                                                         data[0] += rdata[idx] * g;
00516                                                         data[1] += rdata[idx + 1] * g;
00517                                                 }
00518 
00519                                         }
00520                                 }
00521                         }
00522                         data[0] /= n;
00523                         data[1] /= n;
00524                         return true;
00525                 }
00526                 return false;
00527 
00528         }
00529 //      throw InvalidParameterException("Error, unkown mode");
00530         return false;
00531 }
00532 
00533 void PawelProjector::prepcubes(int nx, int ny, int nz, int ri, Vec3i origin,
00534                                        int& nn, IPCube* ipcube) const {
00535         const float r = float(ri)*float(ri);
00536         const int ldpx = origin[0];
00537         const int ldpy = origin[1];
00538         const int ldpz = origin[2];
00539         //cout<<"  ldpx  "<<ldpx<<"  i2  "<<ldpy<<"  i1 "<<ldpz<<endl;
00540         float t;
00541         nn = -1;
00542         for (int i1 = 0; i1 < nz; i1++) {
00543                 t = float(i1 - ldpz);
00544                 const float xx = t*t;
00545                 for (int i2 = 0; i2 < ny; i2++) {
00546                         t = float(i2 - ldpy);
00547                         const float yy = t*t + xx;
00548                         bool first = true;
00549                         for (int i3 = 0; i3 < nx; i3++) {
00550                                 t = float(i3 - ldpx);
00551                                 const float rc = t*t + yy;
00552                                 if (first) {
00553                                         // first pixel on this line
00554                                         if (rc > r) continue;
00555                                         first = false;
00556                                         nn++;
00557                                         if (ipcube != NULL) {
00558                                                 ipcube[nn].start = i3;
00559                                                 ipcube[nn].end   = i3;
00560                                                 ipcube[nn].loc[0] = i3 - ldpx;
00561                                                 ipcube[nn].loc[1] = i2 - ldpy;
00562                                                 ipcube[nn].loc[2] = i1 - ldpz;
00563                                                 //cout<<"  start  "<<i3<<"  i2  "<<i2<<"  i1 "<<i1<<endl;
00564                                         }
00565                                 } else {
00566                                         // second or later pixel on this line
00567                                         if (ipcube != NULL) {
00568                                                 if (rc <= r) ipcube[nn].end = i3;
00569                                         }
00570                                 }
00571                         }
00572                 }
00573         }
00574 }
00575 
00576 
00577 EMData *PawelProjector::project3d(EMData * image) const
00578 {
00579         if (!image) {
00580                 return 0;
00581         }
00582         int ri;
00583         int nx = image->get_xsize();
00584         int ny = image->get_ysize();
00585         int nz = image->get_zsize();
00586         int dim = Util::get_max(nx,ny,nz);
00587         if (nz == 1) {
00588                 LOGERR("The PawelProjector needs a volume!");
00589                 return 0;
00590         }
00591 
00592         Vec3i origin(0,0,0);
00593         // If a sensible origin isn't passed in, choose the middle of
00594         // the cube.
00595         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00596         else {origin[0] = nx/2;}
00597         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00598         else {origin[1] = ny/2;}
00599         if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00600         else {origin[2] = nz/2;}
00601 
00602         if (params.has_key("radius")) {ri = params["radius"];}
00603         else {ri = dim/2 - 1;}
00604 
00605         // Determine the number of rows (x-lines) within the radius
00606         int nn = -1;
00607         prepcubes(nx, ny, nz, ri, origin, nn);
00608         // nn is now the number of rows-1 within the radius
00609         // so we can create and fill the ipcubes
00610         IPCube* ipcube = new IPCube[nn+1];
00611         prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00612 
00613         Transform* rotation = params["transform"];
00614         int nangles = 0;
00615         vector<float> anglelist;
00616         // Do we have a list of angles?
00617         /*
00618         if (params.has_key("anglelist")) {
00619                 anglelist = params["anglelist"];
00620                 nangles = anglelist.size() / 3;
00621         } else {*/
00622 
00623                 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00624                 /*
00625                 Dict p = t3d->get_rotation("spider");
00626 
00627                 string angletype = "SPIDER";
00628                 float phi = p["phi"];
00629                 float theta = p["theta"];
00630                 float psi = p["psi"];
00631                 anglelist.push_back(phi);
00632                 anglelist.push_back(theta);
00633                 anglelist.push_back(psi);
00634                 */
00635                 nangles = 1;
00636         //}
00637 
00638         // initialize return object
00639         EMData* ret = new EMData();
00640         ret->set_size(nx, ny, nangles);
00641         ret->to_zero();
00642         //for (int i = 0 ; i <= nn; i++)  cout<<" IPCUBE "<<ipcube[i].start<<"  "<<ipcube[i].end<<"  "<<ipcube[i].loc[0]<<"  "<<ipcube[i].loc[1]<<"  "<<ipcube[i].loc[2]<<endl;
00643         // loop over sets of angles
00644         for (int ia = 0; ia < nangles; ia++) {
00645                 //int indx = 3*ia;
00646                 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
00647                 //Transform rotation(d);
00648                 if (2*(ri+1)+1 > dim) {
00649                         // Must check x and y boundaries
00650                         for (int i = 0 ; i <= nn; i++) {
00651                                 int k = ipcube[i].loc[1] + origin[1];
00652                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00653                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00654                                         // check for pixels out-of-bounds
00655                                         int iox = int(vb[0]);
00656                                         if ((iox >= 0) && (iox < nx-1)) {
00657                                                 int ioy = int(vb[1]);
00658                                                 if ((ioy >= 0) && (ioy < ny-1)) {
00659                                                         int ioz = int(vb[2]);
00660                                                         if ((ioz >= 0) && (ioz < nz-1)) {
00661                                                                 // real work for pixels in bounds
00662                                         //cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<endl;
00663                                         //cout<<"  TAKE"<<endl;
00664                                                                 float dx = vb[0] - iox;
00665                                                                 float dy = vb[1] - ioy;
00666                                                                 float dz = vb[2] - ioz;
00667                                                                 float a1 = (*image)(iox,ioy,ioz);
00668                                                                 float a2 = (*image)(iox+1,ioy,ioz) - a1;
00669                                                                 float a3 = (*image)(iox,ioy+1,ioz) - a1;
00670                                                                 float a4 = (*image)(iox,ioy,ioz+1) - a1;
00671                                                                 float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00672                                                                                 + (*image)(iox+1,ioy+1,ioz);
00673                                                                 float a61 = -(*image)(iox,ioy,ioz+1)
00674                                                                                 + (*image)(iox+1,ioy,ioz+1);
00675                                                                 float a6 = -a2 + a61;
00676                                                                 float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00677                                                                                 + (*image)(iox,ioy+1,ioz+1);
00678                                                                 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00679                                                                                 + (*image)(iox+1,ioy+1,ioz+1);
00680                                                                 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00681                                                                                 + (a7 + a8*dx)*dy)
00682                                                                                 + a3*dy + dx*(a2 + a5*dy);
00683                                                         } //else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLAATED Z "<<endl; }
00684                                                 }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED Y "<<endl; }
00685                                         }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED X "<<endl; }
00686                                         vb += rotation->get_matrix3_row(0);
00687                                 }
00688                         }
00689 
00690                 } else {
00691                         // No need to check x and y boundaries
00692                         for (int i = 0 ; i <= nn; i++) {
00693                                 int k = ipcube[i].loc[1] + origin[1];
00694                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00695                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00696                                         int iox = int(vb[0]);
00697                                         int ioy = int(vb[1]);
00698                                         int ioz = int(vb[2]);
00699                                         float dx = vb[0] - iox;
00700                                         float dy = vb[1] - ioy;
00701                                         float dz = vb[2] - ioz;
00702                                         float a1 = (*image)(iox,ioy,ioz);
00703                                         float a2 = (*image)(iox+1,ioy,ioz) - a1;
00704                                         float a3 = (*image)(iox,ioy+1,ioz) - a1;
00705                                         float a4 = (*image)(iox,ioy,ioz+1) - a1;
00706                                         float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00707                                                         + (*image)(iox+1,ioy+1,ioz);
00708                                         float a61 = -(*image)(iox,ioy,ioz+1)
00709                                                         + (*image)(iox+1,ioy,ioz+1);
00710                                         float a6 = -a2 + a61;
00711                                         float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00712                                                         + (*image)(iox,ioy+1,ioz+1);
00713                                         float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00714                                                         + (*image)(iox+1,ioy+1,ioz+1);
00715                                         (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00716                                                         + (a7 + a8*dx)*dy)
00717                                                         + a3*dy + dx*(a2 + a5*dy);
00718                                         vb += rotation->get_matrix3_row(0);
00719                                 }
00720                         }
00721                 }
00722         }
00723         ret->update();
00724         EMDeleteArray(ipcube);
00725         if(rotation) {delete rotation; rotation=0;}
00726         return ret;
00727 }
00728 
00729 // EMData *PawelProjector::project3d(EMData * image) const
00730 // {
00731 //      if (!image) {
00732 //              return 0;
00733 //      }
00734 //      int ri;
00735 //      int nx = image->get_xsize();
00736 //      int ny = image->get_ysize();
00737 //      int nz = image->get_zsize();
00738 //      int dim = Util::get_min(nx,ny,nz);
00739 //      if (nz == 1) {
00740 //              LOGERR("The PawelProjector needs a volume!");
00741 //              return 0;
00742 //      }
00743 //      Vec3i origin(0,0,0);
00744 //      // If a sensible origin isn't passed in, choose the middle of
00745 //      // the cube.
00746 //      if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00747 //      else {origin[0] = nx/2;}
00748 //      if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00749 //      else {origin[1] = ny/2;}
00750 //      if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00751 //      else {origin[2] = nz/2;}
00752 //
00753 //      if (params.has_key("radius")) {ri = params["radius"];}
00754 //      else {ri = dim/2 - 1;}
00755 //
00756 //      // Determine the number of rows (x-lines) within the radius
00757 //      int nn = -1;
00758 //      prepcubes(nx, ny, nz, ri, origin, nn);
00759 //      // nn is now the number of rows-1 within the radius
00760 //      // so we can create and fill the ipcubes
00761 //      IPCube* ipcube = new IPCube[nn+1];
00762 //      prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00763 //
00764 //      int nangles = 0;
00765 //      vector<float> anglelist;
00766 //      // Do we have a list of angles?
00767 //      if (params.has_key("anglelist")) {
00768 //              anglelist = params["anglelist"];
00769 //              nangles = anglelist.size() / 3;
00770 //      } else {
00771 //              Transform3D* t3d = params["t3d"];
00772 //              if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
00773 //
00774 //              // This part was modified by David Woolford -
00775 //              // Before this the code worked only for SPIDER and EMAN angles,
00776 //              // but the framework of the Transf      orm3D allows for a generic implementation
00777 //              // as specified here.
00778 //              Dict p = t3d->get_rotation(Transform3D::SPIDER);
00779 //
00780 //              string angletype = "SPIDER";
00781 //              float phi = p["phi"];
00782 //              float theta = p["theta"];
00783 //              float psi = p["psi"];
00784 //              anglelist.push_back(phi);
00785 //              anglelist.push_back(theta);
00786 //              anglelist.push_back(psi);
00787 //              nangles = 1;
00788 //      }
00789 //
00790 //      // End David Woolford modifications
00791 //
00792 //      // initialize return object
00793 //      EMData* ret = new EMData();
00794 //      ret->set_size(nx, ny, nangles);
00795 //      ret->to_zero();
00796 //
00797 //      // loop over sets of angles
00798 //      for (int ia = 0; ia < nangles; ia++) {
00799 //              int indx = 3*ia;
00800 //              Transform3D rotation(Transform3D::SPIDER, anglelist[indx],anglelist[indx+1], anglelist[indx+2]);
00801 //              if (2*(ri+1)+1 > dim) {
00802 //                      // Must check x and y boundaries
00803 //                      for (int i = 0 ; i <= nn; i++) {
00804 //                              int k = ipcube[i].loc[1] + origin[1];
00805 //                              Vec3f vb = ipcube[i].loc*rotation + origin;
00806 //                              for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00807 //                                      // check for pixels out-of-bounds
00808 //                                      int iox = int(vb[0]);
00809 //                                      if ((iox >= 0) && (iox < nx-1)) {
00810 //                                              int ioy = int(vb[1]);
00811 //                                              if ((ioy >= 0) && (ioy < ny-1)) {
00812 //                                                      int ioz = int(vb[2]);
00813 //                                                      if ((ioz >= 0) && (ioz < nz-1)) {
00814 //                                                              // real work for pixels in bounds
00815 //                                                              float dx = vb[0] - iox;
00816 //                                                              float dy = vb[1] - ioy;
00817 //                                                              float dz = vb[2] - ioz;
00818 //                                                              float a1 = (*image)(iox,ioy,ioz);
00819 //                                                              float a2 = (*image)(iox+1,ioy,ioz) - a1;
00820 //                                                              float a3 = (*image)(iox,ioy+1,ioz) - a1;
00821 //                                                              float a4 = (*image)(iox,ioy,ioz+1) - a1;
00822 //                                                              float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00823 //                                                                      + (*image)(iox+1,ioy+1,ioz);
00824 //                                                              float a61 = -(*image)(iox,ioy,ioz+1)
00825 //                                                                      + (*image)(iox+1,ioy,ioz+1);
00826 //                                                              float a6 = -a2 + a61;
00827 //                                                              float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00828 //                                                                      + (*image)(iox,ioy+1,ioz+1);
00829 //                                                              float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00830 //                                                                      + (*image)(iox+1,ioy+1,ioz+1);
00831 //                                                              (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00832 //                                                                      + (a7 + a8*dx)*dy)
00833 //                                                                      + a3*dy + dx*(a2 + a5*dy);
00834 //                                                      }
00835 //                                              }
00836 //                                      }
00837 //                                      vb += rotation.get_matrix3_row(0);
00838 //                              }
00839 //                      }
00840 //
00841 //              } else {
00842 //                      // No need to check x and y boundaries
00843 //                      for (int i = 0 ; i <= nn; i++) {
00844 //                              int k = ipcube[i].loc[1] + origin[1];
00845 //                              Vec3f vb = ipcube[i].loc*rotation + origin;
00846 //                              for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00847 //                                      int iox = int(vb[0]);
00848 //                                      int ioy = int(vb[1]);
00849 //                                      int ioz = int(vb[2]);
00850 //                                      float dx = vb[0] - iox;
00851 //                                      float dy = vb[1] - ioy;
00852 //                                      float dz = vb[2] - ioz;
00853 //                                      float a1 = (*image)(iox,ioy,ioz);
00854 //                                      float a2 = (*image)(iox+1,ioy,ioz) - a1;
00855 //                                      float a3 = (*image)(iox,ioy+1,ioz) - a1;
00856 //                                      float a4 = (*image)(iox,ioy,ioz+1) - a1;
00857 //                                      float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00858 //                                              + (*image)(iox+1,ioy+1,ioz);
00859 //                                      float a61 = -(*image)(iox,ioy,ioz+1)
00860 //                                              + (*image)(iox+1,ioy,ioz+1);
00861 //                                      float a6 = -a2 + a61;
00862 //                                      float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00863 //                                              + (*image)(iox,ioy+1,ioz+1);
00864 //                                      float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00865 //                                              + (*image)(iox+1,ioy+1,ioz+1);
00866 //                                      (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00867 //                                                      + (a7 + a8*dx)*dy)
00868 //                                              + a3*dy + dx*(a2 + a5*dy);
00869 //                                      vb += rotation.get_matrix3_row(0);
00870 //                              }
00871 //                      }
00872 //              }
00873 //      }
00874 //      ret->update();
00875 //      EMDeleteArray(ipcube);
00876 //      return ret;
00877 // }
00878 
00879 EMData *StandardProjector::project3d(EMData * image) const
00880 {
00881         Transform* t3d = params["transform"];
00882         if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00883 //      Dict p = t3d->get_rotation();
00884         if ( image->get_ndim() == 3 )
00885         {
00886 
00887 #ifdef EMAN2_USING_CUDA
00888                 if(EMData::usecuda == 1) {
00889                         if(!image->isrodataongpu()) image->copy_to_cudaro();
00890                         //cout << "CUDA PROJ" << endl;
00891                         Transform* t3d = params["transform"];
00892                         if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00893                         float * m = new float[12];
00894                         t3d->copy_matrix_into_array(m);
00895                         image->bindcudaarrayA(true);
00896                         //EMData* e = new EMData(0,0,image->get_xsize(),image->get_ysize(),1);
00897                         EMData *e = new EMData();
00898                         e->set_size_cuda(image->get_xsize(), image->get_ysize(), 1);
00899                         e->rw_alloc();
00900                         standard_project(m,e->getcudarwdata(), e->get_xsize(), e->get_ysize());
00901                         image->unbindcudaarryA();
00902                         delete [] m;
00903                 
00904                         e->update();
00905                         e->set_attr("xform.projection",t3d);
00906                         e->set_attr("apix_x",(float)image->get_attr("apix_x"));
00907                         e->set_attr("apix_y",(float)image->get_attr("apix_y"));
00908                         e->set_attr("apix_z",(float)image->get_attr("apix_z"));
00909                         //e_>copy_from_device();
00910                         if(t3d) {delete t3d; t3d=0;}
00911                         return e;
00912                 }
00913 #endif
00914                 int nx = image->get_xsize();
00915                 int ny = image->get_ysize();
00916                 int nz = image->get_zsize();
00917 
00918 //              Transform3D r(Transform3D::EMAN, az, alt, phi);
00919                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
00920                 int xy = nx * ny;
00921 
00922                 EMData *proj = new EMData();
00923                 proj->set_size(nx, ny, 1);
00924 
00925                 Vec3i offset(nx/2,ny/2,nz/2);
00926 
00927                 float *sdata = image->get_data();
00928                 float *ddata = proj->get_data();
00929                 for (int k = -nz / 2; k < nz - nz / 2; k++) {
00930                         int l = 0;
00931                         for (int j = -ny / 2; j < ny - ny / 2; j++) {
00932                                 ddata[l]=0;
00933                                 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00934 
00935                                         Vec3f coord(i,j,k);
00936                                         Vec3f soln = r*coord;
00937                                         soln += offset;
00938 
00942 //                                      printf(" ");
00943 
00944                                         float x2 = soln[0];
00945                                         float y2 = soln[1];
00946                                         float z2 = soln[2];
00947 
00948                                         float x = (float)Util::fast_floor(x2);
00949                                         float y = (float)Util::fast_floor(y2);
00950                                         float z = (float)Util::fast_floor(z2);
00951 
00952                                         float t = x2 - x;
00953                                         float u = y2 - y;
00954                                         float v = z2 - z;
00955 
00956                                         size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
00957 // 
00958                                         if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
00959                                         if      (x2 > (nx-1) || y2  > (ny-1) || z2 > (nz-1) ) continue;
00960 
00961                                         if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
00962                                                 ddata[l] +=
00963                                                                 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
00964                                                                 sdata[ii + nx + 1], sdata[ii + xy],     sdata[ii + xy + 1], sdata[ii + xy + nx],
00965                                                                 sdata[ii + xy + nx + 1], t, u, v);
00966                                         }
00967                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
00968                                                 ddata[l] += sdata[ii];
00969                                         }
00970                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
00971                                                 ddata[l] +=     Util::linear_interpolate(sdata[ii], sdata[ii + xy],v);
00972                                         }
00973                                         else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
00974                                                 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + nx],u);
00975                                         }
00976                                         else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
00977                                                 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + 1],t);
00978                                         }
00979                                         else if ( x2 == (nx - 1) ) {
00980                                                 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v);
00981                                         }
00982                                         else if ( y2 == (ny - 1) ) {
00983                                                 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v);
00984                                         }
00985                                         else if ( z2 == (nz - 1) ) {
00986                                                 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u);
00987                                         }
00988                                 }
00989                         }
00990                 }
00991                 proj->update();
00992                 proj->set_attr("xform.projection",t3d);
00993                 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
00994                 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
00995                 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
00996                 
00997                 if(t3d) {delete t3d; t3d=0;}
00998                 return proj;
00999         }
01000         else if ( image->get_ndim() == 2 ) {
01001 
01002                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
01003 
01004                 int nx = image->get_xsize();
01005                 int ny = image->get_ysize();
01006 
01007                 EMData *proj = new EMData();
01008                 proj->set_size(nx, 1, 1);
01009                 proj->to_zero();
01010 
01011                 float *sdata = image->get_data();
01012                 float *ddata = proj->get_data();
01013 
01014                 Vec2f offset(nx/2,ny/2);
01015                 for (int j = -ny / 2; j < ny - ny / 2; j++) { // j represents a column of pixels in the direction of the angle
01016                         int l = 0;
01017                         for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
01018 
01019                                 Vec2f coord(i,j);
01020                                 Vec2f soln = r*coord;
01021                                 soln += offset;
01022 
01023                                 float x2 = soln[0];
01024                                 float y2 = soln[1];
01025 
01026                                 float x = (float)Util::fast_floor(x2);
01027                                 float y = (float)Util::fast_floor(y2);
01028 
01029                                 int ii = (int) (x + y * nx);
01030                                 float u = x2 - x;
01031                                 float v = y2 - y;
01032 
01033                                 if (x2 < 0 || y2 < 0 ) continue;
01034                                 if      (x2 > (nx-1) || y2  > (ny-1) ) continue;
01035 
01036                                 if (  x2 < (nx - 1) && y2 < (ny - 1) ) {
01037                                         ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v);
01038                                 }
01039                                 else if (x2 == (nx-1) && y2 == (ny-1) ) {
01040                                         ddata[l] += sdata[ii];
01041                                 }
01042                                 else if (x2 == (nx-1) ) {
01043                                         ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + nx], v);
01044                                 }
01045                                 else if (y2 == (ny-1) ) {
01046                                         ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + 1], u);
01047                                 }
01048                         }
01049                 }
01050                 proj->set_attr("xform.projection",t3d);
01051                 proj->update();
01052                 if(t3d) {delete t3d; t3d=0;}
01053                 return proj;
01054         }
01055         else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
01056 }
01057 
01058 EMData *MaxValProjector::project3d(EMData * image) const
01059 {
01060         Transform* t3d = params["transform"];
01061         if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
01062 //      Dict p = t3d->get_rotation();
01063         if ( image->get_ndim() == 3 )
01064         {
01065 
01066                 int nx = image->get_xsize();
01067                 int ny = image->get_ysize();
01068                 int nz = image->get_zsize();
01069 
01070 //              Transform3D r(Transform3D::EMAN, az, alt, phi);
01071                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
01072                 int xy = nx * ny;
01073 
01074                 EMData *proj = new EMData();
01075                 proj->set_size(nx, ny, 1);
01076 
01077                 Vec3i offset(nx/2,ny/2,nz/2);
01078 
01079                 float *sdata = image->get_data();
01080                 float *ddata = proj->get_data();
01081                 for (int k = -nz / 2; k < nz - nz / 2; k++) {
01082                         int l = 0;
01083                         for (int j = -ny / 2; j < ny - ny / 2; j++) {
01084                                 ddata[l]=0;
01085                                 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
01086 
01087                                         Vec3f coord(i,j,k);
01088                                         Vec3f soln = r*coord;
01089                                         soln += offset;
01090 
01094 //                                      printf(" ");
01095 
01096                                         float x2 = soln[0];
01097                                         float y2 = soln[1];
01098                                         float z2 = soln[2];
01099 
01100                                         float x = (float)Util::fast_floor(x2);
01101                                         float y = (float)Util::fast_floor(y2);
01102                                         float z = (float)Util::fast_floor(z2);
01103 
01104                                         float t = x2 - x;
01105                                         float u = y2 - y;
01106                                         float v = z2 - z;
01107 
01108                                         size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
01109 // 
01110                                         if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
01111                                         if      (x2 > (nx-1) || y2  > (ny-1) || z2 > (nz-1) ) continue;
01112 
01113                                         if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
01114                                                 ddata[l] = Util::get_max(ddata[l],
01115                                                                 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
01116                                                                 sdata[ii + nx + 1], sdata[ii + xy],     sdata[ii + xy + 1], sdata[ii + xy + nx],
01117                                                                 sdata[ii + xy + nx + 1], t, u, v));
01118                                         }
01119                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
01120                                                 ddata[l] = Util::get_max(ddata[l],sdata[ii]);
01121                                         }
01122                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
01123                                                 ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + xy],v));
01124                                         }
01125                                         else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
01126                                                 ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + nx],u));
01127                                         }
01128                                         else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
01129                                                 ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + 1],t));
01130                                         }
01131                                         else if ( x2 == (nx - 1) ) {
01132                                                 ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v));
01133                                         }
01134                                         else if ( y2 == (ny - 1) ) {
01135                                                 ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v));
01136                                         }
01137                                         else if ( z2 == (nz - 1) ) {
01138                                                 ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u));
01139                                         }
01140                                 }
01141                         }
01142                 }
01143                 proj->update();
01144                 proj->set_attr("xform.projection",t3d);
01145                 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
01146                 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
01147                 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
01148                 
01149                 if(t3d) {delete t3d; t3d=0;}
01150                 return proj;
01151         }
01152         else if ( image->get_ndim() == 2 ) {
01153 
01154                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
01155 
01156                 int nx = image->get_xsize();
01157                 int ny = image->get_ysize();
01158 
01159                 EMData *proj = new EMData();
01160                 proj->set_size(nx, 1, 1);
01161                 proj->to_zero();
01162 
01163                 float *sdata = image->get_data();
01164                 float *ddata = proj->get_data();
01165 
01166                 Vec2f offset(nx/2,ny/2);
01167                 for (int j = -ny / 2; j < ny - ny / 2; j++) { // j represents a column of pixels in the direction of the angle
01168                         int l = 0;
01169                         for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
01170 
01171                                 Vec2f coord(i,j);
01172                                 Vec2f soln = r*coord;
01173                                 soln += offset;
01174 
01175                                 float x2 = soln[0];
01176                                 float y2 = soln[1];
01177 
01178                                 float x = (float)Util::fast_floor(x2);
01179                                 float y = (float)Util::fast_floor(y2);
01180 
01181                                 int ii = (int) (x + y * nx);
01182                                 float u = x2 - x;
01183                                 float v = y2 - y;
01184 
01185                                 if (x2 < 0 || y2 < 0 ) continue;
01186                                 if      (x2 > (nx-1) || y2  > (ny-1) ) continue;
01187 
01188                                 if (  x2 < (nx - 1) && y2 < (ny - 1) ) {
01189                                         ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v));
01190                                 }
01191                                 else if (x2 == (nx-1) && y2 == (ny-1) ) {
01192                                         ddata[l] =      Util::get_max(ddata[l],sdata[ii]);
01193                                 }
01194                                 else if (x2 == (nx-1) ) {
01195                                         ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii],sdata[ii + nx], v));
01196                                 }
01197                                 else if (y2 == (ny-1) ) {
01198                                         ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii],sdata[ii + 1], u));
01199                                 }
01200                         }
01201                 }
01202                 proj->set_attr("xform.projection",t3d);
01203                 proj->update();
01204                 if(t3d) {delete t3d; t3d=0;}
01205                 return proj;
01206         }
01207         else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
01208 }
01209 
01210 
01211 // EMData *FourierGriddingProjector::project3d(EMData * image) const
01212 // {
01213 //      if (!image) {
01214 //              return 0;
01215 //      }
01216 //      if (3 != image->get_ndim())
01217 //              throw ImageDimensionException(
01218 //                              "FourierGriddingProjector needs a 3-D volume");
01219 //      if (image->is_complex())
01220 //              throw ImageFormatException(
01221 //                              "FourierGriddingProjector requires a real volume");
01222 //      const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
01223 //      const int nx = image->get_xsize();
01224 //      const int ny = image->get_ysize();
01225 //      const int nz = image->get_zsize();
01226 //      if (nx != ny || nx != nz)
01227 //              throw ImageDimensionException(
01228 //                              "FourierGriddingProjector requires nx==ny==nz");
01229 //      const int m = Util::get_min(nx,ny,nz);
01230 //      const int n = m*npad;
01231 //
01232 //      int K = params["kb_K"];
01233 //      if ( K == 0 ) K = 6;
01234 //      float alpha = params["kb_alpha"];
01235 //      if ( alpha == 0 ) alpha = 1.25;
01236 //      Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
01237 //
01238 //      // divide out gridding weights
01239 //      EMData* tmpImage = image->copy();
01240 //      tmpImage->divkbsinh(kb);
01241 //      // pad and center volume, then FFT and multiply by (-1)**(i+j+k)
01242 //      //EMData* imgft = tmpImage->pad_fft(npad);
01243 //      //imgft->center_padded();
01244 //      EMData* imgft = tmpImage->norm_pad(false, npad);
01245 //      imgft->do_fft_inplace();
01246 //      imgft->center_origin_fft();
01247 //      imgft->fft_shuffle();
01248 //      delete tmpImage;
01249 //
01250 //      // Do we have a list of angles?
01251 //      int nangles = 0;
01252 //      vector<float> anglelist;
01253 //      // Do we have a list of angles?
01254 //      if (params.has_key("anglelist")) {
01255 //              anglelist = params["anglelist"];
01256 //              nangles = anglelist.size() / 3;
01257 //      } else {
01258 //              // This part was modified by David Woolford -
01259 //              // Before this the code worked only for SPIDER and EMAN angles,
01260 //              // but the framework of the Transform3D allows for a generic implementation
01261 //              // as specified here.
01262 //              Transform3D* t3d = params["t3d"];
01263 //              if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
01264 //              Dict p = t3d->get_rotation(Transform3D::SPIDER);
01265 //
01266 //              string angletype = "SPIDER";
01267 //              float phi = p["phi"];
01268 //              float theta = p["theta"];
01269 //              float psi = p["psi"];
01270 //              anglelist.push_back(phi);
01271 //              anglelist.push_back(theta);
01272 //              anglelist.push_back(psi);
01273 //              nangles = 1;
01274 //      }
01275 //
01276 //      // End David Woolford modifications
01277 //
01278 //      // initialize return object
01279 //      EMData* ret = new EMData();
01280 //      ret->set_size(nx, ny, nangles);
01281 //      ret->to_zero();
01282 //      // loop over sets of angles
01283 //      for (int ia = 0; ia < nangles; ia++) {
01284 //              int indx = 3*ia;
01285 //              Transform3D tf(Transform3D::SPIDER, anglelist[indx],anglelist[indx+1],anglelist[indx+2]);
01286 //              EMData* proj = imgft->extractplane(tf, kb);
01287 //              if (proj->is_shuffled()) proj->fft_shuffle();
01288 //              proj->center_origin_fft();
01289 //              proj->do_ift_inplace();
01290 //              EMData* winproj = proj->window_center(m);
01291 //              delete proj;
01292 //              for (int iy=0; iy < ny; iy++)
01293 //                      for (int ix=0; ix < nx; ix++)
01294 //                              (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
01295 //              delete winproj;
01296 //      }
01297 //      delete imgft;
01298 //      ret->update();
01299 //
01300 //      return ret;
01301 // }
01302 
01303 
01304 EMData *FourierGriddingProjector::project3d(EMData * image) const
01305 {
01306         if (!image) {
01307                 return 0;
01308         }
01309         if (3 != image->get_ndim())
01310                 throw ImageDimensionException(
01311                                                                           "FourierGriddingProjector needs a 3-D volume");
01312         if (image->is_complex())
01313                 throw ImageFormatException(
01314                                                                    "FourierGriddingProjector requires a real volume");
01315         const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
01316         const int nx = image->get_xsize();
01317         const int ny = image->get_ysize();
01318         const int nz = image->get_zsize();
01319         if (nx != ny || nx != nz)
01320                 throw ImageDimensionException(
01321                                                                           "FourierGriddingProjector requires nx==ny==nz");
01322         const int m = Util::get_min(nx,ny,nz);
01323         const int n = m*npad;
01324 
01325         int K = params["kb_K"];
01326         if ( K == 0 ) K = 6;
01327         float alpha = params["kb_alpha"];
01328         if ( alpha == 0 ) alpha = 1.25;
01329         Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
01330 
01331         // divide out gridding weights
01332         EMData* tmpImage = image->copy();
01333         tmpImage->divkbsinh(kb);
01334         // pad and center volume, then FFT and multiply by (-1)**(i+j+k)
01335         //EMData* imgft = tmpImage->pad_fft(npad);
01336         //imgft->center_padded();
01337         EMData* imgft = tmpImage->norm_pad(false, npad);
01338         imgft->do_fft_inplace();
01339         imgft->center_origin_fft();
01340         imgft->fft_shuffle();
01341         delete tmpImage;
01342 
01343         // Do we have a list of angles?
01344         int nangles = 0;
01345         vector<float> anglelist;
01346         // Do we have a list of angles?
01347         if (params.has_key("anglelist")) {
01348                 anglelist = params["anglelist"];
01349                 nangles = anglelist.size() / 3;
01350         } else {
01351                 // This part was modified by David Woolford -
01352                 // Before this the code worked only for SPIDER and EMAN angles,
01353                 // but the framework of the Transform3D allows for a generic implementation
01354                 // as specified here.
01355                 Transform* t3d = params["transform"];
01356                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01357                 Dict p = t3d->get_rotation("spider");
01358 
01359                 string angletype = "SPIDER";
01360                 float phi = p["phi"];
01361                 float theta = p["theta"];
01362                 float psi = p["psi"];
01363                 anglelist.push_back(phi);
01364                 anglelist.push_back(theta);
01365                 anglelist.push_back(psi);
01366                 nangles = 1;
01367                 if(t3d) {delete t3d; t3d=0;}
01368         }
01369 
01370         // End David Woolford modifications
01371 
01372         // initialize return object
01373         EMData* ret = new EMData();
01374         ret->set_size(nx, ny, nangles);
01375         ret->to_zero();
01376         // loop over sets of angles
01377         for (int ia = 0; ia < nangles; ia++) {
01378                 int indx = 3*ia;
01379                 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
01380                 Transform tf(d);
01381                 EMData* proj = imgft->extract_plane(tf, kb);
01382                 if (proj->is_shuffled()) proj->fft_shuffle();
01383                 proj->center_origin_fft();
01384                 proj->do_ift_inplace();
01385                 EMData* winproj = proj->window_center(m);
01386                 delete proj;
01387                 for (int iy=0; iy < ny; iy++)
01388                         for (int ix=0; ix < nx; ix++)
01389                                 (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
01390                 delete winproj;
01391         }
01392         delete imgft;
01393 
01394         if (!params.has_key("anglelist")) {
01395                 Transform* t3d = params["transform"];
01396                 ret->set_attr("xform.projection",t3d);
01397                 if(t3d) {delete t3d; t3d=0;}
01398         }
01399         ret->update();
01400         return ret;
01401 }
01402 
01403 // BEGIN Chao projectors and backprojector addition (04/25/06)
01404 int ChaoProjector::getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) const
01405 /*
01406    purpose: count the number of voxels within a sphere centered
01407             at origin and with a radius ri.
01408 
01409      input:
01410      volsize contains the size information (nx,ny,nz) about the volume
01411      ri      radius of the object embedded in the cube.
01412      origin  coordinates for the center of the volume
01413 
01414      output:
01415      nnz    total number of voxels within the sphere (of radius ri)
01416      nrays  number of rays in z-direction.
01417 */
01418 {
01419         int  ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
01420         int  ftm=0, status = 0;
01421 
01422         r2    = ri*ri;
01423         *nnz  = 0;
01424         *nrays = 0;
01425         int nx = (int)volsize[0];
01426         int ny = (int)volsize[1];
01427         int nz = (int)volsize[2];
01428 
01429         int xcent = (int)origin[0];
01430         int ycent = (int)origin[1];
01431         int zcent = (int)origin[2];
01432 
01433         // need to add some error checking
01434         for (ix = 1; ix <=nx; ix++) {
01435             xs  = ix-xcent;
01436             xx  = xs*xs;
01437             for (iy = 1; iy <= ny; iy++) {
01438                 ys = iy-ycent;
01439                 yy = ys*ys;
01440                 ftm = 1;
01441                 for (iz = 1; iz <= nz; iz++) {
01442                     zs = iz-zcent;
01443                     zz = zs*zs;
01444                     rs = xx + yy + zz;
01445                     if (rs <= r2) {
01446                         (*nnz)++;
01447                         if (ftm) {
01448                            (*nrays)++;
01449                            ftm = 0;
01450                         }
01451                     }
01452                 }
01453             } // end for iy
01454         } // end for ix
01455         return status;
01456 }
01457 
01458 #define cube(i,j,k) cube[ ((k-1)*ny + j-1)*nx + i-1 ]
01459 #define sphere(i)   sphere[(i)-1]
01460 #define cord(i,j)   cord[((j)-1)*3 + (i) -1]
01461 #define ptrs(i)     ptrs[(i)-1]
01462 #define dm(i)       dm[(i)-1]
01463 
01464 int ChaoProjector:: cb2sph(float *cube, Vec3i volsize, int    ri, Vec3i origin,
01465                            int    nnz0, int     *ptrs, int *cord, float *sphere) const
01466 {
01467     int    xs, ys, zs, xx, yy, zz, rs, r2;
01468     int    ix, iy, iz, jnz, nnz, nrays;
01469     int    ftm = 0, status = 0;
01470 
01471     int xcent = (int)origin[0];
01472     int ycent = (int)origin[1];
01473     int zcent = (int)origin[2];
01474 
01475     int nx = (int)volsize[0];
01476     int ny = (int)volsize[1];
01477     int nz = (int)volsize[2];
01478 
01479     r2      = ri*ri;
01480     nnz     = 0;
01481     nrays    = 0;
01482     ptrs(1) = 1;
01483 
01484     for (ix = 1; ix <= nx; ix++) {
01485        xs  = ix-xcent;
01486        xx  = xs*xs;
01487        for ( iy = 1; iy <= ny; iy++ ) {
01488            ys = iy-ycent;
01489            yy = ys*ys;
01490            jnz = 0;
01491 
01492            ftm = 1;
01493            // not the most efficient implementation
01494            for (iz = 1; iz <= nz; iz++) {
01495                zs = iz-zcent;
01496                zz = zs*zs;
01497                rs = xx + yy + zz;
01498                if (rs <= r2) {
01499                   jnz++;
01500                   nnz++;
01501                   sphere(nnz) = cube(iz, iy, ix);
01502 
01503                   //  record the coordinates of the first nonzero ===
01504                   if (ftm) {
01505                      nrays++;
01506                      cord(1,nrays) = iz;
01507                      cord(2,nrays) = iy;
01508                      cord(3,nrays) = ix;
01509                      ftm = 0;
01510                   }
01511                }
01512             } // end for (iz..)
01513             if (jnz > 0) {
01514                 ptrs(nrays+1) = ptrs(nrays) + jnz;
01515             }  // endif (jnz)
01516        } // end for iy
01517     } // end for ix
01518     if (nnz != nnz0) status = -1;
01519     return status;
01520 }
01521 
01522 // decompress sphere into a cube
01523 int ChaoProjector::sph2cb(float *sphere, Vec3i volsize, int  nrays, int    ri,
01524                           int      nnz0, int     *ptrs, int  *cord, float *cube) const
01525 {
01526     int       status=0;
01527     int       r2, i, j, ix, iy, iz,  nnz;
01528 
01529     int nx = (int)volsize[0];
01530     int ny = (int)volsize[1];
01531     // int nz = (int)volsize[2];
01532 
01533     r2      = ri*ri;
01534     nnz     = 0;
01535     ptrs(1) = 1;
01536 
01537     // no need to initialize
01538     // for (i = 0; i<nx*ny*nz; i++) cube[i]=0.0;
01539 
01540     nnz = 0;
01541     for (j = 1; j <= nrays; j++) {
01542        iz = cord(1,j);
01543        iy = cord(2,j);
01544        ix = cord(3,j);
01545        for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
01546            nnz++;
01547            cube(iz,iy,ix) = sphere(nnz);
01548        }
01549     }
01550     if (nnz != nnz0) status = -1;
01551     return status;
01552 }
01553 
01554 #define x(i)        x[(i)-1]
01555 #define y(i,j)      y[(j-1)*nx + i - 1]
01556 
01557 // project from 3D to 2D (single image)
01558 int ChaoProjector::fwdpj3(Vec3i volsize, int nrays, int      , float *dm,
01559                           Vec3i  origin, int    ri, int *ptrs, int *cord,
01560                           float      *x, float  *y) const
01561 {
01562     /*
01563         purpose:  y <--- proj(x)
01564         input  :  volsize  the size (nx,ny,nz) of the volume
01565                   nrays    number of rays within the compact spherical
01566                            representation
01567                   nnz      number of voxels within the sphere
01568                   dm       an array of size 9 storing transformation
01569                            associated with the projection direction
01570                   origin   coordinates of the center of the volume
01571                   ri       radius of the sphere
01572                   ptrs     the beginning address of each ray
01573                   cord     the coordinates of the first point in each ray
01574                   x        3d input volume
01575                   y        2d output image
01576     */
01577 
01578     int    iqx, iqy, i, j, xc, yc, zc;
01579     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
01580     int    status = 0;
01581 
01582     int xcent = origin[0];
01583     int ycent = origin[1];
01584     int zcent = origin[2];
01585 
01586     int nx = volsize[0];
01587 
01588     dm1 = dm(1);
01589     dm4 = dm(4);
01590 
01591     if ( nx > 2*ri ) {
01592         for (i = 1; i <= nrays; i++) {
01593             zc = cord(1,i)-zcent;
01594             yc = cord(2,i)-ycent;
01595             xc = cord(3,i)-xcent;
01596 
01597             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01598             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01599 
01600             for (j = ptrs(i); j< ptrs(i+1); j++) {
01601                iqx = ifix(xb);
01602                iqy = ifix(yb);
01603 
01604                ct   = x(j);
01605                dipx =  xb - (float)(iqx);
01606                dipy = (yb - (float)(iqy)) * ct;
01607 
01608                dipy1m = ct - dipy;
01609                dipx1m = 1.0f - dipx;
01610 
01611                y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
01612                y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m;
01613                y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;
01614                y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
01615 
01616                xb += dm1;
01617                yb += dm4;
01618            }
01619         }
01620     }
01621     else {
01622         fprintf(stderr, " nx must be greater than 2*ri\n");
01623         exit(1);
01624     }
01625     return status;
01626 }
01627 #undef x
01628 #undef y
01629 
01630 #define y(i)        y[(i)-1]
01631 #define x(i,j)      x[((j)-1)*nx + (i) - 1]
01632 
01633 // backproject from 2D to 3D for a single image
01634 int ChaoProjector::bckpj3(Vec3i volsize, int nrays, int      , float *dm,
01635                           Vec3i  origin, int    ri, int *ptrs, int *cord,
01636                           float      *x, float *y) const
01637 {
01638     int       i, j, iqx,iqy, xc, yc, zc;
01639     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
01640     int       status = 0;
01641 
01642     int xcent = origin[0];
01643     int ycent = origin[1];
01644     int zcent = origin[2];
01645 
01646     int nx = volsize[0];
01647 
01648     if ( nx > 2*ri) {
01649         for (i = 1; i <= nrays; i++) {
01650             zc = cord(1,i) - zcent;
01651             yc = cord(2,i) - ycent;
01652             xc = cord(3,i) - xcent;
01653 
01654             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01655             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01656 
01657             for (j = ptrs(i); j <ptrs(i+1); j++) {
01658                 iqx = ifix((float)(xb));
01659                 iqy = ifix((float)(yb));
01660 
01661                 dx = xb - (float)(iqx);
01662                 dy = yb - (float)(iqy);
01663                 dx1m = 1.0f - dx;
01664                 dy1m = 1.0f - dy;
01665                 dxdy = dx*dy;
01666 /*
01667 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
01668 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
01669 c     &                     + dx  *dy1m*x(iqx+1, iqy)
01670 c     &                     + dx  *dy  *x(iqx+1, iqy+1)
01671 c
01672 c              --- faster version of the above commented out
01673 c                  code (derived by summing the following table
01674 c                  of coefficients along  the colunms) ---
01675 c
01676 c                        1         dx        dy      dxdy
01677 c                     ------   --------  --------  -------
01678 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)
01679 c                                        x(i,j+1) -x(i,j+1)
01680 c                              x(i+1,j)           -x(i+1,j)
01681 c                                                x(i+1,j+1)
01682 c
01683 */
01684                y(j) += x(iqx,iqy)
01685                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
01686                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
01687                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1)
01688                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
01689 
01690                xb += dm(1);
01691                yb += dm(4);
01692             } // end for j
01693         } // end for i
01694      }
01695     else {
01696         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
01697     }
01698 
01699     return status;
01700 }
01701 
01702 #undef x
01703 #undef y
01704 #undef dm
01705 
01706 // funny F90 style strange rounding function
01707 int ChaoProjector::ifix(float a) const
01708 {
01709     int ia;
01710 
01711     if (a>=0) {
01712        ia = (int)floor(a);
01713     }
01714     else {
01715        ia = (int)ceil(a);
01716     }
01717     return ia;
01718 }
01719 
01720 #define dm(i,j)          dm[((j)-1)*9 + (i) -1]
01721 #define anglelist(i,j)   anglelist[((j)-1)*3 + (i) - 1]
01722 
01723 // SPIDER stype transformation
01724 void ChaoProjector::setdm(vector<float> anglelist, string const , float *dm) const
01725 { // convert Euler angles to transformations, dm is an 9 by nangles array
01726 
01727         float  psi, theta, phi;
01728         double cthe, sthe, cpsi, spsi, cphi, sphi;
01729         int    j;
01730 
01731         int nangles = anglelist.size() / 3;
01732 
01733         // now convert all angles
01734         for (j = 1; j <= nangles; j++) {
01735                 phi   = static_cast<float>(anglelist(1,j)*dgr_to_rad);
01736                 theta = static_cast<float>(anglelist(2,j)*dgr_to_rad);
01737                 psi   = static_cast<float>(anglelist(3,j)*dgr_to_rad);
01738 
01739                 //              cout << phi << " " << theta << " " << psi << endl;
01740                 cthe  = cos(theta);
01741                 sthe  = sin(theta);
01742                 cpsi  = cos(psi);
01743                 spsi  = sin(psi);
01744                 cphi  = cos(phi);
01745                 sphi  = sin(phi);
01746 
01747                 dm(1,j)=static_cast<float>(cphi*cthe*cpsi-sphi*spsi);
01748                 dm(2,j)=static_cast<float>(sphi*cthe*cpsi+cphi*spsi);
01749                 dm(3,j)=static_cast<float>(-sthe*cpsi);
01750                 dm(4,j)=static_cast<float>(-cphi*cthe*spsi-sphi*cpsi);
01751                 dm(5,j)=static_cast<float>(-sphi*cthe*spsi+cphi*cpsi);
01752                 dm(6,j)=static_cast<float>(sthe*spsi);
01753                 dm(7,j)=static_cast<float>(sthe*cphi);
01754                 dm(8,j)=static_cast<float>(sthe*sphi);
01755                 dm(9,j)=static_cast<float>(cthe);
01756         }
01757 }
01758 #undef anglelist
01759 
01760 #define images(i,j,k) images[ ((k-1)*nyvol + j-1)*nxvol + i-1 ]
01761 
01762 EMData *ChaoProjector::project3d(EMData * vol) const
01763 {
01764 
01765         int nrays, nnz, status, j;
01766         float *dm;
01767         int   *ptrs, *cord;
01768         float *sphere, *images;
01769 
01770         int nxvol = vol->get_xsize();
01771         int nyvol = vol->get_ysize();
01772         int nzvol = vol->get_zsize();
01773         Vec3i volsize(nxvol,nyvol,nzvol);
01774 
01775         int dim = Util::get_min(nxvol,nyvol,nzvol);
01776         if (nzvol == 1) {
01777                 LOGERR("The ChaoProjector needs a volume!");
01778                 return 0;
01779         }
01780         Vec3i origin(0,0,0);
01781         // If a sensible origin isn't passed in, choose the middle of
01782         // the cube.
01783         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01784         else {origin[0] = nxvol/2+1;}
01785         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01786         else {origin[1] = nyvol/2+1;}
01787         if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
01788         else {origin[2] = nzvol/2+1;}
01789 
01790         int ri;
01791         if (params.has_key("radius")) {ri = params["radius"];}
01792         else {ri = dim/2 - 1;}
01793 
01794         // retrieve the voxel values
01795         float *cube = vol->get_data();
01796 
01797         // count the number of voxels within a sphere centered at icent,
01798         // with radius ri
01799         status = getnnz(volsize, ri, origin, &nrays, &nnz);
01800         // need to check status...
01801 
01802         // convert from cube to sphere
01803         sphere = new float[nnz];
01804         ptrs   = new int[nrays+1];
01805         cord   = new int[3*nrays];
01806         if (sphere == NULL || ptrs == NULL || cord == NULL) {
01807                 fprintf(stderr,"ChaoProjector::project3d, failed to allocate!\n");
01808                 exit(1);
01809         }
01810         for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01811         for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01812         for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01813 
01814         status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01815         // check status
01816 
01817         int nangles = 0;
01818         vector<float> anglelist;
01819         string angletype = "SPIDER";
01820         // Do we have a list of angles?
01821         if (params.has_key("anglelist")) {
01822                 anglelist = params["anglelist"];
01823                 nangles = anglelist.size() / 3;
01824         } else {
01825                 Transform* t3d = params["transform"];
01826                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01827                 // This part was modified by David Woolford -
01828                 // Before this the code worked only for SPIDER and EMAN angles,
01829                 // but the framework of the Transform3D allows for a generic implementation
01830                 // as specified here.
01831                 Dict p = t3d->get_rotation("spider");
01832                 if(t3d) {delete t3d; t3d=0;}
01833 
01834                 float phi   = p["phi"];
01835                 float theta = p["theta"];
01836                 float psi   = p["psi"];
01837                 anglelist.push_back(phi);
01838                 anglelist.push_back(theta);
01839                 anglelist.push_back(psi);
01840                 nangles = 1;
01841         }
01842         // End David Woolford modifications
01843 
01844         dm = new float[nangles*9];
01845         setdm(anglelist, angletype, dm);
01846 
01847                 // return images
01848         EMData *ret = new EMData();
01849         ret->set_size(nxvol, nyvol, nangles);
01850         ret->set_complex(false);
01851         ret->set_ri(true);
01852 
01853         images = ret->get_data();
01854 
01855         for (j = 1; j <= nangles; j++) {
01856                 status = fwdpj3(volsize, nrays, nnz   , &dm(1,j), origin, ri,
01857                                                 ptrs   ,  cord, sphere, &images(1,1,j));
01858         // check status?
01859         }
01860 
01861         // deallocate all temporary work space
01862         EMDeleteArray(dm);
01863         EMDeleteArray(ptrs);
01864         EMDeleteArray(cord);
01865         EMDeleteArray(sphere);
01866 
01867         if (!params.has_key("anglelist")) {
01868                 Transform* t3d = params["transform"];
01869                 ret->set_attr("xform.projection",t3d);
01870                 if(t3d) {delete t3d; t3d=0;}
01871         }
01872         ret->update();
01873         return ret;
01874 }
01875 
01876 
01877 #undef images
01878 
01879 #define images(i,j,k) images[ ((k)-1)*nximg*nyimg + ((j)-1)*nximg + (i)-1 ]
01880 // backproject from 2D to 3D (multiple images)
01881 EMData *ChaoProjector::backproject3d(EMData * imagestack) const
01882 {
01883         int nrays, nnz, status, j;
01884         float *dm;
01885         int   *ptrs, *cord;
01886         float *sphere, *images, *cube;
01887 
01888         int nximg   = imagestack->get_xsize();
01889         int nyimg   = imagestack->get_ysize();
01890         int nslices = imagestack->get_zsize();
01891 
01892         int dim = Util::get_min(nximg,nyimg);
01893         Vec3i volsize(nximg,nyimg,dim);
01894 
01895         Vec3i origin(0,0,0);
01896         // If a sensible origin isn't passed in, choose the middle of
01897         // the cube.
01898         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01899         else {origin[0] = nximg/2+1;}
01900         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01901         else {origin[1] = nyimg/2+1;}
01902         if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01903         else {origin[2] = dim/2+1;}
01904 
01905         int ri;
01906         if (params.has_key("radius")) {ri = params["radius"];}
01907         else {ri = dim/2 - 1;}
01908 
01909         // retrieve the voxel values
01910         images = imagestack->get_data();
01911 
01912         // count the number of voxels within a sphere centered at icent,
01913         // with radius ri
01914         status = getnnz(volsize, ri, origin, &nrays, &nnz);
01915         // need to check status...
01916 
01917         // convert from cube to sphere
01918         sphere = new float[nnz];
01919         ptrs   = new int[nrays+1];
01920         cord   = new int[3*nrays];
01921         if (sphere == NULL || ptrs == NULL || cord == NULL) {
01922                 fprintf(stderr,"ChaoProjector::backproject3d, failed to allocate!\n");
01923                 exit(1);
01924         }
01925         for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01926         for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01927         for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01928 
01929         int nangles = 0;
01930         vector<float> anglelist;
01931         string angletype = "SPIDER";
01932         // Do we have a list of angles?
01933         if (params.has_key("anglelist")) {
01934                 anglelist = params["anglelist"];
01935                 nangles = anglelist.size() / 3;
01936         } else {
01937                 Transform* t3d = params["transform"];
01938                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01939                 // This part was modified by David Woolford -
01940                 // Before this the code worked only for SPIDER and EMAN angles,
01941                 // but the framework of the Transform3D allows for a generic implementation
01942                 // as specified here.
01943                 //  This was broken by david.  we need here a loop over all projections and put all angles on stack  PAP 06/28/09
01944                 Dict p = t3d->get_rotation("spider");
01945                 if(t3d) {delete t3d; t3d=0;}
01946 
01947                 float phi = p["phi"];
01948                 float theta = p["theta"];
01949                 float psi = p["psi"];
01950                 anglelist.push_back(phi);
01951                 anglelist.push_back(theta);
01952                 anglelist.push_back(psi);
01953                 nangles = 1;
01954         }
01955 
01956         // End David Woolford modifications
01957 
01958         if (nslices != nangles) {
01959                 LOGERR("the number of images does not match the number of angles");
01960                 return 0;
01961         }
01962 
01963         dm = new float[nangles*9];
01964         setdm(anglelist, angletype, dm);
01965 
01966         // return volume
01967         EMData *ret = new EMData();
01968         ret->set_size(nximg, nyimg, dim);
01969         ret->set_complex(false);
01970         ret->set_ri(true);
01971         ret->to_zero();
01972 
01973         cube = ret->get_data();
01974         // cb2sph should be replaced by something that touches only ptrs and cord
01975         status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01976         // check status
01977 
01978         for (j = 1; j <= nangles; j++) {
01979                 status = bckpj3(volsize, nrays, nnz, &dm(1,j), origin, ri,
01980                          ptrs   , cord , &images(1,1,j), sphere);
01981         // check status?
01982         }
01983 
01984         status = sph2cb(sphere, volsize, nrays, ri, nnz, ptrs, cord, cube);
01985         // check status?
01986 
01987         // deallocate all temporary work space
01988         EMDeleteArray(dm);
01989         EMDeleteArray(ptrs);
01990         EMDeleteArray(cord);
01991         EMDeleteArray(sphere);
01992 
01993         ret->update();
01994         return ret;
01995 }
01996 
01997 #undef images
01998 #undef cube
01999 #undef sphere
02000 #undef cord
02001 #undef ptrs
02002 #undef dm
02003 
02004 EMData *GaussFFTProjector::backproject3d(EMData * ) const
02005 {
02006     // no implementation yet
02007     EMData *ret = new EMData();
02008     return ret;
02009 }
02010 
02011 #define images(i,j,k) images[ (k)*nx*ny + ((j)-1)*nx + (i)-1 ]
02012 
02013 // EMData *PawelProjector::backproject3d(EMData * imagestack) const
02014 // {
02015 //
02016 //     float *images;
02017 //
02018 //     if (!imagestack) {
02019 //      return 0;
02020 //     }
02021 //     int ri;
02022 //     int nx      = imagestack->get_xsize();
02023 //     int ny      = imagestack->get_ysize();
02024 // //     int nslices = imagestack->get_zsize();
02025 //     int dim = Util::get_min(nx,ny);
02026 //     images  = imagestack->get_data();
02027 //
02028 //     Vec3i origin(0,0,0);
02029 //     // If a sensible origin isn't passed in, choose the middle of
02030 //     // the cube.
02031 //     if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
02032 //     else {origin[0] = nx/2;}
02033 //     if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
02034 //     else {origin[1] = ny/2;}
02035 //     if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
02036 //     else {origin[2] = dim/2;}
02037 //
02038 //     if (params.has_key("radius")) {ri = params["radius"];}
02039 //     else {ri = dim/2 - 1;}
02040 //
02041 //     // Determine the number of rows (x-lines) within the radius
02042 //     int nn = -1;
02043 //     prepcubes(nx, ny, dim, ri, origin, nn);
02044 //     // nn is now the number of rows-1 within the radius
02045 //     // so we can create and fill the ipcubes
02046 //     IPCube* ipcube = new IPCube[nn+1];
02047 //     prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02048 //
02049 //      int nangles = 0;
02050 //      vector<float> anglelist;
02051 //      // Do we have a list of angles?
02052 //      if (params.has_key("anglelist")) {
02053 //              anglelist = params["anglelist"];
02054 //              nangles = anglelist.size() / 3;
02055 //      } else {
02056 //              Transform3D* t3d = params["t3d"];
02057 //              if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
02058 //              // This part was modified by David Woolford -
02059 //              // Before this the code worked only for SPIDER and EMAN angles,
02060 //              // but the framework of the Transform3D allows for a generic implementation
02061 //              // as specified here.
02062 //              Dict p = t3d->get_rotation(Transform3D::SPIDER);
02063 //
02064 //              string angletype = "SPIDER";
02065 //              float phi = p["phi"];
02066 //              float theta = p["theta"];
02067 //              float psi = p["psi"];
02068 //              anglelist.push_back(phi);
02069 //              anglelist.push_back(theta);
02070 //              anglelist.push_back(psi);
02071 //              nangles = 1;
02072 //      }
02073 //
02074 //      // End David Woolford modifications
02075 //
02076 //     // initialize return object
02077 //     EMData* ret = new EMData();
02078 //     ret->set_size(nx, ny, dim);
02079 //     ret->to_zero();
02080 //
02081 //     // loop over sets of angles
02082 //     for (int ia = 0; ia < nangles; ia++) {
02083 //        int indx = 3*ia;
02084 //         Transform3D rotation(Transform3D::SPIDER, float(anglelist[indx]),
02085 //                             float(anglelist[indx+1]),
02086 //                             float(anglelist[indx+2]));
02087 //        float dm1 = rotation.at(0,0);
02088 //        float dm4 = rotation.at(1,0);
02089 //
02090 //        if (2*(ri+1)+1 > dim) {
02091 //           // Must check x and y boundaries
02092 //           LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02093 //           return 0;
02094 //        } else {
02095 //           // No need to check x and y boundaries
02096 //           for (int i = 0 ; i <= nn; i++) {
02097 //              int iox = (int)ipcube[i].loc[0]+origin[0];
02098 //              int ioy = (int)ipcube[i].loc[1]+origin[1];
02099 //              int ioz = (int)ipcube[i].loc[2]+origin[2];
02100 //
02101 //              Vec3f vb = rotation*ipcube[i].loc + origin;
02102 //              for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02103 //                 float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02104 //                 int   iqx = (int)floor(xbb);
02105 //
02106 //                 float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02107 //                 int   iqy = (int)floor(ybb);
02108 //
02109 //                 float dipx = xbb - iqx;
02110 //                 float dipy = ybb - iqy;
02111 //
02112 //                 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02113 //                     + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02114 //                     + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02115 //                     + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02116 //                     - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02117 //                 iox++;
02118 //              } // end for j
02119 //        } // end for i
02120 //        } // end if
02121 //     } // end for ia
02122 //
02123 //     ret->update();
02124 //     EMDeleteArray(ipcube);
02125 //     return ret;
02126 // }
02127 EMData *PawelProjector::backproject3d(EMData * imagestack) const
02128 {
02129 
02130         float *images;
02131 
02132         if (!imagestack) {
02133                 return 0;
02134         }
02135         int ri;
02136         int nx      = imagestack->get_xsize();
02137         int ny      = imagestack->get_ysize();
02138 //     int nslices = imagestack->get_zsize();
02139         int dim = Util::get_min(nx,ny);
02140         images  = imagestack->get_data();
02141 
02142         Vec3i origin(0,0,0);
02143     // If a sensible origin isn't passed in, choose the middle of
02144     // the cube.
02145         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
02146         else {origin[0] = nx/2;}
02147         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
02148         else {origin[1] = ny/2;}
02149         if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
02150         else {origin[2] = dim/2;}
02151 
02152         if (params.has_key("radius")) {ri = params["radius"];}
02153         else {ri = dim/2 - 1;}
02154 
02155     // Determine the number of rows (x-lines) within the radius
02156         int nn = -1;
02157         prepcubes(nx, ny, dim, ri, origin, nn);
02158     // nn is now the number of rows-1 within the radius
02159     // so we can create and fill the ipcubes
02160         IPCube* ipcube = new IPCube[nn+1];
02161         prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02162 
02163         int nangles = 0;
02164         vector<float> anglelist;
02165         // Do we have a list of angles?
02166         if (params.has_key("anglelist")) {
02167                 anglelist = params["anglelist"];
02168                 nangles = anglelist.size() / 3;
02169         } else {
02170                 Transform* t3d = params["transform"];
02171                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
02172                 // This part was modified by David Woolford -
02173                 // Before this the code worked only for SPIDER and EMAN angles,
02174                 // but the framework of the Transform3D allows for a generic implementation
02175                 // as specified here.
02176                 Dict p = t3d->get_rotation("spider");
02177                 if(t3d) {delete t3d; t3d=0;}
02178 
02179                 string angletype = "SPIDER";
02180                 float phi = p["phi"];
02181                 float theta = p["theta"];
02182                 float psi = p["psi"];
02183                 anglelist.push_back(phi);
02184                 anglelist.push_back(theta);
02185                 anglelist.push_back(psi);
02186                 nangles = 1;
02187         }
02188 
02189         // End David Woolford modifications
02190 
02191     // initialize return object
02192         EMData* ret = new EMData();
02193         ret->set_size(nx, ny, dim);
02194         ret->to_zero();
02195 
02196     // loop over sets of angles
02197         for (int ia = 0; ia < nangles; ia++) {
02198                 int indx = 3*ia;
02199                 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
02200                 Transform rotation(d);
02201                 float dm1 = rotation.at(0,0);
02202                 float dm4 = rotation.at(1,0);
02203 
02204                 if (2*(ri+1)+1 > dim) {
02205           // Must check x and y boundaries
02206                         LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02207                         return 0;
02208                 } else {
02209           // No need to check x and y boundaries
02210                         for (int i = 0 ; i <= nn; i++) {
02211                                 int iox = (int)ipcube[i].loc[0]+origin[0];
02212                                 int ioy = (int)ipcube[i].loc[1]+origin[1];
02213                                 int ioz = (int)ipcube[i].loc[2]+origin[2];
02214 
02215                                 Vec3f vb = rotation*ipcube[i].loc + origin;
02216                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02217                                         float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02218                                         int   iqx = (int)floor(xbb);
02219 
02220                                         float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02221                                         int   iqy = (int)floor(ybb);
02222 
02223                                         float dipx = xbb - iqx;
02224                                         float dipy = ybb - iqy;
02225 
02226                                         (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02227                                                         + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02228                                                         + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02229                                                                         + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02230                                                                                         - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02231                                         iox++;
02232                                 } // end for j
02233                         } // end for i
02234                 } // end if
02235         } // end for ia
02236 
02237         ret->update();
02238         EMDeleteArray(ipcube);
02239         return ret;
02240 }
02241 #undef images
02242 
02243 EMData *StandardProjector::backproject3d(EMData * ) const
02244 {
02245    // no implementation yet
02246    EMData *ret = new EMData();
02247    return ret;
02248 }
02249 
02250 EMData *MaxValProjector::backproject3d(EMData * ) const
02251 {
02252    // no implementation yet
02253    EMData *ret = new EMData();
02254    return ret;
02255 }
02256 
02257 
02258 EMData *FourierGriddingProjector::backproject3d(EMData * ) const
02259 {
02260    // no implementation yet
02261    EMData *ret = new EMData();
02262    return ret;
02263 }
02264 
02265 // End Chao's projector addition 4/25/06
02266 
02267 void EMAN::dump_projectors()
02268 {
02269         dump_factory < Projector > ();
02270 }
02271 
02272 map<string, vector<string> > EMAN::dump_projectors_list()
02273 {
02274         return dump_factory_list < Projector > ();
02275 }
02276 
02277 /* vim: set ts=4 noet nospell: */

Generated on Tue Jun 11 13:46:18 2013 for EMAN2 by  doxygen 1.3.9.1