EMAN::PawelProjector Class Reference
Pawel Penczek's optimized projection routine. More...
#include <projector.h>
Inheritance diagram for EMAN::PawelProjector:
Collaboration diagram for EMAN::PawelProjector:
Public Member Functions | |
| EMData * | project3d (EMData *image) const |
| Project an 3D image into a 2D image. | |
| EMData * | backproject3d (EMData *image) const |
| Back-project a 2D image into a 3D image. | |
| string | get_name () const |
| Get the projector's name. | |
| string | get_desc () const |
| TypeDict | get_param_types () const |
| Get processor parameter information in a dictionary. | |
Static Public Member Functions | |
| static Projector * | NEW () |
Private Member Functions | |
| void | prepcubes (int nx, int ny, int nz, int ri, Vec3i origin, int &nn, IPCube *ipcube=NULL) const |
Classes | |
| struct | IPCube |
Detailed Description
Pawel Penczek's optimized projection routine.Definition at line 244 of file projector.h.
Member Function Documentation
Project an 3D image into a 2D image.
- Returns:
- A 2D image from the projection.
Implements EMAN::Projector.
Definition at line 712 of file projector.cpp.
References anglelist, EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Transform::get_matrix3_row(), EMAN::Util::get_min(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::Dict::has_key(), EMAN::PawelProjector::IPCube::loc, LOGERR, nn(), NullPointerException, nx, ny, EMAN::Projector::params, prepcubes(), EMAN::EMData::set_size(), EMAN::EMData::to_zero(), and EMAN::EMData::update().
00713 { 00714 if (!image) { 00715 return 0; 00716 } 00717 int ri; 00718 int nx = image->get_xsize(); 00719 int ny = image->get_ysize(); 00720 int nz = image->get_zsize(); 00721 int dim = Util::get_min(nx,ny,nz); 00722 if (nz == 1) { 00723 LOGERR("The PawelProjector needs a volume!"); 00724 return 0; 00725 } 00726 Vec3i origin(0,0,0); 00727 // If a sensible origin isn't passed in, choose the middle of 00728 // the cube. 00729 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];} 00730 else {origin[0] = nx/2;} 00731 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];} 00732 else {origin[1] = ny/2;} 00733 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];} 00734 else {origin[2] = nz/2;} 00735 00736 if (params.has_key("radius")) {ri = params["radius"];} 00737 else {ri = dim/2 - 1;} 00738 00739 // Determine the number of rows (x-lines) within the radius 00740 int nn = -1; 00741 prepcubes(nx, ny, nz, ri, origin, nn); 00742 // nn is now the number of rows-1 within the radius 00743 // so we can create and fill the ipcubes 00744 IPCube* ipcube = new IPCube[nn+1]; 00745 prepcubes(nx, ny, nz, ri, origin, nn, ipcube); 00746 00747 Transform* rotation = params["transform"]; 00748 int nangles = 0; 00749 vector<float> anglelist; 00750 // Do we have a list of angles? 00751 /* 00752 if (params.has_key("anglelist")) { 00753 anglelist = params["anglelist"]; 00754 nangles = anglelist.size() / 3; 00755 } else {*/ 00756 00757 if ( rotation == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified"); 00758 /* 00759 Dict p = t3d->get_rotation("spider"); 00760 00761 string angletype = "SPIDER"; 00762 float phi = p["phi"]; 00763 float theta = p["theta"]; 00764 float psi = p["psi"]; 00765 anglelist.push_back(phi); 00766 anglelist.push_back(theta); 00767 anglelist.push_back(psi); 00768 */ 00769 nangles = 1; 00770 //} 00771 00772 // initialize return object 00773 EMData* ret = new EMData(); 00774 ret->set_size(nx, ny, nangles); 00775 ret->to_zero(); 00776 00777 // loop over sets of angles 00778 for (int ia = 0; ia < nangles; ia++) { 00779 //int indx = 3*ia; 00780 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]); 00781 //Transform rotation(d); 00782 if (2*(ri+1)+1 > dim) { 00783 // Must check x and y boundaries 00784 for (int i = 0 ; i <= nn; i++) { 00785 int k = ipcube[i].loc[1] + origin[1]; 00786 Vec3f vb = ipcube[i].loc*(*rotation) + origin; 00787 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 00788 // check for pixels out-of-bounds 00789 int iox = int(vb[0]); 00790 if ((iox >= 0) && (iox < nx-1)) { 00791 int ioy = int(vb[1]); 00792 if ((ioy >= 0) && (ioy < ny-1)) { 00793 int ioz = int(vb[2]); 00794 if ((ioz >= 0) && (ioz < nz-1)) { 00795 // real work for pixels in bounds 00796 float dx = vb[0] - iox; 00797 float dy = vb[1] - ioy; 00798 float dz = vb[2] - ioz; 00799 float a1 = (*image)(iox,ioy,ioz); 00800 float a2 = (*image)(iox+1,ioy,ioz) - a1; 00801 float a3 = (*image)(iox,ioy+1,ioz) - a1; 00802 float a4 = (*image)(iox,ioy,ioz+1) - a1; 00803 float a5 = -a2 -(*image)(iox,ioy+1,ioz) 00804 + (*image)(iox+1,ioy+1,ioz); 00805 float a61 = -(*image)(iox,ioy,ioz+1) 00806 + (*image)(iox+1,ioy,ioz+1); 00807 float a6 = -a2 + a61; 00808 float a7 = -a3 - (*image)(iox,ioy,ioz+1) 00809 + (*image)(iox,ioy+1,ioz+1); 00810 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1) 00811 + (*image)(iox+1,ioy+1,ioz+1); 00812 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx 00813 + (a7 + a8*dx)*dy) 00814 + a3*dy + dx*(a2 + a5*dy); 00815 } 00816 } 00817 } 00818 vb += rotation->get_matrix3_row(0); 00819 } 00820 } 00821 00822 } else { 00823 // No need to check x and y boundaries 00824 for (int i = 0 ; i <= nn; i++) { 00825 int k = ipcube[i].loc[1] + origin[1]; 00826 Vec3f vb = ipcube[i].loc*(*rotation) + origin; 00827 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 00828 int iox = int(vb[0]); 00829 int ioy = int(vb[1]); 00830 int ioz = int(vb[2]); 00831 float dx = vb[0] - iox; 00832 float dy = vb[1] - ioy; 00833 float dz = vb[2] - ioz; 00834 float a1 = (*image)(iox,ioy,ioz); 00835 float a2 = (*image)(iox+1,ioy,ioz) - a1; 00836 float a3 = (*image)(iox,ioy+1,ioz) - a1; 00837 float a4 = (*image)(iox,ioy,ioz+1) - a1; 00838 float a5 = -a2 -(*image)(iox,ioy+1,ioz) 00839 + (*image)(iox+1,ioy+1,ioz); 00840 float a61 = -(*image)(iox,ioy,ioz+1) 00841 + (*image)(iox+1,ioy,ioz+1); 00842 float a6 = -a2 + a61; 00843 float a7 = -a3 - (*image)(iox,ioy,ioz+1) 00844 + (*image)(iox,ioy+1,ioz+1); 00845 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1) 00846 + (*image)(iox+1,ioy+1,ioz+1); 00847 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx 00848 + (a7 + a8*dx)*dy) 00849 + a3*dy + dx*(a2 + a5*dy); 00850 vb += rotation->get_matrix3_row(0); 00851 } 00852 } 00853 } 00854 } 00855 ret->update(); 00856 EMDeleteArray(ipcube); 00857 if(rotation) {delete rotation; rotation=0;} 00858 return ret; 00859 }
Back-project a 2D image into a 3D image.
- Returns:
- A 3D image from the backprojection.
Implements EMAN::Projector.
Definition at line 1950 of file projector.cpp.
References anglelist, EMAN::Transform::at(), EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::EMData::get_data(), EMAN::Util::get_min(), EMAN::Transform::get_rotation(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::Dict::has_key(), images, EMAN::PawelProjector::IPCube::loc, LOGERR, nn(), NullPointerException, nx, ny, EMAN::Projector::params, phi, prepcubes(), EMAN::EMData::set_size(), EMAN::Dict::size(), EMAN::PawelProjector::IPCube::start, theta, EMAN::EMData::to_zero(), and EMAN::EMData::update().
01951 { 01952 01953 float *images; 01954 01955 if (!imagestack) { 01956 return 0; 01957 } 01958 int ri; 01959 int nx = imagestack->get_xsize(); 01960 int ny = imagestack->get_ysize(); 01961 // int nslices = imagestack->get_zsize(); 01962 int dim = Util::get_min(nx,ny); 01963 images = imagestack->get_data(); 01964 01965 Vec3i origin(0,0,0); 01966 // If a sensible origin isn't passed in, choose the middle of 01967 // the cube. 01968 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];} 01969 else {origin[0] = nx/2;} 01970 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];} 01971 else {origin[1] = ny/2;} 01972 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];} 01973 else {origin[2] = dim/2;} 01974 01975 if (params.has_key("radius")) {ri = params["radius"];} 01976 else {ri = dim/2 - 1;} 01977 01978 // Determine the number of rows (x-lines) within the radius 01979 int nn = -1; 01980 prepcubes(nx, ny, dim, ri, origin, nn); 01981 // nn is now the number of rows-1 within the radius 01982 // so we can create and fill the ipcubes 01983 IPCube* ipcube = new IPCube[nn+1]; 01984 prepcubes(nx, ny, dim, ri, origin, nn, ipcube); 01985 01986 int nangles = 0; 01987 vector<float> anglelist; 01988 // Do we have a list of angles? 01989 if (params.has_key("anglelist")) { 01990 anglelist = params["anglelist"]; 01991 nangles = anglelist.size() / 3; 01992 } else { 01993 Transform* t3d = params["transform"]; 01994 if ( t3d == NULL ) throw NullPointerException("The transform3d object (required for projection), was not specified"); 01995 // This part was modified by David Woolford - 01996 // Before this the code worked only for SPIDER and EMAN angles, 01997 // but the framework of the Transform3D allows for a generic implementation 01998 // as specified here. 01999 Dict p = t3d->get_rotation("spider"); 02000 if(t3d) {delete t3d; t3d=0;} 02001 02002 string angletype = "SPIDER"; 02003 float phi = p["phi"]; 02004 float theta = p["theta"]; 02005 float psi = p["psi"]; 02006 anglelist.push_back(phi); 02007 anglelist.push_back(theta); 02008 anglelist.push_back(psi); 02009 nangles = 1; 02010 } 02011 02012 // End David Woolford modifications 02013 02014 // initialize return object 02015 EMData* ret = new EMData(); 02016 ret->set_size(nx, ny, dim); 02017 ret->to_zero(); 02018 02019 // loop over sets of angles 02020 for (int ia = 0; ia < nangles; ia++) { 02021 int indx = 3*ia; 02022 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]); 02023 Transform rotation(d); 02024 float dm1 = rotation.at(0,0); 02025 float dm4 = rotation.at(1,0); 02026 02027 if (2*(ri+1)+1 > dim) { 02028 // Must check x and y boundaries 02029 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n"); 02030 return 0; 02031 } else { 02032 // No need to check x and y boundaries 02033 for (int i = 0 ; i <= nn; i++) { 02034 int iox = (int)ipcube[i].loc[0]+origin[0]; 02035 int ioy = (int)ipcube[i].loc[1]+origin[1]; 02036 int ioz = (int)ipcube[i].loc[2]+origin[2]; 02037 02038 Vec3f vb = rotation*ipcube[i].loc + origin; 02039 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 02040 float xbb = (j-ipcube[i].start)*dm1 + vb[0]; 02041 int iqx = (int)floor(xbb); 02042 02043 float ybb = (j-ipcube[i].start)*dm4 + vb[1]; 02044 int iqy = (int)floor(ybb); 02045 02046 float dipx = xbb - iqx; 02047 float dipy = ybb - iqy; 02048 02049 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia) 02050 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia)) 02051 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia) 02052 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia) 02053 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia))); 02054 iox++; 02055 } // end for j 02056 } // end for i 02057 } // end if 02058 } // end for ia 02059 02060 ret->update(); 02061 EMDeleteArray(ipcube); 02062 return ret; 02063 }
| string EMAN::PawelProjector::get_name | ( | ) | const [inline, virtual] |
Get the projector's name.
Each projector is indentified by unique name.
- Returns:
- The projector's name.
Implements EMAN::Projector.
Definition at line 250 of file projector.h.
| string EMAN::PawelProjector::get_desc | ( | ) | const [inline, virtual] |
Implements EMAN::Projector.
Definition at line 255 of file projector.h.
00256 { 00257 return "Pawel Penczek's optimized real-space projection generation. Minor interpolation artifacts."; 00258 }
| static Projector* EMAN::PawelProjector::NEW | ( | ) | [inline, static] |
| TypeDict EMAN::PawelProjector::get_param_types | ( | ) | const [inline, virtual] |
Get processor parameter information in a dictionary.
Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.
- Returns:
- A dictionary containing the parameter info.
Reimplemented from EMAN::Projector.
Definition at line 265 of file projector.h.
References EMAN::EMObject::FLOAT, EMAN::EMObject::FLOATARRAY, EMAN::EMObject::INT, EMAN::TypeDict::put(), EMAN::EMObject::STRING, and EMAN::EMObject::TRANSFORM.
00266 { 00267 TypeDict d; 00268 d.put("transform", EMObject::TRANSFORM); 00269 d.put("origin_x", EMObject::INT); 00270 d.put("origin_y", EMObject::INT); 00271 d.put("origin_z", EMObject::INT); 00272 d.put("radius", EMObject::INT); 00273 d.put("anglelist", EMObject::FLOATARRAY); 00274 d.put("angletype", EMObject::STRING); 00275 d.put("theta", EMObject::FLOAT); 00276 d.put("psi", EMObject::FLOAT); 00277 return d; 00278 }
| void PawelProjector::prepcubes | ( | int | nx, | |
| int | ny, | |||
| int | nz, | |||
| int | ri, | |||
| Vec3i | origin, | |||
| int & | nn, | |||
| IPCube * | ipcube = NULL | |||
| ) | const [private] |
Definition at line 520 of file projector.cpp.
References EMAN::PawelProjector::IPCube::end, EMAN::PawelProjector::IPCube::loc, EMAN::PawelProjector::IPCube::start, and t.
Referenced by backproject3d(), and project3d().
00521 { 00522 const float r = float(ri*ri); 00523 const int ldpx = origin[0]; 00524 const int ldpy = origin[1]; 00525 const int ldpz = origin[2]; 00526 float t; 00527 nn = -1; 00528 for (int i1 = 0; i1 < nz; i1++) { 00529 t = float(i1 - ldpz); 00530 const float xx = t*t; 00531 for (int i2 = 0; i2 < ny; i2++) { 00532 t = float(i2 - ldpy); 00533 const float yy = t*t + xx; 00534 bool first = true; 00535 for (int i3 = 0; i3 < nz; i3++) { 00536 t = float(i3 - ldpx); 00537 const float rc = t*t + yy; 00538 if (first) { 00539 // first pixel on this line 00540 if (rc > r) continue; 00541 first = false; 00542 nn++; 00543 if (ipcube != NULL) { 00544 ipcube[nn].start = i3; 00545 ipcube[nn].end = i3; 00546 ipcube[nn].loc[0] = i3 - ldpx; 00547 ipcube[nn].loc[1] = i2 - ldpy; 00548 ipcube[nn].loc[2] = i1 - ldpz; 00549 } 00550 } else { 00551 // second or later pixel on this line 00552 if (ipcube != NULL) { 00553 if (rc <= r) ipcube[nn].end = i3; 00554 } 00555 } 00556 } 00557 } 00558 } 00559 }
The documentation for this class was generated from the following files:
