Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
B
BioEM
Project overview
Project overview
Details
Activity
Releases
Repository
Repository
Files
Commits
Branches
Tags
Contributors
Graph
Compare
Issues
0
Issues
0
List
Boards
Labels
Service Desk
Milestones
Merge Requests
0
Merge Requests
0
CI / CD
CI / CD
Pipelines
Jobs
Schedules
Operations
Operations
Incidents
Environments
Analytics
Analytics
CI / CD
Repository
Value Stream
Wiki
Wiki
Members
Members
Collapse sidebar
Close sidebar
Activity
Graph
Create a new issue
Jobs
Commits
Issue Boards
Open sidebar
MPIBP-Hummer
BioEM
Commits
8639adb7
Commit
8639adb7
authored
Apr 20, 2014
by
David Rohr
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
fix file endings
parent
4a123b50
Changes
4
Hide whitespace changes
Inline
Side-by-side
Showing
4 changed files
with
469 additions
and
469 deletions
+469
-469
bioem.cpp
bioem.cpp
+30
-30
include/map.h
include/map.h
+7
-7
map.cpp
map.cpp
+312
-312
param.cpp
param.cpp
+120
-120
No files found.
bioem.cpp
View file @
8639adb7
...
...
@@ -57,8 +57,8 @@ int bioem::configure(int ac, char* av[])
param
.
writeAngles
=
false
;
param
.
dumpMap
=
false
;
param
.
loadMap
=
false
;
RefMap
.
readMRC
=
false
;
RefMap
.
readMultMRC
=
false
;
RefMap
.
readMRC
=
false
;
RefMap
.
readMultMRC
=
false
;
// *************************************************************************************
cout
<<
" ++++++++++++ FROM COMMAND LINE +++++++++++
\n\n
"
;
...
...
@@ -68,17 +68,17 @@ int bioem::configure(int ac, char* av[])
try
{
po
::
options_description
desc
(
"Command line inputs"
);
desc
.
add_options
()
(
"Inputfile"
,
po
::
value
<
std
::
string
>
(),
"(Mandatory) Name of input parameter file"
)
(
"Modelfile"
,
po
::
value
<
std
::
string
>
()
,
"(Mandatory) Name of model file"
)
(
"Particlesfile"
,
po
::
value
<
std
::
string
>
(),
"(Mandatory) Name of paricles file"
)
(
"ReadPDB"
,
"(Optional) If reading model file in PDB format"
)
(
"ReadMRC"
,
"(Optional) If reading particle file in MRC format"
)
(
"ReadMultipleMRC"
,
"(Optional) If reading Multiple MRCs"
)
(
"DumpMaps"
,
"(Optional) Dump maps after they were red from maps file"
)
(
"LoadMapDump"
,
"(Optional) Read Maps from dump instead of maps file"
)
(
"help"
,
"(Optional) Produce help message"
)
;
desc
.
add_options
()
(
"Inputfile"
,
po
::
value
<
std
::
string
>
(),
"(Mandatory) Name of input parameter file"
)
(
"Modelfile"
,
po
::
value
<
std
::
string
>
()
,
"(Mandatory) Name of model file"
)
(
"Particlesfile"
,
po
::
value
<
std
::
string
>
(),
"(Mandatory) Name of paricles file"
)
(
"ReadPDB"
,
"(Optional) If reading model file in PDB format"
)
(
"ReadMRC"
,
"(Optional) If reading particle file in MRC format"
)
(
"ReadMultipleMRC"
,
"(Optional) If reading Multiple MRCs"
)
(
"DumpMaps"
,
"(Optional) Dump maps after they were red from maps file"
)
(
"LoadMapDump"
,
"(Optional) Read Maps from dump instead of maps file"
)
(
"help"
,
"(Optional) Produce help message"
)
;
po
::
positional_options_description
p
;
...
...
@@ -86,8 +86,8 @@ int bioem::configure(int ac, char* av[])
p
.
add
(
"Modelfile"
,
-
1
);
p
.
add
(
"Particlesfile"
,
-
1
);
p
.
add
(
"ReadPDB"
,
-
1
);
p
.
add
(
"ReadMRC"
,
-
1
);
p
.
add
(
"ReadMultipleMRC"
,
-
1
);
p
.
add
(
"ReadMRC"
,
-
1
);
p
.
add
(
"ReadMultipleMRC"
,
-
1
);
p
.
add
(
"DumpMaps"
,
-
1
);
p
.
add
(
"LoadMapDump"
,
-
1
);
...
...
@@ -125,17 +125,17 @@ int bioem::configure(int ac, char* av[])
Model
.
readPDB
=
true
;
}
if
(
vm
.
count
(
"ReadMRC"
))
{
cout
<<
"Reading particle file in MRC format.
\n
"
;
RefMap
.
readMRC
=
true
;
}
if
(
vm
.
count
(
"ReadMultipleMRC"
))
{
cout
<<
"Reading Multiple MRCs.
\n
"
;
RefMap
.
readMultMRC
=
true
;
}
if
(
vm
.
count
(
"ReadMRC"
))
{
cout
<<
"Reading particle file in MRC format.
\n
"
;
RefMap
.
readMRC
=
true
;
}
if
(
vm
.
count
(
"ReadMultipleMRC"
))
{
cout
<<
"Reading Multiple MRCs.
\n
"
;
RefMap
.
readMultMRC
=
true
;
}
if
(
vm
.
count
(
"DumpMaps"
))
{
...
...
@@ -379,8 +379,8 @@ int bioem::run()
int
bioem
::
compareRefMaps
(
int
iOrient
,
int
iConv
,
const
myfloat_t
*
conv_map
,
mycomplex_t
*
localmultFFT
,
myfloat_t
sumC
,
myfloat_t
sumsquareC
,
const
int
startMap
)
{
//***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps *****
//***************************************************************************************
//***** BioEM routine for comparing reference maps to convoluted maps *****
if
(
FFTAlgo
)
{
//With FFT Algorithm
...
...
@@ -405,8 +405,8 @@ int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, myc
inline
void
bioem
::
calculateCCFFT
(
int
iRefMap
,
int
iOrient
,
int
iConv
,
myfloat_t
sumC
,
myfloat_t
sumsquareC
,
mycomplex_t
*
localConvFFT
,
mycomplex_t
*
localCCT
,
myfloat_t
*
lCC
)
{
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
const
mycomplex_t
*
RefMapFFT
=
&
RefMap
.
RefMapsFFT
[
iRefMap
*
param
.
FFTMapSize
];
for
(
int
i
=
0
;
i
<
param
.
param_device
.
NumberPixels
*
param
.
param_device
.
NumberFFTPixels1D
;
i
++
)
...
...
include/map.h
View file @
8639adb7
...
...
@@ -34,16 +34,16 @@ public:
int
precalculate
(
bioem_param
&
param
,
bioem
&
bio
);
int
PreCalculateMapsFFT
(
bioem_param
&
param
);
int
read_int
(
int
*
currlong
,
FILE
*
fin
,
int
swap
);
int
read_float
(
float
*
currfloat
,
FILE
*
fin
,
int
swap
);
int
read_float_empty
(
FILE
*
fin
);
int
read_char_float
(
float
*
currfloat
,
FILE
*
fin
)
;
int
test_mrc
(
const
char
*
vol_file
,
int
swap
);
int
read_MRC
(
const
char
*
filename
,
bioem_param
&
param
);
int
read_int
(
int
*
currlong
,
FILE
*
fin
,
int
swap
);
int
read_float
(
float
*
currfloat
,
FILE
*
fin
,
int
swap
);
int
read_float_empty
(
FILE
*
fin
);
int
read_char_float
(
float
*
currfloat
,
FILE
*
fin
)
;
int
test_mrc
(
const
char
*
vol_file
,
int
swap
);
int
read_MRC
(
const
char
*
filename
,
bioem_param
&
param
);
mycomplex_t
*
RefMapsFFT
;
bool
readMRC
,
readMultMRC
;
bool
readMRC
,
readMultMRC
;
int
ntotRefMap
;
int
numPixels
;
...
...
map.cpp
View file @
8639adb7
...
...
@@ -35,54 +35,54 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
cout
<<
"Particle Maps read from Map Dump
\n
"
;
}
else
if
(
readMRC
)
{
ntotRefMap
=
0
;
if
(
readMultMRC
)
{
ifstream
input
(
filemap
);
if
(
!
input
.
good
())
{
cout
<<
"Failed to open file contaning MRC names: "
<<
filemap
<<
"
\n
"
;
exit
(
0
);
}
char
line
[
512
]
=
{
0
};
char
mapname
[
100
];
char
tmpm
[
10
]
=
{
0
};
const
char
*
indifile
;
while
(
!
input
.
eof
())
{
input
.
getline
(
line
,
512
);
char
tmpVals
[
100
]
=
{
' '
};
strncpy
(
tmpVals
,
line
,
100
);
sscanf
(
tmpVals
,
"%99c"
,
mapname
);
// Check for last line
strncpy
(
tmpm
,
mapname
,
3
);
if
(
strcmp
(
tmpm
,
"XXX"
)
!=
0
)
{
indifile
=
mapname
;
//Reading Multiple MRC
read_MRC
(
indifile
,
param
);
}
for
(
int
i
=
0
;
i
<
3
;
i
++
)
mapname
[
i
]
=
'X'
;
for
(
int
i
=
3
;
i
<
100
;
i
++
)
mapname
[
i
]
=
{
0
};
}
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from MULTIPLE MRC Files in: "
<<
filemap
<<
"
\n
"
;
}
else
{
read_MRC
(
filemap
,
param
);
cout
<<
"
\n
++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from ONE MRC File: "
<<
filemap
<<
"
\n
"
;
}
}
else
if
(
readMRC
)
{
ntotRefMap
=
0
;
if
(
readMultMRC
)
{
ifstream
input
(
filemap
);
if
(
!
input
.
good
())
{
cout
<<
"Failed to open file contaning MRC names: "
<<
filemap
<<
"
\n
"
;
exit
(
0
);
}
char
line
[
512
]
=
{
0
};
char
mapname
[
100
];
char
tmpm
[
10
]
=
{
0
};
const
char
*
indifile
;
while
(
!
input
.
eof
())
{
input
.
getline
(
line
,
512
);
char
tmpVals
[
100
]
=
{
' '
};
strncpy
(
tmpVals
,
line
,
100
);
sscanf
(
tmpVals
,
"%99c"
,
mapname
);
// Check for last line
strncpy
(
tmpm
,
mapname
,
3
);
if
(
strcmp
(
tmpm
,
"XXX"
)
!=
0
)
{
indifile
=
mapname
;
//Reading Multiple MRC
read_MRC
(
indifile
,
param
);
}
for
(
int
i
=
0
;
i
<
3
;
i
++
)
mapname
[
i
]
=
'X'
;
for
(
int
i
=
3
;
i
<
100
;
i
++
)
mapname
[
i
]
=
{
0
};
}
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from MULTIPLE MRC Files in: "
<<
filemap
<<
"
\n
"
;
}
else
{
read_MRC
(
filemap
,
param
);
cout
<<
"
\n
++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"Particle Maps read from ONE MRC File: "
<<
filemap
<<
"
\n
"
;
}
}
else
{
int
nummap
=
-
1
;
...
...
@@ -177,7 +177,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
}
cout
<<
"Total Number of particles: "
<<
ntotRefMap
;
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
;
cout
<<
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
;
sum_RefMap
=
(
myfloat_t
*
)
mallocchk
(
sizeof
(
myfloat_t
)
*
ntotRefMap
);
sumsquare_RefMap
=
(
myfloat_t
*
)
mallocchk
(
sizeof
(
myfloat_t
)
*
ntotRefMap
);
...
...
@@ -265,266 +265,266 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio)
set_pointers
();
}
int
bioem_RefMap
::
read_MRC
(
const
char
*
filename
,
bioem_param
&
param
)
{
myfloat_t
st
,
st2
;
unsigned
long
count
;
FILE
*
fin
;
float
currfloat
;
int
nc
,
nr
,
ns
,
swap
,
header_ok
=
1
;
float
xlen
,
ylen
,
zlen
;
int
mode
,
ncstart
,
nrstart
,
nsstart
,
ispg
,
nsymbt
,
lskflg
;
float
a_tmp
,
b_tmp
,
g_tmp
;
int
mx
,
my
,
mz
,
mapc
,
mapr
,
maps
;
float
dmin
,
dmax
,
dmean
;
int
n_range_viol0
,
n_range_viol1
;
fin
=
fopen
(
filename
,
"rb"
);
if
(
fin
==
NULL
)
{
cout
<<
"ERROR opening MRC: "
<<
filename
;
exit
(
1
);
}
n_range_viol0
=
test_mrc
(
filename
,
0
);
n_range_viol1
=
test_mrc
(
filename
,
1
);
if
(
n_range_viol0
<
n_range_viol1
)
{
//* guess endianism
swap
=
0
;
if
(
n_range_viol0
>
0
)
{
printf
(
" Warning: %i header field range violations detected in file %s
\n
"
,
n_range_viol0
,
filename
);
}
}
else
{
swap
=
1
;
if
(
n_range_viol1
>
0
)
{
printf
(
"Warning: %i header field range violations detected in file %s
\n
"
,
n_range_viol1
,
filename
);
}
}
printf
(
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
);
printf
(
"Reading Information from MRC: %s
\n
"
,
filename
);
header_ok
*=
read_int
(
&
nc
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nr
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ns
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mode
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ncstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nrstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nsstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mx
,
fin
,
swap
);
header_ok
*=
read_int
(
&
my
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mz
,
fin
,
swap
);
header_ok
*=
read_float
(
&
xlen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
ylen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
zlen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
a_tmp
,
fin
,
swap
);
header_ok
*=
read_float
(
&
b_tmp
,
fin
,
swap
);
header_ok
*=
read_float
(
&
g_tmp
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mapc
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mapr
,
fin
,
swap
);
header_ok
*=
read_int
(
&
maps
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmin
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmax
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmean
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ispg
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nsymbt
,
fin
,
swap
);
header_ok
*=
read_int
(
&
lskflg
,
fin
,
swap
);
printf
(
"Number Columns = %8d
\n
"
,
nc
);
printf
(
"Number Rows = %8d
\n
"
,
nr
);
printf
(
"Number Sections = %8d
\n
"
,
ns
);
printf
(
"MODE = %4d (only data type mode 2: 32-bit)
\n
"
,
mode
);
printf
(
"NSYMBT = %4d (# bytes symmetry operators)
\n
"
,
nsymbt
);
/* printf(" NCSTART = %8d (index of first column, counting from 0)\n",ncstart);
printf("> NRSTART = %8d (index of first row, counting from 0)\n",nrstart);
printf(" NSSTART = %8d (index of first section, counting from 0)\n",nsstart);
printf(" MX = %8d (# of X intervals in unit cell)\n",mx);
printf(" MY = %8d (# of Y intervals in unit cell)\n",my);
printf(" MZ = %8d (# of Z intervals in unit cell)\n",mz);
printf(" X length = %8.3f (unit cell dimension)\n",xlen);
printf(" Y length = %8.3f (unit cell dimension)\n",ylen);
printf(" Z length = %8.3f (unit cell dimension)\n",zlen);
printf(" Alpha = %8.3f (unit cell angle)\n",a_tmp);
printf(" Beta = %8.3f (unit cell angle)\n",b_tmp);
printf(" Gamma = %8.3f (unit cell angle)\n",g_tmp);
printf(" MAPC = %8d (columns axis: 1=X,2=Y,3=Z)\n",mapc);
printf(" MAPR = %8d (rows axis: 1=X,2=Y,3=Z)\n",mapr);
printf(" MAPS = %8d (sections axis: 1=X,2=Y,3=Z)\n",maps);
printf(" DMIN = %8.3f (minimum density value - ignored)\n",dmin);
printf(" DMAX = %8.3f (maximum density value - ignored)\n",dmax);
printf(" DMEAN = %8.3f (mean density value - ignored)\n",dmean);
printf(" ISPG = %8d (space group number - ignored)\n",ispg);
printf(" NSYMBT = %8d (# bytes storing symmetry operators)\n",nsymbt);
printf(" LSKFLG = %8d (skew matrix flag: 0:none, 1:follows)\n",lskflg);*/
if
(
header_ok
==
0
)
{
cout
<<
"ERROR reading MRC header: "
<<
filename
;
exit
(
1
);
}
if
(
nr
!=
param
.
param_device
.
NumberPixels
||
nc
!=
param
.
param_device
.
NumberPixels
)
{
cout
<<
"PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( "
<<
param
.
param_device
.
NumberPixels
<<
", i "
<<
nc
<<
", j "
<<
nr
<<
")"
<<
"
\n
"
;
exit
(
1
);
}
if
(
mode
!=
2
)
{
cout
<<
"ERROR with MRC mode "
<<
mode
<<
"
\n
"
;
cout
<<
"Currently mode 2 is the only one allowed"
<<
"
\n
"
;
exit
(
1
);
}
else
{
rewind
(
fin
);
for
(
count
=
0
;
count
<
256
;
++
count
)
if
(
read_float_empty
(
fin
)
==
0
)
{
cout
<<
"ERROR Converting Data: "
<<
filename
;
exit
(
1
);
}
for
(
count
=
0
;
count
<
(
unsigned
long
)
nsymbt
;
++
count
)
if
(
read_char_float
(
&
currfloat
,
fin
)
==
0
)
{
cout
<<
"ERROR Converting Data: "
<<
filename
;
exit
(
1
);
}
for
(
int
nmap
=
0
;
nmap
<
ns
;
nmap
++
)
{
st
=
0.0
;
st2
=
0.0
;
for
(
int
j
=
0
;
j
<
nr
;
j
++
)
for
(
int
i
=
0
;
i
<
nc
;
i
++
)
{
if
(
read_float
(
&
currfloat
,
fin
,
swap
)
==
0
)
{
cout
<<
"ERROR Converting Data: "
<<
filename
;
exit
(
1
);
}
else
{
Ref
[
nmap
+
ntotRefMap
].
points
[
i
][
j
]
=
currfloat
;
st
+=
currfloat
;
st2
+=
currfloat
*
currfloat
;
}
}
//Normaling maps to zero mean and unit standard deviation
st
/=
float
(
nr
*
nc
);
st2
=
sqrt
(
st2
/
float
(
nr
*
nc
)
-
st
*
st
);
for
(
int
j
=
0
;
j
<
nr
;
j
++
)
for
(
int
i
=
0
;
i
<
nc
;
i
++
){
Ref
[
nmap
+
ntotRefMap
].
points
[
i
][
j
]
=
Ref
[
nmap
+
ntotRefMap
].
points
[
i
][
j
]
/
st2
-
st
/
st2
;
// if(nmap+ntotRefMap==300)cout << i << " " << j << " " << nmap+ntotRefMap << " " << Ref[nmap+ntotRefMap].points[i][j] << "\n";
}
}
ntotRefMap
+=
ns
;
// cout << ntotRefMap << "\n";
}
fclose
(
fin
);
return
(
0
);
}
int
bioem_RefMap
::
read_float
(
float
*
currfloat
,
FILE
*
fin
,
int
swap
)
{
unsigned
char
*
cptr
,
tmp
;
if
(
fread
(
currfloat
,
4
,
1
,
fin
)
!=
1
)
return
0
;
if
(
swap
==
1
)
{
cptr
=
(
unsigned
char
*
)
currfloat
;
tmp
=
cptr
[
0
];
cptr
[
0
]
=
cptr
[
3
];
cptr
[
3
]
=
tmp
;
tmp
=
cptr
[
1
];
cptr
[
1
]
=
cptr
[
2
];
cptr
[
2
]
=
tmp
;
}
return
1
;
}
int
bioem_RefMap
::
read_int
(
int
*
currlong
,
FILE
*
fin
,
int
swap
)
{
unsigned
char
*
cptr
,
tmp
;
if
(
fread
(
currlong
,
4
,
1
,
fin
)
!=
1
)
return
0
;
if
(
swap
==
1
)
{
cptr
=
(
unsigned
char
*
)
currlong
;
tmp
=
cptr
[
0
];
cptr
[
0
]
=
cptr
[
3
];
cptr
[
3
]
=
tmp
;
tmp
=
cptr
[
1
];
cptr
[
1
]
=
cptr
[
2
];
cptr
[
2
]
=
tmp
;
}
return
1
;
}
int
bioem_RefMap
::
read_float_empty
(
FILE
*
fin
)
{
float
currfloat
;
if
(
fread
(
&
currfloat
,
4
,
1
,
fin
)
!=
1
)
return
0
;
return
1
;
}
int
bioem_RefMap
::
read_char_float
(
float
*
currfloat
,
FILE
*
fin
)
{
char
currchar
;
if
(
fread
(
&
currchar
,
1
,
1
,
fin
)
!=
1
)
return
0
;
*
currfloat
=
(
float
)
currchar
;
return
1
;
}
int
bioem_RefMap
::
test_mrc
(
const
char
*
vol_file
,
int
swap
)
{
FILE
*
fin
;
int
nc
,
nr
,
ns
,
mx
,
my
,
mz
;
int
mode
,
ncstart
,
nrstart
,
nsstart
;
float
xlen
,
ylen
,
zlen
;
int
i
,
header_ok
=
1
,
n_range_viols
=
0
;
int
mapc
,
mapr
,
maps
;
float
alpha
,
beta
,
gamma
;
float
dmin
,
dmax
,
dmean
,
dummy
,
xorigin
,
yorigin
,
zorigin
;
fin
=
fopen
(
vol_file
,
"rb"
);
if
(
fin
==
NULL
)
{
cout
<<
"ERROR opening MRC: "
<<
vol_file
;
exit
(
1
);
}
//* read header info
header_ok
*=
read_int
(
&
nc
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nr
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ns
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mode
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ncstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nrstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nsstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mx
,
fin
,
swap
);
header_ok
*=
read_int
(
&
my
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mz
,
fin
,
swap
);
header_ok
*=
read_float
(
&
xlen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
ylen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
zlen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
alpha
,
fin
,
swap
);
header_ok
*=
read_float
(
&
beta
,
fin
,
swap
);
header_ok
*=
read_float
(
&
gamma
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mapc
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mapr
,
fin
,
swap
);
header_ok
*=
read_int
(
&
maps
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmin
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmax
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmean
,
fin
,
swap
);
for
(
i
=
23
;
i
<
50
;
++
i
)
header_ok
*=
read_float
(
&
dummy
,
fin
,
swap
);
header_ok
*=
read_float
(
&
xorigin
,
fin
,
swap
);
header_ok
*=
read_float
(
&
yorigin
,
fin
,
swap
);
header_ok
*=
read_float
(
&
zorigin
,
fin
,
swap
);
fclose
(
fin
);
if
(
header_ok
==
0
)
{
cout
<<
"ERROR reading MRC header: "
<<
vol_file
;
exit
(
1
);
}
n_range_viols
+=
(
nc
>
5000
);
n_range_viols
+=
(
nc
<
0
);
n_range_viols
+=
(
nr
>
5000
);
n_range_viols
+=
(
nr
<
0
);
n_range_viols
+=
(
ns
>
5000
);
n_range_viols
+=
(
ns
<
0
);
n_range_viols
+=
(
ncstart
>
5000
);
n_range_viols
+=
(
ncstart
<-
5000
);
n_range_viols
+=
(
nrstart
>
5000
);
n_range_viols
+=
(
nrstart
<-
5000
);
n_range_viols
+=
(
nsstart
>
5000
);
n_range_viols
+=
(
nsstart
<-
5000
);
n_range_viols
+=
(
mx
>
5000
);
n_range_viols
+=
(
mx
<
0
);
n_range_viols
+=
(
my
>
5000
);
n_range_viols
+=
(
my
<
0
);
n_range_viols
+=
(
mz
>
5000
);
n_range_viols
+=
(
mz
<
0
);
n_range_viols
+=
(
alpha
>
360.0
f
);
n_range_viols
+=
(
alpha
<-
360.0
f
);
n_range_viols
+=
(
beta
>
360.0
f
);
n_range_viols
+=
(
beta
<-
360.0
f
);
n_range_viols
+=
(
gamma
>
360.0
f
);
n_range_viols
+=
(
gamma
<-
360.0
f
);
return
n_range_viols
;
}
int
bioem_RefMap
::
read_MRC
(
const
char
*
filename
,
bioem_param
&
param
)
{
myfloat_t
st
,
st2
;
unsigned
long
count
;
FILE
*
fin
;
float
currfloat
;
int
nc
,
nr
,
ns
,
swap
,
header_ok
=
1
;
float
xlen
,
ylen
,
zlen
;
int
mode
,
ncstart
,
nrstart
,
nsstart
,
ispg
,
nsymbt
,
lskflg
;
float
a_tmp
,
b_tmp
,
g_tmp
;
int
mx
,
my
,
mz
,
mapc
,
mapr
,
maps
;
float
dmin
,
dmax
,
dmean
;
int
n_range_viol0
,
n_range_viol1
;
fin
=
fopen
(
filename
,
"rb"
);
if
(
fin
==
NULL
)
{
cout
<<
"ERROR opening MRC: "
<<
filename
;
exit
(
1
);
}
n_range_viol0
=
test_mrc
(
filename
,
0
);
n_range_viol1
=
test_mrc
(
filename
,
1
);
if
(
n_range_viol0
<
n_range_viol1
)
{
//* guess endianism
swap
=
0
;
if
(
n_range_viol0
>
0
)
{
printf
(
" Warning: %i header field range violations detected in file %s
\n
"
,
n_range_viol0
,
filename
);
}
}
else
{
swap
=
1
;
if
(
n_range_viol1
>
0
)
{
printf
(
"Warning: %i header field range violations detected in file %s
\n
"
,
n_range_viol1
,
filename
);
}
}
printf
(
"
\n
+++++++++++++++++++++++++++++++++++++++++++
\n
"
);
printf
(
"Reading Information from MRC: %s
\n
"
,
filename
);
header_ok
*=
read_int
(
&
nc
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nr
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ns
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mode
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ncstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nrstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nsstart
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mx
,
fin
,
swap
);
header_ok
*=
read_int
(
&
my
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mz
,
fin
,
swap
);
header_ok
*=
read_float
(
&
xlen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
ylen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
zlen
,
fin
,
swap
);
header_ok
*=
read_float
(
&
a_tmp
,
fin
,
swap
);
header_ok
*=
read_float
(
&
b_tmp
,
fin
,
swap
);
header_ok
*=
read_float
(
&
g_tmp
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mapc
,
fin
,
swap
);
header_ok
*=
read_int
(
&
mapr
,
fin
,
swap
);
header_ok
*=
read_int
(
&
maps
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmin
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmax
,
fin
,
swap
);
header_ok
*=
read_float
(
&
dmean
,
fin
,
swap
);
header_ok
*=
read_int
(
&
ispg
,
fin
,
swap
);
header_ok
*=
read_int
(
&
nsymbt
,
fin
,
swap
);
header_ok
*=
read_int
(
&
lskflg
,
fin
,
swap
);
printf
(
"Number Columns = %8d
\n
"
,
nc
);
printf
(
"Number Rows = %8d
\n
"
,
nr
);
printf
(
"Number Sections = %8d
\n
"
,
ns
);
printf
(
"MODE = %4d (only data type mode 2: 32-bit)
\n
"
,
mode
);
printf
(
"NSYMBT = %4d (# bytes symmetry operators)
\n
"
,
nsymbt
);
/* printf(" NCSTART = %8d (index of first column, counting from 0)\n",ncstart);
printf("> NRSTART = %8d (index of first row, counting from 0)\n",nrstart);
printf(" NSSTART = %8d (index of first section, counting from 0)\n",nsstart);
printf(" MX = %8d (# of X intervals in unit cell)\n",mx);
printf(" MY = %8d (# of Y intervals in unit cell)\n",my);
printf(" MZ = %8d (# of Z intervals in unit cell)\n",mz);
printf(" X length = %8.3f (unit cell dimension)\n",xlen);
printf(" Y length = %8.3f (unit cell dimension)\n",ylen);
printf(" Z length = %8.3f (unit cell dimension)\n",zlen);
printf(" Alpha = %8.3f (unit cell angle)\n",a_tmp);
printf(" Beta = %8.3f (unit cell angle)\n",b_tmp);
printf(" Gamma = %8.3f (unit cell angle)\n",g_tmp);
printf(" MAPC = %8d (columns axis: 1=X,2=Y,3=Z)\n",mapc);
printf(" MAPR = %8d (rows axis: 1=X,2=Y,3=Z)\n",mapr);
printf(" MAPS = %8d (sections axis: 1=X,2=Y,3=Z)\n",maps);
printf(" DMIN = %8.3f (minimum density value - ignored)\n",dmin);
printf(" DMAX = %8.3f (maximum density value - ignored)\n",dmax);
printf(" DMEAN = %8.3f (mean density value - ignored)\n",dmean);
printf(" ISPG = %8d (space group number - ignored)\n",ispg);
printf(" NSYMBT = %8d (# bytes storing symmetry operators)\n",nsymbt);
printf(" LSKFLG = %8d (skew matrix flag: 0:none, 1:follows)\n",lskflg);*/
if
(
header_ok
==
0
)
{
cout
<<
"ERROR reading MRC header: "
<<
filename
;
exit
(
1
);
}
if
(
nr
!=
param
.
param_device
.
NumberPixels
||
nc
!=
param
.
param_device
.
NumberPixels
)
{
cout
<<
"PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( "
<<
param
.
param_device
.
NumberPixels
<<
", i "
<<
nc
<<
", j "
<<
nr
<<
")"
<<
"
\n
"
;
exit
(
1
);
}
if
(
mode
!=
2
)
{
cout
<<
"ERROR with MRC mode "
<<
mode
<<
"
\n
"
;
cout
<<
"Currently mode 2 is the only one allowed"
<<
"
\n
"
;
exit
(
1
);
}
else
{
rewind
(
fin
);
for
(
count
=
0
;
count
<
256
;
++
count
)
if
(
read_float_empty
(
fin
)
==
0
)
{
cout
<<
"ERROR Converting Data: "
<<
filename
;
exit
(
1
);
}
for
(
count
=
0
;
count
<
(
unsigned
long
)
nsymbt
;
++
count
)
if
(
read_char_float
(
&
currfloat
,
fin
)
==
0
)
{
cout
<<
"ERROR Converting Data: "
<<
filename
;
exit
(
1
);
}
for
(
int
nmap
=
0
;
nmap
<
ns
;
nmap
++
)
{
st
=
0.0
;
st2
=
0.0
;
for
(
int
j
=
0
;
j
<
nr
;
j
++
)
for
(
int
i
=
0
;
i
<
nc
;
i
++
)