EMAN2
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 EMData *PawelProjector::project3d(EMData * image) const
00577 {
00578         if (!image)  return 0;
00579         int ri;
00580         int nx = image->get_xsize();
00581         int ny = image->get_ysize();
00582         int nz = image->get_zsize();
00583         int dim = Util::get_max(nx,ny,nz);
00584         int dimc = dim/2;
00585         if (nz == 1) {
00586                 LOGERR("The PawelProjector needs a volume!");
00587                 return 0;
00588         }
00589 
00590         Vec3i origin(0,0,0);
00591 
00592         if (params.has_key("radius")) {ri = params["radius"];}
00593         else {ri = dim/2 - 1;}
00594 
00595         Transform* rotation = params["transform"];
00596         //int nangles = 0;
00597         //vector<float> anglelist;
00598         // Do we have a list of angles?
00599         /*
00600         if (params.has_key("anglelist")) {
00601                 anglelist = params["anglelist"];
00602                 nangles = anglelist.size() / 3;
00603         } else {*/
00604 
00605                 //if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
00606                 /*
00607                 Dict p = t3d->get_rotation("spider");
00608 
00609                 string angletype = "SPIDER";
00610                 float phi = p["phi"];
00611                 float theta = p["theta"];
00612                 float psi = p["psi"];
00613                 anglelist.push_back(phi);
00614                 anglelist.push_back(theta);
00615                 anglelist.push_back(psi);
00616                 */
00617                 int nangles = 1;
00618         //}
00619 
00620         //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;
00621         // loop over sets of angles
00622         //for (int ia = 0; ia < nangles; ia++) {
00623         EMData* ret = new EMData();
00624         int ia = 0;
00625                 //int indx = 3*ia;
00626                 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
00627                 //Transform rotation(d);
00628                 if (2*(ri+1)+1 > dim) {
00629                         // initialize return object
00630                         ret->set_size(dim, dim, nangles);
00631                         ret->to_zero();
00632                         origin[0] = dimc;
00633                         origin[1] = dimc;
00634                         origin[2] = dimc;
00635                         Vec3i loc(-dimc,0,0);
00636                         Vec3i vorg(nx/2,ny/2,nz/2);
00637                         // This code is for arbitrary dimensions, so must check x and y boundaries
00638                         for (int l = 0 ; l < dim; l++) {  // Z
00639                                 loc[2] = l - dimc;
00640                                 for (int k = 0 ; k < dim; k++) {  // Y
00641                                         loc[1] = k - dimc;
00642                                         Vec3f vb = loc*(*rotation) + vorg;
00643                                         for (int j = 0; j < dim; j++) {  //X
00644                                                 // check for pixels out-of-bounds
00645                                                 //cout<<j<<"  j  "<<k<<"  k  "<<"  iox  "<<vb[0]<<"  ioy  "<<vb[1]<<"  ioz "<<vb[2]<<endl;
00646                                                 int iox = int(vb[0]);
00647                                                 if ((iox >= 0) && (iox < nx-1)) {
00648                                                         int ioy = int(vb[1]);
00649                                                         if ((ioy >= 0) && (ioy < ny-1)) {
00650                                                                 int ioz = int(vb[2]);
00651                                                                 if ((ioz >= 0) && (ioz < nz-1)) {
00652                                                                         // real work for pixels in bounds
00653                                                                 //cout<<j<<"  j  "<<k<<"  k  "<<"  iox  "<<iox<<"  ioy  "<<ioy<<"  ioz "<<ioz<<endl;
00654                                                                 //cout<<"  TAKE"<<endl;
00655                                                                         float dx = vb[0] - iox;
00656                                                                         float dy = vb[1] - ioy;
00657                                                                         float dz = vb[2] - ioz;
00658                                                                         float a1 = (*image)(iox,ioy,ioz);
00659                                                                         float a2 = (*image)(iox+1,ioy,ioz) - a1;
00660                                                                         float a3 = (*image)(iox,ioy+1,ioz) - a1;
00661                                                                         float a4 = (*image)(iox,ioy,ioz+1) - a1;
00662                                                                         float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00663                                                                                                 + (*image)(iox+1,ioy+1,ioz);
00664                                                                         float a61 = -(*image)(iox,ioy,ioz+1)
00665                                                                                                 + (*image)(iox+1,ioy,ioz+1);
00666                                                                         float a6 = -a2 + a61;
00667                                                                         float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00668                                                                                                 + (*image)(iox,ioy+1,ioz+1);
00669                                                                         float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00670                                                                                                 + (*image)(iox+1,ioy+1,ioz+1);
00671                                                                         (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00672                                                                                                 + (a7 + a8*dx)*dy)
00673                                                                                                 + a3*dy + dx*(a2 + a5*dy);
00674                                                                 } //else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLAATED Z "<<endl; }
00675                                                         }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED Y "<<endl; }
00676                                                 }//else {cout<<"  iox  "<<iox<<"  ioy  "<<int(vb[1])<<"  ioz "<<int(vb[2])<<"  VIOLATED X "<<endl; }
00677                                                 vb += rotation->get_matrix3_row(0);
00678                                         }
00679                                 }
00680                         }
00681 
00682                 } else {
00683                         // If a sensible origin isn't passed in, choose the middle of
00684                         // the cube.
00685                         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
00686                         else {origin[0] = nx/2;}
00687                         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
00688                         else {origin[1] = ny/2;}
00689                         if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
00690                         else {origin[2] = nz/2;}
00691                         // Determine the number of rows (x-lines) within the radius
00692                         int nn = -1;
00693                         prepcubes(nx, ny, nz, ri, origin, nn);
00694                         // nn is now the number of rows-1 within the radius
00695                         // so we can create and fill the ipcubes
00696                         IPCube* ipcube = new IPCube[nn+1];
00697                         prepcubes(nx, ny, nz, ri, origin, nn, ipcube);
00698                         // initialize return object
00699                         ret->set_size(nx, ny, nangles);
00700                         ret->to_zero();
00701                         // No need to check x and y boundaries
00702                         for (int i = 0 ; i <= nn; i++) {
00703                                 int k = ipcube[i].loc[1] + origin[1];
00704                                 Vec3f vb = ipcube[i].loc*(*rotation) + origin;
00705                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
00706                                         int iox = int(vb[0]);
00707                                         int ioy = int(vb[1]);
00708                                         int ioz = int(vb[2]);
00709                                         float dx = vb[0] - iox;
00710                                         float dy = vb[1] - ioy;
00711                                         float dz = vb[2] - ioz;
00712                                         float a1 = (*image)(iox,ioy,ioz);
00713                                         float a2 = (*image)(iox+1,ioy,ioz) - a1;
00714                                         float a3 = (*image)(iox,ioy+1,ioz) - a1;
00715                                         float a4 = (*image)(iox,ioy,ioz+1) - a1;
00716                                         float a5 = -a2 -(*image)(iox,ioy+1,ioz)
00717                                                                 + (*image)(iox+1,ioy+1,ioz);
00718                                         float a61 = -(*image)(iox,ioy,ioz+1)
00719                                                                 + (*image)(iox+1,ioy,ioz+1);
00720                                         float a6 = -a2 + a61;
00721                                         float a7 = -a3 - (*image)(iox,ioy,ioz+1)
00722                                                                 + (*image)(iox,ioy+1,ioz+1);
00723                                         float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1)
00724                                                                 + (*image)(iox+1,ioy+1,ioz+1);
00725                                         (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx
00726                                                                 + (a7 + a8*dx)*dy)
00727                                                                 + a3*dy + dx*(a2 + a5*dy);
00728                                         vb += rotation->get_matrix3_row(0);
00729                                 }
00730                         }
00731                         EMDeleteArray(ipcube);
00732                 }
00733         //}
00734         ret->update();
00735         if(rotation) {delete rotation; rotation=0;}
00736         return ret;
00737 }
00738 
00739 EMData *StandardProjector::project3d(EMData * image) const
00740 {
00741         Transform* t3d = params["transform"];
00742         if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00743 //      Dict p = t3d->get_rotation();
00744         if ( image->get_ndim() == 3 )
00745         {
00746 
00747 #ifdef EMAN2_USING_CUDA
00748                 if(EMData::usecuda == 1) {
00749                         if(!image->isrodataongpu()) image->copy_to_cudaro();
00750                         //cout << "CUDA PROJ" << endl;
00751                         Transform* t3d = params["transform"];
00752                         if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00753                         float * m = new float[12];
00754                         t3d->copy_matrix_into_array(m);
00755                         image->bindcudaarrayA(true);
00756                         //EMData* e = new EMData(0,0,image->get_xsize(),image->get_ysize(),1);
00757                         EMData *e = new EMData();
00758                         e->set_size_cuda(image->get_xsize(), image->get_ysize(), 1);
00759                         e->rw_alloc();
00760                         standard_project(m,e->getcudarwdata(), e->get_xsize(), e->get_ysize());
00761                         image->unbindcudaarryA();
00762                         delete [] m;
00763                 
00764                         e->update();
00765                         e->set_attr("xform.projection",t3d);
00766                         e->set_attr("apix_x",(float)image->get_attr("apix_x"));
00767                         e->set_attr("apix_y",(float)image->get_attr("apix_y"));
00768                         e->set_attr("apix_z",(float)image->get_attr("apix_z"));
00769                         //e_>copy_from_device();
00770                         if(t3d) {delete t3d; t3d=0;}
00771                         return e;
00772                 }
00773 #endif
00774                 int nx = image->get_xsize();
00775                 int ny = image->get_ysize();
00776                 int nz = image->get_zsize();
00777 
00778 //              Transform3D r(Transform3D::EMAN, az, alt, phi);
00779                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
00780                 int xy = nx * ny;
00781 
00782                 EMData *proj = new EMData();
00783                 proj->set_size(nx, ny, 1);
00784 
00785                 Vec3i offset(nx/2,ny/2,nz/2);
00786 
00787                 float *sdata = image->get_data();
00788                 float *ddata = proj->get_data();
00789                 for (int k = -nz / 2; k < nz - nz / 2; k++) {
00790                         int l = 0;
00791                         for (int j = -ny / 2; j < ny - ny / 2; j++) {
00792                                 ddata[l]=0;
00793                                 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00794 
00795                                         Vec3f coord(i,j,k);
00796                                         Vec3f soln = r*coord;
00797                                         soln += offset;
00798 
00802 //                                      printf(" ");
00803 
00804                                         float x2 = soln[0];
00805                                         float y2 = soln[1];
00806                                         float z2 = soln[2];
00807 
00808                                         float x = (float)Util::fast_floor(x2);
00809                                         float y = (float)Util::fast_floor(y2);
00810                                         float z = (float)Util::fast_floor(z2);
00811 
00812                                         float t = x2 - x;
00813                                         float u = y2 - y;
00814                                         float v = z2 - z;
00815 
00816                                         size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
00817 // 
00818                                         if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
00819                                         if      (x2 > (nx-1) || y2  > (ny-1) || z2 > (nz-1) ) continue;
00820 
00821                                         if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
00822                                                 ddata[l] +=
00823                                                                 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
00824                                                                 sdata[ii + nx + 1], sdata[ii + xy],     sdata[ii + xy + 1], sdata[ii + xy + nx],
00825                                                                 sdata[ii + xy + nx + 1], t, u, v);
00826                                         }
00827                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
00828                                                 ddata[l] += sdata[ii];
00829                                         }
00830                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
00831                                                 ddata[l] +=     Util::linear_interpolate(sdata[ii], sdata[ii + xy],v);
00832                                         }
00833                                         else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
00834                                                 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + nx],u);
00835                                         }
00836                                         else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
00837                                                 ddata[l] += Util::linear_interpolate(sdata[ii], sdata[ii + 1],t);
00838                                         }
00839                                         else if ( x2 == (nx - 1) ) {
00840                                                 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v);
00841                                         }
00842                                         else if ( y2 == (ny - 1) ) {
00843                                                 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v);
00844                                         }
00845                                         else if ( z2 == (nz - 1) ) {
00846                                                 ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u);
00847                                         }
00848                                 }
00849                         }
00850                 }
00851                 proj->update();
00852                 proj->set_attr("xform.projection",t3d);
00853                 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
00854                 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
00855                 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
00856                 
00857                 if(t3d) {delete t3d; t3d=0;}
00858                 return proj;
00859         }
00860         else if ( image->get_ndim() == 2 ) {
00861 
00862                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
00863 
00864                 int nx = image->get_xsize();
00865                 int ny = image->get_ysize();
00866 
00867                 EMData *proj = new EMData();
00868                 proj->set_size(nx, 1, 1);
00869                 proj->to_zero();
00870 
00871                 float *sdata = image->get_data();
00872                 float *ddata = proj->get_data();
00873 
00874                 Vec2f offset(nx/2,ny/2);
00875                 for (int j = -ny / 2; j < ny - ny / 2; j++) { // j represents a column of pixels in the direction of the angle
00876                         int l = 0;
00877                         for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00878 
00879                                 Vec2f coord(i,j);
00880                                 Vec2f soln = r*coord;
00881                                 soln += offset;
00882 
00883                                 float x2 = soln[0];
00884                                 float y2 = soln[1];
00885 
00886                                 float x = (float)Util::fast_floor(x2);
00887                                 float y = (float)Util::fast_floor(y2);
00888 
00889                                 int ii = (int) (x + y * nx);
00890                                 float u = x2 - x;
00891                                 float v = y2 - y;
00892 
00893                                 if (x2 < 0 || y2 < 0 ) continue;
00894                                 if      (x2 > (nx-1) || y2  > (ny-1) ) continue;
00895 
00896                                 if (  x2 < (nx - 1) && y2 < (ny - 1) ) {
00897                                         ddata[l] += Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v);
00898                                 }
00899                                 else if (x2 == (nx-1) && y2 == (ny-1) ) {
00900                                         ddata[l] += sdata[ii];
00901                                 }
00902                                 else if (x2 == (nx-1) ) {
00903                                         ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + nx], v);
00904                                 }
00905                                 else if (y2 == (ny-1) ) {
00906                                         ddata[l] += Util::linear_interpolate(sdata[ii],sdata[ii + 1], u);
00907                                 }
00908                         }
00909                 }
00910                 proj->set_attr("xform.projection",t3d);
00911                 proj->update();
00912                 if(t3d) {delete t3d; t3d=0;}
00913                 return proj;
00914         }
00915         else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
00916 }
00917 
00918 EMData *MaxValProjector::project3d(EMData * image) const
00919 {
00920         Transform* t3d = params["transform"];
00921         if ( t3d == NULL ) throw NullPointerException("The transform object containing the angles(required for projection), was not specified");
00922 //      Dict p = t3d->get_rotation();
00923         if ( image->get_ndim() == 3 )
00924         {
00925 
00926                 int nx = image->get_xsize();
00927                 int ny = image->get_ysize();
00928                 int nz = image->get_zsize();
00929 
00930 //              Transform3D r(Transform3D::EMAN, az, alt, phi);
00931                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
00932                 int xy = nx * ny;
00933 
00934                 EMData *proj = new EMData();
00935                 proj->set_size(nx, ny, 1);
00936 
00937                 Vec3i offset(nx/2,ny/2,nz/2);
00938 
00939                 float *sdata = image->get_data();
00940                 float *ddata = proj->get_data();
00941                 for (int k = -nz / 2; k < nz - nz / 2; k++) {
00942                         int l = 0;
00943                         for (int j = -ny / 2; j < ny - ny / 2; j++) {
00944                                 ddata[l]=0;
00945                                 for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
00946 
00947                                         Vec3f coord(i,j,k);
00948                                         Vec3f soln = r*coord;
00949                                         soln += offset;
00950 
00954 //                                      printf(" ");
00955 
00956                                         float x2 = soln[0];
00957                                         float y2 = soln[1];
00958                                         float z2 = soln[2];
00959 
00960                                         float x = (float)Util::fast_floor(x2);
00961                                         float y = (float)Util::fast_floor(y2);
00962                                         float z = (float)Util::fast_floor(z2);
00963 
00964                                         float t = x2 - x;
00965                                         float u = y2 - y;
00966                                         float v = z2 - z;
00967 
00968                                         size_t ii = (size_t) ((size_t)x + (size_t)y * nx + (size_t)z * xy);
00969 // 
00970                                         if (x2 < 0 || y2 < 0 || z2 < 0 ) continue;
00971                                         if      (x2 > (nx-1) || y2  > (ny-1) || z2 > (nz-1) ) continue;
00972 
00973                                         if (x2 < (nx - 1) && y2 < (ny - 1) && z2 < (nz - 1)) {
00974                                                 ddata[l] = Util::get_max(ddata[l],
00975                                                                 Util::trilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],
00976                                                                 sdata[ii + nx + 1], sdata[ii + xy],     sdata[ii + xy + 1], sdata[ii + xy + nx],
00977                                                                 sdata[ii + xy + nx + 1], t, u, v));
00978                                         }
00979                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) && z2 == (nz - 1) ) {
00980                                                 ddata[l] = Util::get_max(ddata[l],sdata[ii]);
00981                                         }
00982                                         else if ( x2 == (nx - 1) && y2 == (ny - 1) ) {
00983                                                 ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + xy],v));
00984                                         }
00985                                         else if ( x2 == (nx - 1) && z2 == (nz - 1) ) {
00986                                                 ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + nx],u));
00987                                         }
00988                                         else if ( y2 == (ny - 1) && z2 == (nz - 1) ) {
00989                                                 ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii], sdata[ii + 1],t));
00990                                         }
00991                                         else if ( x2 == (nx - 1) ) {
00992                                                 ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + nx], sdata[ii + xy], sdata[ii + xy + nx],u,v));
00993                                         }
00994                                         else if ( y2 == (ny - 1) ) {
00995                                                 ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + xy], sdata[ii + xy + 1],t,v));
00996                                         }
00997                                         else if ( z2 == (nz - 1) ) {
00998                                                 ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx], sdata[ii + nx + 1],t,u));
00999                                         }
01000                                 }
01001                         }
01002                 }
01003                 proj->update();
01004                 proj->set_attr("xform.projection",t3d);
01005                 proj->set_attr("apix_x",(float)image->get_attr("apix_x"));
01006                 proj->set_attr("apix_y",(float)image->get_attr("apix_y"));
01007                 proj->set_attr("apix_z",(float)image->get_attr("apix_z"));
01008                 
01009                 if(t3d) {delete t3d; t3d=0;}
01010                 return proj;
01011         }
01012         else if ( image->get_ndim() == 2 ) {
01013 
01014                 Transform r = t3d->inverse(); // The inverse is taken here because we are rotating the coordinate system, not the image
01015 
01016                 int nx = image->get_xsize();
01017                 int ny = image->get_ysize();
01018 
01019                 EMData *proj = new EMData();
01020                 proj->set_size(nx, 1, 1);
01021                 proj->to_zero();
01022 
01023                 float *sdata = image->get_data();
01024                 float *ddata = proj->get_data();
01025 
01026                 Vec2f offset(nx/2,ny/2);
01027                 for (int j = -ny / 2; j < ny - ny / 2; j++) { // j represents a column of pixels in the direction of the angle
01028                         int l = 0;
01029                         for (int i = -nx / 2; i < nx - nx / 2; i++,l++) {
01030 
01031                                 Vec2f coord(i,j);
01032                                 Vec2f soln = r*coord;
01033                                 soln += offset;
01034 
01035                                 float x2 = soln[0];
01036                                 float y2 = soln[1];
01037 
01038                                 float x = (float)Util::fast_floor(x2);
01039                                 float y = (float)Util::fast_floor(y2);
01040 
01041                                 int ii = (int) (x + y * nx);
01042                                 float u = x2 - x;
01043                                 float v = y2 - y;
01044 
01045                                 if (x2 < 0 || y2 < 0 ) continue;
01046                                 if      (x2 > (nx-1) || y2  > (ny-1) ) continue;
01047 
01048                                 if (  x2 < (nx - 1) && y2 < (ny - 1) ) {
01049                                         ddata[l] =      Util::get_max(ddata[l],Util::bilinear_interpolate(sdata[ii], sdata[ii + 1], sdata[ii + nx],sdata[ii + nx + 1], u, v));
01050                                 }
01051                                 else if (x2 == (nx-1) && y2 == (ny-1) ) {
01052                                         ddata[l] =      Util::get_max(ddata[l],sdata[ii]);
01053                                 }
01054                                 else if (x2 == (nx-1) ) {
01055                                         ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii],sdata[ii + nx], v));
01056                                 }
01057                                 else if (y2 == (ny-1) ) {
01058                                         ddata[l] =      Util::get_max(ddata[l],Util::linear_interpolate(sdata[ii],sdata[ii + 1], u));
01059                                 }
01060                         }
01061                 }
01062                 proj->set_attr("xform.projection",t3d);
01063                 proj->update();
01064                 if(t3d) {delete t3d; t3d=0;}
01065                 return proj;
01066         }
01067         else throw ImageDimensionException("Standard projection works only for 2D and 3D images");
01068 }
01069 
01070 
01071 // EMData *FourierGriddingProjector::project3d(EMData * image) const
01072 // {
01073 //      if (!image) {
01074 //              return 0;
01075 //      }
01076 //      if (3 != image->get_ndim())
01077 //              throw ImageDimensionException(
01078 //                              "FourierGriddingProjector needs a 3-D volume");
01079 //      if (image->is_complex())
01080 //              throw ImageFormatException(
01081 //                              "FourierGriddingProjector requires a real volume");
01082 //      const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
01083 //      const int nx = image->get_xsize();
01084 //      const int ny = image->get_ysize();
01085 //      const int nz = image->get_zsize();
01086 //      if (nx != ny || nx != nz)
01087 //              throw ImageDimensionException(
01088 //                              "FourierGriddingProjector requires nx==ny==nz");
01089 //      const int m = Util::get_min(nx,ny,nz);
01090 //      const int n = m*npad;
01091 //
01092 //      int K = params["kb_K"];
01093 //      if ( K == 0 ) K = 6;
01094 //      float alpha = params["kb_alpha"];
01095 //      if ( alpha == 0 ) alpha = 1.25;
01096 //      Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
01097 //
01098 //      // divide out gridding weights
01099 //      EMData* tmpImage = image->copy();
01100 //      tmpImage->divkbsinh(kb);
01101 //      // pad and center volume, then FFT and multiply by (-1)**(i+j+k)
01102 //      //EMData* imgft = tmpImage->pad_fft(npad);
01103 //      //imgft->center_padded();
01104 //      EMData* imgft = tmpImage->norm_pad(false, npad);
01105 //      imgft->do_fft_inplace();
01106 //      imgft->center_origin_fft();
01107 //      imgft->fft_shuffle();
01108 //      delete tmpImage;
01109 //
01110 //      // Do we have a list of angles?
01111 //      int nangles = 0;
01112 //      vector<float> anglelist;
01113 //      // Do we have a list of angles?
01114 //      if (params.has_key("anglelist")) {
01115 //              anglelist = params["anglelist"];
01116 //              nangles = anglelist.size() / 3;
01117 //      } else {
01118 //              // This part was modified by David Woolford -
01119 //              // Before this the code worked only for SPIDER and EMAN angles,
01120 //              // but the framework of the Transform3D allows for a generic implementation
01121 //              // as specified here.
01122 //              Transform3D* t3d = params["t3d"];
01123 //              if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
01124 //              Dict p = t3d->get_rotation(Transform3D::SPIDER);
01125 //
01126 //              string angletype = "SPIDER";
01127 //              float phi = p["phi"];
01128 //              float theta = p["theta"];
01129 //              float psi = p["psi"];
01130 //              anglelist.push_back(phi);
01131 //              anglelist.push_back(theta);
01132 //              anglelist.push_back(psi);
01133 //              nangles = 1;
01134 //      }
01135 //
01136 //      // End David Woolford modifications
01137 //
01138 //      // initialize return object
01139 //      EMData* ret = new EMData();
01140 //      ret->set_size(nx, ny, nangles);
01141 //      ret->to_zero();
01142 //      // loop over sets of angles
01143 //      for (int ia = 0; ia < nangles; ia++) {
01144 //              int indx = 3*ia;
01145 //              Transform3D tf(Transform3D::SPIDER, anglelist[indx],anglelist[indx+1],anglelist[indx+2]);
01146 //              EMData* proj = imgft->extractplane(tf, kb);
01147 //              if (proj->is_shuffled()) proj->fft_shuffle();
01148 //              proj->center_origin_fft();
01149 //              proj->do_ift_inplace();
01150 //              EMData* winproj = proj->window_center(m);
01151 //              delete proj;
01152 //              for (int iy=0; iy < ny; iy++)
01153 //                      for (int ix=0; ix < nx; ix++)
01154 //                              (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
01155 //              delete winproj;
01156 //      }
01157 //      delete imgft;
01158 //      ret->update();
01159 //
01160 //      return ret;
01161 // }
01162 
01163 
01164 EMData *FourierGriddingProjector::project3d(EMData * image) const
01165 {
01166         if (!image) {
01167                 return 0;
01168         }
01169         if (3 != image->get_ndim())
01170                 throw ImageDimensionException(
01171                                                                           "FourierGriddingProjector needs a 3-D volume");
01172         if (image->is_complex())
01173                 throw ImageFormatException(
01174                                                                    "FourierGriddingProjector requires a real volume");
01175         const int npad = params.has_key("npad") ? int(params["npad"]) : 2;
01176         const int nx = image->get_xsize();
01177         const int ny = image->get_ysize();
01178         const int nz = image->get_zsize();
01179         if (nx != ny || nx != nz)
01180                 throw ImageDimensionException(
01181                                                                           "FourierGriddingProjector requires nx==ny==nz");
01182         const int m = Util::get_min(nx,ny,nz);
01183         const int n = m*npad;
01184 
01185         int K = params["kb_K"];
01186         if ( K == 0 ) K = 6;
01187         float alpha = params["kb_alpha"];
01188         if ( alpha == 0 ) alpha = 1.25;
01189         Util::KaiserBessel kb(alpha, K, (float)(m/2), K/(2.0f*n), n);
01190 
01191         // divide out gridding weights
01192         EMData* tmpImage = image->copy();
01193         tmpImage->divkbsinh(kb);
01194         // pad and center volume, then FFT and multiply by (-1)**(i+j+k)
01195         //EMData* imgft = tmpImage->pad_fft(npad);
01196         //imgft->center_padded();
01197         EMData* imgft = tmpImage->norm_pad(false, npad);
01198         imgft->do_fft_inplace();
01199         imgft->center_origin_fft();
01200         imgft->fft_shuffle();
01201         delete tmpImage;
01202 
01203         // Do we have a list of angles?
01204         int nangles = 0;
01205         vector<float> anglelist;
01206         // Do we have a list of angles?
01207         if (params.has_key("anglelist")) {
01208                 anglelist = params["anglelist"];
01209                 nangles = anglelist.size() / 3;
01210         } else {
01211                 // This part was modified by David Woolford -
01212                 // Before this the code worked only for SPIDER and EMAN angles,
01213                 // but the framework of the Transform3D allows for a generic implementation
01214                 // as specified here.
01215                 Transform* t3d = params["transform"];
01216                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01217                 Dict p = t3d->get_rotation("spider");
01218 
01219                 string angletype = "SPIDER";
01220                 float phi = p["phi"];
01221                 float theta = p["theta"];
01222                 float psi = p["psi"];
01223                 anglelist.push_back(phi);
01224                 anglelist.push_back(theta);
01225                 anglelist.push_back(psi);
01226                 nangles = 1;
01227                 if(t3d) {delete t3d; t3d=0;}
01228         }
01229 
01230         // End David Woolford modifications
01231 
01232         // initialize return object
01233         EMData* ret = new EMData();
01234         ret->set_size(nx, ny, nangles);
01235         ret->to_zero();
01236         // loop over sets of angles
01237         for (int ia = 0; ia < nangles; ia++) {
01238                 int indx = 3*ia;
01239                 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
01240                 Transform tf(d);
01241                 EMData* proj = imgft->extract_plane(tf, kb);
01242                 if (proj->is_shuffled()) proj->fft_shuffle();
01243                 proj->center_origin_fft();
01244                 proj->do_ift_inplace();
01245                 EMData* winproj = proj->window_center(m);
01246                 delete proj;
01247                 for (int iy=0; iy < ny; iy++)
01248                         for (int ix=0; ix < nx; ix++)
01249                                 (*ret)(ix,iy,ia) = (*winproj)(ix,iy);
01250                 delete winproj;
01251         }
01252         delete imgft;
01253 
01254         if (!params.has_key("anglelist")) {
01255                 Transform* t3d = params["transform"];
01256                 ret->set_attr("xform.projection",t3d);
01257                 if(t3d) {delete t3d; t3d=0;}
01258         }
01259         ret->update();
01260         return ret;
01261 }
01262 
01263 // BEGIN Chao projectors and backprojector addition (04/25/06)
01264 int ChaoProjector::getnnz(Vec3i volsize, int ri, Vec3i origin, int *nrays, int *nnz) const
01265 /*
01266    purpose: count the number of voxels within a sphere centered
01267             at origin and with a radius ri.
01268 
01269      input:
01270      volsize contains the size information (nx,ny,nz) about the volume
01271      ri      radius of the object embedded in the cube.
01272      origin  coordinates for the center of the volume
01273 
01274      output:
01275      nnz    total number of voxels within the sphere (of radius ri)
01276      nrays  number of rays in z-direction.
01277 */
01278 {
01279         int  ix, iy, iz, rs, r2, xs, ys, zs, xx, yy, zz;
01280         int  ftm=0, status = 0;
01281 
01282         r2    = ri*ri;
01283         *nnz  = 0;
01284         *nrays = 0;
01285         int nx = (int)volsize[0];
01286         int ny = (int)volsize[1];
01287         int nz = (int)volsize[2];
01288 
01289         int xcent = (int)origin[0];
01290         int ycent = (int)origin[1];
01291         int zcent = (int)origin[2];
01292 
01293         // need to add some error checking
01294         for (ix = 1; ix <=nx; ix++) {
01295             xs  = ix-xcent;
01296             xx  = xs*xs;
01297             for (iy = 1; iy <= ny; iy++) {
01298                 ys = iy-ycent;
01299                 yy = ys*ys;
01300                 ftm = 1;
01301                 for (iz = 1; iz <= nz; iz++) {
01302                     zs = iz-zcent;
01303                     zz = zs*zs;
01304                     rs = xx + yy + zz;
01305                     if (rs <= r2) {
01306                         (*nnz)++;
01307                         if (ftm) {
01308                            (*nrays)++;
01309                            ftm = 0;
01310                         }
01311                     }
01312                 }
01313             } // end for iy
01314         } // end for ix
01315         return status;
01316 }
01317 
01318 #define cube(i,j,k) cube[ ((k-1)*ny + j-1)*nx + i-1 ]
01319 #define sphere(i)   sphere[(i)-1]
01320 #define cord(i,j)   cord[((j)-1)*3 + (i) -1]
01321 #define ptrs(i)     ptrs[(i)-1]
01322 #define dm(i)       dm[(i)-1]
01323 
01324 int ChaoProjector:: cb2sph(float *cube, Vec3i volsize, int    ri, Vec3i origin,
01325                            int    nnz0, int     *ptrs, int *cord, float *sphere) const
01326 {
01327     int    xs, ys, zs, xx, yy, zz, rs, r2;
01328     int    ix, iy, iz, jnz, nnz, nrays;
01329     int    ftm = 0, status = 0;
01330 
01331     int xcent = (int)origin[0];
01332     int ycent = (int)origin[1];
01333     int zcent = (int)origin[2];
01334 
01335     int nx = (int)volsize[0];
01336     int ny = (int)volsize[1];
01337     int nz = (int)volsize[2];
01338 
01339     r2      = ri*ri;
01340     nnz     = 0;
01341     nrays    = 0;
01342     ptrs(1) = 1;
01343 
01344     for (ix = 1; ix <= nx; ix++) {
01345        xs  = ix-xcent;
01346        xx  = xs*xs;
01347        for ( iy = 1; iy <= ny; iy++ ) {
01348            ys = iy-ycent;
01349            yy = ys*ys;
01350            jnz = 0;
01351 
01352            ftm = 1;
01353            // not the most efficient implementation
01354            for (iz = 1; iz <= nz; iz++) {
01355                zs = iz-zcent;
01356                zz = zs*zs;
01357                rs = xx + yy + zz;
01358                if (rs <= r2) {
01359                   jnz++;
01360                   nnz++;
01361                   sphere(nnz) = cube(iz, iy, ix);
01362 
01363                   //  record the coordinates of the first nonzero ===
01364                   if (ftm) {
01365                      nrays++;
01366                      cord(1,nrays) = iz;
01367                      cord(2,nrays) = iy;
01368                      cord(3,nrays) = ix;
01369                      ftm = 0;
01370                   }
01371                }
01372             } // end for (iz..)
01373             if (jnz > 0) {
01374                 ptrs(nrays+1) = ptrs(nrays) + jnz;
01375             }  // endif (jnz)
01376        } // end for iy
01377     } // end for ix
01378     if (nnz != nnz0) status = -1;
01379     return status;
01380 }
01381 
01382 // decompress sphere into a cube
01383 int ChaoProjector::sph2cb(float *sphere, Vec3i volsize, int  nrays, int    ri,
01384                           int      nnz0, int     *ptrs, int  *cord, float *cube) const
01385 {
01386     int       status=0;
01387     int       r2, i, j, ix, iy, iz,  nnz;
01388 
01389     int nx = (int)volsize[0];
01390     int ny = (int)volsize[1];
01391     // int nz = (int)volsize[2];
01392 
01393     r2      = ri*ri;
01394     nnz     = 0;
01395     ptrs(1) = 1;
01396 
01397     // no need to initialize
01398     // for (i = 0; i<nx*ny*nz; i++) cube[i]=0.0;
01399 
01400     nnz = 0;
01401     for (j = 1; j <= nrays; j++) {
01402        iz = cord(1,j);
01403        iy = cord(2,j);
01404        ix = cord(3,j);
01405        for (i = ptrs(j); i<=ptrs(j+1)-1; i++, iz++) {
01406            nnz++;
01407            cube(iz,iy,ix) = sphere(nnz);
01408        }
01409     }
01410     if (nnz != nnz0) status = -1;
01411     return status;
01412 }
01413 
01414 #define x(i)        x[(i)-1]
01415 #define y(i,j)      y[(j-1)*nx + i - 1]
01416 
01417 // project from 3D to 2D (single image)
01418 int ChaoProjector::fwdpj3(Vec3i volsize, int nrays, int      , float *dm,
01419                           Vec3i  origin, int    ri, int *ptrs, int *cord,
01420                           float      *x, float  *y) const
01421 {
01422     /*
01423         purpose:  y <--- proj(x)
01424         input  :  volsize  the size (nx,ny,nz) of the volume
01425                   nrays    number of rays within the compact spherical
01426                            representation
01427                   nnz      number of voxels within the sphere
01428                   dm       an array of size 9 storing transformation
01429                            associated with the projection direction
01430                   origin   coordinates of the center of the volume
01431                   ri       radius of the sphere
01432                   ptrs     the beginning address of each ray
01433                   cord     the coordinates of the first point in each ray
01434                   x        3d input volume
01435                   y        2d output image
01436     */
01437 
01438     int    iqx, iqy, i, j, xc, yc, zc;
01439     float  ct, dipx, dipy, dipx1m, dipy1m, xb, yb, dm1, dm4;
01440     int    status = 0;
01441 
01442     int xcent = origin[0];
01443     int ycent = origin[1];
01444     int zcent = origin[2];
01445 
01446     int nx = volsize[0];
01447 
01448     dm1 = dm(1);
01449     dm4 = dm(4);
01450 
01451     if ( nx > 2*ri ) {
01452         for (i = 1; i <= nrays; i++) {
01453             zc = cord(1,i)-zcent;
01454             yc = cord(2,i)-ycent;
01455             xc = cord(3,i)-xcent;
01456 
01457             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01458             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01459 
01460             for (j = ptrs(i); j< ptrs(i+1); j++) {
01461                iqx = ifix(xb);
01462                iqy = ifix(yb);
01463 
01464                ct   = x(j);
01465                dipx =  xb - (float)(iqx);
01466                dipy = (yb - (float)(iqy)) * ct;
01467 
01468                dipy1m = ct - dipy;
01469                dipx1m = 1.0f - dipx;
01470 
01471                y(iqx  ,iqy)   = y(iqx  ,iqy)   + dipx1m*dipy1m;
01472                y(iqx+1,iqy)   = y(iqx+1,iqy)   + dipx*dipy1m;
01473                y(iqx+1,iqy+1) = y(iqx+1,iqy+1) + dipx*dipy;
01474                y(iqx  ,iqy+1) = y(iqx  ,iqy+1) + dipx1m*dipy;
01475 
01476                xb += dm1;
01477                yb += dm4;
01478            }
01479         }
01480     }
01481     else {
01482         fprintf(stderr, " nx must be greater than 2*ri\n");
01483         exit(1);
01484     }
01485     return status;
01486 }
01487 #undef x
01488 #undef y
01489 
01490 #define y(i)        y[(i)-1]
01491 #define x(i,j)      x[((j)-1)*nx + (i) - 1]
01492 
01493 // backproject from 2D to 3D for a single image
01494 int ChaoProjector::bckpj3(Vec3i volsize, int nrays, int      , float *dm,
01495                           Vec3i  origin, int    ri, int *ptrs, int *cord,
01496                           float      *x, float *y) const
01497 {
01498     int       i, j, iqx,iqy, xc, yc, zc;
01499     float     xb, yb, dx, dy, dx1m, dy1m, dxdy;
01500     int       status = 0;
01501 
01502     int xcent = origin[0];
01503     int ycent = origin[1];
01504     int zcent = origin[2];
01505 
01506     int nx = volsize[0];
01507 
01508     if ( nx > 2*ri) {
01509         for (i = 1; i <= nrays; i++) {
01510             zc = cord(1,i) - zcent;
01511             yc = cord(2,i) - ycent;
01512             xc = cord(3,i) - xcent;
01513 
01514             xb = zc*dm(1)+yc*dm(2)+xc*dm(3) + xcent;
01515             yb = zc*dm(4)+yc*dm(5)+xc*dm(6) + ycent;
01516 
01517             for (j = ptrs(i); j <ptrs(i+1); j++) {
01518                 iqx = ifix((float)(xb));
01519                 iqy = ifix((float)(yb));
01520 
01521                 dx = xb - (float)(iqx);
01522                 dy = yb - (float)(iqy);
01523                 dx1m = 1.0f - dx;
01524                 dy1m = 1.0f - dy;
01525                 dxdy = dx*dy;
01526 /*
01527 c               y(j) = y(j) + dx1m*dy1m*x(iqx  , iqy)
01528 c     &                     + dx1m*dy  *x(iqx  , iqy+1)
01529 c     &                     + dx  *dy1m*x(iqx+1, iqy)
01530 c     &                     + dx  *dy  *x(iqx+1, iqy+1)
01531 c
01532 c              --- faster version of the above commented out
01533 c                  code (derived by summing the following table
01534 c                  of coefficients along  the colunms) ---
01535 c
01536 c                        1         dx        dy      dxdy
01537 c                     ------   --------  --------  -------
01538 c                      x(i,j)   -x(i,j)   -x(i,j)    x(i,j)
01539 c                                        x(i,j+1) -x(i,j+1)
01540 c                              x(i+1,j)           -x(i+1,j)
01541 c                                                x(i+1,j+1)
01542 c
01543 */
01544                y(j) += x(iqx,iqy)
01545                     +  dx*(-x(iqx,iqy)+x(iqx+1,iqy))
01546                     +  dy*(-x(iqx,iqy)+x(iqx,iqy+1))
01547                     +  dxdy*( x(iqx,iqy) - x(iqx,iqy+1)
01548                              -x(iqx+1,iqy) + x(iqx+1,iqy+1) );
01549 
01550                xb += dm(1);
01551                yb += dm(4);
01552             } // end for j
01553         } // end for i
01554      }
01555     else {
01556         fprintf(stderr, "bckpj3: nx must be greater than 2*ri\n");
01557     }
01558 
01559     return status;
01560 }
01561 
01562 #undef x
01563 #undef y
01564 #undef dm
01565 
01566 // funny F90 style strange rounding function
01567 int ChaoProjector::ifix(float a) const
01568 {
01569     int ia;
01570 
01571     if (a>=0) {
01572        ia = (int)floor(a);
01573     }
01574     else {
01575        ia = (int)ceil(a);
01576     }
01577     return ia;
01578 }
01579 
01580 #define dm(i,j)          dm[((j)-1)*9 + (i) -1]
01581 #define anglelist(i,j)   anglelist[((j)-1)*3 + (i) - 1]
01582 
01583 // SPIDER stype transformation
01584 void ChaoProjector::setdm(vector<float> anglelist, string const , float *dm) const
01585 { // convert Euler angles to transformations, dm is an 9 by nangles array
01586 
01587         float  psi, theta, phi;
01588         double cthe, sthe, cpsi, spsi, cphi, sphi;
01589         int    j;
01590 
01591         int nangles = anglelist.size() / 3;
01592 
01593         // now convert all angles
01594         for (j = 1; j <= nangles; j++) {
01595                 phi   = static_cast<float>(anglelist(1,j)*dgr_to_rad);
01596                 theta = static_cast<float>(anglelist(2,j)*dgr_to_rad);
01597                 psi   = static_cast<float>(anglelist(3,j)*dgr_to_rad);
01598 
01599                 //              cout << phi << " " << theta << " " << psi << endl;
01600                 cthe  = cos(theta);
01601                 sthe  = sin(theta);
01602                 cpsi  = cos(psi);
01603                 spsi  = sin(psi);
01604                 cphi  = cos(phi);
01605                 sphi  = sin(phi);
01606 
01607                 dm(1,j)=static_cast<float>(cphi*cthe*cpsi-sphi*spsi);
01608                 dm(2,j)=static_cast<float>(sphi*cthe*cpsi+cphi*spsi);
01609                 dm(3,j)=static_cast<float>(-sthe*cpsi);
01610                 dm(4,j)=static_cast<float>(-cphi*cthe*spsi-sphi*cpsi);
01611                 dm(5,j)=static_cast<float>(-sphi*cthe*spsi+cphi*cpsi);
01612                 dm(6,j)=static_cast<float>(sthe*spsi);
01613                 dm(7,j)=static_cast<float>(sthe*cphi);
01614                 dm(8,j)=static_cast<float>(sthe*sphi);
01615                 dm(9,j)=static_cast<float>(cthe);
01616         }
01617 }
01618 #undef anglelist
01619 
01620 #define images(i,j,k) images[ ((k-1)*nyvol + j-1)*nxvol + i-1 ]
01621 
01622 EMData *ChaoProjector::project3d(EMData * vol) const
01623 {
01624 
01625         int nrays, nnz, status, j;
01626         float *dm;
01627         int   *ptrs, *cord;
01628         float *sphere, *images;
01629 
01630         int nxvol = vol->get_xsize();
01631         int nyvol = vol->get_ysize();
01632         int nzvol = vol->get_zsize();
01633         Vec3i volsize(nxvol,nyvol,nzvol);
01634 
01635         int dim = Util::get_min(nxvol,nyvol,nzvol);
01636         if (nzvol == 1) {
01637                 LOGERR("The ChaoProjector needs a volume!");
01638                 return 0;
01639         }
01640         Vec3i origin(0,0,0);
01641         // If a sensible origin isn't passed in, choose the middle of
01642         // the cube.
01643         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01644         else {origin[0] = nxvol/2+1;}
01645         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01646         else {origin[1] = nyvol/2+1;}
01647         if (params.has_key("origin_z")) {origin[2] = params["origin_z"];}
01648         else {origin[2] = nzvol/2+1;}
01649 
01650         int ri;
01651         if (params.has_key("radius")) {ri = params["radius"];}
01652         else {ri = dim/2 - 1;}
01653 
01654         // retrieve the voxel values
01655         float *cube = vol->get_data();
01656 
01657         // count the number of voxels within a sphere centered at icent,
01658         // with radius ri
01659         status = getnnz(volsize, ri, origin, &nrays, &nnz);
01660         // need to check status...
01661 
01662         // convert from cube to sphere
01663         sphere = new float[nnz];
01664         ptrs   = new int[nrays+1];
01665         cord   = new int[3*nrays];
01666         if (sphere == NULL || ptrs == NULL || cord == NULL) {
01667                 fprintf(stderr,"ChaoProjector::project3d, failed to allocate!\n");
01668                 exit(1);
01669         }
01670         for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01671         for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01672         for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01673 
01674         status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01675         // check status
01676 
01677         int nangles = 0;
01678         vector<float> anglelist;
01679         string angletype = "SPIDER";
01680         // Do we have a list of angles?
01681         if (params.has_key("anglelist")) {
01682                 anglelist = params["anglelist"];
01683                 nangles = anglelist.size() / 3;
01684         } else {
01685                 Transform* t3d = params["transform"];
01686                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01687                 // This part was modified by David Woolford -
01688                 // Before this the code worked only for SPIDER and EMAN angles,
01689                 // but the framework of the Transform3D allows for a generic implementation
01690                 // as specified here.
01691                 Dict p = t3d->get_rotation("spider");
01692                 if(t3d) {delete t3d; t3d=0;}
01693 
01694                 float phi   = p["phi"];
01695                 float theta = p["theta"];
01696                 float psi   = p["psi"];
01697                 anglelist.push_back(phi);
01698                 anglelist.push_back(theta);
01699                 anglelist.push_back(psi);
01700                 nangles = 1;
01701         }
01702         // End David Woolford modifications
01703 
01704         dm = new float[nangles*9];
01705         setdm(anglelist, angletype, dm);
01706 
01707                 // return images
01708         EMData *ret = new EMData();
01709         ret->set_size(nxvol, nyvol, nangles);
01710         ret->set_complex(false);
01711         ret->set_ri(true);
01712 
01713         images = ret->get_data();
01714 
01715         for (j = 1; j <= nangles; j++) {
01716                 status = fwdpj3(volsize, nrays, nnz   , &dm(1,j), origin, ri,
01717                                                 ptrs   ,  cord, sphere, &images(1,1,j));
01718         // check status?
01719         }
01720 
01721         // deallocate all temporary work space
01722         EMDeleteArray(dm);
01723         EMDeleteArray(ptrs);
01724         EMDeleteArray(cord);
01725         EMDeleteArray(sphere);
01726 
01727         if (!params.has_key("anglelist")) {
01728                 Transform* t3d = params["transform"];
01729                 ret->set_attr("xform.projection",t3d);
01730                 if(t3d) {delete t3d; t3d=0;}
01731         }
01732         ret->update();
01733         return ret;
01734 }
01735 
01736 
01737 #undef images
01738 
01739 #define images(i,j,k) images[ ((k)-1)*nximg*nyimg + ((j)-1)*nximg + (i)-1 ]
01740 // backproject from 2D to 3D (multiple images)
01741 EMData *ChaoProjector::backproject3d(EMData * imagestack) const
01742 {
01743         int nrays, nnz, status, j;
01744         float *dm;
01745         int   *ptrs, *cord;
01746         float *sphere, *images, *cube;
01747 
01748         int nximg   = imagestack->get_xsize();
01749         int nyimg   = imagestack->get_ysize();
01750         int nslices = imagestack->get_zsize();
01751 
01752         int dim = Util::get_min(nximg,nyimg);
01753         Vec3i volsize(nximg,nyimg,dim);
01754 
01755         Vec3i origin(0,0,0);
01756         // If a sensible origin isn't passed in, choose the middle of
01757         // the cube.
01758         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01759         else {origin[0] = nximg/2+1;}
01760         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01761         else {origin[1] = nyimg/2+1;}
01762         if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01763         else {origin[2] = dim/2+1;}
01764 
01765         int ri;
01766         if (params.has_key("radius")) {ri = params["radius"];}
01767         else {ri = dim/2 - 1;}
01768 
01769         // retrieve the voxel values
01770         images = imagestack->get_data();
01771 
01772         // count the number of voxels within a sphere centered at icent,
01773         // with radius ri
01774         status = getnnz(volsize, ri, origin, &nrays, &nnz);
01775         // need to check status...
01776 
01777         // convert from cube to sphere
01778         sphere = new float[nnz];
01779         ptrs   = new int[nrays+1];
01780         cord   = new int[3*nrays];
01781         if (sphere == NULL || ptrs == NULL || cord == NULL) {
01782                 fprintf(stderr,"ChaoProjector::backproject3d, failed to allocate!\n");
01783                 exit(1);
01784         }
01785         for (int i = 0; i<nnz; i++) sphere[i] = 0.0;
01786         for (int i = 0; i<nrays+1; i++) ptrs[i] = 0;
01787         for (int i = 0; i<3*nrays; i++) cord[i] = 0;
01788 
01789         int nangles = 0;
01790         vector<float> anglelist;
01791         string angletype = "SPIDER";
01792         // Do we have a list of angles?
01793         if (params.has_key("anglelist")) {
01794                 anglelist = params["anglelist"];
01795                 nangles = anglelist.size() / 3;
01796         } else {
01797                 Transform* t3d = params["transform"];
01798                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
01799                 // This part was modified by David Woolford -
01800                 // Before this the code worked only for SPIDER and EMAN angles,
01801                 // but the framework of the Transform3D allows for a generic implementation
01802                 // as specified here.
01803                 //  This was broken by david.  we need here a loop over all projections and put all angles on stack  PAP 06/28/09
01804                 Dict p = t3d->get_rotation("spider");
01805                 if(t3d) {delete t3d; t3d=0;}
01806 
01807                 float phi = p["phi"];
01808                 float theta = p["theta"];
01809                 float psi = p["psi"];
01810                 anglelist.push_back(phi);
01811                 anglelist.push_back(theta);
01812                 anglelist.push_back(psi);
01813                 nangles = 1;
01814         }
01815 
01816         // End David Woolford modifications
01817 
01818         if (nslices != nangles) {
01819                 LOGERR("the number of images does not match the number of angles");
01820                 return 0;
01821         }
01822 
01823         dm = new float[nangles*9];
01824         setdm(anglelist, angletype, dm);
01825 
01826         // return volume
01827         EMData *ret = new EMData();
01828         ret->set_size(nximg, nyimg, dim);
01829         ret->set_complex(false);
01830         ret->set_ri(true);
01831         ret->to_zero();
01832 
01833         cube = ret->get_data();
01834         // cb2sph should be replaced by something that touches only ptrs and cord
01835         status = cb2sph(cube, volsize, ri, origin, nnz, ptrs, cord, sphere);
01836         // check status
01837 
01838         for (j = 1; j <= nangles; j++) {
01839                 status = bckpj3(volsize, nrays, nnz, &dm(1,j), origin, ri,
01840                          ptrs   , cord , &images(1,1,j), sphere);
01841         // check status?
01842         }
01843 
01844         status = sph2cb(sphere, volsize, nrays, ri, nnz, ptrs, cord, cube);
01845         // check status?
01846 
01847         // deallocate all temporary work space
01848         EMDeleteArray(dm);
01849         EMDeleteArray(ptrs);
01850         EMDeleteArray(cord);
01851         EMDeleteArray(sphere);
01852 
01853         ret->update();
01854         return ret;
01855 }
01856 
01857 #undef images
01858 #undef cube
01859 #undef sphere
01860 #undef cord
01861 #undef ptrs
01862 #undef dm
01863 
01864 EMData *GaussFFTProjector::backproject3d(EMData * ) const
01865 {
01866     // no implementation yet
01867     EMData *ret = new EMData();
01868     return ret;
01869 }
01870 
01871 #define images(i,j,k) images[ (k)*nx*ny + ((j)-1)*nx + (i)-1 ]
01872 
01873 // EMData *PawelProjector::backproject3d(EMData * imagestack) const
01874 // {
01875 //
01876 //     float *images;
01877 //
01878 //     if (!imagestack) {
01879 //      return 0;
01880 //     }
01881 //     int ri;
01882 //     int nx      = imagestack->get_xsize();
01883 //     int ny      = imagestack->get_ysize();
01884 // //     int nslices = imagestack->get_zsize();
01885 //     int dim = Util::get_min(nx,ny);
01886 //     images  = imagestack->get_data();
01887 //
01888 //     Vec3i origin(0,0,0);
01889 //     // If a sensible origin isn't passed in, choose the middle of
01890 //     // the cube.
01891 //     if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
01892 //     else {origin[0] = nx/2;}
01893 //     if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
01894 //     else {origin[1] = ny/2;}
01895 //     if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
01896 //     else {origin[2] = dim/2;}
01897 //
01898 //     if (params.has_key("radius")) {ri = params["radius"];}
01899 //     else {ri = dim/2 - 1;}
01900 //
01901 //     // Determine the number of rows (x-lines) within the radius
01902 //     int nn = -1;
01903 //     prepcubes(nx, ny, dim, ri, origin, nn);
01904 //     // nn is now the number of rows-1 within the radius
01905 //     // so we can create and fill the ipcubes
01906 //     IPCube* ipcube = new IPCube[nn+1];
01907 //     prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
01908 //
01909 //      int nangles = 0;
01910 //      vector<float> anglelist;
01911 //      // Do we have a list of angles?
01912 //      if (params.has_key("anglelist")) {
01913 //              anglelist = params["anglelist"];
01914 //              nangles = anglelist.size() / 3;
01915 //      } else {
01916 //              Transform3D* t3d = params["t3d"];
01917 //              if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified");
01918 //              // This part was modified by David Woolford -
01919 //              // Before this the code worked only for SPIDER and EMAN angles,
01920 //              // but the framework of the Transform3D allows for a generic implementation
01921 //              // as specified here.
01922 //              Dict p = t3d->get_rotation(Transform3D::SPIDER);
01923 //
01924 //              string angletype = "SPIDER";
01925 //              float phi = p["phi"];
01926 //              float theta = p["theta"];
01927 //              float psi = p["psi"];
01928 //              anglelist.push_back(phi);
01929 //              anglelist.push_back(theta);
01930 //              anglelist.push_back(psi);
01931 //              nangles = 1;
01932 //      }
01933 //
01934 //      // End David Woolford modifications
01935 //
01936 //     // initialize return object
01937 //     EMData* ret = new EMData();
01938 //     ret->set_size(nx, ny, dim);
01939 //     ret->to_zero();
01940 //
01941 //     // loop over sets of angles
01942 //     for (int ia = 0; ia < nangles; ia++) {
01943 //        int indx = 3*ia;
01944 //         Transform3D rotation(Transform3D::SPIDER, float(anglelist[indx]),
01945 //                             float(anglelist[indx+1]),
01946 //                             float(anglelist[indx+2]));
01947 //        float dm1 = rotation.at(0,0);
01948 //        float dm4 = rotation.at(1,0);
01949 //
01950 //        if (2*(ri+1)+1 > dim) {
01951 //           // Must check x and y boundaries
01952 //           LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
01953 //           return 0;
01954 //        } else {
01955 //           // No need to check x and y boundaries
01956 //           for (int i = 0 ; i <= nn; i++) {
01957 //              int iox = (int)ipcube[i].loc[0]+origin[0];
01958 //              int ioy = (int)ipcube[i].loc[1]+origin[1];
01959 //              int ioz = (int)ipcube[i].loc[2]+origin[2];
01960 //
01961 //              Vec3f vb = rotation*ipcube[i].loc + origin;
01962 //              for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
01963 //                 float xbb = (j-ipcube[i].start)*dm1 + vb[0];
01964 //                 int   iqx = (int)floor(xbb);
01965 //
01966 //                 float ybb = (j-ipcube[i].start)*dm4 + vb[1];
01967 //                 int   iqy = (int)floor(ybb);
01968 //
01969 //                 float dipx = xbb - iqx;
01970 //                 float dipy = ybb - iqy;
01971 //
01972 //                 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
01973 //                     + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
01974 //                     + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
01975 //                     + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
01976 //                     - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
01977 //                 iox++;
01978 //              } // end for j
01979 //        } // end for i
01980 //        } // end if
01981 //     } // end for ia
01982 //
01983 //     ret->update();
01984 //     EMDeleteArray(ipcube);
01985 //     return ret;
01986 // }
01987 EMData *PawelProjector::backproject3d(EMData * imagestack) const
01988 {
01989 
01990         float *images;
01991 
01992         if (!imagestack) {
01993                 return 0;
01994         }
01995         int ri;
01996         int nx      = imagestack->get_xsize();
01997         int ny      = imagestack->get_ysize();
01998 //     int nslices = imagestack->get_zsize();
01999         int dim = Util::get_min(nx,ny);
02000         images  = imagestack->get_data();
02001 
02002         Vec3i origin(0,0,0);
02003     // If a sensible origin isn't passed in, choose the middle of
02004     // the cube.
02005         if (params.has_key("origin_x")) {origin[0] = params["origin_x"];}
02006         else {origin[0] = nx/2;}
02007         if (params.has_key("origin_y")) {origin[1] = params["origin_y"];}
02008         else {origin[1] = ny/2;}
02009         if (params.has_key("origin_z")) {origin[1] = params["origin_z"];}
02010         else {origin[2] = dim/2;}
02011 
02012         if (params.has_key("radius")) {ri = params["radius"];}
02013         else {ri = dim/2 - 1;}
02014 
02015     // Determine the number of rows (x-lines) within the radius
02016         int nn = -1;
02017         prepcubes(nx, ny, dim, ri, origin, nn);
02018     // nn is now the number of rows-1 within the radius
02019     // so we can create and fill the ipcubes
02020         IPCube* ipcube = new IPCube[nn+1];
02021         prepcubes(nx, ny, dim, ri, origin, nn, ipcube);
02022 
02023         int nangles = 0;
02024         vector<float> anglelist;
02025         // Do we have a list of angles?
02026         if (params.has_key("anglelist")) {
02027                 anglelist = params["anglelist"];
02028                 nangles = anglelist.size() / 3;
02029         } else {
02030                 Transform* t3d = params["transform"];
02031                 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified");
02032                 // This part was modified by David Woolford -
02033                 // Before this the code worked only for SPIDER and EMAN angles,
02034                 // but the framework of the Transform3D allows for a generic implementation
02035                 // as specified here.
02036                 Dict p = t3d->get_rotation("spider");
02037                 if(t3d) {delete t3d; t3d=0;}
02038 
02039                 string angletype = "SPIDER";
02040                 float phi = p["phi"];
02041                 float theta = p["theta"];
02042                 float psi = p["psi"];
02043                 anglelist.push_back(phi);
02044                 anglelist.push_back(theta);
02045                 anglelist.push_back(psi);
02046                 nangles = 1;
02047         }
02048 
02049         // End David Woolford modifications
02050 
02051     // initialize return object
02052         EMData* ret = new EMData();
02053         ret->set_size(nx, ny, dim);
02054         ret->to_zero();
02055 
02056     // loop over sets of angles
02057         for (int ia = 0; ia < nangles; ia++) {
02058                 int indx = 3*ia;
02059                 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]);
02060                 Transform rotation(d);
02061                 float dm1 = rotation.at(0,0);
02062                 float dm4 = rotation.at(1,0);
02063 
02064                 if (2*(ri+1)+1 > dim) {
02065           // Must check x and y boundaries
02066                         LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n");
02067                         return 0;
02068                 } else {
02069           // No need to check x and y boundaries
02070                         for (int i = 0 ; i <= nn; i++) {
02071                                 int iox = (int)ipcube[i].loc[0]+origin[0];
02072                                 int ioy = (int)ipcube[i].loc[1]+origin[1];
02073                                 int ioz = (int)ipcube[i].loc[2]+origin[2];
02074 
02075                                 Vec3f vb = rotation*ipcube[i].loc + origin;
02076                                 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) {
02077                                         float xbb = (j-ipcube[i].start)*dm1 + vb[0];
02078                                         int   iqx = (int)floor(xbb);
02079 
02080                                         float ybb = (j-ipcube[i].start)*dm4 + vb[1];
02081                                         int   iqy = (int)floor(ybb);
02082 
02083                                         float dipx = xbb - iqx;
02084                                         float dipy = ybb - iqy;
02085 
02086                                         (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia)
02087                                                         + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia))
02088                                                         + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia)
02089                                                                         + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia)
02090                                                                                         - images(iqx,iqy+1,ia)+images(iqx,iqy,ia)));
02091                                         iox++;
02092                                 } // end for j
02093                         } // end for i
02094                 } // end if
02095         } // end for ia
02096 
02097         ret->update();
02098         EMDeleteArray(ipcube);
02099         return ret;
02100 }
02101 #undef images
02102 
02103 EMData *StandardProjector::backproject3d(EMData * ) const
02104 {
02105    // no implementation yet
02106    EMData *ret = new EMData();
02107    return ret;
02108 }
02109 
02110 EMData *MaxValProjector::backproject3d(EMData * ) const
02111 {
02112    // no implementation yet
02113    EMData *ret = new EMData();
02114    return ret;
02115 }
02116 
02117 
02118 EMData *FourierGriddingProjector::backproject3d(EMData * ) const
02119 {
02120    // no implementation yet
02121    EMData *ret = new EMData();
02122    return ret;
02123 }
02124 
02125 // End Chao's projector addition 4/25/06
02126 
02127 void EMAN::dump_projectors()
02128 {
02129         dump_factory < Projector > ();
02130 }
02131 
02132 map<string, vector<string> > EMAN::dump_projectors_list()
02133 {
02134         return dump_factory_list < Projector > ();
02135 }
02136 
02137 /* vim: set ts=4 noet nospell: */