EMAN2
Functions
emdata_transform.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

EMData * do_fft () const
 This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata.h", NEVER directly include this file. More...
 
void do_fft_inplace ()
 Do FFT inplace. More...
 
EMData * do_ift ()
 return the inverse fourier transform (IFT) image of the current image. More...
 
void do_ift_inplace ()
 
std::string render_amp8 (int x, int y, int xsize, int ysize, int bpl, float scale, int min_gray, int max_gray, float min_render, float max_render, float gamma, int flags)
 Render the image into an 8-bit image. More...
 
EMBytes render_ap24 (int x, int y, int xsize, int ysize, int bpl, float scale, int min_gray, int max_gray, float min_render, float max_render, float gamma, int flags)
 Render the image into an 8-bit image. More...
 
void render_amp24 (int x, int y, int xsize, int ysize, int bpl, float scale, int min_gray, int max_gray, float min_render, float max_render, void *ref, void cmap(void *, int coord, unsigned char *tri))
 Render the image into a 24-bit image. More...
 
void ri2ap ()
 convert the complex image from real/imaginary to amplitude/phase More...
 
void ap2ri ()
 convert the complex image from amplitude/phase to real/imaginary More...
 
void ri2inten ()
 convert the complex image from real/imaginary to Intensity/0. More...
 
EMData * bispecRotTransInvN (int N, int NK)
 This computes the rotational and translational bispectral invariants of an image. More...
 
EMData * bispecRotTransInvDirect (int type=0)
 This computes the rotational and translational bispectral invariants of an image. More...
 
void insert_clip (const EMData *const block, const IntPoint &origin)
 Insert a clip into this image. More...
 
void insert_scaled_sum (EMData *block, const FloatPoint &center, float scale=1.0, float mult_factor=1.0)
 Add a scaled image into another image at a specified location. More...
 

Function Documentation

◆ ap2ri()

void ap2ri ( )

convert the complex image from amplitude/phase to real/imaginary

Referenced by EMAN::EMData::apply_radial_func(), do_ift(), do_ift_inplace(), and ri2inten().

◆ bispecRotTransInvDirect()

EMData * EMData::bispecRotTransInvDirect ( int  type = 0)

This computes the rotational and translational bispectral invariants of an image.

the output is a single 3d Volume whose x,y labels are lengths, corresponding to the two lengths of sides of a triangle the z label is for the angle

Definition at line 1332 of file emdata_transform.cpp.

