#include <projector.h>
Inheritance 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 () |
Static Public Attributes | |
static const string | NAME = "pawel" |
Private Member Functions | |
void | prepcubes (int nx, int ny, int nz, int ri, Vec3i origin, int &nn, IPCube *ipcube=NULL) const |
Classes | |
struct | IPCube |
Definition at line 248 of file projector.h.
Back-project a 2D image into a 3D image.
Implements EMAN::Projector.
Definition at line 1972 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().
01973 { 01974 01975 float *images; 01976 01977 if (!imagestack) { 01978 return 0; 01979 } 01980 int ri; 01981 int nx = imagestack->get_xsize(); 01982 int ny = imagestack->get_ysize(); 01983 // int nslices = imagestack->get_zsize(); 01984 int dim = Util::get_min(nx,ny); 01985 images = imagestack->get_data(); 01986 01987 Vec3i origin(0,0,0); 01988 // If a sensible origin isn't passed in, choose the middle of 01989 // the cube. 01990 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];} 01991 else {origin[0] = nx/2;} 01992 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];} 01993 else {origin[1] = ny/2;} 01994 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];} 01995 else {origin[2] = dim/2;} 01996 01997 if (params.has_key("radius")) {ri = params["radius"];} 01998 else {ri = dim/2 - 1;} 01999 02000 // Determine the number of rows (x-lines) within the radius 02001 int nn = -1; 02002 prepcubes(nx, ny, dim, ri, origin, nn); 02003 // nn is now the number of rows-1 within the radius 02004 // so we can create and fill the ipcubes 02005 IPCube* ipcube = new IPCube[nn+1]; 02006 prepcubes(nx, ny, dim, ri, origin, nn, ipcube); 02007 02008 int nangles = 0; 02009 vector<float> anglelist; 02010 // Do we have a list of angles? 02011 if (params.has_key("anglelist")) { 02012 anglelist = params["anglelist"]; 02013 nangles = anglelist.size() / 3; 02014 } else { 02015 Transform* t3d = params["transform"]; 02016 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified"); 02017 // This part was modified by David Woolford - 02018 // Before this the code worked only for SPIDER and EMAN angles, 02019 // but the framework of the Transform3D allows for a generic implementation 02020 // as specified here. 02021 Dict p = t3d->get_rotation("spider"); 02022 if(t3d) {delete t3d; t3d=0;} 02023 02024 string angletype = "SPIDER"; 02025 float phi = p["phi"]; 02026 float theta = p["theta"]; 02027 float psi = p["psi"]; 02028 anglelist.push_back(phi); 02029 anglelist.push_back(theta); 02030 anglelist.push_back(psi); 02031 nangles = 1; 02032 } 02033 02034 // End David Woolford modifications 02035 02036 // initialize return object 02037 EMData* ret = new EMData(); 02038 ret->set_size(nx, ny, dim); 02039 ret->to_zero(); 02040 02041 // loop over sets of angles 02042 for (int ia = 0; ia < nangles; ia++) { 02043 int indx = 3*ia; 02044 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]); 02045 Transform rotation(d); 02046 float dm1 = rotation.at(0,0); 02047 float dm4 = rotation.at(1,0); 02048 02049 if (2*(ri+1)+1 > dim) { 02050 // Must check x and y boundaries 02051 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n"); 02052 return 0; 02053 } else { 02054 // No need to check x and y boundaries 02055 for (int i = 0 ; i <= nn; i++) { 02056 int iox = (int)ipcube[i].loc[0]+origin[0]; 02057 int ioy = (int)ipcube[i].loc[1]+origin[1]; 02058 int ioz = (int)ipcube[i].loc[2]+origin[2]; 02059 02060 Vec3f vb = rotation*ipcube[i].loc + origin; 02061 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 02062 float xbb = (j-ipcube[i].start)*dm1 + vb[0]; 02063 int iqx = (int)floor(xbb); 02064 02065 float ybb = (j-ipcube[i].start)*dm4 + vb[1]; 02066 int iqy = (int)floor(ybb); 02067 02068 float dipx = xbb - iqx; 02069 float dipy = ybb - iqy; 02070 02071 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia) 02072 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia)) 02073 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia) 02074 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia) 02075 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia))); 02076 iox++; 02077 } // end for j 02078 } // end for i 02079 } // end if 02080 } // end for ia 02081 02082 ret->update(); 02083 EMDeleteArray(ipcube); 02084 return ret; 02085 }
string EMAN::PawelProjector::get_desc | ( | ) | const [inline, virtual] |
Implements EMAN::Projector.
Definition at line 259 of file projector.h.
00260 { 00261 return "Pawel Penczek's optimized real-space projection generation. Minor interpolation artifacts."; 00262 }
string EMAN::PawelProjector::get_name | ( | ) | const [inline, virtual] |
Get the projector's name.
Each projector is indentified by unique name.
Implements EMAN::Projector.
Definition at line 254 of file projector.h.
References NAME.
00255 { 00256 return NAME; 00257 }
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.
Reimplemented from EMAN::Projector.
Definition at line 269 of file projector.h.
References EMAN::EMObject::FLOAT, EMAN::EMObject::FLOATARRAY, EMAN::EMObject::INT, EMAN::TypeDict::put(), EMAN::EMObject::STRING, and EMAN::EMObject::TRANSFORM.
00270 { 00271 TypeDict d; 00272 d.put("transform", EMObject::TRANSFORM); 00273 d.put("origin_x", EMObject::INT); 00274 d.put("origin_y", EMObject::INT); 00275 d.put("origin_z", EMObject::INT); 00276 d.put("radius", EMObject::INT); 00277 d.put("anglelist", EMObject::FLOATARRAY); 00278 d.put("angletype", EMObject::STRING); 00279 d.put("theta", EMObject::FLOAT); 00280 d.put("psi", EMObject::FLOAT); 00281 return d; 00282 }
static Projector* EMAN::PawelProjector::NEW | ( | ) | [inline, static] |
void PawelProjector::prepcubes | ( | int | nx, | |
int | ny, | |||
int | nz, | |||
int | ri, | |||
Vec3i | origin, | |||
int & | nn, | |||
IPCube * | ipcube = NULL | |||
) | const [private] |
Definition at line 531 of file projector.cpp.
References EMAN::PawelProjector::IPCube::end, EMAN::PawelProjector::IPCube::loc, and t.
Referenced by backproject3d(), and project3d().
00532 { 00533 const float r = float(ri)*float(ri); 00534 const int ldpx = origin[0]; 00535 const int ldpy = origin[1]; 00536 const int ldpz = origin[2]; 00537 //cout<<" ldpx "<<ldpx<<" i2 "<<ldpy<<" i1 "<<ldpz<<endl; 00538 float t; 00539 nn = -1; 00540 for (int i1 = 0; i1 < nz; i1++) { 00541 t = float(i1 - ldpz); 00542 const float xx = t*t; 00543 for (int i2 = 0; i2 < ny; i2++) { 00544 t = float(i2 - ldpy); 00545 const float yy = t*t + xx; 00546 bool first = true; 00547 for (int i3 = 0; i3 < nx; i3++) { 00548 t = float(i3 - ldpx); 00549 const float rc = t*t + yy; 00550 if (first) { 00551 // first pixel on this line 00552 if (rc > r) continue; 00553 first = false; 00554 nn++; 00555 if (ipcube != NULL) { 00556 ipcube[nn].start = i3; 00557 ipcube[nn].end = i3; 00558 ipcube[nn].loc[0] = i3 - ldpx; 00559 ipcube[nn].loc[1] = i2 - ldpy; 00560 ipcube[nn].loc[2] = i1 - ldpz; 00561 //cout<<" start "<<i3<<" i2 "<<i2<<" i1 "<<i1<<endl; 00562 } 00563 } else { 00564 // second or later pixel on this line 00565 if (ipcube != NULL) { 00566 if (rc <= r) ipcube[nn].end = i3; 00567 } 00568 } 00569 } 00570 } 00571 } 00572 }
Project an 3D image into a 2D image.
Implements EMAN::Projector.
Definition at line 575 of file projector.cpp.
References anglelist, EMDeleteArray(), EMAN::PawelProjector::IPCube::end, EMAN::Transform::get_matrix3_row(), EMAN::Util::get_max(), 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().
00576 { 00577 if (!image) { 00578 return 0; 00579 } 00580 int ri; 00581 int nx = image->get_xsize(); 00582 int ny = image->get_ysize(); 00583 int nz = image->get_zsize(); 00584 int dim = Util::get_max(nx,ny,nz); 00585 if (nz == 1) { 00586 LOGERR("The PawelProjector needs a volume!"); 00587 return 0; 00588 } 00589 00590 Vec3i origin(0,0,0); 00591 // If a sensible origin isn't passed in, choose the middle of 00592 // the cube. 00593 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];} 00594 else {origin[0] = nx/2;} 00595 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];} 00596 else {origin[1] = ny/2;} 00597 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];} 00598 else {origin[2] = nz/2;} 00599 00600 if (params.has_key("radius")) {ri = params["radius"];} 00601 else {ri = dim/2 - 1;} 00602 00603 // Determine the number of rows (x-lines) within the radius 00604 int nn = -1; 00605 prepcubes(nx, ny, nz, ri, origin, nn); 00606 // nn is now the number of rows-1 within the radius 00607 // so we can create and fill the ipcubes 00608 IPCube* ipcube = new IPCube[nn+1]; 00609 prepcubes(nx, ny, nz, ri, origin, nn, ipcube); 00610 00611 Transform* rotation = params["transform"]; 00612 int nangles = 0; 00613 vector<float> anglelist; 00614 // Do we have a list of angles? 00615 /* 00616 if (params.has_key("anglelist")) { 00617 anglelist = params["anglelist"]; 00618 nangles = anglelist.size() / 3; 00619 } else {*/ 00620 00621 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified"); 00622 /* 00623 Dict p = t3d->get_rotation("spider"); 00624 00625 string angletype = "SPIDER"; 00626 float phi = p["phi"]; 00627 float theta = p["theta"]; 00628 float psi = p["psi"]; 00629 anglelist.push_back(phi); 00630 anglelist.push_back(theta); 00631 anglelist.push_back(psi); 00632 */ 00633 nangles = 1; 00634 //} 00635 00636 // initialize return object 00637 EMData* ret = new EMData(); 00638 ret->set_size(nx, ny, nangles); 00639 ret->to_zero(); 00640 //for (int i = 0 ; i <= nn; i++) cout<<" IPCUBE "<<ipcube[i].start<<" "<<ipcube[i].end<<" "<<ipcube[i].loc[0]<<" "<<ipcube[i].loc[1]<<" "<<ipcube[i].loc[2]<<endl; 00641 // loop over sets of angles 00642 for (int ia = 0; ia < nangles; ia++) { 00643 //int indx = 3*ia; 00644 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]); 00645 //Transform rotation(d); 00646 if (2*(ri+1)+1 > dim) { 00647 // Must check x and y boundaries 00648 for (int i = 0 ; i <= nn; i++) { 00649 int k = ipcube[i].loc[1] + origin[1]; 00650 Vec3f vb = ipcube[i].loc*(*rotation) + origin; 00651 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 00652 // check for pixels out-of-bounds 00653 int iox = int(vb[0]); 00654 if ((iox >= 0) && (iox < nx-1)) { 00655 int ioy = int(vb[1]); 00656 if ((ioy >= 0) && (ioy < ny-1)) { 00657 int ioz = int(vb[2]); 00658 if ((ioz >= 0) && (ioz < nz-1)) { 00659 // real work for pixels in bounds 00660 //cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<endl; 00661 //cout<<" TAKE"<<endl; 00662 float dx = vb[0] - iox; 00663 float dy = vb[1] - ioy; 00664 float dz = vb[2] - ioz; 00665 float a1 = (*image)(iox,ioy,ioz); 00666 float a2 = (*image)(iox+1,ioy,ioz) - a1; 00667 float a3 = (*image)(iox,ioy+1,ioz) - a1; 00668 float a4 = (*image)(iox,ioy,ioz+1) - a1; 00669 float a5 = -a2 -(*image)(iox,ioy+1,ioz) 00670 + (*image)(iox+1,ioy+1,ioz); 00671 float a61 = -(*image)(iox,ioy,ioz+1) 00672 + (*image)(iox+1,ioy,ioz+1); 00673 float a6 = -a2 + a61; 00674 float a7 = -a3 - (*image)(iox,ioy,ioz+1) 00675 + (*image)(iox,ioy+1,ioz+1); 00676 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1) 00677 + (*image)(iox+1,ioy+1,ioz+1); 00678 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx 00679 + (a7 + a8*dx)*dy) 00680 + a3*dy + dx*(a2 + a5*dy); 00681 } //else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLAATED Z "<<endl; } 00682 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED Y "<<endl; } 00683 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED X "<<endl; } 00684 vb += rotation->get_matrix3_row(0); 00685 } 00686 } 00687 00688 } else { 00689 // No need to check x and y boundaries 00690 for (int i = 0 ; i <= nn; i++) { 00691 int k = ipcube[i].loc[1] + origin[1]; 00692 Vec3f vb = ipcube[i].loc*(*rotation) + origin; 00693 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 00694 int iox = int(vb[0]); 00695 int ioy = int(vb[1]); 00696 int ioz = int(vb[2]); 00697 float dx = vb[0] - iox; 00698 float dy = vb[1] - ioy; 00699 float dz = vb[2] - ioz; 00700 float a1 = (*image)(iox,ioy,ioz); 00701 float a2 = (*image)(iox+1,ioy,ioz) - a1; 00702 float a3 = (*image)(iox,ioy+1,ioz) - a1; 00703 float a4 = (*image)(iox,ioy,ioz+1) - a1; 00704 float a5 = -a2 -(*image)(iox,ioy+1,ioz) 00705 + (*image)(iox+1,ioy+1,ioz); 00706 float a61 = -(*image)(iox,ioy,ioz+1) 00707 + (*image)(iox+1,ioy,ioz+1); 00708 float a6 = -a2 + a61; 00709 float a7 = -a3 - (*image)(iox,ioy,ioz+1) 00710 + (*image)(iox,ioy+1,ioz+1); 00711 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1) 00712 + (*image)(iox+1,ioy+1,ioz+1); 00713 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx 00714 + (a7 + a8*dx)*dy) 00715 + a3*dy + dx*(a2 + a5*dy); 00716 vb += rotation->get_matrix3_row(0); 00717 } 00718 } 00719 } 00720 } 00721 ret->update(); 00722 EMDeleteArray(ipcube); 00723 if(rotation) {delete rotation; rotation=0;} 00724 return ret; 00725 }
const string PawelProjector::NAME = "pawel" [static] |