#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 | |
Projector * | NEW () |
Static Public Attributes | |
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 |
Definition at line 248 of file projector.h.
|
Back-project a 2D image into a 3D image.
Implements EMAN::Projector. Definition at line 2127 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, phi, prepcubes(), EMAN::EMData::set_size(), EMAN::Dict::size(), EMAN::PawelProjector::IPCube::start, theta, EMAN::EMData::to_zero(), EMAN::EMData::update(), EMAN::Vec3f, and EMAN::Vec3i. 02128 { 02129 02130 float *images; 02131 02132 if (!imagestack) { 02133 return 0; 02134 } 02135 int ri; 02136 int nx = imagestack->get_xsize(); 02137 int ny = imagestack->get_ysize(); 02138 // int nslices = imagestack->get_zsize(); 02139 int dim = Util::get_min(nx,ny); 02140 images = imagestack->get_data(); 02141 02142 Vec3i origin(0,0,0); 02143 // If a sensible origin isn't passed in, choose the middle of 02144 // the cube. 02145 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];} 02146 else {origin[0] = nx/2;} 02147 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];} 02148 else {origin[1] = ny/2;} 02149 if (params.has_key("origin_z")) {origin[1] = params["origin_z"];} 02150 else {origin[2] = dim/2;} 02151 02152 if (params.has_key("radius")) {ri = params["radius"];} 02153 else {ri = dim/2 - 1;} 02154 02155 // Determine the number of rows (x-lines) within the radius 02156 int nn = -1; 02157 prepcubes(nx, ny, dim, ri, origin, nn); 02158 // nn is now the number of rows-1 within the radius 02159 // so we can create and fill the ipcubes 02160 IPCube* ipcube = new IPCube[nn+1]; 02161 prepcubes(nx, ny, dim, ri, origin, nn, ipcube); 02162 02163 int nangles = 0; 02164 vector<float> anglelist; 02165 // Do we have a list of angles? 02166 if (params.has_key("anglelist")) { 02167 anglelist = params["anglelist"]; 02168 nangles = anglelist.size() / 3; 02169 } else { 02170 Transform* t3d = params["transform"]; 02171 if ( t3d == NULL ) throw NullPointerException("The transform object (required for projection), was not specified"); 02172 // This part was modified by David Woolford - 02173 // Before this the code worked only for SPIDER and EMAN angles, 02174 // but the framework of the Transform3D allows for a generic implementation 02175 // as specified here. 02176 Dict p = t3d->get_rotation("spider"); 02177 if(t3d) {delete t3d; t3d=0;} 02178 02179 string angletype = "SPIDER"; 02180 float phi = p["phi"]; 02181 float theta = p["theta"]; 02182 float psi = p["psi"]; 02183 anglelist.push_back(phi); 02184 anglelist.push_back(theta); 02185 anglelist.push_back(psi); 02186 nangles = 1; 02187 } 02188 02189 // End David Woolford modifications 02190 02191 // initialize return object 02192 EMData* ret = new EMData(); 02193 ret->set_size(nx, ny, dim); 02194 ret->to_zero(); 02195 02196 // loop over sets of angles 02197 for (int ia = 0; ia < nangles; ia++) { 02198 int indx = 3*ia; 02199 Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]); 02200 Transform rotation(d); 02201 float dm1 = rotation.at(0,0); 02202 float dm4 = rotation.at(1,0); 02203 02204 if (2*(ri+1)+1 > dim) { 02205 // Must check x and y boundaries 02206 LOGERR("backproject3d, pawel, 2*(ri+1)+1 > dim\n"); 02207 return 0; 02208 } else { 02209 // No need to check x and y boundaries 02210 for (int i = 0 ; i <= nn; i++) { 02211 int iox = (int)ipcube[i].loc[0]+origin[0]; 02212 int ioy = (int)ipcube[i].loc[1]+origin[1]; 02213 int ioz = (int)ipcube[i].loc[2]+origin[2]; 02214 02215 Vec3f vb = rotation*ipcube[i].loc + origin; 02216 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 02217 float xbb = (j-ipcube[i].start)*dm1 + vb[0]; 02218 int iqx = (int)floor(xbb); 02219 02220 float ybb = (j-ipcube[i].start)*dm4 + vb[1]; 02221 int iqy = (int)floor(ybb); 02222 02223 float dipx = xbb - iqx; 02224 float dipy = ybb - iqy; 02225 02226 (*ret)(iox,ioy,ioz) += images(iqx,iqy,ia) 02227 + dipy*(images(iqx,iqy+1,ia)-images(iqx,iqy,ia)) 02228 + dipx*(images(iqx+1,iqy,ia)-images(iqx,iqy,ia) 02229 + dipy*(images(iqx+1,iqy+1,ia)-images(iqx+1,iqy,ia) 02230 - images(iqx,iqy+1,ia)+images(iqx,iqy,ia))); 02231 iox++; 02232 } // end for j 02233 } // end for i 02234 } // end if 02235 } // end for ia 02236 02237 ret->update(); 02238 EMDeleteArray(ipcube); 02239 return ret; 02240 }
|
|
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 }
|
|
Get the projector's name. Each projector is indentified by unique name.
Implements EMAN::Projector. Definition at line 254 of file projector.h. 00255 {
00256 return NAME;
00257 }
|
|
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::TypeDict::put(). 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 }
|
|
Definition at line 264 of file projector.h. 00265 { 00266 return new PawelProjector(); 00267 }
|
|
Definition at line 533 of file projector.cpp. References EMAN::PawelProjector::IPCube::end, EMAN::PawelProjector::IPCube::loc, nn(), EMAN::PawelProjector::IPCube::start, t, and EMAN::Vec3i. Referenced by backproject3d(), and project3d(). 00534 { 00535 const float r = float(ri)*float(ri); 00536 const int ldpx = origin[0]; 00537 const int ldpy = origin[1]; 00538 const int ldpz = origin[2]; 00539 //cout<<" ldpx "<<ldpx<<" i2 "<<ldpy<<" i1 "<<ldpz<<endl; 00540 float t; 00541 nn = -1; 00542 for (int i1 = 0; i1 < nz; i1++) { 00543 t = float(i1 - ldpz); 00544 const float xx = t*t; 00545 for (int i2 = 0; i2 < ny; i2++) { 00546 t = float(i2 - ldpy); 00547 const float yy = t*t + xx; 00548 bool first = true; 00549 for (int i3 = 0; i3 < nx; i3++) { 00550 t = float(i3 - ldpx); 00551 const float rc = t*t + yy; 00552 if (first) { 00553 // first pixel on this line 00554 if (rc > r) continue; 00555 first = false; 00556 nn++; 00557 if (ipcube != NULL) { 00558 ipcube[nn].start = i3; 00559 ipcube[nn].end = i3; 00560 ipcube[nn].loc[0] = i3 - ldpx; 00561 ipcube[nn].loc[1] = i2 - ldpy; 00562 ipcube[nn].loc[2] = i1 - ldpz; 00563 //cout<<" start "<<i3<<" i2 "<<i2<<" i1 "<<i1<<endl; 00564 } 00565 } else { 00566 // second or later pixel on this line 00567 if (ipcube != NULL) { 00568 if (rc <= r) ipcube[nn].end = i3; 00569 } 00570 } 00571 } 00572 } 00573 } 00574 }
|
|
Project an 3D image into a 2D image.
Implements EMAN::Projector. Definition at line 577 of file projector.cpp. References 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, prepcubes(), EMAN::EMData::set_size(), EMAN::PawelProjector::IPCube::start, EMAN::EMData::to_zero(), EMAN::EMData::update(), EMAN::Vec3f, and EMAN::Vec3i. 00578 { 00579 if (!image) { 00580 return 0; 00581 } 00582 int ri; 00583 int nx = image->get_xsize(); 00584 int ny = image->get_ysize(); 00585 int nz = image->get_zsize(); 00586 int dim = Util::get_max(nx,ny,nz); 00587 if (nz == 1) { 00588 LOGERR("The PawelProjector needs a volume!"); 00589 return 0; 00590 } 00591 00592 Vec3i origin(0,0,0); 00593 // If a sensible origin isn't passed in, choose the middle of 00594 // the cube. 00595 if (params.has_key("origin_x")) {origin[0] = params["origin_x"];} 00596 else {origin[0] = nx/2;} 00597 if (params.has_key("origin_y")) {origin[1] = params["origin_y"];} 00598 else {origin[1] = ny/2;} 00599 if (params.has_key("origin_z")) {origin[2] = params["origin_z"];} 00600 else {origin[2] = nz/2;} 00601 00602 if (params.has_key("radius")) {ri = params["radius"];} 00603 else {ri = dim/2 - 1;} 00604 00605 // Determine the number of rows (x-lines) within the radius 00606 int nn = -1; 00607 prepcubes(nx, ny, nz, ri, origin, nn); 00608 // nn is now the number of rows-1 within the radius 00609 // so we can create and fill the ipcubes 00610 IPCube* ipcube = new IPCube[nn+1]; 00611 prepcubes(nx, ny, nz, ri, origin, nn, ipcube); 00612 00613 Transform* rotation = params["transform"]; 00614 int nangles = 0; 00615 vector<float> anglelist; 00616 // Do we have a list of angles? 00617 /* 00618 if (params.has_key("anglelist")) { 00619 anglelist = params["anglelist"]; 00620 nangles = anglelist.size() / 3; 00621 } else {*/ 00622 00623 if ( rotation == NULL ) throw NullPointerException("The transform object (required for projection), was not specified"); 00624 /* 00625 Dict p = t3d->get_rotation("spider"); 00626 00627 string angletype = "SPIDER"; 00628 float phi = p["phi"]; 00629 float theta = p["theta"]; 00630 float psi = p["psi"]; 00631 anglelist.push_back(phi); 00632 anglelist.push_back(theta); 00633 anglelist.push_back(psi); 00634 */ 00635 nangles = 1; 00636 //} 00637 00638 // initialize return object 00639 EMData* ret = new EMData(); 00640 ret->set_size(nx, ny, nangles); 00641 ret->to_zero(); 00642 //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; 00643 // loop over sets of angles 00644 for (int ia = 0; ia < nangles; ia++) { 00645 //int indx = 3*ia; 00646 //Dict d("type","spider","phi",anglelist[indx],"theta",anglelist[indx+1],"psi",anglelist[indx+2]); 00647 //Transform rotation(d); 00648 if (2*(ri+1)+1 > dim) { 00649 // Must check x and y boundaries 00650 for (int i = 0 ; i <= nn; i++) { 00651 int k = ipcube[i].loc[1] + origin[1]; 00652 Vec3f vb = ipcube[i].loc*(*rotation) + origin; 00653 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 00654 // check for pixels out-of-bounds 00655 int iox = int(vb[0]); 00656 if ((iox >= 0) && (iox < nx-1)) { 00657 int ioy = int(vb[1]); 00658 if ((ioy >= 0) && (ioy < ny-1)) { 00659 int ioz = int(vb[2]); 00660 if ((ioz >= 0) && (ioz < nz-1)) { 00661 // real work for pixels in bounds 00662 //cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<endl; 00663 //cout<<" TAKE"<<endl; 00664 float dx = vb[0] - iox; 00665 float dy = vb[1] - ioy; 00666 float dz = vb[2] - ioz; 00667 float a1 = (*image)(iox,ioy,ioz); 00668 float a2 = (*image)(iox+1,ioy,ioz) - a1; 00669 float a3 = (*image)(iox,ioy+1,ioz) - a1; 00670 float a4 = (*image)(iox,ioy,ioz+1) - a1; 00671 float a5 = -a2 -(*image)(iox,ioy+1,ioz) 00672 + (*image)(iox+1,ioy+1,ioz); 00673 float a61 = -(*image)(iox,ioy,ioz+1) 00674 + (*image)(iox+1,ioy,ioz+1); 00675 float a6 = -a2 + a61; 00676 float a7 = -a3 - (*image)(iox,ioy,ioz+1) 00677 + (*image)(iox,ioy+1,ioz+1); 00678 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1) 00679 + (*image)(iox+1,ioy+1,ioz+1); 00680 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx 00681 + (a7 + a8*dx)*dy) 00682 + a3*dy + dx*(a2 + a5*dy); 00683 } //else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLAATED Z "<<endl; } 00684 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED Y "<<endl; } 00685 }//else {cout<<" iox "<<iox<<" ioy "<<int(vb[1])<<" ioz "<<int(vb[2])<<" VIOLATED X "<<endl; } 00686 vb += rotation->get_matrix3_row(0); 00687 } 00688 } 00689 00690 } else { 00691 // No need to check x and y boundaries 00692 for (int i = 0 ; i <= nn; i++) { 00693 int k = ipcube[i].loc[1] + origin[1]; 00694 Vec3f vb = ipcube[i].loc*(*rotation) + origin; 00695 for (int j = ipcube[i].start; j <= ipcube[i].end; j++) { 00696 int iox = int(vb[0]); 00697 int ioy = int(vb[1]); 00698 int ioz = int(vb[2]); 00699 float dx = vb[0] - iox; 00700 float dy = vb[1] - ioy; 00701 float dz = vb[2] - ioz; 00702 float a1 = (*image)(iox,ioy,ioz); 00703 float a2 = (*image)(iox+1,ioy,ioz) - a1; 00704 float a3 = (*image)(iox,ioy+1,ioz) - a1; 00705 float a4 = (*image)(iox,ioy,ioz+1) - a1; 00706 float a5 = -a2 -(*image)(iox,ioy+1,ioz) 00707 + (*image)(iox+1,ioy+1,ioz); 00708 float a61 = -(*image)(iox,ioy,ioz+1) 00709 + (*image)(iox+1,ioy,ioz+1); 00710 float a6 = -a2 + a61; 00711 float a7 = -a3 - (*image)(iox,ioy,ioz+1) 00712 + (*image)(iox,ioy+1,ioz+1); 00713 float a8 = -a5 - a61 - (*image)(iox,ioy+1,ioz+1) 00714 + (*image)(iox+1,ioy+1,ioz+1); 00715 (*ret)(j,k,ia) += a1 + dz*(a4 + a6*dx 00716 + (a7 + a8*dx)*dy) 00717 + a3*dy + dx*(a2 + a5*dy); 00718 vb += rotation->get_matrix3_row(0); 00719 } 00720 } 00721 } 00722 } 00723 ret->update(); 00724 EMDeleteArray(ipcube); 00725 if(rotation) {delete rotation; rotation=0;} 00726 return ret; 00727 }
|
|
Definition at line 56 of file projector.cpp. |