1333{
1334
1335 int EndP = this -> get_xsize(); // length(fTrueVec);
1336 int Mid = (int) ((1+EndP)/2);
1337 int End = 2*Mid-1;
1338
1339 int CountxyMax = End*End;
1340
1341// int *SortfkInds = new int[CountxyMax];
1342 int *kVecX = new int[CountxyMax];
1343 int *kVecY = new int[CountxyMax];
1344 float *fkVecR = new float[CountxyMax];
1345 float *fkVecI = new float[CountxyMax];
1346 float *absD1fkVec = new float[CountxyMax];
1347// float *absD1fkVecSorted = new float[CountxyMax];
1348
1349
1350 float *jxjyatan2 = new float[End*End];
1351
1352
1353 EMData * ThisCopy = new EMData(End,End);
1354
1355 for (int jx=0; jx <End ; jx++) { // create jxjyatan2
1356 for (int jy=0; jy <End ; jy++) {
1357 float ValNow = this -> get_value_at(jx,jy);
1358 ThisCopy -> set_value_at(jx,jy,ValNow);
1359 jxjyatan2[jy*End + jx] = atan2((float)(jy+1-Mid) , (float)(jx +1 -Mid));
1360// cout<< " jxM= " << jx+1<<" jyM= " << jy+1<< "ValNow" << ValNow << endl; // Works
1361 }}
1362
1363
1364 EMData* fk = ThisCopy -> do_fft();
1365 fk ->process_inplace("xform.fourierorigin.tocenter");
1366
1367// Create kVecX , kVecy etc
1368
1369 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) { // kEx twice the value of the Fourier
1370 // x variable: EMAN index for real, imag
1371 int kx = kEx/2; // kx is the value of the Fourier variable
1372 int kIx = kx+Mid-1; // This is the value of the index for a matlab image (-1)
1373 int kCx = -kx ;
1374 int kCIx = kCx+ Mid-1 ;
1375 for (int kEy= 0 ; kEy<End; kEy++) { // This is the value of the EMAN index
1376 int kIy = kEy ; // This is the value of the index for a matlab image (-1)
1377 int ky = kEy+1-Mid; // (kEy+ Mid-1)%End - Mid+1 ; // This is the actual value of the Fourier variable
1378 float realVal = fk -> get_value_at(kEx ,kEy) ;
1379 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
1380 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
1381 float fkAng = atan2(imagVal,realVal);
1382
1383 float NewRealVal ;
1384 float NewImagVal ;
1385 float AngMatlab ;
1386
1387 if (kIx==Mid-1) {
1388// AngMatlab = -fkAng - 2.*M_PI*(kIy+ 1-Mid)*(Mid)/End;
1389 }
1390
1391 if (kIx>Mid-1){
1392// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << AngMatlab << endl;
1393 }
1394
1395 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
1396 NewRealVal = absVal*cos(AngMatlab);
1397 NewImagVal = absVal*sin(AngMatlab);
1398
1399
1400 fkVecR[ kIy +kIx *End] = NewRealVal ;
1401 fkVecR[(End-1-kIy)+kCIx*End] = NewRealVal ;
1402 fkVecI[ kIy +kIx *End] = NewImagVal ;
1403 fkVecI[(End-1-kIy)+kCIx*End] = -NewImagVal ;
1404 absD1fkVec[(End-1-kIy) + kIx *End] = absVal;
1405 absD1fkVec[(End-1-kIy) + kCIx *End] = absVal;
1406 kVecX[kIy+kIx *End] = kx ;
1407 kVecX[kIy+kCIx *End] = kCx ;
1408 kVecY[kIy+kIx *End] = ky ;
1409 kVecY[kIy+kCIx *End] = ky ;
1410
1411 // cout << " kIxM= " << kIx+1 << " kIy=" << kIy+1 << " fkVecR[i] =" << NewRealVal << " fkVecI[i]=" << NewImagVal <<" angle[i]= " << AngMatlab << " Total Index" << kIy+kIx *End << endl;
1412
1413// printf("kx=%d,ky=%d,tempVal =%f+ i %4.2f \n",kx,ky,realVal,imagVal );
1414// cout << "kx = " << kx << "; ky = "<< ky << "; val is" << realVal<<"+ i "<<imagVal<< endl;
1415
1416// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; fkAng*9/ 2pi is " << fkAng*9/2/M_PI<< endl;
1417// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; absval is " << absVal<< "; realval is " << NewRealVal<< "; imagval is " << NewImagVal<< endl;
1418 }
1419 }
1420
1421
1422// for (int TotalInd = 0 ; TotalInd < CountxyMax ; TotalInd++){
1423// int kx = kVecX[TotalInd]; // This is the value of the index for a matlab image (-1)
1424// int kIx = kx+Mid-1; // This is the value of the index for a matlab image (-1)
1425// int ky = kVecY[TotalInd];
1426// int kIy = ky+Mid-1; // This is the value of the index for a matlab image (-1)
1427 //float fkR = fkVecR[kIy+kIx *End] ;
1428 //float fkI = fkVecI[kIy+kIx *End] ;
1429// }
1430
1431 float frR= 3.0/4.0;
1432 frR= 1;
1433 int LradRange= (int) (1+floor(Mid/frR -.1)) ;
1434
1435 float *radRange = new float[LradRange]; //= 0:.75:(Mid-1);
1436 for (int irad=0; irad < LradRange; irad++){
1437 radRange[irad] = frR*irad;
1438// cout << " irad = " << irad << " radRange[irad]= " << radRange[irad] << " LradRange= " << LradRange << endl;
1439 }
1440
1441 cout << "Starting the calculation of invariants" << endl;
1442
1443
1444 if (type==0) {
1445 int LthetaRange = 59;
1446 float ftR = (2.0f*M_PI/LthetaRange );
1447 float *thetaRange = new float[LthetaRange]; //= 0:.75:(Mid-1);
1448
1449 for (int ith=0; ith < LthetaRange; ith++){
1450 thetaRange[ith] = ftR*ith; }
1451
1452 int TotalVol = LradRange*LradRange*LthetaRange;
1453
1454 float *RotTransInv = new float[TotalVol];
1455 float *WeightInv = new float[TotalVol];
1456
1457 for (int jW=0; jW<TotalVol; jW++) {
1458 RotTransInv[jW] = 0;
1459 WeightInv[jW] = 0;
1460 }
1461
1462 for (int jW=0; jW<TotalVol; jW++) {
1463 RotTransInv[jW] = 0;
1464 WeightInv[jW] = 0;
1465 }
1466 // float *RotTransInv = new float[LradRange*LradRange ] ;
1467 // float *RotTransInvN = new float[LradRange*LradRange*(NMax+1) ] ;
1468
1469 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){ // Main Section for type 0
1470 int kx = kVecX[Countkxy] ;
1471 int ky = kVecY[Countkxy] ;
1472 float k2 = ::sqrt((float)(kx*kx+ky*ky));
1473 float phiK =0;
1474 if (k2>0) phiK = jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; // phiK=atan2(ky,kx);
1475 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
1476 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
1477 // printf("Countkxy=%d,\t kx=%d, ky=%d, fkR=%3.2f,fkI=%3.2f \n", Countkxy, kx, ky, fkR, fkI);
1478
1479 if ((k2==0)|| (k2>Mid) ) { continue;}
1480
1481 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){ // This is the innermost loop
1482 int qx = kVecX[Countqxy] ;
1483 int qy = kVecY[Countqxy] ;
1484 float q2 = ::sqrt((float)(qx*qx+qy*qy));
1485 if ((q2==0)|| (q2>Mid) ) {continue;}
1486 float phiQ =0;
1487 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; // phiQ=atan2(qy,qx);
1488 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
1489 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
1490 int kCx = (-kx-qx);
1491 int kCy = (-ky-qy);
1492 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
1493 int kCIy = ((kCy+Mid+2*End)%End);
1494 kCx = ((kCIx+End-1)%End)+1-Mid; // correct
1495 kCy = ((kCIy+End-1)%End)+1-Mid ; // correct
1496
1497// float C2 = ::sqrt((float)(kCx*kCx+ kCy*kCy));
1498 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
1499 float fCR = fkVecR[CountCxy];
1500 float fCI = fkVecI[CountCxy];
1501 /* if (Countkxy==1) {
1502 printf(" Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", Countqxy, absD1fkVec[Countqxy],qx, qy);
1503 printf(" CountCxy=%d, absD1fkVec[CountCxy]=%f,kCx=%d,kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
1504 }*/
1505// float phiC = jxjyatan2[ (kCy+Mid-1)*End + kCx+Mid-1];
1506 float phiQK = (4*M_PI+phiQ-phiK);
1507 while (phiQK> (2*M_PI)) phiQK -= (2*M_PI);
1508
1509
1510
1511 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
1512
1513 if ((q2<k2) ) continue;
1514// if ((q2<k2) || (C2<k2) || (C2<q2)) continue;
1515
1516 // printf(" CountCxy=%d, absD1fkVec[CountCxy]=%f,kCx=%d,kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
1517
1518 // up to here, matched perfectly with Matlab
1519
1520 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
1521 int k2IndHi = k2IndLo;
1522 float k2Lo= radRange[k2IndLo];
1523 if (k2IndLo+1< LradRange) {
1524 k2IndHi = k2IndLo+1;
1525 }
1526// float k2Hi= radRange[k2IndHi];
1527
1528 float kCof =k2-k2Lo;
1529
1530 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
1531 int q2IndHi=q2IndLo;
1532 float q2Lo= radRange[q2IndLo];
1533 if (q2IndLo+1 < LradRange) {
1534 q2IndHi = q2IndLo+1 ;
1535 }
1536 float qCof = q2-q2Lo;
1537
1538 if ((qCof<0) || (qCof >1) ) {
1539 cout<< "Weird! qCof="<< qCof << " q2="<< q2 << " q2IndLo="<< q2IndLo << endl ;
1540 int x ;
1541 cin >> x ;
1542 }
1543
1544 int thetaIndLo = 0; while ((phiQK>=thetaRange[thetaIndLo+1])&& (thetaIndLo+1<LthetaRange)) thetaIndLo +=1;
1545 int thetaIndHi = thetaIndLo;
1546
1547 float thetaLo = thetaRange[thetaIndLo];
1548 float thetaHi = thetaLo;
1549 float thetaCof = 0;
1550
1551 if (thetaIndLo+1< LthetaRange) {
1552 thetaIndHi = thetaIndLo +1;
1553 }else{
1554 thetaIndHi=0;
1555 }
1556
1557 thetaHi = thetaRange[thetaIndHi];
1558
1559 if (thetaHi==thetaLo) {
1560 thetaCof =0 ;
1561 } else {
1562 thetaCof = (phiQK-thetaLo)/(thetaHi-thetaLo);
1563 }
1564
1565 if ((thetaCof>2*M_PI) ) {
1566 cout<< "Weird! thetaCof="<< thetaCof <<endl ;
1567 thetaCof=0;
1568 }
1569
1570
1571 // if ((thetaIndLo>=58) || (k2IndLo >= LradRange-1) || (q2IndLo >= LradRange-1) ) {
1572
1573
1574 for (int jk =1; jk<=2; jk++){
1575 for (int jq =1; jq<=2; jq++){
1576 for (int jtheta =1; jtheta<=2; jtheta++){
1577
1578 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1))
1579 * (thetaCof+(1-2*thetaCof)*(jtheta==1));
1580
1581
1582 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
1583 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
1584 int thetaInd = thetaIndLo*(jtheta==1) + thetaIndHi*(jtheta ==2);
1585 int TotalInd = thetaInd*LradRange*LradRange+q2Ind*LradRange+k2Ind;
1586 /* if (TotalInd+1 >= LthetaRange*LradRange*LradRange) {
1587 cout << "Weird!!! TotalInd="<< TotalInd << " IndMax" << LthetaRange*LradRange*LradRange << " LradRange=" << LradRange << endl;
1588 cout << "k2Ind= "<< k2Ind << " q2Ind="<< q2Ind << " thetaInd="<< thetaInd << " q2IndLo="<< q2IndLo << " q2IndHi="<< q2IndHi << endl;
1589 cout << "k2=" << k2 << "q2=" << q2 << " phiQK=" << phiQK*180.0/M_PI<< endl;
1590 }*/
1591
1592 RotTransInv[TotalInd] += Weight*bispectemp;
1593 WeightInv[TotalInd] += Weight;
1594 // cout << "k2Ind= "<< k2Ind << " q2Ind="<< q2Ind << "Weight=" << Weight << endl;
1595 }}}
1596 } // Countqxy
1597 } // Countkxy
1598
1599 cout << "Finished Main Section " << endl;
1600
1601 /* RotTransInvN[jr1 + LradRange*jr2+LradRange*LradRange*N] = RotTransInvTemp ;*/
1602
1603 cout << " LradRange " <<LradRange <<" LthetaRange " << LthetaRange << endl;
1604 EMData *RotTransInvF = new EMData(LradRange,LradRange,LthetaRange);
1605 EMData *WeightImage = new EMData(LradRange,LradRange,LthetaRange);
1606
1607 // cout << "FFFFFFF" << endl;
1608 //
1609 // RotTransInvF -> set_size(LradRange,LradRange,LthetaRange);
1610 //
1611 // cout << "GGGG" << endl;
1612
1613 for (int jtheta =0; jtheta < LthetaRange; jtheta++){ // write out to RotTransInvF
1614 for (int jq =0; jq<LradRange; jq++){ // LradRange
1615 for (int jk =0; jk<LradRange ; jk++){// LradRange
1616 // cout << "Hi There" << endl;
1617 int TotalInd = jtheta*LradRange*LradRange+jq*LradRange+jk;
1618 float Weight = WeightInv[TotalInd];
1619 WeightImage -> set_value_at(jk,jq,jtheta,Weight);
1620 RotTransInvF -> set_value_at(jk,jq,jtheta,0);
1621 if (Weight <= 0) continue;
1622 RotTransInvF -> set_value_at(jk,jq,jtheta,RotTransInv[TotalInd] / Weight);// include /Weight
1623 int newjtheta = (LthetaRange- jtheta)%LthetaRange;
1624 RotTransInvF -> set_value_at(jq,jk,newjtheta,RotTransInv[TotalInd]/Weight );// include /Weight
1625 }
1626 }
1627 }
1628
1629 cout << " Almost Done " << endl;
1630#ifdef _WIN32
1631 _unlink("WeightImage.???");
1632#else
1633 int rslt;
1634 rslt = system("rm -f WeightImage.???"); rslt++;
1635#endif //_WIN32
1636 WeightImage -> write_image("WeightImage.img");
1637
1638 return RotTransInvF ;
1639 } // Finish type 0
1640
1641 if (type==1) {
1642 int TotalVol = LradRange*LradRange;
1643
1644 float *RotTransInv = new float[TotalVol];
1645 float *WeightInv = new float[TotalVol];
1646
1647 for (int jW=0; jW<TotalVol; jW++) {
1648 RotTransInv[jW] = 0;
1649 WeightInv[jW] = 0;
1650 }
1651
1652
1653 for (int Countkxy =0; Countkxy<CountxyMax; Countkxy++){
1654 int kx = kVecX[Countkxy] ;
1655 int ky = kVecY[Countkxy] ;
1656 float k2 = ::sqrt((float)(kx*kx+ky*ky));
1657 float fkR = fkVecR[(ky+Mid-1) + (kx+Mid-1) *End] ;
1658 float fkI = fkVecI[(ky+Mid-1) + (kx+Mid-1) *End] ;
1659 // printf("Countkxy=%d,\t kx=%d, ky=%d, fkR=%3.2f,fkI=%3.2f \n", Countkxy, kx, ky, fkR, fkI);
1660
1661 if ((k2==0)|| (k2>Mid) ) { continue;}
1662
1663 for (int Countqxy =0; Countqxy<CountxyMax; Countqxy++){ // This is the innermost loop
1664
1665// up to here, matched perfectly with Matlab
1666 int qx = kVecX[Countqxy] ;
1667 int qy = kVecY[Countqxy] ;
1668 float q2 = ::sqrt((float)(qx*qx+qy*qy));
1669 if ((q2==0)|| (q2>Mid) ) {continue;}
1670 if ((q2<k2) ) continue;
1671
1672 float fqR = fkVecR[(qy+Mid-1) + (qx+Mid-1) *End] ;
1673 float fqI = fkVecI[(qy+Mid-1) + (qx+Mid-1) *End] ;
1674
1675 int kCx = (-kx-qx);
1676 int kCy = (-ky-qy);
1677 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
1678 int kCIy = ((kCy+Mid+2*End)%End);
1679 kCx = ((kCIx+End-1)%End)+1-Mid; // correct
1680 kCy = ((kCIy+End-1)%End)+1-Mid ; // correct
1681
1682// float C2 = ::sqrt((float)(kCx*kCx+ kCy*kCy));
1683 int CountCxy = (kCx+Mid-1)*End+(kCy+Mid-1);
1684 float fCR = fkVecR[CountCxy];
1685 float fCI = fkVecI[CountCxy];
1686
1687
1688 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI));
1689
1690
1691 int k2IndLo = 0; while ((k2>=radRange[k2IndLo+1]) && (k2IndLo+1 < LradRange ) ) k2IndLo +=1;
1692 int k2IndHi = k2IndLo;
1693 float k2Lo= radRange[k2IndLo];
1694 if (k2IndLo+1< LradRange) {
1695 k2IndHi = k2IndLo+1;
1696 }
1697// float k2Hi= radRange[k2IndHi];
1698
1699 float kCof =k2-k2Lo;
1700
1701 int q2IndLo = 0; while ((q2>=radRange[q2IndLo+1]) && (q2IndLo+1 < LradRange ) ) q2IndLo +=1;
1702 int q2IndHi=q2IndLo;
1703 float q2Lo= radRange[q2IndLo];
1704 if (q2IndLo+1 < LradRange) {
1705 q2IndHi = q2IndLo+1 ;
1706 }
1707 float qCof = q2-q2Lo;
1708
1709
1710 for (int jk =1; jk<=2; jk++){
1711 for (int jq =1; jq<=2; jq++){
1712
1713 float Weight = (kCof+(1-2*kCof)*(jk==1))*(qCof+(1-2*qCof)*(jq==1));
1714
1715 int k2Ind = k2IndLo*(jk==1) + k2IndHi*(jk==2);
1716 int q2Ind = q2IndLo*(jq==1) + q2IndHi*(jq==2);
1717 int TotalInd = q2Ind*LradRange+k2Ind;
1718 RotTransInv[TotalInd] += Weight*bispectemp;
1719 WeightInv[TotalInd] += Weight;
1720 // cout << "k2Ind= "<< k2Ind << " q2Ind="<< q2Ind << "Weight=" << Weight << endl;
1721 }}
1722 } // Countqxy
1723 } // Countkxy
1724
1725// cout << "Finished Main Section " << endl;
1726// cout << " LradRange " <<LradRange << endl;
1727
1728
1729 EMData *RotTransInvF = new EMData(LradRange,LradRange);
1730 EMData *WeightImage = new EMData(LradRange,LradRange);
1731
1732 for (int jk =0; jk<LradRange ; jk++){// LradRange
1733 for (int jq =jk; jq<LradRange; jq++){ // LradRange
1734 int TotalInd = jq*LradRange+jk;
1735 int TotalIndBar = jq*LradRange+jk;
1736 float Weight = WeightInv[TotalInd] + WeightInv[TotalIndBar];
1737 if (Weight <=0) continue;
1738 WeightImage -> set_value_at(jk,jq,Weight);
1739 WeightImage -> set_value_at(jq,jk,Weight);
1740 float ValNow = cbrt( (RotTransInv[TotalInd] + RotTransInv[TotalIndBar]) / Weight ) ;
1741 RotTransInvF -> set_value_at(jk,jq,ValNow);// include /Weight
1742 RotTransInvF -> set_value_at(jq,jk,ValNow );// include /Weight
1743 }}
1744
1745#ifdef _WIN32
1746 _unlink("WeightImage.???");
1747#else
1748 int rslt;
1749 rslt = system("rm -f WeightImage.???"); rslt++;
1750#endif //_WIN32
1751 WeightImage -> write_image("WeightImage.img");
1752
1753 return RotTransInvF ;
1754 }
1755 return 0;
1756}
EMAN::EMData
EMData stores an image's data and defines core image processing routines.
Definition: emdata.h:82
set_value_at
void set_value_at(Vec3i loc, float val)
set_value_at with Vec3i
Definition: emdata_core.h:552
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
get_value_at
float get_value_at(int x, int y, int z) const
Get the pixel density value at coordinates (x,y,z).
Definition: emdata_core.h:221
write_image
void write_image(const string &filename, int img_index=0, EMUtil::ImageType imgtype=EMUtil::IMAGE_UNKNOWN, bool header_only=false, const Region *region=0, EMUtil::EMDataType filestoragetype=EMUtil::EM_FLOAT, bool use_host_endian=true)
write the header and data out to an image.
get_xsize
int get_xsize() const
Get the image x-dimensional size.
Definition: emdata_metadata.h:544
do_fft
EMData * do_fft() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
Definition: emdata_transform.cpp:55
x
#define x(i)
Definition: projector.cpp:1517

