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

Generated on Thu May 3 10:06:27 2012 for EMAN2 by  doxygen 1.4.7