00001 #include "vuFourierVolume_IO.h"
00002 #include "Volume/Fourier/Unimodal/3d/3d.h"
00003
00004 #define MAX_FILTER_SIZE 6
00005
00006 template <int S> bool vuFourierVolume_IO<S>::
00007 preprocessSpatialInput(byte *data, dword XSize, dword YSize, dword ZSize,
00008 float s, float mult_pad, dword add_pad,
00009 bool doWrapAndInit)
00010 {
00011 dword dim = computeDimensions(XSize, YSize, ZSize, mult_pad, add_pad);
00012
00013 m_XSize = m_YSize = m_ZSize = dim;
00014 m_Volume = new float[m_XSize * m_YSize * m_ZSize * 2 * S];
00015
00016
00017 readSpatial(data, m_Volume, m_XSize, m_YSize, m_ZSize, XSize, YSize, ZSize);
00018
00019
00020 preprocess();
00021
00022 if (doWrapAndInit) {
00023 setWrap(MAX_FILTER_SIZE/2);
00024
00025 XSize = m_XSize;
00026 YSize = m_YSize;
00027 ZSize = m_ZSize;
00028
00029 m_XSize = XSize + m_Wrap * 2;
00030 m_YSize = YSize + m_Wrap * 2;
00031 m_ZSize = ZSize + m_Wrap * 2;
00032
00033 float *dest = new float[m_XSize * m_YSize * m_ZSize * 2 * S];
00034
00035
00036 padFourier(m_Volume, dest, m_XSize, m_YSize, m_ZSize, XSize, YSize, ZSize);
00037 CHECKNDELETE(m_Volume);
00038 m_Volume = dest;
00039
00040 wrapAndInitialize(s);
00041 }
00042 else
00043 m_Wrap = 0;
00044
00045 return true;
00046 }
00047
00048 template <int S>
00049 void vuFourierVolume_IO<S>::preprocess()
00050 {
00051 cout << "shifting volume. " << flush;
00052 shift3D(m_Volume, m_XSize, m_YSize, m_ZSize);
00053
00054 cout << "transforming volume. " << flush;
00055 transform3D(m_Volume, m_XSize, m_YSize, m_ZSize);
00056
00057 cout << "shifting volume. " << flush;
00058 shift3D(m_Volume, m_XSize, m_YSize, m_ZSize);
00059 cout << "done\n" << flush;
00060 }
00061
00062 template <int S> void vuFourierVolume_IO<S>::
00063 transform3D(float* vol, dword xx, dword yy, dword zz)
00064 {
00065 fftwnd_plan plan;
00066
00067 plan = fftw3d_create_plan(xx, yy, zz, FFTW_FORWARD,
00068 FFTW_ESTIMATE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00069
00070 fftwnd(plan, S, (fftw_complex*)vol, S, 1, NULL, 0, 0);
00071
00072 fftwnd_destroy_plan(plan);
00073 dword count = xx * yy * zz * 2 * S;
00074 float size_f = (float) xx* yy * zz;
00075 float *ptr = vol;
00076 for (dword i = 0; i < count; ++i) *(ptr++) /= size_f;
00077 }
00078
00079 template <int S> void vuFourierVolume_IO<S>::
00080 inverseTransform3D(float* vol, dword xx, dword yy, dword zz)
00081 {
00082 fftwnd_plan plan;
00083
00084 plan = fftw3d_create_plan(xx, yy, zz, FFTW_BACKWARD,
00085 FFTW_ESTIMATE | FFTW_IN_PLACE | FFTW_USE_WISDOM);
00086
00087 fftwnd(plan, S, (fftw_complex*)vol, S, 1, NULL, 0, 0);
00088
00089 fftwnd_destroy_plan(plan);
00090 }
00091
00092
00093 template <int S> void vuFourierVolume_IO<S>::
00094 shift3D(float* data,dword XSize,dword YSize,dword ZSize)
00095 {
00096 dword xstep = XSize * 2 * S;
00097 dword ystep = YSize * XSize * 2 * S;
00098 dword index1 = 0;
00099 int start1 = 0;
00100 int val1 = -1;
00101
00102 dword index2;
00103 int val2, start2;
00104
00105 for(dword k=0; k<ZSize; k++) {
00106 index2 = index1;
00107 start2 = start1;
00108 val2 = val1;
00109 for(dword j=0; j<YSize; j++) {
00110 for(dword i=start2; i<xstep; i+=4*S) {
00111 float *ptr = &data[index2+i];
00112 for (int l=0; l<S; l++) {
00113 *(ptr++) *= -1.0f;
00114 *(ptr++) *= -1.0f;
00115 }
00116 }
00117 val2 *= -1;
00118 start2 += val2 * 2 * S;
00119 index2 += xstep;
00120 }
00121 val1 *= -1;
00122 start1 += val1*2*S;
00123 index1 += ystep;
00124 }
00125 }
00126
00127 template <int S> bool vuFourierVolume_IO<S>::
00128 scaleAndWriteToFourierFile(const char *fileName, float scale)
00129 {
00130 if (scale !=1.0 && scale > 0 || true) {
00131 dword count = m_XSize * m_YSize * m_ZSize * 2 * S;
00132 float *ptr = m_Volume;
00133 for (dword i=0; i<count; i++) *(ptr++) *= scale;
00134 }
00135
00136 FILE *file = fopen(fileName, "wb");
00137 if (file == NULL) return false;
00138
00139 dword XSize = m_XSize - 2 * m_Wrap;
00140 dword YSize = m_YSize - 2 * m_Wrap;
00141 dword ZSize = m_ZSize - 2 * m_Wrap;
00142
00143 bool status = vu1712<S>::writeHeader(file, XSize, YSize, ZSize);
00144 if (status) {
00145 status = vu1712<S>::writeData(file, m_Volume,
00146 m_XSize, m_YSize, m_ZSize,
00147 XSize, YSize, ZSize);
00148 }
00149 fclose(file);
00150 return status;
00151 }
00152
00153
00154
00155
00156
00157 template <int S>
00158 void vuFourierVolume_IO<S>::readSpatial(byte *in, float* out,
00159 dword XX, dword YY, dword ZZ,
00160 dword XXsmall, dword YYsmall, dword ZZsmall)
00161 {
00162 byte *i_ptr = in;
00163 float *v_ptr = out;
00164
00165 dword lowX = (XX - XXsmall) / 2;
00166 dword lowY = (YY - YYsmall) / 2;
00167 dword lowZ = (ZZ - ZZsmall) / 2;
00168
00169 dword hiX = XX - lowX - XXsmall;
00170 dword hiY = YY - lowY - YYsmall;
00171 dword hiZ = ZZ - lowZ - ZZsmall;
00172
00173 v_ptr = pad(v_ptr, lowZ * XX * YY * 2 * S);
00174 for (dword k = 0; k < ZZsmall; k++) {
00175 v_ptr = pad(v_ptr, lowY * XX * 2 * S);
00176 for (dword j = 0; j < YYsmall; j++) {
00177 v_ptr = pad(v_ptr, lowX * 2 * S);
00178 for (dword i = 0; i < XXsmall; i++) {
00179 for (int l=0; l<S; l++) {
00180 *(v_ptr++) = *(i_ptr++);
00181 *(v_ptr++) = 0.0f;
00182 }
00183 }
00184 v_ptr = pad(v_ptr, hiX * 2 * S);
00185 }
00186 v_ptr = pad(v_ptr, hiY * XX * 2 * S);
00187 }
00188 v_ptr = pad(v_ptr, hiZ * XX * YY * 2 * S);
00189 }
00190
00191
00192 template <int S>
00193 float* vuFourierVolume_IO<S>::pad(float* v, dword n)
00194 {
00195 for (dword i = 0; i < n; ++i) *(v++) = 0.0f;
00196 return v;
00197 }
00198
00199 template <int S>
00200 void vuFourierVolume_IO<S>::padFourier(float *in, float* out,
00201 dword XX, dword YY, dword ZZ,
00202 dword XXsmall, dword YYsmall, dword ZZsmall)
00203 {
00204 float *i_ptr = in;
00205 float *v_ptr = out;
00206
00207 dword lowX = (XX - XXsmall) / 2;
00208 dword lowY = (YY - YYsmall) / 2;
00209 dword lowZ = (ZZ - ZZsmall) / 2;
00210
00211 dword hiX = XX - lowX - XXsmall;
00212 dword hiY = YY - lowY - YYsmall;
00213 dword hiZ = ZZ - lowZ - ZZsmall;
00214
00215 v_ptr = pad(v_ptr, lowZ * XX * YY * 2 * S);
00216 for (dword k = 0; k < ZZsmall; k++) {
00217 v_ptr = pad(v_ptr, lowY * XX * 2 * S);
00218 for (dword j = 0; j < YYsmall; j++) {
00219 v_ptr = pad(v_ptr, lowX * 2 * S);
00220 for (dword i = 0; i < XXsmall; i++) {
00221 for (int l=0; l<S; l++) {
00222 *(v_ptr++) = *(i_ptr++);
00223 *(v_ptr++) = *(i_ptr++);
00224 }
00225 }
00226 v_ptr = pad(v_ptr, hiX * 2 * S);
00227 }
00228 v_ptr = pad(v_ptr, hiY * XX * 2 * S);
00229 }
00230 v_ptr = pad(v_ptr, hiZ * XX * YY * 2 * S);
00231 }
00232
00233 template <int S>
00234 bool vuFourierVolume_IO<S>::getSpatialDataFromVUF(byte *&spatVol,
00235 dword &XSize,
00236 dword &YSize,
00237 dword &ZSize,
00238 float scale,
00239 vuString fileName)
00240 {
00241 CHECKNDELETE(spatVol);
00242
00243 if (S != 1) {
00244 cerr << "vuFourierVolume_IO::getSpatialDataFromVUF(): ";
00245 cerr << "only S=1 is supported" << endl;
00246 return false;
00247 }
00248
00249 float *freqVol = NULL;
00250 bool status = true;
00251 if (fileName.isEmpty()) return false;
00252 FILE *file = fopen(fileName, "rb");
00253 if (file == NULL) return false;
00254
00255 status = vu1712<S>::readHeader(file, XSize, YSize, ZSize);
00256
00257 if (status) {
00258 freqVol = new float[XSize * YSize * ZSize * S * 2];
00259 status = vu1712<S>::readData(file, freqVol,
00260 XSize, YSize, ZSize,
00261 XSize, YSize, ZSize);
00262 }
00263 fclose(file);
00264
00265 if (status) {
00266 vuFourierVolume_IO<S>::inverseTransform3D(freqVol, XSize, YSize, ZSize);
00267
00268 dword count = XSize * YSize * ZSize;
00269 float *s_ptr = freqVol;
00270 spatVol = new byte[count];
00271 byte *d_ptr = spatVol;
00272
00273 if (scale == 0) scale = 1.0;
00274
00275 for (dword i=0; i<count; i++) {
00276 register float tmp = *s_ptr * scale;
00277
00278 if (tmp > 255)
00279 tmp = 255.0f;
00280 else if (tmp < 0)
00281 tmp = 0.0f;
00282
00283 *d_ptr = (byte)tmp;
00284
00285 d_ptr += 1;
00286 s_ptr += 2;
00287 }
00288 }
00289
00290 CHECKNDELETE(freqVol);
00291 return status;
00292 }