References do_fft(), EMAN::EMData::EMData(), get_value_at(), get_xsize(), set_value_at(), sqrt(), write_image(), and x.

◆ bispecRotTransInvN()

EMData * EMData::bispecRotTransInvN ( int  N,
int  NK 
)

This computes the rotational and translational bispectral invariants of an image.

The invariants are labelled by the Fourier Harmonic label given by N. fVec is the real input image NK is the number of Fourier components one wishes to use in calculating this bispectrum the output is a single 2D image whose x,y labels are lengths, corresponding to the two lengths of sides of a triangle

Definition at line 1034 of file emdata_transform.cpp.

1035{
1036
1037 int EndP = this -> get_xsize(); // length(fTrueVec);
1038 int Mid = (int) ((1+EndP)/2);
1039 int End = 2*Mid-1;
1040
1041 int CountxyMax = End*End;
1042
1043 int *SortfkInds = new int[CountxyMax];
1044 int *kVecX = new int[CountxyMax];
1045 int *kVecY = new int[CountxyMax];
1046 float *fkVecR = new float[CountxyMax];
1047 float *fkVecI = new float[CountxyMax];
1048 float *absD1fkVec = new float[CountxyMax];
1049 float *absD1fkVecSorted = new float[CountxyMax];
1050
1051 float *jxjyatan2 = new float[End*End]; // jxjyatan2[jy*End + jx] = atan2(jy+1-Mid , jx +1 -Mid);
1052
1053 EMData * ThisCopy = new EMData(End,End);
1054
1055 for (int jx=0; jx <End ; jx++) {
1056 for (int jy=0; jy <End ; jy++) {
1057 float ValNow = this -> get_value_at(jx,jy);
1058 ThisCopy -> set_value_at(jx,jy,ValNow);
1059// cout<< " jxM= " << jx+1<<" jyM= " << jy+1<< "ValNow" << ValNow << endl; // Works
1060 }}
1061
1062
1063 EMData* fk = ThisCopy -> do_fft();
1064 fk ->process_inplace("xform.fourierorigin.tocenter");
1065
1066// EMData* fk
1067 EMData* fkRCopy = new EMData(End,End);
1068 EMData* fkICopy = new EMData(End,End);
1069 EMData* fkCopy = new EMData(End,End);
1070
1071
1072 for (int jCount= 0; jCount<End*End; jCount++) {
1073// jCount = jy*End + jx;
1074 int jx = jCount%End ;
1075 int jy = (jCount-jx)/End ;
1076 jxjyatan2[jCount] = atan2((float)(jy+1-Mid) , (float)(jx +1-Mid));
1077 }
1078
1079
1080 for (int kEx= 0; kEx<2*Mid; kEx=kEx+2) { // kEx twice the value of the Fourier
1081 // x variable: EMAN index for real, imag
1082 int kx = kEx/2; // kx is the value of the Fourier variable
1083 int kIx = kx+Mid-1; // This is the value of the index for a matlab image (-1)
1084 int kCx = -kx ;
1085 int kCIx = kCx+ Mid-1 ;
1086 for (int kEy= 0 ; kEy<End; kEy++) { // This is the value of the EMAN index
1087 int kIy = kEy ; // This is the value of the index for a matlab image (-1)
1088 int ky = kEy+1-Mid; // (kEy+ Mid-1)%End - Mid+1 ; // This is the actual value of the Fourier variable
1089 float realVal = fk -> get_value_at(kEx ,kEy) ;
1090 float imagVal = fk -> get_value_at(kEx+1,kEy) ;
1091 float absVal = ::sqrt(realVal*realVal+imagVal*imagVal);
1092 float fkAng = atan2(imagVal,realVal);
1093
1094 float NewRealVal ;
1095 float NewImagVal ;
1096 float AngMatlab ;
1097
1098 if (kIx==Mid-1) {
1099// AngMatlab = -fkAng - 2.*M_PI*(kIy+ 1-Mid)*(Mid)/End;
1100// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << AngMatlab << endl;
1101 }
1102
1103 if (kIx>Mid-1){
1104// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << AngMatlab << endl;
1105 }
1106
1107 AngMatlab = fkAng - 2.0f*M_PI*(kx +ky)*(Mid)/End;
1108 NewRealVal = absVal*cos(AngMatlab);
1109 NewImagVal = absVal*sin(AngMatlab);
1110
1111
1112 fkVecR[kIy+kIx *End] = NewRealVal ;
1113 fkVecR[kIy+kCIx*End] = NewRealVal ;
1114 fkVecI[kIy+kIx *End] = NewImagVal ;
1115 fkVecI[kIy+kCIx*End] = -NewImagVal ;
1116 absD1fkVec[kIy + kIx *End] = absVal;
1117 absD1fkVec[kIy + kCIx *End] = absVal;
1118 kVecX[kIy+kIx *End] = kx ;
1119 kVecX[kIy+kCIx *End] = kCx ;
1120 kVecY[kIy+kIx *End] = ky ;
1121 kVecY[kIy+kCIx *End] = ky ;
1122// printf("kx=%d,ky=%d,tempVal =%f+ i %4.2f \n",kx,ky,realVal,imagVal );
1123// cout << "kx = " << kx << "; ky = "<< ky << "; val is" << realVal<<"+ i "<<imagVal<< endl;
1124
1125// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; fkAng*9/ 2pi is " << fkAng*9/2/M_PI<< endl;
1126// cout << "kIMx = "<< kIx+1 << "; kIMy = "<< kIy+1 <<"; absval is " << absVal<< "; realval is " << NewRealVal<< "; imagval is " << NewImagVal<< endl;
1127 fkCopy -> set_value_at(kIx ,kIy, absVal);
1128 fkCopy -> set_value_at(kCIx,kIy, absVal);
1129 fkRCopy -> set_value_at(kIx, kIy, NewRealVal);
1130 fkRCopy -> set_value_at(kCIx,kIy, NewRealVal);
1131 fkICopy -> set_value_at(kIx, kIy, NewImagVal);
1132 fkICopy -> set_value_at(kCIx,kIy,-NewImagVal);
1133
1134 }
1135 }
1136#ifdef _WIN32
1137 _unlink("fkCopy.???");
1138 _unlink("fk?Copy.???");
1139#else
1140 int rslt;
1141 rslt = system("rm -f fkCopy.???"); rslt++;
1142 rslt = system("rm -f fk?Copy.???"); rslt++;
1143#endif //_WIN32
1144 fkCopy -> write_image("fkCopy.img");
1145 fkRCopy -> write_image("fkRCopy.img");
1146 fkICopy -> write_image("fkICopy.img");
1147
1148 cout << "Starting the sort "<< endl;
1149
1150 vector< pair<float, int> > absInds;
1151 for(int i = 0; i < CountxyMax; ++i ) {
1152 pair<float,int> p;
1153 p = make_pair(absD1fkVec[i],i); // p = make_pair(rand(),i);
1154 absInds.push_back( p);
1155 }
1156
1157 std::sort(absInds.begin(),absInds.end());
1158
1159 for(int i = 0; i < CountxyMax; ++i ) {
1160 pair<float,int> p ;
1161 p = absInds[i] ;
1162 absD1fkVecSorted[CountxyMax-1-i] = p.first ;
1163 SortfkInds[CountxyMax-1-i] = p.second ;
1164 }
1165
1166 cout << "Ending the sort "<< endl;
1167
1168// float AngsMat[] ={2.8448, -0.3677,-0.2801,-1.0494,-1.7836,-2.5179, 2.9959, 3.0835,-0.1290,-0.8876,2.1829, 2.2705,1.5011,0.7669,0.0327,-0.7366,-0.6489,2.4215,-1.6029,1.4676,1.5552,0.7859,0.0517,-0.6825,-1.4518,-1.3642,1.7063,-1.7845,1.2859,1.3736,0.6043,-0.1299,-0.8642,-1.6335,-1.5459,1.5247,-1.6546,1.4159,1.5036,0.7342,0,-0.7342,-1.5036,-1.4159,1.6546,-1.5247,1.5459,1.6335,0.8642,0.1299,-0.6043,-1.3736,-1.286,1.7846,-1.7063,1.3642,1.4519,0.6825,-0.0517,-0.7859,-1.5553,-1.4676,1.6029,-2.4216,0.649,0.7366,-0.0327,-0.767,-1.5012,-2.2705,-2.1829,0.8877,0.1291,-3.0836,-2.9959,2.5179,1.7837,1.0495,0.2801,0.3677,-2.8449};
1169
1170
1171 for(int i = 0; i < CountxyMax; ++i ) { // creates a new fkVec
1172// int Si = SortfkInds[i];
1173// int kIx = (int) Si/End; kIx = (int) i/End; // i = kIx*End+kIy
1174// int kIy = Si - kIx*End; kIy = i - kIx*End;
1175// int iC = (End-1-kIx)*End + (End-1-kIy);
1176// if (i<30) { cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " valAft=" << absD1fkVecSorted[i]<< " valBef=" << absD1fkVec[Si] << " SortfkInds = " << Si <<endl; }// This worked
1177// cout<< "i= " << i << " kIx= " << kIx << " kIy=" << kIy << " fkVecR[i] =" << fkVecR[i]<< " fkVecI[i]=" << fkVecI[i] <<" angle[i]= " << fkAng << endl;
1178 }
1179 cout<< "Ratio of Last Amplitude to First Amplitude= " << absD1fkVecSorted[NK] /absD1fkVecSorted[0] << endl;
1180
1181// pause;
1182
1183// for(int i = 0; i < NK; ++i ) { // Prints out the new fkVec , CountxyMax
1184// int Si= SortfkInds[i];
1185// int kIx = (int) Si/End; // i = kIx*End+kIy
1186// int kIy = Si - kIx*End;
1187// cout << " kIxM= " << kIx+1 << " kIyM=" << kIy+1 << " fkVecAbs=" << ::sqrt(fkVecR[Si]*fkVecR[Si] + fkVecI[Si]* fkVecI[Si]) << " fkVecAbs=" << absD1fkVecSorted[i] << " kx= " << kVecX[Si] << " ky=" << kVecY[Si] << endl;
1188// }
1189
1190// angEMAN+angMat+angDiff =0 mod 2 pi
1191
1192// angDiff= 2*pi*((-4):4)*(Mid)/End; angEMAN+angMat+angDiff= integer*2 *pi
1193// [ absD1fkVecSorted, SortfkInds] =sort( absD1fkVec,'descend') ;
1194// Util::sort_mat(&absD1fkVec[0],&absD1fkVec[Countxy],&SortfkInds[0],&SortfkInds[Countxy]);
1195
1196
1197// Let radial sampling be 0:0.5:(Mid-1)
1198
1199 // int NK= min(12,CountxyMax) ;
1200
1201
1202
1203 cout << "NK = " << NK << endl;
1204 float frR= 3.0/4.0;
1205 int LradRange= (int) (floor(Mid/frR)) ;
1206
1207 float *radRange = new float[LradRange]; //= 0:.75:(Mid-1);
1208 radRange[0]=0;
1209 for (int irad=1; irad < LradRange; irad++){
1210 radRange[irad] = radRange[irad-1] + frR; }
1211
1212
1213
1214 // should equal to (2*Mid-1)
1215 cout << "Starting the calculation of invariants for N= " << N << endl;
1216
1217/* int NMax=5; */
1218
1219 EMData* RotTransInv = new EMData();
1220 RotTransInv -> set_size(LradRange,LradRange);
1221
1222
1223// float *RotTransInv = new float[LradRange*LradRange ] ;
1224// float *RotTransInvN = new float[LradRange*LradRange*(NMax+1) ] ;
1225
1226// for (int N=0 ; N<NMax; N++) {
1227
1228 for (int jr1=0; jr1 < LradRange ; jr1++ ) { // LradRange
1229 float r1= radRange[jr1];
1230// cout << "Pre jr2 "<< endl;
1231 for (int jr2=0; jr2<LradRange; jr2++ ) { //LradRange
1232 float r2= radRange[jr2];
1233 float RotTransInvTemp=0;
1234 for (int jCountkxy =0; jCountkxy<NK; jCountkxy++){
1235 int Countkxy =SortfkInds[jCountkxy] ; // 1: CountxyMax
1236 int kx = kVecX[Countkxy] ;
1237 int ky = kVecY[Countkxy] ;
1238 float k2 = (float)(kx*kx+ky*ky);
1239 if (k2==0) { continue;}
1240 float phiK =0;
1241 if (k2>0) phiK= jxjyatan2[ (ky+Mid-1)*End + kx+Mid-1]; phiK= atan2((float)ky,(float)kx);
1242
1243 float fkR = fkVecR[Countkxy] ;
1244 float fkI = fkVecI[Countkxy] ;
1245/* printf("jCountkxy=%d, Countkxy=%d,absD1fkVec(Countkxy)=%f,\t\t kx=%d, ky=%d \n", jCountkxy, Countkxy, absD1fkVec[Countkxy], kx, ky);*/
1246
1247 for (int jCountqxy =0; jCountqxy<NK; jCountqxy++){
1248 int Countqxy =SortfkInds[jCountqxy] ; // Countqxy is the index for absD1fkVec
1249 int qx = kVecX[Countqxy] ;
1250 int qy = kVecY[Countqxy] ;
1251 int q2 = qx*qx+qy*qy;
1252 if (q2==0) {continue;}
1253 float phiQ =0;
1254 if (q2>0) phiQ = jxjyatan2[ (qy+Mid-1)*End + qx+Mid-1]; phiQ=atan2((float)qy,(float)qx);
1255 float fqR = fkVecR[Countqxy] ;
1256 float fqI = fkVecI[Countqxy] ;
1257 int kCx = (-kx-qx);
1258 int kCy = (-ky-qy);
1259 int kCIx = ((kCx+Mid+2*End)%End);// labels of the image in C
1260 int kCIy = ((kCy+Mid+2*End)%End);
1261 kCx = kCIx-Mid; // correct
1262 kCy = kCIy-Mid; // correct
1263 int CountCxy = kCIx*End+kCIy;
1264 float fCR = fkVecR[CountCxy];
1265 float fCI = fkVecI[CountCxy];
1266 if (jr1+jr2==-1) {
1267 printf("jCountqxy=%d , Countqxy=%d, absD1fkVec(Countqxy)=%f,qx=%d, qy=%d \n", jCountqxy, Countqxy, absD1fkVec[Countqxy],qx, qy);
1268 printf(" CountCxy=%d,absD1fkVec[CountCxy]=%f, kCx=%d, kCy=%d \n",CountCxy, absD1fkVec[CountCxy], kCx, kCy );
1269 }
1270 for (int p=0; p<NK; p++){
1271// printf("p=%d, SortfkInds[p]=%d, CountCxy =%d \n", p,SortfkInds[p], CountCxy);
1272 if (SortfkInds[p]==CountCxy){
1273 float Arg1 = 2.0f*M_PI*r1*::sqrt((float) q2)/End;
1274 float Arg2 = 2.0f*M_PI*r2*::sqrt((float) k2)/End;
1275// printf("Arg1=%4.2f, Arg2=%4.2f, \n",Arg1, Arg2 );
1276// if (Arg1+ Arg2<15) {
1277 float bispectemp = (fkR*(fqR*fCR -fqI*fCI) -fkI*(fqI*fCR +fqR*fCI))
1278 * cos(N*(phiK-phiQ+M_PI));
1279 bispectemp -= (fkR*(fqR*fCI + fqI*fCR) +fkI*(fqR*fCR - fqI*fCI))
1280 * sin(N*(phiK-phiQ+M_PI));
1281 float bess1 = calc_bessel(N, Arg1 );
1282 float bess2 = calc_bessel(N, Arg2 );
1283// printf("fkr=%4.2f, fqr=%4.2f, bess1=%4.2f,bess2=%4.2f \n",fkR, fqR, bess1, bess2);
1284/* printf("p =%d, SortfkInds[p]=%d, CountCxy=%d, Arg1 =%4.2f, bess1=%4.2f, \n",
1285 p, SortfkInds[p],CountCxy, Arg1, bess1);*/
1286 RotTransInvTemp = RotTransInvTemp + bispectemp * bess1*bess2 ;
1287// }
1288 }
1289 }
1290 } // jCountqxy
1291 } // jCountkxy
1292 RotTransInv -> set_value_at(jr1,jr2, RotTransInvTemp) ;
1293/* RotTransInvN[jr1 + LradRange*jr2+LradRange*LradRange*N] = RotTransInvTemp ;*/
1294 } //jr2
1295 } //jr1
1296// }//N
1297
1298 return RotTransInv ;
1299
1300
1301}
set_size
void set_size(int nx, int ny=1, int nz=1, bool noalloc=false)
Resize this EMData's main board memory pointer.
calc_bessel
float calc_bessel(const int n, const float &x)
Definition: emdata_transform.cpp:1027

