00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "mrilib.h"
00034 #include "r_misc.h"
00035 #include "r_new_resam_dset.h"
00036
00037
00038
00039 static int apply_dataxes ( THD_3dim_dataset * dset, THD_dataxes * dax );
00040 static int apply_orientation ( THD_3dim_dataset * dset, FD_brick * fdb,
00041 char orient[] );
00042 static int valid_resam_inputs( THD_3dim_dataset * , THD_3dim_dataset *,
00043 double, double, double, char [], int );
00044
00045 #define LR_RESAM_IN_FAIL -1
00046 #define LR_RESAM_IN_NONE 0
00047 #define LR_RESAM_IN_REORIENT 1
00048 #define LR_RESAM_IN_RESAM 2
00049
00050 static char * this_file = "r_new_resam_dset.c";
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073 THD_3dim_dataset * r_new_resam_dset
00074 (
00075 THD_3dim_dataset * din,
00076 THD_3dim_dataset * min,
00077 double dx,
00078 double dy,
00079 double dz,
00080 char orient [],
00081 int resam,
00082 int * sublist
00083 )
00084 {
00085 THD_3dim_dataset * dout;
00086 THD_3dim_dataset * dtmp;
00087 THD_dataxes new_daxes;
00088 FD_brick * fdb = NULL;
00089 int work;
00090 char * l_orient, l_orient_s[4] = "";
00091
00092
00093 if ( min ) l_orient = l_orient_s;
00094 else l_orient = orient;
00095
00096 work = valid_resam_inputs( din, min, dx, dy, dz, l_orient, resam );
00097
00098 if ( work == LR_RESAM_IN_FAIL )
00099 return NULL;
00100
00101 if ( sublist )
00102 dtmp = THD_copy_dset_subs(din, sublist);
00103 else
00104 dtmp = din;
00105
00106 dout = EDIT_empty_copy( dtmp );
00107 if( ! ISVALID_3DIM_DATASET( dout ) )
00108 {
00109 fprintf( stderr, "ERROR: <%s> - failed to duplicate datset at %p\n",
00110 this_file, din );
00111 return NULL;
00112 }
00113
00114
00115 EDIT_dset_items(dout, ADN_prefix, "junk_resample", ADN_none );
00116
00117
00118 if ( min &&
00119 (dout->view_type != min->view_type) &&
00120 (dout->dblk && dout->dblk->diskptr) )
00121 {
00122 dout->view_type = min->view_type;
00123 THD_init_diskptr_names( dout->dblk->diskptr, NULL, NULL, NULL,
00124 min->view_type, True );
00125 }
00126
00127
00128
00129 if ( work & LR_RESAM_IN_REORIENT )
00130 {
00131 if ( ( fdb = THD_oriented_brick( dout, l_orient ) ) == NULL )
00132 {
00133 fprintf( stderr, "<%s>: failed to create THD_brick with orient "
00134 "string <%.6s>\n", this_file, l_orient );
00135 THD_delete_3dim_dataset( dout, FALSE );
00136 return NULL;
00137 }
00138
00139 if ( apply_orientation( dout, fdb, l_orient ) != 0 )
00140 {
00141 THD_delete_3dim_dataset( dout, FALSE );
00142 return NULL;
00143 }
00144 }
00145
00146 new_daxes = *dout->daxes;
00147
00148
00149 if ( work & LR_RESAM_IN_RESAM )
00150 {
00151
00152 new_daxes.type = DATAXES_TYPE;
00153
00154
00155 if ( min && dx == 0.0 ) new_daxes = *min->daxes;
00156 else
00157 {
00158
00159 if ( min ) *dout->daxes = *min->daxes;
00160
00161 if ( r_dxyz_mod_dataxes(dx, dy, dz, dout->daxes, &new_daxes) != 0 )
00162 {
00163 THD_delete_3dim_dataset( dout, FALSE );
00164 return NULL;
00165 }
00166 }
00167
00168 if ( apply_dataxes( dout, &new_daxes ) != 0 )
00169 {
00170 THD_delete_3dim_dataset( dout, FALSE );
00171 return NULL;
00172 }
00173 }
00174
00175
00176 dout->warp_parent = dtmp;
00177 dout->warp_parent_idcode = dtmp->idcode;
00178
00179
00180 dout->wod_flag = True;
00181 dout->vox_warp->type = ILLEGAL_TYPE;
00182 *dout->wod_daxes = new_daxes;
00183
00184 dout->daxes->parent = (XtPointer)dout;
00185
00186 if ( r_fill_resampled_data_brick( dout, resam ) != 0 )
00187 {
00188 THD_delete_3dim_dataset( dout, FALSE );
00189 dout = NULL;
00190 }
00191
00192 if ( sublist ) DSET_delete(dtmp);
00193
00194 return dout;
00195 }
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 int r_fill_resampled_data_brick( THD_3dim_dataset * dset, int resam )
00206 {
00207 THD_diskptr * diskptr = dset->dblk->diskptr;
00208 MRI_IMAGE * im;
00209 char * newdata, * dptr;
00210 float bfac;
00211 int ival, dsize;
00212 int nx, ny, nz, nxy, nxyz, nv;
00213 int slice;
00214
00215 if ( DSET_LOADED(dset) )
00216 {
00217 fprintf( stderr, "error <%s>: trying to fill pre-loaded dataset\n",
00218 this_file );
00219 return FAIL;
00220 }
00221
00222 DSET_lock(dset);
00223
00224
00225 dset->warp = myXtNew( THD_warp );
00226 *dset->warp = IDENTITY_WARP;
00227 ADDTO_KILL( dset->kl, dset->warp );
00228
00229 diskptr->byte_order = mri_short_order();
00230 diskptr->storage_mode = STORAGE_BY_BRICK;
00231
00232
00233 nx = dset->daxes->nxx; ny = dset->daxes->nyy; nz = dset->daxes->nzz;
00234 nv = diskptr->nvals;
00235
00236 nxy = nx * ny;
00237 nxyz = nx * ny * nz;
00238
00239
00240 for ( ival = 0; ival < nv; ival++ )
00241 {
00242
00243 dsize = mri_datum_size( DSET_BRICK_TYPE(dset, ival) );
00244
00245 if ( (newdata = (char *)malloc( nxyz * dsize )) == NULL )
00246 {
00247 fprintf( stderr, "r frdb: alloc failure: %d bytes!\n",
00248 nxyz * dsize );
00249 return FAIL;
00250 }
00251
00252 dptr = newdata;
00253
00254
00255 bfac = DBLK_BRICK_FACTOR(dset->dblk,ival);
00256 DBLK_BRICK_FACTOR(dset->dblk,ival) = 0.0;
00257
00258
00259 for ( slice = 0; slice < nz; slice++)
00260 {
00261 im = AFNI_dataset_slice( dset, 3, slice, ival, resam );
00262 if ( im == NULL )
00263 {
00264 fprintf( stderr, "r_fill_resampled_data_brick: failure to "
00265 "compute dataset slice %d\n", slice );
00266 free( newdata );
00267 return FAIL;
00268 }
00269
00270 memcpy( (void *)dptr, mri_data_pointer(im), nxy*dsize );
00271 mri_free( im );
00272 dptr += nxy*dsize;
00273 }
00274
00275 DBLK_BRICK_FACTOR(dset->dblk,ival) = bfac;
00276
00277 EDIT_substitute_brick(dset, ival, DSET_BRICK_TYPE(dset,ival),
00278 (void *)newdata);
00279 }
00280
00281 dset->dblk->malloc_type = DATABLOCK_MEM_MALLOC;
00282 dset->wod_flag = False;
00283
00284
00285 THD_load_statistics( dset );
00286
00287 return 0;
00288 }
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298 int r_dxyz_mod_dataxes( double dx, double dy, double dz,
00299 THD_dataxes * daxin, THD_dataxes * daxout
00300 )
00301 {
00302 double rex, rey, rez;
00303 double lxx, lyy, lzz;
00304 int ret_val;
00305
00306 if ( ! ISVALID_DATAXES( daxin ) || ! ISVALID_DATAXES( daxout ) )
00307 return -1;
00308
00309 *daxout = *daxin;
00310
00311 if ( dx <= 0.0 || dy <= 0.0 || dz <= 0.0 )
00312 return -1;
00313
00314 rex = (daxout->xxdel > 0) ? dx : -dx;
00315 rey = (daxout->yydel > 0) ? dy : -dy;
00316 rez = (daxout->zzdel > 0) ? dz : -dz;
00317
00318 lxx = daxin->nxx * daxin->xxdel;
00319 lyy = daxin->nyy * daxin->yydel;
00320 lzz = daxin->nzz * daxin->zzdel;
00321
00322 daxout->nxx = (int)( lxx/rex + 0.499 );
00323 daxout->nyy = (int)( lyy/rey + 0.499 );
00324 daxout->nzz = (int)( lzz/rez + 0.499 );
00325
00326
00327 daxout->xxorg = daxin->xxorg + 0.5*(lxx - daxin->xxdel)
00328 - 0.5*(daxout->nxx - 1)*rex;
00329
00330 daxout->yyorg = daxin->yyorg + 0.5*(lyy - daxin->yydel)
00331 - 0.5*(daxout->nyy - 1)*rey;
00332
00333 daxout->zzorg = daxin->zzorg + 0.5*(lzz - daxin->zzdel)
00334 - 0.5*(daxout->nzz - 1)*rez;
00335
00336
00337 daxout->xxdel = rex;
00338 daxout->yydel = rey;
00339 daxout->zzdel = rez;
00340
00341
00342
00343 daxout->xxmin = daxout->xxorg;
00344 daxout->xxmax = daxout->xxorg + (daxout->nxx-1)*daxout->xxdel;
00345 if ( daxout->xxmin > daxout->xxmax )
00346 {
00347 double tmp = daxout->xxmin;
00348 daxout->xxmin = daxout->xxmax;
00349 daxout->xxmax = tmp;
00350 }
00351
00352 daxout->yymin = daxout->yyorg;
00353 daxout->yymax = daxout->yyorg + (daxout->nyy-1)*daxout->yydel;
00354 if ( daxout->yymin > daxout->yymax )
00355 {
00356 double tmp = daxout->yymin;
00357 daxout->yymin = daxout->yymax;
00358 daxout->yymax = tmp;
00359 }
00360
00361 daxout->zzmin = daxout->zzorg;
00362 daxout->zzmax = daxout->zzorg + (daxout->nzz-1)*daxout->zzdel;
00363 if ( daxout->zzmin > daxout->zzmax )
00364 {
00365 double tmp = daxout->zzmin;
00366 daxout->zzmin = daxout->zzmax;
00367 daxout->zzmax = tmp;
00368 }
00369
00370 #ifdef EXTEND_BBOX
00371 daxout->xxmin -= 0.5 * daxout->xxdel;
00372 daxout->xxmax += 0.5 * daxout->xxdel;
00373 daxout->yymin -= 0.5 * daxout->yydel;
00374 daxout->yymax += 0.5 * daxout->yydel;
00375 daxout->zzmin -= 0.5 * daxout->zzdel;
00376 daxout->zzmax += 0.5 * daxout->zzdel;
00377 #endif
00378
00379 return ret_val = 0;
00380 }
00381
00382
00383
00384
00385
00386
00387
00388
00389 Boolean r_is_valid_orient_str ( char ostr [] )
00390 {
00391 int o1, o2, o3;
00392
00393 if ( ostr == NULL )
00394 return False;
00395
00396 o1 = ORCODE(toupper(ostr[0]));
00397 o2 = ORCODE(toupper(ostr[1]));
00398 o3 = ORCODE(toupper(ostr[2]));
00399
00400 if ( ( o1 != ILLEGAL_TYPE ) &&
00401 ( o2 != ILLEGAL_TYPE ) &&
00402 ( o3 != ILLEGAL_TYPE ) &&
00403 OR3OK(o1,o2,o3) )
00404 return True;
00405 else
00406 return False;
00407 }
00408
00409
00410
00411
00412
00413
00414
00415
00416 int r_orient_str2vec ( char ostr [], THD_ivec3 * ovec )
00417 {
00418 int o1, o2, o3;
00419
00420 if ( !ostr || !ovec )
00421 {
00422 fprintf( stderr, "%s: r_orient_str2vec - invalid parameter pair "
00423 "(%p,%p)\n", this_file, ostr, ovec );
00424 return -1;
00425 }
00426
00427 ovec->ijk[0] = o1 = ORCODE(toupper(ostr[0]));
00428 ovec->ijk[1] = o2 = ORCODE(toupper(ostr[1]));
00429 ovec->ijk[2] = o3 = ORCODE(toupper(ostr[2]));
00430
00431 if ( ( o1 == ILLEGAL_TYPE ) ||
00432 ( o2 == ILLEGAL_TYPE ) ||
00433 ( o3 == ILLEGAL_TYPE ) ||
00434 ( !OR3OK(o1,o2,o3) ) )
00435 {
00436 fprintf( stderr, "%s: r_orient_str2vec - bad ostr <%.4s>\n",
00437 this_file, ostr );
00438 return -2;
00439 }
00440
00441 return 0;
00442 }
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456 static int valid_resam_inputs( THD_3dim_dataset * dset, THD_3dim_dataset * mset,
00457 double dx, double dy, double dz,
00458 char orient [], int resam )
00459 {
00460 int ret_val = LR_RESAM_IN_NONE;
00461
00462 if( ! ISVALID_3DIM_DATASET(dset) )
00463 {
00464 fprintf( stderr, "ERROR: <%s> - invalid input dataset\n", this_file );
00465 return LR_RESAM_IN_FAIL;
00466 }
00467
00468
00469 if ( resam < 0 || resam > LAST_RESAM_TYPE )
00470 {
00471 fprintf( stderr, "ERROR: <%s> - invalid resample mode of %d\n",
00472 this_file, resam );
00473 return LR_RESAM_IN_FAIL;
00474 }
00475
00476 if( mset )
00477 {
00478
00479 if ( ! ISVALID_3DIM_DATASET(mset) || !ISVALID_DATAXES(mset->daxes) )
00480 {
00481 fprintf( stderr, "ERROR: <%s> - invalid master dataset\n",
00482 this_file );
00483 return LR_RESAM_IN_FAIL;
00484 }
00485
00486
00487 if ( ! orient )
00488 {
00489 fprintf( stderr, "ERROR: <%s> - orientation memory should "
00490 "come with master\n", this_file );
00491 return LR_RESAM_IN_FAIL;
00492 }
00493
00494
00495 orient[0] = ORIENT_typestr[mset->daxes->xxorient][0];
00496 orient[1] = ORIENT_typestr[mset->daxes->yyorient][0];
00497 orient[2] = ORIENT_typestr[mset->daxes->zzorient][0];
00498 orient[3] = '\0';
00499
00500
00501 if ( dx != 0.0 && (dx < 0.0 || dy <= 0.0 || dz <= 0.0) )
00502 {
00503 fprintf( stderr, "ERROR: <%s> - invalid (dx,dy,dz) = (%f,%f,%f)\n",
00504 this_file, dx, dy, dz );
00505 return LR_RESAM_IN_FAIL;
00506 }
00507
00508 return LR_RESAM_IN_RESAM | LR_RESAM_IN_REORIENT;
00509 }
00510
00511
00512
00513
00514 if ( dx == 0.0 && dy == 0.0 && dz == 0.0 )
00515 ret_val &= ~LR_RESAM_IN_RESAM;
00516 else if ( dx <= 0.0 || dy <= 0.0 || dz <= 0.0 )
00517 {
00518 fprintf( stderr, "ERROR: <%s> - invalid (dx,dy,dz) = (%f,%f,%f)\n",
00519 this_file, dx, dy, dz );
00520 return LR_RESAM_IN_FAIL;
00521 }
00522 else
00523 ret_val |= LR_RESAM_IN_RESAM;
00524
00525
00526 if ( orient == NULL )
00527 ret_val &= ~LR_RESAM_IN_REORIENT;
00528 else if ( r_is_valid_orient_str ( orient ) == False )
00529 {
00530 fprintf( stderr, "ERROR: <%s> - invalid orientation string <%.6s>\n",
00531 this_file, orient );
00532 return LR_RESAM_IN_FAIL;
00533 }
00534 else
00535 ret_val |= LR_RESAM_IN_REORIENT;
00536
00537 return ret_val;
00538 }
00539
00540
00541
00542
00543
00544
00545
00546
00547 static int apply_orientation( THD_3dim_dataset * dset, FD_brick * fdb,
00548 char orient[] )
00549 {
00550 THD_dataxes * daxp = dset->daxes;
00551 THD_ivec3 ivnxyz, ivorient;
00552 THD_fvec3 fvdel, fvorg;
00553 float org4[4], del4[4];
00554 int aa1, aa2, aa3;
00555 int a1, a2, a3;
00556 int ret_val;
00557
00558 LOAD_IVEC3(ivnxyz, fdb->n1, fdb->n2, fdb->n3);
00559
00560 if ( (ret_val = r_orient_str2vec(orient, &ivorient)) != 0 )
00561 return ret_val;
00562
00563 LOAD_FVEC3( fvdel,
00564 ORIENT_sign[ivorient.ijk[0]]=='+' ? fdb->del1 : -fdb->del1,
00565 ORIENT_sign[ivorient.ijk[1]]=='+' ? fdb->del2 : -fdb->del2,
00566 ORIENT_sign[ivorient.ijk[2]]=='+' ? fdb->del3 : -fdb->del3 );
00567
00568
00569 org4[1] = daxp->xxorg; org4[2] = daxp->yyorg; org4[3] = daxp->zzorg;
00570 del4[1] = daxp->xxdel; del4[2] = daxp->yydel; del4[3] = daxp->zzdel;
00571
00572 a1 = fdb->a123.ijk[0]; a2 = fdb->a123.ijk[1]; a3 = fdb->a123.ijk[2];
00573 aa1 = abs(a1); aa2 = abs(a2); aa3 = abs(a3);
00574
00575 LOAD_FVEC3( fvorg,
00576 (a1 > 0) ? org4[aa1] : org4[aa1]+(fdb->n1-1)*del4[aa1],
00577 (a2 > 0) ? org4[aa2] : org4[aa2]+(fdb->n2-1)*del4[aa2],
00578 (a3 > 0) ? org4[aa3] : org4[aa3]+(fdb->n3-1)*del4[aa3] );
00579
00580
00581 EDIT_dset_items( dset,
00582 ADN_nxyz, ivnxyz,
00583 ADN_xyzdel, fvdel,
00584 ADN_xyzorg, fvorg,
00585 ADN_xyzorient, ivorient,
00586 ADN_none );
00587 return 0;
00588 }
00589
00590
00591
00592
00593
00594
00595
00596
00597 static int apply_dataxes( THD_3dim_dataset * dset, THD_dataxes * dax )
00598 {
00599 THD_ivec3 ivnxyz;
00600 THD_fvec3 fvdel, fvorg;
00601
00602 LOAD_IVEC3( ivnxyz, dax->nxx, dax->nyy, dax->nzz );
00603 LOAD_FVEC3( fvorg, dax->xxorg, dax->yyorg, dax->zzorg );
00604 LOAD_FVEC3( fvdel, dax->xxdel, dax->yydel, dax->zzdel );
00605
00606 EDIT_dset_items( dset,
00607 ADN_nxyz, ivnxyz,
00608 ADN_xyzdel, fvdel,
00609 ADN_xyzorg, fvorg,
00610 ADN_none );
00611
00612
00613 *dset->wod_daxes = *dax;
00614 *dset->daxes = *dax;
00615
00616 return 0;
00617 }
00618