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

projector.cpp

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

Generated on Fri Apr 30 15:38:56 2010 for EMAN2 by  doxygen 1.3.9.1