References calc_bessel(), do_fft(), EMAN::EMData::EMData(), get_value_at(), get_xsize(), set_size(), set_value_at(), sqrt(), and write_image().

◆ do_fft()

EMData * EMData::do_fft ( ) const

This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata.h", NEVER directly include this file.

return the fast fourier transform (FFT) image of the current image. the current image is not changed. The result is in real/imaginary format.

Returns
The FFT of the current image in real/imaginary format.

Definition at line 55 of file emdata_transform.cpp.

56{
57 ENTERFUNC;
58#ifdef FFT_CACHING
59 if (fftcache!=0) {
60 return fftcache->copy();
61 }
62#endif //FFT_CACHING
63
64 if (is_complex() ) { // ming add 08/17/2010
65#ifdef NATIVE_FFT
66 LOGERR(" NATIVE_FFT does not support complex to complex."); // PAP
67 throw ImageFormatException("real image expected. Input image is complex image.");
68#else
69 EMData *temp_in=copy();
70 EMData *dat= copy_head();
71 int offset;
72 if(is_fftpadded()) {
73 offset = is_fftodd() ? 1 : 2;
74 }
75 else offset=0;
76 //printf("offset=%d\n",offset);
77 EMfft::complex_to_complex_nd(temp_in->get_data(),dat->get_data(),nx-offset,ny,nz);
78
79 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
80
81 dat->update();
82 delete temp_in;
83 EXITFUNC;
84 return dat;
85#endif // NATIVE_FFT
86 } else {
87 int nxreal = nx;
88 int offset = 2 - nx%2;
89 int nx2 = nx + offset;
90 EMData* dat = copy_head();
91 dat->set_size(nx2, ny, nz);
92 //dat->to_zero(); // do not need it, real_to_complex will do it right anyway
93 if (offset == 1) dat->set_fftodd(true);
94 else dat->set_fftodd(false);
95
96 float *d = dat->get_data();
97 //std::cout<<" do_fft "<<rdata[5]<<" "<<d[5]<<std::endl;
98 EMfft::real_to_complex_nd(get_data(), d, nxreal, ny, nz);
99
100 dat->update();
101 dat->set_fftpad(true);
102 dat->set_complex(true);
103 dat->set_attr("is_intensity",false);
104 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(true);
105 dat->set_ri(true);
106
107 EXITFUNC;
108#ifdef FFT_CACHING
109// printf("%p %d\n",this,nxyz);
110 if (nxyz<80000000) {
111 fftcache=dat->copy();
112 }
113#endif //FFT_CACHING
114 return dat;
115 }
116}
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
copy_head
EMData * copy_head() const
Make an image with a copy of the current image's header.
Definition: emdata_core.cpp:100
is_fftodd
bool is_fftodd() const
Does this image correspond to a (real-space) odd nx?
Definition: emdata_metadata.h:864
is_fftpadded
bool is_fftpadded() const
Is this image already extended along x for ffts?
Definition: emdata_metadata.h:831
is_complex
bool is_complex() const
Is this a complex image?
Definition: emdata_metadata.h:648
get_data
float * get_data() const
Get the image pixel density data in a 1D float array.
Definition: emdata_metadata.h:78
ImageFormatException
#define ImageFormatException(desc)
Definition: exception.h:147
LOGERR
#define LOGERR
Definition: log.h:51
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49

References copy(), copy_head(), ENTERFUNC, EXITFUNC, get_data(), ImageFormatException, is_complex(), is_fftodd(), is_fftpadded(), LOGERR, EMAN::EMData::nx, EMAN::EMData::nxyz, EMAN::EMData::ny, and EMAN::EMData::nz.

Referenced by bispecRotTransInvDirect(), bispecRotTransInvN(), EMAN::EMData::calc_mutual_correlation(), EMAN::EMData::convolute(), and EMAN::EMData::make_footprint().

◆ do_fft_inplace()

void EMData::do_fft_inplace ( )

Do FFT inplace.

And return the FFT image.

Returns
The FFT of the current image in real/imaginary format.

Definition at line 118 of file emdata_transform.cpp.

119{
120 ENTERFUNC;
121
122 if ( is_complex() ) {
123 LOGERR("real image expected. Input image is complex image.");
124 throw ImageFormatException("real image expected. Input image is complex image.");
125 }
126
127 size_t offset;
128 int nxreal;
129 get_data(); // Required call if GPU caching is being used. Otherwise harmless
130 if (!is_fftpadded()) {
131 // need to extend the matrix along x
132 // meaning nx is the un-fftpadded size
133 nxreal = nx;
134 offset = 2 - nx%2;
135 if (1 == offset) set_fftodd(true);
136 else set_fftodd(false);
137 int nxnew = nx + offset;
138 set_size(nxnew, ny, nz);
139 for (int iz = nz-1; iz >= 0; iz--) {
140 for (int iy = ny-1; iy >= 0; iy--) {
141 for (int ix = nxreal-1; ix >= 0; ix--) {
142 size_t oldxpos = ix + (iy + iz*ny)*(size_t)nxreal;
143 size_t newxpos = ix + (iy + iz*ny)*(size_t)nxnew;
144 (*this)(newxpos) = (*this)(oldxpos);
145 }
146 }
147 }
148 set_fftpad(true);
149 } else {
150 offset = is_fftodd() ? 1 : 2;
151 nxreal = nx - offset;
152 }
153 EMfft::real_to_complex_nd(rdata, rdata, nxreal, ny, nz);
154
155 set_complex(true);
156 if(ny==1 && nz==1) set_complex_x(true);
157 set_ri(true);
158
159 update();
160
161 EXITFUNC;
162 //return this;
163}
rdata
#define rdata(i)
Definition: analyzer.cpp:592
set_fftodd
void set_fftodd(bool is_fftodd)
Mark this image as having (real-space) odd nx.
Definition: emdata_metadata.h:882
set_fftpad
void set_fftpad(bool is_fftpadded)
Mark this image as already extended along x for ffts.
Definition: emdata_metadata.h:850
set_ri
void set_ri(bool is_ri)
Mark this image as a real/imaginary format complex image.
Definition: emdata_metadata.h:817
update
void update()
Mark EMData as changed, statistics, etc will be updated at need.
Definition: emdata_metadata.h:138
set_complex_x
void set_complex_x(bool is_complex_x)
Marks this image a 1D FFT image in X direction.
Definition: emdata_metadata.h:747
set_complex
void set_complex(bool is_complex)
Mark this image as a complex image.
Definition: emdata_metadata.h:712

References ENTERFUNC, EXITFUNC, get_data(), ImageFormatException, is_complex(), is_fftodd(), is_fftpadded(), LOGERR, EMAN::EMData::nx, EMAN::EMData::ny, EMAN::EMData::nz, EMAN::EMData::rdata, set_complex(), set_complex_x(), set_fftodd(), set_fftpad(), set_ri(), set_size(), and update().

◆ do_ift()

EMData * EMData::do_ift ( )

return the inverse fourier transform (IFT) image of the current image.

the current image may be changed if it is in amplitude/phase format as opposed to real/imaginary format - if this change is performed it is not undone.

Exceptions
ImageFormatExceptionIf the image is not a complex image.
Returns
The current image's inverse fourier transform image.

Definition at line 330 of file emdata_transform.cpp.

331{
332 ENTERFUNC;
333
334 if (!is_complex()) {
335 LOGERR("complex image expected. Input image is real image.");
336 throw ImageFormatException("complex image expected. Input image is real image.");
337 }
338
339 if (!is_ri()) {
340 LOGWARN("run IFT on AP data, only RI should be used. Converting.");
341 }
342
343 get_data(); // Required call if GPU caching is being used. Otherwise harmless
344 EMData* dat = copy_head();
345 dat->set_size(nx, ny, nz);
346 ap2ri();
347
348 float *d = dat->get_data();
349 int ndim = get_ndim();
350
351 /* Do inplace IFT on a image copy, because the complex to real transform of
352 * nd will destroy its input array even for out-of-place transforms.
353 */
354 memcpy((char *) d, (char *) rdata, (size_t)nx * ny * nz * sizeof(float));
355
356 int offset = is_fftodd() ? 1 : 2;
357 //cout << "Sending offset " << offset << " " << nx-offset << endl;
358 if (ndim == 1) {
359 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
360 } else {
361 EMfft::complex_to_real_nd(d, d, nx - offset, ny, nz);
362
363 size_t row_size = (nx - offset) * sizeof(float);
364 for (size_t i = 1; i < (size_t)ny * nz; i++) {
365 memmove((char *) &d[i * (nx - offset)], (char *) &d[i * nx], row_size);
366 }
367 }
368
369 dat->set_size(nx - offset, ny, nz); //remove the padding
370#if defined USE_FFTW3 //native fft and ACML already done normalization
371 // SCALE the inverse FFT
372 float scale = 1.0f / ((nx - offset) * ny * nz);
373 dat->mult(scale);
374#endif //USE_FFTW3
375 dat->set_fftodd(false);
376 dat->set_fftpad(false);
377 dat->set_complex(false);
378 if(dat->get_ysize()==1 && dat->get_zsize()==1) dat->set_complex_x(false);
379 dat->set_ri(false);
380 dat->set_attr("is_intensity",false);
381 dat->update();
382
383
384 EXITFUNC;
385 return dat;
386}
is_ri
bool is_ri() const
Is this image a real/imaginary format complex image?
Definition: emdata_metadata.h:797
get_ndim
int get_ndim() const
Get image dimension.
Definition: emdata_metadata.h:596
ap2ri
void ap2ri()
convert the complex image from amplitude/phase to real/imaginary
LOGWARN
#define LOGWARN
Definition: log.h:53

References ap2ri(), copy_head(), ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageFormatException, is_complex(), is_fftodd(), is_ri(), LOGERR, LOGWARN, EMAN::EMData::nx, EMAN::EMData::ny, EMAN::EMData::nz, EMAN::EMData::rdata, and EMAN::EMData::scale().

◆ do_ift_inplace()

void EMData::do_ift_inplace ( )

Definition at line 392 of file emdata_transform.cpp.

393{
394 ENTERFUNC;
395
396 if (!is_complex()) {
397 LOGERR("complex image expected. Input image is real image.");
398 throw ImageFormatException("complex image expected. Input image is real image.");
399 }
400
401 if (!is_ri()) {
402 LOGWARN("run IFT on AP data, only RI should be used. ");
403 }
404 ap2ri();
405
406 int offset = is_fftodd() ? 1 : 2;
407 float* data = get_data();
408 EMfft::complex_to_real_nd(data, data, nx - offset, ny, nz);
409
410#if defined USE_FFTW3 //native fft and ACML already done normalization
411 // SCALE the inverse FFT
412 int nxo = nx - offset;
413 float scale = 1.0f / ((size_t)nxo * ny * nz);
414 mult(scale);
415#endif //USE_FFTW3
416
417 set_fftpad(true);
418 set_complex(false);
419 if(ny==1 && nz==1) set_complex_x(false);
420 set_ri(false);
421 update();
422
423 EXITFUNC;
424// return this;
425}
mult
void mult(int n)
multiply an integer number to each pixel value of the image.
Definition: emdata_core.h:97

References ap2ri(), ENTERFUNC, EXITFUNC, get_data(), ImageFormatException, is_complex(), is_fftodd(), is_ri(), LOGERR, LOGWARN, mult(), EMAN::EMData::nx, EMAN::EMData::ny, EMAN::EMData::nz, EMAN::EMData::scale(), set_complex(), set_complex_x(), set_fftpad(), set_ri(), and update().

◆ insert_scaled_sum()

void insert_scaled_sum ( EMData *  block,
const FloatPoint &  center,
float  scale = 1.0,
float  mult_factor = 1.0 
)

Add a scaled image into another image at a specified location.

This is used, for example, to accumulate gaussians in programs like pdb2mrc.py. The center of 'block' will be positioned at 'center' with scale factor 'scale'. Densities will be interpolated in 'block' and multiplied by 'mult'.

Parameters
blockThe image to inserted.
centerThe center of the inserted block in 'this'.
scaleScale factor.
mult_factorNumber used to multiply the block's densities.
Exceptions
ImageDimensionExceptionIf 'this' image is not 2D/3D.

◆ render_amp24()

void EMData::render_amp24 ( int  x,
int  y,
int  xsize,
int  ysize,
int  bpl,
float  scale,
int  min_gray,
int  max_gray,
float  min_render,
float  max_render,
void *  ref,
void   cmapvoid *, int coord, unsigned char *tri 
)

Render the image into a 24-bit image.

2D image only.

Parameters
x
y
xsize
ysize
bpl
scale
min_gray
max_gray
min_render
max_render
ref
cmap
Exceptions
ImageDimensionExceptionIf the image is not 2D.

Definition at line 705 of file emdata_transform.cpp.

709{
710 ENTERFUNC;
711
712 if (get_ndim() != 2) {
713 throw ImageDimensionException("2D only");
714 }
715
716 if (is_complex()) {
717 ri2ap();
718 }
719
720 if (render_max <= render_min) {
721 render_max = render_min + 0.01f;
722 }
723
724 std::string ret=std::string();
725 ret.resize(iysize*bpl);
726 unsigned char *data=(unsigned char *)ret.data();
727
728 float rm = render_min;
729 float inv_scale = 1.0f / scale;
730 int ysize = iysize;
731 int xsize = ixsize;
732 const int scale_n = 100000;
733
734 int ymin = 0;
735 if ( iysize * inv_scale > ny) {
736 ymin = (int) (iysize - ny / inv_scale);
737 }
738 float gs = (maxgray - mingray) / (render_max - render_min);
739 if (render_max < render_min) {
740 gs = 0;
741 rm = FLT_MAX;
742 }
743 int dsx = -1;
744 int dsy = 0;
745 if (inv_scale == floor(inv_scale)) {
746 dsx = (int) inv_scale;
747 dsy = (int) (inv_scale * nx);
748 }
749 int addi = 0;
750 int addr = 0;
751
752 if (dsx == -1) {
753 addi = (int) floor(inv_scale);
754 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
755 }
756
757 int remx = 0;
758 int remy = 0;
759 int xmin = 0;
760 if (x0 < 0) {
761 xmin = (int) (-x0 / inv_scale);
762 xsize -= (int) floor(x0 / inv_scale);
763 x0 = 0;
764 }
765
766 if ((xsize - xmin) * inv_scale > (nx - x0)) {
767 xsize = (int) ((nx - x0) / inv_scale + xmin);
768 }
769 int ymax = ysize - 1;
770 if (y0 < 0) {
771 ymax = (int) (ysize + y0 / inv_scale - 1);
772 ymin += (int) floor(y0 / inv_scale);
773 y0 = 0;
774 }
775
776
777 if (xmin < 0) {
778 xmin = 0;
779 }
780
781 if (ymin < 0) {
782 ymin = 0;
783 }
784 if (xsize > ixsize) {
785 xsize = ixsize;
786 }
787 if (ymax > iysize) {
788 ymax = iysize;
789 }
790
791 int lmax = nx * ny - 1;
792 unsigned char tri[3];
793 float* image_data = get_data();
794 if (is_complex()) {
795 if (dsx != -1) {
796 int l = y0 * nx;
797 for (int j = ymax; j >= ymin; j--) {
798 int ll = x0;
799 for (int i = xmin; i < xsize; i++, ll += dsx) {
800 if (l + ll > lmax || ll >= nx - 2) {
801 break;
802 }
803 int kk = 0;
804 if (ll >= nx / 2) {
805 if (l >= (ny - inv_scale) * nx) {
806 kk = 2 * (ll - nx / 2) + 2;
807 }
808 else {
809 kk = 2 * (ll - nx / 2) + l + 2 + nx;
810 }
811 }
812 else {
813 kk = nx * ny - (l + 2 * ll) - 2;
814 }
815 int k = 0;
816 float t = image_data[kk];
817 if (t <= rm) {
818 k = mingray;
819 }
820 else if (t >= render_max) {
821 k = maxgray;
822 }
823 else {
824 k = (int) (gs * (t - render_min));
825 k += mingray;
826 }
827 tri[0] = static_cast < unsigned char >(k);
828 cmap(ref, kk, tri);
829 data[i * 3 + j * bpl] = tri[0];
830 data[i * 3 + 1 + j * bpl] = tri[1];
831 data[i * 3 + 2 + j * bpl] = tri[2];
832 }
833 l += dsy;
834 }
835 }
836 else {
837 remy = 10;
838 for (int j = ymax, l = y0 * nx; j >= ymin; j--) {
839 int br = l;
840 remx = 10;
841 for (int i = xmin, ll = x0; i < xsize - 1; i++) {
842 if (l + ll > lmax || ll >= nx - 2) {
843 break;
844 }
845 int kk = 0;
846 if (ll >= nx / 2) {
847 if (l >= (ny * nx - nx)) {
848 kk = 2 * (ll - nx / 2) + 2;
849 }
850 else {
851 kk = 2 * (ll - nx / 2) + l + 2 + nx;
852 }
853 }
854 else {
855 kk = nx * ny - (l + 2 * ll) - 2;
856 }
857 int k = 0;
858 float t = image_data[kk];
859 if (t <= rm) {
860 k = mingray;
861 }
862 else if (t >= render_max) {
863 k = maxgray;
864 }
865 else {
866 k = (int) (gs * (t - render_min));
867 k += mingray;
868 }
869 tri[0] = static_cast < unsigned char >(k);
870 cmap(ref, kk, tri);
871 data[i * 3 + j * bpl] = tri[0];
872 data[i * 3 + 1 + j * bpl] = tri[1];
873 data[i * 3 + 2 + j * bpl] = tri[2];
874 ll += addi;
875 remx += addr;
876 if (remx > scale_n) {
877 remx -= scale_n;
878 ll++;
879 }
880 }
881 l = br + addi * nx;
882 remy += addr;
883 if (remy > scale_n) {
884 remy -= scale_n;
885 l += nx;
886 }
887 }
888 }
889 }
890 else {
891 if (dsx != -1) {
892 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
893 int br = l;
894 for (int i = xmin; i < xsize; i++, l += dsx) {
895 if (l > lmax) {
896 break;
897 }
898 float t = image_data[l];
899 int k = 0;
900 if (t <= rm) {
901 k = mingray;
902 }
903 else if (t >= render_max) {
904 k = maxgray;
905 }
906 else {
907 k = (int) (gs * (t - render_min));
908 k += mingray;
909 }
910 tri[0] = static_cast < unsigned char >(k);
911 cmap(ref, l, tri);
912 data[i * 3 + j * bpl] = tri[0];
913 data[i * 3 + 1 + j * bpl] = tri[1];
914 data[i * 3 + 2 + j * bpl] = tri[2];
915 }
916 l = br + dsy;
917 }
918 }
919 else {
920 remy = 10;
921 for (int j = ymax, l = x0 + y0 * nx; j >= ymin; j--) {
922 int br = l;
923 remx = 10;
924 for (int i = xmin; i < xsize; i++) {
925 if (l > lmax) {
926 break;
927 }
928 float t = image_data[l];
929 int k = 0;
930 if (t <= rm) {
931 k = mingray;
932 }
933 else if (t >= render_max) {
934 k = maxgray;
935 }
936 else {
937 k = (int) (gs * (t - render_min));
938 k += mingray;
939 }
940 tri[0] = static_cast < unsigned char >(k);
941 cmap(ref, l, tri);
942 data[i * 3 + j * bpl] = tri[0];
943 data[i * 3 + 1 + j * bpl] = tri[1];
944 data[i * 3 + 2 + j * bpl] = tri[2];
945 l += addi;
946 remx += addr;
947 if (remx > scale_n) {
948 remx -= scale_n;
949 l++;
950 }
951 }
952 l = br + addi * nx;
953 remy += addr;
954 if (remy > scale_n) {
955 remy -= scale_n;
956 l += nx;
957 }
958 }
959 }
960 }
961
962 EXITFUNC;
963}
ri2ap
void ri2ap()
convert the complex image from real/imaginary to amplitude/phase
Definition: emdata_transform.cpp:999
ImageDimensionException
#define ImageDimensionException(desc)
Definition: exception.h:166

References ENTERFUNC, EXITFUNC, get_data(), get_ndim(), ImageDimensionException, is_complex(), EMAN::EMData::nx, EMAN::EMData::ny, ri2ap(), and EMAN::EMData::scale().

◆ render_amp8()

std::string render_amp8 ( int  x,
int  y,
int  xsize,
int  ysize,
int  bpl,
float  scale,
int  min_gray,
int  max_gray,
float  min_render,
float  max_render,
float  gamma,
int  flags 
)

Render the image into an 8-bit image.

2D images only. flags provide a way to do unusual things with this function, such as calculating a histogram of the rendered area.

Parameters
xorigin of the area to render
y
xsizesize of the area to render in output pixels
ysize
bplbytes per line, if asrgb remember *3
scalescale factor for rendering
min_grayminimum gray value to render (0-255)
max_graymaximum gray value to render (0-255)
min_renderfloat image density corresponding to min_gray
max_renderfloat image density corresponding to max_gray
gamma
flags1-duplicate each output pixel 3x for RGB rendering,2-add a 256 int greyscale histogram to the end of the image array,4-invert y axis,8-render 32 bit 0xffRRGGBB
Exceptions
ImageDimensionExceptionIf the image is not 2D.

◆ render_ap24()

EMBytes EMData::render_ap24 ( int  x,
int  y,
int  xsize,
int  ysize,
int  bpl,
float  scale,
int  min_gray,
int  max_gray,
float  min_render,
float  max_render,
float  gamma,
int  flags 
)

Render the image into an 8-bit image.

2D images only. flags provide a way to do unusual things with this function, such as calculating a histogram of the rendered area.

Parameters
xorigin of the area to render
y
xsizesize of the area to render in output pixels
ysize
bplbytes per line, if asrgb remember *3
scalescale factor for rendering
min_grayminimum gray value to render (0-255)
max_graymaximum gray value to render (0-255)
min_renderfloat image density corresponding to min_gray
max_renderfloat image density corresponding to max_gray
gamma
flags1-duplicate each output pixel 3x for RGB rendering,2-add a 256 int greyscale histogram to the end of the image array,4-invert y axis,8-render 32 bit 0xffRRGGBB
Exceptions
ImageDimensionExceptionIf the image is not 2D.

Definition at line 429 of file emdata_transform.cpp.

432{
433 ENTERFUNC;
434
435 int asrgb;
436 int hist=(flags&2)/2;
437 int invy=(flags&4)?1:0;
438
439 if (!is_complex()) throw ImageDimensionException("complex only");
440
441 if (get_ndim() != 2) {
442 throw ImageDimensionException("2D only");
443 }
444
445 if (is_complex()) ri2ap();
446
447 if (render_max <= render_min) {
448 render_max = render_min + 0.01f;
449 }
450
451 if (gamma<=0) gamma=1.0;
452
453 // Calculating a full floating point gamma for
454 // each pixel in the image slows rendering unacceptably
455 // however, applying a gamma-mapping to an 8 bit colorspace
456 // has unaccepable accuracy. So, we oversample the 8 bit colorspace
457 // as a 12 bit colorspace and apply the gamma mapping to that
458 // This should produce good accuracy for gamma values
459 // larger than 0.5 (and a high upper limit)
460 static int smg0=0,smg1=0; // while this destroys threadsafety in the rendering process
461 static float sgam=0; // it is necessary for speed when rendering large numbers of small images
462 static unsigned char gammamap[4096];
463 if (gamma!=1.0 && (smg0!=mingray || smg1!=maxgray || sgam!=gamma)) {
464 for (int i=0; i<4096; i++) {
465 if (mingray<maxgray) gammamap[i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
466 else gammamap[4095-i]=(unsigned char)(mingray+(maxgray-mingray+0.999)*pow(((float)i/4096.0f),gamma));
467 }
468 }
469 smg0=mingray; // so we don't recompute the map unless something changes
470 smg1=maxgray;
471 sgam=gamma;
472
473 if (flags&8) asrgb=4;
474 else if (flags&1) asrgb=3;
475 else throw ImageDimensionException("must set flag 1 or 8");
476
477 EMBytes ret=EMBytes();
478// ret.resize(iysize*bpl);
479 ret.assign(iysize*bpl+hist*1024,char(mingray));
480 unsigned char *data=(unsigned char *)ret.data();
481 unsigned int *histd=(unsigned int *)(data+iysize*bpl);
482 if (hist) {
483 for (int i=0; i<256; i++) histd[i]=0;
484 }
485
486 float rm = render_min;
487 float inv_scale = 1.0f / scale;
488 int ysize = iysize;
489 int xsize = ixsize;
490
491 int ymin = 0;
492 if (iysize * inv_scale > ny) {
493 ymin = (int) (iysize - ny / inv_scale);
494 }
495
496 float gs = (maxgray - mingray) / (render_max - render_min);
497 float gs2 = 4095.999f / (render_max - render_min);
498// float gs2 = 1.0 / (render_max - render_min);
499 if (render_max < render_min) {
500 gs = 0;
501 rm = FLT_MAX;
502 }
503
504 int dsx = -1;
505 int dsy = 0;
506 int remx = 0;
507 int remy = 0;
508 const int scale_n = 100000;
509
510 int addi = 0;
511 int addr = 0;
512 if (inv_scale == floor(inv_scale)) {
513 dsx = (int) inv_scale;
514 dsy = (int) (inv_scale * nx);
515 }
516 else {
517 addi = (int) floor(inv_scale);
518 addr = (int) (scale_n * (inv_scale - floor(inv_scale)));
519 }
520
521 int xmin = 0;
522 if (x0 < 0) {
523 xmin = (int) (-x0 / inv_scale);
524 xsize -= (int) floor(x0 / inv_scale);
525 x0 = 0;
526 }
527
528 if ((xsize - xmin) * inv_scale > (nx - x0)) {
529 xsize = (int) ((nx - x0) / inv_scale + xmin);
530 }
531 int ymax = ysize - 1;
532 if (y0 < 0) {
533 ymax = (int) (ysize + y0 / inv_scale - 1);
534 ymin += (int) floor(y0 / inv_scale);
535 y0 = 0;
536 }
537
538 if (xmin < 0) xmin = 0;
539 if (ymin < 0) ymin = 0;
540 if (xsize > ixsize) xsize = ixsize;
541 if (ymax > iysize) ymax = iysize;
542
543 int lmax = nx * ny - 1;
544
545 int mid=nx*ny/2;
546 float* image_data = get_data();
547 if (dsx != -1) {
548 int l = y0 * nx;
549 for (int j = ymax; j >= ymin; j--) {
550 int ll = x0;
551 for (int i = xmin; i < xsize; i++) {
552 if (l + ll > lmax || ll >= nx - 2) break;
553
554 int k = 0;
555 unsigned char p;
556 int ph;
557 if (ll >= nx / 2) {
558 if (l >= (ny - inv_scale) * nx) k = 2 * (ll - nx / 2) + 2;
559 else k = 2 * (ll - nx / 2) + l + 2 + nx;
560 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
561 else k+=mid;
562 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
563 }
564 else {
565 k = nx * ny - (l + 2 * ll) - 2;
566 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
567 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
568 else k+=mid;
569 }
570 float t = image_data[k];
571 if (t <= rm) p = mingray;
572 else if (t >= render_max) p = maxgray;
573 else if (gamma!=1.0) {
574 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
575 p = gammamap[k]; // apply gamma using precomputed gamma map
576 }
577 else {
578 p = (unsigned char) (gs * (t - render_min));
579 p += mingray;
580 }
581 if (ph<256) {
582 data[i * asrgb + j * bpl] = p*(255-ph)/256;
583 data[i * asrgb + j * bpl+1] = p*ph/256;
584 data[i * asrgb + j * bpl+2] = 0;
585 }
586 else if (ph<512) {
587 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
588 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
589 data[i * asrgb + j * bpl] = 0;
590 }
591 else {
592 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
593 data[i * asrgb + j * bpl] = p*(ph-512)/256;
594 data[i * asrgb + j * bpl+1] = 0;
595 }
596 if (hist) histd[p]++;
597 ll += dsx;
598 }
599 l += dsy;
600 }
601 }
602 else {
603 remy = 10;
604 int l = y0 * nx;
605 for (int j = ymax; j >= ymin; j--) {
606 int br = l;
607 remx = 10;
608 int ll = x0;
609 for (int i = xmin; i < xsize - 1; i++) {
610 if (l + ll > lmax || ll >= nx - 2) {
611 break;
612 }
613 int k = 0;
614 unsigned char p;
615 int ph;
616 if (ll >= nx / 2) {
617 if (l >= (ny * nx - nx)) k = 2 * (ll - nx / 2) + 2;
618 else k = 2 * (ll - nx / 2) + l + 2 + nx;
619 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
620 else k+=mid;
621 ph = (int)(image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
622 }
623 else {
624 k = nx * ny - (l + 2 * ll) - 2;
625 if (k>=mid) k-=mid; // These 2 lines handle the Fourier origin being in the corner, not the middle
626 else k+=mid;
627 ph = (int)(-image_data[k+1]*768/(2.0*M_PI))+384; // complex phase as integer 0-767
628 }
629
630 float t = image_data[k];
631 if (t <= rm)
632 p = mingray;
633 else if (t >= render_max) {
634 p = maxgray;
635 }
636 else if (gamma!=1.0) {
637 k=(int)(gs2 * (t-render_min)); // map float value to 0-4096 range
638 p = gammamap[k]; // apply gamma using precomputed gamma map
639 }
640 else {
641 p = (unsigned char) (gs * (t - render_min));
642 p += mingray;
643 }
644 if (ph<256) {
645 data[i * asrgb + j * bpl] = p*(255-ph)/256;
646 data[i * asrgb + j * bpl+1] = p*ph/256;
647 data[i * asrgb + j * bpl+2] = 0;
648 }
649 else if (ph<512) {
650 data[i * asrgb + j * bpl+1] = p*(511-ph)/256;
651 data[i * asrgb + j * bpl+2] = p*(ph-256)/256;
652 data[i * asrgb + j * bpl] = 0;
653 }
654 else {
655 data[i * asrgb + j * bpl+2] = p*(767-ph)/256;
656 data[i * asrgb + j * bpl] = p*(ph-512)/256;
657 data[i * asrgb + j * bpl+1] = 0;
658 }
659 if (hist) histd[p]++;
660 ll += addi;
661 remx += addr;
662 if (remx > scale_n) {
663 remx -= scale_n;
664 ll++;
665 }
666 }
667 l = br + addi * nx;
668 remy += addr;
669 if (remy > scale_n) {
670 remy -= scale_n;
671 l += nx;
672 }
673 }
674 }
675
676 // this replicates r -> g,b
677 if (asrgb==4) {
678 for (int j=ymin*bpl; j<=ymax*bpl; j+=bpl) {
679 for (int i=xmin; i<xsize*4; i+=4) {
680 data[i+j+3]=255;
681 }
682 }
683 }
684
685 EXITFUNC;
686
687 // ok, ok, not the most efficient place to do this, but it works
688 if (invy) {
689 int x,y;
690 char swp;
691 for (y=0; y<iysize/2; y++) {
692 for (x=0; x<ixsize; x++) {
693 swp=ret[y*bpl+x];
694 ret[y*bpl+x]=ret[(iysize-y-1)*bpl+x];
695 ret[(iysize-y-1)*bpl+x]=swp;
696 }
697 }
698 }
699
700 // return PyString_FromStringAndSize((const char*) data,iysize*bpl);
701 return ret;
702}
y
#define y(i, j)
Definition: projector.cpp:1516

References ENTERFUNC, EXITFUNC, EMAN::EMData::flags, get_data(), get_ndim(), ImageDimensionException, is_complex(), EMAN::EMData::nx, EMAN::EMData::ny, ri2ap(), EMAN::EMData::scale(), x, and y.

◆ ri2ap()

void EMData::ri2ap ( )

convert the complex image from real/imaginary to amplitude/phase

Definition at line 999 of file emdata_transform.cpp.

1000{
1001 ENTERFUNC;
1002
1003 if (!is_complex() || !is_ri()) {
1004 return;
1005 }
1006
1007 float * data = get_data();
1008
1009 size_t size = (size_t)nx * ny * nz;
1010 for (size_t i = 0; i < size; i += 2) {
1011 float f = (float)hypot(data[i], data[i + 1]);
1012 if (data[i] == 0 && data[i + 1] == 0) {
1013 data[i + 1] = 0;
1014 }
1015 else {
1016 data[i + 1] = atan2(data[i + 1], data[i]);
1017 }
1018 data[i] = f;
1019 }
1020
1021 set_ri(false);
1022 update();
1023 EXITFUNC;
1024}

References ENTERFUNC, EXITFUNC, get_data(), is_complex(), is_ri(), EMAN::EMData::nx, EMAN::EMData::ny, EMAN::EMData::nz, set_ri(), and update().

Referenced by EMAN::EMData::add_incoherent(), get_fft_amplitude(), get_fft_phase(), render_amp24(), and render_ap24().

◆ ri2inten()

void EMData::ri2inten ( )

convert the complex image from real/imaginary to Intensity/0.

This conversion cannot be reversed, and the image remains marked as R/I. Also sets the is_intensity flag, which is used by routines like calc_radial_dist to avoid returning intensity squared.

Definition at line 979 of file emdata_transform.cpp.

980{
981 ENTERFUNC;
982
983 if (!is_complex()) return;
984 if (!is_ri()) ap2ri();
985
986 float * data = get_data();
987 size_t size = (size_t)nx * ny * nz;
988 for (size_t i = 0; i < size; i += 2) {
989 data[i]=data[i]*data[i]+data[i+1]*data[i+1];
990 data[i+1]=0;
991 }
992
993 set_attr("is_intensity", int(1));
994 update();
995 EXITFUNC;
996}
set_attr
void set_attr(const string &key, EMObject val)
Set a header attribute's value.

References ap2ri(), ENTERFUNC, EXITFUNC, get_data(), is_complex(), is_ri(), EMAN::EMData::nx, EMAN::EMData::ny, EMAN::EMData::nz, set_attr(), and update().

Generated by doxygen 1.9.2